# Tri-axial deformation in nuclei with realistic NN interactions

E. N. E. van Dalen and H. Müther  Institut für Theoretische Physik,
Universität Tübingen, D-72076 Tübingen, Germany

###### Abstract

The structure of finite nuclei is investigated by employing an interaction
model which is based on the low-momentum interaction $V_{lowk}$. It is
supplemented by a density-dependent contact interaction fitted to reproduce
the saturation properties of infinite nuclear matter within the Hartree-Fock
approach. The calculations of finite nuclei are performed in a basis of plane
waves discretized in a cartesian box of appropriate size. As a first example
the structure of Ne isotopes is considered ranging from 18Ne to the neutron
drip line. Rather good agreement is obtained for the bulk properties of these
nuclei without any free parameter. The basis is also appropriate to describe
other deformed nuclei and the transition from discrete nuclei to homogeneous
matter which is supposed to occur in the crust of neutron stars.

finite nuclei, neutron star crust, neutron drip line, pasta phase, realistic
NN interaction

###### pacs:

21.60.Jz, 21.30.Fe, 21.65.-f, 26.60.Gj

##  I Introduction

The availability of radioactive ion beam facilities in several laboratories -
like SPIRAL2 at Ganil and the future GSI facility FAIR in Germany - has made
it possible to measure properties of many neutron-rich nuclei far away from
the valley of $\beta$-stability. The study of these exotic nuclear systems is
an interesting field of physics. It is driven by the expectation of observing
nuclear properties, which are different from those encountered so far, i.e.
near the valley of stability, such as special deformations and neutron halos.

The role of deformations has been considered to some extent. Some attention
has been devoted to the consequences for the shell closures far from
stability. Deformations can influence the double $\beta$-decay rate. However,
almost all calculations of the double $\beta$-decay matrix elements performed
so far have been made assuming spherical symmetric nuclei, whereas we know
that in many cases, e.g. 76Ge, they are not. A study of the double
$\beta$-decay of 48Ti reaches the conclusion that a large mismatch of
deformation between the parent and the grand-daughter nuclei can reduce the
value of nuclear matrix elements by factors as large as 2-3 caurier:2008 .
Therefore, an accurate description of the geometrical structure of nuclei is
important, in particular for neutron-rich nuclei far away from the line of
$\beta$-stability. The geometrical structures of these nuclei can be quite
complex compared to deeply-bound, double magic nuclei such as 4He, 16O, and
40Ca. Furthermore, neutron-rich deformed nuclei can be regarded as the
starting point of the transition to infinite homogeneous matter, which should
occur in the outer crust of a neutron star.

Most of the studies on nuclear system far away from the valley of stability
have been performed using density functional theories (DFT). Non-
relativisticbender:2003  as well as covariantniksic:2011  density functional
theories have been developed during the last twenty years and successfully
been applied to describe a large variety of nuclear phenomena with high
precision. The basis of this success, however, are the careful
parameterizations of the underlying density functionals. The variety of
developments and the number of different parameterizations is growing so
rapidly that only a few specialists are able to keep an overview in specific
subfields of DFT.

The alternative would be ab initio calculations, which use a realistic model
of the nucleon-nucleon (NN) interaction, based on the description of the NN
scattering data. However, there are at least four obstacles on the way to
derive properties of finite nuclei in such an ab initio approach.

The first problem arises from the fact that also the realistic NN interaction
is not unique, the number of models, which all fit the NN scattering data is
large, ranging from rather stiff models for local NN interactionsarg1  to
softer One-Boson-Exchange modelsmach1  or the chiral modelschirar ; chira1 ;
chira2 , which have received a lot of attention, recently. The second obstacle
is due to the fact that even the softer version of realistic NN interaction
models contain strong central and tensor components, which make it inevitable
to employ many-body theories well beyond the mean-field approximation (see
e.g. polls00  and references cited there). However, obstacle number three,
even employing sophisticated many-body theories one is in general not able to
reproduce the saturation point of nuclear matter, the simplest nuclear
structure one can imagine. Finally, the effective NN interactions resulting
from such a microscopic approach are rather complicate and are typically
determined in terms of matrix elements using an appropriate basis of single-
particle states. So problem number four is the choice of an appropriate set of
basis states to describe weakly bound nuclei allowing for a general
deformation.

In the following we will address these obstacles for a nuclear structure
calculation a little bit more in detail and present an approach, which is a
kind of compromise: A nuclear structure calculation based on the realistic NN
interaction supplemented by a DFT correction with a minimal number of
adjustable parameters. First applications of this approach in nuclear matter
and deeply bound nuclei have been presented in vandalen:2009b . It is the aim
of the work presented here, to test this approach in the description of weakly
bound nuclei.

A possible solution of the problems one and two of an ab initio approach,
outlined above, is the use of a low-momentum interaction, which can be derived
from a realistic NN interaction by means of renormalization techniques
bogner:2001 ; bogner:2005 ; bogner:2007 ; bozek:2006 ; vandalen:2010 . It is
obtained by separating the low-momentum and high-momentum components of a
realistic NN interaction. This means that all short-range components of the NN
interaction are integrated out. Therefore, a non-perturbative treatment of NN
correlations is not necessary anymore. If the cut-off $\Lambda$ is
appropriately chosen, i.e. around $\Lambda$ = 2 fm-1, the resulting low-
momentum interaction $V_{lowk}$ turns out to be independent of the underlying
realistic interaction. So using $V_{lowk}$ we have a unique NN interaction,
which describes the NN scattering data up to the pion threshold and can be
used in a mean field approach.

However, using the low-momentum interaction $V_{lowk}$ a new problem arises.
Hartree Fock calculations employing $V_{lowk}$ in infinite nuclear matter
produce a binding energy per nucleon increasing with density in a monotonic
way bozek:2006 ; kuckei:2003 . Hence, the emergence of a saturation point is
prevented, the problem number three mentioned even gets worse. While in
Brueckner-Hartree-Fock calculations using realistic NN interactions the
calculated saturation point does not agree with the empirical data, one does
not even get one using $V_{lowk}$.

A way out of this problem is to add three- or many-nucleon forces in the
Hamiltonian. There are various possible sources for such many-nucleon
interactions. The renormalization procedure, which in the case of the two-
nucleon problem leads to the $V_{lowk}$ interaction, gives rise to many-
nucleon forces when applied to systems with three or more nucleons. Such
effective three-nucleon forces, which arise from the renormalization with
respect to high momentum components within a pure nucleonic model, should
provide similar effects as e.g. the density and energy dependence of the
Brueckner G-matrix and therefore support the saturation mechanism. However,
also sub-nucleonic degrees of freedom or relativistic effects may be
represented in terms of a medium dependent effective NN force or three-nucleon
forces.

In the present study we do not care about the origin of such many-nucleon
forces but just accept the fact that a Hamiltonian based on $V_{lowk}$ must be
supplemented to be used in calculations of nuclear structure. This supplement
has been chosen to contain basic contact interaction terms of a non-
relativistic DFT like the Skyrme modelsk1 ; bv81  and adjust the parameters to
reproduce the empirical saturation point of symmetric nuclear matter. This has
been done in vandalen:2009b  and there it could also be demonstrated that this
approach yields very reasonable results for the bulk properties of closed
shell nuclei like 16O and 40Ca.

These double magic nuclei can be described with a spherical geometry using
spherical coordinates coraggio:2005 ; coraggio:2006 . However, non-spherical
nuclei require a geometry allowing for the description of more complex
structures. Therefore, some recent investigations martini:2011 ; ebran:2011
applied cylindrical coordinates to describe axially-symmetric-deformed nuclei.
Although one can describe simple prolate and oblate deformations with the
cylindrical coordinates, the description of more complex tri-axial structures
is not so easy. Therefore, we have chosen a cartesian geometry to allow the
description of tri-axial geometrical structures in these nuclei.

Furthermore, nuclear structure calculations of finite nuclei with realistic
forces employ a set of appropriate single-particle wave functions as a set of
basis states, typically those of a harmonic oscillator (HO) potential. This
choice of the HO basis seems to be quite reasonable if one wants to describe
the structure of deeply-bound closed-shell nuclei. However, the HO basis may
be questionable when weakly-bound nuclei close to the proton or neutron drip
line are to be considered, since it is not appropriate to describe the tail of
the single-particle wave functions for the weakly-bound valence nucleons. For
this kind of studies a basis of plane wave (PW) states may be more suitable
vandalen:2009b .

One aim of this study is to establish a tool, which allows a microscopic
description of tri-axial geometrical structures in nuclei and in the
transition from isolated nuclei to homogeneous matter based on a realistic NN
interaction. Therefore, we will apply a cubic Wigner-Seitz cell using
cartesian coordinates and a tri-axial PW basis. This geometry will allow us to
extend our investigation to tri-axial geometrical structures of nuclei. The
nuclear structure calculations will be performed in the framework of a
Hartree-Fock (HF) calculation using a $V_{lowk}$ potential supplemented with a
contact interaction.

As an example we apply this approach to investigate the density distributions
of neon isotopes. An issue of interest is the situation of N = 20 shell
closure in neon, which lies far away from the valley of $\beta$-stability.
Compared to recent investigations martini:2011 ; ebran:2011 , our calculations
are performed using a tri-axial PW basis and cartesian coordinates in a cubic
Wigner-Seitz cell. This geometry will allow us to extend our investigation to
tri-axial geometrical structures instead of only axially-symmetric
deformations. Furthermore, we use a realistic interaction, a $V_{lowk}$
potential, instead of a phenomenological interaction. In addition, we will go
beyond the description of stable neon isotopes to the situation with unbound
neutrons, which will have some relevance for the transition from discrete
nuclei to homogeneous matter which is supposed to occur, e.g. in the crust of
neutron stars.

The plan of this paper is as follows. The method is treated in Sec. II. The
procedure to renormalize the low-momentum interaction and obtain the matrix
elements of $V_{lowk}$ is reviewed in Sec. II.1. In this section we also
describe our interaction model supplementing $V_{lowk}$ by an appropriate
contact interaction. In Sec. II.2, we work out how to solve the HF equation in
a basis of tri-axial plane wave states in cartesian coordinates. The results
for density distributions of the neon isotopes are presented and discussed in
Sec. III. Finally, Sec. IV contains a summary and the conclusions of our work.

##  II Method

###  II.1 Model of the NN interaction

The main idea of $V_{lowk}$ interaction is to separate the long-range part
from the short-range part of a realistic NN interaction and restrict the
nuclear structure calculation to the low-momentum components, i.e. the long-
range part. This long-range or low-momentum part of the NN interaction is
fairly well described in terms of meson-exchange, whereas the quark degrees of
freedom are important in the description of the short-range or high-momentum
components of the NN interaction. In other words, a model space, which
accounts for the low momentum degrees of freedom, is defined and the effective
Hamiltonian is renormalized for this low-momentum regime to account for the
effects of the high-momentum components, which are integrated out. In our
case, the matrix elements for $V_{lowk}$ are determined by using the unitary-
model-operator approach (UMOA), a model space technique to disentangle low and
high momentum components. This approach has been described in the literature
suzuki:1982 . Therefore, we will restrict ourselves purely to the basics.

First, we define the projection operator $P$, which projects onto the low-
momentum subspace. In our case, i.e. the effective two-nucleon problem, the
operator $P$ projects onto the model space of two-nucleon wave functions with
relative momenta $k$ smaller than a chosen cut-off $\Lambda$. Next, the
projection operator $Q$ is defined. The projection operator $Q$ projects onto
the complement of this subspace, i.e. the high-momentum subspace. In the UMOA
approach, the aim is to define a unitary transformation $U$ in such a way that
the transformed Hamiltonian does not couple $P$ and $Q$ spaces. The technique
to obtain this unitary transformation has been described in Ref. bozek:2006 ;
fuji:2004 . With the help of this unitary transformation, the low momentum
interaction can be obtained. It is defined as

$V_{eff}=U^{-1}(h_{0}+v_{12})U-h_{0},$ (1)

where $v_{12}$ represents the bare NN interaction and the starting Hamiltonian
$h_{0}$ denotes the one-body part of the two-body system. The phase shifts
obtained from solving the Lipmann-Schwinger equation for NN scattering using
this $V_{eff}$ with a cut-off $\Lambda$ are identical to those obtained from
solving the Lipmann-Schwinger equation for the original realistic interaction
$v_{12}$ without a cut-off. If the original interaction is realistic, this
means that this $V_{eff}$ will also be a realistic interaction and reproduces
the NN scattering phase shifts up to a cut-off $\Lambda$. This $V_{eff}$ will
from now on be called $V_{lowk}$.

An interesting feature is that the resulting $V_{lowk}$ is found to be
essentially model independent. In other words, this $V_{lowk}$ potential will
be independent of the underlying realistic interaction $v_{12}$, if the cut-
off $\Lambda$ is chosen around $\Lambda$ = 2 fm-1. In the calculations
discussed in section III below, we actually consider the proton-neutron part
of the CD Bonn potential mach1  for the bare interaction $v_{12}$ using a cut-
off parameter $\Lambda$ = 2 to obtain $V_{lowk}$.

The fact that the $V_{lowk}$ does not induce any short-range correlations into
the nuclear wave function, seems to lead to problems in the description of the
saturation behavior of nuclear matter. Due to this lack of short-range
correlation effects, the emergence of a saturation point is prevented in
calculations of symmetric nuclear matter bogner:2005 ; kuckei:2003 . Hence,
three-body interaction terms or density-dependent two-nucleon interactions are
needed to obtain saturation in isospin symmetric nuclear matter. Therefore,
the effective interaction $V_{lowk}$ is supplemented by a contact interaction,
which we previously presented in Ref. vandalen:2009b . This contact
interaction has the form

$\Delta\mathcal{V}=\Delta\mathcal{V}_{0}+\Delta\mathcal{V}_{3},$ (2)

with

$\Delta\mathcal{V}_{0}=\frac{1}{4}t_{0}[(2+x_{0})\rho^{2}-(2x_{0}+1)(\rho_{n}^{2}+\rho_{p}^{2})]$
(3)

and

$\Delta\mathcal{V}_{3}=\frac{1}{24}t_{3}\rho^{\alpha}[(2+x_{3})\rho^{2}-(2x_{3}+1)(\rho_{n}^{2}+\rho_{p}^{2})],$
(4)

where $\rho_{p}$ and $\rho_{n}$ refer to the local proton and neutron
densities. Furthermore, the matter density is denoted as
$\rho=\rho_{p}+\rho_{n}$. The parameters of the contact interaction are
obtained in the following way: the parameters $\alpha=0.5$ and $x_{0}=0.0$
have been fixed. The three remaining parameters $t_{0}$, $t_{3}$, and $x_{3}$
have then been fitted to reproduce the empirical saturation point of nuclear
matter and the symmetry energy at saturation density. The results for these
fitting parameters are listed in Table 1.

$t_{0}$ [MeV fm3] |  $t_{3}$ [MeV fm3+3α] |  $x_{3}$
---|---|---
584.1 | 8330.7 | -0.5
Table 1:  The fitted parameters $t_{0}$, $t_{3}$ and $x_{3}$ defining the
contact interaction of Eq. (2).

###  II.2 Hartree-Fock Calculations in a tri-axial Plane Wave Basis

While the Hartree-Fock or mean-field contribution arising from the contact
interaction can easily be evaluated in configuration space, the effective
interaction $V_{lowk}$ is non-local and is typically defined in terms of
antisymmetrized two-body matrixelements

$<\alpha\gamma|V_{lowk}|\beta\delta>\,,$ (5)

using a set of single-particle basis states $|\alpha>,\,|\beta>\dots$. This
set of basis states should be appropriate to describe the self-consistent
single-particle states emerging from the Hartree-Fock calculation for deformed
nuclei, but also allow the evaluation of the matrix elements of $V_{lowk}$ in
a straightforward and simple manor.

As we are interested in the description of nuclei with weakly bound valence
nucleons and close to the neutron drip line, we did not use a basis of
harmonic oscillator states but considered a set of plane-wave states with
periodic boundary conditions on the surface of a cubic box with a box size of
$L$ in each cartesian direction, i.e. a volume of $L^{3}$. This means that the
cartesian components of the momenta are discretized and we obtain e.g. for the
$x$-direction

$k_{xj}=\frac{2\pi\,j}{L}\quad\mbox{for}\quad j=0,1,\dots$

and corresponding values for the $y$ and $z$ direction. Therefore each basis
state $\alpha$ is identified by the 3 cartesian components forming the
momentum $\vec{k}$ and the quantum number $s=\pm 1/2$ defining the projection
of the nucleon spin in the directions of the z-axis. Therefore the orthonormal
single-particle basis states can be written in configuration space

$\Phi_{\alpha}(\vec{r})=\frac{1}{\sqrt{L^{3}}}\exp(i\vec{k}\cdot\vec{r})\,\chi_{s}\,,$
(6)

with $\chi_{s}$ denoting the two-component Pauli spinor. The Hartree-Fock
single-particle states are then expanded in this basis

$|\Psi_{n}>=\sum_{\alpha=1}^{N}|\alpha><\alpha|\Psi_{n}>=\sum_{\alpha=1}^{N}c_{n\alpha}|\alpha>\,,$
(7)

and the $V_{lowk}$ part of the Hartree-Fock Hamiltonian can be written in
terms of the two-body matrix elements, of (5)

$<\alpha|H_{lowk}^{HF}|\beta>=\sum_{\gamma,\delta}<\alpha\gamma|V_{lowk}|\beta\delta>\rho_{\gamma\delta}\,,$
(8)

where $\rho_{\gamma\delta}$ represents the single-particle density matrix

$\rho_{\gamma\delta}=\sum_{n=1}^{A}c_{n\gamma}^{*}c_{n\delta}=\sum_{n=1}^{A}<\Psi_{n}|\gamma><\delta|\Psi_{n}>\,,$
(9)

with a sum including all occupied HF states $|\Psi_{n}>$. It is
straightforward to supplement this $V_{lowk}$ part of the Hartree-Fock
Hamiltonian by the corresponding terms for the kinetic energy, the contact
interaction of (2) including the so-called rearrangement terms and the Coulomb
potential for the protons pgreinh ; Bonche:2005 ; Ryssens:2014

$<\alpha|H^{HF}|\beta>=<\alpha|H_{kin}+H_{lowk}^{HF}+\Delta\mathcal{V}^{HF}+V_{Coul}^{HF}|\beta>\,.$
(10)

This HF Hamiltonian is diagonalized

$\sum_{\beta=1}^{N}<\alpha|H^{HF}|\beta>c_{n\beta}=\varepsilon_{n}c_{n\alpha}\,,$
(11)

to obtain the expansion coefficients $c_{n\alpha}$ and the HF single-particle
energies $\varepsilon_{n}$. Since the definition of HF Hamiltonian depends on
the expansion coefficients (see e.g. eq.(8)) the evaluation of $H^{HF}$ and
its diagonalization must be iterated in the usual way until a self-consistent
solution is obtained. We initiate this iteration by approximating the HF
states by appropriate eigenstates of a tri-axial deformed Woods Saxon
potential. We note that the Coulomb contribution in (10) has been calculated
using the so-called Slater approximation for the exchange partpgreinh ;
Ryssens:2014  and that the energy for the spurious center of mass motion has
been subtracted from the energy after the variational calculation.

We still have to specify some details concerning the definition and the
truncation of the PW basis defined in (6) and used in (7). Most of the
calculations discussed below have been performed assuming for the size of the
cartesian box a length $L$ of 18 fm. It turned out that stable results were
obtained if the expansion in (7) was restricted to momenta $\vec{k}$ with
$|\vec{k}|\leq$ 2 fm-1, which leads to a typical number of 739 different
momenta $\vec{k}$ to be included. Including the spin degree of freedom this
implies for the number of basis states $N=1478$ in (7) for neutron states and
the same number of coefficients representing the proton states. The only
symmetry, which has been imposed, was to assume the HF states to be states of
good parity. This symmetry ensures that the center of the density distribution
is in the center of the cartesian box.

For the description of spherical nuclei it is of course easier to use a
spherical basis like e.g. the eigenstates of a free particle moving in a
spherical cavity with radius $R$ as it has been done e.g. in vandalen:2009b .
In this case the wave functions of the orthogonal basis can be separated in a
radial part and an angular part

$\Phi_{iljm}(\vec{r})=R_{il}(r){\cal Y}_{ljm}(\vartheta,\varphi)\,.$ (12)

The angular part ${\cal Y}_{ljm}$ represent the spherical harmonics including
the spin degrees of freedom by coupling the orbital angular momentum $l$ with
the spin to a single-particle angular momentum $j$. For the case of spherical
nuclei one obtains single-particle states with definite angular momentum $j$
and a degeneracy of $2j+1$. The expansion of the HF states is restricted to an
expansion of the radial wave function and typically it is sufficient to
consider of the order of 10 expansion coefficients.

The spherical basis of (12) can also be used to describe deformed nuclei. In
this case one must allow the HF states to include admixtures of different
angular momenta $ljm$ and would obtain again larger dimensions for the
expansion of the HF states. A disadvantage of the spherical basis is the fact
that it does not really allow for a good description of the limit of
homogeneous nuclear matter. This is obvious if one uses a basis of radial
functions (i.e. spherical Bessel functions) with a node at the border of the
cavity $r=R$. A Slater determinant built from such basis states for free
particles leads to a non-uniform density distribution. But also using a basis
with mixed boundary conditions for different partial waves does not lead to a
homogeneous density distribution.

The plane wave basis of (6) constrained to a cartesian box contains the limit
of homogeneous matter in a trivial way. In fact, this limit has been used as
one check of the computer code and we were able to reproduce results of
infinite matter. The HF calculation of 16O and 40Ca starting from a spherical
initial wave function has been another test. In this case the HF iteration
preserves the symmetry of the initial state and also the calculation in the
cartesian basis yield spherical solutions, which are within the numerical
accuracy identical to the solutions obtained in the spherical basis (12).

##  III Results and discussion

The model for the NN interaction discussed above has been developed in
vandalen:2009b  and used to describe the properties of closed shell nuclei 16O
and 40Ca. It is the aim of the investigations presented here to test this
approach for the evaluation of properties of open shell nuclei. Special
attention will be paid to the determination of the neutron drip line and the
formation of non-spherical density distributions. As an example we consider
the structure of the Ne isotopes with even nucleon numbers $A$ ranging from
$A=18$ to $A=32$. We also consider the case of $A=34$, which is beyond the
neutron drip line and an example of $A=42$ with 32 neutrons, which should
partly be evaporated to the continuum.

isotope |  cal. E/A [MeV] |  exp. E/A [MeV]
---|---|---
18Ne | -7.20 | -7.34
20Ne | -7.96 | -8.03
22Ne | -8.03 | -8.08
24Ne | -8.02 | -7.99
26Ne | -7.63 | -7.75
28Ne | -7.31 | -7.39
30Ne | -7.01 | -7.07
32Ne | -6.72 | -6.67
34Ne | -6.17 | -
Table 2: Results for the energies per nucleon using the $V_{low-k}$ plus
contact term interaction model in the tri-axial plane wave basis. The
calculated energies have been corrected by the one-body center-of-mass energy.
The experimental values are taken from Ref. audi:1993 .

The results for the binding energy per nucleon are displayed in table 2 and
compared to the experimental data taken from Ref. audi:1993 . The difference
between the calculated values and the experimental data is for all cases below
3 MeV for the total binding energy and therefore comparable to the deviations
typically obtained in pure phenomenological DFT calculationsebran:2011 . This
is very encouraging noting that the interaction model considered here is based
on a realistic NN interaction supplemented by a phenomenological part with
only three parameters adjusted in the description of nuclear matter.

Also the position of the neutron drip-line is reproduced, as one can see from
a comparison of the total binding energies deduced from this table. The total
binding energy evaluated for 34Ne (34 times 6.17 MeV) is smaller than the
corresponding binding energy of 32Ne.

The position of the drip-line can also be extracted from the spectra of the
single-particle states displayed in the two panels of Fig. 1 for protons and
neutrons, respectively. The HF calculations have been performed without a
constraint to time reversal symmetry. Nevertheless, the resulting single-
particle energies form pairs of almost degenerate states exhibiting
differences up to a few hundred keV. The plot shows the position of the lower
of such two single-particle energies.

Figure 1: (Color online) Single-particle energies for protons (left panel) and
neutrons (right panel) obtained from Hartree-Fock calculations for various Ne
isotopes with even nucleon number $A$. Note that each line represents the
position of two degenerate states. The energies of all occupied and the lowest
unoccupied level are presented. The dashed lines represent the position of the
Fermi energy.

The single-particle energies for the protons decrease in a rather monotonic
way with increasing nucleon number $A$. This reflects the strong attraction of
the proton-neutron interaction leading to a stronger binding for the protons
with an increasing number of the neutrons.

The single-particle energies for the neutrons, displayed in the right panel of
Fig. 1 do not show such a monotonic behavior. Most of the single-particle
energies change only very little with neutron number. This is true in
particular for $A\geq 22$, i.e. the isotopes with more neutrons than the 22Ne,
which is the isotope with largest binding energy per nucleon.

Going from $A$ to $A+2$ only the position of the states, which are not
occupied in the isotope with $A$ nucleons but get occupied in the $A+2$
isotope receive a stronger binding. Inspecting the neutron single-particle
energies, which are close to the Fermi energy in the cases of $A$=32 and
$A$=34, the position of the neutron drip-line discussed above gets confirmed.
While all neutron states occupied in 32Ne have single-particle energies below
0, one pair of states to be occupied in 34Ne has an energy larger than 0.
These states are almost degenerate with the lowest particle states above the
Fermi energy.

Some of these features are also represented by the Fermi energies, which we
define as the arithmetic mean value of the highest occupied state and the
lowest unoccupied state and show in Fig. 2.

Figure 2: (Color online) The Fermi energies for the neutrons and the protons
in neon isotopes.

The Fermi energy for the protons steadily decrease with increasing nucleon
number, while the Fermi energy for the neutrons increase very weakly for
$A\geq$22 and cross the zero at 34Ne.

Another interesting quantity to be deduced from the single-particle spectra is
the energy gap, i.e. the energy difference between the highest occupied state
and the lowest unoccupied state. The corresponding results are shown in Fig.
3. The absolute values for the gaps in the single-particle spectra for protons
and neutrons are large than 6 MeV in the case of 20Ne. For all other isotopes
considered at least one of those values for the gap are smaller than 4 MeV,
which is an indication that for these nuclei $T=1$ pairing correlations may be
relevant.

Figure 3: (Color online) The energy gap between the last occupied state and
first unoccupied state for various neon isotopes.

The saturation point of nuclear matter is identified by the energy and the
corresponding value for the baryon density. A many-body calculation is
considered successful only, if both quantities are reproduced. In close
analogy we should be satisfied with a microscopic theory describing the bulk
properties of finite nuclei, if the binding energies of the nuclei as well as
the radii of the charge distribution are described.

Corresponding results for the are displayed in Fig. 4. The agreement between
the radii calculated for the proton distribution and the values for the charge
distributionsmarinova  are remarkable. This demonstrates that the contact
interaction of (2) is appropriate to stabilize the density of infinite nuclear
matter as well as the surface of light nuclei at the empirical value.

Fig. 4 also displays the root mean square radii calculated for the density
distributions for the neutrons. The radius increases with the number of
neutrons, representing the development of neuron skin for $N>Z$ with a rather
monotonic slope. A small deviation from this monotonic slope is observed only
at $A$=30, which corresponds to a neutron number $N=20$ a value at which the
$sd$-shell in a simple shell-model would be filled.

Figure 4: (Color online) The neutron and proton root mean square radius of
neon isotopes. The experimental data for the charge radii have been obtained
from measurements of the optical isotope shifts in marinova

The radius, however, is only one observable to be deduced from the density
distribution. In Fig. 5, the neutron and proton density distributions of 24Ne
are shown in terms of contour plots projected to the three cartesian planes.

Figure 5: (Color online) Contour plots of the neutron and proton density
distributions of 24Ne in the xy, xz, and yz plane.

These plots show that the nucleus is elongated in one direction, which is
characterized by a vector pointing to the corner of the cartesian box, in
which the calculation is performed, with positive values for the $y$ and $z$
coordinate and negative value for the $x$ coordinate. The elongation axes for
protons and neutrons are identical, which reflects the attraction between
protons and neutrons.

To visualize more details about the symmetry of these density distributions,
Fig. 6 displays the density distribution in a plane, which contains the axes
of prolate elongation (left panel) and on perpendicular to the elongation
axis. These plots demonstrate that the resulting density distributions are
indeed triaxial and could not be described assuming a cylindrical symmetry, as
has been tried in the work of Ref. ebran:2011 . The same can be observed for
the other isotopes.

Figure 6: (Color online) Contour plots of the neutron density distributions of
24Ne in a plane parallel (left panel) and perpendicular (right panel) to the
elongation axis. The contour lines represent 8 densities ranging in equal
step-sizes between zero and maximal density.

What is the gain in energy, which is related to the deformation of the nuclei?
To answer this question, we have supplemented the calculations in a cartesian
box, which we discussed so far with corresponding calculations in a spherical
cavity (see vandalen:2009b  and discussion around eq.(12) above) assuming
spherical symmetry and, if necessary, a partial occupation of the degenerate
states close to the Fermi energy.

The tri-axial calculations predict stable isotopes up to 32Ne. The HF
calculations allowing only for spherical structures gives stable isotopes up
to 30Ne. Therefore, the differences in the binding energies between tri-axial
calculations and spherical calculations are plotted in Fig. 7 only for
isotopes up to 30Ne. The total energy gain allowing for tri-axial deformation
is largest in the case of 24Ne (0.8 MeV per nucleon, which corresponds to 19.2
MeV in total). The energy gain is reduced for the heavier isotopes and is
close to zero in the case of 30Ne, in which all neutron states of a spherical
$sd$-shell can be filled (see also discussion of the radii of the neutron
distributions above).

Figure 7: (Color online) The energy difference between calculations allowing
for tri-axial structures and calculations allowing only spherical structures.
The binding energy results considered are those without cm correction in both
cases.

In the calculation it is easy to go beyond the description of stable neon
isotopes to the situation of nuclear structures with unbound neutrons. This
should have some relevance for the transition from discrete nuclei to
homogeneous matter which is supposed to occur, e.g. in the crust of neutron
stars. As an example of a system with unbound neutrons we consider a system
with 10 protons and 32 neutrons. Results for the density distributions for the
neutrons are plotted in Fig. 8. In this system, the density profile in z
direction deviates clearly from that in the x and y direction. One has a
quasi-nucleus at the center of the cubic box with a region of a low-density
neutron gas extending in the $z$ direction to the boundary of the box. Since
the box displayed must be considered as one member of a lattice of boxes with
periodic boundary conditions, this region of a low-density neutron gas
continues to the boxes neighbored in $z$-direction, forming a rod of neutrons
with quasi-nuclear structures located along this rod.

This observation is of interest in the sense that it demonstrates that the gas
of evaporated neutrons may show structures like the one discussed here.
However, the present calculation is not realistic in the sense that the
cartesian boxes considered in our calculation are to small to simulate the
structure of quasi-nuclei embedded in a gas of quasi-free neutrons in the
crust of neutron stars. Furthermore, to apply such calculations to the neutron
star crust the contribution of leptons has to be included and the
$\beta$-equilibrium condition has to be imposed.

,

Figure 8: (Color online) Contour plots of the neutron density distributions of
a quasi-nuclear structure with 10 protons and 32 neutrons in the xz plane
(left panel) and yz plane (right panel) of the box considered in the
calculation. The contour lines represent the densities 0.005, 0.02 0.045,
0.07, and 0.09 nucleons fm-3, respectively.

##  IV Summary and conclusion

The aim of this study is to establish a general tool, which allows a Hartree-
Fock (HF) calculation of weakly bound open shell nuclei based on a general
realistic NN interaction. For that purpose the HF equations have been solved
in basis of plane wave states with periodic boundary conditions in a cartesian
box. Stable results are obtained for light nuclei employing a box with a box
size of around 20 fm and basis size of around 1500 states for protons and
neutrons, respectively. This allows the description of structures with general
symmetry. The variational procedure is restricted to single-particle states
with good parity. This stabilizes the position of the center of mass of the
nucleus to the center of the box.

The model used for the NN interaction is a renormalized low momentum
potential, $V_{lowk}$, derived from a realistic NN interaction which fits the
NN scattering data. This $V_{lowk}$ interaction has been supplemented with a
contact interaction with three parameters adjusted to reproduce the empirical
saturation point of symmetric nuclear matter and the symmetry energy in
nuclear matter.

Calculations have been performed for Ne isotopes with even number of neutrons,
ranging from nucleon number $A$=18 to $A$=34. Results for the bulk properties
of these nuclei are very encouraging. The total energy is reproduced with an
accuracy of better than 2 percent and also the radii of the charge
distributions of the stable isotopes are described with high accuracy. The
calculations reproduce the position of the neutron drip-line for the Ne
isotopes.

The deformation degrees of freedom are very important providing an extra
binding energy up to almost 20 MeV as compared to a calculation restricted to
spherical symmetry with partial occupation probabilities for degenerate states
at the Fermi energy. The shapes of the nuclear density distributions are tri-
axial with a prolate elongation axis.

In a next step one may perform Hartree-Fock-Bogoliubov calculations, imposing
the time-reversal symmetry in the HF variation and using the realistic
$V_{lowk}$ for a consistent description of pairing properties. This would also
enable us to investigate the role of pairing in weakly bound nuclei. Also the
study of isotope chains different from Ne may be of interest.

##  V Acknowledgments

This work has been supported by a grant (Mu 705/7-1) of the Deutsche
Forschungsgemeinschaft DFG.

## References

  * (1) E. Caurier, J. Men$\acute{e}$ndez, F. Nowacki, and A. Poves, Phys. Rev. Lett. 100, 052503 (2008).
  * (2) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * (3) T. Niksic, D. Vretenar, and P. Ring, Prog. Part. Nucl. Phys. 66, 519 (2011).
  * (4) R.B. Wiringa, V.G.J. Stoks, and R. Schiavilla, Phys. Rev. C 51, 38 (1995).
  * (5) R. Machleidt, F. Sammarruca, and Y. Song, Phys. Rev. C 53, R1483 (1996).
  * (6) R. Machleidt and D.R. Entem, Phys. Rep. 503, 1 (2011).
  * (7) E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev. Mod. Phys. 81, 1773 (2009).
  * (8) N. Kaiser, S. Fritsch, and W. Weise, Nucl. Phys. A 697, 255 (2002).
  * (9) H. Müther and A. Polls, Prog. Part. Nucl. Phys. 45, 243 (2000).
  * (10) E. N. E. van Dalen, P. Gögelein, and H. Müther, Phys. Rev. C 80, 044312 (2009).
  * (11) S. K. Bogner, T.T.S. Kuo, and L. Coraggio, Nucl. Phys. A684, 432c (2001).
  * (12) S.K. Bogner, A. Schwenk, R.J. Furnstahl, and A. Nogga, Nucl. Phys. A763, 59 (2005).
  * (13) S. K. Bogner, R. J. Furnstahl, S. Ramanan, and A. Schwenk, Nucl. Phys. A784, 79 (2007).
  * (14) P. Bożek, D.J. Dean, and H. Müther, Phys. Rev. C 74, 014303 (2006).
  * (15) E. N. E. van Dalen and H. M¨uther, IJMPE 19, 2077 (2010).
  * (16) J. Kuckei, F. Montani, H. Müther, and A. Sedrakian, Nucl. Phys. A 723, 32 (2003).
  * (17) T.H.R. Skyrme, Nucl. Phys. 9, 615 (1959).
  * (18) P. Bonche and D. Vautherin, Nucl Phys. A 372, 496 (1981).
  * (19) L. Coraggio, A. Covello, A. Gargano, N. Itaco, T. T. S. Kuo, and R. Machleidt, Phys. Rev. C 71, 014307 (2005).
  * (20) L. Coraggio, A. Covello, A. Gargano, N. Itaco, and T. T. S. Kuo, Phys. Rev. C 75, 057303 (2007); Phys. Rev. C 73, 014304 (2006).
  * (21) M. Martini, S. Péru, and M. Dupuis, Phys. Rev. C 83, 034309 (2011).
  * (22) J.-P. Ebran, E. Khan, D. Peña Arteaga, and D. Vretenar, Phys. Rev. C 83, 064323 (2011).
  * (23) K. Suzuki, Prog. Theoret. Phys. 68, 246 (1986).
  * (24) S. Fujii, R. Okamoto, and K. Suzuki, Phys. Rev. C 69, 034328 (2004).
  * (25) P.-G. Reinhard in Computational Nuclear Physics I, Eds.: K. Langanke, J.A. Maruhn, and S.E. Koonin, Springer Verlag, 28 (1991).
  * (26) P. Bonche, H. Flocard, P.-H. Heenen, Comp. Phys. Comm. 171, 49 (2005).
  * (27) W. Ryssens, V. Hellemans, B. Bender, and P.-H. Heenen, arXiv:1405.1897v1 (2014).
  * (28) G. Audi and W. H. Wapstra, Nucl. Phys. A565, 1 (1993).
  * (29) K. Marinova, W. Geithner, M. Kowalska, K. Blaum, S. Kappertz, M. Keim, S. Kloos, G. Kotrotsios, P. Lievens, R. Neugart, H. Simon, and S. Wilbert, Phys. Rev. C 84, 034313 (2011).


# Exclusive production of heavy quarkonia as a probe of the low $x$ and low
scale gluon PDF

a, S. P. Jonesb, A. D. Martinc, M. G. Ryskinc,d and T. Teubnera
a Department of Mathematical Sciences, University of Liverpool, Liverpool L69
3BX, U.K.
b Theoretical Physics Department, CERN, Geneva, Switzerland
c Institute for Particle Physics Phenomenology, Durham University, Durham DH1
3LE, U.K.
d Petersburg Nuclear Physics Institute, NRC Kurchatov Institute, Gatchina, St.
Petersburg, 188300, Russia
E-mail

###### Abstract:

We discuss the exclusive $J/\psi$ photoproduction process, as measured
recently at LHCb, as a means of constraining and ultimately determining the
low $x$ and low $Q$ gluon PDF. The scale dependence of the theoretical
prediction for this process is shown to be systematically improved via a
taming of the known $\overline{\text{MS}}$ result, this amounts to resumming
logarithmically enhanced small-$x$ terms and implementing a small-$Q$ power
correction. The cross section level predictions allow the behaviour of the
gluon PDF in the low $(x,Q)$ domain to be determined.

##  1 Introduction

The recent measurement by LHCb at $\sqrt{s}=\text{13 TeV}$ [1] of central
exclusive production of $J/\psi$ mesons in $pp\rightarrow
p\,+\,\,J/\psi\,+\,p$ events allows constraints on the low $x$ and low $Q$
scale behaviour of the global gluon parton distribution function (PDF).
Together with data from HERA [2] and the LHCb at a lower collision energy of
$\sqrt{s}=\text{7 TeV}$ [1], the kinematic coverage is extended down to about
$x=(M_{J/\psi}/\sqrt{s})e^{-Y}\sim 3\times 10^{-6}$ with $Q\sim\text{1.5
GeV},$ and where $2<Y(J/\psi)<4.5$ is the rapidity interval of the LHCb
forward detector.

We study exclusive $J/\psi$ photoproduction to NLO within collinear
factorisation and show that a systematic taming of the $\overline{\text{MS}}$
approach leads to a refined theoretical result that, until recently, was
plagued by large scale dependence and large logarithms at high energy (small
$x$). As we will see, this is overcome through implementation of a $Q_{0}$ cut
and resummation of a class of double logarithms that, by virtue of the low
scale process, are necessary to ensure a reliable and stable theoretical
prediction.

In Section 2, we briefly recall our model construction and emphasise the
challenges one faces in incorporating the exclusive $J/\psi$ framework into
global PDF analyses. In Section 3, we discuss the new and improved theoretical
result, before mapping this to the cross section level which allows
preliminary predictions about the behaviour of the low $x$ and low scale gluon
PDF to be made.

##  2 Model framework and challenges

We describe the ultraperipheral event $pp\rightarrow p\,\,J/\psi\,\,p$ within
collinear factorisation, with the set up for the quasi-elastic (semi) hard
scattering subprocess $\gamma p\rightarrow J/\psi\,p$ following [3].
Explicitly,

$\text{Im $A$}\sim\sum_{i=q,g}\text{Im}(F_{i}\otimes
C_{i})\otimes\phi_{Q\bar{Q}}^{V},$ (1)

where the $C_{i}$ are the perturbatively calculable subprocess kernels, and
the $F_{i}$ denote generalised parton distribution functions (GPDs). These are
a generalisation of the conventional collinear PDFs used in the global
analyses, required by the off forward kinematics and the presence of skewing
parameter, $\xi$, see e.g. [4]. The projection of the open charm quarks onto
the $J/\psi$ wave function, $\phi_{Q\bar{Q}}^{V}$, is made within LO NRQCD.
Hoodbhoy [5] showed that relativistic corrections to the $c\bar{c}\,\,J/\psi$
transition vertex do not amount to a sizeable correction to the cross section
and would mostly affect the normalisation only which, as we will see, is
adequate given the description of the data in the HERA region (where the
global gluon PDF is known better). Relativistic corrections do not affect the
behaviour as a function of energy, which is completely driven by the $x$
dependence of the gluon PDF.

Our construction of the GPD grids is via the Shuvaev integral transform. It is
a well established means of obtaining GPDs, with the non perturbative
collinear PDFs as input, assuming we are at a sufficiently small $x$ and that
the transform is used in the correct region of parameter space to avoid Regge
pole singularities; $\xi<|X|\ll 1$, see Fig. 1 for our set up. At NLO, the
accuracy is $\mathcal{O}(x)$ so, even for the slightly larger $x$ probed at
HERA, the error is marginal. The imaginary part of the coefficient function
describes the hard matrix element and we include the real part of the
amplitude through the dispersion relation,

$\frac{\text{Re $A$}}{\text{Im
$A$}}=\frac{\pi}{2}\frac{\partial\ln\text{Im}A/W^{2}}{\partial\ln W^{2}},$

where $W$ is the $\gamma p$ subprocess centre of mass energy.

Figure 1: The LO gluon initiated (left panel) and NLO quark initiated (right
panel) $J/\psi$ production diagrams. Here the $F_{i}$ are the GPDs, the dashed
boxed encloses the perturbatively calculable five leg coefficient function and
$\phi^{V}_{Q\bar{Q}}$ describes the formation of the $J/\psi$.

The LHCb collaboration extract the cross section for $\gamma p\rightarrow
J/\psi\,p$ using

$\frac{\text{d}\sigma^{\text{th}}(pp)}{\text{d}Y}=S^{2}(W_{+})N_{+}\sigma^{\text{th}}_{+}(\gamma
p)+S^{2}(W_{-})N_{-}\sigma^{\text{th}}_{-}(\gamma p),$ (2)

with $S^{2}(W_{\pm})$ and $N_{\pm}$ rapidity gap survival factors and photon
fluxes, respectively, for $\gamma p$ centre of mass energies
$W^{2}_{\pm}=M_{J/\psi}\sqrt{s}e^{\pm|Y|}.$ The survival factors, estimated
from general rescattering principles, account for factorisation breaking
corrections and describe the probability that the rapidity gap will not become
populated with additional emissions that would otherwise destroy the
exclusivity of the event. The $W_{-}$ component samples $x\sim 10^{-2}$ and so
$\sigma_{-}$ may be provided by the HERA data while the $W_{+}$ component
samples $x\sim 10^{-5}$ or less and so $\sigma_{+}$ can be extracted from eqn.
(2.2) using $\text{d}\sigma^{\text{th}}(pp)/\text{d}Y$ provided by LHCb.

##  3 Stability of amplitudes at NLO and cross section predictions

The NLO contribution for exclusive $J/\psi$ photoproduction in the
$\overline{\text{MS}}$ collinear factorisation scheme has been known for some
time [3]. However, the result in this scheme has been brought under control,
only recently, by means of eliminating a double counting effect and resumming
a class of large logarithms, see [6] for full details. One may analyse the
high energy asymptotics of the NLO amplitude and observe that the result
contains terms
$\sim\alpha_{s}(\mu_{R})^{2}\ln(\mu_{F}^{2}/m_{c}^{2})\ln(1/\xi)$, which upset
the perturbative convergence through a resulting large scale dependence at
small $\xi$. The renormalisation and factorisation scales are denoted by
$\mu_{R}$ and $\mu_{F}$ respectively. However, by setting $\mu_{F}=m_{c}$,
these terms are completely absorbed providing a resummation for this class of
large logarithmically enhanced terms.

The upshot is a shifting of terms from the NLO coefficient function to the LO
GPD with a remaining scale dependence, $\mu_{F}$, in the NLO coefficient
function and a residual scale dependence, $\mu_{f}$, in the NLO GPD. The
resulting $\gamma p\rightarrow J/\psi\,\,p$ amplitude is therefore recast in
the form

$A(\mu_{f})=C_{g}^{\text{LO}}\otimes
F_{g}(\mu_{F})+\sum_{i=q,g}C^{\text{NLO}}_{\text{rem},i}(\mu_{F})\otimes
F_{i}(\mu_{f}).$ (3)

With this choice, we have absorbed large contributions arising from a specific
scale and momentum fraction hierarchy into the parametrisation of the input
GPD.

We must, in addition, also investigate the effect of a subtraction that is
important at low scales. Recall the DGLAP evolution starts at some low, yet
still perturbative, scale $Q_{0}$. At LO, all contributions below $Q_{0}$ are
included in the input PDFs at $Q_{0}$. At NLO, however, we need to subtract
from the evolution the contribution of $t$-channel loop momentum
$|l|^{2}<Q_{0}^{2}$, otherwise we double count. In this way, we restrict the
virtuality of the four momentum circulating in the gluon ladder diagrams to be
above $Q_{0}$. Note that this has never been a ubiquitous feature of an
$\overline{\text{MS}}$ calculation but is important as the subtraction amounts
to a power correction of $\mathcal{O}(Q_{0}^{2}/\mu_{F}^{2})$, which is
sizeable here because the process sits at a low scale, of the order of
$m_{c}$111Actually, it can be shown that the $Q_{0}$ subtraction effectively
absorbs the quark contribution so that the exclusive $J/\psi$ process is
predominantly a probe of the gluon PDF.. Fig. 2 shows the vast suppression of
scale variation achieved compared to the $\overline{\text{MS}}$ approach.

Figure 2: $\overline{\text{MS}}$ scale variations of $\text{Im }A/W^{2}$ vs.
$W$ at LO and NLO generated using CTEQ6.6 partons at
$\mu_{F}^{2}=\mu_{R}^{2}=Q_{0}^{2},\,m_{c}^{2},\,2m_{c}^{2}$ (left panel). The
corresponding plot with now the scale fixing and $Q_{0}$ cut implemented, with
variations now made w.r.t $\mu_{f}=\mu_{R}$ and
$\mu_{F}^{2}=m_{c}^{2}=2.4\,\text{GeV}^{2}$ fixed (right panel).
$Q_{0}=1.3\,\text{GeV}$ is the starting PDF scale for CTEQ6.6.

Might a fully fledged BFKL resummation improve the situation further? It was
suggested in [7] to resum instead terms like $\sim\alpha_{s}\ln(1/\xi)$
present in the coefficient function to obtain better stability with respect to
factorisation scale variations. However, we do not do so here as the $Q_{0}$
subtraction would need to be implemented in a manner consistent with the use
of the LO BFKL kernel and, furthermore, the LO BFKL gives rise to a too hard
gluon PDF, which is inconsistent with the LHCb data.

Now, with the NLO result sufficiently stable, the natural next step is to
compare the predictions with the data, that is, at the level of the cross
section. In Fig. 3 (left panel), we show predictions using three sets of
global PDF fits [8, 9, 10], evaluated at
$Q_{0}=\mu_{F}=\mu_{f}=\mu_{R}=m_{c}$. The results agree favourably in the
HERA regime, where the global partons are better constrained, whilst in the
LHCb regime we observe huge differences between the various global
PDFs222Exclusive $J/\psi$ photoproduction was also studied within the $k_{t}$
factorisation framework, see [11]. Here, a DLL modified input gluon ansatz
mimicking the resummation of double logarithms through one step of DGLAP
evolution was used and fitted to the experimental data. In this approach,
however, only a subset of the full NLO corrections are incorporated. As the
tamed collinear factorisation result includes all the effects within $k_{t}$
factorisation (to NLO accuracy) but with the full NLO contribution realised,
it is in this sense an improved theoretical prediction. Moreover, the gluon
extracted in this alternative approach is not readily comparable to the
$\overline{\text{MS}}$ partons of the global analyses, without a non-trivial
scheme conversion. Hence, we do not compare at a quantitative level here..
This, together with the right panel in which we propagate the individual
global fit uncertainties to the cross section, shown as a shaded band for
MMHT14 and boundary curves for NNPDF3.0, demonstrates the utility of the
exclusive $J/\psi$ data and supports the claim that their inclusion into the
global PDF analyses would provide a serious constraint for the low $x$ and low
$Q$ behaviour of the gluon PDF.

Figure 3:  Left panel - Cross section predictions using three distinct sets of
global partons [8, 9, 10] with the scales $\mu_{f}^{2}=\mu_{R}^{2}=m_{c}^{2}$
(solid lines). Also shown for CT14 is the prediction with scales
$\mu_{f}^{2}=\mu_{R}^{2}=2m_{c}^{2}$ (dashed line), which demonstrates the
stability of the cross section prediction with respect to scale variations.
Right panel - Cross section predictions using two sets of global partons [8,
10] also displaying the global PDF 68% c. l. uncertainty, which greatly
exceeds the experimental uncertainty. The data are from [2] and the LHCb
$W_{+}$ solutions are constructed from [1], via eqn. (2.2).

## References

  * [1] R. Aaij et al. [LHCb Collaboration], JHEP 1810 (2018) 167; J. Phys. G 41, 055002 (2014).
  * [2] ZEUS Collaboration (S. Chekanov et al.) Eur. Phys. J. C 24 (2002) 345; Nucl. Phys. B 695 (2004) 3; H1 Collaboration (A. Aktas et al.), Eur. Phys. J. C 46 (2006) 585; (C. Alexa et al.) Eur. Phys. J. C 73 (2013) 2466.
  * [3] D. Y. Ivanov, A. Schafer, L. Szymanowski and G. Krasnikov, Eur. Phys. J. C 75 (2015) 75.
  * [4] M. Diehl, Phys. Rept. 388 (2003) 41.
  * [5] P. Hoodbhoy, Phys. Rev. D 56 (1997) 388.
  * [6] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 43 (2016) 035002; Eur. Phys. J. C 76 (2016) 633.
  * [7] D. Y. Ivanov, B. Pire, L. Szymanowski and J. Wagner, [1510.06710 [hep-ph]].
  * [8] R. D. Ball et al. [NNPDF Collaboration], JHEP 1504 (2015) 040.
  * [9] S. Dulat et al., Phys. Rev. D 93 (2016) 033006.
  * [10] L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, Eur. Phys. J. C 75 (2015) 204.
  * [11] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 44, no. 3, 03LT01 (2017).


# A fully hybrid integrated Erbium-based laser

Yang Liu1,2,∗,†, Zheru Qiu1,2,∗, Xinru Ji1,2, Andrea Bancora1,2, Grigory
Lihachev 1,2, Johann Riemensberger1,2, Rui Ning Wang1,2, Andrey Voloshin1,2,
and Tobias J. Kippenberg1,2,‡ 1Institute of Physics, Swiss Federal Institute
of Technology Lausanne (EPFL), CH-1015 Lausanne, Switzerland
2Center for Quantum Science and Engineering,Swiss Federal Institute of
Technology Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Erbium-doped fiber lasers [1, 2, 3] exhibit high coherence and low noise as
required for applications in fiber optic sensing [4], gyroscopes, LiDAR, and
optical frequency metrology [5]. Endowing Erbium-based gain in photonic
integrated circuits can provide a basis for miniaturizing low-noise fiber
lasers to chip-scale form factor, and enable large-volume applications. Yet,
while major progress has been made in the last decade on integrated lasers
based on silicon photonics with III-V gain media [6, 7, 8, 9, 10, 11, 12, 13,
14], the integration of Erbium lasers on chip has been compounded by large
laser linewidth. Recent advances in photonic integrated circuit-based high-
power Erbium-doped amplifiers make a new class of rare-earth-ion-based lasers
possible [15]. Here, we demonstrate a fully integrated chip-scale Erbium laser
that achieves high power, narrow linewidth, frequency agility and the
integration of a III-V pump laser. The laser circuit is based on an Erbium-
implanted ultralow-loss silicon nitride (Si3N4) photonic integrated circuit
[15]. This device achieves single-mode lasing with a free-running intrinsic
linewidth of 50 Hz, a relative intensity noise of $<$-150 dBc/Hz at $>$10 MHz
offset, and an output power up to 17 mW, approaching the performance of fiber
lasers [16] and state-of-the-art semiconductor extended cavity lasers [12, 17,
8, 13]. An intra-cavity microring-based Vernier filter enables wavelength
tunability of $>$ 40 nm within the C- and L-bands while attaining side mode
suppression ratio (SMSR) of $>$ 70 dB, surpassing legacy fiber lasers in
tuning and SMRS performance. This new class of low-noise, tuneable Erbium
waveguide laser could find applications in LiDAR [18], microwave photonics
[19, 20], optical frequency synthesis [21], and free-space communications. Our
approach is extendable to other wavelengths where rare-earth ions can provide
gain, and more broadly, constitutes a novel way to photonic integrated
circuit-based rare-earth-ion-doped lasers.

Erbium-doped fiber lasers (EDFLs) [1, 2, 3] have become indispensable sources
of high coherence laser light for distributed acoustic sensing [22, 4],
optical gyroscopes, free-space optical transmission [23], optical frequency
metrology[5], and high-power laser machining [24, 25] and are considered the
’gold standard’ of laser phase noise. EDFLs exhibit many advantages such as
all-fiberized cavities, alignment-free components, and benefit from the
advantageous Erbium-based gain properties including slow gain dynamics,
temperature insensitivity, low amplification related noise figure, lower
spontaneous emission power coupled to oscillating modes than short
semiconductor gain media [26, 27], and excellent confinement of laser
radiation for high beam quality. These properties along with low phase noise
have led to wide proliferation of Erbium-based fiber lasers in industrial
applications. Erbium ions can provide equally a basis for compact photonic
integrated circuit-based lasers [28] that can benefit from manufacturing at
lower cost, smaller form factor and reduced susceptibility to environmental
vibrations compared to fiber lasers. Prior efforts have been made to implement
chip-based waveguide lasers using Erbium-doped materials such as Al2O3[29,
30], TeO2[31], LiNbO3[32], and Erbium silicate compounds [33] as waveguide
claddings or cores, but the demonstrated laser intrinsic linewidth remained at
the level of MHz [34, 35, 36, 37], far above the sub-100-Hz linewidth achieved
in commercial fiber lasers and state-of-the-art heterogeneously or hybrid
integrated semiconductor-based lasers (Supplementary Note 1). One major
obstacle to realizing narrow-linewidth Erbium waveguide lasers is the
challenge of integrating long and low-loss active waveguides ranging from tens
of centimeters to meters—the lengths routinely deployed in fiber lasers to
ensure low phase noise, single-frequency operation, and sufficient round-trip
gain [16].

Here, we overcome this challenge and demonstrate hybrid integrated Erbium-
doped waveguide lasers (EDWLs) using Si3N4 photonic integrated circuits that
achieve narrow linewidth, frequency agility, high power, and the integration
with pump lasers. Crucial to this advance are meter-scale-long Erbium-
implanted silicon nitride (Er:Si3N4) photonic integrated circuits that can
provide $>$ 30 dB net gain [15] with $>$100 mW output power. The Si3N4
photonic integrated circuit moreover exhibits absence of two-photon absorption
in telecommunication bands[38], radiation hardness for space compatibility,
high power handling of up to tens of watts[39], a lower temperature
sensitivity than silicon, and low Brillouin scattering (a power-limiting
factor in silica-based fiber lasers)[40].

Figure 1:  A hybrid integrated Er:Si3N4 laser. (A) Schematic of a hybrid
integrated Vernier laser consisting of an Erbium-implanted silicon nitride
Er:Si3N4 photonic integrated circuit and an edge-coupled III-V semiconductor
pump laser diode. The intra-cavity microring-based Vernier filter enables
single-mode lasing operation within the Erbium-based gain bandwidth. (B)
Optical image of an Er:Si3N4 laser circuit integrated with micro heaters for
wavelength and phase tuning. The green dashed circle indicates the Erbium-
implanted gain spiral. (C) Optical images of the Erbium-implanted spiral
waveguides and (D) the coupling regime of the Vernier filter indicated by
coloured boxes in (B). (E) Fabrication process flow of the Er:Si3N4 photonic
integrated circuit based on selective Erbium ion implantation.

## Results

### Hybrid integrated Erbium-based Vernier lasers

The laser device is structured as a linear optical cavity with a spiral
Erbium-doped gain waveguide and two reflectors formed by Sagnac loop mirrors
at both ends (Fig.1A). One dichroic loop mirror that consists of a dichroic
directional coupler allows for laser reflection near 1550 nm and optical pump
transmission near 1480 nm, and the other reflector deploys a short waveguide
splitter for broadband reflection. The optical pump can also be injected via a
waveguide taper connected to a microring bus waveguide. The laser device
(Fig.1B) exhibits a compact footprint of only 2 × 3 mm2 with a densely-packed
0.2-m-long Erbium-doped Si3N4 spiral waveguide (Fig.1C) with a cross section
of $0.7\times 2.1\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$ . A narrow-band
intra-cavity Vernier filter designed to achieve sub-GHz 3 dB bandwidth and 5
THz FSR using two cascaded add-drop microring resonators (100 GHz FSRs with 2
GHz difference) (Fig.1D) is deployed to ensure single-mode lasing operation
with a small laser cavity mode spacing of ca. 200 MHz (Supplementary Note 2).
Integrated microheaters are used to align the Vernier filter peak transmission
wavelength to a cavity longitudinal mode. This integrated laser circuit was
fabricated using the photonic Damascene process [41], followed by selective
Erbium ion implantation, post annealing, and heater fabrication (Fig.1E; see
Methods and Supplementary Note 3).

Figure 2:  A hybrid integrated Er:Si3N4 vernier laser operated at single-mode
lasing.  (A) Optical image of a hybrid integrated Er:Si3N4 Vernier laser edge-
coupled with a pump laser diode chip (3SP Technologies, 1943 LCv1). Green
luminescence was observed, stemming from the transition from higher-lying
levels of excited Erbium ions to the ground state. (B) Measured optical
spectrum of single-mode lasing. The inset shows the output power as a function
of the pump power. (C) Measured time-frequency spectrogram of the heterodyne
beatnote between the packaged EDWL and a fully-stabilized frequency comb
(FC1500, Menlo Systems GmbH) over 4 hours. (D) Experimental setup for Vernier
filter characterization (device ID: D85_04_F04_C15_V1). (E) Illustration of
the Vernier effect by measuring the superposed resonances through the
intermediate bus waveguide of the Vernier filter. (F) Zoomed-in range of the
measured transmission. Colored circles indicate the resonances of each
microring. (G) Wide-range transmission overlaid with the Erbium ion gain
spectrum. (H) Frequency spacing variation between adjacent ring resonances,
yielding a Vernier spacing of 4.65 THz, corresponding to 37.1 nm. (I) The
curve fitting of the measured through port transmission of the resonance
indicated in (F), and the calculated filtering response at the drop port.

To demonstrate a fully integrated EDWL, we performed photonic packaging via
hybrid integration in a custom 14-pin butterfly package. The 1480 nm InP
Fabry-Pérot (FP) laser diode (LD) was edge coupled to one of the laser
cavities on an Er:Si3N4 photonic integrated circuit (Fig.2A), with simulated
coupling loss of < 3 dB. The laser output waveguide was end-coupled and glued
with a cleaved UHNA-7 optical fiber spliced to a SMF-28 optical fiber pigtail,
exhibiting 2.7 dB coupling loss at 1550 nm. The pump LD, a Peltier element, a
thermistor, and all microheaters are connected to butterfly pins using wire
bonding. The integrated micro-heaters were used for the temperature control of
the Vernier filter and the phase-shifter section to configure single-mode
lasing and wavelength tuning. The Erbium ions can be optically excited by the
pump light emitted from the multi-longitudinal-mode pump LD ($>$4 nm spectral
linewidth near 1480 nm), providing 1.9 dB/cm of measured net gain coefficient
[15]. The optical spectrum of the collected laser output shows a single-mode
lasing operation with $>$ 70 dB of side mode suppression ratio (SMSR) at 0.1
nm resolution bandwidth (Fig.2B). This high 72-dB SMSR was made possible using
the drop port of the the narrow passband intra-cavity Vernier filter, which
can select the lasing mode and reject the broadband amplified spontaneous
emission noise. This record high SMSR surpasses what has been reported in
integrated Erbium lasers, fiber lasers, and integrated semiconductor-based
lasers (Supplementary Note 1), typically below 60 dB that is usually limited
by intra-cavity filtering performance. Conversely, this is challenging to
implement in legacy fiber-based Erbium lasers where the filtering components
based on long Bragg gratings can only offer several GHz wide passband with
grating side lobes and lack of broadband wavelength tuning capability. We
observed an off-chip lasing threshold pump power of ca. 20 mW and an on-chip
slope efficiency of 6.7 $\%$ when sweeping the pump power (Fig.2B inset),
which can be further optimized by reducing the coupling loss and the cavity
loss. The fully packaged laser showed a frequency drift of $<$ 20 MHz over 4
hours (Fig.2C) when performing a heterodyne beatnote measurement with a fully-
stabilized optical frequency comb indicating a good frequency stability due
the monolithic nature of the laser comprised of both cavity and gain medium.
During a 24-hour test, this laser showed a frequency drift of $<$ 140 MHz
without mode hops (Supplementary Note 4), representing a comparable long-term
frequency stability as a commercial diode laser (Toptica CTL).

### Single-mode lasing and wavelength tuning

The use of photonic integrated circuits and Vernier structures (Fig.2) enables
to endow the integrated Erbium laser with broad wavelength tuning, a
capability that bulk fiber lasers lack. We investigated the intra-cavity
filtering properties by characterizing the optical transmission of the middle
bus waveguide (Fig.2D). The measured transmission of the individual resonators
used for the Vernier filter is shown in Fig.2F and the designed 2 GHz FSR was
experimentally attained (98 GHz and 100 GHz, respectively), leading to a
measured Vernier filter FSR of 4.65 THz that corresponds to 37.1 nm span near
1550 nm wavelength (Fig.2G,H). Such a large Vernier FSR ensures the single-
wavelength lasing within the Erbium emission wavelength range (Fig.2G). By
overlapping the resonances from the two resonators, i.e. vanishing the
frequency spacing (Fig.2H), the lasing wavelength is determined. By fitting
the resonance linewidth near 194.8 THz (Fig.2I), we obtain an external
coupling rate $\kappa_{\mathrm{ex,0}}/2\pi=$ 411 MHz (between the microring
and the bus waveguide) and an intrinsic loss rate $\kappa_{\mathrm{0}}/2\pi=$
42.5 MHz. This strong over-coupled configuration
($\frac{\kappa_{\mathrm{ex}}}{\kappa_{\mathrm{0}}}>10$) can ensure that the
Vernier filter _simultaneously_ achieves a narrow 3-dB passband bandwidth of
636 MHz and in principle a low insertion loss. Such strong overcoupling can
allow for low loss operation of the Vernier filter, which however in the
current device was not attained. The Vernier filter exhibits an insertion loss
of -3.2 dB due to the parasitic loss induced by the coupling from the
fundamental waveguide mode to higher order modes, which leads to a suboptimal
coupling ideality [42] of $I=0.87$ (Supplementary Note 5).

Figure 3:  Demonstration of wideband tuning of the laser wavelength.  (A)
Experimental setup for laser wavelength tuning demonstration (device ID:
D85_04_F8_C508_VL2.0) (B) Operating principle of the wavelength tuning of the
vernier laser. The traces indicate the transmission of each over-coupled
microresonator and the entire vernier filter, respectively. (C) Two
dimensional laser wavelength tuning map, showing the wavelength of the
predominant lasing mode as a function of the electrical power applied to the
two micro heaters. The dashed schematically indicates the approach to coarse
wavelength tuning. The white regions indicate that the expected laser emission
was missing due to the microring resonance misalignment or a competing lasing
mode when approaching the edge of the WDM filter transmission band. (D)
Measured optical spectra of single-mode lasing tuned over 40 nm wavelength
range. The optical spectrum analyzer’s resolution bandwidth is set to 0.1 nm.

Next, we demonstrate the wavelength tunability (Fig. 3A). The coarse tuning of
laser wavelength was carried out by switching the aligned resonance of two
microresonators (Fig. 3B). The step size of ca. 0.8 nm was determined by the
microring FSR. Fine tuning of the wavelength can be achieved by simultaneously
shifting the two resonators in the same direction and adjusting the phase
shifter to align the corresponding cavity longitudinal mode (164 MHz spacing)
to the Vernier filter passband. Figure 3C shows the 2-dimensional (2D) laser
wavelength tuning map when varying the electrical power applied to the
microheaters. From the recorded entire 2D map of wavelengths we selected the
settings marked in Fig.3C. This allowed for continuous and deterministic
tuning over the entire wavelength band from 1548.1 nm to 1585.8 nm,
maintaining power of $>4\leavevmode\nobreak\ \mathrm{mW}$ and SMSR of
$>70\leavevmode\nobreak\ \mathrm{dB}$ (Fig.3D). Such wavelength tunability
cannot be achieved in conventional rare-earth-ion-doped fiber lasers without
the use of free space etalon filters. The wavelength tuning range was limited
by the Vernier filter FSR and the wavelength-division multiplexing coupler
transmission band (Supplementary Note 6). During heater power scanning, we
note that a few of wavelength tuning steps were missed due to the misalignment
of microring resonances of the Vernier filter . During tuning, the phase
shifter was adjusted to maximize the output power at the desired mode. A
maximum fiber-coupled output power of ca. $17$ mW were measured at 1585 nm
with 219 mW pump power. Other competing lasing modes apart from the
predominant lasing mode were observed when using high pump power, due to the
fact that the large Si3N4 waveguide cross section allows for multiple
transversal optical modes that can coincidently satisfy the lasing condition
(Supplementary Note 8).

Figure 4:  Laser noise properties and the fully hybrid integration of an EDWL.
(A) Experimental setups for the measurement of laser frequency noise, relative
intensity noise, and intrinsic laser linewidth. (B) Optical image of a fully
hybrid integrated EDWL assembly. (C) Measured laser frequency noise based on
heterodyne detection with reference lasers. EDWL1 device ID:
D85_04_F8_C508_VL2.0; EDWL2 device ID: D85_04_F8_C508_VL2.1; EDWL3 device ID:
D85_04_F8_C508_VL2.1. (D) Measured and fitted spectra of delayed self-
heterodyne interferometric measurement for intrinsic laser linewidth
investigation. (E) Measured laser relative intensity noise (RIN) based on
direct photodetetion. (F) Relaxation oscillation peaks under varied pump
power.

### Frequency and intensity noise measurements

To demonstrate the low noise features of the free-running EDWLs, we
characterized the frequency noise, the intrinsic laser linewidth, and the
relative intensity noise (RIN), respectively (Fig. 4A). Firstly, a reference
external cavity diode laser (free running Toptica CTL) was tuned close to the
lasing wavelength near 1560 nm of an EDWL (not packaged) with ca. 3 mW output
power for heterodyne photodetection. The in-phase and quadrature components of
the sampled beatnote time trace was processed using Welch’s method [43] to
obtain the single-side power spectral density (PSD) of frequency noise
$S_{\delta v}(f)$. The frequency noise PSD (red line) reached a plateau of
$h_{0}=62.0\leavevmode\nobreak\ \mathrm{Hz^{2}/Hz}$ at the offset frequency of
6 MHz, corresponding to a Lorentzian linewidth of $\pi
h_{0}=194.8\leavevmode\nobreak\ \mathrm{Hz}$; this measured white noise floor
was masked by the ECDL’s white noise floor (Fig. 4C). We also applied the
delayed self-heterodyne interferometric measurement [44] to validate the
intrinsic linewidth (Fig. 4D), which generates a power spectrum of the
autocorrelation of the laser line under sub-coherence condition (Supplementary
Note 9). In the offset frequency range from 10 kHz to 2.5 MHz where a
relaxation oscillation peak was observed, the Erbium laser shows a higher
frequency noise due to the laser cavity fluctuation caused by the pump laser
noise transduction and the thermorefractive noise in the microresonator [45].
The measured frequency noise at offset frequencies of $<$10 kHz was dominated
by ECDL characteristic noise features [46]. We achieved a record low intrinsic
linewidth (blue line) of $\pi h_{0}=50.1\leavevmode\nobreak\ \mathrm{Hz}$
($h_{0}=15.9\leavevmode\nobreak\ \mathrm{Hz^{2}/Hz}$ ) in an Erbium waveguide
laser with a higher output power of 10 mW, when beating against a low-noise
Erbium fiber laser (Koheras Adjustik). The fully packaged EDWL (Fig. 4B) with
2.8 mW output power shows a comparable intrinsic linewidth (purple line) and a
lower frequency noise at the mid-range offset frequencies. Using laser cavity
designs with reduced cold cavity losses and increased mode area, hertz-
linewidth EDWL can be feasibly achieved (Supplementary Note 10).

The full width at half maximum (FWHM) of the integral linewidth associated
with Gaussian contribution was obtained by integrating the frequency noise PSD
from the inverse of measurement time (1/$T_{0}$) up to the frequency where
$S_{\delta v}(f)$ intersects with the $\beta$-separation line $S_{\delta
v}(f)=8\ln(2)f/\pi^{2}$ (dashed line) [47]. With the integrated surface $A$,
we obtained a minimum FWHM linewidth ($8\ln(2)A^{1/2}$) of the free-running
EDWL is 82.2 kHz at 1 ms measurement time, which does not yet supersede a
fiber laser, but is lower than 166.6 kHz of an ECDL (Toptica CTL)
characterized as a reference laser for comparison. For comparison, the
commercial stabilized fiber-based laser shows 2.4 kHz of the FWHM linewidth at
1 ms measurement time.

Next, we show that the Erbium waveguide laser features a lower RIN compared to
a commercial fiber laser (Koheras Adjustik, KOH45) (Supplementary Note 11).
The waveguide laser shows a RIN down to $-130\leavevmode\nobreak\
\mathrm{dBc/Hz}$ (yellow& purple) at mid-range offset frequencies between 10
kHz and 1 MHz, lower than the fiber laser RIN (grey) that has a PSD pole
induced by relaxation oscillation (Fig. 4E). The mid-range RIN was mainly
limited by the pump laser RIN transduction which even contributed to an
increased RIN by 5 dB for the un-packaged EDWL. The pump RIN noise
transduction at frequency above $20\leavevmode\nobreak\ \mathrm{MHz}$ was
suppressed due to the slow dynamics of Erbium ions. The relaxation oscillation
frequency can be calculated by
$f_{\mathrm{r}}\approx\frac{1}{2\pi}\sqrt{\frac{P_{\mathrm{cav}}\kappa}{P_{\mathrm{sat}}\tau}}$
where $P_{\mathrm{sat}}$ is saturation power of the gain medium,
$P_{\mathrm{cav}}$ is the laser cavity power,$\kappa$ is the cold cavity loss
rate, and $\tau$ is the Erbium ion upper-state lifetime (Supplementary Note
12). The waveguide laser RIN reduced to $<-155\leavevmode\nobreak\
\mathrm{dBc/Hz}$ at offset frequencies of $>10$ MHz. We observed that the
relaxation oscillation frequency the waveguide laser varied from 0.3 MHz to
2.4 MHz when increasing the optical pump power (Fig. 4F), which is higher than
the case in the fiber laser (typically $<100\leavevmode\nobreak\ \mathrm{kHz}$
). This higher relaxation oscillation frequency originates from the smaller
saturation power and the shorter Erbium upper-state lifetime of 3.4 ms [15].

### Summary

In summary, we have demonstrated a photonic integrated circuit-based Erbium
laser that achievesd sub-100 Hz intrinsic linewidth, low RIN noise, $>$ 72 dB
SMSR, and 40 nm wide wavelength tunability with power exceeding 10 mW. The
Erbium-doped waveguide lasers use foundry compatible silicon nitride
waveguides, and have the potential to combine fiber-laser coherence with low
size, weight, power and cost of integrated photonics. Such a laser may find
application in existing applications such as coherent sensing, and may equally
provide a disruptive solution for emerging applications that require high
volumes, such as lasers for coherent FMCW LiDAR, or for coherent optical
communications where iTLA (integrated tunable laser assembly) have been widely
deployed, but fiber lasers’ high coherence is increasingly demanded for
advanced high-speed modulation formats while their use has been impeded by the
high cost and large size. Co-doping other rare-earth ions such as ytterbium
(emission at $1.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$) and thulium
($0.8\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$1.45\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$2.0\text{\,}\mathrm{\SIUnitSymbolMicro m}$) will moreover allow access to
other wavelengths. Looking to the future, the compatibility of silicon nitride
with heterogeneously integrated thin-film lithium niobate[48], as well as
piezoelectric thin films [49, 11], and Erbium waveguide amplifiers [15]
provides the capability to create fully-integrated high-speed, low-noise,
high-power optical engines for LiDAR, long-haul optical coherent
communications, and analog optical links.

Funding Information: This work was supported by the Air Force Office of
Scientific Research (AFOSR) under Award No. FA9550-19-1-0250, and by contract
W911NF2120248 (NINJA) from the Defense Advanced Research Projects Agency
(DARPA), Microsystems Technology Office (MTO). This work further supported by
the EU H2020 research and innovation programme under grant No. 965124
(FEMTOCHIP), by the SNSF under grant no. 201923 (Ambizione), and by the Marie
Sklodowska-Curie IF grant No. 898594 (CompADC) and grant No. 101033663
(RaMSoM). Acknowledgments: Silicon nitride samples were fabricated in the EPFL
Center of MicroNanoTechnology (CMi).

Author contributions: Y.L. and Z.Q. performed the experiments. Y.L. carried
out data analysis and simulations. Y.L. and Z.Q. designed Si3N4 waveguide
laser chips. X.J., G. L. and J.R. provided experimental supports. A. B. and A.
V. designed and performed the device packaging. R.N.W., Z.Q. and X.J.
fabricated the passive Si3N4 samples. Y.L. wrote the manuscript with the
assistance from Z.Q. and the input from all co-authors. T.J.K supervised the
project.

Data Availability Statement: The code and data used to produce the plots
within this work will be released on the repository Zenodo upon publication of
this preprint.

Competing interests T.J.K. is a cofounder and shareholder of LiGenTec SA, a
start-up company offering Si3N4 photonic integrated circuits as a foundry
service.

## References

  * Barnes  _et al._ [1989] W. L. Barnes, P. R. Morkel, L. Reekie, and D. N. Payne, Optics Letters 14, 1002 (1989).
  * Suzuki  _et al._ [1989] K. Suzuki, Y. Kimura, and M. Nakazawa, Japanese Journal of Applied Physics 28, L1000 (1989).
  * Chow  _et al._ [1996] J. Chow, G. Town, B. Eggleton, M. Ibsen, K. Sugden, and I. Bennion, IEEE Photonics Technology Letters 8, 60 (1996).
  * Zadok  _et al._ [2012] A. Zadok, Y. Antman, N. Primerov, A. Denisov, J. Sancho, and L. Thevenaz, Laser & Photonics Reviews 6, L1 (2012).
  * Xu and Wise [2013] C. Xu and F. W. Wise, Nature Photonics 7, 875 (2013).
  * Tran  _et al._ [2019] M. A. Tran, D. Huang, J. Guo, T. Komljenovic, P. A. Morton, and J. E. Bowers, IEEE Journal of Selected Topics in Quantum Electronics 26, 1 (2019).
  * de Beeck  _et al._ [2021] C. O. de Beeck, F. M. Mayor, S. Cuyvers, S. Poelman, J. F. Herrmann, O. Atalar, T. P. McKenna, B. Haq, W. Jiang, J. D. Witmer, _et al._ , Optica 8, 1288 (2021).
  * Fan  _et al._ [2020] Y. Fan, A. van Rees, P. J. Van der Slot, J. Mak, R. M. Oldenbeuving, M. Hoekman, D. Geskus, C. G. Roeloffzen, and K.-J. Boller, Optics express 28, 21713 (2020).
  * Morton  _et al._ [2022] P. A. Morton, C. Xiang, J. B. Khurgin, C. D. Morton, M. Tran, J. Peters, J. Guo, M. J. Morton, and J. E. Bowers, Journal of Lightwave Technology 40, 1802 (2022).
  * Guo  _et al._ [2022] Y. Guo, X. Li, M. Jin, L. Lu, J. Xie, J. Chen, and L. Zhou, APL Photonics (2022).
  * Lihachev  _et al._ [2022] G. Lihachev, J. Riemensberger, W. Weng, J. Liu, H. Tian, A. Siddharth, V. Snigirev, V. Shadymov, A. Voloshin, R. N. Wang, J. He, S. A. Bhave, and T. J. Kippenberg, Nature Communications 13, 3522 (2022).
  * Xiang  _et al._ [2021] C. Xiang, J. Guo, W. Jin, L. Wu, J. Peters, W. Xie, L. Chang, B. Shen, H. Wang, Q.-F. Yang, D. Kinghorn, M. Paniccia, K. J. Vahala, P. A. Morton, and J. E. Bowers, Nature Communications 12, 6650 (2021).
  * Li  _et al._ [2021a] B. Li, W. Jin, L. Wu, L. Chang, H. Wang, B. Shen, Z. Yuan, A. Feshali, M. Paniccia, K. J. Vahala, and J. E. Bowers, Optics Letters 46, 5201 (2021a).
  * Xiang  _et al._ [2023] C. Xiang, W. Jin, O. Terra, B. Dong, H. Wang, L. Wu, J. Guo, T. J. Morin, E. Hughes, J. Peters, Q.-X. Ji, A. Feshali, M. Paniccia, K. J. Vahala, and J. E. Bowers, Three-dimensional integration enables ultra-low-noise, isolator-free Si photonics (2023), arXiv:2301.09989 [physics].
  * Liu  _et al._ [2022] Y. Liu, Z. Qiu, X. Ji, A. Lukashchuk, J. He, J. Riemensberger, M. Hafermann, R. N. Wang, J. Liu, C. Ronning, and T. J. Kippenberg, Science 376, 1309 (2022).
  * Fu  _et al._ [2017] S. Fu, W. Shi, Y. Feng, L. Zhang, Z. Yang, S. Xu, X. Zhu, R. A. Norwood, and N. Peyghambarian, Journal of the Optical Society of America B 34, A49 (2017).
  * Xiang  _et al._ [2020] C. Xiang, W. Jin, J. Guo, J. D. Peters, M. J. Kennedy, J. Selvidge, P. A. Morton, and J. E. Bowers, Optica 7, 20 (2020).
  * Kim  _et al._ [2021] I. Kim, R. J. Martins, J. Jang, T. Badloe, S. Khadir, H.-Y. Jung, H. Kim, J. Kim, P. Genevet, and J. Rho, Nature Nanotechnology 16, 508 (2021).
  * Liu  _et al._ [2020a] Y. Liu, A. Choudhary, D. Marpaung, and B. J. Eggleton, Advances in Optics and Photonics 12, 485 (2020a).
  * Marpaung  _et al._ [2019] D. Marpaung, J. Yao, and J. Capmany, Nature Photonics 13, 80 (2019).
  * Spencer  _et al._ [2018] D. T. Spencer, T. Drake, T. C. Briles, J. Stone, L. C. Sinclair, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, A. Bluestone, N. Volet, T. Komljenovic, L. Chang, S. H. Lee, D. Y. Oh, M.-G. Suh, K. Y. Yang, M. H. P. Pfeiffer, T. J. Kippenberg, E. Norberg, L. Theogarajan, K. Vahala, N. R. Newbury, K. Srinivasan, J. E. Bowers, S. A. Diddams, and S. B. Papp, Nature 557, 81 (2018).
  * Soriano-Amat  _et al._ [2021] M. Soriano-Amat, H. F. Martins, V. Durán, L. Costa, S. Martin-Lopez, M. Gonzalez-Herraez, and M. R. Fernández-Ruiz, Light: Science & Applications 10, 51 (2021).
  * Dix-Matthews  _et al._ [2021] B. P. Dix-Matthews, S. W. Schediwy, D. R. Gozzard, E. Savalle, F.-X. Esnault, T. Lévèque, C. Gravestock, D. D’Mello, S. Karpathakis, M. Tobar, and P. Wolf, Nature Communications 12, 515 (2021).
  * Nilsson [2010] J. Nilsson, Optics InfoBase Conference Papers 20, 10.1364/cleo.2010.ctuc1 (2010), iSBN: 9781557528896 Publisher: IEEE.
  * Jackson [2012] S. D. Jackson, Nature Photonics 6, 423 (2012).
  * Siegman [1986] A. E. Siegman,  _Lasers_ (Univ. Science books, Mill Valley, Calif, 1986).
  * Streifer  _et al._ [1982] W. Streifer, D. Scifres, and R. Burnham, IEEE Journal of Quantum Electronics 18, 1918 (1982), conference Name: IEEE Journal of Quantum Electronics.
  * Bradley and Pollnau [2011] J. Bradley and M. Pollnau, Laser & Photonics Reviews 5, 368 (2011), zSCC: 0000274.
  * Rönn  _et al._ [2019] J. Rönn, W. Zhang, A. Autere, X. Leroux, L. Pakarinen, C. Alonso-Ramos, A. Säynätjoki, H. Lipsanen, L. Vivien, E. Cassan, and Z. Sun, Nature Communications 10, 432 (2019).
  * Mu  _et al._ [2020] J. Mu, M. Dijkstra, J. Korterik, H. Offerhaus, and S. M. García-Blanco, Photonics Research 8, 1634 (2020).
  * Frankis  _et al._ [2020] H. C. Frankis, H. M. Mbonde, D. B. Bonneville, C. Zhang, R. Mateman, A. Leinse, and J. D. B. Bradley, Photonics Research 8, 127 (2020).
  * Cai  _et al._ [2022] M. Cai, K. Wu, J. Xiang, Z. Xiao, T. Li, C. Li, and J. Chen, IEEE Journal of Selected Topics in Quantum Electronics 28, 1 (2022).
  * Wang  _et al._ [2018] X. Wang, P. Zhou, Y. He, and Z. Zhou, Optical Materials Express 8, 2970 (2018).
  * Li  _et al._ [2018] N. Li, E. S. Magden, Z. Su, N. Singh, A. Ruocco, M. Xin, M. Byrd, P. T. Callahan, J. D. Bradley, D. Vermeulen, and M. R. Watts, 2018 Conference on Lasers and Electro-Optics, CLEO 2018 - Proceedings 26, 4560 (2018), iSBN: 9781943580422\.
  * Belt  _et al._ [2013] M. Belt, T. Huffman, M. L. Davenport, W. Li, J. S. Barton, and D. J. Blumenthal, Optics Letters 38, 4825 (2013).
  * Purnawirman  _et al._ [2017] Purnawirman, N. Li, E. S. Magden, G. Singh, N. Singh, A. Baldycheva, E. S. Hosseini, J. Sun, M. Moresco, T. N. Adam, G. Leake, D. Coolbaugh, J. D. B. Bradley, and M. R. Watts, Optics Express 25, 13705 (2017), zSCC: NoCitationData[s0].
  * Li  _et al._ [2021b] T. Li, K. Wu, M. Cai, Z. Xiao, H. Zhang, C. Li, J. Xiang, Y. Huang, and J. Chen, APL Photonics 6, 101301 (2021b).
  * Liu  _et al._ [2021] J. Liu, G. Huang, R. N. Wang, J. He, A. S. Raja, T. Liu, N. J. Engelsen, and T. J. Kippenberg, Nature Communications 12, 2236 (2021).
  * Brasch  _et al._ [2016] V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, Science 351, 357 (2016).
  * Gyger  _et al._ [2020] F. Gyger, J. Liu, F. Yang, J. He, A. S. Raja, R. N. Wang, S. A. Bhave, T. J. Kippenberg, and L. Thévenaz, Physical Review Letters 124, 013902 (2020).
  * Pfeiffer  _et al._ [2016] M. H. P. Pfeiffer, A. Kordts, V. Brasch, M. Zervas, M. Geiselmann, J. D. Jost, and T. J. Kippenberg, Optica 3, 20 (2016), publisher: Optica Publishing Group.
  * Pfeiffer  _et al._ [2017] M. H. Pfeiffer, J. Liu, M. Geiselmann, and T. J. Kippenberg, Physical Review Applied 7, 024026 (2017).
  * Welch [1967] P. Welch, IEEE Transactions on Audio and Electroacoustics 15, 70 (1967).
  * Richter  _et al._ [1986] L. Richter, H. Mandelberg, M. Kruger, and P. McGrath, IEEE Journal of Quantum Electronics 22, 2070 (1986).
  * Huang  _et al._ [2019] G. Huang, E. Lucas, J. Liu, A. S. Raja, G. Lihachev, M. L. Gorodetsky, N. J. Engelsen, and T. J. Kippenberg, Physical Review A 99, 061801 (2019).
  * Liu  _et al._ [2020b] J. Liu, E. Lucas, A. S. Raja, J. He, J. Riemensberger, R. N. Wang, M. Karpov, H. Guo, R. Bouchand, and T. J. Kippenberg, Nature Photonics 14, 486 (2020b).
  * Di Domenico  _et al._ [2010] G. Di Domenico, S. Schilt, and P. Thomann, Applied Optics 49, 4801 (2010).
  * Churaev  _et al._ [2021] M. Churaev, A. Riedhauser, R. N. Wang, C. Möhl, T. Blésin, M. A. Anderson, V. Snigirev, A. Siddharth, Y. Popoff, D. Caimi, S. Hönl, J. Riemensberger, J. Liu, P. Seidler, and T. J. Kippenberg,  , 1 (2021), _eprint: 2112.02018.
  * Liu  _et al._ [2020c] J. Liu, H. Tian, E. Lucas, A. S. Raja, G. Lihachev, R. N. Wang, J. He, T. Liu, M. H. Anderson, W. Weng, S. A. Bhave, and T. J. Kippenberg, Nature 583, 385 (2020c).
  * Polman  _et al._ [1991] A. Polman, D. C. Jacobson, D. J. Eaglesham, R. C. Kistler, and J. M. Poate, Journal of Applied Physics 70, 3778 (1991).

## Methods

### Device fabrication and ion implantation

This ultralow loss Si3N4 photonic integrated laser circuit is fabricated using
the photonic Damascene process [41]. We applied selective Erbium ion
implantation [50] (a total fluence of $1\times 10^{16}\leavevmode\nobreak\
\text{ions}\leavevmode\nobreak\ \text{cm}^{-2}$ at a maximum beam energy of 2
MeV) to the pre-fabricated passive Si3N4 photonic integrated circuits to endow
the spiral waveguide with Erbium-based optical gain, while we kept other
passive components un-doped by selectively masking a portion of the chip with
photoresist (Fig.1E and Supplementary Note 1). We achieved a high doping
concentration of $3.25\times 10^{20}\leavevmode\nobreak\
\text{ions}\leavevmode\nobreak\ \text{cm}^{-3}$, more than one order of
magnitude higher than that of conventional Erbium-doped fibers. This allows
for high roundtrip net gain of 1.9 dB/cm (characterized from a 4.5-mm-long
Erbium-doped waveguide) [15]. After ion implantation, we annealed the sample
at $1000\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for one hour to activate the
Erbium ions and heal implantation defects. The optical gain is provided by the
stimulated emission of erbium ions excited by an optical pump (Fig.1A inset).
Micro-heaters were subsequently added atop the silica upper cladding after the
ion implantation and post annealing processes.