# On the abnormal temperature dependent elastic properties of fused silica
irradiated by ultrafast lasers

Pieter Vlugter<EMAIL_ADDRESS>Yves Bellouard  Galatea Lab, STI/IEM,
Ecole Polytechnique Fédérale de Lausanne (EPFL), 2002 Neuchâtel, Switzerland

###### Abstract

Materials with thermal-invariant elastic properties of materials are of
interest for resonant device frequency and dimensional stability of precision
devices. Here, we demonstrate that the temperature coefficient of elasticity
(TCE) of amorphous silica can be locally reduced using femtosecond laser
irradiation. Notably, a self-organized laser-induced modification shows a
decrease of the TCE beyond 50%. The origin of this dramatic decrease is
attributed to the presence of sub-micron permanent densified zones within the
fine multi-layers structure of the laser-modified zone.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

Amorphous silica a-SiO2 has unique behavior among materials in general. Unlike
many other glasses and the vast majority of materials in general, it exhibits
a positive temperature coefficient of elasticity (TCE), i.e. its Young’s
modulus increases with increasing temperature as reported by R.B. Sosman in
1927 Uhlmann (2012). The cause for this unusual behavior has since then been
disputed. One recurring idea, proposed by Babcock in 1954 Babcock  _et al._
(1954), is that a-SiO2 contains two types of structure, coexisting within the
glass, each of which behaves normally locally, where their relative amounts
change continuously with temperature and pressure. According this model, the
material undergoes a reversible displacive transformation from an elastic weak
state ($\alpha$-form) to a stiffer one ($\beta$-form), similar to
cristobalite, one of the cristalline polymorphs of SiO2. While for the
crystalline phases, the transformation, $\alpha\rightleftharpoons\beta$,
occurs in a relative small temperature window, it would transform over a broad
temperature range in the case of a-SiO2. It is widely acknowledged that the
structure of a-SiO2 is not fully random and possess a short-range order (¡
0.5nm), while the medium-range structure (0.5 nm-1.5 nm) is consisting in a
network of rings of various sizes formed by SiO2 tetrahedra Zachariasen
(1932); King (1967); Galeener (1982), where the six-membered ring structures
is the most abundant member. According to molecular dynamics simulations done
by L. Huang _et al._ , these six-number ring-structures are specifically
responsible for the stiffening effect, by conformation changes of the ring
structures Huang and Kieffer (2004), see Figure 1. The four coordinated
tetrahedra SiO2 act as a rigid body unit Hammonds  _et al._ (1996). During the
short-range thermal expansion, the distance between Si and O atom of the four
coordinated tetrahedra SiO2 molecule Shi  _et al._ (2020), causes the six-
membered ring structures to transform from an asymmetric ($\alpha$-form) to a
more symmetric, and stiffer, configuration ($\beta$-form).

Figure 1: An impression of the six-membered ring’s conformation changes into
an $\alpha$ to the $\beta$ form due to temperature increase, adapted from
Huang and Kieffer (2004); Deng  _et al._ (2020). Where the red and white atom
represents silicon and oxygen, respectively.

a-SiO2’s TCE is altered while undergoing irreversible transformations towards
different metastable amorphous states, a transformation mechanism coined by M.
Grimsditch as ’amorphous polyamorphism’ Grimsditch (1984). Two amorphous state
families are primarily identified, the low (LDA) and high density amorphous
states (HDA), respectively. Interestingly, like other abnormal properties (the
boson peak Nakayama (2002), softening upon compression Bridgman (1925)),
permanently densified a-SiO2 reduces the magnitude of the TCE. Furthermore,
unlike amorphous ice, a-SiO2 possesses a broad range of states, within the
categories of LDA and HDA states, that can be reached through various
thermodynamic paths McMillan2007-jy or using various intense irradiation
methods that stimulate bond-breaking processes and localized reordering
Schenker and Oldham (1997); Reghioua  _et al._ (2019). As there is no unique
path towards a given state, different densification methods, such as fast
neutron exposure Strakna  _et al._ (1963), thermal quenching Le Parc  _et al._
(2006) and cold- and hot-compression Guerette  _et al._ (2018), will
inherently lead to distinct densification state, amorphous structure and TCE
reduction level.

Here, we focus our attention on the local changes of a-SiO2’ TCE resulting
from femtosecond laser pulses bulk-exposure. Thanks to their high peak-power
that triggers non-linear absorption processes, and despite having a wavelength
for which the material is transparent, these lasers are capable of modifying
locally transparent materials only in the focal volume Stuart  _et al._
(1995); Glezer  _et al._ (1997). In a-SiO2, the modified volume is associated
with microscale material property changes, e.g. refractive index Davis  _et
al._ (1996), chemical susceptibility Marcinkevi Ius  _et al._ (2001), thermal
conductivity Bellouard  _et al._ (2008), thermal expansion coefficient Vlugter
_et al._ (2019) and Young’s modulus Bellouard  _et al._ (2006); Vlugter and
Bellouard (2020). These laser-induced physical properties changes are
particularly useful in the context of monolithic micro-devices and led to the
creation of monolithic devices, e.g. complex integrated photonics circuits
Minoshima  _et al._ (2002), tuned resonators Yang and Bellouard (2017);
Vlugter and Bellouard (2020), optomechanical devices Bellouard  _et al._
(2005) and optofluidics He  _et al._ (2014); Haque  _et al._ (2014); Schaap
_et al._ (2011), made from a single piece of material and functionalized
through laser exposure sequences were reported. Specifically, a-SiO2 is
attractive for its remarkable material properties, such as its low thermal
expansion, its low optical losses over a broad spectral range as well as its
low mechanical losses and high thermal stability. Unfortunately, its TCE
remains not negligible, $\sim$210 ppm Spinner (1956); Bao and Chen (2016),
especially in contrast with nickel-based alloys commonly referred as Elinvar
or Nivarox - that have a TCE around and below ppm levels Kikuchi  _et al._
(1987). To further enhance a-SiO2 applicability, we explore how femtosecond
laser-irradiated zones can also be used to lower the intrinsic TCE value of
the material. Specifically, we investigate two types of laser-induced
modifications: one, commonly referred as ’type I’, consisting of an
homogeneous modification throughout the exposed zone, and a second one,
usually called ’type II’, consisting of self-organized nanoplanes parallel one
to another.

## Experimental methods

TCE properties are retrieved using a method based on measuring the mechanical
resonance of laser modified cantilevers as a function of temperature. We take
advantage of a process based on femtosecond laser exposure and chemical
etching to fabricate micro-cantilevers. The volume inside those cantilevers
are exposed to ultrafast laser as depicted in Figure 2a. The fabrication of
these structures, dimensions and exposure strategy are further detailed in
Vlugter and Bellouard (2020).

The resonance frequency dependency on temperature, also referred as
temperature coefficient of frequency TCF, expressed as function of material
properties Melamud  _et al._ (2009) is;

$\displaystyle\textrm{TCF}=\frac{1}{f_{0}}\frac{\partial f}{\partial T}$
$\displaystyle=$ $\displaystyle(\textrm{CTE}+\textrm{TCE})/2$ $\displaystyle=$
$\displaystyle\left(\textrm{CTE}+\frac{1}{E_{0}}\frac{\partial E}{\partial
T}\right)/2$

Where $f$ is the resonance frequency and $E$ is the Young’s modulus - the
subscript 0 refers to measured property at room temperature. Additionally, the
TCF depends on the thermal expansion coefficient (CTE). Nevertheless, the low
CTE of a-SiO2 ($\sim$ 0.5 ppm/∘C) makes this dependence negligible compared to
the TCE ($\sim$ 220 ppm/∘C) which is three orders of magnitude higher Spinner
(1956); Bao and Chen (2016). As we will see later, the laser induced CTE
changes, reported in Vlugter  _et al._ (2019), are likewise negligible
compared to the laser induced TCE changes.

Figure 2: a) a schematic of the mechanical resonator with in red the laser
affected zones (LAZ). The insets above show two scanning electron microscope
(SEM) images, typical cross-section of type I and II modifications. A
schematic further depicts a set of nanograting modifications spaced by
$l_{s}$, and consisting of alternating porous (TCEpl) and non-porous layers
(TCEnl) planes. $\vec{v}_{s}$ is the laser writing direction, $\vec{k}$ is the
optical propagation axis and $\vec{E}_{el}$ is the direction of the electric
field defining the polarization of the laser. The measured TCF and related TCE
are shown in graphs b-f). b refers to the case of a type I modification. c and
e, show the result in the case of a type II modification for nanogratings
oriented along the cantilever long-axis while, d and f) for a type II
modification, but here with the nanogratings oriented in the transverse
direction to the cantilever long-axis.

The resonator volume is not entirely and not homogeneously exposed to the
laser. Instead, the modified zones consists of a set of lines, spaced by a
distance $l_{s}$, and located just below the surface (Figure 2a). The TCE
extracted from the measured TCF expressed in equation Experimental methods is
the sum of pristine (TCEs) and laser affected zones (TCEl). The change of TCE
in the modified volume TCEl is estimated based on a mathematical model,
detailed in the Appendix A.

The temperature is controlled and cycled between 10 ∘C to 40 ∘C, using a
dedicated environmental chamber described in Vlugter  _et al._ (2019). The
cantilever is mounted on a piezoelectric-actuator that excites it on its
resonant frequency. The resonator’s response is measured using a triangulation
sensor (Keyence LKH022) and a phase-lock-loop scheme is applied to track the
resonance frequency changes during temperature cycles. Type I modifications
are produced using a femtosecond laser system emitting 50 fs-pulses at 850 nm
(Satsuma seeding an OPA from Amplitude), while type II modifications are
obtained with an Yb-doped fiber-amplified system emitting 380 fs-pulses at
1030 nm (Yuzu from Amplitude). In the sequel, the pulse energies used to
obtain the modified zone are kept constant (250 nJ for type I and 240 nJ for
type II), while the exposure dose is varied by changing the scanning speed. As
for the substrates, we consider two types of a-SiO2: one with a high OH
content (Corning, ref. 7980) and one with a low OH content (Heraeus, Suprasil
3002).

## Results

Figure 2b shows the TCF and the corresponding TCEl variation as a function of
the deposited energy for type I modifications. With these experimental
conditions, we observe a TCF decrease of 7 ppm/∘C and estimate a TCE decrease
of about 20%. Figure 2c and d, considers the TCF and estimated TCE versus the
deposited energy are depicted for a type II modifications for two opposite
nanograting orientations. The change of TCE shows a distinct behavior in both
cases. A saturating behavior for the case displayed in Figure 2d and f,
corresponding to a transverse nanograting orientation with respect to the
cantilever long-axis is found, while a linear decay for the case in Figure 2c
and e is observed and corresponds to nanogratings aligned along the cantilever
long-axis. This observation is useful to understand better what the origin of
the reduction is and will be discussed more in depth in the following section.
For type II modifications, the TCE reductions are of larger magnitude compared
to type I modifications. For nanograting aligned along the cantilever long-
axis (Figure 2c and e) the reduction is eventually reaching values
$\Delta$TCE$\sim$60$\%$ for a deposited energy of 160 J/mm2. The LAZ’s CTE
changes are negligible compared with the TCE changes, based on data reported
in Vlugter  _et al._ (2019), which confirms our earlier stated hypothesis.

## Hypothesis towards the origin of the TCE reduction

Considering the nanogratings morphology Lancry  _et al._ (2013) and the laser
patterns considered here, several features could be responsible for the
observed TCE reduction. Namely, the porous layers (TCEpl) and/or the non-
porous layers (TCEnl), permanent densifications, elastically deformed material
surrounding the LAZ (TCEϵ) and laser induced point defects are possible
candidates to explain the observed behavior. Figure 3a illustrates the
characteristic laser induced material changes and are the possible origins for
the TCE reduction.

Figure 3: Cross-section of two adjacent laser-written tracks with different
regions, characteristics of a type II modification (so called ’nanograting’
structure). TCEnl and TCEpl are the temperature coefficient of elasticity of
the non-porous and porous layers, respectively. TCEs and TCEϵ are the
temperature-dependence coefficient of elasticity of the pristine material and
the strained region, respectively. The red contour map illustrates the strain
induced by the laser induced volume expansion Champion  _et al._ (2013).
Further, the increase of point defect density is illustrated and the increase
of the three- and four-membered rings.

### Elastic deformation

The porous character is responsible for a net volume expansion of the modified
volume Champion and Bellouard (2012). Consequently, the surrounding unexposed
material is elastically densified. From compression studies in the elastic
regime, it is known that the athermal Young’s modulus softens Deschamps  _et
al._ (2014); Bridgman (1925). This could indicate that the symmetric ring
structures undergo a displacive transformation towards a weaker asymmetric
form, due to an external pressure Huang and Kieffer (2004). Following this
hypothesis, one can expect that this permanent external pressure prevents this
weaker asymmetric form to transform back into its stiffer symmetric
configuration. As a consequence, the TCE would reduce as the stiffening
mechanism (described in the introduction section) is obstructed due to the
presence of post-laser exposure residual permanent stress. To support this
hypothesis, studies on other abnormal material properties of SiO2 (i.e. boson
peak Deschamps  _et al._ (2011, 2009); Hemley  _et al._ (1986)) show that
elastic densification does reduce the magnitude of these abnormal behaviors.
Unfortunately, to the best of our knowledge, no experimental result has been
reported related to the correlation between elastic deformation and TCE.

Figure 4: Stress ($\hat{\sigma}_{l,0}$) and TCF versus laser-written line-
spacing.The exposure conditions are; pulse energy 300 nJ, repetition rate 750
kHz, scanning velocity 5 mm/s, pulse width 300 fs and laser-wavelength 1030
nm.

In the sequel, we make a comparative study between an experiment for probing
the stress present in LAZs and the one to measure the TCE as described
earlier. To probe the stress, a bimorph deflection experiment is used, like in
Champion and Bellouard (2012), where the tip displacement of an laser-exposed
cantilever is related to the stress using a modified Stoney’s equation as
detailed in Appendix B.

Unlike the previously discussed experiments, the exposure conditions (i.e.
pulse energy and exposure dose) are kept constant, but the line spacing
between individual laser-written lines is varied ($l_{s}$ in Figure 2a). This
parameter effectively controls the average stress present in laser-written
patterns, see Figure 4. On the one hand, the stress seems to saturate for
line-spacing below 4 $\mu$m (green curve in Figure 4), which is in agreement
with our previous observations Vlugter  _et al._ (2019). On the other hand,
the TCF does not follow the same trend, but rather, seems to be linearly
proportional with line-spacing. There is therefore no obvious correlation
between the elastic deformation of the unexposed material in between laser-
written line and the reduction of TCE. This agrees well with the MD-
simulations reported by Huang _et al._ Huang and Kieffer (2004), where they
concluded that the $\alpha$-$\beta$ transitions are saturated beyond the
pressure window ranging from -9 to 6 GPa. Therefore, for the pressure levels
obtained in our experiments, a sufficient amount of asymmetric ring-structures
can be thermally transformed to his stiffer counterpart.

### Porous and non-porous layer

At this point of the discussion, it is excluded that the elastic deformation
has an effect on the observed TCE reduction. To quantify whether porous or
non-porous layers are contributing to the reduction of the TCE, we revisit the
results regarding the two orthogonal nanograting orientations, depicted in
Figure 2. For the nanogratings aligned along the loading direction (Figure 2,
the TCE of the non-porous layers (TCEnl) are dominant, since it is equivalent
to a mechanical parallel configuration. Furthermore, the non-porous layer is
nearly a factor ten stiffer than porous layer and occupies nearly 90% of the
volume Vlugter and Bellouard (2020). For the orthogonal arrangements (Figure
2d and f), the porous layers contribution is this time prominent, as this
arrangement of layers is equivalent to in-series connected stiffness. This
reasoning is supported by estimates made in equation LABEL:equa_BetaPol0 that
considers nanogratings aligned along the loading direction and equation Porous
and non-porous layer, for which they are aligned transversely. For
$V_{nl}\sim$ 0.9,$V_{pl}\sim$ 0.1 (Bricchi and Kazansky (2006)) and
$E_{pl}\sim$ $E_{nl}$0.1 (Vlugter and Bellouard (2020)), equation
LABEL:equa_BetaPol0 reveals that the behavior of $\textrm{TCE}^{\perp}_{l}$ in
Figure 2c and e is mostly governed by non-porous layers contributions to the
TCE, while for an orthogonal nanograting orientation, equation Porous and non-
porous layer, both non- and porous- layers contribute nearly equally,
resulting in the behavior observed in Figure 2d and f.

$\displaystyle\textrm{TCE}^{\perp}_{l}=$
$\displaystyle\frac{\textrm{TCE}_{nl}V_{nl}E_{nl}+\textrm{TCE}_{pl}V_{pl}E_{pl}}{E_{nl}V_{nl}+E_{pl}V_{pl}}\hskip
14.22636pt\rightarrow$ $\displaystyle\textrm{TCE}^{\perp}_{l}\sim$
$\displaystyle\textrm{TCE}_{nl}0.9+\textrm{TCE}_{pl}0.01$

$\displaystyle\textrm{TCE}^{\Neswline}_{l}=$
$\displaystyle\frac{\textrm{TCE}_{nl}V_{nl}E_{pl}+\textrm{TCE}_{pl}V_{pl}E_{nl}}{V_{nl}E_{pl}+V_{pl}E_{nl}}\hskip
14.22636pt\rightarrow$ $\displaystyle\hskip
14.22636pt\textrm{TCE}^{\Neswline}_{l}\sim$
$\displaystyle\textrm{TCE}_{nl}0.9+\textrm{TCE}_{pl}$

### Annealing

To further test the hypothesis that permanently densified zones (i.e. zones
within the laser modified volume and in between nano-porous layers) are the
main source for TCE reduction, we investigated the behavior of the cantilever
after various thermal annealing steps. Depending on the temperature it is
performed, thermal annealing of laser affected zones have different effects.
For instance, for low annealing temperatures ($T_{ann}<\textrm{500}$ ∘C),
laser induced point defects (e.g., dangling bonds, non-bridging oxygen hole
centers (NBOHC), etc.) are specifically altered Skuja  _et al._ (2005); Sun
_et al._ (2000), while for intermediate temperatures (500 ∘C
$<T_{ann}<\textrm{900}$ ∘C), permanent densification is quenched Guerette  _et
al._ (2018); Cornet  _et al._ (2019). Finally, the porous layers are removed
for annealing temperatures beyond the so called ”glass annealing temperature”
Cavillon  _et al._ (2020), which for the Corning 7980 glass is found at
$T_{ann}\sim$1042 ∘C. Therefore, we logically expect changes in the TCE
depending on the annealing temperature that may provide further hints about
the origin of the TCE reduction.

In this set of experiments, we characterized the athermal Young’s modulus and
the TCE of post-annealed resonators under various conditions. In practice, the
resonators are placed in a nitrogen-purged furnace and the temperature is
elevated to a set of characteristics levels $T_{ann}$. Heating and cooling
rates are set below 1 ∘C/min, to avoid additional residual stress and
undesired transient effects. To study the fine changes of the TCE as function
of the annealing temperature, we reuse the specimens with close-spaced
exposure patterns discussed in Figure 4. For investigating the changes in
athermal Young’s moduli, we reuse the samples discussed in Figure 2d. The
athermal modulus is known to decay as function of the deposited in an
inversely proportional manner Vlugter and Bellouard (2020). This reduction of
the elastic properties is attributed to the porous planes Vlugter and
Bellouard (2020). Figure 5a and b depicts the athermal Young’s modulus (E) and
TCE normalized relative changes, respectively. The references for
normalization are the non-annealed laser-exposed material, meaning that 100%
corresponds to the non-annealed case and 0% to the unexposed one.

Overall, the laser-induced TCE and Young modulus changes remain almost
unaffected ($<$2%) for increasing annealing $T_{ann}$ up to 600 ∘C. This
indicates that point defects like E’ centers, NBOHC, ODC(I) and ODC(II) do not
have a noticeable effect on the reduction of the TCE as their concentration
gradually deplete at these annealing temperatures and durations Sun  _et al._
(2000); Zoubir  _et al._ (2006); Witcher  _et al._ (2013). Likewise,
interstitial oxygen, found in the Raman spectrum Lancry  _et al._ (2013);
Bellouard  _et al._ (2016), disappears at these annealing temperatures
Zimmermann (2017). In Sun  _et al._ (2000), it is suggested that interstitial
oxygen and E’ centers interact and form peroxy radicals (for $T_{ann}<$400
∘C). Although these chemical rearrangements highly likely have an effect on
the thermal elastic behavior, its significance is too small to measure due to
a relative low concentration.

The annealing experiments to $\Delta$E (Figure 5a) indicate that porous planes
survive at least up to a temperature of $T_{ann}<$900 ∘C. This high thermal-
stability of the porous planes have been also observed elsewhere using a
different methodology Zhang  _et al._ (2014); Wei  _et al._ (2020) and is
further extensively studied and modeled in Cavillon  _et al._ (2020); Wang
_et al._ (2021). Contrarily to the Young’s modulus, the TCE recovers at
$T_{ann}$=900 ∘C. As the porous layers are not recovered at these annealing
temperatures, it suggests that the non-porous layers are the main origin of
TCE reduction.

From Tann = 700 ∘C, the amplitude of $\Delta$TCE reduces significantly. At
$T_{ann}$ = 900 ∘C $\Delta$TCE is approximately -25%. This could either mean
that the recovered non-porous layers have not a similar molecular arrangement
as the pristine material, resulting in a negative $\Delta$TCE. Another
explanation is that the porous layers possess a positive $\Delta$TCE, since
the porous layers still exists for these annealing conditions. Effectively
this would mean that the TCE reduction of the non-porous layers is higher than
what we effectively measured, hence TCE${}_{nl}>$TCEl.

Raman spectroscopy studies by Wang et al., Witcher et al. and Zhang et al.
monitored the D1 and D2 peaks as function of the annealing temperature, these
peaks are usually attributed to a three- and four-members rings structures
Witcher  _et al._ (2013); Zhang  _et al._ (2014); Wang  _et al._ (2021).
Typically, they observed $T_{relax}$ of the D1 and D2 peaks is in the range
from 600 to 700 ∘C. The annealing behavior of these densification signatures
in the Raman spectrum correlates well with the annealing behavior of the TCE
found in this work. In a previous work, based on a dedicated study related to
the athermal Young’s modulus, we suggested that the non-porous layers in the
nanograting modification are densified Vlugter and Bellouard (2020). This fact
was also reported in earlier work by Bricchi et al., but this time, based on
birefringence measurements Bricchi and Kazansky (2006).

Figure 5: a) The change in athermal Young’s modulus ($\Delta$E) versus the
annealing temperature. b) The change in LAZ’s temperature-elastic coefficient
($\Delta TCE$) versus annealing temperatures.

Let us now consider the annealing behavior for specimens densified with other
non laser-based methods, namely, cold- Grimsditch (1986) and hot-mechanical
compression Guerette  _et al._ (2018); Cornet  _et al._ (2019), fast neutron
radiation Primak (1958) or shockwave induced densification Okuno  _et al._
(1999). Hot-compressed and shockwave densification reveal a similar trends for
densification versus $T_{ann}$, a plateau for the lower temperatures and a
decrease at higher temperatures. However, for the cold-compressed methods
(where SiO2 is compressed beyond its ’plastic’ limit of 9 GPa at room
temperature Grimsditch (1986)) and for fast-neutron radiation scenario, a
monotonic decrease is revealed Primak (1958). Shockwave-induced densification
reveals a relaxation temperature at T${}_{relax}\sim 600^{\circ}$C Okuno  _et
al._ (1999), while for hot-compressed specimens, it depends on the magnitude
of the applied compression pressure Guerette  _et al._ (2018) and temperature
Cornet  _et al._ (2019). Cornet et al. studied the relaxation of hot-densified
SiO2 using small- and wide- X-ray diffraction methods Cornet  _et al._ (2019).
They revealed that for higher compression temperatures, the densified
structure is more thermally stable (i.e. a higher Trelax). In their case,
Trelax was ranging between 400 ∘C to 850 ∘C for compression temperatures
ranging between 400 ∘C to 1020 ∘C.

The annealing behavior in this work correlates well with the annealing
behavior of the densification signatures found in the Raman spectrum, as well
with other non laser-based densification methods as can be seen in Figure 6
that provides a visual summary. Summarizing these observations, the evidences
are strong that the TCE reduction measured in this work is due to permanent
densification where parts of the material is transformed to HDA. The latter is
also mentioned in Wang  _et al._ (2021) .

Figure 6: An overview of the relaxation temperatures during annealing for
laser induced material changes and for non-laser induced densifications in
a-SiO2. Great care should be taken when extracting precise values from this
graph, as the exact values depend on the exact material composition (e.g. OH
or Cl content) and annealing conditions (e.g. duration and heat-up/cool-down
rate). Nevertheless, it provides a useful global overview.

To summarize the preceding subsections, we have provided evidences that the
non-porous layer is the main origin for the reduction of the TCE and not the
porous layer. In this scenario, the elastically-deformed material due to
volume changes do not play a role in the reduction of the TCE, neither do the
point defects. We gave further evidences that the TCE reduction correlates
well with permanent densification signatures, and therefore, we suggest that
the nanograting modification includes volumes that are in a high density
amorphous state. Additionally, we highlight that possible densification around
nano-pores in the porous layer do not contribute to a measurable reduction of
TCE as observed in this work.

## Discussion on the nature of the TCE reduction

Various works have demonstrated that densifying SiO2 reduces the TCE. Le Parc
obtained a TCE reduction of thermally quenching a-SiO2 at atmospheric pressure
Le Parc  _et al._ (2006). Based on Brillouin scattering measurement, the
author estimated a TCE reduction of 5% maximum. Strakna et al. studied the
reduction of TCE for fast-neutron irradiated a-SiO2 samples in the temperature
range from 100K to 300K Strakna  _et al._ (1963). They showed that the TCE is
reduced by approximately 70%. In recent work, Guerrete et al. investigated the
TCE for hot-compressed a-SiO2 Guerette  _et al._ (2018) and found a TCE
reduction of more than 80%, for compression-pressures of 6 GPa at a
temperature of 1100 ∘C. Figure 7 summarises the reduction of TCE obtained with
various densification processes, compared with the result of this study (in
the particular case of type II modification, a deposited energy of 160 J/mm2
and a polarisation perpendicular to the laser-scanning direction).

Figure 7: A summary of $\Delta$TCE obtained by different densification process
as reported in literature (Guerrete et al. 2018 Guerette  _et al._ (2018), Le
Parc et al. 2006 Le Parc  _et al._ (2006), Fraser et al. 1968 Fraser (1968),
Strakna et al. 1963, Strakna  _et al._ (1963)). The typical value obtained in
this study is also shown.

Based on this overview, the evidences are strong that compacting a-SiO2
results in a reduction of the TCE. In a densified material, the compaction
hinders conformation changes of the six-membered rings structures, due to the
lack of space Guerette  _et al._ (2018); Huang  _et al._ (2017). Furthermore,
the ring size distribution seems to have an effect on the magnitude of the TCE
reduction. For instance, densification induced by ionizing radiation, three
and four ring-members are generated at the expense of a decay of the number of
six-membered ring structures Vukelić  _et al._ (2010); Shcheblanov and
Povarnitsyn (2016). As the three- and four- membered rings are stiffer
structures, they are less prone to conformation changes Huang  _et al._ (2017)
and consequently, no stiffening effect upon a temperature increase is expected
for lower number ring-members, as this would be the case for six-membered
rings in fused silica. Works on intermediate glassRabia  _et al._ (2016);
Huang  _et al._ (2017) suggest a correlation between the magnitude of the TCE
and the amount of the six-membered rings Rabia  _et al._ (2016); Huang  _et
al._ (2017).

The TCE reduction obtained in this study seems remarkably high, when comparing
it with results obtained for the homogeneous modification (Type I), Figure 2a.
This indirectly suggests that the densification of the non-porous planes in
the nanograting modification is higher than homogeneous modification. This is
in other works also indirectly suggested, based on the study to the athermal
elastic properties Vlugter and Bellouard (2020) and refractive index Bricchi
and Kazansky (2006) of the non-porous planes. Where both the athermal elastic
properties and refractive index is significantly higher compared to the one of
the Type I modification.

The onset for densification of the Type I modification is generally attributed
to fast cooling after laser exposure Wang  _et al._ (2021) or/and due to a
bond-breaking mechanism, induced by ionization radiation, creating smaller
ring-members Shcheblanov and Povarnitsyn (2016). As reviewed above, both
mechanisms, fast thermal quenching and ionizing radiation, produce similar TCE
reductions as obtained for the homeogenous modification. The result correlates
well with similar densification mechanisms. Simultaneously, this raises the
question what other mechanism is taking part in reducing the TCE of the
nanograting regime, and this with a magnitude of approximately 60%. In Vlugter
and Bellouard (2020), we envisioned that the nanograting modification’s non-
porous planes densification mechanism is different to the one found in type I
exposure conditions. Accordingly, we further suspected that pressures
generated at the porous-plane location further densify the surrounding
material, i.e. the non-porous planes, at elevated temperatures as the volume
change of the porous planes exert a static pressure. Based on the new insights
provided here, these hypotheses might not hold. From the annealing
experiments, in this work, we can state that this static pressure of the
porous planes does not play a role. The TCE reduction, and therefore the
densification of these zones, recovers for $T_{ann}>$600 ∘C, while the static
pressure of the porous planes is still present. This observation suggests that
the onset of the densification is of temporal nature and appears during the
exposure or shortly after.

A precise justification for the onset for higher densification in the non-
porous planes is beyond the scope of this work and deserves a more dedicated
study. We would like nevertheless to point-out a possible line of thoughts.
The nano-pores found in the porous layers are a product of glass decomposition
due to the locally enhanced plasma density Bhardwaj  _et al._ (2006); Lancry
_et al._ (2013). The localized enhanced electron density Rudenko  _et al._
(2016) leads to a rapid, sub-nanosecond Bulgakova  _et al._ (2015) increase of
the lattice temperature. A shockwave is generated due to fast and local
thermal expansion of the material. One one hand, the negative pressure caused
by the rarefaction wave, results in a tensile force onto the material and
initiates a cavitation process, leading to the formation of nanopores Rudenko
_et al._ (2018). On the other hand, at the front of the shockwave, compressive
stress exists. This compressive stress might be the cause of this high TCE
reduction found in non-porous layer.

## Conclusion

With our earlier works, Vlugter  _et al._ (2019) that described the behavior
of the coefficient of thermal expansion and Vlugter and Bellouard (2020) that
addressed the change in Young’s modulus at room temperature, this work
finalizes a thermomechanical description of femtosecond irradiation effects on
a-SiO2, by providing insights related to the laser-induced changes induced to
the temperature coefficient of elasticity (TCE).

Here, we observed that the local decrease of the temperature coefficient of
elasticity of a-SiO2 exposed to femtosecond laser irradiation can be
significant, eventually reaching 60% and among the highest reported effects
for the various methods reported to decrease the TCE of silica.

These findings demonstrate the opportunity of tailoring the thermomechanical
behavior of dielectrics in a non-contact and localised manner. This is useful
for e.g. watchmaking, frequency references and precision devices which
requires high dimensional stability. Furthermore, the applied methodology in
this work can be used for exploring the TCE reductions in other glass formers
using ultrafast laser.

## Acknowledgments

The Galatea Lab acknowledges the sponsorship of Richemont International SA.

##  Appendix A Derivation of resonance frequency - material property relation

Figure A.1: A schematic drawing of the resonator with the LAZ in red. The
dimensions, second moment of area $I_{\cdot}$ and material properties
($E_{\cdot}$, TCE⋅, CTE⋅).

In this Appendix we express the resonator’s resonance frequency in terms of
material properties of the laser exposed and pristine material, since it is
exposed partly, see Figure A.1. The lowest resonant frequency of a clamped-
free cantilever is expressed in equation A.1. $\rho$ and $V$ are the mass
density and the volume of the cantilever. E. and I. are the Young’s modulus
and the second moment of inertia, respectively, with the subscripts .l and .s
referring to the laser affected zones (LAZ) and the non-exposed ones,
respectively. The circumflex as diacritical mark indicates that the variable
is a volume average and the subscript .0 indicates a temperature invariant
variable. $\lambda$ is the wavenumber, which corresponds to a resonance mode.

$f=\frac{\lambda^{2}}{2\pi}\sqrt{\frac{\hat{E}_{l}I_{l}+E_{s}I_{s}}{\rho
Vl^{3}}}$ (A.1)

with $\rho=\rho_{0}/(1+\textrm{CTE}\Delta T)$, $V=V_{0}(1+\textrm{CTE}\Delta
T)$, $l=l_{0}(1+\hat{\textrm{CTE}}_{l}\Delta T)$,
$I_{l}=I_{l,0}(1+\hat{\textrm{CTE}}_{l}\Delta T)^{4}$ and
$I_{s}=I_{s,0}(1+\textrm{CTE}_{s}\Delta T)^{4}$, equation A.2 becomes

$f=\frac{\lambda^{2}}{2\pi}\sqrt{\frac{\hat{E}_{l}I_{l,0}(1+\hat{\textrm{CTE}_{l}}\Delta
T)+E_{s}I_{s,0}(1+\textrm{CTE}_{s}\Delta T)}{\rho_{0}V_{0}l_{0}^{3}}}$ (A.2)

The temperature-dependent elastic constants are linearized to the first order,
i.e. $E_{s}=E_{s,0}+dE_{s}/dT\Delta T$. The thermal derivative of the
frequency, evaluated at $\Delta T=0$, is;

$\displaystyle\frac{\partial f}{\partial T}\bigg{\rvert}_{\Delta
T=0}=\frac{\lambda}{4\pi}\frac{1}{\sqrt{\hat{E}_{l,0}I_{l,0}+E_{s,0}I_{s,0}}}\frac{1}{\sqrt{l_{0}^{3}\rho_{0}V_{0}}}\cdots$
$\displaystyle\left[I_{l,0}\left(\frac{d\hat{E}_{l}}{dT}+\hat{E}_{l,0}\hat{\textrm{CTE}}_{l}\right)+I_{s,0}\left(\frac{d\hat{E}_{l}}{dT}+E_{s,0}\textrm{CTE}_{s}\right)\right]$
(A.3)

Simplifying further the equation, with the assumption that the terms
containing the coefficient of thermal expansion ($\textrm{CTE}_{.}$) are much
smaller (by an order of magnitude $\sim 4$) than with the ones containing
elastic constant thermal derivatives ($\partial E_{\cdot}/\partial T$), gives

$\frac{\partial f}{\partial T}\bigg{\rvert}_{\Delta
T=0}=\frac{f_{0}}{2}\frac{1}{\hat{E}_{l,0}I_{l,0}+E_{s,0}I_{s,0}}\left[I_{l,0}\left(\frac{d\hat{E}_{l}}{dT}\right)+I_{s,0}\left(\frac{dE_{s}}{dT}\right)\right]$
(A.4)

rewriting the equation in terms of $\textrm{TCE}_{.}$, i.e.
$dE_{l}/dT=E_{l,0}\textrm{TCE}_{l}$, gives,

$\frac{\partial f}{\partial T}\bigg{\rvert}_{\Delta
T=0}=\frac{f_{0}}{2}\left[\frac{I_{l,0}\hat{E}_{l,0}\hat{\textrm{TCE}}_{l}+I_{s,0}E_{s,0}\textrm{TCE}_{s}}{\hat{E}_{l,0}I_{l,0}+E_{s,0}I_{s,0}}\right]$
(A.5)

$\hat{\textrm{TCE}}_{l}$ and $\hat{E}_{l,0}$ are the average properties of the
modified volume, which are further specified in equations A.6 and A.7, in
terms of the volume fractions $V_{l}$ and $V_{s}$. Where LAZ and pristine
volume fraction is expressed as $V_{l}=w_{0}/l_{s}$ and $V_{s}=1-V_{l}$ (see
Figure A.1).

$\hat{\textrm{TCE}}_{l}=\frac{\textrm{TCE}_{s}V_{s}E_{l,0}+\textrm{TCE}_{l}V_{l}E_{s,0}}{V_{s}E_{l,0}+V_{l}E_{s,0}}$
(A.6)

$\hat{E}_{l,0}=\frac{E_{l,0}E_{s,0}}{E_{l,0}V_{s}+E_{s,0}V_{l}}$ (A.7)

##  Appendix B Derivation of bimorph stress relation

In this annex, we derive a mathematical expression for the stress in the
modified layer $\hat{\sigma}_{l,0}$ of the bimorph structure. In Figure B.1 a
schematic is depicted of a cantilever being exposed just below the top-layer
and consequently bend with a bending radius $r$ due to volume expansion of the
laser affected zone. The stress and strain profile is illustrate in the inset
A-A Figure B.1. Where the stress in the laser affected layer,
$\hat{\sigma}_{l,0}$, is opposite in sign compared with the stress in the
pristine bulk. For this system, the stress and strain in the pristine bulk is
zero at a distance $t_{s}$/3 from the bottom surface, defining the neutral-
axis Hsueh (2002). We define the neutral-axis as the origin of the
z-coordinate, i.e. z=0.

Figure B.1: A schematic of a laser-exposed cantilever, with the dimensions
indicated as well as the stress and strain profiles across the thickness.
$\hat{\sigma}_{l,0}$ is the stress in laser exposed layer.

The relationship between bending radius and strain is;

$\epsilon_{xx}(z)=\frac{z}{r}$ (B.1)

Newton’s third law is used to analytical determine $\hat{\sigma}_{l,0}$, see
equation B.2. The first two terms represents the force in pristine bulk, with
$E_{s}$ as the Young’s modulus of this volume. Where the last term represent
the force in the LAZ.

$0=\frac{\partial\epsilon_{xx}(z)}{\partial
z}E_{s}\left(\int_{-t_{s}/3}^{t_{s}2/3}z\cdot\textrm{d}z+\int_{t_{s}2/3+t_{l}}^{t_{s}2/3+t_{l}+t_{s,t}}z\cdot\textrm{d}z\right)+\hat{\sigma}_{l,0}t_{l}$
(B.2)

Rewriting equation B.2 gives the expression for $\hat{\sigma}_{l,0}$;

$\hat{\sigma}_{l,0}=-\frac{E_{s}t^{2}_{s}}{6t_{l}r}-\frac{E_{s}}{t_{l}r_{0}}\left[\frac{1}{2}\left(t_{s}\frac{2}{3}+t_{l}+t_{s,t}\right)^{2}-\frac{1}{2}\left(\frac{t_{s}}{3}+t_{l}\right)^{2}\right]$
(B.3)

Note, that the first term is the Stoney’s equation.

## References

  * Uhlmann (2012) D. Uhlmann,  _Elasticity and Strength in Glasses: Glass: Science and Technology_ (Elsevier, 2012).
  * Babcock  _et al._ (1954) C. L. Babcock, S. W. Barber, and K. Fajans, Ind. Eng. Chem. 46, 161 (1954).
  * Zachariasen (1932) W. H. Zachariasen, J. Am. Chem. Soc. 54, 3841 (1932).
  * King (1967) S. V. King, Nature 213, 1112 (1967).
  * Galeener (1982) F. L. Galeener, Solid State Commun. 44, 1037 (1982).
  * Huang and Kieffer (2004) L. Huang and J. Kieffer, Phys. Rev. B Condens. Matter 69, 224203 (2004).
  * Hammonds  _et al._ (1996) K. D. Hammonds, . M. T. Do VE, and . A. P. Giddy, pages 1057, 1079 (1996).
  * Shi  _et al._ (2020) Y. Shi, D. Ma, A. P. Song, B. Wheaton, M. Bauchy, and S. R. Elliott, J. Non-Cryst. Solids 528, 119760 (2020).
  * Deng  _et al._ (2020) B. Deng, Y. Shi, and F. Yuan, Materialia 12, 100752 (2020).
  * Grimsditch (1984) M. Grimsditch, Phys. Rev. Lett. 52, 2379 (1984).
  * Nakayama (2002) T. Nakayama, Rep. Prog. Phys. 65, 1195 (2002).
  * Bridgman (1925) P. W. Bridgman, . Am.J.Sci. 10 (1925).
  * Schenker and Oldham (1997) R. E. Schenker and W. G. Oldham, J. Appl. Phys. 82, 1065 (1997).
  * Reghioua  _et al._ (2019) I. Reghioua, M. Lancry, O. Cavani, S. L. Floch, D. R. Neuville, and N. Ollier, Appl. Phys. Lett. 115, 251101 (2019).
  * Strakna  _et al._ (1963) R. E. Strakna, A. E. Clark, D. L. Bradley, and W. M. Slie, J. Appl. Phys. 34, 1439 (1963).
  * Le Parc  _et al._ (2006) R. Le Parc, C. Levelut, J. Pelous, V. Martinez, and B. Champagnon, J. Phys. Condens. Matter 18, 7507 (2006).
  * Guerette  _et al._ (2018) M. Guerette, M. R. Ackerson, J. Thomas, E. B. Watson, and L. Huang, J. Chem. Phys. 148, 194501 (2018).
  * Stuart  _et al._ (1995) B. C. Stuart, M. D. Feit, A. M. Rubenchik, B. W. Shore, and M. D. Perry, Phys. Rev. Lett. 74, 2248 (1995).
  * Glezer  _et al._ (1997) E. N. Glezer, M. Milosavljevic, L. Huang, R. J. Finlay, T. H. Her, J. P. Callan, and E. Mazur, Opt. Lett. 22, 422 (1997).
  * Davis  _et al._ (1996) K. M. Davis, K. Miura, N. Sugimoto, and K. Hirao, Opt. Lett. 21, 1729 (1996).
  * Marcinkevi Ius  _et al._ (2001) A. Marcinkevi Ius, S. Juodkazis, M. Watanabe, M. Miwa, S. Matsuo, H. Misawa, and J. Nishii, Opt. Lett. 26, 277 (2001).
  * Bellouard  _et al._ (2008) Y. Bellouard, E. Barthel, A. A. Said, M. Dugan, and P. Bado, Opt. Express 16, 19520 (2008).
  * Vlugter  _et al._ (2019) P. Vlugter, E. Block, and Y. Bellouard, Phys. Rev. Materials 3, 053802 (2019).
  * Bellouard  _et al._ (2006) Y. Bellouard, T. Colomb, C. Depeursinge, M. Dugan, A. A. Said, and P. Bado, Opt. Express 14, 8360 (2006).
  * Vlugter and Bellouard (2020) P. Vlugter and Y. Bellouard, Phys. Rev. Materials 4, 023607 (2020).
  * Minoshima  _et al._ (2002) K. Minoshima, A. Kowalevicz, E. Ippen, and J. Fujimoto, Opt. Express 10, 645 (2002).
  * Yang and Bellouard (2017) T. Yang and Y. Bellouard, Phys. Rev. Applied 7, 064002 (2017).
  * Bellouard  _et al._ (2005) Y. Bellouard, A. Said, and P. Bado, Opt. Express 13, 6635 (2005).
  * He  _et al._ (2014) F. He, Y. Liao, J. Lin, J. Song, L. Qiao, Y. Cheng, and K. Sugioka, Sensors 14, 19402 (2014).
  * Haque  _et al._ (2014) M. Haque, K. K. C. Lee, S. Ho, L. A. Fernandes, and P. R. Herman, Lab Chip 14, 3817 (2014).
  * Schaap  _et al._ (2011) A. Schaap, Y. Bellouard, and T. Rohrlack, Biomed. Opt. Express 2, 658 (2011).
  * Spinner (1956) S. Spinner, J. Am. Ceram. Soc. 39, 113 (1956).
  * Bao and Chen (2016) Y. Bao and G. Chen, Meas. Sci. Technol. 27, 065101 (2016).
  * Kikuchi  _et al._ (1987) M. Kikuchi, H. Masumoto, and T. Nakayama, Transactions of the Japan Institute of Metals 28, 585 (1987).
  * Melamud  _et al._ (2009) R. Melamud, S. A. Chandorkar, B. Kim, H. K. Lee, J. C. Salvia, G. Bahl, M. A. Hopcroft, and T. W. Kenny, J. Microelectromech. Syst. 18, 1409 (2009).
  * Lancry  _et al._ (2013) M. Lancry, B. Poumellec, J. Canning, K. Cook, J.-C. Poulin, and F. Brisset, Laser Photonics Rev. 7, 953 (2013).
  * Champion  _et al._ (2013) A. Champion, M. Beresna, P. Kazansky, and Y. Bellouard, Opt. Express 21, 24942 (2013).
  * Champion and Bellouard (2012) A. Champion and Y. Bellouard, Opt. Mater. Express, OME 2, 789 (2012).
  * Deschamps  _et al._ (2014) T. Deschamps, J. Margueritat, C. Martinet, A. Mermet, and B. Champagnon, Sci. Rep. 4, 7193 (2014).
  * Deschamps  _et al._ (2011) T. Deschamps, C. Martinet, D. de Ligny, J. L. Bruneel, and B. Champagnon, J. Chem. Phys. 134, 234503 (2011).
  * Deschamps  _et al._ (2009) T. Deschamps, C. Martinet, D. R. Neuville, D. de Ligny, C. Coussa-Simon, and B. Champagnon, J. Non-Cryst. Solids 355, 2422 (2009).
  * Hemley  _et al._ (1986) R. J. Hemley, H. K. Mao, P. M. Bell, and B. O. Mysen, Phys. Rev. Lett. 57, 747 (1986).
  * Bricchi and Kazansky (2006) E. Bricchi and P. G. Kazansky, Appl. Phys. Lett. 88, 111119 (2006).
  * Skuja  _et al._ (2005) L. Skuja, M. Hirano, H. Hosono, and K. Kajihara, Phys. Status Solidi C 2, 15 (2005).
  * Sun  _et al._ (2000) H.-B. Sun, S. Juodkazis, M. Watanabe, S. Matsuo, H. Misawa, and J. Nishii, J. Phys. Chem. B 104, 3450 (2000).
  * Cornet  _et al._ (2019) A. Cornet, C. Martinet, V. Martinez, and D. de Ligny, J. Chem. Phys. 151, 164502 (2019).
  * Cavillon  _et al._ (2020) M. Cavillon, Y. Wang, B. Poumellec, F. Brisset, and M. Lancry, Appl. Phys. A: Mater. Sci. Process. 126, 876 (2020).
  * Zoubir  _et al._ (2006) A. Zoubir, C. Rivero, R. Grodsky, K. Richardson, M. Richardson, T. Cardinal, and M. Couzi, Phys. Rev. B Condens. Matter 73, 224117 (2006).
  * Witcher  _et al._ (2013) J. J. Witcher, W. J. Reichman, L. B. Fletcher, N. W. Troy, and D. M. Krol, OSA Optical materials express 3 (2013).
  * Bellouard  _et al._ (2016) Y. Bellouard, A. Champion, B. McMillen, S. Mukherjee, R. R. Thomson, C. Pépin, P. Gillet, and Y. Cheng, Optica 3, 1285 (2016).
  * Zimmermann (2017) F. Zimmermann,  _Ultrashort pulse induced nanostructures in transparent materials_ , Ph.D. thesis, Friedrich-Schiller-Universitat Jena (2017).
  * Zhang  _et al._ (2014) F. Zhang, Y. Yu, C. Cheng, Y. Dai, H. Zhang, and J. Qiu, Appl. Phys. B 117, 53 (2014).
  * Wei  _et al._ (2020) S.-E. Wei, Y. Wang, H. Yao, M. Cavillon, B. Poumellec, G.-D. Peng, and M. Lancry, Sensors Fiber Optic Sensors and Applications XI, 20, 762 (2020).
  * Wang  _et al._ (2021) Y. Wang, M. Cavillon, N. Ollier, B. Poumellec, and M. Lancry, Phys. Status Solidi , 2100023 (2021).
  * Grimsditch (1986) M. Grimsditch, Phys. Rev. B Condens. Matter 34, 4372 (1986).
  * Primak (1958) W. Primak, Phys. Rev. 110, 1240 (1958).
  * Okuno  _et al._ (1999) M. Okuno, B. Reynard, Y. Shimada, Y. Syono, and C. Willaime, Phys. Chem. Miner. 26, 304 (1999).
  * Fraser (1968) D. B. Fraser, J. Appl. Phys. 39, 5868 (1968).
  * Huang  _et al._ (2017) L. Huang, F. Yuan, M. Guerette, Q. Zhao, and S. Sundararaman, J. Mater. Res. 32, 174 (2017).
  * Vukelić  _et al._ (2010) S. Vukelić, P. Kongsuwan, S. Ryu, and Y. L. Yao, J. Manuf. Sci. Eng. 132 (2010).
  * Shcheblanov and Povarnitsyn (2016) N. S. Shcheblanov and M. E. Povarnitsyn, EPL 114, 26004 (2016).
  * Rabia  _et al._ (2016) M. K. Rabia, S. Degioanni, C. Martinet, J. Le Brusq, B. Champagnon, and D. Vouagner, J. Phys. Condens. Matter 28, 075402 (2016).
  * Bhardwaj  _et al._ (2006) V. R. Bhardwaj, E. Simova, P. P. Rajeev, C. Hnatovsky, R. S. Taylor, D. M. Rayner, and P. B. Corkum, Phys. Rev. Lett. 96, 057404 (2006).
  * Rudenko  _et al._ (2016) A. Rudenko, J.-P. Colombier, and T. E. Itina, Phys. Rev. B Condens. Matter 93, 075427 (2016).
  * Bulgakova  _et al._ (2015) N. M. Bulgakova, V. P. Zhukov, S. V. Sonina, and Y. P. Meshcheryakov, J. Appl. Phys. 118, 233108 (2015).
  * Rudenko  _et al._ (2018) A. Rudenko, J.-P. Colombier, and T. E. Itina, Phys. Chem. Chem. Phys. 20, 5887 (2018).
  * Hsueh (2002) C.-H. Hsueh, Appl. PHYSICS 91, 9652 (2002).


University of Technology Nuremberg, Germany and
https://stoianmihail.github.iomihail.stoian@utn.dehttps://orcid.org/0000-0002-8843-3374
Mihail Stoian <ccs2012> <concept>
<concept_id>10003752.10003809.10011254.10011258</concept_id>
<concept_desc>Theory of computation Dynamic programming</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
</ccs2012> [500]Theory of computation Dynamic programming

# TSP Escapes the $O(2^{n}n^{2})$ Curse

Mihail Stoian

###### Abstract

The dynamic programming solution to the traveling salesman problem due to
Bellman, and independently Held and Karp, runs in time $O(2^{n}n^{2})$, with
no improvement in the last sixty years. We break this barrier for the first
time by designing an algorithm that runs in deterministic time
$2^{n}n^{2}/2^{\Omega(\sqrt{\log n})}$. We achieve this by strategically
remodeling the dynamic programming recursion as a min-plus matrix product, for
which faster-than-naïve algorithms exist.

###### keywords:

traveling salesman problem, polynomial method, dynamic programming

###### category:

##  1 Introduction

The Traveling Salesman Problem (TSP) needs no introduction: It is probably one
of the most studied problems in computer science since its establishment
itself as a field (and not only; we recommend Schrijver’s [38] survey). Due to
its simplicity, namely finding the shortest tour through $n$ cities, it is
also one of the introductory problems in dynamic programming. Bellman [4] and,
in an independent work Held and Karp [31], proposed a dynamic programming
solution that runs in $O(2^{n}n^{2})$-time and takes $O(2^{n}n)$-space, using
a table indexed by vertex subsets. This has been _the_ textbook solution for
decades.

The unweighted case, namely deciding whether the graph has a Hamiltonian
cycle, can indeed be solved faster: In his breakthrough, Björklund [6] showed
a Monte Carlo algorithm that runs in $O^{*}(1.657^{n})$-time. Using a standard
technique, it can be adapted to solve our particular problem, yet with an
$O(w)$ dependence in the running time, where $w$ is the sum of all edge
weights. In this context, Nederlof [37] provided an $O(1.9999^{n})$-time
randomized algorithm with constant error probability for the _bipartite_ TSP,
_assuming the matrix-multiplication exponent, $\omega$, is equal to 2_. He
also outlines how his methods can be useful for solving TSP in non-bipartite
graphs. Degree-boundedness is also a fertile ground for such improvements [8,
18]. In the Euclidean setting, the $2^{O(1-1/d)}$ lower-bound [20] has already
been attained [19]. However, despite these advances, the standard dynamic
programming solution in time $O(2^{n}n^{2})$ is still the best algorithm for
the general setting. This was also the status quo when Cook [16] wrote his
compendium on TSP ten years ago:

“ Unfortunately for TSP fans, no good algorithm is known for the problem. The
best result thus far is a solution method, discovered in 1962, that runs in
time proportional to $n^{2}2^{n}$. — William Cook [16] ”

This leads us to our driving research question:

> Is there an $o(2^{n}n^{2})$-time algorithm for TSP?

A positive answer would break another _psychological_ barrier, as Koutis and
Williams [34] refer to it, and, if simple enough, can be used as the new
textbook algorithm for TSP.

##  2 Our Results

We answer our research question in the positive. We introduce an algorithm
that, to our best knowledge, is the first to break the longstanding
$O(2^{n}n^{2})$ time-barrier for TSP. We remodel the standard dynamic
programming recursion as a min-plus matrix product, namely:

######  Theorem 2.1 (Backbone).

If the Min-Plus Matrix Product of $n\times n$ matrices can be solved in time
$T(n)$, then TSP can be solved in time $O(2^{n}T(n)/n)$.

If we plug in the fastest algorithm for the min-plus matrix product due to
Williams [46], which was been later derandomized by Chan and Williams [12], we
obtain our claimed result:

######  Corollary 2.2 (Main).

TSP can be solved in deterministic time $2^{n}n^{2}/2^{\Omega(\sqrt{\log
n})}$.

The factor $\frac{1}{2^{\Omega(\sqrt{\log n})}}$ may seem a bit unusual to the
untrained reader. In particular, this factor is better than
$\frac{1}{\log^{k}n}$ time, for every $k$. We next outline our algorithm.

##  3 Algorithm

Let us first recall the standard dynamic programming (DP) solution:

$\displaystyle\mathrm{dp}(S,k)$ $\displaystyle=\displaystyle\min_{j\in
S\setminus\\{k\\}}\left(\mathrm{dp}(S\setminus\\{k\\},j)+c_{j,k}\right),k\in
S$ (1) $\displaystyle\mathrm{dp}(\\{1\\},1)$ $\displaystyle=0,$ (2)

where $S\subseteq[n]\vcentcolon=\\{1,\ldots,n\\}$ and $c$ is the $n\times n$
cost matrix of the instance. At the end, the final solution is computed as
$\displaystyle\min_{j}\left(\mathrm{dp}([n],j)+c_{j,1}\right)$. The algorithm
runs in time $O(2^{n}n^{2})$, since for each state $(S,k)$ we have to loop
over all $j$’s in $S\setminus\\{k\\}$. A convenient way to compute the
$\mathrm{dp}$-table is to fix a cardinality $\ell\in\\{2,\ldots,n\\}$ and
compute the $\mathrm{dp}$-values for all sets of cardinality $\ell$. In other
words, the $\mathrm{dp}$-table is computed _layer-wise_.

Key Insight. The next observation is rather crucial for understanding our
algorithm: To some extent, Eq. (1) _is indeed an inner product in the min-plus
semi-ring_ between the vector $\mathrm{dp}(S\setminus\\{k\\},\vcentcolon)$ and
the $j$th column of the cost matrix $c$ (for a fixed $k\in S$). Thus, the
computation of the entire row $\mathrm{dp}(S,\vcentcolon)$ consists of several
min-plus inner products. Next, we exploit this new perspective.

Remodeling the DP. To this end, we remark that if we represented these inner
products as a matrix product, the output could also be shared with other rows
of the current layer. With this in mind, we can _group_ the rows of the
$\mathrm{dp}$-table of the previous layer, i.e., the $(\ell-1)$th layer, into
batches of size $n$ (so that the dimension matches).111The rows of the final
batch may be padded with $\infty$ in case the batch does not have exactly size
$n$. We perform the min-plus product of the current batch with the cost matrix
and then update the rows for the current cardinality $\ell$. After all
$\lceil{n\choose\ell-1}/n\rceil$ batches have been iterated, the $\ell$th
layer of the $\mathrm{dp}$-table has been fully computed.

Pseudocode. We outline the pseudocode in Alg. 1. We first initialize the table
with $\infty$ and set up the base case, Eq. (2). We then iterate over all
cardinalities $\ell\in\\{2,\ldots,n\\}$ and generate the corresponding
batches. Specifically, the function $\textsc{Batches}(\ell-1)$ generates
internally the set $\\{S\subseteq[n]\mid|S|=\ell-1\\}$ containing all sets of
cardinality $\ell-1$ and then chunks it into batches of size $n$ (the final
batch is handled individually in the min-plus product), returning these one by
one in the variable $\mathcal{B}$. For each batch, we perform the min-plus
product between the corresponding rows of the $\mathrm{dp}$-table, i.e.,
$\mathrm{dp}(\mathcal{B},\vcentcolon)$, and the cost matrix $c$. This results
in a temporary matrix $p$ which is then iterated to fill in the values of the
current layer. Namely, for each set $\mathcal{B}[i]$, we take the $k$’s that
are not present within it and update the value of the set
$\mathcal{B}[i]\cup\\{k\\}$ with $p_{i,k}$. In particular, $p_{i,k}$ stores
the min-plus inner product between $\mathrm{dp}(\mathcal{B}[i],\vcentcolon)$
and $c_{\vcentcolon,k}$. Hence, the operation in Eq. (1), which once took
linear time, can now be done in $O(1)$, given the min-plus product $p$ was
computed before. Finally, we return the solution.

Algorithm 1 $\textsc{TSPviaMinPlusProduct}(n,c)$

1:$\mathrm{dp}(\vcentcolon,\vcentcolon)\leftarrow\infty$

2:$\mathrm{dp}(\\{1\\},1)\leftarrow 0$

3:for $\ell=2,\ldots,n$ do

4: for $\mathcal{B}\in\textsc{Batches}(\ell-1)$ do

5:
$p\leftarrow\textsc{MinPlusProduct}(\mathrm{dp}(\mathcal{B},\vcentcolon),c)$

6: for $i=1,...,|\mathcal{B}|$ do

7:
$\mathrm{dp}(\mathcal{B}[i]\cup\\{k\\},k)\leftarrow\min(\mathrm{dp}(\mathcal{B}[i]\cup\\{k\\},k),p_{i,k})$,
$\forall k\in[n]\setminus\mathcal{B}[i]$

8: end for

9: end for

10:end for

11:return $\displaystyle\min_{k}\left(\mathrm{dp}([n],k)+c_{k,1}\right)$

We show that the algorithm returns the optimal solution in the claimed running
time:

See 2.1

######  Proof 3.1.

The main observation is that Alg. 1 slightly changes the computation flow:
While Eq. (1) _pulls_ the values from the previous layer to compute the
current layer, our algorithm _pushes_ them instead. Indeed, once the min-plus
product of the current batch has been computed as matrix $p$, we can gradually
update the values of the current layer $\ell$. Since all batches of the
previous layer are iterated, the correctness of the previous algorithm is
preserved.

Let us now analyze the running time. For each layer $\ell$, there are
$\lceil{n\choose\ell-1}\rceil$ batches to iterate. In each batch, we perform a
min-plus product between the $\mathrm{dp}$-table _restricted_ to the batch
$\mathcal{B}$ (this has $n\times n$ size) and the $n\times n$ cost matrix in
time $T(n)$. We can iterate the $n\times n$ intermediate matrix $p$ to update
the values of the current layer. Thus, the total running time reads

$O\left(\sum_{\ell=2}^{n}\frac{\binom{n}{{\ell-1}}}{n}\cdot\left(T(n)+n^{2}\right)\right)=O\left(\frac{2^{n}}{n}T(n)\right).$

Let us now take a closer look at $T(n)$, the running time of the min-plus
matrix product. Fortunately, there has been a rich line of research on
computing it, culminating in Williams’ algorithm [46], which runs in time
$n^{3}/2^{\Omega(\sqrt{\log n})}$ and was later derandomized by Chan and
Williams [12]. Thus, our Alg. 1 can run in time
$2^{n}n^{2}/2^{\Omega(\sqrt{\log n})}$.

##  4 Related Work

Traveling Salesman Problem. The research on TSP is rather vast and it would
take an entire survey to cover it. Apart from those mentioned in the
introduction, approximation algorithms for TSP are also a beloved research
area: Christofides’ 1.5-approximation algorithm [15] (independently discovered
by Serdyukov [39]) is probably one of the classics of approximation
algorithms. In a breakthrough result, Karlin, Klein, and Gharan [32] showed
that for some $\varepsilon>10^{-36}$, there is an
$(1.5-\varepsilon)$-approximation algorithm. Remarkable results have also been
obtained for asymmetric TSP [3, 40] and in other settings as well [44, 49].
Beyond that, there is also a PTAS for the Euclidean setting [2, 36]. Another
favorite research area is the $k$-OPT heuristic, which is a local search
algorithm that is allowed to replace $k$ edges of a given solution to obtain a
better cost [35, 5, 30, 17]. A further prominent area is that of polynomial-
space algorithms, under which setting TSP has an $O(4^{n}n^{\log n})$-time
algorithm [26, 7]. Alternatively, one can trade off time and space, as shown
by Koivisto and Parviainen [33]. TSP also enjoys a quantum speedup, as shown
by Ambainis et al. [1].

Min-Plus Matrix Product. While our improved algorithm for TSP directly uses
the fastest deterministic algorithm for the min-plus product, which uses the
polynomial method from circuit complexity [46, 12], the line of research on
this problem has been vast, with incremental improvements over the years,
mainly for its main application, the all-pairs shortest paths (APSP) problem
[23, 45, 24, 21, 41, 27, 42, 50, 43, 11, 28, 10, 29]. Notably, there is an
open problem in fine-grained complexity whether there exists an algorithm for
the min-plus product which runs in truly subcubic time, i.e.,
$O(n^{3-\varepsilon})$ for $\varepsilon>0$ [47]. To this end, recent work has
focused on understanding for which instances this product can be solved in
truly subcubic time. Such instances are bounded-difference matrices [9, 13],
later generalized to the case where one of the matrices is of
$O(1)$-approximate rank [48], matrices with monotone rows or columns [14, 25,
22], or geometrically weighted matrices [10].

##  5 Discussion

While modest, our improvement shows that the well-known dynamic programming
solution due to Bellman [4] and, independently Held and Karp [31], designed
more than sixty years ago and taught in every algorithms course, is _not_ the
best we can hope for. This breaks another _psychological_ barrier, as Koutis
and Williams [34] refer to it.

Open Problems. A natural question is whether our algorithm can be sped up.
Given the APSP hypothesis [47], it seems that we may indeed need to resort to
special instances. For example, note that the cost matrix of the TSP instance
remains constant throughout the algorithm. This opens up the possibility of
preprocessing it, even in $O(2^{n}n)$-time, so that upcoming matrix products
can be performed faster.

In this light, it is natural to ask whether an
$O(2^{n}n^{2-\varepsilon})$-time algorithm exists for $\varepsilon>0$. Even
more interesting is the (still) open question of whether there is an
$O(1.9999^{n})$-time algorithm, for which the work of Nederlof [37] has
already paved the way.

Accessibility. For the “TSP fans” of Cook [16], we have prepared a proof-of-
concept implementation of Alg. 1 on the author’s Github
page.222github.com/stoianmihail

## References

  * [1] Andris Ambainis, Kaspars Balodis, Jānis Iraids, Martins Kokainis, Krišjānis Prūsis, and Jevgēnijs Vihrovs.  Quantum speedups for exponential-time dynamic programming algorithms.  In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1783–1793. SIAM, 2019.
  * [2] Sanjeev Arora.  Polynomial time approximation schemes for euclidean traveling salesman and other geometric problems.  Journal of the ACM (JACM), 45(5):753–782, 1998.
  * [3] Arash Asadpour, Michel X Goemans, Aleksander Mądry, Shayan Oveis Gharan, and Amin Saberi.  An o (log n/log log n)-approximation algorithm for the asymmetric traveling salesman problem.  Operations Research, 65(4):1043–1061, 2017.
  * [4] Richard Bellman.  Dynamic programming treatment of the travelling salesman problem.  Journal of the ACM (JACM), 9(1):61–63, 1962.
  * [5] Mark de Berg, Kevin Buchin, Bart MP Jansen, and Gerhard Woeginger.  Fine-grained complexity analysis of two classic tsp variants.  ACM Transactions on Algorithms (TALG), 17(1):1–29, 2020.
  * [6] Andreas Bjorklund.  Determinant sums for undirected hamiltonicity.  SIAM Journal on Computing, 43(1):280–299, 2014.
  * [7] Andreas Björklund and Thore Husfeldt.  Exact algorithms for exact satisfiability and number of perfect matchings.  Algorithmica, 52:226–249, 2008.
  * [8] Andreas Björklund, Thore Husfeldt, Petteri Kaski, and Mikko Koivisto.  The travelling salesman problem in bounded degree graphs.  In International Colloquium on Automata, Languages, and Programming, pages 198–209. Springer, 2008.
  * [9] Karl Bringmann, Fabrizio Grandoni, Barna Saha, and Virginia Vassilevska Williams.  Truly subcubic algorithms for language edit distance and rna folding via fast bounded-difference min-plus product.  SIAM Journal on Computing, 48(2):481–512, 2019.
  * [10] Timothy M Chan.  More algorithms for all-pairs shortest paths in weighted graphs.  In Proceedings of the thirty-ninth annual ACM symposium on Theory of computing, pages 590–598, 2007.
  * [11] Timothy M Chan.  All-pairs shortest paths with real weights in o (n 3/log n) time.  Algorithmica, 50:236–243, 2008.
  * [12] Timothy M Chan and Ryan Williams.  Deterministic apsp, orthogonal vectors, and more: Quickly derandomizing razborov-smolensky.  In Proceedings of the twenty-seventh annual ACM-SIAM symposium on Discrete algorithms, pages 1246–1255. SIAM, 2016.
  * [13] Shucheng Chi, Ran Duan, and Tianle Xie.  Faster algorithms for bounded-difference min-plus product.  In Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1435–1447. SIAM, 2022.
  * [14] Shucheng Chi, Ran Duan, Tianle Xie, and Tianyi Zhang.  Faster min-plus product for monotone instances.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 1529–1542, 2022.
  * [15] Nicos Christofides.  Worst-case analysis of a new heuristic for the travelling salesman problem.  In Operations Research Forum, volume 3, page 20. Springer, 2022.
  * [16] William J Cook.  In pursuit of the traveling salesman: mathematics at the limits of computation.  Princeton University Press, 2015.
  * [17] Marek Cygan, Łukasz Kowalik, and Arkadiusz Socała.  Improving tsp tours using dynamic programming over tree decompositions.  ACM Transactions on Algorithms (TALG), 15(4):1–19, 2019.
  * [18] Marek Cygan and Marcin Pilipczuk.  Faster exponential-time algorithms in graphs of bounded average degree.  Information and Computation, 243:75–85, 2015.
  * [19] Mark De Berg, Hans L Bodlaender, Sándor Kisfaludi-Bak, and Sudeshna Kolay.  An eth-tight exact algorithm for euclidean tsp.  In 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), pages 450–461. IEEE, 2018.
  * [20] Mark De Berg, Hans L Bodlaender, Sándor Kisfaludi-Bak, Dániel Marx, and Tom C Van Der Zanden.  A framework for exponential-time-hypothesis–tight algorithms and lower bounds in geometric intersection graphs.  SIAM Journal on Computing, 49(6):1291–1331, 2020.
  * [21] Wlodzimierz Dobosiewicz.  A more efficient algorithm for the min-plus multiplication.  International journal of computer mathematics, 32(1-2):49–60, 1990.
  * [22] Anita Dürr.  Improved bounds for rectangular monotone min-plus product and applications.  Information Processing Letters, 181:106358, 2023.
  * [23] Robert W Floyd.  Algorithm 97: shortest path.  Communications of the ACM, 5(6):345–345, 1962.
  * [24] Michael L Fredman.  New bounds on the complexity of the shortest path problem.  SIAM Journal on Computing, 5(1):83–89, 1976.
  * [25] Yuzhou Gu, Adam Polak, Virginia Vassilevska Williams, and Yinzhan Xu.  Faster monotone min-plus product, range mode, and single source replacement paths.  arXiv preprint arXiv:2105.02806, 2021.
  * [26] Yuri Gurevich and Saharon Shelah.  Expected computation time for hamiltonian path problem.  SIAM Journal on Computing, 16(3):486–502, 1987.
  * [27] Yijie Han.  Improved algorithm for all pairs shortest paths.  Information Processing Letters, 91(5):245–250, 2004.
  * [28] Yijie Han.  An o (n 3 (log log n/log n) 5/4) time algorithm for all pairs shortest path.  Algorithmica, 51:428–434, 2008.
  * [29] Yijie Han and Tadao Takaoka.  An o(n3 loglogn/log2n) time algorithm for all pairs shortest paths.  Journal of Discrete Algorithms, 38-41:9–19, 2016.  URL: https://www.sciencedirect.com/science/article/pii/S1570866716300296, doi:10.1016/j.jda.2016.09.001.
  * [30] Sophia Heimann, Hung P Hoang, and Stefan Hougardy.  The $k$-opt algorithm for the traveling salesman problem has exponential running time for $k\geq 5$.  arXiv preprint arXiv:2402.07061, 2024.
  * [31] Michael Held and Richard M Karp.  A dynamic programming approach to sequencing problems.  Journal of the Society for Industrial and Applied mathematics, 10(1):196–210, 1962.
  * [32] Anna R Karlin, Nathan Klein, and Shayan Oveis Gharan.  A (slightly) improved approximation algorithm for metric tsp.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 32–45, 2021.
  * [33] Mikko Koivisto and Pekka Parviainen.  A space–time tradeoff for permutation problems.  In Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms, pages 484–492. SIAM, 2010.
  * [34] Ioannis Koutis and Ryan Williams.  Algebraic fingerprints for faster algorithms.  Communications of the ACM, 59(1):98–105, 2015.
  * [35] Mark W Krentel.  Structure in locally optimal solutions.  In 30th Annual Symposium on Foundations of Computer Science, pages 216–221. IEEE Computer Society, 1989.
  * [36] Joseph SB Mitchell.  Guillotine subdivisions approximate polygonal subdivisions: A simple polynomial-time approximation scheme for geometric tsp, k-mst, and related problems.  SIAM Journal on computing, 28(4):1298–1309, 1999.
  * [37] Jesper Nederlof.  Bipartite tsp in $o(1.9999^{n})$ time, assuming quadratic time matrix multiplication.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 40–53, 2020.
  * [38] Alexander Schrijver.  On the history of combinatorial optimization (till 1960).  Handbooks in operations research and management science, 12:1–68, 2005.
  * [39] Anatoliy I Serdyukov.  On some extremal tours in graphs.  Upravlyaemye systemy, 17:76–79, 1978.
  * [40] Ola Svensson, Jakub Tarnawski, and László A. Végh.  A constant-factor approximation algorithm for the asymmetric traveling salesman problem.  J. ACM, 67(6), nov 2020.  doi:10.1145/3424306.
  * [41] Tadao Takaoka.  A new upper bound on the complexity of the all pairs shortest path problem.  Information Processing Letters, 43(4):195–199, 1992.
  * [42] Tadao Takaoka.  A faster algorithm for the all-pairs shortest path problem and its application.  In Computing and Combinatorics: 10th Annual International Conference, COCOON 2004, Jeju Island, Korea, August 17-20, 2004. Proceedings 10, pages 278–289. Springer, 2004.
  * [43] Tadao Takaoka.  An o (n3loglogn/logn) time algorithm for the all-pairs shortest path problem.  Information Processing Letters, 96(5):155–161, 2005.
  * [44] Vera Traub and Jens Vygen.  Approaching 3/2 for the s-t-path tsp.  Journal of the ACM (JACM), 66(2):1–17, 2019.
  * [45] Stephen Warshall.  A theorem on boolean matrices.  Journal of the ACM (JACM), 9(1):11–12, 1962.
  * [46] Ryan Williams.  Faster all-pairs shortest paths via circuit complexity.  In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pages 664–673, 2014.
  * [47] Virginia Vassilevska Williams.  On some fine-grained questions in algorithms and complexity.  In Proceedings of the international congress of mathematicians: Rio de janeiro 2018, pages 3447–3487. World Scientific, 2018.
  * [48] Virginia Vassilevska Williams and Yinzhan Xu.  Truly subcubic min-plus product for less structured matrices, with applications.  In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 12–29. SIAM, 2020.
  * [49] Rico Zenklusen.  A 1.5-approximation for path tsp.  In Proceedings of the thirtieth annual ACM-SIAM symposium on discrete algorithms, pages 1539–1549. SIAM, 2019.
  * [50] Uri Zwick.  A slightly improved sub-cubic algorithm for the all pairs shortest paths problem with real edge lengths.  In International Symposium on Algorithms and Computation, pages 921–932. Springer, 2004.


* García Pérez et al. (2016) García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144, doi: 10.3847/0004-6256/151/6/144
  * Gratton et al. (2004) Gratton, R., Sneden, C., & Carretta, E. 2004, ARA&A, 42, 385, doi: 10.1146/annurev.astro.42.053102.133945
  * Gratton et al. (2012) Gratton, R. G., Carretta, E., & Bragaglia, A. 2012, A&A Rev., 20, 50, doi: 10.1007/s00159-012-0050-3
  * Grevesse et al. (2007) Grevesse, N., Asplund, M., & Sauval, A. J. 2007, Space Sci. Rev., 130, 105, doi: 10.1007/s11214-007-9173-7
  * Griffith et al. (2019) Griffith, E., Johnson, J. A., & Weinberg, D. H. 2019, ApJ, 886, 84, doi: 10.3847/1538-4357/ab4b5d
  * Griffith et al. (2021) Griffith, E., Weinberg, D. H., Johnson, J. A., et al. 2021, ApJ, 909, 77, doi: 10.3847/1538-4357/abd6be
  * Griffith et al. (2023) Griffith, E. J., Hogg, D. W., Dalcanton, J. J., et al. 2023, arXiv e-prints, arXiv:2307.05691, doi: 10.48550/arXiv.2307.05691
  * Griffith et al. (2022) Griffith, E. J., Weinberg, D. H., Buder, S., et al. 2022, ApJ, 931, 23, doi: 10.3847/1538-4357/ac5826
  * Gunn et al. (2006) Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332, doi: 10.1086/500975
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Hansen et al. (2018) Hansen, C. J., El-Souri, M., Monaco, L., et al. 2018, ApJ, 855, 83, doi: 10.3847/1538-4357/aa978f
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harris (1996) Harris, W. E. 1996, AJ, 112, 1487, doi: 10.1086/118116
  * Hasselquist et al. (2016) Hasselquist, S., Shetrone, M., Cunha, K., et al. 2016, ApJ, 833, 81, doi: 10.3847/1538-4357/833/1/81
  * Hasselquist et al. (2017) Hasselquist, S., Shetrone, M., Smith, V., et al. 2017, ApJ, 845, 162, doi: 10.3847/1538-4357/aa7ddc
  * Hasselquist et al. (2021) Hasselquist, S., Hayes, C. R., Lian, J., et al. 2021, ApJ, 923, 172, doi: 10.3847/1538-4357/ac25f9
  * Hayes et al. (2022) Hayes, C. R., Masseron, T., Sobeck, J., et al. 2022, ApJS, 262, 34, doi: 10.3847/1538-4365/ac839f
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85, doi: 10.1038/s41586-018-0625-x
  * Holtzman et al. (2015) Holtzman, J. A., Shetrone, M., Johnson, J. A., et al. 2015, AJ, 150, 148, doi: 10.1088/0004-6256/150/5/148
  * Hubeny & Lanz (2011) Hubeny, I., & Lanz, T. 2011, Synspec: General Spectrum Synthesis Program, Astrophysics Source Code Library, record ascl:1109.022.  http://ascl.net/1109.022
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iben (1965) Iben, Icko, J. 1965, ApJ, 142, 1447, doi: 10.1086/148429
  * Johnson & Pilachowski (2010) Johnson, C. I., & Pilachowski, C. A. 2010, ApJ, 722, 1373, doi: 10.1088/0004-637X/722/2/1373
  * Jönsson et al. (2020) Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120, doi: 10.3847/1538-3881/aba592
  * Kirby et al. (2019) Kirby, E. N., Xie, J. L., Guo, R., et al. 2019, ApJ, 881, 45, doi: 10.3847/1538-4357/ab2c02
  * Kollmeier et al. (2017) Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv e-prints, arXiv:1711.03234, doi: 10.48550/arXiv.1711.03234
  * Leung & Bovy (2019a) Leung, H. W., & Bovy, J. 2019a, MNRAS, 483, 3255, doi: 10.1093/mnras/sty3217
  * Leung & Bovy (2019b) —. 2019b, MNRAS, 489, 2079, doi: 10.1093/mnras/stz2245
  * Liu et al. (2016) Liu, F., Yong, D., Asplund, M., Ramírez, I., & Meléndez, J. 2016, MNRAS, 457, 3934, doi: 10.1093/mnras/stw247
  * Loaiza-Tacuri et al. (2023) Loaiza-Tacuri, V., Cunha, K., Souto, D., et al. 2023, arXiv e-prints, arXiv:2309.10778, doi: 10.48550/arXiv.2309.10778
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Martell et al. (2016) Martell, S. L., Shetrone, M. D., Lucatello, S., et al. 2016, ApJ, 825, 146, doi: 10.3847/0004-637X/825/2/146
  * Masseron et al. (2016) Masseron, T., Merle, T., & Hawkins, K. 2016, BACCHUS: Brussels Automatic Code for Characterizing High accUracy Spectra, Astrophysics Source Code Library, record ascl:1605.004.  http://ascl.net/1605.004
  * McWilliam et al. (2018) McWilliam, A., Piro, A. L., Badenes, C., & Bravo, E. 2018, ApJ, 857, 97, doi: 10.3847/1538-4357/aab772
  * McWilliam et al. (2013) McWilliam, A., Wallerstein, G., & Mottini, M. 2013, ApJ, 778, 149, doi: 10.1088/0004-637X/778/2/149
  * Mészáros et al. (2012) Mészáros, S., Allende Prieto, C., Edvardsson, B., et al. 2012, AJ, 144, 120, doi: 10.1088/0004-6256/144/4/120
  * Mészáros et al. (2020) Mészáros, S., Masseron, T., García-Hernández, D. A., et al. 2020, MNRAS, 492, 1641, doi: 10.1093/mnras/stz3496
  * Mészáros et al. (2021) Mészáros, S., Masseron, T., Fernández-Trincado, J. G., et al. 2021, MNRAS, 505, 1645, doi: 10.1093/mnras/stab1208
  * Mishenina et al. (2015) Mishenina, T., Pignatari, M., Carraro, G., et al. 2015, MNRAS, 446, 3651, doi: 10.1093/mnras/stu2337
  * Myers et al. (2022) Myers, N., Donor, J., Spoo, T., et al. 2022, AJ, 164, 85, doi: 10.3847/1538-3881/ac7ce5
  * Ness et al. (2019) Ness, M. K., Johnston, K. V., Blancato, K., et al. 2019, ApJ, 883, 177, doi: 10.3847/1538-4357/ab3e3c
  * Ness et al. (2022) Ness, M. K., Wheeler, A. J., McKinnon, K., et al. 2022, ApJ, 926, 144, doi: 10.3847/1538-4357/ac4754
  * Nidever et al. (2015) Nidever, D. L., Holtzman, J. A., Allende Prieto, C., et al. 2015, AJ, 150, 173, doi: 10.1088/0004-6256/150/6/173
  * Nidever et al. (2020) Nidever, D. L., Hasselquist, S., Hayes, C. R., et al. 2020, ApJ, 895, 88, doi: 10.3847/1538-4357/ab7305
  * Nomoto et al. (2006) Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424, doi: 10.1016/j.nuclphysa.2006.05.008
  * Osorio et al. (2020) Osorio, Y., Allende Prieto, C., Hubeny, I., Mészáros, S., & Shetrone, M. 2020, A&A, 637, A80, doi: 10.1051/0004-6361/201937054
  * Pinsonneault et al. (2018) Pinsonneault, M. H., Elsworth, Y. P., Tayar, J., et al. 2018, ApJS, 239, 32, doi: 10.3847/1538-4365/aaebfd
  * Plez (2012) Plez, B. 2012, Turbospectrum: Code for spectral synthesis, Astrophysics Source Code Library, record ascl:1205.004
  * Poovelil et al. (2020) Poovelil, V. J., Zasowski, G., Hasselquist, S., et al. 2020, ApJ, 903, 55, doi: 10.3847/1538-4357/abb93e
  * Price-Whelan et al. (2021) Price-Whelan, A. M., Hogg, D. W., Johnston, K. V., et al. 2021, ApJ, 910, 17, doi: 10.3847/1538-4357/abe1b7
  * Ratcliffe & Ness (2023) Ratcliffe, B. L., & Ness, M. K. 2023, ApJ, 943, 92, doi: 10.3847/1538-4357/aca8a1
  * Rybizki et al. (2017) Rybizki, J., Just, A., & Rix, H.-W. 2017, A&A, 605, A59, doi: 10.1051/0004-6361/201730522
  * Sales-Silva et al. (2022) Sales-Silva, J. V., Daflon, S., Cunha, K., et al. 2022, ApJ, 926, 154, doi: 10.3847/1538-4357/ac4254
  * Santana et al. (2021) Santana, F. A., Beaton, R. L., Covey, K. R., et al. 2021, AJ, 162, 303, doi: 10.3847/1538-3881/ac2cbc
  * Sbordone et al. (2007) Sbordone, L., Bonifacio, P., Buonanno, R., et al. 2007, A&A, 465, 815, doi: 10.1051/0004-6361:20066385
  * Schiavon et al. (2017) Schiavon, R. P., Zamora, O., Carrera, R., et al. 2017, MNRAS, 465, 501, doi: 10.1093/mnras/stw2162
  * Schiavon et al. (2023) Schiavon, R. P., Phillips, S. G., Myers, N., et al. 2023, MNRAS, doi: 10.1093/mnras/stad3020
  * Shetrone et al. (2015) Shetrone, M., Bizyaev, D., Lawler, J. E., et al. 2015, ApJS, 221, 24, doi: 10.1088/0067-0049/221/2/24
  * Shetrone et al. (2019) Shetrone, M., Tayar, J., Johnson, J. A., et al. 2019, ApJ, 872, 137, doi: 10.3847/1538-4357/aaff66
  * Silva et al. (2006) Silva, G. M. D., Sneden, C., Paulson, D. B., et al. 2006, AJ, 131, 455, doi: 10.1086/497968
  * Skúladóttir et al. (2019) Skúladóttir, Á., Hansen, C. J., Salvadori, S., & Choplin, A. 2019, A&A, 631, A171, doi: 10.1051/0004-6361/201936125
  * Smiljanic et al. (2016) Smiljanic, R., Romano, D., Bragaglia, A., et al. 2016, A&A, 589, A115, doi: 10.1051/0004-6361/201528014
  * Smith et al. (2000) Smith, V. V., Suntzeff, N. B., Cunha, K., et al. 2000, AJ, 119, 1239, doi: 10.1086/301276
  * Smith et al. (2021) Smith, V. V., Bizyaev, D., Cunha, K., et al. 2021, AJ, 161, 254, doi: 10.3847/1538-3881/abefdc
  * Spina et al. (2021) Spina, L., Ting, Y.-S., De Silva, G. M., et al. 2021, MNRAS, 503, 3279, doi: 10.1093/mnras/stab471
  * Ting & Weinberg (2022) Ting, Y.-S., & Weinberg, D. H. 2022, ApJ, 927, 209, doi: 10.3847/1538-4357/ac5023
  * Van der Swaelmen et al. (2013) Van der Swaelmen, M., Hill, V., Primas, F., & Cole, A. A. 2013, A&A, 560, A44, doi: 10.1051/0004-6361/201321109
  * Ventura et al. (2001) Ventura, P., D’Antona, F., Mazzitelli, I., & Gratton, R. 2001, ApJ, 550, L65, doi: 10.1086/319496
  * Vincenzo et al. (2019) Vincenzo, F., Spitoni, E., Calura, F., et al. 2019, MNRAS, 487, L47, doi: 10.1093/mnrasl/slz070
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Weinberg et al. (2019) Weinberg, D. H., Holtzman, J. A., Hasselquist, S., et al. 2019, ApJ, 874, 102, doi: 10.3847/1538-4357/ab07c7
  * Weinberg et al. (2022) Weinberg, D. H., Holtzman, J. A., Johnson, J. A., et al. 2022, ApJS, 260, 32, doi: 10.3847/1538-4365/ac6028
  * Wheeler et al. (2023a) Wheeler, A. J., Abruzzo, M. W., Casey, A. R., & Ness, M. K. 2023a, AJ, 165, 68, doi: 10.3847/1538-3881/acaaad
  * Wheeler et al. (2023b) Wheeler, A. J., Casey, A. R., & Abruzzo, M. W. 2023b, arXiv e-prints, arXiv:2310.19823, doi: 10.48550/arXiv.2310.19823
  * Wilson et al. (2019) Wilson, J. C., Hearty, F. R., Skrutskie, M. F., et al. 2019, PASP, 131, 055001, doi: 10.1088/1538-3873/ab0075
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237
  * Zamora et al. (2015) Zamora, O., García-Hernández, D. A., Allende Prieto, C., et al. 2015, AJ, 149, 181, doi: 10.1088/0004-6256/149/6/181
  * Zasowski et al. (2013) Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81, doi: 10.1088/0004-6256/146/4/81
  * Zasowski et al. (2017) Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198, doi: 10.3847/1538-3881/aa8df9