11institutetext: INAF-OATS, Via G.B. Tiepolo 11, 34143 Trieste, Italy
11email<EMAIL_ADDRESS>22institutetext: Dipartimento di Fisica ”Enrico
Fermi”, Università di Pisa and INFN – Sezione di Pisa, Largo B. Pontecorvo 3,
56127 Pisa, Italy
22email<EMAIL_ADDRESS>
# Revisiting V1309 Sco 2008 outburst spectra

Observational evidence for theoretical modeling of stellar mergers

Elena Mason 11 Steven N. Shore 22

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. V1309 Sco is the only certain noncompact stellar merger, due to its
indisputable preoutburst light curve matching that of a contact binary of
almost equal mass stars. Therefore, anything that can be deduced from the
existing observations serves as benchmark constraints for models.

Aims. We present some observational evidences to guide future hydrodynamical
simulations and common envelope studies.

Methods. Using archive spectra taken at high and mid spectral resolution
during the V1309 Sco outburst and late decline, together with the inferential
methods we developed to study nova ejecta through panchromatic high resolution
spectroscopic follow ups, we constrain the physical state, structure, dynamics
and geometry of the transient originated in the stellar merger.

Results. We found that the emitted spectra arise from two distinct
contributions: matter expelled during the 2008 outburst and circumbinary gas
produced during historic mass loss episodes. These two components likely have
orthogonal geometry with the 2008 mass loss displaying a dust-laden bipolar
ejecta produced by a time limited rapidly accelerating wind and the
circumbinary gas having a donut-like shape. A central source powers them both,
having produced a fluorescent light pulse, but we cannot precisely determine
the time it started or its spectral energy distribution. We can, however,
place its upper energy cutoff at about 54 eV and the bulk of its emission at
$<$20 eV. We also know that the central source turned off within months from
the outburst and before the ejecta turned optically thin.

###### Key Words.:

##  1 Introduction

Discovered as a transient in 2008 (e.g. Nakano et al. 2008), V1309 Sco was
first thought to be a microlensing event (Bensby, private communication). It
was soon classified as a classical nova (Naito & Fujii 2008). Monitored during
its early decline through high resolution spectroscopy at the VLT+UVES, V1309
Sco was suggested to be a so called ”red-nova” (Mason et al. 2010) and
eventually shown to be a non-compact stellar merger (Tylenda et al. 2011).
Here, we re-visit the UVES spectra published in Mason et al. (2010), together
with the late decline sequence of VLT+XShooter spectra published in Kaminski
et al. (2015).

Over the past several years we have studied panchromatic high resolution
spectra of classical novae unveiling the dynamic of their ejecta. Our intent
is to use that experience to discern the dynamics and the physics of the
stellar merger, and to produce empirical constraints to guide numerical
simulations.

##  2 The ESO archive data

The V1309 Sco spectra analyzed herein are in the European Southern Observatory
(ESO) archives. These correspond to six UltraViolet Echelle Spectrograph
(UVES) epochs taken with the DIC1(346+580) and DIC2(437+860) settings which,
together, cover the wavelength range $\sim$3000-10300 Å in each epoch; plus 3
XShooter epochs covering the wavelength range $\sim$3000 Å - 2.4 $\mu$m each.
The UVES data were obtained during regular service mode operations, where
downloaded from the archive111http://archive.eso.org/eso/eso_archive_main.html
and reduced on a local computer using the esorex pipeline (version 6.1.3) with
gasgano interface (version 2.4.8). The data were flux calibrated using the
”master response” function produced by the quality control group, rather than
the ”instrument response” function that can be generated with the
spectrophotometric standard star observation of the night (if taken), because
experience showed that the spectra are better calibrated in the whole spectral
range when the master response is used. The XShooter data were, instead, taken
during commissioning and science
verification222https://www.eso.org/sci/facilities/paranal/instruments/xshooter/news.html
and links therein (bottom page).. The reduced spectra are those already
published in Kaminski et al. (2015) that provide the details of the data
reduction process. Here we add that the extraction of the third XShooter
spectrum was done using a reduced extraction window, smaller than the spectrum
spatial full width half maximum (FWHM), to avoid the part contaminated by a
nearby star (the seeing allowed such a crude deblending). Adopting the time of
the transient’s discovery, 2008 Sep. 2.46 UT, as $t_{0}$, the UVES spectra
sample the first decline phase, the plateau and the subsequent steep decline
phase displayed by the V band light curve (Figure 1 in Mason et al. 2010). The
XShooter spectra, taken about one year later, monitor the late decline (e.g.
Figure 2 in Kaminski et al. 2015). The observing date and the age of the
transient at each epoch are in Table 1.

Table 1: The date of the UVES and XShooter observations, together with the age of the transient at each epoch (considering 2008 Sep. 2.46, i.e. MJD 54711.46, the time of discovery, the $t_{0}$ epoch). Obs. Date (UT) | Epoch # | Age (d) | Instrument
---|---|---|---
2008-09-13 | 1 | 10.6 | UVES
2008-09-18 | 2 | 15.6 | UVES
2008-09-20 | 3 | 18.5 | UVES
2008-09-28 | 4 | 25.6 | UVES
2008-10-08 | 5 | 35.6 | UVES
2008-10-20 | 6 | 47.6 | UVES
2009-05-04 | 7 | 243.7 | XShooter
2009-08-13 | 8 | 344.7 | XShooter
2009-09-29 | 9 | 391.6 | XShooter

##  3 Analysis of the UVES data set

Figure 1: GASS/LAB H I 21 cm temperature profile (solid blue line and red dots
for the GASS and LAB measurements, respectively), compared with the K I
$\lambda$7699 (black line) and Na I D (gray line) interstellar absorptions.
The units of the right vertical axis are erg/cm2/s/Å. Telluric absorptions are
marked with the $\oplus$ symbol. Velocities are in the Local Standard of Rest
(LSR) reference frame.

###  3.1 Reddening and distance

The reddening for V1309 Sco cannot be inferred from the Na I or Ca II
interstellar absorption lines since these are saturated and severely blended
with low velocity absorption structures from V1309 Sco itself. However, this
is not the case for the K I absorption doublet which is clean and traces the
same interstellar gas as Na I. Fig.1 compares the Galactic-All-Sky-
Survey/Leiden-Argentine-Bonn (GASS/LAB) 21 cm temperature profile (Kalberla et
al. 2005) with the absorption line profile from the K I $\lambda$7698.96 line
in the first UVES spectrum. The comparison suggests that the bulk of the
neutral H mapped by GASS/LAB is between V1309 Sco and the observer and that we
can integrate the whole emission profile between $-40$ and $+$65 km/s for a
reliable estimate of the reddening. We derive a H I column density of
$\sim$3.7$\times$1021 cm-2 and E(B-V)=0.65$\pm$0.01 mag.

The K I and Na I profiles show evidence for an additional absorption component
centered at v${}_{LSR}\sim-$60 km/s which does not have a corresponding
feature in the GASS/LAB profile. We favor the interpretation that they are of
interstellar origin since we observe them in all the ground state transitions
(e.g. Na I, K I, etc) but in none of the transitions from ionized atoms (e.g.
Fe II, Sc II, etc) that certainly originate in the transient’s gas. The
difference between the metal atomic absorptions and the H I profile can be
explained by the patchy line of sight toward the Galactic center and the
coordinate offset between the V1309 Sco the H I surveys pointings, that can
exclude or include intervening clouds. The LAB nearest pointing is offset by
$\geq$0.21o, while the effective beamsize is 0.266o. Expanding the range for
the integral of the GASS/LAB profile to $-80$ to $+65$ km/s to include such
additional absorptions, does not change the resulting H I column density and
reddening.

Using different diagnostics, however, the reddening changes depending on
whether the $-60$ km/s absorptions are taken into account. In particular,
measuring the K I equivalent width, EW, and adopting the Munari & Zwitter
(1997) empirical relation, we get E(B$-$V)$\simeq$0.60$\pm$0.02 mag from the
average of the measurements from the six UVES spectra for the absorptions
centered at $\sim 0$ km/s alone, and E(B$-$V)$\simeq$0.84$\pm$0.04 when
including the $\sim-60$ km/s absorptions.

Finally, we can obtain a reddening estimate using the EW of diffuse
interstellar bands (DIBs). Although most of the strongest DIBs are either
severely blended with absorption features from the transient or in the
wavelength gaps of the UVES setup, we can measure the EW of the DIB at 6613.3
Å and 6283.8 Å, albeit with a lower confidence because of the broad profile
and the imperfect telluric correction. Averaging the measurements from all six
UVES spectra we get EW$\simeq$0.15$\pm$0.03 Å for the 6613 Å DIB and
EW$\simeq$0.64$\pm$0.14 for the 6284 Å DIB. They correspond to E(B-V) in the
range 0.67-0.73 and 0.54-1.03, respectively, depending on the calibration
(e.g. Jenniskens & Desert 1994, Friedman et al. 2011). The two DIBs match, in
velocity, the positions of the $\sim 0$ km/s clouds in the LSR, possibly
indicating a smaller column density for the cloud at $\sim-60$ km/s.

Noting that these three methods are completely independent, we adopt
0.65$\leq$E(B-V)$\leq$0.84 mag and work with these two extremes, whenever the
dereddening correction is needed to the analysis.

The additional interstellar gas cloud at LSR velocity $\sim-60$ km/s, has an
interesting implication for the distance of V1309 Sco. The V1309 Sco
coordinates point toward the Galactic center ($l\sim 360$ and $b\sim-3$ deg),
at which the Dame et al. (2001) CO maps show the ”3 kpc ring” molecular gas
velocities compatible with the observed $-60$ km/s. This implies that the
transient is behind the ”3 kpc molecular ring” and therefore at about $\sim$ 5
kpc from the Sun. This is somewhat farther than the 3-3.5 kpc previously
estimated (e.g. Tylenda et al. 2011).

Figure 2: Examples of line profiles and their evolution across the UVES
spectra. Each subpanel is for a different transition (specified in the
subpanel title); early epochs are on top, and later ones in the bottom of each
subpanel. The vertical shaded area marks the velocity range (-200,-115) km/s
where the persistent absorption is located. Velocities are in the heliocentric
rest frame. Y-axis units are in erg/cm2/s/Å.

Figure 3: Line profiles of the forbidden doublets [O I] and [Ca II] (the
dashed line is for the weakest of the transitions in the doublet). Times
increase from top to bottom within each subpanel; while the shaded area marks
the velocity range (-200,-115) km/s where persistent absorptions are observed.
Velocities are in the heliocentric rest frame. Y-axis units are in
erg/cm2/s/Å.

###  3.2 The line profiles, their evolution and the inferred structure

Figures 2 and 3 show a selection of line profiles and their evolution across
the six UVES epochs. We distinguish two different portions of each line in
velocity space: 1) the range ($-$360,$-$200) km/s, which we call the ”high
velocity component” (HVC); and 2) the range (–200,–50) km/s which we call the
”low velocity component” (LVC) and which displays what looks like a persistent
absorption in the velocity range (–200,–115) km/s (the gray band in Figs. 2
and 3). These velocity limits are approximate, since they may vary slightly
with the transition. The two ranges manifest different evolution. Their line
ratios are also substantially different, suggesting different opacity and gas
conditions.

The HVC is constantly evolving in our sequence of spectra. It appears as a
weak, shallow and asymmetric absorption that progressively increases in
strength from epoch 1 to 4 and eventually turns into a strong emission profile
(from epoch 5 on), which peaks at $\sim-$210 km/s and has an extended blue
wing. The strength of the emission relative to the continuum, also increases
with time. This evolution is particularly evident in relatively isolated
transitions (i.e. free from confounding blends) that are either resonance
lines (e.g. Na I and K I panels in Fig. 2) or have metastable lower energy
levels (e.g. Sc II panel in Fig. 2). Recombination lines, such as the early H
Balmer and He I transitions, are always in emission (or emission dominated),
with the line strength steadily increasing from epoch 1 to 6 and later
(XShooter spectra), relative to the continuum. Because of its velocity, we
associate the HVC with gas expelled from the system during the outburst. Its
evolution is also consistent with gas whose physical conditions of opacity,
ionization, and possibly density, change because of expansion. The shallow
absorption initially increases in strength because of the of a propagating
recombination wave (Shore et al. 2011, Shore 2012); it then turns into
emission once the optical depth has substantially decreased. Similarly, the
strength of the H and He I high velocity components grows constantly as self
absorption effects decrease.

Figure 4: Selected ranges in the U-band spectrum showing the ”curtain”. This
is mostly produced by ground transitions from neutral and singly ionized
metals and completely absorbing the continuum blue-ward of the Balmer jump.
The top panel shows the whole wavelength range 3200-4100 Å; the middle and the
bottom panels show zoomed in portions of the curtain. Black, blue, yellow,
light-blue, red and dotted-gray are for epoch 1, 2, 3, 4, 5 and 6,
respectively. Figure 5: Comparison between the observations (black line, epoch
1 spectrum, smoothed with a running box of 5) and the curtain simulation
produced as described in the text (light blue line). Units of the y axis are
erg/s/cm2/Å. The simulation has been arbitrarily scaled.

The LVC consists of a ubiquitous absorption component and a relatively weak
emission that together resemble a P Cyg profile (where blending with other
transitions is not affecting significantly the whole profile). On close
inspection, the absorption component is persistent and almost invariant in
width, profile and, to a less extent, depth, especially in the resonance
transitions in the first four epochs. Moreover, the 3200-3700 Å wavelength
range displays a plethora of neutral and singly ionized metal ground and low
excitation energy transitions that form a curtain absorbing the continuum
radiation during all six UVES epochs (see Fig.9). Those absorptions do not
change in time, what changes is the continuum. This is why we dubbed them
”persistent absorptions”. They are also observed in the [O I] and [Ca II]
doublets and this has interesting physical implications that we discuss in
Sect. 3.4. To test the effective presence of a curtain we performed a simple
simulation to compare with our data: we selected transitions from the neutral
and singly ionized metals Fe, Ti, and Cr, using the Kurucz’s line list (Kurucz
& Bell 1995), and computed their opacity assuming local thermodynamic
equilibrium (LTE) conditions, an excitation temperature Tex=104 K, and solar
abundances. No thermal broadening was assumed and a uniform dynamical
broadening of 8 km/s was applied to all, in agreement with the observed lines
width. Each line was weighted only by its energy level and degeneracy, and the
elemental abundance. By including only 3 elements and assuming LTE conditions,
when in fact the we are probably observing the relics of a wind (see also
Sect. 5), we are far from model fitting the data. Instead, we demonstrate that
the LVC is consistent with a curtain, i.e. passive gas absorbing incident
radiation. In particular, Fig. 5 shows that even such a crude simulation is
capable of reproducing the observed pattern of absorptions in the data.

The curtain and its persistence suggest that the absorbing gas is not part of
the 2008 mass loss but, rather, circumbinary material. This picture is also
supported by the subsequent evolution of the whole LVC. Although the detailed
evolution of the profiles appears different for different ions, the emission
component generally increases over the continuum while the absorption
decreases and appears either in front of the continuum or in front of the
emission components (including the HVC), or both. This evolution is consistent
with a change in the emissivity. There are no signatures of dynamical effects
such as expanding gas, whose change in opacity would produce recombination
waves, or shocking gas clouds, whose interaction would cause substantial line
profile changes. Hence, we propose that the gas responsible for the LVC is
material produced in past mass loss episodes, and that the increase in
emissivity in this material results from a light pulse (the 2008 outburst) and
its time delayed fluorescence and diffusion across the medium. The
circumbinary gas is ionization bounded, at least in the time interval spanned
by the six UVES epochs (as shown by the persistent absorptions), and probably
beyond. Degrading the UVES spectra to the XShooter resolution and inspecting
the Na I D line profile, it appears that partially neutral circumbinary
material is also present at the time of the first two XShooter observations
(see Fig. 6). This further indicates that such a circumbinary gas does not
interact dynamically with the low velocity components at least until day 345.
Moreover, since the XShooter epoch 9 spectrum shows no new transitions and is
a fainter replica of the epoch 7 and 8 spectra (see Section 4), we conclude
that there is no dynamical interaction between the LVC and the HVC until at
least day 392.

Figure 6: Na I D emission lines evolution across time (time increase from top
to bottom). The Black color is for the original spectra (UVES or XShooter),
while the light blue is for the UVES spectra degraded to the XShooter
resolution. Comparing epochs 5 and 6 (UVES) with epochs 7 and 8, it is evident
that what appear as a ”red wing plateau” in low resolution spectra, could well
be a signature of the non-resolved persistent absorption. Velocities are in
the heliocentric rest frame. Y-axis units are erg/cm2/s/Å.  Figure 7:
Comparison between the UVES epoch 6 Na I D and the forbidden transitions ([O
I] and [Ca II]) profiles. The color code for the identification of each line
is in the figure panel itself. The spectrum in the figure has been corrected
for reddening and telluric absorptions. Velocities are in the heliocentric
rest frame.

A final, remarkable aspect of the observed line profiles is that the HVC
completely lacks emission with positive velocities or better symmetric with
respect to the systemic velocity (which we estimate to be $\sim-$120 km/s in
Section 3.3). Since we regard it unlikely that there was an asymmetric mass
loss from just one stellar hemisphere, we suggest that the receding component
of the outflow is obscured or rather self-obscured by dust within the ejecta.
Shore et al. (2018 and references therein) showed how the dust can suppress
the receding part of an emission line profile in an axis-symmetric ejecta
inclined to the line of sight. In this case, in order to hide the receding
portion of the emission line, the dust must be well mixed in the ejecta and
have a large opacity. Using the He I $\lambda$5876 emission line in the epoch
5 spectrum, since it is HVC dominated (i.e. unbiased by the LVC), and the
intensity ratio between the line peak and the point in the red wing which is
symmetric around the peak with respect the systemic velocity, we find almost 4
mag (3.8 mag, precisely, and $\tau_{d}$=3.5) of dust extinction across the
line. Using the He I $\lambda$7065 emission to perform the same estimate we
find 3 mag of dust extinction (or $\tau_{d}$=2.7) between the high velocity
peak emission and its corresponding symmetric point with respect to the
systemic velocity. Assuming silicate dust, as suggested by the 2010
observations of Nicholls et al. (2013), our observed $\tau_{d}$ ratio is
consistent with the grains being $\gtrsim$0.1 $\mu$m already 35 days after
discovery (using Draine & Lee, 1984, absorption efficiency plot in their
Figure 5). Repeating the exercise with the epoch 6 spectrum, we obtain
$\tau_{d}(5875)$=3.6 and $\tau_{d}(7065)$=3.2 for a ratio of 1.13 which may
indicate grain growth to sizes of order of $\sim$0.3 $\mu$m. We could not
repeat the measurements and calculation in the first four spectra, since,
although the He I lines are present, the $\lambda$7065 red wing is biased by
partial blending with adjacent emission lines, while the He I $\lambda$5876
shows also the contribution of the LVC emission in its red wing.

The global picture, therefore, is that the HVC is mass lost during the 2008
outburst, and the LVC is surrounding circumbinary material expelled sometime
in the past. Their exact geometry is not constrained at this stage. However,
there are two possibilities: either the two structures are very distant to
each other so that the HVC does not reach the low velocity gas within the time
of the spectroscopic observations333For an expansion velocity of the HVC of
$\sim$400 km/s (see Section 5) the very lower limit for the LVC’s distance is
400$\times$392$\times$86400 km $\simeq$90 AU; or they have orthogonal
geometries so that they cannot interact regardless of their time of expulsion
and relative distances. We will return to this in Section 5.

###  3.3 The systemic velocity

Having established the circumbinary nature of the LVC, some quantitative
information can be derived if it has a minimal axisymmetric distribution
around the central source (e.g. a spherical shell, an equatorial disk, etc.).
Indeed, for optically thin lines, such as [O I]$\lambda\lambda$6300,6364 and
[Ca II]$\lambda\lambda$7291,7323 that are fully in emission in epoch 5 and 6,
the whole emitting volume is visible. Thus, the line centroids match the
volume geometrical center where the merger is sited. Fig. 7 compares the
profiles of the [O I], [Ca II] and Na I transitions in the epoch 6 spectrum.
It shows that the blue wing of the emission, i.e. the outermost layer of the
low velocity gas, might still be unexcited (undetected), implying a bias
toward lower values if the line centroid is used. To overcome such a bias we
measured the [O I]$\lambda$5577 line, that is not a ground state transition
and is free from absorption components, and compared its flux centroid with
that of the peak of the other forbidden transitions as well as the bisector
between the blue wing of the early epoch absorption and the red wing of the
late epoch emission. The resulting radial velocity of the system, $\gamma$, is
between $-$120 and $-$125 km/s. Hence, the projected terminal velocity of -180
km/s observed in the Na I D persistent absorption reduces to $\sim-$60 km/s.
This value of $\gamma$ revises the $-$80 km/s estimated by Mason et al.
(2010), and the $-$175 km/s reported by Kaminski et al. (2015), who, however,
concluded that the derived velocity was unlikely representative of the system
radial velocity. Our $\sim-$120 km/s heliocentric systemic velocity is in good
agreement with the v${}_{LSR}=-$100 km/s, v${}_{helio}\sim-$110 km/s, measured
by Kaminski et al. (2018).

###  3.4 The column density of the low velocity component

As mentioned in Section 3.2, the low velocity persistent absorption is also
evident in the forbidden transitions [O I]$\lambda\lambda$6300,6364, [Ca
I]$\lambda$6527, and [Ca II]$\lambda\lambda$7291,7323. The presence of an
absorption component in forbidden transitions indicates a very large column
density of the intervening gas. After removing the telluric absorptions around
the [Ca II] and [O I] lines using molecfit (Smette et al 2015, Kausch et al.
2015) and IRAF/telluric, we followed Savage & Sembach (1991) to compute the O
I, Ca I and Ca II column density. Estimated column density for each absorption
and each epoch are listed in Table 2. Uncertainties are of the order of
$few\times 10^{14}$ and $few\times 10^{16}$ for the [Ca I] and [Ca II]
transitions, respectively, and $few\times 10^{19-20}$ for the [O I]. The
larger uncertainties in the [O I] are explained by the ambiguity on the local
continuum level position (given the weak emission component) and the suspected
presence of blending features (e.g. [Ni I], see Allende Prieto et al. 2001).
Using the equivalent width relation for optically thin gas (e.g. Spitzer
1968), we obtained similar values for the Ca0,+ and the O0 column density.
Adopting Asplund et al. (2021) solar abundances for Ca and assuming that
calcium is either neutral or singly ionized, we derive N${}_{H}\sim$2.5$\times
10^{23}$ cm-2. Similarly, we derive N${}_{H}\sim$5.3$\times 10^{23}$ cm-2 from
the average of the NO measurements. A factor of 2 uncertainty is compatible
with the method and the error in the measurements.

Table 2: The Ca, Ca+ and O column density in cm-2. Epoch/Age |  [Ca II]$\lambda$7323 |  [Ca II]$\lambda$7291 |  [O I]$\lambda$6364 |  [O I]$\lambda$6300 |  [Ca I]$\lambda$6573
---|---|---|---|---|---
1/10 | 5.0e17 | 5.0e17 | 3.4e20 | 3.8e20 | 8.0e13
2/15 | 5.5e17 | 4.0e17 | 2.0e20 | 2.2e20 | 1.3e14
3/17 | 5.5e17 | 3.7e17 | 1.9e20 | 2.3e20 | 1.4e14
4/25 | 5.5e17 | 5.0e17 | 1.6e20 | 3.4e20 | 2.4e14

###  3.5 The light pulse energy distribution

The low velocity gas absorbs in the neutral and singly ionized metals such as
Fe, Sc, Ti, V, Cr, Sr, Ba and most importantly Ca, the neutral (alkali metals)
Na and K, and neutral (non metals) O and H. All these transitions are either
ground transitions or produced by pumping mechanisms, or fluorescence.

The H Balmer line absorptions indicate a very opaque Ly$\alpha$ capable of
substantially populating the $n=2$ energy level. The emission part, instead,
traces recombination in addition to decay from the excited levels, but such an
emission is minimal in the LVC. Fe, Sc, etc, because of their low ionization
potential energies are mainly singly ionized. They belong to the so called
iron-curtain elements whose resonance transitions in the UV range pump optical
emissions with metastable low energy levels (Shore & Aufdenberg 1993). The
emission component of these transitions is relatively strong over the
continuum just because of such a pumping mechanism. The absorption is also
quite strong because of the overpopulation of the metastable low energy
levels. As already noted by Mason et al. (2010), a pumping mechanism also
occurs in the case of the Ca II transitions and explains the relatively strong
emission component of the [Ca II] doublet compared to [O I]. The Ca II H&K
transitions are extremely opaque (saturated) and overpopulate their upper
level to the point that a detectable number of de-excitation occurs for the
alternative channel of the NIR Ca II triplet and the resonance [Ca II]
doublet.

In contrast, the evolution of the oxygen transitions illustrates how the
fluorescence strengthens with time. The [O I] doublet upper level is populated
either by collisions or by the emission of the intercombination line O I]
$\lambda$1641, pumped directly by the ground transition $\lambda$1302 or
indirectly by the fluorescent match of the Lyβ and the O I $\lambda$1025 lines
(see Shore & Wahlgren 2010). The minimal emission component of the [O I]
doublet in the first four UVES epochs, together with the lack of the [O
I]$\lambda$5577 line444The [O I]$\lambda$5577 emission appears in the epoch 5
UVES spectrum and persists across the XShooter spectra. and the relative
weakness of the emission component of the triplets at $\sim$8446 Å
(fluorescent) and $\sim$7774 Å (recombination), suggest that fluorescence and
recombination are only weakly occurring. They become important in the epoch 5
and 6 spectra, once the incoming radiation has progressed through and affected
a substantial fraction of the circumbinary material. That some recombination
is, however, occurring since at least epoch 1, is shown by the very weak LVC
emission at the He I $\lambda$5876 in the first four UVES epochs555The line
strengthens until epoch 3 and declines afterward. In epoch 5 and 6 if present,
the LVC He I emission is at the level of or below the red wing of the HVC
emission.. Its weakness, together with the lack of other LVC He I emissions
indicate that very little He has been ionized and is recombining. This implies
either that the incident He-ionizing radiation is minimal, perhaps diluted by
a distance effect, or that the hydrogen is extremely optically thick and
absorbs the vast majority of the photons with high enough energy to ionize
helium. Either way, the set of observed lines and their profile suggest that
the radiation reaching the low velocity gas has a high energy tail up to
$\gtrsim$24 eV, but the bulk of its distribution peaks well below 20 eV. This
is further supported by the lack of C I and N I LVC emissions that, if present
would indicate recombination (IP$\sim$11 and 15 eV for the C and N,
respectively) and of [N II] absorption that suggest there is no ionized
nitrogen (the [N II] $\log(gf)$ is close to that of [O I], while the N
abundance is much greater than that of Ca), although for the nitrogen also
apply the same reasoning than for the helium, i.e. it is little to no ionized
because of the extreme opacity of the hydrogen.

We also note that the He I $\lambda$5876 reaches maximum strength in the epoch
3 spectrum and that in epoch 4 the line flux has already decreased by a factor
$\exp(-1)$ indicating a recombination time of $\sim$ 8 days. This, using
Hummer & Storey (1998) He I recombination coefficients, indicates electron
densities of the order of n${}_{e}\sim$7$\times$107 cm-3.

The HVC is more ionized since it displays both He I and He II emissions. The
He I lines are already present in the epoch 1 spectrum and strengthen with
time until epoch 3, declining afterward. They eventually disappear in the
XShooter spectra sequence taken about one year after outburst. The He II
$\lambda$4686 emission, instead, is weak, but present, only in the epochs 5
and 6 spectra. This underlines once more the different conditions in the high
velocity gas with respect to the LVC. The former is initially opaque with no
emission except for the H Balmer series and the He I lines. In particular, it
is too opaque to produce He+2. Only once the HVC gas has expanded sufficiently
does the He II emission line appear. Also, the velocity profile of the
$\lambda$4686 does not match that of the He I lines but is weighted toward the
low velocity range of the HVC (see Fig.8). In contrast to an explosive ejecta
where the higher ionization emissions appear (first) in the ejecta outskirts
where the density is lower, the ionization structure here is similar to that
of the Stromgren sphere and the density gradient, if any, is not as steep.
Hence, the HVC is also irradiated by an underlying source and either receives
more flux or see a ”harder” source than the LVC gas. In the first case the low
and the high velocity component are far apart, in the second case the low
velocity component does not directly see the central source. We have no way to
tell the two scenarios apart.

Figure 8: Comparison of the He I (black) and He II (light blue) line profiles
in the epoch 6 spectrum. The He II line has been scaled by a factor 4. The
vertical axis is flux in unit of erg/cm2/s/Å. Velocities are in the
heliocentric rest frame. See text for details.

The detection of the He II $\lambda$4686 line increases the upper energy
cutoff of the central source to $\sim$55 eV, in the extreme UV. Such high
energies are sufficient to produce also C+3, N+3 and O+3 but we do not see C
III, N III or O III emission lines (e.g. $\lambda$4650, $\lambda$4640, and
$\lambda$3047 and $\lambda$3759, respectively), as expected from
recombination. Their nondetection could simply result from a low emission
measure - after all the He II $\lambda$4686 is less than 2% of the He I
$\lambda$5876\. This suggests that the central source spectral energy
distribution extends into the EUV range, but peaks at $\sim$20 eV or less.
Either way, the disappearance of the He I and He II emission lines in the
XShooter spectra suggests that the central ionizing source has turned off
sometime between day +48 and day +244.

Figure 9: U to V band of the XShooter spectra. Black, light blue, and gray
lines are for epoch 7, 8 and 9, respectively. Y-axis units are erg/cm2/s/Å; X
axis units are Å. The epoch 9 spectrum was cut below $\sim$3800 Å, because no
signal was evident there in the 2D frame. The epoch 8 and 9 spectra are almost
identical in the continuum (where present), with the epoch 9 spectrum having
stronger molecular bands in emission. Bracketed numbers refer to the Moore’s
multiplet identification (Moore 1945).

##  4 Analysis of the XShooter data set

The XShooter data set consists of three epochs about one year after the
outburst, at days +244, +345 and +392. By this time the spectra appearance has
changed dramatically. The continuum in the UVB arm (3000-5500 Å) gradually
fades below the detection threshold of the exposures. It is still detected in
the epoch 7 spectrum but consistent with zero in the subsequent two epochs.
Hence, we cannot verify whether the curtain produced by the low ionization
potential metals that we observed in the UVES spectra is still present. In
general, however, the lower spectral resolution of XShooter compared to UVES
prevents firm conclusions about what, if any, narrow absorption that might
still be present and absorbing against the continuum. Our analysis is,
therefore, limited to the emission line spectrum. But, as discussed in Sect.
3.2 and shown in Fig. 6, it is quite possible that the LVC is still present in
the epoch 7 and 8 spectra. In epoch 9 the Na I D emission has almost
disappeared, so the presence the low velocity persistent absorption cannot be
determined.

The emission line spectrum is produced mainly by the high velocity component.
It is difficult to identify all the emission lines with certainty since at any
given wavelength there are a number of candidate low ionization potential
metal identification (van Hoof 2018, NIST database), and we have already seen
numerous blends in the higher resolution UVES spectra. However, we can safely
say that the majority of the observed emission are Fe II and [Fe II]
transitions. The Fe II lines were already present in the UVES spectra with the
high velocity emission first appearing in the epoch 4 spectrum, reaching the
observed maximum 10 days after and subsequently gradually fading or
disappearing. The [Fe II] emissions, absent in the UVES spectra, are observed
in the epoch 7 XShooter spectrum when they have maximum strength. They
slightly weaken in the subsequent epochs. All the [Fe II] emission are ground
or very low lying level transitions (the energy of the lower level is in the
range 0-0.23 eV), with an energy jump of $<$3 eV; while the Fe II lines have
upper energy levels $<$6 eV and, de-exciting, feed the upper levels of the
forbidden transitions. That is, the [Fe II] lines arise both from collisional
excitation and from cascade down from the low levels of the Fe II transitions,
and, therefore, seems consistent with just the dilution of the high velocity
component in the circumbinary space.

Other forbidden transitions that appear in the HVC and strengthen in the
XShooter epochs are Mg I] $\lambda$4571, [O I], [Ca II], and [S II]. The [S
II] is particularly interesting because the sulfur ionization potential
($\sim$10 eV) is lower than for the hydrogen and oxygen, but much higher than
for iron and related metallic species. Hence, it sets a stringent limit on the
hardness of the ionizing source since we do not see any forbidden [O II]. Of
the two [S II] doublets $\lambda\lambda$4068,4076 and
$\lambda\lambda$6716,6730 only $\lambda$4068 is fairly strong and with a good
S/N; the others are weak and obvious only in the epoch 7 and 8 spectra.
Measurements of the doublet at $\sim$6700 Å are compromised by blending with a
molecular band. Nonetheless, using the five level models from Keenan et al.
(1996), our line ratios $\log(F(4068)/F(6716+6730))$
$\log(F(4076)/F(6716+6730))$ are beyond their density range for any adopted
reddening, suggesting electron densities in excess of 105 cm-3. The large
F(4068)/F(4076) ratio, $\sim$2.5 to 4, depending on the epoch, suggests
densities close or above the critical value ($n_{cri}\sim$2-3$\times$106
cm-3)666The expected ratio in the low and high density limit is 2 and $\sim$5,
respectively. The fact that the ratio increases between epoch 7 (2.5) and 8
(4), suggests in addition that there is a density gradient in the high
velocity component, with the emissivity of lowest density portion dropping in
favor of that of the higher density part, although the mid-resolution of
XShooter does not suffice to identify sub-components within the line. The
oxygen F(5577)/F(6300) measured line ratio ($\sim$0.2), compared with
theoretical multilevel balance predictions assuming that the ionization
equilibrium is dominated by electron collisions and charge exchange reactions
and not an external radiation field (e.g. Hamann 1994), also indicates
$n_{e}\gtrapprox$106 cm-3. Finally, solving the equation of the detailed
balance in the three level approximation for the [O I]
$\lambda\lambda$6300,6364, $\lambda$5577, $\lambda$2972 in the high density
limit, we also obtain densities in the range 2$\times$10${}^{6}\leq
n_{e}\leq$4$\times$106 cm-3, with an obvious addition by the low velocity
component in the red wing of the $\lambda$5577 line from epoch 8 on. The
estimated densities increase to 4.6$\times$106 cm-3 adopting E(B-V)=0.84 mag.
Again the apparent increase of the density with time indicate the presence of
a density gradient.

##  5 Discussion and conclusions

In the spirit of Sherlock Holmes ”We balance probabilities and choose the most
likely. It is the scientific use of the imagination.” (Conan Doyle, 1902), the
picture that emerges from the whole analysis is that there are two distinct
line forming regions contributing to the observed spectrum at any time: the
low and the high velocity component.

We identify the HVC as the mass outflow produced with the 2008 outburst. Its
expansion velocities are likely larger than observed because we are viewing
the binary at some intermediate inclination (see below), although the measured
radial velocity (up to $\sim-$400 km/s) is augmented by the systemic radial
velocity ($\gamma\sim-$120 km/s)777For $i\sim$45o and taking into account the
systemic radial velocity, the outflow velocity could be as large as $-$400
km/s.. The extremely peaked, narrow profile that is offset with respect to
$\gamma$, of the HVC once in emission, is consistent with sudden, time limited
mass loss in a rapidly accelerating, possibly somewhat collimated, wind (for
instance, as modeled for the outbursting symbiotic MWC560 by Shore et al.,
1994); the peak of the emission line corresponding to where the velocity
approaches saturation, its width resulting from the velocity gradient. The
missing ”red side” indicates the presence of dust obscuring the receding part
of the outflow in an intermediate inclination system: dust in bipolar ejecta
viewed edge-on cannot asymmetrically affect the line profile; conversely, the
inclination does not play a role in the observed asymmetry only in case of
spherical ejecta. Ground truth about the dust effect comes from the inspection
of the sub-mm SiO and CO line profiles in Kaminski et al. (2018, their figure
8 in particular), which are immune to extinction and show the whole emission
region with its receding part. Finally, part of the emission just redward of
about $-$200 km/s, is also blocked by the LVC and argues some geometrical
thickness of the latter. The HVC shows a ionization gradient with the highest
energy transitions in its inner/low velocity layers, and the low energy
emission more in its outskirt. This ionization structure is similar to the
nested ionized regions in a planetary nebula and is also consistent with a
density distribution produced by rapidly accelerating wind of relatively short
duration; it also demonstrates the existence of a central ionizing source.
Initially very opaque, the HVC turns transparent only after the central source
has ”turned off”, so that nothing but forbidden transitions from low
ionization potential energy metals can develop. After about one year from the
outburst the HVC shows a density of several 106 cm-3.

The LVC matches circumbinary gas with a substantial column density and which
absorbs and is excited in its inner layers by the outburst (i.e. the gas is
ionization bounded). Its fluorescence, observed particularly well in the [O I]
and [Ca II] lines, is similar to that of the SN 1987A ring (Dwek & Felten
1992). The time required to reach the fluorescence maximum depends on the
distance of the gas from the central source, the duration of the central
source pulse, and the physical conditions of the gas and its local emissivity.
This implies that we cannot disentangle the light delay time (to reach the
gas) and the diffusion time (within the gas), but the maximum flux is observed
once the light from the whole emitting volume has reached the observer. This
happens in epoch 5 indicating a representative time scale for the combined
delay+diffusion time of $\sim$35 days. Without the delay time we cannot
precisely determine the distance of the LVC from the central source or when it
left the system, but an upper limit for its distance can be estimated assuming
no diffusion time, so that 35/2 light days, i.e. $\sim$3000 AU, is the LVC
standoff distance. We know, however, that the diffusion time is not negligible
since the [O I]$\lambda\lambda$6300,6364 absorptions were already present in
the first spectrum indicating that the light pulse had already reached the LVC
in less than 10 days. Without knowing the delay time it is not possible to
accurately compute the LVC density, but the mere presence of the [O I]
emission lines and the $\lambda$5577 in particular, argues for values not much
larger than the critical density, $n_{e}\sim$108 cm-3. This value is
consistent with that derived in Section 3.5 from the recombination timescale
of the He I, implying that the LVC is stationary in the time interval of our
observations and further supporting the scenario in which the LVC is distinct
from the transient (i.e. central object plus ”ejecta”) and unaffected by it
other than by fluorescence. The large column density (3-5$\times$1023 cm-2)
suggest a substantial mass loss from past times. When combined with the O I
$\lambda$5577 critical density ($\sim$4$\times 10^{23}$/108) the radial
extension of the LVC is as large as $\sim$260 AU, suggesting long formation
time scales. Taken together, this is in line with McCollum et al. (2014)
analysis of the Spitzer preoutburst data whose color matches those of ”highly
evolved stars with large amount of chemically complex mass loss, such as AGB
stars”. But the LVC is even more complex since it shows, in addition, (see Na
I profile in Fig.7) non-saturated narrow absorption components at lower
velocity (in the range $-140$,$-100$ km/s) which might arise from mass lost
from the Lagrangian point L2 and the spiral structure predicted by theoretical
models (e.g. Webbink 1976, Pejcha 2016a, 2016b, 2017, MacLeod & Loeb 2020).
The fact that they are observed only in the Na I transitions indicate a
significantly lower column density than in the main LVC absorption.

There are no signs of dynamical interaction between the two components. All
the observed emission processes are consistent with photon processes and not
with collisions: from the absorbed blue ($<$4000 Å) continuum suggestive of
Rayleigh scattering by the HVC, to the HVC ionization structure, to the
fluoresclent emission lines in the LVC. This lack of dynamical interaction
signatures suggests that the LVC and the HVC are orthogonal to each other. The
alternate scenario of two distant concentric shells seems less likely, for,
otherwise, we should not be able to observe the HVC while the low velocity gas
is opaque (which it always is), nor would we see its excited inner layers –
unless it is extremely fragmented and patchy. The line profiles do not
indicate this. The LVC, however, is geometrical thick, since it intercepts
part of the light from the HVC. We, therefore, propose that the HVC is a
somewhat bipolar rapidly accelerating wind, and the LVC has a torus-like
geometry. Although V1309 Sco does not yet appear convincingly elongated in
submm ALMA images (and we suggest this might be due to its distance in
addition to its age), bipolar geometries have been inferred for similar
objects from submm observations (Kaminski et al 2018).

In summary, we have observationally obtained a description of the structure
and geometry of V1309 Sco during the 2008 outburst, and inferred approximate
densities for the emitting volumes and qualitative dynamics. We also provided
evidence for a high energy cutoff of the emission from the central source, but
we cannot further constrain its bolometric luminosity or its precise spectral
energy distribution (SED). Nonetheless, whatever the energy generation
mechanism(s) was and however it varied, the central source certainly turned
off on time scales of months. It emitted in the UV range but spectroscopic
observations were not taken in this energy band during V1309 Sco event. This
spectral range is especially important during the rise time, since it would
have provided a stronger constraint of the SED, the start time of the event,
and when the ejecta turned UV optically thick. In light of the column density
we inferred for the V1309 Sco circumbinary material, the curtain would have
sufficient opacity to completely redistribute the UV flux to the optical if
the source was completely covered. Hence, prompt FUV and NUV observations
should be attempted for the next candidate merger since the relative
emissivity of the UV and optical flux depends on the source coverage, if
independent reddening estimates can be obtained. The same argument goes for
the IR and dust redistribution and is independent of the distance. Also, V1309
Sco optical light curve is photometric redistribution by the mass expelled
with the outburst and by the circumbinary gas as they recombine. Since the
reprocessing material is not spherical, the observed optical light is only a
fraction of the central source bolometric emission.

V1309 Sco, might have been a peculiar event, but it is the only certain
noncompact binary star merger we know. It is therefore important to consider
all of the above constraints when producing theoretical models of noncompact
stars mergers. Stellar mergers resulting from a common envelope event have
been modeled by numerous authors both as general case (see Ivanova et al. 2013
for a comprehensive review; see also Pejcha et 2016a, 2016b, Metzger & Pejcha
2017, MacLeod & Loeb 2020 for more recent developments) and simulated
specifically for V1309 Sco (e.g. Nandez et al. 2014, Pejcha et al. 2017 and
references therein). Direct comparison between those models and our results is
not straightforward since theoretical codes are based on a number of
simplifications that might not apply to V1309 Sco specific case (e.g. the
point mass plunging in the extended envelope of a giant star), and seek to
reproduce photometric behaviors which might be incomplete and are notoriously
ambiguous. Still, the following points are readily noticeable.

1) All the models predict a substantial mass loss during the phases preceding
the dynamical merger and resulting in the presence of circumbinary material
with a toroidal or disk-like geometry. While we confirm the presence of such a
circumbinary material, its extension and formation timescale are not quite in
agreement with the theoretical predictions. Our circumbinary material is well
detached from the central binary as shown by the time delay between the
appearance of the blue shifted LVC absorption and its emission counterpart.
The simulated medium is much more compact, relatively close or even connected
to the binary with scale of tens of giant-star radii or orbital separations,
and piled up on time scales of years at the very most. Observationally, for a
standoff distance of $\leq$3000 AU and a drift velocity of $\sim$60 km/s, the
gas left the system hundreds of years before the outburst. The independent
estimate of the gas extension along the orbital plane of $\simeq$260 AU, also
argue for long formation time scales and is orders of magnitude larger than
that predicted by the simulations (e.g. compare with the $\sim$60 binary
separation, $a$, or stellar donor radii, $R_{1}$, in Figures 9 of Pejcha et
al. 2017 and 5 of MacLeod & Loeb 2020). The observed column density is also
order of magnitude smaller than predicted by Pejcha et al. (2017) simulations.

2) In the attempt to reproduce the observed light curves of V1309 Sco and
alike transients, Metzger & Pejcha (2017) proposed shocks between the slowly
moving outflow from L2 and the fast ejecta launched during the dynamical phase
of the merger. In particular, second light curve peaks or plateau result
depending on the shock energetic capable (or less) to keep hydrogen ionized.
Our non detection of any signature of dynamical interaction between, and
within, the HVC and the LVC in the time span of our observations (392 days)
dismiss the proposed scenario for V1309 Sco. It will be important to monitor
through sufficient high resolution spectroscopy future candidate stellar
mergers to collect kinematic and dynamical information of the expanding
gas(es) during their major photometric phases. Here, we recall the light echo
effect produced by circumstellar dust that is off the line of sight and that
might have same, bluer or redder color depending on the differential
extinction between the scattering dust and the transients. This is fairly well
known for supernovae (e.g. Patat 2005) and has also been modeled for classical
novae (Shore 2013). Ultimately only (spectro)polarimetry and line profile
evolution studies can serve as a diagnostic for the mechanisms behind a given
light curve. Imagery and bolometry alone are insufficient.

We conclude by noting the importance of considering the observed spectral
signatures when defining the parameter space to be explored in simulations.
Models and simulations should predict spectroscopic observables. Ideally,
detailed treatment of the relevant radiative processes should be included in
the common envelope studies to predict the observed spectrum and its time
evolution.

###### Acknowledgements.

Based on observations collected at the European Organisation for Astronomical
Research in the Southern Hemisphere under ESO programmes 281.D-5055(A) and
60.A-9445(A). This research has made use of NASA’s Astrophysics Data System
Bibliographic Services. This work also made extensive use of the NIST atomic
data base and online tool and the van Hoof’ Atomic line list v3.00b4 online
tool. It also used the Partial Grotrian diagram of astrophysical interest by
Moore & Merill (1968, NSRDS - NBS 23). We thank the referee, Tomasz Kaminski,
for comments and suggestions that improved the paper. EM also thanks the
Lodovico Coccato and Christian Hummel from ESO USD and DMO for their help with
molecfit.

## References

  * (1)
  * prieto (2001) Allende Prieto, C.; Lambert, D.L. ; Asplund, M., 2001, ApJ, 556L, 63
  * Asplund (2021) Asplund, M.; Amarsi, A.M.; Grevesse, N., 2021, A&A, 653, 141
  * Dame (2001) Dame, T.M.; Hartmann, D.; Thaddeus, P., 2001, ApJ, 547, 792
  * doyle (1902) Conan Doyle, A., 1902, in The Hound of the Baskervilles, ch. 4 any edition.
  * Draine+ (1984) Draine, B. T. & Lee, H. M., 1984, ApJ, 285, 89
  * dwek (1992) Dewk, E.; Felten, J.E., 1992, ApJ, 387, 551
  * friedman (2011) Friedman, S.D.; York, D.G.; McCall, B.J.; et al. 2011, ApJ, 727, 33
  * hamann (1994) Hamann, F., 1994, ApJSS, 93, 485
  * hummer and storey (1998) Hummer, D.G.; Storey, P.J., 1998, MNRAS, 297, 1073
  * ivanova (2013) Ivanova, N.; Justham, S.; Chen, X.; et al., 2013, A&ARv, 21, 59
  * Jenniskens (1994) Jenniskens, P.; Desert, F.X., 1994, A&ASS, 106, 39
  * kalberla (2005) Kalberla, P.M.W., Burton, W. B., & Hartmann, D., et al., 2005, A&A, 440, 775
  * kaminski (2015) Kaminski, T.; Mason, E.; Tylenda, R.; et al., 2015, A&A 580, A34
  * tomas (2018) Kamiński, T.; Steffen, W.; Tylenda, R., et al., 2018, A&A, 617, 129
  * kausch (2015) Kausch, W.; Noll, S.; Smette, A.; et al., 2015, A&A 576, A78
  * keenan (1996) Keenan, F. P.; Aller, L.H.; Bell, K.L.; et al., 1996, MNRAS, 281, 1073
  * kuruz (1988) Kurucz, R.L.; Bell, B., 1995 Atomic Line Data, Kurucz CD-ROM No. 23. Cambridge, Mass.
  * macleod (2020) MacLeod, M.; Loeb, A., 2020, ApJ, 895, 29
  * Mason (2010) Mason, E.; Diaz, M.; Williams, R. E. ; et al., 2010, A&A, 516, 108
  * mccollum (2014) McCollum, B.; Laine, S.; Väisänen, P.; et al., 2014, AJ, 147, 11
  * metzger (2017) Metzger, B.D.; Pejcha, O., MNRAS, 471, 3200
  * moore (1945) Moore, C.E., 1945, in A multiplet table of astrophysical interest, national standard reference data system, NBS 40
  * munari (1997) Munari, U., & Zwitter, T., 1997, A&A, 318, 269
  * iauc2 (2008) Naito, H.; Fujii, M., 2008, IAUC 8972, 2
  * iauc1 (2008) Nakano, S. ; Nishiyama, K. ; Kabashima, F. ; et al., 2008, IAUC 8972, 1
  * nandez (2014) Nandez, J.L.A.; Ivanova, N.; Lombardi, J.C. JR, 2014, ApJ, 786, 39
  * nicholls (2013) Nicholls, C. P.; Melis, C.; Soszynski, I., et al., 2013, MNRAS, 431, L33
  * patat (2005) Patat, F., 2005, MNRAS, 357, 1161
  * pejcha (2016a) Pejcha, O.; Metzger, B.D.; Tomida, K., 2016a, MNRAS, 455, 4351
  * pejcha (2016b) Pejcha, O.; Metzger, B.D.; Tomida, K., 2016b, MNRAS, 461, 2527
  * pejcha (2017) Pejcha, O.; Metzger, B.D.; Tyles, J.G.; et al., 2017, ApJ, 850, 59
  * savage (1991) Savage, B.D.; Semabch, K.R., 1991, ApJ, 379, 245
  * shore (1993) Shore, S.N.; Aufdenberg, J.P.; 1993, ApJ, 416, 355
  * shore (1994) Shore, S.N.; Aufdenberg, J.P.; Michalitsianos, A.G., 1994, ApJ, 108, 671
  * shore (2010) Shore, S.N.; Wahlgren, G.M., 2010, A&A, 515A, 108
  * shore (2011) Shore, S. N.; Augusteijn, T.; Ederoclite, A., et al., 2011, A&A, 533L, 8
  * shore (2012) Shore, S.N., 2012, BASI, 40, 185
  * shore (2013) Shore, S.N., 2013, A&A, 559, L7
  * shore (2018) Shore, S.N.; Kuin, N.P.; Mason, E., et al., 2018, A&A, 619A, 104
  * smette (2015) Smette, A.; Sana, H., Noll, S.; et al., 2015, A&A 576, A77
  * spitzer (1968) Spitzer, L., 1968, Diffuse matter in space, New York: Interscience Publication
  * Tylenda (2011) Tylenda, R.; Hajduk, M.; Kaminski, T.; et al., 2011, A&A 528, A114
  * peter (2018) van Hoof, P.A.M., 2018, Galaxies, 6, 63
  * webbink (1976) Webbink, R.F., 1976, ApJ, 209, 829
  * (46)


# VC Theoretical Explanation of Double Descent

Eng Hock Lee
Department of Electrical and Computer Engineering
University of Minnesota, Twin Cities
Minneapolis, MN 55455
<EMAIL_ADDRESS>
&Vladimir Cherkassky
Department of Electrical and Computer Engineering
University of Minnesota, Twin Cities
Minneapolis, MN 55455
<EMAIL_ADDRESS>

###### Abstract

There has been growing interest in generalization performance of large
multilayer neural networks that can be trained to achieve zero training error,
while generalizing well on test data. This regime is known as ‘second descent’
and it appears to contradict conventional view that optimal model complexity
should reflect optimal balance between underfitting and overfitting, aka the
bias-variance trade-off. This paper presents VC-theoretical analysis of double
descent and shows that it can be fully explained by classical VC
generalization bounds. We illustrate an application of analytic VC-bounds for
modeling double descent for classification problems, using empirical results
for several learning methods, such as SVM, Least Squares, and Multilayer
Perceptron classifiers. In addition, we discuss several possible reasons for
misinterpretation of VC-theoretical results in the machine learning community.

##  1 Introduction

There have been many recent successful applications of Deep Learning (DL).
However, at present time, various DL methods are driven mainly by heuristic
improvements, while theoretical and conceptual understanding of this
technology remains limited. For example, large neural networks can be trained
to fit available data (achieving zero training error) and still achieve good
generalization for test data. This contradicts conventional statistical wisdom
that overfitting leads to poor generalization. This phenomenon has been
systematically described by Belkin, et al. [1] who introduced the appropriate
terminology (‘double descent’) and pointed out the difference between the
classical regime (first descent) and the modern one (second descent). This
disagreement between the classical statistical view and modern machine
learning practice provides motivation for new theoretical explanations of the
generalization ability of DL networks and other over-parameterized estimators.
Several different explanations include: special properties of multilayer
network parameterization [2], choosing proper inductive bias during second
descent [1], effect of Stochastic Gradient Descent (SGD) training [3, 4, 5],
the effect of various heuristics (used for training) on generalization [6],
and the effect of margin on generalization [7]. The current consensus view on
the ‘generalization paradox’ in DL networks is summarized below:

  * –

Existing indices for model complexity (or capacity), such as VC-dimension,
cannot explain generalization performance of DL networks.

  * –

‘Classical’ theories developed in ML and statistics cannot explain
generalization performance of DL networks. In particular, classical VC
generalization bounds cannot be used to explain double descent. Specifically,
the ability of large DL networks to achieve zero training error (during second
descent mode) effectively ‘rules out all of the VC-dimension arguments as a
possible explanation for the generalization performance of state-of-the-art
neural networks’ [3].

This paper demonstrates that these assertions are incorrect, and that
classical VC-theoretical results can fully explain generalization performance
of DL networks, including ‘double descent’, for classification problems. In
particular, we show that proper application of VC-bounds using correct
estimates of VC-dimension provides accurate modeling of double descent curves,
for various classifiers trained using stochastic gradient descent (SGD), least
squares loss and standard SVM loss. The proposed VC-theoretical explanation
provides many additional insights on generalization performance during first
descent vs. second descent, and on the effect of statistical properties of the
data on the shape of double descent curves.

Next, we briefly review VC-theoretical concepts and results necessary for
understanding generalization performance of all learning methods based on
minimization of training error [8, 9, 10, 11]:

  1. 1.

Finite VC dimension provides necessary and sufficient conditions for good
generalization.

  2. 2.

VC theory provides analytic bounds on (unknown) test error, as a function of
training error, VC dimension and the number of training samples.

Clearly, these VC-theoretic results contradict an existing consensus view that
VC-theory cannot account for generalization performance of large DL networks.
This disagreement results from misinterpretation of basic VC-theoretical
concepts in DL research. These are a few examples of such misunderstanding:

  * –

A common view that VC-dimension grows with the number of parameters (weights),
and therefore, ‘traditional measures of model complexity struggle to explain
the generalization ability of large artificial neural networks’ [3]. In fact,
it is well known that VC-dimension can be equal, or larger, or smaller, than
the number of parameters [8, 11].

  * –

Another common view is that ‘VC dimension depends only on the model family and
data distribution, and not on the training procedure used to find models’
[12]. In fact, VC dimension does not depend on data distribution [8, 11, 13].
Furthermore, VC dimension certainly depends on SGD algorithm [8, 11].

For classification problems, VC theory provides analytic generalization bounds
for (unknown) Prediction Risk (or test error), as a function of Empirical Risk
(or training error) and VC-dimension (h) of a set of admissible models, aka
approximating functions. That is, for a given training data set (of size n),
VC generalization bound has the following form [8, 9, 10, 11]:

$R_{tst}\leq
R_{trn}+\frac{\varepsilon}{2}\left(1+\sqrt{1+\frac{4R_{trn}}{\varepsilon}}\right)$
(1) $\text{where
}\varepsilon=\frac{a_{1}}{n}\left(h\left(\ln\left(\frac{a_{2}n}{h}\right)+1\right)-\ln{\frac{\eta}{4}}\right),\eta=\min\left(\frac{4}{\sqrt{n}},1\right)$

This VC bound (1) holds with probability $1-\eta$ (aka confidence level) for
all possible models (functions) including the one minimizing the training
error ($R_{trn}$). The second additive term in (1), called the confidence
interval, aka excess risk, depends on both the empirical risk (training error)
and VC dimension (h). This bound describes the relationship between training
error, test error and VC-dimension, and it is often used for conceptual
understanding of model complexity control, i.e. understanding the effect of
VC-dimension on test error. Application of this bound for accurate modeling of
double descent curves requires:

  * –

Selecting proper values of positive constants $a_{1}$ and $a_{2}$. The worst-
case values $a_{1}=4$ and $a_{2}=2$, provided in VC-theory [8, 9] result in
VC-bounds that are too loose for real-life data sets [11]. For classification
problems, we suggest using the values $a_{1}=a_{2}=1$, used for all empirical
results presented in this paper.

  * –

Analytic estimate of VC-dimension. For many learning methods, including DL,
analytic estimate of VC-dimension are not known. For example, for SGD style
algorithms, the effect of various heuristics (e.g., initialization of weights,
etc.) on VC-dimension is difficult (or impossible?) to quantify analytically.

Note that VC bound (1) provides conceptual explanation of both first and
second descent. That is, first descent corresponds to minimizing this bound
when training error is non-zero [8, 9, 11]. Whereas second descent corresponds
to minimizing this bound when training error is kept at zero, using a set of
models having small VC-dimension. This can be shown by setting the training
error in bound (1) to zero, resulting in the following simplified bound for
test error during the second descent:

$R_{tst}\leq\varepsilon,\qquad\text{where
}\varepsilon=\frac{h}{n}\left(\ln\left(\frac{n}{h}\right)+1\right)$ (2)

More formally, since VC bound (1) depends only on two factors, training error
and VC-dimension, there are two different strategies for minimizing this bound
[11]:

  * –

For a set of functions (models) with fixed VC-dimension, reduce the training
error. This strategy leads to well-known classical bias-variance trade-off aka
‘first descent’;

  * –

For small (fixed) training error, minimize the VC-dimension. This strategy
corresponds to second descent, when training error is zero.

These two strategies correspond to different methods for controlling VC-
dimension, that have been known long before DL. For example, the second
strategy corresponds to margin maximization in SVM. Traditional learning
methods typically implement a single strategy, whereas practitioners in DL
observed the effect of both strategies when varying a single hyper parameter,
such as network size or the number of epochs.

The main technical reason for misapplication of VC-theory, besides
misinterpretation of VC-dimension, is that VC bound (1) remains virtually
unknown in the DL community. That is, all technical arguments suggesting that
VC-theory is unable to explain double descent, are based on analysis of
uniform convergence bounds [1, 3, 14, 15, 16, 17]. In such bounds, the
confidence interval term, aka the excess error, is of the order
$\mathcal{O}\left(\sqrt{h/n}\right)$, i.e. it does not depend on empirical
risk (training error). However, VC-theory also provides more accurate uniform
relative convergence bounds, such as VC-bound (1), presented in [8, 9, 10,
11], where the confidence interval term also depends on the training error.
So, whereas it is true that uniform convergence bounds cannot explain double
descent, it can be fully explained by uniform relative convergence bounds.

Training of DL networks is based on stochastic gradient descent (SGD), which
incorporates several heuristic rules to ensure that the norm of weights
remains small. These rules include: initialization of weights to small random
values and re-normalization during training. Consequently, for large DL
networks, the model complexity is determined by the norm of weights, rather
than the number of weights (parameters). Further, this dependence of VC-
dimension on the training algorithm helps explain why theoretical estimates of
VC-dimension based only on network topology have found little practical use
[11].

##  2 Application of VC Bounds for Modeling Double Descent

This section presents a VC-theoretical explanation of double descent for
classification, for a single-layer network shown in Figure 1. The same network
setting was used for analysis of double descent in recent papers [1, 18, 19,
20]. In this network, a classifier is estimated in two steps:

  * –

First, input vector x is encoded using N nonlinear features. Commonly, random
features (weak features) are used, such as random ReLU or Random Fourier
Features (RFF);

  * –

Second, a linear model is estimated in this N-dimensional feature space.

This simplified setting enables VC theoretical analysis of double descent,
because the analytic estimates of VC-dimension are known. That is, since the
network output is formed as a linear combination of N features, analytic
estimate of VC-dimension for linear hyperplanes
$f(\boldsymbol{z},\boldsymbol{w})=(\boldsymbol{wz})+b$ is known [8, 9, 10]:

$h\leq\min\left(||\boldsymbol{w}||^{2},N\right)+1$ (3)

This bound holds under the assumption that all training samples are enclosed
within a sphere of radius 1, in Z-space. In summary, VC-dimension can be
bounded by the input dimensionality (N), or by the norm of weights. These are
two different mechanisms for controlling VC-dimension.

The double descent phenomenon can be observed for various learning methods
used to estimate weights $\boldsymbol{w}$ for the network structure in Figure
1. Next, we present empirical results showing the application of VC bounds (1)
and (3) for modeling double descent when network weights are estimated using
SVM or Least-Squares (LS) classifiers. For large networks trained using LS,
when N is larger than sample size (n), minimization of squared error is
performed using pseudo-inverse, which finds a solution corresponding to the
minimization of the norm squared $||\boldsymbol{w}||^{2}$.

Figure 1: Single hidden layer network estimating a linear classifier in
nonlinear feature space (Z-space).

In all experimental results in this section, double descent is observed when
the network size (N) is gradually increased. Specifically, according to
analytic bound (3):

  * –

When network size (N) is small, the VC-dimension initially grows linearly with
N. This corresponds to the first descent, or traditional bias-variance trade-
off.

  * –

For overparameterized networks (large N), VC-dimension is controlled by the
norm squared of weights, leading to second descent.

We use two types of random nonlinear features [1, 21], ReLU and RFF. Random
ReLU features are formed as:

$\boldsymbol{Z}_{i}=\max\left(\langle\boldsymbol{v}_{i},\boldsymbol{X}\rangle,0\right),\quad
i=1,...,N$

where random vectors $\boldsymbol{v}_{1},…,\boldsymbol{v}_{N}$ are sampled
uniformly from the range [-1,1]. Random Fourier Features (RFF) are formed as:

$\boldsymbol{Z}_{i}=\exp\left(\sqrt{-1}\langle\boldsymbol{v}_{i},\boldsymbol{X}\rangle\right),\quad
i=1,...,N$

Where random $\boldsymbol{v}_{1},…,\boldsymbol{v}_{N}$ are sampled from
Gaussian distribution with standard deviation $\sigma=0.05$. In all
experiments, input (x) values were pre-scaled to [0, 1] range, for both
training and test data.

Following the nonlinear mapping X $\rightarrow$ Z, all z-values are re-scaled
to [-1, 1] range. Such re-scaling is performed to satisfy the condition for
bound (3), stating that all training samples in Z-space should be enclosed
within a sphere of radius 1.

Training samples (z, y) are used to estimate a decision boundary in Z-space.
Two different methods (LS and SVM classifiers) are used for estimating linear
decision function $f(\boldsymbol{z},\boldsymbol{w})=(\boldsymbol{wz})+b$ from
training data, in order to show double descent curves for two different loss
functions, LS and SVM loss. For SVM modeling, the regularization parameter C
is set to 64 in all experiments. Empirical test error is estimated using an
independent test set.

Figure 2: Application of VC-bounds for MNIST digit 5 vs 8 data set using
random ReLU features. Figure 3: Application of VC-bounds for MNIST digit 5 vs
8 data set using RFF.

Most empirical results in this paper were obtained for MNIST digits adapted
for binary classification (digit 5 vs 8), where digits are grey-scale images
of size 28x28. The training set size n = 800, and test set size is 2,000. We
have also used other data sets for modeling and observed similar results. See
Appendix for additional results.

Empirical results in Figures 2 and 3 show application of VC bounds to modeling
MNIST data. They show:

  * –

Empirical training and test error curves, as a function of N ( the number of
nonlinear features), along with VC-theoretical estimate of test error obtained
via bounds (1) and (3). These curves show that analytic VC-bounds can explain
(and predict) double descent;

  * –

The ‘norm squared’ of estimated linear model, as a function of the number of
features N;

  * –

For SVM, we also show the number of support vectors for trained SVM model.

Modeling results for random ReLU and RFF features are similar, so we only
comment on results in Figure 2:

  * –

For small N, VC-dimension grows linearly with N for SVM method. Empirical
results show that first descent error curves can be explained by VC-bound (1),
because the minimum of VC-bound closely corresponds to the minimum of test
error. This can be clearly seen for LS classifier, and less obvious for SVM.

  * –

For large N, VC-dimension is controlled by the norm squared, according to
bound (3). These results show that second descent can be explained by VC-bound
(1), for both SVM and LS learning methods.

Whereas empirical results for both SVM and LS in Figure 2 are qualitatively
similar, their double descent curves show different values of interpolation
threshold N* (where the training error reaches zero). For SVM, the value N*
$\approx$ 100 is achieved when the number of features equals the number of
support vectors. For LS classifier, the interpolation point N* $\approx$ 800
is achieved when the number of features equals the number of training samples.

The dependence of test error on the norm of weights in large networks has been
known to practitioners, and some limited theoretical explanation is provided
in [1, 22, 23]. For example, [1, 22] suggest that minimum norm provides
inductive bias by favoring models with higher degree of smoothness. However,
these papers do not mention VC bounds that clearly relate the VC-dimension to
the norm of weights, and explain generalization performance for linear
classifiers.

The dependence of interpolation threshold N* on training sample size for LS
classifiers has also been observed in the DL literature. However, in the
absence of sound theoretical framework for double descent, interpretation of
this empirical dependency leads to convoluted explanations. For example,
Nakkiran, et al. [12] investigated the effect of varying the number of
training samples on test error, for a fixed-size DL network. In particular,
they observed two double descent error curves for the same network trained
using smaller and larger size training data, showing that during second
descent, near interpolation threshold, the test error for a network trained
with larger data set is worse than for the same network trained on smaller
data set. This phenomenon was called ‘sample-wise non-monotonicity’, and a new
theory was proposed for explaining regimes where ‘increasing the number of
training samples actually hurts test performance’. However, this phenomenon
has a simple VC-theoretical explanation, as explained next. Note that for LS
classifiers during second descent the shape of the ‘norm squared’ closely
follows the shape of test error, according to VC bound (2), as evident in
Figures 2 and 3. Since for LS classifiers the interpolation threshold is given
by training sample size, there is a region near interpolation threshold, where
VC-dimension for a smaller training size is smaller than for a larger training
size. According to VC-bound (2), in this region of second descent we can
expect a smaller (better) test error for a smaller training size.

VC theoretical framework can also help to understand the effect of statistical
characteristics of training data on generalization curves. Next, we present
empirical results demonstrating the effect on noisy data on the shape of
double descent curves, along with their VC theoretical explanation. For these
experiments, we use a single-layer network trained using SVM and LS classifier
using random ReLU features. We use digits data with corrupted class labels.
The training set size is 800 (400 per class), and the test set size is 2,000.
Figure 4 shows the effect of noise level on the shape of double descent
curves, for SVM and LS classifiers. Results for both SVM and LS models show
double descent curves, but their shape is different. For the SVM model
estimated using ‘clean’ data (0% label noise), there is no visible first
descent at all, but for noisy data we observe both first and second descent.
For the LS model, we clearly observe first and second descent for both clean
and noisy training data. For LS curves, the interpolation threshold ($\approx$
training size 800) is the same for different noise levels, but for SVM the
value of interpolation threshold increases with noise level in the data. This
can be explained by noting that for SVM, the interpolation threshold is
reached when the training data becomes linearly separable (in nonlinear
feature space, or Z-space in Figure 1). Therefore, for SVM the interpolation
threshold is given by the number of support vectors needed to separate
training data. For noisy data, a SVM model requires a larger number of support
vectors, resulting in larger interpolation threshold.

The ability to generalize for noisy data can be explained by noting that
during second descent VC bound (2) depends only on VC-dimension (the norm of
weights). With increasing noise (in the data), the norm of weights increases,
resulting in degradation of test error and flattening of second descent test
error curve (as evident in Figure 4, for both SVM and LS).

Empirical results in Figure 4 also show that for noisy data, generalization
performance during second descent degrades, relative to optimal first descent
model. This is contrary to the popular view that DL networks usually provide
superior generalization performance during second descent [1, 3, 4, 12].

We suspect that superior performance during second descent, reported in the DL
community, can be explained by using large and ‘clean’ data sets (common in
Big Data). For such training data sets (of size n), generalization performance
during second descent is likely to be good, because the VC bound (2) on test
error depends only on the ratio of VC-dimension to sample size $(h/n)$.

Figure 4: Effect of noise on double descent for MNIST data, for SVM and LS
classifiers using random ReLU features.

##  3 Modeling Double Descent for Fully Connected Multilayer Networks

Figure 5: Modeling results for second descent as a function of network width
(N). Figure 6: Second descent as a function of the number of epochs, for N =
10 and N = 100.

Empirical results for a simplified network setting in Section 2 provide
insights for generalization performance of over-parameterized multilayer
networks. Such general DL networks use SGD training that keeps the norm of
weights small, so that the model complexity is determined by the norm of
weights, rather than the number of weights (parameters). However, direct
application of analytic VC bounds to modeling double descent may be tricky,
because we need to address two challenging research issues:

  1. 1.

How to estimate VC-dimension for general DL networks, where analytic estimates
do not exist;

  2. 2.

Understanding design choices for setting multiple ‘tuning’ parameters, such as
network width, number of training epochs, weight initialization, etc. All of
these hyperparameters can be used to control the VC dimension of DL networks.
Double descent curves show dependence of training and test error on a single
complexity parameter, when all other tuning parameters are preset to ‘good’
values.

For these reasons, the application of analytic VC bounds to general DL
networks is difficult (or impossible?). However, it can still be done for
restricted and well-defined network settings. In this section, we consider a
fully connected network with a single hidden layer, as in Figure 1, where the
network weights in both layers are estimated during training via SGD. In this
case, z-features are adaptively estimated from training data, in contrast to
fixed random features used earlier in Section 2.

Let us consider two factors (hyperparameters) controlling complexity of such
networks trained via SGD: the number of hidden units (N), and the number of
training epochs. Empirical results showing double descent curves as a function
of these two factors have been extensively reported in DL literature [1, 12].
However, our purpose here is not to replicate such double descent curves, but
to explain them using VC bounds (1) and (3). In order to make it possible, we
have to specify network setting where VC-dimension can be approximately
estimated. Therefore, we only consider modeling for second descent, where
training error is kept very small (or zero), and test error is bounded by (2).
This can be achieved when the number of hidden units N is large, or the number
of training epochs is large. So, our experiments intend to show application of
VC bounds only in such restricted settings, where varying one complexity
factor (for example, N) has an effect on VC-dimension and test error,
according to bound (2).

Further, in order to estimate VC dimension, we hypothesize that under such
restricted settings, during second descent, the norm of weights in the output
layer can be used to approximate the VC-dimension of a neural network. The
reasoning (behind this hypothesis) is that for such restricted settings the
training error is zero, so a linear decision boundary in Z-space should have a
large separation margin, i.e. VC dimension is controlled mainly by the norm of
weights in the output layer. This assertion appears to be supported by
experiments for fully connected networks, trained on several different data
sets.

Next, we present empirical results for MNIST digits data, under the following
experimental setting:

  * –

200 training and 2,000 test examples (of digits 5 and 8);

  * –

Fully connected network using random ReLU activation function in hidden units
[24];

  * –

Training using SGD with learning rate 0.001 and momentum 0.95. The learning
rate is reduced by 10% for every 500 epochs. Batch normalization is used
during training.

  * –

Weights initialized, prior to training, using Xavier uniform distribution,
following [25].

Our design choices for SGD implementation mainly follows earlier studies [1,
25].

Figure 5 shows modeling results for second descent mode, as a function of the
number of hidden units N (the number of epochs is set to 6000 in all
experiments). The top part shows empirical training and test error curves, and
the VC-bound on test error, that closely approximates empirical test error.
The bottom part of the figure shows the VC-dimension, estimated as the norm
squared of the output layer weights (for trained network). Figure 6 shows
modeling results for second descent mode, as a function of the number of
epochs (for networks with N = 10 and 100 hidden units). Note that VC-bound and
VC-dimension can be reliably estimated only in second descent mode, when
training error is close to zero. This region, where training error is smaller
than 1%, is indicated by dotted vertical line. These results show that in the
region where training error is very small, increasing the number of epochs
results is a small increase in VC-dimension and a slight decrease of training
error. This is a particular form of memorization-complexity trade-off,
implicit in VC bound (1), when training error is very small (close to zero).

These results demonstrate applicability of VC bounds for modeling second
descent in fully connected multilayer networks. In addition, we can see the
effect of each complexity parameter (network size N and the number of epochs)
on VC-dimension. This can be used for ranking tuning parameters, according to
their ability to control VC-dimension of DL networks during second descent.

##  4 Summary

This paper provides a VC-theoretical explanation of ‘double descent’ in
multilayer networks. We show that for simplified network structures where VC-
dimension can be analytically estimated, VC generalization bounds can be
applied directly to the model and predict a double descent phenomenon.

VC-theoretical framework is helpful for improved understanding of empirical
results observed in DL, such as: the effect of various heuristics on
generalization, relative performance of first and second descent, etc. Another
important VC-theoretical insight is that during second descent VC-dimension
depends on the norm of weights. According to VC-theoretical explanation,
second descent occurs when zero training error is achieved using an estimator
having small VC-dimension, i.e. small norm of weights. This phenomenon is
general, and it does not depend on a particular training algorithm or on a
chosen parameterization (such as multilayer network). Therefore, double
descent can be observed for other learning methods, such as SVM estimators,
and not only for DL networks trained by SGD algorithm.

Empirical results presented in this paper contradict the consensus opinion
that VC-theory cannot explain prediction performance of neural networks.
Possible future work in this area may investigate VC theoretical modeling of
double descent for low-dimensional data, and also for regression problems.
Note that for regression problems VC-theoretical bounds have a different form
[8, 9, 10, 11], and these bounds have not been previously used for modeling
second descent.

## Acknowledgments and Disclosure of Funding

The authors would like to thank Dr. Vladimir Vapnik for providing valuable
suggestions and insights on application of VC-bounds for modeling double
descent.

## References

  * [1] Belkin, M., Hsu, D., Ma, S. & Mandal, S. (2019) Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences, 116(32):15849-15854.
  * [2] Bengio, Y. (2009) Learning deep architectures for AI. Foundations and Trends in Machine Learning, 2(1):1-127.
  * [3] Zhang, C., Bengio, S., Hardt, M., Recht, B. & Vinyals, O. (2021) Understanding deep learning requires rethinking generalization. Communications of the ACM, 64(3):107-115.
  * [4] Neyshabur, B., Tomioka, R. & Srebro, N. (2014) In search of the real inductive bias: On the role of implicit regularization in deep learning. arXiv preprint arXiv:1412.6614.
  * [5] Dinh, L., Pascanu, R., Bengio, S. & Bengio, Y. (2017) Sharp minima can generalize for deep nets. In International Conference on Machine Learning, pages 1019-1028.
  * [6] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. & Salakhutdinov, R. (2014) Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15(56):1929-1958.
  * [7] Bartlett, P., Foster, D. & Telgarsky, M. (2017) Spectrally-normalized margin bounds for neural networks. In Advances in Neural Information Processing Systems, pages 6241-6250.
  * [8] Vapnik, V. (1998) Statistical Learning Theory. John Wiley & Sons, New York.
  * [9] Vapnik, V. (1999) The Nature of Statistical Learning Theory, Springer.
  * [10] Vapnik, V. (2006) Estimation of Dependencies Based on Empirical Data, Second Edition, Springer.
  * [11] Cherkassky, V. & Mulier, F. (2007) Learning from Data, Second Edition, Wiley-Interscience.
  * [12] Nakkiran, P., Kaplun, G., Bansal, Y., Yang, T., Barak, B. & Sutskever, I. (2021) Deep double descent: Where bigger models and more data hurt. Journal of Statistical Mechanics: Theory and Experiment, 2021(12).
  * [13] Schölkopf, B., Smola, A.J. & Bach, F. (2002) Learning with kernels: support vector machines, regularization, optimization, and beyond, MIT press.
  * [14] Tewari, A. & Bartlett, P. (2014) Learning Theory. Academic Press Library in Signal Processing.
  * [15] Bartlett. P., Montanari, A. & Rakhlin, A. (2021) Deep learning: a statistical viewpoint, Acta Numerica.
  * [16] Koltchinskii, V. (2001) Rademacher penalties and structural risk minimization, IEEE Transactions on Information Theory, 47(5):1902-1914.
  * [17] Sokolić, J., Giryes, R., Sapiro, G. & Rodrigues, M. (2017) Robust large margin deep neural networks, IEEE Transactions on Signal Processing, 65(16):4265-4280.
  * [18] Belkin, M., Hsu, D. & J. Xu. (2020) Two Models of Double Descent for Weak Features, SIAM Journal on Mathematics of Data Science, 2(4).
  * [19] Neal, B., Mittal, S., Baratin, A., Tantia, V., Scicluna, M., Lacoste-Julien, S. & Mitliagkas, I. (2018) A modern take on the bias-variance tradeoff in neural networks. arXiv preprint arXiv:1810.08591.
  * [20] Hastie, T., Montanari, A., Rosset, S. & Ryan J. (2019) Surprises in high-dimensional ridgeless least squares interpolation. arXiv preprint arXiv:1903.08560v5.
  * [21] Rahimi, A. & Recht, B. (2008) Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, pages 1177-1184.
  * [22] Belkin, M., Hsu, D., Ma, S. & Mandal, S. (2018) To understand deep learning we need to understand kernel learning. In International Conference on Machine Learning, pages 541- 549.
  * [23] Gunasekar, S., Woodworth, B., Bhojanapalli, S., Neyshabur, B. & N. Srebro. (2017) Implicit regularization in matrix factorization. In Advances in Neural Information Processing Systems, pages 6151-6159.
  * [24] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L. & Desmaison, A. (2019) Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing systems, pages 8026-8037.
  * [25] Glorot, X. & Bengio, Y. (2010) Understanding the difficulty of training deep feedforward neural networks. In Proceedings of Machine Learning Reseach, pages 249-256.

## Appendix

This Appendix includes additional empirical results for modeling double
descent in a single-layer network shown in Figure 1 in the main paper. First
set of results investigates the choice of theoretical constants $a_{1}$ and
$a_{2}$ in VC bound (1) reproduced below:

$R_{tst}\leq
R_{trn}+\frac{\varepsilon}{2}\left(1+\sqrt{1+\frac{4R_{trn}}{\varepsilon}}\right)$
$\text{where
}\varepsilon=\frac{a_{1}}{n}\left(h\left(\ln\left(\frac{a_{2}n}{h}\right)+1\right)-\ln{\frac{\eta}{4}}\right),\eta=\min\left(\frac{4}{\sqrt{n}},1\right)$

VC theory [8, 9] specifies their range and provides the values corresponding
to pessimistic assumptions (about unknown data distributions):

  * –

the range [0, 4] for $a_{1}$ and [0,2] for $a_{2}$.

  * –

worst-case correspond to values $a_{1}=4$ and $a_{2}=2$.

These worst-case values result in upper bounds that are too crude for real-
life data sets. Therefore, for low-noise data sets in the main paper, we used
values $a_{1}=1$ and $a_{2}=1$. However, for noisy data we should use larger
values. The choice of proper values for these theoretical constants for noisy
data is discussed next, using LS classifiers with ReLU features, for digits
data (the same 5 vs. 8 data set as in the main paper). For this data set, we
introduce noise by using corrupted class labels. That is, we consider 3 data
sets, with 0% noise, 5% noise and 10% noise, where 0% noise refers to original
‘clean’ data set (used in the main paper). Figure 7 shows modeling results for
data set with 5% noise, using VC bounds with $a_{1}=1$ and $a_{2}=1$. These
results show that VC-bounds underestimate empirical curves for both first and
especially second descent. However, using values $a_{1}=3$ and $a_{2}=1$
results in practical VC bounds that provide accurate modeling of double
descent for this data, at various noise levels. See empirical results in
Figure 8.

Figure 7: Modeling double descent for digits data with 5% label noise, using
values $a_{1}=1$ and $a_{2}=1$. Figure 8: Modeling double descent for digits
data with corrupted class labels, using values $a_{1}=3$ and $a_{2}=1$.

Next, we show experimental results for the same digits data set, when images
are corrupted by random Gaussian noise. Here, the noise level is given by
standard deviation of the Gaussian noise $\hat{\sigma}=0,0.1,0.2$. Empirical
results in Figure 9 show that VC bounds (with values $a_{1}=3$ and $a_{2}=1$)
provide accurate modeling of double descent, at various noise levels.

Figure 9: Modeling double descent for digits data with corrupted pixel, using
values $a_{1}=3$ and $a_{2}=1$.

We can conclude that these values $a_{1}=3$ and $a_{2}=1$ provide robust VC
theoretical modeling of double descent for noisy data. For example, Figure 10,
shows modeling double descent for CIFAR10 data set (cat vs automobile),
extracted from CIFAR10 data base. This data set has 800 training samples and
2,000 test samples. Modeling results are obtained for the network with random
ReLU features, trained using LS classifier.

Figure 10: Modeling double descent for CIFAR data (Cat vs Automobile), using
values $a_{1}=3$ and $a_{2}=1$.

Last set of results shows the effect of varying the number of training samples
on test error, for a fixed-size network. This setting was used in [12].
Results in Figure 11 show modeling double descent for digits data set, using
fixed-size network with N=500 and N=1500. These results are obtained for the
network with random ReLU features and are trained using LS classifier with
$a_{1}=1$ and $a_{2}=1$. They show very accurate modeling of double descent
using VC-bounds, under this setting.

Figure 11: Modeling the effect of varying the number of training samples on
test error, for a fixed-size network.