# Enhanced laser-driven proton acceleration via improved fast electron heating
in a controlled pre-plasma

Leonida A. Gizzi  Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy
INFN, Sez. Pisa, Italy  correspondence to<EMAIL_ADDRESS>Elisabetta Boella  Physics Department, Lancaster University, Bailrigg,
Lancaster LA1 4YW, UK  Cockcroft Institute, Sci-Tech Daresbury, Keckwick Lane,
Warrington WA4 4AD, UK<EMAIL_ADDRESS>Luca Labate  Intense Laser
Irradiation Laboratory, INO-CNR, Pisa, Italy  INFN, Sez. Pisa, Italy
<EMAIL_ADDRESS>Federica Baffigi  Intense Laser Irradiation
Laboratory, INO-CNR, Pisa, Italy  Pablo J. Bilbao  Physics Department,
Lancaster University, Bailrigg, Lancaster LA1 4YW, UK  Fernando Brandi
Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy  Gabriele
Cristoforetti  Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy
Alberto Fazzi  Dipartimento di Energia, Politecnico di Milano, Italy  INFN,
Sezione di Milano, Italy  Lorenzo Fulgentini  Intense Laser Irradiation
Laboratory, INO-CNR, Pisa, Italy  Dario Giove  INFN-LASA, Segrate, Italy
Petra Koester  Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy
Daniele Palla  Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy
Paolo Tomassini  Intense Laser Irradiation Laboratory, INO-CNR, Pisa, Italy

###### Abstract

The interaction of ultraintense laser pulses with solids is largely affected
by the plasma gradient at the vacuum-solid interface, which modifies the
absorption and ultimately, controls the energy distribution function of heated
electrons. A micrometer scale-length plasma has been predicted to yield a
significant enhancement of the energy and weight of the fast electron
population and to play a major role in laser-driven proton acceleration with
thin foils. We report on recent experimental results on proton acceleration
from laser interaction with foil targets at ultra-relativistic intensities. We
show a three-fold increase of the proton cut-off energy when a micrometer
scale-length pre-plasma is introduced by irradiation with a low energy
femtosecond pre-pulse. Our realistic numerical simulations agree with the
observed gain of the proton cut-off energy and confirm the role of stochastic
heating of fast electrons in the enhancement of the accelerating sheath field.

## Introduction

Ion acceleration driven by ultraintense lasers using Target Normal Sheath
Acceleration (TNSA) [1] is establishing itself as a powerful technique to
access relatively high energy ion beams in a compact and affordable layout.
Since the original investigations, effort has been dedicated [2, 3] to enhance
the cut-off energy and the flux of the accelerated ions via a number of
advances of laser and target specifications using plasma mirror for ultra-high
laser contrast with ultra-thin targets [4], cryogenic [5] or nano-structured
targets [6]. To date, however, practical exploitation of laser-driven ion
acceleration relies heavily on the original TNSA configuration based on thin
foil targets and optimized contrast without plasma mirror, possibly operating
at the repetition rate required for applications like radiobiology [7, 8],
where high dose irradiation is needed for meaningful studies. In this context,
great attention is being dedicated to the control of accelerated ions,
including energy cut-off, beam divergence, charge and emittance. Target
optimisation and engineering, looking at different properties of surface,
geometry and conductivity, are becoming crucial in this effort to ensure
reliable and reproducible performance of beamlines based on laser-driven
acceleration.

From a more fundamental viewpoint, recent investigations to further improve
ion acceleration performance focus on the interaction parameters affecting the
energy density of laser-heated fast electrons, namely their divergence, peak
current and temperature. It is well known that fast electron generation is
strongly modified by the size of the laser focal spot [9, 10], and more
effectively by the scale-length of the plasma at the vacuum-target interface.
Indeed, several modelling studies based on Particle-In-Cell (PIC) simulations
predict a significant increase of the proton cut-off energy with optimized
pre-plasma conditions. In _Sentoku et al._[11], simulations were used to
investigate fast electron heating in the standing wave occurring in the pre-
plasma from superposition of incident and reflected wave, showing a linear
increase of the ion energy cut-off with the scale-length of the pre-plasma.
Similar results were shown in a more recent modelling study [12] where the
role of temporal pulse shape was also explored. The role of standing waves was
then further investigated in a colliding pulses configuration with two pulses
incident on the target from opposite incident angles [13], showing enhancement
in the fast electron temperature and ion energy cut-off for both ultra-high
contrast laser systems using double plasma mirror. In the same study, a strong
linear increase of the cut-off energy with the scale-length of the pre-plasma
is shown for both colliding pulse and single pulse at constant energy, with
the highest proton cut-off energy occurring for a micrometer scale-length
plasma. More in details, in _Sentoku et al._ [11], a moderate increase of the
proton energy up to 1.5 times was found at mildly relativistic intensities of
$\mathrm{10^{18}}\>\mathrm{W/cm^{2}}$ at 1.053 $\mathrm{\mu m}$ wavelength for
a micrometer-scale pre-plasma. The increase was explained as due to the
stochastic heating of electrons in the standing wave structure responsible for
an enhanced ponderomotive acceleration of electrons. In the same study it is
shown that this enhancement of proton energy can be expected to grow to a
factor of 4 or more for longer scale-length and eventually saturate above 6-8
$\mathrm{\mu m}$ scale-length. On the other hand, this enhancement factor may
be partially reduced at ultra-relativistic intensities due to the
$\propto\sqrt{I}$ scaling compared to the $\propto I$ of moderately or sub-
relativistic intensities. In _Ferri et al._ [13] simulations at a higher
intensity of $\mathrm{7\times 10^{19}}\>\mathrm{W/cm^{2}}$ at 0.8 $\mathrm{\mu
m}$ wavelength ($\mathrm{\mathrm{a_{0}=5.7}}$) leads approximately to a three-
fold increase of the proton cut-off for single pulse irradiation of a 3
$\mathrm{\mu m}$ thick Al target and for a density scale-length of 0.6
$\mathrm{\mu m}$. In _Kumar et al._ [12] simulations at ultra-relativistic
intensities of $\mathrm{10^{20}}\>\mathrm{W/cm^{2}}$ and
$\mathrm{\mathrm{a_{0}=10.6}}$, very similar to the value of our experiment,
were investigated. In this intensity regime PIC simulations show that in the
presence of a pre-plasma consisting of a 10 $\mathrm{\mu m}$ linear density
ramp, the proton cut-off energy is expected to be approximately 8.5 MeV, while
without the pre-plasma, the cut-off energy is 5 MeV, showing again an
enhancement of approximately a factor 1.7.

In spite of these theoretical predictions, to date, a clear experimental
evidence of stochastic heating of fast electrons in proton acceleration
measurements is still lacking. This is mainly due to the complex interplay
typically occurring in real laser-target interactions, between plasma
formation at the front side and the deformation and disruption of the target
rear side, where ion acceleration mainly occurs, making changes to the fast
electron distribution function impossible or difficult to detect. In fact, in
most previous experiments, the generation of the required plasma gradient at
the vacuum-solid interface is achieved with uncontrolled precursor laser
irradiation due to Amplified Spontaneous Emission (ASE) [14] or nanosecond
irradiation [15] that have unpredictable, detrimental effects on the target
integrity.

The use of low energy, femtosecond pre-pulse for plasma generation allows a
controlled and predictable pre-plasma formation (see Methods). Previous
attempts to use this technique to show the effect of pre-plasma on fast
electron heating, were hindered by the relatively low intensity [16] or by the
presence of the ASE [17]. More recently [18], the same technique was used with
ASE-free interaction at ultra-relativistic intensities. A modest increase of
the proton cut-off energy was shown for 2 $\mu$m thick plastic targets with
numerical simulations suggesting self-focusing and pulse-front steepening.

We consider a similar framework to investigate ultra-relativistic interaction
with thin foil targets characterized by a vacuum-solid interface consisting of
a micrometer scale-length pre-plasma (see Methods). The interaction
configuration is modelled using realistic numerical PIC simulations that
account for the presence of the pre-plasma at the vacuum-target interface.
Simulations point at the presence of a strong stochastic heating of fast
electrons which is responsible for a robust enhancement of the TNSA field that
boosts the proton energy. Our measurements confirm this scenario, showing an
unprecedented three-fold increase of the cut-off energy of accelerated protons
compared to the measurements obtained without pre-plasma. While unfolding the
effect of a novel physical mechanism, our results also provide a recipe for
the optimum target thickness for a given pre-pulse timing, opening a novel
route to the optimization and control of efficient and robust TNSA ion
sources.

## Results

Experiment. The experiment was carried out as part of the development of a
laser driven particle acceleration programme for the development of a proton
beamline [19]. Details of the full experimental set up are given in the
Methods and a full characterization of the interaction conditions is given
elsewhere [19, 20, 21]. Here we focus on a set of proton spectra obtained from
irradiation of thin foil targets at a peak intensity of $\mathrm{2.4\times
10^{20}}\>\mathrm{W/cm^{2}}$, with and without a pre-plasma forming fs pre-
pulse. Proton spectra were obtained simultaneously with two different
detectors. The plots of Fig.1(left) shows two sets of proton spectra obtained
with the Thomson Parabola Spectrometer (TPS) (see Methods) on 10 $\mathrm{\mu
m}$ Al targets and on 25 $\mathrm{\mu m}$ Ti targets, respectively. The three
spectra with 10 $\mathrm{\mu m}$ Al targets were obtained without the pre-
plasma and with the standard contrast of our laser[19], and show a cut-off
energy around 6 MeV. The four spectra with 25 $\mathrm{\mu m}$ Ti targets were
obtained with a pre-plasma generated by the fs pre-pulse arriving 10.4 ns
before the main pulse, having a scale-length of 4.5 $\mu$m at the relativistic
critical density at the time of interaction (see Methods).

Figure 1: (left) Comparison of experimental proton spectra obtained from the
TPS for 10 $\mathrm{\mu m}$ thick Al targets without pre-plasma and 25
$\mathrm{\mu m}$ thick Ti target with the pre-plasma. Sample error bars on the
measured proton energy are shown for two spectra and are indicative of the
error-bars for all the curves. (right) Comparison of experimental spectra
obtained from deconvolution of the temporal evolution of the TOF signals from
a set of Titanium foil targets of different thickness of 2, 5, 12.5 and 25
$\mathrm{\mu m}$. All data were taken without the pre-plasma, except the
purple line for 25 $\mathrm{\mu m}$ with cut-off energy at 10 MeV. Sample
error bars on the measured proton energy are shown on the two spectra with 25
$\mu$m Ti targets.

The proton spectra of Fig.1 (left) obtained with the 25 $\mathrm{\mu m}$ Ti,
in the presence of the pre-plasma, show a cut-off energy in the range between
$\sim$9 and 12 MeV, namely up to a 2-fold increase with respect to the 10
$\mathrm{\mu m}$ Al target. These results clearly show that, in spite of the
significant increase of target thickness compared to the 10 $\mathrm{\mu m}$
Al, the presence of the pre-plasma significantly modifies the proton spectrum
obtained with the 25 $\mathrm{\mu m}$ Ti, with a strong enhancement of the
maximum energy of accelerated protons. A systematic increase of the signal at
low energy is also visible for the Ti data compared with Al targets. The 25
$\mathrm{\mu m}$ Ti plots of Fig.1(left) also demonstrate that the observed
increase of the cut-off energy in the presence of a pre-plasma is repeatable
and statistically significant.

To rule out the the possibility that this difference could be explained by the
different target material, shots were taken with Ti targets of different
thickness with and without the pre-plasma. Fig.1 (right) shows the proton
spectra obtained with the Time Of Flight (TOF) detector (see Methods) from the
irradiation of Ti targets with thickness ranging from 2 to 25 $\mathrm{\mu
m}$. All the data were taken without the pre-plasma, except the purple line
with cut-off at 10 MeV that was obtained with the pre-plasma. According to
these plots, the cut-off energy without the pre-plasma is maximum around 7 MeV
for a 5 $\mathrm{\mu m}$ Ti target and drops slightly to 6 MeV for 12.5 and to
4 MeV for 25 $\mathrm{\mu m}$ Ti.

This is consistent with the expected behaviour for TNSA accelerated protons
where cut-off energy is expected to decrease with increasing target
thickness[2, 3]. As already observed in our experimental conditions[19], the
cut-off energy also decreases for the thinnest target of 2 $\mathrm{\mu m}$
Ti. This is expected due to the finite temporal contrast of our laser system
at the ps level as previously discussed in details[19]. Similar cut-off
energies were obtained for 2.4 $\mathrm{\mu m}$ Al targets that showed a cut-
off energy of <4 MeV [19].

Remarkably, Fig.1 (right) shows that the cut-off energy of the 25 $\mathrm{\mu
m}$ Ti target with the pre-plasma exhibits the highest value of 10 MeV, a
value approximately 2.5 times the value obtained for the same target thickness
without the pre-plasma. These results confirm those obtained with the TPS, as
shown in Fig.1 (left) where a cut-off energy as high as 12 MeV is shown,
corresponding to a three-fold increase of the cut-off energy without the pre-
pulse for the same 25 $\mathrm{\mu m}$ Ti target thickness. We observe that
the highest cut-off energy with the pre-plasma was found only for the 25
$\mathrm{\mu m}$ thick Titanium target. In all other explored cases of thinner
targets, a lower or similar cut-off energy was found in the pre-plasma case,
compared with the no pre-plasma case, when the same target thickness was used.
As discussed below, this is due to the effect of the shock launched by the
pre-pulse used to generate the pre-plasma that, for thinner targets, reaches
the rear surface before the arrival of the main pulse, disrupting the sheath
field.

Modelling. We carried out extensive PIC simulations (see Methods) to model our
interaction conditions and ion acceleration measurements and to identify the
mechanism responsible for the observed three-fold increase of the proton cut-
off energy in the presence of a pre-plasma. We considered a pre-plasma scale-
length ranging from a minimum of 0.13 up to 5 $\mathrm{\mu m}$. The spectrum
of accelerated protons obtained from 3D PIC simulations for a 10 $\mathrm{\mu
m}$ Al for the case of no pre-plasma and for values of the plasma scale-length
of $L_{g}=0.13\>\mathrm{\mu m}$, $L_{g}=0.42\>\mathrm{\mu m}$ and
$L_{g}=1.27\>\mathrm{\mu m}$ are shown in Fig.2(a), while Fig.2(b) reports the
proton cut-off energy for the same cases. Here, the spectra are plotted when
the proton energy has reached the saturation value in the simulation.

According to these results, simulations predict a cut-off energy just below 5
MeV for the pre-plasma free case. This value increases up to 7.5 MeV for the
420 nm plasma scale-length and exceeds 10 MeV for the 1.27 $\mathrm{\mu m}$
scale-length. These simulations show a linear dependence of the cut-off energy
with the pre-plasma scale-length, predicting a cut-off energy increasing
rapidly with the pre-plasma scale-length, with a two-fold increase for a
scale-length of 1.27 $\mathrm{\mu m}$.

A comparison with our experimental results of Fig.1(left) for the case of 10
$\mathrm{\mu m}$ Al without pre-plasma, yielding a measured cut-off energy of
6.0 MeV, shows a similar value for a pre-plasma scale-length of approximately
100 nm. As anticipated, this scale-length is consistent with the presence of
the ps pedestal of our laser pulse. In fact, taking into account the ps
contrast our laser system [19], we expect the formation of a very short scale-
length pre-plasma a few ps before the laser peak. Considering a pre-plasma
expanding at the ion-acoustic wave speed, we estimate a density scale-length
at the critical density of the order of 100 nm, a value that is fully
consistent with the scale-length expected from the numerical simulations
results discussed above. Fig.2 confirms that in the case of a 10 $\mathrm{\mu
m}$ Al target the proton cut-off energy could exceed 10 MeV already for a
modest pre-plasma scale-length just above 1.27 $\mathrm{\mu m}$. However, we
stress here that in these and other similar PIC simulations, the presence of
the pre-plasma in front of the target is included ad hoc. In reality, the
formation of any pre-plasma generated on the target surface by either the ps
pedestal or by a custom pre-pulse, also involves the bulk of the target. The
shock generated at the target surface by any pre-pulse interaction eventually
reaches the rear surface of the target, disrupting the formation of the TNSA
field. This problem can only be neglected in the case of sufficiently thick
targets, where shock break-out at the target rear surface occurs after the
main pulse interaction, so that TNSA and ion acceleration take place before
any disruption of the target rear surface has occurred.

Figure 2: 3D PIC simulation results: (a) Proton energy spectra for a
$10\>\mathrm{\mu m}$ Al target for different pre-plasma scalength:
$L_{g}=0.13$ (blue), $0.42$ (red) and $1.27\>\mathrm{\mu m}$ (green). The
black curve corresponds to a sharp plasma-vacuum transition. (b) Cut-off
energy vs pre-plasma scalength for the cases in (a).

Shock break-out. Here we refer to the collisional shock generated by the
interaction of the pre-pulse at the target surface. We point out that shock
generation inside the target, driven by fast electrons, has also been
considered elsewhere [22] for ultrarelativistic intensities. In our
experimental conditions, this effect can be neglected for the pre-pulse, due
to the low intensity. As for the main pulse, break-out at the target rear face
of any shock produced by fast electrons would anyway take place tens of
picoseconds after proton acceleration [22], so this effect is not expected to
play a role in our measurements. The role of collisional shock generation
induced by pre-pulse irradiation at the target surface and propagation in the
bulk target was investigated with our hydrodynamic simulations (see Methods).
We found a shock wave propagation speed of 1.7 $\mathrm{\mu m}$/ns for the Al
target. Consequently, in the case of the 10 $\mathrm{\mu m}$ Al target, shock
break-out occurs 5.9 ns after pre-pulse irradiation, that is 4.5 ns before
main pulse irradiation that in our experiment occurs 10.4 ns after pre-pulse
irradiation. Therefore, in our experimental conditions, the enhancement seen
in PIC simulations for the 10 $\mathrm{\mu m}$ Al target with 1.27
$\mathrm{\mu m}$ scale-length could not be observed because disruption of the
rear surface would occur before main pulse irradiation.

In contrast, for Ti targets the hydrodynamic code predicts a shock propagation
speed of 1.5 $\mathrm{\mu m}$/ns that, for the 25 $\mathrm{\mu m}$ Ti target,
gives a shock break-out time at the rear surface of 16.6 ns after pre-pulse
irradiation. In our experimental conditions this corresponds to 6.2 ns after
the main pulse irradiation, thus allowing TNSA to take place without
perturbation. These considerations indeed show that the presence of the pre-
plasma for Al targets below 19 $\mathrm{\mu m}$ and Ti targets below
approximately 17.5 $\mathrm{\mu m}$ are subject to the effect of the shock
break-out at the rear face and disruption of the TNSA field. This also
explains why higher cut-off energy was only found for the 25 $\mathrm{\mu m}$
thickness Ti among all explored targets. Incidentally we observe that the TNSA
proton cut-off energy for sharp vacuum-target interface is indeed expected to
decrease with increasing target thickness, thus making such relatively thick
targets, in general, not ideal for efficient proton/ion acceleration. On the
other hand, both simulations and experimental data clearly show that for a
sufficiently large scale-length pre-plasma, the gain in cut-off energy due to
interaction with the pre-plasma can largely overcome such losses and lead to
record cut-off values.

Figure 3: 2D PIC simulations: (a, b) Laser electric field and longitudinal
density profile (black solid line), (c, d) electron density, (e, f)
longitudinal electric field and (g, h) electron kinetic energy density at
$t=0.11\>\mathrm{ps}$ for a Ti plasma with $L_{g}=5$ (left column) and
$0.13\>\mathrm{\mu m}$ (right column). The black dashed lines in (a, b) denote
the position of the critical density and the relativistic critical density.
The inset in (c) is an enlargement showing details of the electron cavities on
the target surface. The black solid lines in (e, f) represent the electric
field lineout in the middle of the simulation box.

Scaling. In order to investigate the interaction with a more extended pre-
plasma, we also carried out 2D PIC simulations (see Methods) with a scale-
length beyond the value of 1.27 $\mathrm{\mu m}$ explored above with 3D
simulations. Being computationally much more affordable that the 3D
simulations, 2D simulations allowed us to explore the role of longer pre-
plasma gradient in front of the target with scale-lengths of the order of
those predicted by hydrodynamic simulations (see Methods).

In Fig.3 the results from a simulation with a long pre-plasma
($L_{g}=5\>\mathrm{\mu m}$) are compared with results obtained with
$L_{g}=0.13\>\mathrm{\mu m}$ . When a long plasma gradient is present in front
of the target, the laser penetrates through the underdense plasma (Fig.3(a)),
where it undergoes self-focusing and steepening. When the laser reaches the
plasma region where $n_{c}\leq n\leq a_{0}n_{c}$, $n$ being the electron
density, $n_{c}=m_{e}\omega^{2}/(4\pi e^{2})$ being the critical density and
$a_{0}n_{c}$ relativistic critical density for $a_{0}>>1$[23], it expels
electrons forming semi-circular cavities (see Fig.3(c)) [24, 25]. It is
interesting to notice that nothing of this is observed when the laser reaches
the target surface when $L_{g}=0.13\>\mathrm{\mu m}$, corresponding to a
nearly sharp, almost unperturbed, vacuum-plasma transition (see Fig.3(b)). In
the latter case the ponderomotive force associated with the laser is not large
enough to push the target surface inward and bore a hole (Fig.3(d)) and the
laser is reflected back by the opaque target.

## Discussion

The complex laser dynamics in the plasma section across the critical density
for the case $L_{g}=5\>\mathrm{\mu m}$ leads to the formation of a standing
wave due to a superposition between the incident and reflected electromagnetic
wave [11, 12]. Here the role of the angle of incidence (15° in our case) has
been discussed recently [26] showing effective contribution of stochastic
heating for relatively long scale-length and even larger angles of incidence.
Signatures of the standing wave are clearly visible in Fig.3(e) while such
signatures are absent in Fig.3(f), where the longitudinal electric field is
plotted for $L_{g}=5$ and $0.13\>\mathrm{\mu m}$, respectively.

Figure 4: 2D PIC simulations results: (a) Electron distribution for density
scale-length $L_{g}=0.13$ (orange), $0.25$ (blue), $1.27$ (green), $3.38$
(red) and $5.07\>\mathrm{\mu m}$ (black). (b) Proton cut-off energy vs pre-
plasma scale length and (c) proton angular distribution for the cases in (a).

Electrons undergo stochastic motion in the standing wave [11, 27] resulting in
an overall more efficient electron heating, with a larger number of electrons
reaching higher energies (Figs.3(g) and (h)). The improved coupling between
the laser pulse and the electrons in the presence of an extended pre-plasma is
clearly visible in Fig.4(a), which shows the electron distribution in the
whole box after electrons have recirculated through the target, for different
values of $L_{g}$. All the spectra are made up by a cold and a hot electron
component. While the cold electrons have similar temperatures, the hot
electrons exhibit higher temperatures for longer plasma gradients, with
temperatures varying between $0.31\>\mathrm{MeV}$ for $L_{g}=0.13\>\mathrm{\mu
m}$ and $1.49\>\mathrm{MeV}$ for $L_{g}=5.07\>\mathrm{\mu m}$. As a
consequence, in the latter cases, the sheath field is enhanced and protons
attain higher energies (Fig.4(b)). As noted in Fig.2, the proton cut-off
energy increases linearly with $L_{g}$. However, from Fig.4(b) we observe that
for $L_{g}\geq 1.27\>\mathrm{\mu m}$, the maximum proton energy increases at a
slower rate compared to shorter pre-plasma gradients. A similar behaviour was
reported also in _Sentoku et al._ [11], but for higher values of $L_{g}$. On
one side this difference might be due to the use of a lower laser intensity in
their case. On the other side it is an indication that the plasma scale-length
used in our experimental conditions may already be in the saturation regime.

A comparison with the 3D simulations of Fig.2 shows that for the same scale-
length of $1.27\>\mathrm{\mu m}$, 2D simulations predict a significantly
higher cut-off energy of approximately 18.5 MeV, instead of 11 MeV. This is
expected due to the slower decay [28] of the sheath field in 2D compared to 3D
simulations, as also discussed below (see Methods). Therefore a quantitative
comparison with experimental values is not straightforward. However, it is
clear from these 2D simulation results that the presence of an extended pre-
plasma enables a more effective energy coupling to fast electrons that is at
the origin of the higher proton cut-off energy recorded for the
$25\>\mathrm{\mu m}$ Ti target, with $L_{g}\geq 5\>\mathrm{\mu m}$ and
yielding the observed three-fold increase of the maximum proton energy.
Finally, simulations seem to indicate that beside an increase of the cut-off
energy, the presence of a longer plasma gradient leads to a more collimated
proton beam (Fig.4(c)). This could not be verified in our experiment and will
be the subject of future investigations.

As a final remark, we observe that further proton cut-off enhancement can be
obtained optimizing the target thickness for a given pre-pulse timing relative
to the main pulse and for a given plasma density scale-length. In the case of
a pre-pulse to main pulse delay of 10.4 ns, as in our case, the minimum Ti
target thickness to avoid the detrimental effect on shock break-out on proton
acceleration is approximately 15 $\mu$m. Considering the target thickness
dependence of Fig.1(b), we can expect a final proton cut-off energy without
pre-plasma close to the 5.8 MeV measured for the 12.5 $\mu$m target. The
demonstrated three-fold increase in the presence of the pre-plasma would
therefore yield approximately 17 MeV of proton cut-off energy for a 15 $\mu$m
thick Ti target. Moreover, considering that, as anticipated, the 4.5 $\mu$m
scale-length used in our experiment is already in the saturation regime of the
fast electron heating mechanism, we can consider a shorter scale-length and a
correspondingly thinner target. In the case of a 5 $\mu$m and a pre-pulse to
main pulse delay of 3.3 ns, according to our model, we can expect a proton
cut-off energy in excess of 20 MeV.

Further investigation is needed to fully understand the role of laser pulse
self-focusing, pulse steepening and stochastic heating of fast electrons in
these interactions and their interplay. Our results demonstrate that with a
reliable technique to modify the vacuum-target interface, the effect of these
mechanisms on TNSA proton acceleration can be controlled to a much larger
extent than indicated by previous studies, with a potential strong impact on
applications of laser accelerated proton beams.

## Methods

Experimental set up. The experiment was carried out at the Intense Laser
Irradiation Laboratory using the ILIL-PW Ti:Sa laser installation [29] based
on Chirped Pulse Amplification (CPA)[30]. The laser pulse energy was up to 4.3
J on target with 27 fs pulse duration, with $\lambda=0.795\>\mathrm{\mu m}$
central laser wavelength. The 100 mm diameter beam was focused by an F/4.5
Off-Axis Parabolic (OAP) mirror with an angle of incidence of 15°. The focal
spot was elliptical, with an average diameter of 4.4 $\mathrm{\mu m}$ (FWHM).
The nominal intensity on target was up to $\mathrm{2.4\times
10^{20}}\>\mathrm{W/cm^{2}}$ with a relativistic dimensionless parameter
$\mathrm{\mathrm{a_{0}=8.55\times 10^{-10}\lambda[{\mathrm{\mu
m}}]I[W/cm^{2}]^{1/2}}=10.6}$ The target consisted of a thin foil mounted on a
remotely controlled motorized support with a sub-micrometer positioning
precision. Titanium and Aluminium foils were used, with thickness ranging from
2 to 25 $\mathrm{\mu m}$.

The main laser features and previous ion acceleration measurements are
described elsewhere [31, 19, 32], including details of the experimental set
up, laser contrast features and ion detectors.

Ion acceleration data were taken using the same set of ion beam diagnostics
already presented in [19], including a TPS and a TOF SiC detector. The SiC
detector was placed at a distance of 136 cm from the target rear side at 1.4
deg from the target normal and was filtered using a 48 $\mathrm{\mu m}$ thick
Al foil. The TPS and the TOF detector were used simultaneously so that a
cross-comparison of the signals obtained from the two devices was possible
[33, 34]. We observe that since the two detectors are not exactly on the same
line of sight, they are sampling slightly different portions of the
accelerated proton beam and, therefore, details of the spectra taken by the
two detectors may be different. We refer to the Supplementary Material for
additional information on the comparison between the two detectors.

In the experimental campaign described here we investigated the role of a
micrometer scale-length preplasma at the vacuum-solid interface. An effective
approach to plasma formation at the target surface is based on the use of a
low energy femtosecond pre-pulse [17], provided no effect of the ASE is
present. The plasma produced in these conditions is reproducible and an
accurate prediction of the density scale-length at the relevant electron
density and at the time of the main pulse irradiation is obtained with
hydrodynamic modelling. The shock generated by the this fs pre-pulse is also
be modelled and the break-out time of the shock front at the target rear face
is predicted with high accuracy. This is a crucial issue in proton
acceleration measurements: if the main pulse interaction occurs before the
shock break-out, proton acceleration measurements will not be affected by the
pre-pulse, and the role of the pre-plasma on the fast electrons can be thus
identified (see Supplemental information).

In our experiment, a femtosecond pre-pulse was set up on purpose by
introducing a small misalignment of the Pockels cell at the exit of the
regenerative amplifier. In this way, a small fraction of the pulse circulating
in the regenerative amplifier leaks through the Pockels cell every 10.4 ns.
The first pulse immediately before the main pulse is partially amplified in
the amplifier chain and, in particular, in the 3rd multipass amplifier, thanks
to the pumping timing scheme of the entire system, while all the other pulses
will be further attenuated. In such a way, a pre-pulse intensity on target,
relative to the main pulse, of $3\times 10^{-4}$ was obtained. The duration of
this pre-pulse was measured to be of $\sim 31\,\mathrm{fs}$ (this value was
taken into account for the hydrodynamic simulations). The intensity of the
pulse leaving the regenerative cavity one round-trip before the pre-pulse,
that is, $20.8\,\mathrm{ns}$ before the main pulse, was measured to be smaller
by at least a factor $10^{-6}$.

The effect of this pre-pulse on the target was modelled using hydrodynamic
simulations to calculate the expected profile of the pre-plasma in the
underdense region, at the time of interaction of the main pulse (see below).

Numerical simulations: Particle-in-cell. Numerical PIC simulations were
carried out using the code OSIRIS [35]. In the first set of 3D PIC simulations
a 10 $\mathrm{\mu m}$ Al foil was considered with a 32 nm thick H+ contaminant
layer at the back of the foil. The Al foil was assumed to be pre-ionised to
Al3+ with an initial temperature of 100 keV. We verified that different level
of target ionization and different initial electron temperature in the range
$10\,\mathrm{-}\,200\>\mathrm{KeV}$ did not lead to significant changes on the
final H+ energy. The density of the Al layer was 34.14 times the critical
density, nc = 1.76 $\times$ 1021 cm-3, while the density of the H+ layer on
the rear side was 1.14 times the critical density. The front of the target was
modelled following a hyperbolic tangent density profile centered at
$x=25.3\>\mathrm{\mu m}$ with scale-length $L_{g}$ to mimic the presence of
the pre-plasma. The laser intensity used for this first set of simulations was
$\mathrm{1.6\times 10^{20}}\>\mathrm{W/cm^{2}}$, corresponding to a normalized
vector potential $a_{0}=8.6$. The p-polarised pulse was focused on target to a
spot size $w_{0}$ of 2.2 $\mathrm{\mu m}$ and had a duration at FWHM of 27 fs.
The simulation box was 200 $\times$ 18 $\times$ 18 $\mathrm{\mu m}^{3}$ with a
cell size of 60 $\times$ 60 $\times$ 60 nm3, equivalent to approximately 13
points per laser wavelength, with 8 particles per cell. Simulations performed
with different resolution and number of particles per cell yielded similar
results.

A second set of 2D PIC simulations was performed to carry out a parametric
investigation of the role of target thickness and materials. For this set of
simulations we employed a laser intensity on target
$I=2.4\times\mathrm{10^{20}}\>\mathrm{W/cm^{2}}$
($\mathrm{\mathrm{a_{0}=10.6}}$) and we modelled the interaction of this pulse
with a 25 $\mathrm{\mu m}$ Ti target having a peak density
$n_{Ti}=56.9\>n_{c}$ and a longitudinal profile similar to that employed in
the previous simulation campaign (all the rest of the parameters are kept the
same). A simulation box of 607 $\times$ 26 $\mathrm{\mu m}^{2}$ with a cell
size of 30 $\times$ 30 nm2 equivalent to about 26 points per laser wavelength
has been used together with 36 particles per cell. We stress that, as it is
well known [36, 28], 2D simulations overestimate the final ion energy due to
the logarithmic growth in time of the proton cutoff energy. In our simulations
the cutoff energies shown in Fig.4b were measured at 1 ps after the beginning
of the laser interaction with the target. While a one-to-one comparison with
experimental results not possible, 2D numerical simulations still offer a
valid insight of the underpinning physics [36].

Numerical simulations: Pre-plasma. For the modelling of the pre-pulse
generated pre-plasma, we used the 2D Eulerian hydrocode POLLUX [37]. The code
models laser absorption via inverse bremsstrahlung and thermal transport via
flux-limited Spitzer-Harm conductivity. Ionization is calculated assuming
local thermodynamic equilibrium, while a perfect gas equation of state is used
for electrons.

According to the simulations, at the time of pre-pulse interaction with the
target, the plasma has a temperature of a few eV that rapidly cools down while
expanding. At the time on main pulse interaction, the scale-length of the pre-
plasma is 7 $\mathrm{\mu m}$ at the classical critical density (1.7 $\times
10^{21}$ cm-3) for the 800 nm laser light and 4.5 $\mathrm{\mu m}$ at the
relativistic critical density (1.8 $\times 10^{22}$ cm-3). Hydrodynamic
simulations also show that the shock launched by the fs pre-pulse propagates
through the target at a velocity of 1.5 $\mu$m/ns in the case of Titanium.

We refer to the Supplementary material for additional information on these
simulations.

## Conclusions

Our study demonstrates experimentally that a micrometer-size scale-length pre-
plasma can strongly enhance proton acceleration driven by ultraintense laser
interaction with not-so-thin foil targets, leading to a three-fold increase of
the cut-off energy with respect to case with no pre-plasma. For targets
sufficiently thick to prevent disruption of the TNSA accelerating field by
shock break-out at the rear side, the observed enhancement largely overcomes
the reduction of the cut-off energy due to the increased thickness. This
behaviour was modelled by detailed PIC simulations and confirmed by our
measurements. Our results open the possibility to significantly extend the
cut-off energy of robust TNSA based proton sources without plasma mirror or
ultra-thin targets and with a minor adjustment of a standard CPA laser front-
end.

## References

  * [1] Snavely, R. A. _et al._ Intense high-energy proton beams from petawatt-laser irradiation of solids.  _Phys. Rev. Lett._ 85, 2945, DOI: 10.1103/PhysRevLett.85.2945 (2000).
  * [2] Daido, H., Nishiuchi, M. & Pirozhkov, A. S.  Review of laser-driven ion sources and their applications.  _Reports Prog. Phys._ 75, 056401, DOI: 10.1088/0034-4885/75/5/056401 (2012).
  * [3] Macchi, A., Borghesi, M. & Passoni, M.  Ion acceleration by superintense laser-plasma interaction.  _Reviews of Modern Physics_ 85, 751–793, DOI: 10.1103/RevModPhys.85.751 (2013).
  * [4] Kaymak, V. _et al._ Boosted acceleration of protons by tailored ultra-thin foil targets.  _Scientific Reports_ 9, DOI: 10.1038/s41598-019-55011-2 (2019).
  * [5] Scott, G. G. _et al._ Dual Ion Species Plasma Expansion from Isotopically Layered Cryogenic Targets.  _Physical Review Letters_ 120, 204801, DOI: 10.1103/PhysRevLett.120.204801 (2018).
  * [6] Ferri, J. _et al._ Enhancement of laser-driven ion acceleration in non-periodic nanostructured targets.  _J. Plasma Phys._ 86, 1–12, DOI: 10.1017/S0022377819000898 (2019).
  * [7] Bayart, E. _et al._ Fast dose fractionation using ultra-short laser accelerated proton pulses can increase cancer cell mortality, which relies on functional PARP1 protein.  _Sci. Rep._ 9, 1–10, DOI: 10.1038/s41598-019-46512-1 (2019).
  * [8] Giulietti, A. (ed.) _Laser-Driven Particle Acceleration Towards Radiobiology and Medicine_ (Springer, 2016).
  * [9] Dover, N. P. _et al._ Effect of small focus on electron heating and proton acceleration in ultra-relativistic laser-solid interactions.  _Physical Review Letters_ 124, 1–6, DOI: 10.1103/PhysRevLett.124.084802 (2020).
  * [10] Steinke, S. _et al._ Acceleration of high charge ion beams with achromatic divergence by petawatt laser pulses.  _Phys. Rev. Accel. Beams_ 23, 1–7, DOI: 10.1103/PhysRevAccelBeams.23.021302 (2020).
  * [11] Sentoku, Y. _et al._ High-energy ion generation in interaction of short laser pulse with high-density plasma.  _Applied Physics B: Lasers and Optics_ 74, 207–215, DOI: 10.1007/s003400200796 (2002).
  * [12] Kumar, S., Gopal, K. & Gupta, D. N.  Proton acceleration from overdense plasma target interacting with shaped laser pulses in the presence of preplasmas.  _Plasma Physics and Controlled Fusion_ 61, DOI: 10.1088/1361-6587/ab216e (2019).
  * [13] Ferri, J., Siminos, E. & Fülöp, T.  Enhanced target normal sheath acceleration using colliding laser pulses.  _Commun. Phys._ 2, 1–8, DOI: 10.1038/s42005-019-0140-x (2019).
  * [14] Kaluza, M. _et al._ Influence of the Laser Prepulse on Proton Acceleration in Thin-Foil Experiments.  _Physical Review Letters_ 93, 045003, DOI: 10.1103/PhysRevLett.93.045003 (2004).
  * [15] McKenna, P. _et al._ Effects of front surface plasma expansion on proton acceleration in ultraintense laser irradiation of foil targets.  _Laser Part. Beams_ 26, 591–596 (2008).
  * [16] Yogo, A. _et al._ Laser prepulse dependency of proton-energy distributions in ultraintense laser-foil interactions with an online time-of-flight technique.  _Phys. Plasmas_ 14, 043104, DOI: 10.1063/1.2721066 (2007).
  * [17] Glinec, Y., Genoud, G., Lundh, O., Persson, A. & Wahlström, C. G.  Evolution of energy spectrum from laser-accelerated protons with a 100 fs intense prepulse.  _Appl. Phys. B Lasers Opt._ 93, 317–321, DOI: 10.1007/s00340-008-3172-6 (2008).
  * [18] Wang, P. _et al._ A simple way to introduce an adjustable femtosecond pre-pulse to enhance laser-driven proton acceleration.  _J. Phys. Conf. Ser._ 1350, 012063, DOI: 10.1088/1742-6596/1350/1/012063 (2019).
  * [19] Gizzi, L. A. _et al._ Light Ion Accelerating Line (L3IA): Test experiment at ILIL-PW.  _Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip._ 909, 160–163, DOI: 10.1016/j.nima.2018.03.016 (2018).
  * [20] Gizzi, L. A. _et al._ Intense proton acceleration in ultrarelativistic interaction with nanochannels.  _Phys. Rev. Res._ 2, 33451, DOI: 10.1103/PhysRevResearch.2.033451 (2020).
  * [21] Cristoforetti, G. _et al._ Laser-driven proton acceleration via excitation of surface plasmon polaritons into TiO2nanotube array targets.  _Plasma Phys. Control. Fusion_ 62, 114001, DOI: 10.1088/1361-6587/abb5e3 (2020).
  * [22] Shaikh, M. _et al._ Tracking ultrafast dynamics of intense shock generation and breakout at target rear.  _Phys. Plasmas_ 25, 2–6, DOI: 10.1063/1.5049815 (2018).
  * [23] Weng, S. M., Mulser, P. & Sheng, Z. M.  Relativistic critical density increase and relaxation and high-power pulse propagation.  _Phys. Plasmas_ 19, 022705, DOI: 10.1063/1.3680638 (2012).
  * [24] Wilks, S. C., Kruer, W. L., Tabak, M. & Langdon, A. B.  Absorption of ultra-intense laser pulses.  _Phys. Rev. Lett._ 69, 1383–1386, DOI: 10.1103/PhysRevLett.69.1383 (1992).
  * [25] Tsung, F., Tochitsky, S. Y., Haberbergher, D. J., Mori, W. B. & Joshi, C.  Co2 laser acceleration of forward directed mev proton beams in a gas target at critical plasma density.  _Journal of Plasma Physics_ 78, 373–382, DOI: 10.1017/S0022377812000189 (2012).
  * [26] Chopineau, L. _et al._ Identification of Coupling Mechanisms between Ultraintense Laser Light and Dense Plasmas.  _Phys. Rev. X_ 9, 011050, DOI: 10.1103/PhysRevX.9.011050 (2019).
  * [27] Nuter, R. _et al._ Influence of a preplasma on electron heating and proton acceleration in ultraintense laser-foil interaction.  _J. Appl. Phys._ 104, DOI: 10.1063/1.3028274 (2008).
  * [28] Babaei, J. _et al._ Rise time of proton cut-off energy in 2D and 3D PIC simulations.  _Phys. Plasmas_ 24, DOI: 10.1063/1.4979901 (2017).
  * [29] Gizzi, L. A. _et al._ A New Line for Laser-Driven Light Ions Acceleration and Related TNSA Studies.  _Appl. Sci._ 7, 984, DOI: 10.3390/app7100984 (2017).
  * [30] Strickland, D. & Mourou, G.  Compression of amplified chirped optical pulses.  _Opt. Commun._ 56, 219–221, DOI: 10.1016/0030-4018(85)90120-8 (1985).
  * [31] Gizzi, L. _et al._ Laser-plasma acceleration with flame and ilil ultraintense lasers.  _Applied Sciences_ 3, 559–580, DOI: 10.3390/app3030559 (2013).
  * [32] Gizzi, L. A. _et al._ Overview and specifications of laser and target areas at the Intense Laser Irradiation Laboratory.  _HIGH POWER LASER Sci. Eng._ DOI: 10.1017/hpl.2020.47 (2020).
  * [33] Gizzi, L. A. _et al._ Role of laser contrast and foil thickness in target normal sheath acceleration.  _Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers, Detect. Assoc. Equip._ 829, 144–148, DOI: 10.1016/j.nima.2016.01.036 (2016).
  * [34] Altana, C. _et al._ Investigation of ion acceleration mechanism through laser-matter interaction in femtosecond domain.  _Nucl. Instruments Methods Phys. Res. Sect. A-Accelerators Spectrometers Detect. Assoc. Equip._ 829, 159–162, DOI: 10.1016/j.nima.2016.02.016 (2016).
  * [35] Fonseca, R. A. _et al._ OSIRIS: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators.  _Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics)_ 2331 LNCS, 342–351, DOI: 10.1007/3-540-47789-6_36 (2002).
  * [36] Sgattoni, A., Londrillo, P., Macchi, A. & Passoni, M.  Laser ion acceleration using a solid target coupled with a low-density layer.  _Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys._ 85, 036405, DOI: 10.1103/PhysRevE.85.036405 (2012).
  * [37] Pert, G. J.  Two-dimensional hydrodynamic models of laser-produced plasmas.  _J. Plasma Phys._ 41, 263–280 (1989).

##  1 Acknowledgements

This project has received funding from the CNR funded Italian research Network
ELI-Italy ((D.M. No.631 08.08.2016) and from the L3IA INFN Experiment of CSN5.
We gratefully acknowledge support from CLF, STFC (UK) for in kind contribution
to the experimental set up described in this experiment. E.B. acknowledges the
OSIRIS Consortium, consisting of UCLA and IST, for the use of the OSIRIS 4.0
framework and the visXD framework. Simulations were performed on the
supercomputers ARCHER (EPCC, UK) under Plasma HEC Consortium EPSRC grants
EP/L000237/1 and EP/R029148/1 and Marconi-Broadwell (CINECA, Italy) under
PRACE allocation.

## Author contributions statement

L.A.G and L.L. conceived the experiment, L.A.G, L.L., F.Ba., F.Br., G.C.,
A.F., L.F., D.G., P.K. and D.P. conducted the experiment, L.A.G., L.L., F.Ba.,
F.Br., G.C., A. F., D.G., P.K., D.P., and P.T. analysed the results. E.B. and
P.J.B. carried out the numerical PIC simulations, L.L. carried out
hydrodynamic simulations, L.A.G. and E.B. prepared the manuscript. All authors
reviewed the manuscript.


# Post-merger evolution of double helium white dwarfs and distribution of
helium-rich hot subdwarfs

Jinlong Yu,1 Xianfei Zhang,2 Guoliang Lü1
1School of Physical Science and Technology, Xinjiang University, Urumqi,
830046, China
2Department of Astronomy, Beijing Normal University, Beijing 100875, China
E-mail<EMAIL_ADDRESS>E-mail<EMAIL_ADDRESS>E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The mergers of double helium white dwarfs are believed to form isolated
helium-rich hot subdwarfs. Observation shows that the helium-rich hot
subdwarfs can be divided into two subgroups based on whether the surface is
carbon-rich or carbon-normal. But it is not clear whether this distribution
directly comes from binary evolution. We adopt the binary population synthesis
(BPS) to obtain the population of single helium-rich hot subdwarfs according
to the channel of double helium white dwarfs merger. We find that the merger
channel can represent the two subgroups in the $T_{\rm{eff}}-\log g$ plane
related to different masses of progenitor helium white dwarfs. For $Z$ = 0.02,
the birth rates and local density of helium-rich hot subdwarf stars by the
mergers of two helium white dwarfs is $\sim 4.82\times 10^{-3}$ $\rm yr^{-1}$
and $\sim$ 290.0 $\rm kpc^{-3}$ at 13.7 Gyr in our Galaxy, respectively. The
proportion of carbon-rich and carbon-normal helium-rich hot subdwarfs are
32$\%$ and 68$\%$, respectively.

###### keywords:

binaries: close — stars: fundamental parameters — stars: white dwarfs

††pubyear: 2015††pagerange: Post-merger evolution of double helium white
dwarfs and distribution of helium-rich hot subdwarfs–Post-merger evolution of
double helium white dwarfs and distribution of helium-rich hot subdwarfs

##  1 Introduction

In general, over 95 percent of stars will end their lives as white dwarf stars
(WDs) and take a long time to cool down; thus, they are abundant and long-
lived. More than half of the stars are in binary systems (Han et al., 2020).
Therefore, after the interaction of binary stars, we can expect that many
binary star systems contain at least one white dwarf, or even both of them are
white dwarfs. According to different chemical compositions, white dwarfs can
be divided into helium-white dwarfs (HeWDs), carbon-oxygen white dwarfs (CO
WDs), and oxygen-neon-magnesium white dwarfs (ONeMg WDs). Almost all the
helium WDs result from binary star interaction because the evolution timescale
for a single star to evolve into a HeWD would exceed Universe’s age by far
(Bergeron et al., 1992; Han, 1998). According to theory and observation, there
are quite a few double HeWDs in the Galaxy, and they are also important
sources of gravitational waves (GW) (Bours et al., 2014). Some of the close
double HeWDs will be merged within a Hubble time due to orbital decay by GW
radiation. The merger of two HeWDs provides a natural mechanism for producing
extremely He-rich stars, e.g., He-rich hot subdwarfs. This channel is
considered to be one of the vital scenarios to explain the formation of single
He-rich hot subdwarfs (Webbink, 1984; Saio & Jeffery, 2000; Han et al., 2002,
2003; Zhang & Jeffery, 2012). Also, this merger channel are supported by the
observation distribution of He-rich hot subdwarfs in the Galactic (Luo et al.,
2019, 2020).

According to the surface effective temperature and spectra, hot subdwarfs are
traditionally grouped into subdwarf B (sdB), subdwarf OB (sdOB), and subdwarf
O (sdO). The classes sdOB and sdB are often merged into sdB owning to their
sdB-like spectra with weak He II (Heber, 2009, 2016). Some hot subdwarfs have
He-dominated spectra, which are known as He-rich hot subdwarfs. According to
the recently complete samples of hot subdwarfs in Gaia, about $10\%$ of the
hot subdwarf stars are He-rich hot subdwarfs (Geier et al., 2017; Geier et
al., 2019; Geier, 2020). He-rich hot subdwarfs are subdivided into carbon-rich
and nitrogen-rich (carbon-normal) according to their surface abundance of
carbon and nitrogen(Jeffery et al., 1997; Drilling et al., 2003; Stroeer et
al., 2007; Jeffery et al., 2021). By analyzing the spectrum, Hirsch (2009)
obtain carbon abundances (log $\beta_{C}$) and nitrogen abundances (log
$\beta_{N}$) that are given by mass fraction, i.e.,$\rm
log\beta_{i}=\frac{m_{i}}{m_{H}+m_{He}+m_{c}+m_{N}}$. They reported 16 carbon-
rich hot subdwarf stars with surface abundances $-$2.60 $\leq$ log $\beta_{C}$
$\leq$ $-$1.56 and 17 carbon-normal hot subdwarf stars with surface abundances
$-$4.52 $\leq$ log $\beta_{C}$ $\leq$ $-$4.14. Besides, the nitrogen
abundances of carbon-rich objects are in the range $-$3.60 $\leq$ log
$\beta_{N}$ $\leq$ $-$2.29, and the nitrogen abundances of carbon-normal stars
are in the range $-$2.81 $\leq$ log $\beta_{N}$ $\leq$ $-$2.20. Naslim et al.
(2010) obtain a carbon-rich single hot subdwarf stars with surface abundances
log $\beta_{C}$ $=$ $-$2.46 and four carbon-normal hot subdwarf stars with
surface abundances $-$4.34 $\leq$ log $\beta_{C}$ $\leq$ $-$3.56. Following
Hirsch (2009), similarly, we define log $\beta_{C}>-2.63$ for C-rich He-rich
hot subdwarfs, and log $\beta_{C}\leq-2.63$ for C-normal He-rich hot
subdwarfs. The values of $-$2.63 corresponed to solar surface carbon
abundance.

Figure 1: Evolutionary tracks of different double HeWDs merger remnants in the
Teff - log g plane with different metallicity, e.g., Z= 0.02, 0.01, and 0.001.
The red and yellow lines represent C-rich mergers and C-normal mergers,
respectively. The red dots and yellow triangles indicate the C-rich and
C-normal hot subdwarfs, which are observed by Hirsch (2009) and Naslim et al.
(2010). Figure 2: Models of double HeWDs merger remnants on the MWD-MWD
planes. Each panel shows the models with different metallicity, e.g., Z= 0.02,
0.01, and 0.001, respectively. The red circles indicate C-rich mergers, and
the yellow triangles indicate C-normal mergers. The dots on the background
show the masses of double HeWD obtained in the binary population synthesis.

In addition to the different abundances of elements, observations show that on
the $T_{\rm eff}-logg$ diagram, the two types of hot subdwarfs, C-rich and
C-normal, are obviously in two different temperature regions, and the dividing
line is roughly 43000 K. An essential question is whether there is a real gap
between C-rich and C-normal He-rich hot subdwarfs on the $T_{\rm eff}-logg$
diagram. The observational bias effects (e. g., lack of statistically
significant sample, errors of effective temperature, etc.) can make the gap,
and theoretical evolutionary channel also can directly form the gap. To
represent the mergers of double HeWDs, Zhang & Jeffery (2012) build models
that combine fast and slow accretion processes and their models are in
agreement with the observation of surface effective temperatures, surface
gravities, and abundances (in particular C and N). They show that the C-rich
and C-normal hot subdwarfs can form from the mergers of different masses of
progenitor HeWDs, e.g., the C-rich stars produced by high-mass HeWDs and
C-normal stars made by low-mass HeWDs, respectively. However, it’s difficult
to explain why these two subgroups exist due to only a few models calculated
in their work.

In this letter, from the perspective of binary star population synthesis
(BPS), we have extended the work of Zhang & Jeffery (2012). The results show
that the distribution of C-rich and C-normal hot subdwarfs is probably a
direct result of the helium white dwarfs’ merger channel. Section §2 describes
the stellar models and the simulation method of BPS. Results calculated from
BPS and conclusions are presented in Sections §3 and §4, respectively.

##  2 MODEL

###  2.1 Merger remnants

We calculate the merger remnants by taking advantage of the detailed stellar
evolution code Modules for Experiments in Stellar Astrophysics (MESA, version
12115, Paxton et al. 2011; Paxton et al. 2013, 2015). Following the methods of
Zhang & Jeffery (2012) and Zhang et al. (2014), we generate the HeWDs models
by evolving a 1.1 $\rm M_{\odot}$ zero-age main-sequence star (ZAMS) with
metallicity Z= 0.02, 0.01, and 0,001, then remove the complete hydrogen
envelope by a high mass-loss rate until the He core reach the required mass.
The naked He core evolves into the WD cooling phase straightly. Then we obtain
the HeWD models when the logarithmic surface luminosity log(L /$L_{\odot}$) =
$-$2\. The masses of HeWDs range from 0.2 $\rm M_{\odot}$ to 0.45 $\rm
M_{\odot}$ in steps of 0.05 $\rm M_{\odot}$ in our calculation.

---
Figure 3: Distribution of the He-rich hot subdwarfs formed via the double
HeWDs merger channel. The red dots and yellow triangles represent the C-rich
and C-normal hot subdwarfs, respectively. The observed samples are from Hirsch
(2009) and Naslim et al. (2010). In the top three panels, red and yellow lines
represent the C-rich and C-normal normalized number density in the effective
temperature range, respectively. The background grayscale distribution
represents the relative number distribution of merged white dwarfs for He-rich
hot subdwarfs.

Due to the SPH simulation of Dan et al. (2014) and Zhang & Jeffery (2012)
calculation, once the secondary (low-mass HeWD) are disrupted, the rest merger
process can decompose into two critical accretion processes before and after.
(1) a hot and fully convection envelope surrounds the primary formed by fast
accretion rate, which contains about half of the secondary’s masses; (2) a
disk created by the rest masses of the secondary, and then quickly disrupted
by expanding envelope, which are still in a relative fast accretion (Zhang et
al., 2014). Considering the convergence, we assume two similar steps to
represent such processes. We use a $10^{-3}$ $\rm M_{\odot}$ $\rm yr^{-1}$
accretion rate for half of the secondary’s masses to create a hot envelope and
make it thoroughly mix, then transfer the rest masses to the primary by the
same accretion rate to represent the disk disruption process. It should be
pointed out that the accretion rate we used is not consistent with the
accretion rate given by the SPH simulation. We use this accretion rate value
mainly to reproduce the similar temperature profile generated in the
simulation. After the accretion processes, we obtain a series of post-merger
models with different masses. Then, we evolve the post-merger models into the
WD phase again.

###  2.2 Evolutionary tracks of post-merger

Figure 1 shows the models with different metallicity, e.g., Z = 0.02, 0.01,
and 0.001, respectively. As shown by Zhang & Jeffery (2012), post-mergers will
undergo a series of inward shell helium flash processes, which will last for
about 0.6 - 2.5 million years. When the shell helium flash ends, the central
helium starts to burn and enters the hot subdwarf star’s main sequence phase.
For post-mergers with different masses, they have different surface carbon
abundances. Generally, the greater the mass, the higher the abundance of
carbon on the surface. As mentioned earlier, we define hot subdwarfs which
surface carbon abundance exceeds the solar abundance as C-rich type, and vice
versa as C-normal type. The evolutionary tracks in red and yellow represent
the C-rich and C-normal post-mergers, respectively. We compared the evolution
tracks with samples of Hirsch (2009) and Naslim et al. (2010). Except for the
evolution model with Z=0.001, the observed samples are in good agreement with
our models. The C-rich type is mainly located in the high-temperature area,
and the C-normal type is located in a relatively low-temperature area.
However, this distribution feature will also be affected by the mass
distribution of the white dwarfs before merging; hence we need to perform
further binary population synthesis calculations.

###  2.3 Binary Population Synthesis

We use BPS to investigate the mass distribution of progenitor double HeWDs.
Following Lü et al. (2012), Zhu et al. (2013), Lü et al. (2017), Zhu et al.
(2017) and Wang et al. (2018), we perform a Monte Carlo simulation of a
population with $10^{7}$ primordial binaries. Starting from zero-age main-
sequence (ZAMS), each binary system evolves using rapid binary star evolution
code BSE (Hurley et al., 2000, 2002). Finally, we select the pairs of double
HeWDs which will merge within a Hubble time to obtain the whole double HeWDs
populations.

In our Monte Carlo simulation, we assumed that all stars are binaries and
circular orbits. We adopt the initial mass function (IMF) of Miller & Scalo
(1979) to generate the primary, and in the mass range of 0.08 to 100 $\rm
M_{\odot}$. We obtain the secondary mass from a constant mass-ratio
distribution (Mazeh et al., 1992), i.e. $n(q^{{}^{\prime}})$ = 1, 0 $\leq$
$q^{{}^{\prime}}$ $\leq$ 1, where $q^{{}^{\prime}}$ = $M_{\rm d}$/$M_{\rm a}$.
We use the distribution of orbital separations given by Han (1998) for compact
binaries:

$an(a)=\begin{cases}0.07(a/a_{0})^{1.2},\qquad a\leq a_{0}\\\
0.07,\qquad\qquad a_{0}\leq a\leq a_{1},\end{cases}$ (1)

where $a_{0}=10\;\rm{R}_{\odot}$, $a_{1}=5.75\times 10^{6}\;\rm{R}_{\odot}$ =
0.13 $\rm pc$. According to the formula, this results in approximately 50
percent of stellar systems having an orbital period of more than 100 yr. These
systems are considered single star systems here(Chen & Han, 2009; Lü et al.,
2011, 2012). We adopt the same value of the IMF parameters, the mass-ratio
distribution, and the initial orbital separation distribution. Such parameters
were used in some previous works, e.g., hot subdwarfs (Han et al., 2002, 2003;
Zhang & Jeffery, 2012); double WD mergers in the Galaxy (Han, 1998; Zhang et
al., 2014); EL CVn stars (Chen et al., 2017), extremely low-mass WDs (ELM WDs)
in DDs (Li et al., 2019) and A-type subdwarfs (Yu et al., 2019). Figure 2
shows the mass distribution of double white dwarfs obtained from BPS and our
post-merger evolution models.

##  3 Result

Once we have the evolutionary tracks of post-mergers and the masses
distribution of progenitor double HeWDs, thus we can combine both results to
represent the He-rich hot subdwarf distribution. About 3000 HeWD+HeWD binaries
will undergo merger to produce hot subdwarf stars of the $10^{7}$ binary
systems. Combining the mass distribution of HeWDs, Figure 2 shows the types of
hot subdwarfs generated by merging HeWDs of different masses and only the
massive HeWDs can produce the C-rich hot subdwarfs. Furthermore, our results
show that to form C-rich hot subdwarfs, the total masses of the pre-mergers of
double helium white dwarfs must be large enough, and the masses of the primary
HeWDs must be large enough either.

Figure 3 shows the distribution of C-rich and C-normal He-rich hot subdwarf
stars formed via the HeWDs mergers with different metallicity in the
$T_{\rm{eff}}-\log g$ plane. The grayscale image shows the distribution of the
number of hot subdwarfs formed by merging helium-white dwarfs. We also show
the distribution of C-rich and C-normal subtypes in the temperature range. To
compare with the observations, Figure 3 shows the observation results of
Hirsch (2009) and Naslim et al. (2010). We can see that C-rich and C-normal
hot subdwarfs are in two different regions in the figure, roughly with 43000K
as the dividing line. Our calculation results reproduce this distribution
feature. It is worth pointing out that from the results of our helium white
dwarf merger model, the distribution characteristics of C-rich and C-normal
hot subdwarfs are not only affected by the masses of progenitors of HeWDs, but
also by their metallicities. As the metallicity decreases, the distribution
intervals of C-rich and C-normal gradually move closer. When the Z=0.001, the
distributions of the two subtypes in the temperature range are almost mixed
together, it isn’t easy to distinguish.

In Han (1998), the birth rates of the mergers of double HeWDs are $\sim
6.00\times 10^{-3}$ $\rm yr^{-1}$ in the Galaxy. We obtain a similar merger
rate of double HeWDs of $\sim 4.82\times 10^{-3}$ $\rm yr^{-1}$. Figure 4
shows the birth rate of C-rich and C-normal He-rich hot subdwarfs. We assume a
constant star formation rate of 5 $\rm M_{\odot}$ $\rm yr^{-1}$ over the past
13.7 Gyr, and the galaxy’s volume is 500 $\rm kpc^{3}$. By our calculation,
the local space density and birthrate of He-rich hot subdwarfs are $\sim$
290.0 $\rm kpc^{-3}$ and $\sim 4.82\times 10^{-3}$ $\rm yr^{-1}$ at 13.7 Gyr
for Z=0.02, respectively. Meanwhile, the proportion of C-rich and C-normal He-
rich hot subdwarfs are 32$\%$ and 68$\%$, respectively. We expect that the
number of C-rich and C-normal He-rich hot subdwarfs we may observe in the
Milky Way are $\sim 4.64\times 10^{4}$ and $\sim 9.86\times 10^{4}$,
respectively. For the metallicity Z=0.01 and 0.001, the birthrates are $\sim
6.23\times 10^{-3}$ $\rm yr^{-1}$ and $\sim 2.80\times 10^{-3}$ $\rm yr^{-1}$
at 13.7 Gyr; the local space densities are $\sim$ 379.0 $\rm kpc^{-3}$ and
$\sim$ 83.0 $\rm kpc^{-3}$ at 13.7 Gyr.

Figure 4: The local space density and birth rate of C-rich and C-normal He-
rich hot subdwarf stars formed via the double HeWDs merger channel. The solid
lines and dashed lines represent the C-normal and C-rich He-rich hot
subdwarfs, respectively.

##  4 Discussion and CONCLUSIONS

It is generally believed that hot subdwarfs are a special kind of stars formed
by the evolution of binary stars, and the details of their formation are still
not exact enough. He-rich hot subdwarfs, which are almost entirely composed of
helium, are a very special subtype. The merger of double helium white dwarfs
is considered an essential channel for forming He-rich hot subdwarfs.
Observedly, He-rich hot subdwarfs can be divided into two groups: C-rich and
C-normal in terms of temperature distribution. It is of great significance for
us to understand the formation of hot subdwarfs to clarify whether this
distribution feature is naturally derived from binary stars’ evolution or
caused by other reasons, such as observation effects. In this paper, we
constructed a series of double HeWD merge models following Zhang & Jeffery
(2012) and performed evolutionary calculations on post-mergers. Combined with
binary population synthesis, we have studied the distribution characteristics
of He-rich hot subdwarfs formed by merging channels.

We found that the He-rich hot subdwarf formed by merging channels is
distributed in two different regions on the $T_{\rm{eff}}-\log g$ plane. Its
interval range is consistent with the corresponding C-rich and C-normal He-
rich hot subdwarf observations. Meanwhile, we also obtain a similar surface
abundance of carbon. The main reason for this difference is that the C-rich
hot subdwarfs are derived from the merger of the larger-masses double helium
white dwarf star than the C-normal hot subdwarfs. Furthermore, as the metal
content increases from our calculation results, the separation between these
two regions will become larger. Therefore, the observed distributions of
C-rich and C-normal He-rich hot subdwarf in different temperature regions are
most likely the results of binary star evolution. Of course, due to the
limited number of samples, we look forward to more observations to confirm
this.

It is worth pointing out that the surface carbon abundance is directly related
to the merging model’s temperature profile and the element mixing process
caused by convection. The merger’s maximum temperature is related to the
accretion rate used in the merger process simulation. The maximum accretion
rate in our calculation can only reach $10^{-3}$ $\rm M_{\odot}$ $\rm
yr^{-1}$, which is lower than the accretion rate used by Zhang & Jeffery
(2012), so the surface carbon element is also lower than Zhang & Jeffery
(2012). However, the distribution of C-rich and C-normal regions is mainly
related to the mass of double helium white dwarfs. The accretion rate only
affects the abundance of carbon. Thus the accretion rate is not affected much
and does not affect our conclusion. In future work, we will further improve
our model calculations and look forward to more accurately reproducing the
formation and evolution of these He-rich hot subdwarfs. Besides, the other
formation channel of He-rich hot subdwarfs, e.g., late hot flasher scenarios,
are also worthy of consideration in future work. Although we can explain the
distribution characteristic of helium-rich hot subdwarf by the evolutionary
channel, we still cannot exclude the possibility of the lacking sufficiently
large observational sample. Therefore, a much larger observational sample
would be helpful to examine our results.

## Acknowledgements

This work received the generous support of the National Natural Science
Foundation of China, project Nos. U2031204, 11763007, 11863005, 12073006 and
U2031203. We would also like to express our gratitude to the Tianshan Youth
Project of Xinjiang No. 2017Q014.

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Bergeron et al. (1992) Bergeron P., Saffer R. A., Liebert J., 1992, ApJ, 394, 228
  * Bours et al. (2014) Bours M. C. P., et al., 2014, MNRAS, 438, 3399
  * Chen & Han (2009) Chen X., Han Z., 2009, MNRAS, 395, 1822
  * Chen et al. (2017) Chen X., Maxted P. F. L., Li J., Han Z., 2017, MNRAS, 467, 1874
  * Dan et al. (2014) Dan M., Rosswog S., Brüggen M., Podsiadlowski P., 2014, MNRAS, 438, 14
  * Drilling et al. (2003) Drilling J. S., Moehler S., Jeffery C. S., Heber U., Napiwotzki R., 2003, in Gray R. O., Corbally C. J., Philip A. G. D., eds, The Garrison Festschrift. p. 27
  * Geier (2020) Geier S., 2020, A&A, 635, A193
  * Geier et al. (2017) Geier S., Østensen R. H., Nemeth P., Gentile Fusillo N. P., Gänsicke B. T., Telting J. H., Green E. M., Schaffenroth J., 2017, A&A, 600, A50
  * Geier et al. (2019) Geier S., Raddi R., Gentile Fusillo N. P., Marsh T. R., 2019, A&A, 621, A38
  * Han (1998) Han Z., 1998, MNRAS, 296, 1019
  * Han et al. (2002) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., Ivanova N., 2002, MNRAS, 336, 449
  * Han et al. (2003) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., 2003, MNRAS, 341, 669
  * Han et al. (2020) Han Z.-W., Ge H.-W., Chen X.-F., Chen H.-L., 2020, Research in Astronomy and Astrophysics, 20, 161
  * Heber (2009) Heber U., 2009, ARA&A, 47, 211
  * Heber (2016) Heber U., 2016, PASP, 128, 082001
  * Hirsch (2009) Hirsch H. A., 2009, PhD thesis, Friedrich-Alexander University Erlangen-Nürnberg
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Jeffery et al. (1997) Jeffery C. S., Drilling J. S., Harrison P. M., Heber U., Moehler S., 1997, A&AS, 125, 501
  * Jeffery et al. (2021) Jeffery C. S., Miszalski B., Snowdon E., 2021, MNRAS, 501, 623
  * Li et al. (2019) Li Z., Chen X., Chen H.-L., Han Z., 2019, ApJ, 871, 148
  * Lü et al. (2011) Lü G., Zhu C., Wang Z., Huo W., Yang Y., 2011, MNRAS, 413, L11
  * Lü et al. (2012) Lü G.-L., Zhu C.-H., Postnov K. A., Yungelson L. R., Kuranov A. G., Wang N., 2012, MNRAS, 424, 2265
  * Lü et al. (2017) Lü G., Zhu C., Wang Z., Iminniyaz H., 2017, ApJ, 847, 62
  * Luo et al. (2019) Luo Y., Németh P., Deng L., Han Z., 2019, ApJ, 881, 7
  * Luo et al. (2020) Luo Y., Németh P., Li Q., 2020, ApJ, 898, 64
  * Mazeh et al. (1992) Mazeh T., Goldberg D., Duquennoy A., Mayor M., 1992, ApJ, 401, 265
  * Miller & Scalo (1979) Miller G. E., Scalo J. M., 1979, ApJS, 41, 513
  * Naslim et al. (2010) Naslim N., Jeffery C. S., Ahmad A., Behara N. T., Şahìn T., 2010, MNRAS, 409, 582
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Saio & Jeffery (2000) Saio H., Jeffery C. S., 2000, MNRAS, 313, 671
  * Stroeer et al. (2007) Stroeer A., Heber U., Lisker T., Napiwotzki R., Dreizler S., Christlieb N., Reimers D., 2007, A&A, 462, 269
  * Wang et al. (2018) Wang T., Li L., Zhu C., Wang Z., Wang A., Wu Q., Liu H.-S., Lü G., 2018, ApJ, 863, 17
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Yu et al. (2019) Yu J., et al., 2019, ApJ, 885, 20
  * Zhang & Jeffery (2012) Zhang X., Jeffery C. S., 2012, MNRAS, 419, 452
  * Zhang et al. (2014) Zhang X., Jeffery C. S., Chen X., Han Z., 2014, MNRAS, 445, 660
  * Zhu et al. (2013) Zhu C., Lü G., Wang Z., 2013, ApJ, 777, 23
  * Zhu et al. (2017) Zhu C., Lü G., Wang Z., 2017, ApJ, 835, 249


In this paper, we proposed a hybrid of PG and inexact projected regularized
Newton method for solving the fused $\ell_{0}$-norms regularization problem
(1). This hybrid framework fully exploits the advantages of PG method and
Newton method, while avoids their disadvantages. We proved a global
convergence by means of KL property, without assuming the uniformly positive
definiteness of the regularized Hessian matrix, and also obtained a
superlinear convergence rate under a Hölderian local error bound on the set of
the second-order stationary points, without assuming the local minimality of
the limit point.

All PGiPN, ZeroFPR and PG have employed the polynomial-time algorithm for
computing a point in the proximal mapping of $g$ when $B=\widehat{B}$, that we
developed in this paper. Numerical tests indicate that our PGiPN not only
produce solutions of better quality, but also requires 2-3 times less running
time than PG and ZeroFPR, where the latter mainly attributes to our subspace
strategy when applying the projected regularized Newton method to solve the
problems.

## References

  * [1] M. Ahookhosh, A. Themelis, and P. Patrinos, A Bregman forward-backward linesearch algorithm for nonconvex composite optimization: superlinear convergence to nonisolated local minima, SIAM Journal on Optimization, 31 (2021), pp. 653–685.
  * [2] A. Aravkin, M. P. Friedlander, F. J. Herrmann, and T. Van Leeuwen, Robust inversion, dimensionality reduction, and randomized sampling, Mathematical Programming, 134 (2012), pp. 101–125.
  * [3] H. Attouch, J. Bolte, P. Redont, and A. Soubeyran, Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality, Mathematics of Operations Research, 35 (2010), pp. 438–457.
  * [4] H. Attouch, J. Bolte, and B. F. Svaiter, Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized Gauss-Seidel methods, Mathematical Programming, 137 (2013), pp. 91–129.
  * [5] G. Bareilles, F. Iutzeler, and J. Malick, Newton acceleration on manifolds identified by proximal gradient methods, Mathematical Programming, 200 (2023), pp. 37–70.
  * [6] H. H. Bauschke, J. M. Borwein, and W. Li, Strong conical hull intersection property, bounded linear regularity, Jameson’s property (g), and error bounds in convex optimization, Mathematical Programming, 86 (1999), pp. 135–160.
  * [7] D. P. Bertsekas, Projected Newton methods for optimization problems with simple constraints, SIAM Journal on Control and Optimization, 20 (1982), pp. 221–246.
  * [8] , Nonlinear programming, Journal of the Operational Research Society, 48 (1997), pp. 334–334.
  * [9] W. Bian and X. Chen, A smoothing proximal gradient algorithm for nonsmooth convex regression with cardinality penalty, SIAM Journal on Numerical Analysis, 58 (2020), pp. 858–883.
  * [10] T. Blumensath and M. E. Davies, Iterative thresholding for sparse approximations, Journal of Fourier Analysis and Applications, 14 (2008), pp. 629–654.
  * [11] , Normalized iterative hard thresholding: Guaranteed stability and performance, IEEE Journal of selected topics in signal processing, 4 (2010), pp. 298–309.
  * [12] J. Bolte, S. Sabach, and M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Mathematical Programming, 146 (2014), pp. 459–494.
  * [13] J. V. Burke and J. J. Moré, On the identification of active constraints, SIAM Journal on Numerical Analysis, 25 (1988), pp. 1197–1211.
  * [14] H. Davenport and A. Schinzel, A combinatorial problem connected with differential equations, American Journal of Mathematics, 87 (1965), pp. 684–694.
  * [15] J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani, Pathwise coordinate optimization, The Annals of Applied Statistics, 1 (2007), pp. 302–332.
  * [16] K. K. Herrity, A. C. Gilbert, and J. A. Tropp, Sparse approximation via iterative thresholding, in 2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings, vol. 3, IEEE, 2006, pp. III–III.
  * [17] S. W. Jewell, T. D. Hocking, P. Fearnhead, and D. M. Witten, Fast nonconvex deconvolution of calcium imaging data, Biostatistics, 21 (2020), pp. 709–726.
  * [18] H. Jiang, S. Luo, and Y. Dong, Simultaneous feature selection and clustering based on square root optimization, European Journal of Operational Research, 289 (2021), pp. 214–231.
  * [19] C. Kanzow and T. Lechner, Efficient regularized proximal quasi-Newton methods for large-scale nonconvex composite optimization problems, arXiv preprint arXiv:2210.07644, (2022).
  * [20] R. Killick, P. Fearnhead, and I. A. Eckley, Optimal detection of changepoints with a linear computational cost, Journal of the American Statistical Association, 107 (2012), pp. 1590–1598.
  * [21] S. R. Land and J. H. Friedman, Variable fusion: A new adaptive signal regression method, tech. rep., Technical Report 656, Department of Statistics, Carnegie Mellon University, 1997.
  * [22] J. D. Lee, Y. Sun, and M. A. Saunders, Proximal Newton-type methods for minimizing composite functions, SIAM Journal on Optimization, 24 (2014), pp. 1420–1443.
  * [23] X. Li, D. Sun, and K.-C. Toh, On efficiently solving the subproblems of a level-set method for fused Lasso problems, SIAM Journal on Optimization, 28 (2018), pp. 1842–1866.
  * [24] J. Liu, S. Ji, and J. Ye, SLEP: Sparse learning with efficient projections, Arizona State University, 6 (2009), p. 7.
  * [25] J. Liu, L. Yuan, and J. Ye, An efficient algorithm for a class of fused Lasso problems, in Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining, 2010, pp. 323–332.
  * [26] R. Liu, S. Pan, Y. Wu, and X. Yang, An inexact regularized proximal Newton method for nonconvex and nonsmooth optimization, arXiv preprint arXiv:2209.09119v5, (2023).
  * [27] Z. Lu, Iterative hard thresholding methods for $\ell_{0}$ regularized convex cone programming, Mathematical Programming, 147 (2014), pp. 125–154.
  * [28] Z. Lu and Y. Zhang, Sparse approximation via penalty decomposition methods, SIAM Journal on Optimization, 23 (2013), pp. 2448–2478.
  * [29] C. Molinari, J. Liang, and J. Fadili, Convergence rates of Forward–Douglas–Rachford splitting method, Journal of Optimization Theory and Applications, 182 (2019), pp. 606–639.
  * [30] B. S. Mordukhovich, X. Yuan, S. Zeng, and J. Zhang, A globally convergent proximal Newton-type method in nonsmooth convex optimization, Mathematical Programming, 198 (2023), pp. 899–936.
  * [31] S. Pan, L. Liang, and Y. Liu, Local optimality for stationary points of group zero-norm regularized problems and equivalent surrogates, Optimization, 72 (2023), pp. 2311–2343.
  * [32] R. A. Poliquin and R. T. Rockafellar, A calculus of prox-regularity, J. Convex Anal, 17 (2010), pp. 203–210.
  * [33] R. T. Rockafellar, Convex analysis, Princeton university press, 1970\.
  * [34] R. T. Rockafellar and R. J.-B. Wets, Variational analysis, vol. 317, Springer, 2009.
  * [35] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268.
  * [36] M. Sharir, Davenport-schinzel sequences and their geometric applications, in Theoretical Foundations of Computer Graphics and CAD, Springer, 1988, pp. 253–278.
  * [37] S. Sra, Scalable nonconvex inexact proximal splitting, Advances in Neural Information Processing Systems, 25 (2012).
  * [38] L. Stella, A. Themelis, and P. Patrinos, Forward–backward quasi-Newton methods for nonsmooth optimization problems, Computational Optimization and Applications, 67 (2017), pp. 443–487.
  * [39] A. Themelis, M. Ahookhosh, and P. Patrinos, On the acceleration of forward-backward splitting via an inexact Newton method, in Splitting Algorithms, Modern Operator Theory, and Applications, Springer, 2019, pp. 363–412.
  * [40] A. Themelis, L. Stella, and P. Patrinos, Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone linesearch algorithms, SIAM Journal on Optimization, 28 (2018), pp. 2274–2303.
  * [41] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight, Sparsity and smoothness via the fused Lasso, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67 (2005), pp. 91–108.
  * [42] K. Ueda and N. Yamashita, Convergence properties of the regularized Newton method for the unconstrained nonconvex optimization, Applied Mathematics and Optimization, 62 (2010), pp. 27–46.
  * [43] S. J. Wright, R. D. Nowak, and M. A. Figueiredo, Sparse reconstruction by separable approximation, IEEE Transactions on Signal Processing, 57 (2009), pp. 2479–2493.
  * [44] F. Wu and W. Bian, Accelerated iterative hard thresholding algorithm for $l_{0}$ regularized regression problem, Journal of Global Optimization, 76 (2020), pp. 819–840.
  * [45] Y. Wu, S. Pan, and X. Yang, A regularized Newton method for $\ell_{q}$-norm composite optimization problems, SIAM Journal on Optimization, 33 (2023), pp. 1676–1706.
  * [46] M.-C. Yue, Z. Zhou, and A. M.-C. So, A family of inexact SQA methods for non-smooth convex minimization with provable convergence guarantees based on the Luo-Tseng error bound property, Mathematical Programming, 174 (2019), pp. 327–358.
  * [47] S. Zhou, L. Pan, and N. Xiu, Newton method for $\ell_{0}$-regularized optimization, Numerical Algorithms, 88 (2021), pp. 1541–1570.