∎

11institutetext: R. Smiljanic, A. R. da Silva, R. E. Giribaldi
22institutetext: Nicolaus Copernicus Astronomical Center, Polish Academy of
Sciences, ul. Bartycka 18, 00-716, Warsaw, Poland
22email<EMAIL_ADDRESS>
# Detecting weak beryllium lines with CUBES††thanks: The authors acknowledge
support by the National Science Centre, Poland, project 2018/31/B/ST9/01469

Rodolfo Smiljanic  André R. da Silva  Riano E. Giribaldi

(Received: date / Accepted: date)

###### Abstract

Beryllium is a light element with one single stable isotope, 9Be, which is a
pure product of cosmic-ray spallation in the interstellar medium. Beryllium
abundances in late-type stars can be used in studies about evolutionary
mixing, Galactic chemical evolution, planet engulfment, and the formation of
globular clusters. Some of these uses of Be abundances figure among the
science cases of the Cassegrain U-Band Efficient Spectrograph (CUBES), a new
near-UV low- and medium-resolution spectrograph under development for the Very
Large Telescope. Here, we report on a study about beryllium abundances in
extremely metal-poor stars in the context of the phase A of CUBES. Our
motivation is to understand the limits for the detection of weak lines in
extremely metal-poor stars of low Be abundances. We analyze simulated CUBES
observations, performed in medium-resolution mode, based on synthetic spectra
for four mock stars with [Fe/H] $\leq$ $-$3.0. We find that detecting the Be
lines is possible in certain cases, but is very challenging and requires high
signal-to-noise ratio. Depending on the atmospheric parameters of the target
stars, and if signal-to-noise per pixel of about 400 can be achieved, it
should be possible to detect Be abundances between $\log$(Be/H) $-$13.1 and
$-$13.6, with a typical uncertainty of $\pm$ 0.15 dex. Using CUBES, the
required data for such studies can be obtained for stars that are fainter by
two magnitudes with respect to what is possible with current instrumentation.

###### Keywords:

Low-mass stars Population II stars Stellar abundances Stellar spectral lines

††journal: Experimental astronomy

##  1 Introduction

The near ultraviolet (near-UV) spectral region (which we define here as the
region between $\sim$ 3000–4000 Å) contains the only absorption lines of
beryllium that can be studied in the spectra of late-type stars, when using
ground-based instrumentation. These are the lines located at 3130.423 and
3131.067 Å which originate from the 2s 2S1/2 $\rightarrow$ 2p 2P1/2,3/2
transitions in singly-ionized beryllium (Kramida, 2005).

Beryllium has atomic number Z = 4 and a single stable isotope, 9Be.
Interestingly, the Be atoms that we observe in stellar photospheres are a pure
product of cosmic-ray spallation, being formed in the interstellar medium from
the break up of CNO nuclei (Tatischeff and Gabici, 2018). As became clear from
early works such as Rebolo et al. (1988) and Gilmore et al. (1992), in old,
metal-poor stars the relation between the abundance of Be and the stellar
metallicity111When referring to metallicity we mean the iron abundance,
[Fe/H], which is the usual rough indicator of the total metal content of a
star. Abundances in the bracket notation for two elements A and B, i.e. [A/B],
mean the difference between the logarithm of the ratio of the abundances by
number in a star and the logarithm of the same ratio in the Sun: [A/B] =
$\log$[N(A)/N(B)]⋆ \- $\log$[N(A)/N(B)]⊙. The abundances by number are given
in a scale where the number of hydrogen atoms is N(H) = 1012. is linear. Such
a linear relation implies that Be is produced as a primary element, i.e. its
production rate is independent of metallicity (Prantzos, 2012). Primary
elements are those produced directly from H and He. The only way to obtain Be
that behaves as a primary element is to assume that the Galactic cosmic rays
that are important for the Be nucleosynthesis are made of accelerated CNO
nuclei which then collide with material of the interstellar medium (Duncan et
al., 1992, 1997). Such a fundamental observation about the evolution of Be in
the Galaxy helps in the understanding of the sources and sites from where
Galactic cosmic-rays are accelerated (see e.g. Lingenfelter, 2019, and
references therein).

The study of Be abundances in old, metal-poor stars ([Fe/H] $<\,-$1.00)
started with the few objects analyzed by Molaro et al. (1984), Molaro and
Beckman (1984), Rebolo et al. (1988), Gilmore et al. (1991), and Ryan et al.
(1992). Near-UV sensitive spectrographs at 8-10m class-telescopes allowed the
determination of Be abundances in larger samples, of the order of 150 metal-
poor stars down to [Fe/H] $\sim\,-$3.00 (Smiljanic et al., 2009; Tan et al.,
2009; Boesgaard et al., 2011). Several science cases can benefit if the sample
of metal-poor stars analyzed for Be abundances increases and is extended to
include stars of even lower metallicities.

In this context, the Cassegrain U-Band Efficient Spectrograph (CUBES, see
Zanutta et al., CUBES Phase-A design overview, ExA, in press as part of the
CUBES Special Issue) is expected to make an important contribution. CUBES is
planned to be installed at the Very Large Telescope (VLT) of the European
Southern Observatory (ESO), in Cerro Paranal, Chile. The CUBES instrument was
first discussed in Bristow et al. (2014) and Barbuy et al. (2014). Future
giant telescopes, like the Extremely Large Telescope (ELT) of ESO, will be
optimized for observations in the red/infrared. Therefore, it was realized
that a UV optimized spectrograph at the VLT would remain the best option for
studies requiring wavelength coverage bluer than 4000 Å(Pasquini, 2014; Evans
et al., 2016). The current CUBES design envisions a high-efficiency
spectrograph able to operate in low- (R $\sim$ 6000) and medium-resolution (R
$\sim$ 24 000) dedicated to the near-UV spectral region (Zanutta et al.). The
two resolution modes are obtained with the use of two different image slicers
(Calcines et al., Design of the VLT-CUBES image slicers, ExA, in press as part
of the CUBES Special Issue) in combination with first-order transmission
gratings (Zeitner et al., ExA, in press as part of the CUBES Special Issue).
The science case for such an instrument has been summarized before in Evans et
al. (2018) and Ernandes et al. (2020), and is further discussed in several
contributions to this Special Issue.

In this work, we present our simulations of the expected performance of CUBES
regarding the detection of weak Be lines in faint, extremely metal-poor stars.
We build upon, and expand, previous simulations reported in Smiljanic (2014)
by incorporating the recent updates in the CUBES design. This work is divided
as follows. In Section 2, we summarize the science case that motivates the
observation of the Be lines in the spectra of metal-poor stars. In Section 3,
we discuss the experimental data and codes used to generate the simulated
stellar spectra. Section 4 presents the analysis that attempts to determine
the Be abundances of the simulated observations and a discussion about the
limits for the detection of Be lines using CUBES spectra. Finally, in Section
5 we summarize our findings and highlight a few important caveats that one
should keep in mind when interpreting our simulations.

##  2 Beryllium in metal-poor stars

The models of inhomogeneous (or stochastic) chemical enrichment of the
Galactic halo built by Suzuki et al. (1999) and Suzuki and Yoshii (2001)
suggested that Be abundances in the early Galaxy should have a smaller scatter
than the abundances of elements produced by stellar nucleosynthesis (such as
oxygen or iron). This led to the suggestion that Be abundances could be used
as a good cosmic clock for the early Galaxy (Beers et al., 2000; Suzuki and
Yoshii, 2001). To test this, Pasquini et al. (2004, 2007) investigated Be
abundances in turn-off stars of two globular clusters. The abundances were
converted into ages using chemical evolution models computed by Valle et al.
(2002). The comparison of these values with the cluster ages derived from
isochrone fitting seemed to support the idea of Be as a cosmochronometer, even
though the uncertainty of the abundances was high. We expect that CUBES will
allow an increase in the number of globular clusters that can be used for such
analysis (see Giribaldi et al., ExA, in press as part of the CUBES Special
Issue).

The possibility to use Be abundances as a time scale is particularly
attractive for studies of early Galaxy formation and evolution. Recently, Gaia
data (Gaia Collaboration et al., 2016, 2018) revealed strong evidence that the
Milky Way suffered a major merger between 8 to 11 Gyr ago (Belokurov et al.,
2018; Helmi et al., 2018; Haywood et al., 2018). Hints of this event had
already been identified in previous data, e.g. Chiba and Beers (2000), but it
was the Gaia data that revealed the full picture of the system involved in the
merger. It is now becoming apparent that the evolution of Be abundances as a
function of metallicity is different in the distinct systems that built up the
Galactic halo (Molaro et al., 2020; Smiljanic et al., 2021). The first hints
that Be could be used to separate accreted and in-situ stars, in a diagram of
[O/Fe] as a function of log(Be/H), were found by Pasquini et al. (2005) and
later confirmed by Smiljanic et al. (2009) and Tan and Zhao (2011).
Independently on whether or not Be is indeed a good cosmic clock, it does
offers a different “nucleosynthetic axis” to study the early chemical
evolution of the Milky Way.

There has also being a push towards attempting to determine Be abundances in
several extremely metal-poor stars with [Fe/H] $\lesssim~{}-$3.00 (Primas et
al., 2000b, a; Boesgaard and Novicki, 2006; Rich and Boesgaard, 2009; Ito et
al., 2013; Placco et al., 2014; Spite et al., 2019; Smiljanic et al., 2021;
Ernandes et al., 2022). One motivation for this was to investigate the
possibility of primordial production of substantial amounts of 9Be. Standard
primordial (or Big Bang) nucleosynthesis predicts the formation of tiny
amounts of Be, at the level of $\log$(Be/H) $\sim$ $-$17.6 (Coc et al., 2014).
Several alternative scenarios could result in an enhanced primordial Be
abundance of up to $\log$(Be/H) $\sim$ $-$15 (Boyd and Kajino, 1989; Orito et
al., 1997; Pospelov et al., 2008; Kusakabe et al., 2014, 2017). Nevertheless,
the Be abundances currently detected in extremely metal-poor stars are at the
level of $\log$(Be/H) $\sim$ $-$13.5 with a few upper limits down to
$\log$(Be/H) $\sim$ $-$14.2 (Placco et al., 2014; Spite et al., 2019;
Smiljanic et al., 2021). This is still one to two orders of magnitude above
what could be produced by non-standard primordial nucleosynthesis, and thus
not at the level that could help to constrain such models.

It is thus important to expand investigations of Be abundances in metal-poor
stars, both for tracing the early star forming conditions of the Galaxy and
for constraining any possible deviation from standard Big Bang
nucleosynthesis. The extremely metal-poor stars (with [Fe/H] $<$ $-$3.0)
useful for these analyses are mostly located at far distances, at the Galactic
halo, and are thus faint. Moreover, the Be lines for these stars are very
small and therefore spectra of high signal-to-noise ratio (SNR) are needed for
a proper analysis. Essentially, all such stars that could be easily analysed
with current instrumentation have already been observed. For example,
Smiljanic et al. (2021) report on new near-UV spectra for star BPS BS
16968-061, with [Fe/H] = $-$2.98 and $V$ = 13.22 mag, obtained with the
Ultraviolet and Visual Echelle Spectrograph (UVES, Dekker et al., 2000) at the
VLT. A total of 20 hours of observation resulted in SNR $\sim$ 200 per pixel
in the region around the Be lines. Only a new generation of highly efficient
near-UV spectrographs, like CUBES, can help to observe increased samples of
such stars using shorter observing times.

##  3 Simulated spectra

To simulate the CUBES observations of metal-poor stars, we computed synthetic
spectra covering not only the near-UV but also wavelengths up to about 7000 Å.
This is needed because we rely on the redder wavelengths to compute the
stellar magnitudes in the $V$ band. The work to compile and test the atomic
and molecular line list needed for these calculations is described in the
companion paper by Giribaldi et al. (in this Special Issue). Here, we describe
only our literature investigation which aimed to compile the best available
data for the two Be spectral lines, including for the calculation of their
hyperfine structure.

###  3.1 Hyperfine structure of the Be lines

Table 1: Adopted data for the Be resonance lines. wavelength (Å) | $\log(gf)$ |  $\sigma$.$\alpha$
---|---|---
3130.423 |  $-$0.178 | 123.212
3131.067 |  $-$0.479 | 123.212

The simulations discussed here, have been performed with an updated line list
with respect to what has been used in our recent work (Smiljanic et al.,
2021). The recent literature was scanned for modern results of the atomic data
needed to compute the Be spectral lines. In particular, our line list now
includes the hyperfine splitting of the Be lines. We thus take the opportunity
to describe these updates.

We adopt absolute transition frequencies measured in Nörtershäuser et al.
(2015) and Krieger et al. (2017), which have uncertainties smaller by almost
two orders of magnitude in comparison with the previous values (Bollinger et
al., 1985). We note that the experimental fine structure splitting is in very
good agreement with the modern theoretical value (Puchalski and Pachucki,
2015), the difference being negligible for our purposes. The transition
wavelengths, in Å, are listed in Table 1.

For the two Be lines, the critical compilation of transition probabilities by
Fuhr and Wiese (2010) recommends results of the accurate theoretical
calculations by Yan et al. (1998). We investigated more recent results of
lifetimes and oscillator strengths (Tang et al., 2009; Safronova and
Safronova, 2013; Shukla et al., 2020) but found no change that could affect
the determination of stellar chemical abundances. The adopted values of
$\log(gf)$ are given in Table 1.

For the broadening due to collisions with hydrogen atoms, calculations based
on the theory developed by Anstee and O’Mara (1991) and Barklem and O’Mara
(1998) are available. The line broadening cross sections are tabulated with
values computed at a velocity of 10000 m s-1, and the cross-sections are
assumed to vary according to a power law on the velocity of the form
$v^{-\alpha}$. The cross section ($\sigma$) and the $\alpha$ parameter
computed by Barklem and O’Mara (2000) are listed in Table 1.

For the hyperfine (HFS) structure, we adopt results from the theoretical
calculations of Yerokhin (2008). The magnetic dipole (A) and electric
quadrupole (B) HFS constants needed to compute the energy shifts of each
sublevel are listed in Table 2. Other calculations reported in the literature
give slightly different values for some of the constants (Puchalski and
Pachucki, 2009, 2014). We have tested these values and found that the
differences are at a level that causes no effect on the HFS needed for
synthesizing stellar spectra.

The hyperfine line splitting is a result of the interaction between the
electronic angular momentum, J, and the nuclear angular momentum, I, resulting
in a total angular momentum F. The nuclear angular momentum for 9Be is I =
3/2. The selection rules allow transitions between levels with $\Delta$F =
$\pm$1,0 but not between two levels with F = 0. This results in the splitting
of the 3130.423 Å line into six components and of the 3131.067 Å line into
four components. The energy shifts and line strengths were computed with
equations that can be found in texts like Woodgate (1970) and Emery (2006),
and with help of the code provided by Mathar (2011). The wavelengths and
$\log(gf)$ of the hyperfine transitions are given in Table 3.

Table 2: Constants needed to compute the energy shifts of the hyperfine sublevels (Yerokhin, 2008). |  2S${}_{\frac{1}{2}}$ |  2P${}_{\frac{1}{2}}$ |  2P${}_{\frac{3}{2}}$
---|---|---|---
A (MHz) |  $-$625.11 |  $-$117.925 |  $-$1.016
B (MHz) | – | – |  $-$2.281
Table 3: Hyperfine structure of the Be lines. wavelength (Å) | $\log(gf)$ | Fl – Fu
---|---|---
3130.425 |  $-$1.382 | 1 – 0
3130.425 |  $-$0.984 | 1 – 1
3130.425 |  $-$0.984 | 1 – 2
3130.421 |  $-$1.683 | 2 – 1
3130.421 |  $-$0.984 | 2 – 2
3130.421 |  $-$0.537 | 2 – 3
3131.069 |  $-$1.683 | 1 – 1
3131.070 |  $-$0.984 | 1 – 2
3131.065 |  $-$0.984 | 2 – 1
3131.066 |  $-$0.984 | 2 – 2

###  3.2 Simulated observations

Synthetic stellar spectra were computed with the code Turbospectrum (Plez,
2012) using the MARCS one-dimensional model atmospheres (Gustafsson et al.,
2008) under the local thermodynamical equilibrium (LTE) approximation. The
spectra were computed between 3000-7000 Å on steps of 0.02 Å.

In total, spectra of four stars were simulated. We chose to simulate spectra
of two main-sequence stars with ($T_{\mathrm{eff}}$, $\log~{}g$) = (6300 K,
4.30) and metallicities [Fe/H] = $-$3.0 and $-$3.5. These are parameters
similar to the stars analysed in Smiljanic et al. (2021). We also simulated
spectra for two subgiants with ($T_{\mathrm{eff}}$, $\log~{}g$) = (5600 K,
3.40) and the same two metallicity values, [Fe/H] = $-$3.0 and $-$3.5. These
temperature and $\log~{}g$ values are the same of the star analyzed in Spite
et al. (2019). All stars were simulated with microturbulence velocity ($\xi$)
of 1.50 km s-1 and using $\alpha$-enhanced models with [$\alpha$/Fe] = +0.4.
With these choices we approximately bracket the parameter space for the
extremely metal-poor stars that one would select to observe.

In addition, for each of the four stars, 20 spectra with different abundances
of Be were computed where one of these 20 was always simulated without the Be
lines. For the stars with [Fe/H] = $-$3.0, the other 19 spectra have log(Be/H)
that varies between $-$13.8 and $-$12.9 in steps of 0.05 dex (i.e. [Be/Fe]
varying from $-$0.12 to +0.78). For the stars with [Fe/H] = $-$3.5, the 19
spectra have log(Be/H) that varies between $-$14.1 and $-$13.2 in steps of
0.05 dex (i.e. [Be/Fe] varying from +0.08 to +0.98). Thus, our initial set of
simulations includes 80 spectra.

To understand the properties of the observed spectra, we make use of the CUBES
end-to-end simulator (E2E) and the CUBES exposure time calculator (ETC)
described in Genoni et al. (The CUBES Instrument model and simulation tools,
ExA, in press as part of the CUBES Special Issue). The E2E simulates the
effects of the atmosphere, the fore-optics, the image slicers, the
spectrograph, and the detector on a real astronomical observation. According
to the E2E, the current CUBES design provides R $\sim$ 23 000 with a sampling
of $\sim$2.35 px (pixels of 0.058 Å) at 3131 Å, the region of the Be lines.

To estimate the SNR of the observations, we used the ETC assuming observations
for a bright ($V$ = 12.5 mag, in the Vega system) and a faint ($V$ = 14 mag)
case. The SNR was estimated for observations with airmass = 1.2 and with a
seeing of 0.9 arcsec at 5500 Å (which corresponds to about 1.0 arcsec at 3100
Å). The exposure time was set to 3000s. This time was chosen to simulate an
observing block of one hour, as is typical for service mode observations in
Paranal, assuming 10 min of overhead, as is the case for observations with
UVES at the VLT. For the main-sequence star, we obtain SNR per pixel $\sim$
200 and 100, for the bright and faint cases, respectively. For the subgiant,
the values are SNR $\sim$ 170 and 85, for the bright and faint cases,
respectively. The lower SNR of the subgiant, reflects the lower flux that it
outputs in the near-UV wavelengths because of its lower $T_{\mathrm{eff}}$. In
addition to these SNR values, we consider a case of higher SNR that simulates
the gain obtained by stacking four observations of each star (i.e., a total of
four observing hours spent in each target). This option of stacked spectra
have SNR = 400, for the main-sequence stars, and SNR = 340, for the subgiant
stars. For comparison, the observations of the same stars were simulated using
the UVES ETC, assuming similar observing conditions and using the blue arm
centered at 346nm with the 1.0 arcsec slit. For the main sequence star, the
obtained SNR values are 30 and 13, for the bright and faint case,
respectively. For the subgiant, the values are SNR = 18 and 7, for the bright
and faint case, respectively.

Gaussian noise equivalent to these SNR levels were added to the spectra using
the task mknoise in IRAF222Used within AstroConda, which is maintained by the
Space Telescope Science Institute (STScI). See
https://astroconda.readthedocs.io/en/latest/index.html.. In total, four
realizations of the noise were created, so that an indication of the random
uncertainties affecting the analysis could be obtained. Four noise simulations
of 80 spectra in three different SNR levels results in a total of 960 spectra
for our analysis.

Figure 1: Comparison of the stellar spectra, with and without noise, for three
different beryllium abundances. Note the different y-axis scale for the panels
showing the subgiant stars. The Be abundances are given in the insets of each
panel. Top row: Normalized spectra for the main-sequence star (left) and the
subgiant star (right) with [Fe/H] = $-$3.0, without noise. Middle row:
Normalized spectra for the bright (left) and faint (right) main-sequence star
with [Fe/H] = $-$3.0. The length of the error bars indicate the $\pm$1$\sigma$
error for SNR = 200 and 100 for the bright and faint cases, respectively.
Bottom row: Normalized spectra for a bright (left) and faint (right) subgiant
star with [Fe/H] = $-$3.0 (SNR = 170 and 85, for the bright and fain stars,
respectively).

##  4 Detecting the weak Be lines

The magnitude of the challenge that we want to address can be appreciated in
Fig. 1. For the main-sequence star, the Be lines are very weak and hard to be
appropriately measured even in the case without noise. It is clear that there
is a change in the detection limit with SNR, when moving from the bright to
the faint case. For the subgiant star, which is cooler and has lower
$\log~{}g$ value, the Be lines are stronger when compared to the main sequence
star, even if the abundance values are the same. The detection of the lines is
easier for this type of star. Besides the effect of the SNR, it is important
to note from these plots that these weak Be lines are spread over only three
to maybe five pixels in the CUBES spectra.

To determine if the Be lines in a given noised spectrum can be detected, we
proceeded as follows. Each noised observed spectrum of a given star is
compared to the 20 synthetic spectra, without noise, of the same star (where
the Be abundances vary by 0.05 dex with one spectrum computed without Be). The
first step is to search for the synthetic spectrum that best fits the
observations. Here, we simply compute the coefficient of determination,
$R^{2}$, which involves the ratio between the residual sum of squares and the
total sum of squares:

$R^{2}=1-\frac{\sum^{n}_{i=1}\,(O_{i}-M_{i})^{2}}{\sum^{n}_{i=1}\,(O_{i}-\overline{O})^{2}},$
(1)

where $n$ is the total number of pixels used in the fitting, $O_{i}$ is the
normalized flux of the simulated observation in pixel $i$, $M_{i}$ is the
normalized flux of the model synthetic spectrum in pixel $i$, and
$\overline{O}$ is the mean of all $O_{i}$ values. We apply this equation in a
region including 5 pixels around the center of each Be line. A good model that
fits the observations will have $R^{2}$ approaching one (as the residual sum
of squares tends to zero). Model selection using $R^{2}$ is usually not
advised, as models can be built with a large number of free parameters to
provide an ever better fit to the data. In our case, however, the only free
parameter is the Be abundance and we simply want to find the model that better
explains the simulated data. In this procedure, each Be line is analyzed
separately.

After finding the best fitting model, we still test if the Be line can be
considered detected or if an upper limit should be reported instead. This is
done through the comparison of two distributions of normalized residuals. One
distribution from the simulated observation with respect to the best-fitting
synthetic spectrum, the other from the simulated observation with respect to
the spectrum without Be. In both cases, the residuals were normalized to the
error expected from the SNR. These two distributions are then compared using a
two-sample Kolmogorov-Smirnov test. The null hypothesis is that there is no
difference between the two distributions. If the null hypothesis can not be
rejected, then the spectrum without Be gives a fit that is statistically as
good as the one with the selected Be abundance. This means that the line could
not be clearly detected. If the null hypothesis is rejected, we consider that
a detection was obtained.

###  4.1 The subgiant star case

As an example, we start looking at the four simulated one hour observations of
a bright ($V$ = 12.5 mag) subgiant star with [Fe/H] = -3.0 and
$\log(\mathrm{Be/H})$ = $-$13.2. This value of Be abundance was the detection
limit for the main-sequence stars observed with UVES analyzed in Smiljanic et
al. (2021). In Figs. 2 and 3 we show examples of fitting made to the 3130 Å
and 3131 lines, respectively. In all four observations, both lines are
detected with best fitting abundance differing by at most 0.05 dex of the
input one.

For the 3130 Å line, three realizations of the observed spectrum are fit with
Be at the same value of the input abundance, $\log(\mathrm{Be/H})$ = $-$13.2.
The fourth observed spectrum is fit with a slightly lower value,
$\log(\mathrm{Be/H})$ = $-$13.25 (this last case is shown at top panel of Fig.
2. A boxplot of the four abundances is shown in the bottom right panel.).
Looking at the $R^{2}$ values of the fits, however, we note that even if there
is one fit that gives a maximum value close to one, there are other fits, with
abundances within $\pm$0.05-0.10 dex of the best fit, that return very similar
$R^{2}$ values (bottom left panel of Fig. 2). This suggests that, in this
case, the fitting uncertainty is of the order of $\pm$0.10 dex. In the fitting
of all four simulated observations, the Kolmogorov-Smirnov test indicates that
the line was detected.

Figure 2: Fitting of the four simulated observations of a bright ($V$ = 12.5
mag) subgiant star with [Fe/H] = -3.0 and $\log(\mathrm{Be/H})$ = $-$13.2,
using five pixels around the 3130 Å line. Top: One selected fit that returns
$\log$(Be/H) = $-$13.25, slightly lower than the input abundance of the
observation. Bottom left: The $R^{2}$ values obtained when fitting the Be line
with spectra computed with a range of abundances. Note that the spectrum
without Be returns a worse fit than all other options that include the Be
lines. Bottom right: Boxplot of the best fitting abundances in each of the
four simulated observations.

For the 3131 Å line the results are very similar, even if this line is weaker.
Two of the noise realizations are fit with Be at the same value of the input
abundance, $-$13.2, but two are fit by a slightly higher value, $-$13.15 (one
of them is shown at top panel of Fig. 3. The boxplot of the four abundances is
shown in the bottom right panel). The $R^{2}$ values here span a wider range.
This happens in part because of the term at the bottom part of Equation 1,
which involves the average of the observed flux values. For the stronger 3130
line, which is also affected by strong blends, this term is larger which has
the effect of reducing the scale of the total $R^{2}$ variation. The $R^{2}$
values indicate that abundances within $\pm$0.05 dex of the best value give
results of similar quality (bottom left panel of Fig. 3). In the fitting of
all four simulated observations, the Kolmogorov-Smirnov test indicates that
the line was detected.

Figure 3: Results of fitting five pixels around the 3131 Å line in the same
observations shown in Fig. 2. Top: One selected fit that returns $\log$(Be/H)
= $-$13.15 while the input abundance of the observation was $-$13.20. Bottom
left: The $R^{2}$ values obtained when fitting the line with spectra of a
range of abundances. The spectrum without Be again returns a worse fit than
all other options. Bottom right: Boxplot of the best fitting abundances in
each of the four simulated observations.

To determine the detection limit, we analyse the other simulated observations
in the same way, looking at the ones with progressively lower Be abundance, in
steps of 0.05 dex. For the weaker line at 3131 Å, we reach a quite different
situation already when the input abundance of the observations is $\log$(Be/H)
= $-$13.25. The best fitting synthetic spectra for the four observations agree
with the input within $\pm$0.05, the $R^{2}$ values indicate a clear peak
close to unity around the same values, but in three cases the Kolmogorov-
Smirnov test indicates that the line was not detected. Exactly the same
happens for $\log$(Be/H) = $-$13.3 and $-$13.35. For the next step, $-$13.4,
all four Kolmogorov-Smirnov tests indicate that a detection was not possible.
Moreover, the $R^{2}$ values now fail to approach unit, remaining between
08-0.9 for three cases, and with one case where the maximum value is below
zero. While this last abundance level seems clearly not detected, the case for
the other values seem more ambiguous. We conclude that the limit for a
detection of the 3131Å Be line in this case is at an abundance of $\log$(Be/H)
= $-$13.25.

Figure 4: Same as Fig. 2, fitting the 3130 Å line for an input abundance of
$\log(\mathrm{Be/H})$ = $-$13.6. Top: One selected best fit with $\log$(Be/H)
= $-$13.6, the same value of the input abundance of the observation. Bottom
left: The $R^{2}$ values obtained when fitting the Be line with spectra with a
range of abundances. Note that the range is different from the previous plots.
Several of the fitted spectra provide $R^{2}$ values above 0.8. Bottom right:
Boxplot of the best fitting abundances in each of the four simulated
observations.

The line at 3130 Å is stronger, so visually it seems that it could be detected
at lower values of the abundance. However, already at the same value of
$\log$(Be/H) = $-$13.25 some ambiguity starts to appear. The best fitting
abundances are within $\pm$0.05 dex of the input, the $R^{2}$ plots indicate
that abundances within $\pm$0.10 of the best value provide similar fits, but
one of the four Kolmogorov-Smirnov tests indicate that a detection was not
possible. Essentially the same happens for observations with $\log$(Be/H) =
$-$13.3 and $-$13.35, although now the $R^{2}$ plots indicate some degeneracy
in the fits of maybe $\pm$0.15. Only for abundances below $\log$(Be/H) =
$-$13.4, two or more Kolmogorov-Smirnov tests indicate that a detection was
not possible. From these considerations, we conclude that the limit for a
clear detection of the 3130Å Be line is most likely around an abundance of
$\log$(Be/H) = $-$13.35, which is the last case before the analysis starts to
return mostly dubious fits.

If we increase the SNR of the observations, for a case of four exposures of
the bright subgiant (which is also equivalent to 16 exposures of the faint
subgiant), the limits improve. For the 3130 Å line, most tests indicate that a
detection is possible down to $\log$(Be/H) = $-$13.6 (see Fig. 4). At this
point, however, the $R^{2}$ plot shows that there is already no difference
between spectra with abundances $\pm$0.15 dex of the best fit. For the weaker
Be line, the tests indicate we can detect the line down to $\log$(Be/H) =
$-$13.5. Finally, for the case of lower SNR (one single exposure of the faint
subgiant), our tests indicate a detection only for the highest level of the
abundances that we simulated, $\log$(Be/H) = $-$12.95.

For the case of a subgiant of lower metallicity, [Fe/H] = $-$3.5, we only
simulated abundances below $\log$(Be/H) = $-$13.2. At this level, for the
cases of SNR = 85, no clear detection seems possible. At SNR = 170, line 3130
Å seems to be detectable at lower Be abundances when compared to the case with
[Fe/H] = $-$3.0. On the other hand, line 3131 Å seems detectable only at a
slightly higher Be abundance. For the line at 3130 Å, this difference might
have to do with the reduced influence of the blends surrounding it. For the
other line, it seems just part of the statistical fluctuation. In the same
way, for the case with SNR = 340, the line at 3130 Å seems to be detectable
down to a lower abundance, $\log$(Be/H) = $-$13.7, with little difference
among the fits with $\pm$0.15 dex.

$V$ (mag) | texp | SNR | [Fe/H] | Line 3130 Å | Line 3131Å
---|---|---|---|---|---
14 | 3000s | 85 |  $-$3.00 |  $\geq$ $-$12.95 |  $\geq$ $-$12.95
14 | 3000s | 85 |  $-$3.50 | – | –
12.5 | 3000s | 170 |  $-$3.00 |  $\geq$ $-$13.35 |  $\geq$ $-$13.25
12.5 | 3000s | 170 |  $-$3.50 |  $\geq$ $-$13.50 |  $\geq$ $-$13.20
12.5 | 4 $\times$ 3000s | 340 |  $-$3.00 |  $\geq$ $-$13.60 |  $\geq$ $-$13.50
12.5 | 4 $\times$ 3000s | 340 |  $-$3.50 |  $\geq$ $-$13.70 |  $\geq$ $-$13.45
Table 4: Limits of beryllium abundance detection for the subgiant case. Note
that the case of SNR = 170 can also be obtained with four 3000s exposures for
the faint, $V$ = 14 mag, star. The case with SNR = 340 would require 16
exposures of 3000s for the same faint star.

Overall, the conclusion here seems to be that, as long as the SNR approaches a
value close to 400, with CUBES we can push the detection limit to $\log$(Be/H)
= $-$13.6 or $-$13.7, with an uncertainty of the order of $\pm$0.15 dex. The
derived detection limits for each level of SNR are given in Table 4.

###  4.2 The main-sequence star case

Finally, let us examine the case for the observations of a main sequence star.
For this case, the Be lines are weaker, so we can anticipate detection limits
at higher abundances. With SNR = 100, neither of the two lines can be
detected, even at the maximum Be abundances that we simulated ($\log$(Be/H) =
$-$12.9 for [Fe/H] = $-$3.0 and $\log$(Be/H) = $-$13.2 for [Fe/H] = $-$3.5).

Even increasing the SNR to 200, the line at 3131 Å remains not detectable. It
becomes marginally possible to detect this line only with SNR = 400 and only
at the highest abundances simulated for the [Fe/H] = $-$3.0 case. For
$\log$(Be/H) = $-$12.9, two of the fits pass the Kolmogorov-Smirnov test and
the $R^{2}$ values are above 0.8.

It seems that it is only the stronger line at 3130 Å that stands a chance of
really being detected. At the same SNR = 200, we determined that the limit is
at $\log$(Be/H) = $-$13.0 (Fig. 5). The abundances of the best fitted spectra
are within $\pm$0.05 dex of the input value, there is a clear peak in the
$R^{2}$ values, and the Kolmogorov-Smirnov test is still consistent with a
detection. For lower abundances in the observed spectra, the fits are still
good, both visually and in the $R^{2}$ plot. However, the Kolmogorov-Smirnov
tests already indicate that the noise level precludes a clear detection. For
the main-sequence star of lower metallicity, [Fe/H] = -3.5, our simulation
with the highest Be abundance, $\log$(Be/H) = $-$13.2, is not detectable.

The final case is where the stacked observations of the main-sequence star
result in SNR = 400. As was the case for the subgiant star, the detection
limit seems to depend on how strong are the surrounding blends. For a star
with [Fe/H] = $-$3.00, the lower abundance that we can detect is $\log$(Be/H)
= $-$13.2 (Fig. 6). In this case, the $R^{2}$ plot indicates that synthetic
spectra with abundances within $\pm$0.10 dex return fits of similar quality.
For a star with [Fe/H] = $-$3.50, the lower abundance we can detect is
smaller, $\log$(Be/H) = $-$13.4 (Fig. 7). In this case, the boxplot shows that
the best fits were obtained with abundances within $\pm$0.15 dex of the input
value. The $R^{2}$ plots indicates that synthetic spectra with abundances
within $\pm$0.10 dex return fits of similar quality. For the main-sequence
star case, the detection limits for each level of SNR are given in Table 5.

Figure 5: Same as Fig. 2, but fitting the 3130 Å line in the case of a bright
main-sequence star with [Fe/H] = $-$3.0 and input abundance of
$\log(\mathrm{Be/H})$ = $-$13.0, observed for 3000s. Top: One selected best
fit with $\log$(Be/H) = $-$13.0, the same value of the input abundance of the
observation. Bottom left: The $R^{2}$ values obtained when fitting the Be line
with spectra with a range of abundances. Bottom right: Boxplot of the best
fitting abundances in each of the four simulated observations.

Figure 6: Same as Fig. 5, but for the case where SNR = 400. Top: One selected
best fit with $\log$(Be/H) = $-$13.15 for the observation that was simulated
with input abundance of $\log$(Be/H) = $-$13.20. Bottom left: The $R^{2}$
values obtained when fitting the Be line with spectra with a range of
abundances. Bottom right: Boxplot of the best fitting abundances in each of
the four simulated observations.

Figure 7: Same as Fig. 6, but for the case of star with [Fe/H] = $-$3.5. Top: One selected best fit with $\log$(Be/H) = $-$13.35 for the observation that was simulated with input abundance of $\log$(Be/H) = $-$13.40. Bottom left: The $R^{2}$ values obtained when fitting the Be line with spectra with a range of abundances. Bottom right: Boxplot of the best fitting abundances in each of the four simulated observations. $V$ (mag) | texp | SNR | [Fe/H] | Line 3130 Å | Line 3131Å
---|---|---|---|---|---
14 | 3000s | 100 |  $-$3.00 | – | –
14 | 3000s | 100 |  $-$3.50 | – | –
12.5 | 3000s | 200 |  $-$3.00 |  $\geq$ $-$13.0 | –
12.5 | 3000s | 200 |  $-$3.50 | – | –
12.5 | 4 $\times$ 3000s | 400 |  $-$3.00 |  $\geq$ $-$13.2 |  $\geq$ $-$12.90
12.5 | 4 $\times$ 3000s | 400 |  $-$3.50 |  $\geq$ $-$13.4 | –
Table 5: Limits of beryllium abundance detection for the mains-sequence star
case. Note that the case of SNR = 200 can also be obtained with four 3000s
exposures for the faint star. The case with SNR = 400 would require 16
exposures of 3000s for the faint star.

##  5 Conclusions

In the context of the phase A study of CUBES, we investigated the limits for
detecting weak Be lines in the spectra of extremely metal-poor stars ([Fe/H]
$\leq$ $-$3.0). The current design of CUBES, described in Zanutta et al. (this
Special Issue), shows that it should be possible to obtain spectra with R
$\sim$ 23 000 and a sampling of 2.35 pixels in the region of the Be lines.

For the simulations, we estimated the SNR that can be obtained using the CUBES
ETC (Genoni et al., this Special Issue). We considered a standard service mode
observing block of one hour with ten minutes of overhead and 50 minutes of
exposure time. To compute the spectra, we adopted the line list described in
the companion paper of Giribaldi et al. (this Special Issue) with the
inclusion of the HFS splitting of the Be lines. We simulated the spectra of
two stars, one subgiant with parameters similar to the star analysed by Spite
et al. (2019) and one main-sequence star with parameters similar to the sample
analysed in Smiljanic et al. (2021).

The simulations indicate that, as long as a SNR of about 400 can be reached,
it should be possible to detect the stronger Be line (at 3130 Å) with
abundances of $\log$(Be/H) $\geq$ $-$13.4 dex and of $\log$(Be/H) $\geq$
$-$13.7 dex, for the main-sequence and subgiant cases, respectively. The
typical uncertainty of the fitting only can be of the order of $\pm$0.10-0.15
dex. The SNR $\sim$ 400 can be reached in 4 hours for a star of $V$ = 12.5 mag
but would still require about 16 hours for stars $V$ $\sim$ 14 mag. In Fig. 8
we indicate the region where CUBES can potentially detect Be abundances in
comparison to the results from Smiljanic et al. (2021).

Figure 8: Relation between the Be abundance and metallicity. The upper limits
for stars 2MASS J18082002-5104378 (Spite et al., 2019) and BD+44 493 (Placco
et al., 2014) are shown as upside-down red and orange triangles, respectively.
The blue circles and blue triangles (for upper limits) are stars analyzed in
Smiljanic et al. (2021). The black circles are stars from Smiljanic et al.
(2009). The solid line is a fit to this last data set. The gray box indicates
the parameter space region with $\log$(Be/H) $\geq$ $-$13.7 for stars with
$-$2.5 $\geq$ [Fe/H] $\geq$ $-$3.5 where CUBES can potentially be used to
detect Be abundances in extremely metal-poor stars and improve the current
statistics.

A search in the PASTEL catalog (Soubiran et al., 2016) and the SAGA database
(Suda et al., 2008) returned 13 known stars with $V$ $\leq$ 14 mag, $\log~{}g$
$\geq$ 3.4, and [Fe/H] $\leq$ $-$3.0, that could be interesting targets for a
Be survey in extremely metal-poor stars (meaning that Be abundances have not
been determined for them before). The number of targets can grow to $\sim$ 21
known objects, if we can push the magnitude limit down to $V$ $\sim$ 14.3 mag.
Populating the region with [Fe/H] $\leq$ $-$3.00 of the $\log$(Be/H) vs.
[Fe/H] diagram can help in understanding how inhomogeneous were the early
stages of the Galactic chemical enrichment (see Smiljanic et al., 2021, for a
discussion). Our results suggest that the lower resolution of the CUBES
spectra limits its capability to improve the very low upper limits around or
below $\log$(Be/H) $\sim$ $-$14, already existing for a few stars in the
literature (e.g., Placco et al., 2014; Spite et al., 2019). The real power of
CUBES will be in increasing the magnitude limit of the sample that can be
potentially analysed down to $\log$(Be/H) $\geq$ $-$13.7.

One should note, however, that our simulations assume a very idealized case
for the analysis. We know a priory, and know perfectly well, the values for
the atmospheric parameters of the stars, for the abundances of all other
elements in the spectrum, and for the total broadening affecting the spectra.
In a real case, all these quantities are unknowns that will be estimated with
some uncertainty. Those uncertainties will affect the final result of the
analysis. Moreover, the only source of noise affecting our simulated
observations is one that is perfectly Gaussian. Finally, deviations from the
local thermodynamic equilibrium in the formation of the Be lines were not
taken into account. Recent results, however, suggest that such effects might
be important in metal-poor stars (Korotin and Kučinskas, 2022). We plan to
improve the simulations, addressing some of these points, in the next phases
of the development of CUBES. It is clear that detecting the weak Be lines will
be a challenging endeavour, but it is also clear that this science case can
only be addressed with high-efficiency new spectrographs like CUBES.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

The line list used to compute our synthetic spectra is available in Github,
https://github.com/RGiribaldi/Master-line-list-for-spectral-synthesis-with-
Turbospectrum. The other datasets generated during and/or analysed during the
current study are available from the corresponding author on reasonable
request.

## References

  * Anstee and O’Mara (1991) Anstee SD, O’Mara BJ (1991) An investigation of Brueckner’s theory of line broadening with application to the sodium D lines. MNRAS, 253:549–560, DOI 10.1093/mnras/253.3.549
  * Barbuy et al. (2014) Barbuy B, Bawden Macanhan V, Bristow P, Castilho B, Dekker H, Delabre B, Diaz M, Gneiding C, Kerber F, Kuntschner H, La Mura G, Maciel W, Meléndez J, Pasquini L, Pereira CB, Petitjean P, Reiss R, Siqueira-Mello C, Smiljanic R, Vernet J (2014) CUBES: cassegrain U-band Brazil-ESO spectrograph. Ap&SS, 354(1):191–204, DOI 10.1007/s10509-014-2039-z
  * Barklem and O’Mara (1998) Barklem PS, O’Mara BJ (1998) The broadening of strong lines of Ca+̂, Mg+̂ and Ba+̂ by collisions with neutral hydrogen atoms. MNRAS, 300(3):863–871, DOI 10.1046/j.1365-8711.1998.01942.x
  * Barklem and O’Mara (2000) Barklem PS, O’Mara BJ (2000) Broadening of lines of Beii, Srii and Baii by collisions with hydrogen atoms and the solar abundance of strontium. MNRAS, 311(3):535–540, DOI 10.1046/j.1365-8711.2000.03090.x
  * Beers et al. (2000) Beers TC, Suzuki TK, Yoshii Y (2000) The Light Elements Be and B as Stellar Chronometers in the Early Galaxy. In: da Silva L, de Medeiros R, Spite M (eds) The Light Elements and their Evolution, vol 198, p 425, astro-ph/0002056
  * Belokurov et al. (2018) Belokurov V, Erkal D, Evans NW, Koposov SE, Deason AJ (2018) Co-formation of the disc and the stellar halo. MNRAS, 478(1):611–619, DOI 10.1093/mnras/sty982, 1802.03414
  * Boesgaard and Novicki (2006) Boesgaard AM, Novicki MC (2006) Beryllium in Disk and Halo Stars: Evidence for a Beryllium Dispersion in Old Stars. ApJ, 641(2):1122–1130, DOI 10.1086/500501, astro-ph/0512317
  * Boesgaard et al. (2011) Boesgaard AM, Rich JA, Levesque EM, Bowler BP (2011) Beryllium and Alpha-element Abundances in a Large Sample of Metal-poor Stars. ApJ, 743(2):140, DOI 10.1088/0004-637X/743/2/140, 1110.2823
  * Bollinger et al. (1985) Bollinger JJ, Wells JS, Wineland DJ, Itano WM (1985) Hyperfine structure of the 2p sup2Psub1/2 state in sup9Besup+. Physical Review A 31(4):2711–2714, DOI 10.1103/PhysRevA.31.2711
  * Boyd and Kajino (1989) Boyd RN, Kajino T (1989) Can 9Be Provide a Test of Cosmological Theories? ApJ Letters, 336:L55, DOI 10.1086/185360
  * Bristow et al. (2014) Bristow P, Barbuy B, Macanhan VB, Castilho B, Dekker H, Delabre B, Diaz M, Gneiding C, Kerber F, Kuntschner H, La Mura G, Reiss R, Vernet J (2014) Introducing CUBES: the Cassegrain U-band Brazil-ESO spectrograph. In: Ramsay SK, McLean IS, Takami H (eds) Ground-based and Airborne Instrumentation for Astronomy V, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 9147, p 914709, DOI 10.1117/12.2054751
  * Chiba and Beers (2000) Chiba M, Beers TC (2000) Kinematics of Metal-poor Stars in the Galaxy. III. Formation of the Stellar Halo and Thick Disk as Revealed from a Large Sample of Nonkinematically Selected Stars. AJ, 119(6):2843–2865, DOI 10.1086/301409, astro-ph/0003087
  * Coc et al. (2014) Coc A, Uzan JP, Vangioni E (2014) Standard big bang nucleosynthesis and primordial CNO abundances after Planck. J. of Cosmol. and Astropart. Phys., 2014(10):050, DOI 10.1088/1475-7516/2014/10/050, 1403.6694
  * Dekker et al. (2000) Dekker H, D’Odorico S, Kaufer A, Delabre B, Kotzlowski H (2000) Design, construction, and performance of UVES, the echelle spectrograph for the UT2 Kueyen Telescope at the ESO Paranal Observatory. In: Iye M, Moorwood AF (eds) Optical and IR Telescope Instrumentation and Detectors, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 4008, pp 534–545, DOI 10.1117/12.395512
  * Duncan et al. (1992) Duncan DK, Lambert DL, Lemke M (1992) The Abundance of Boron in Three Halo Stars. ApJ, 401:584, DOI 10.1086/172088
  * Duncan et al. (1997) Duncan DK, Primas F, Rebull LM, Boesgaard AM, Deliyannis CP, Hobbs LM, King JR, Ryan SG (1997) The Evolution of Galactic Boron and the Production Site of the Light Elements. ApJ, 488(1):338–349, DOI 10.1086/304683
  * Emery (2006) Emery G (2006) Hyperfine Structure in Springer handbook of atomic, molecular, and optical physics, Springer Science+Business Media, Inc., p 253\. DOI 10.1007/978-0-387-26308-3“˙16
  * Ernandes et al. (2020) Ernandes H, Evans CJ, Barbuy B, Castilho B, Cescutti G, Christlieb N, Cristiani S, Di Marcantonio P, Hansen C, Quirrenbach A, Smiljanic R (2020) Stellar astrophysics in the near-UV with VLT-CUBES. In: Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 11447, p 1144760, DOI 10.1117/12.2562497, 2102.02205
  * Ernandes et al. (2022) Ernandes H, Barbuy B, Friaça A, Hill V, Spite M, Spite F, Castilho BV, Evans CJ (2022) Be, V, and Cu in the halo star CS 31082-001 from near-UV spectroscopy. MNRAS, 510(4):5362–5375, DOI 10.1093/mnras/stab3789, 2202.04450
  * Evans et al. (2016) Evans CJ, Puech M, Rodrigues M, Barbuy B, Cuby JG, Dalton G, Fitzsimons E, Hammer F, Jagourel P, Kaper L, Morris SL, Morris TJ (2016) Science requirements and trade-offs for the MOSAIC instrument for the European ELT. In: Evans CJ, Simard L, Takami H (eds) Ground-based and Airborne Instrumentation for Astronomy VI, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 9908, p 99089J, DOI 10.1117/12.2231675, 1608.06542
  * Evans et al. (2018) Evans CJ, Barbuy B, Castilho B, Smiljanic R, Melendez J, Japelj J, Cristiani S, Snodgrass C, Bonifacio P, Puech M, Quirrenbach A (2018) Revisiting the science case for near-UV spectroscopy with the VLT. In: Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 10702, p 107022E, DOI 10.1117/12.2312022, 1806.11173
  * Fuhr and Wiese (2010) Fuhr JR, Wiese WL (2010) Tables of Atomic Transition Probabilities for Beryllium and Boron. Journal of Physical and Chemical Reference Data 39(1):013101–013101, DOI 10.1063/1.3286088
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti T, de Bruijne JHJ, Brown AGA, Vallenari A, et al (2016) The Gaia mission. A&A, 595:A1, DOI 10.1051/0004-6361/201629272, 1609.04153
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Babusiaux C, van Leeuwen F, Barstow MA, Jordi C, Vallenari A, et al (2018) Gaia Data Release 2. Observational Hertzsprung-Russell diagrams. A&A, 616:A10, DOI 10.1051/0004-6361/201832843, 1804.09378
  * Gilmore et al. (1991) Gilmore G, Edvardsson B, Nissen PE (1991) First Detection of Beryllium in a Very Metal Poor Star: A Test of the Standard Big Bang Model. ApJ, 378:17, DOI 10.1086/170402
  * Gilmore et al. (1992) Gilmore G, Gustafsson B, Edvardsson B, Nissen PE (1992) Is beryllium in metal-poor stars of galactic or cosmological origin? Nature, 357(6377):379–384, DOI 10.1038/357379a0
  * Gustafsson et al. (2008) Gustafsson B, Edvardsson B, Eriksson K, Jørgensen UG, Nordlund Å, Plez B (2008) A grid of MARCS model atmospheres for late-type stars. I. Methods and general properties. A&A, 486(3):951–970, DOI 10.1051/0004-6361:200809724, 0805.0554
  * Haywood et al. (2018) Haywood M, Di Matteo P, Lehnert MD, Snaith O, Khoperskov S, Gómez A (2018) In Disguise or Out of Reach: First Clues about In Situ and Accreted Stars in the Stellar Halo of the Milky Way from Gaia DR2. ApJ, 863(2):113, DOI 10.3847/1538-4357/aad235, 1805.02617
  * Helmi et al. (2018) Helmi A, Babusiaux C, Koppelman HH, Massari D, Veljanoski J, Brown AGA (2018) The merger that led to the formation of the Milky Way’s inner stellar halo and thick disk. Nature, 563(7729):85–88, DOI 10.1038/s41586-018-0625-x, 1806.06038
  * Ito et al. (2013) Ito H, Aoki W, Beers TC, Tominaga N, Honda S, Carollo D (2013) Chemical Analysis of the Ninth Magnitude Carbon-enhanced Metal-poor Star BD+44°493. ApJ, 773(1):33, DOI 10.1088/0004-637X/773/1/33, 1306.3614
  * Korotin and Kučinskas (2022) Korotin S, Kučinskas A (2022) Abundance of beryllium in the Sun and stars: The role of non-local thermodynamic equilibrium effects. A&A, 657:L11, DOI 10.1051/0004-6361/202142789, 2201.00532
  * Kramida (2005) Kramida AE (2005) Critical Compilation of Wavelengths and Energy Levels of Singly Ionized Beryllium (Be II). Physica Scripta, 72(4):309–319, DOI 10.1238/Physica.Regular.072a00309
  * Krieger et al. (2017) Krieger A, Nörtershäuser W, Geppert C, Blaum K, Bissell ML, Frömmgen N, Hammen M, Kreim K, Kowalska M, Krämer J, Neugart R, Neyens G, Sánchez R, Tiedemann D, Yordanov DT, Zakova M (2017) Frequency-comb referenced collinear laser spectroscopy of Be+ for nuclear structure investigations and many-body QED tests. Applied Physics B: Lasers and Optics 123(1):15, DOI 10.1007/s00340-016-6579-5, 1609.07655
  * Kusakabe et al. (2014) Kusakabe M, Kim KS, Cheoun MK, Kajino T, Kino Y, Mathews GJ (2014) Revised Big Bang Nucleosynthesis with Long-lived, Negatively Charged Massive Particles: Updated Recombination Rates, Primordial 9Be Nucleosynthesis, and Impact of New 6Li Limits. ApJS, 214(1):5, DOI 10.1088/0067-0049/214/1/5, 1403.4156
  * Kusakabe et al. (2017) Kusakabe M, Mathews GJ, Kajino T, Cheoun MK (2017) Review on effects of long-lived negatively charged massive particles on Big Bang Nucleosynthesis. International Journal of Modern Physics E 26(8):1741004-64, DOI 10.1142/S021830131741004X, 1706.03143
  * Lingenfelter (2019) Lingenfelter RE (2019) The Origin of Cosmic Rays: How Their Composition Defines Their Sources and Sites and the Processes of Their Mixing, Injection, and Acceleration. ApJS, 245(2):30, DOI 10.3847/1538-4365/ab4b58, 1903.06330
  * Mathar (2011) Mathar RJ (2011) Corrigendum to “Universal factorization of 3n-j (j¿2) symbols…” [J. Phys. A: Math. Gen. 37 (2004) 3259]. arXiv e-prints arXiv:1102.5125, 1102.5125
  * Molaro and Beckman (1984) Molaro P, Beckman J (1984) An upper limit to the abundance of 9Be in the population II star HD 76932 from a high resolution spectrum with the IUE. A&A, 139:394–400
  * Molaro et al. (1984) Molaro P, Beckman JE, Castelli F (1984) An upper limit to the beryllium abundance in the population II star HD 140283 from a high resolution IUE spectrum. In: Rolfe E (ed) Fourth European IUE Conference, ESA Special Publication, vol 218, pp 197–201
  * Molaro et al. (2020) Molaro P, Cescutti G, Fu X (2020) Lithium and beryllium in the Gaia-Enceladus galaxy. MNRAS, 496(3):2902–2909, DOI 10.1093/mnras/staa1653, 2006.00787
  * Nörtershäuser et al. (2015) Nörtershäuser W, Geppert C, Krieger A, Pachucki K, Puchalski M, Blaum K, Bissell ML, Frömmgen N, Hammen M, Kowalska M, Krämer J, Kreim K, Neugart R, Neyens G, Sánchez R, Yordanov DT (2015) Precision Test of Many-Body QED in the Be+ 2 p Fine Structure Doublet Using Short-Lived Isotopes. Phys. Rev. Lett. 115(3):033002, DOI 10.1103/PhysRevLett.115.033002, 1507.03830
  * Orito et al. (1997) Orito M, Kajino T, Boyd RN, Mathews GJ (1997) Geometrical Effects of Baryon Density Inhomogeneities on Primordial Nucleosynthesis. ApJ, 488(2):515–523, DOI 10.1086/304716, astro-ph/9609130
  * Pasquini (2014) Pasquini L (2014) UV opportunities at ESO. Ap&SS, 354(1):121–124, DOI 10.1007/s10509-014-2049-x
  * Pasquini et al. (2004) Pasquini L, Bonifacio P, Randich S, Galli D, Gratton RG (2004) Beryllium in turnoff stars of NGC 6397: Early Galaxy spallation, cosmochronology and cluster formation. A&A, 426:651–657, DOI 10.1051/0004-6361:20041254, astro-ph/0407524
  * Pasquini et al. (2005) Pasquini L, Galli D, Gratton RG, Bonifacio P, Randich S, Valle G (2005) Early star formation in the Galaxy from beryllium and oxygen abundances. A&A, 436(3):L57–L60, DOI 10.1051/0004-6361:200500124, astro-ph/0505396
  * Pasquini et al. (2007) Pasquini L, Bonifacio P, Randich S, Galli D, Gratton RG, Wolff B (2007) Beryllium abundance in turn-off stars of NGC 6752. A&A, 464(2):601–607, DOI 10.1051/0004-6361:20066260, astro-ph/0612710
  * Placco et al. (2014) Placco VM, Beers TC, Roederer IU, Cowan JJ, Frebel A, Filler D, Ivans II, Lawler JE, Schatz H, Sneden C, Sobeck JS, Aoki W, Smith VV (2014) Hubble Space Telescope Near-ultraviolet Spectroscopy of the Bright CEMP-no Star BD+44°493. ApJ, 790(1):34, DOI 10.1088/0004-637X/790/1/34, 1406.0538
  * Plez (2012) Plez B (2012) Turbospectrum: Code for spectral synthesis. 1205.004
  * Pospelov et al. (2008) Pospelov M, Pradler J, Steffen FD (2008) Constraints on supersymmetric models from catalytic primordial nucleosynthesis of beryllium. J. of Cosmol. and Astropart. Phys., 2008(11):020, DOI 10.1088/1475-7516/2008/11/020, 0807.4287
  * Prantzos (2012) Prantzos N (2012) Production and evolution of Li, Be, and B isotopes in the Galaxy. A&A, 542:A67, DOI 10.1051/0004-6361/201219043, 1203.5662
  * Primas et al. (2000a) Primas F, Asplund M, Nissen PE, Hill V (2000a) The beryllium abundance in the very metal-poor halo star G 64-12 from VLT/UVES observations. A&A, 364:L42–L46, astro-ph/0009482
  * Primas et al. (2000b) Primas F, Molaro P, Bonifacio P, Hill V (2000b) First UVES observations of beryllium in very metal-poor stars. A&A, 362:666–672, astro-ph/0008402
  * Puchalski and Pachucki (2009) Puchalski M, Pachucki K (2009) Fine and hyperfine splitting of the 2P state in Li and Be+. Physical Review A 79(3):032510, DOI 10.1103/PhysRevA.79.032510, 0901.2633
  * Puchalski and Pachucki (2014) Puchalski M, Pachucki K (2014) Ground-state hyperfine splitting in the Be+ ion. Physical Review A 89(3):032510, DOI 10.1103/PhysRevA.89.032510, 1402.4573
  * Puchalski and Pachucki (2015) Puchalski M, Pachucki K (2015) Quantum electrodynamics m $\alpha$6 and m $\alpha$7ln$\alpha$ corrections to the fine splitting in Li and Be+. Physical Review A 92(1):012513, DOI 10.1103/PhysRevA.92.012513, 1506.02462
  * Rebolo et al. (1988) Rebolo R, Molaro P, Abia C, Beckman JE (1988) Abundances of 9Be in a sample of highly metal-deficient dwarfs : implications for early galactic nucleosynthesis and primordial lithium. A&A, 193:193–201
  * Rich and Boesgaard (2009) Rich JA, Boesgaard AM (2009) Beryllium, Oxygen, and Iron Abundances in Extremely Metal-Deficient Stars. ApJ, 701(2):1519–1533, DOI 10.1088/0004-637X/701/2/1519, 0906.3296
  * Ryan et al. (1992) Ryan SG, Norris JE, Bessell MS, Deliyannis C (1992) Evolution of Beryllium Abundances in the Galactic Halo. ApJ, 388:184, DOI 10.1086/171141
  * Safronova and Safronova (2013) Safronova UI, Safronova MS (2013) Relativistic many-body calculation of energies, lifetimes, polarizabilities, and hyperpolarizabilities in Li-like Be+. Physical Review A 87(3):032502, DOI 10.1103/PhysRevA.87.032502
  * Shukla et al. (2020) Shukla N, Arora B, Sharma L, Srivastava R (2020) Two-dipole and three-dipole dispersion coefficients for interaction of alkaline-earth-metal atoms with alkaline-earth-metal atoms and alkaline-earth-metal ions. Physical Review A 102(2):022817, DOI 10.1103/PhysRevA.102.022817, 2008.04341
  * Smiljanic (2014) Smiljanic R (2014) Stellar abundances of beryllium and CUBES. Ap&SS, 354(1):55–64, DOI 10.1007/s10509-014-1916-9, 1403.6276
  * Smiljanic et al. (2009) Smiljanic R, Pasquini L, Bonifacio P, Galli D, Gratton RG, Randich S, Wolff B (2009) Beryllium abundances and star formation in the halo and in the thick disk. A&A, 499(1):103–119, DOI 10.1051/0004-6361/200810592, 0902.0483
  * Smiljanic et al. (2021) Smiljanic R, Zych MG, Pasquini L (2021) Inhomogeneity in the early Galactic chemical enrichment exposed by beryllium abundances in extremely metal-poor stars. A&A, 646:A70, DOI 10.1051/0004-6361/202039101, 2012.07438
  * Soubiran et al. (2016) Soubiran C, Le Campion JF, Brouillet N, Chemin L (2016) The PASTEL catalogue: 2016 version. A&A, 591:A118, DOI 10.1051/0004-6361/201628497, 1605.07384
  * Spite et al. (2019) Spite M, Bonifacio P, Spite F, Caffau E, Sbordone L, Gallagher AJ (2019) Be and O in the ultra metal-poor dwarf 2MASS J18082002-5104378: the Be-O correlation. A&A, 624:A44, DOI 10.1051/0004-6361/201834741, 1902.11048
  * Suda et al. (2008) Suda T, Katsuta Y, Yamada S, Suwa T, Ishizuka C, Komiya Y, Sorai K, Aikawa M, Fujimoto MY (2008) Stellar Abundances for the Galactic Archeology (SAGA) Database — Compilation of the Characteristics of Known Extremely Metal-Poor Stars. PASJ, 60:1159, DOI 10.1093/pasj/60.5.1159, 0806.3697
  * Suzuki and Yoshii (2001) Suzuki TK, Yoshii Y (2001) A New Model for the Evolution of Light Elements in an Inhomogeneous Galactic Halo. ApJ, 549(1):303–319, DOI 10.1086/319049, astro-ph/0010108
  * Suzuki et al. (1999) Suzuki TK, Yoshii Y, Kajino T (1999) Evolution of Beryllium and Boron in the Inhomogeneous Early Galaxy. ApJ Letters, 522(2):L125–L128, DOI 10.1086/312233, astro-ph/9907182
  * Tan and Zhao (2011) Tan K, Zhao G (2011) A Possible Signature of Non-uniform Be-$\alpha$ Relationships for the Galaxy. ApJ Letters, 738(2):L33, DOI 10.1088/2041-8205/738/2/L33, 1108.2074
  * Tan et al. (2009) Tan KF, Shi JR, Zhao G (2009) Beryllium abundances in metal-poor stars. MNRAS, 392(1):205–215, DOI 10.1111/j.1365-2966.2008.14027.x, 0810.2600
  * Tang et al. (2009) Tang LY, Zhang JY, Yan ZC, Shi TY, Babb JF, Mitroy J (2009) Calculations of polarizabilities and hyperpolarizabilities for the Be+ ion. Physical Review A 80(4):042511, DOI 10.1103/PhysRevA.80.042511, 0908.4060
  * Tatischeff and Gabici (2018) Tatischeff V, Gabici S (2018) Particle Acceleration by Supernova Shocks and Spallogenic Nucleosynthesis of Light Elements. Annual Review of Nuclear and Particle Science 68(1):377–404, DOI 10.1146/annurev-nucl-101917-021151, 1803.01794
  * Valle et al. (2002) Valle G, Ferrini F, Galli D, Shore SN (2002) Evolution of Li, Be, and B in the Galaxy. ApJ, 566(1):252–260, DOI 10.1086/338036, astro-ph/0110327
  * Woodgate (1970) Woodgate GK (1970) Elementary Atomic Structure. Oxford University Press
  * Yan et al. (1998) Yan ZC, Tambasco M, Drake GWF (1998) Energies and oscillator strengths for lithiumlike ions. Physical Review A 57(3):1652–1661, DOI 10.1103/PhysRevA.57.1652
  * Yerokhin (2008) Yerokhin VA (2008) Hyperfine structure of Li and Be+. Physical Review A 78(1):012513, DOI 10.1103/PhysRevA.78.012513, 0805.0677


# Data-driven learning of nonlocal models:
from high-fidelity simulations to constitutive laws

Huaiqian You,1 Yue Yu,1 Stewart Silling,2 Marta D’Elia3

###### Abstract

We show that machine learning can improve the accuracy of simulations of
stress waves in one-dimensional composite materials. We propose a data-driven
technique to learn nonlocal constitutive laws for stress wave propagation
models. The method is an optimization-based technique in which the nonlocal
kernel function is approximated via Bernstein polynomials. The kernel,
including both its functional form and parameters, is derived so that when
used in a nonlocal solver, it generates solutions that closely match high-
fidelity data. The optimal kernel therefore acts as a homogenized nonlocal
continuum model that accurately reproduces wave motion in a smaller-scale,
more detailed model that can include multiple materials. We apply this
technique to wave propagation within a heterogeneous bar with a periodic
microstructure. Several one-dimensional numerical tests illustrate the
accuracy of our algorithm. The optimal kernel is demonstrated to reproduce
high-fidelity data for a composite material in applications that are
substantially different from the problems used as training data.

## Introduction

Nonlocal models use integral operators acting on a lengthscale $\delta$, known
as horizon. This feature allows nonlocal models to capture long-range forces
at small scales and multiscale behavior, and to reduce regularity requirements
on the solutions, which are allowed to be discontinuous or even singular. In
recent decades, nonlocal equations have been successfully used to model
several engineering and scientific applications, including fracture mechanics
(Silling 2000; Ha and Bobaru 2011; Trask et al. 2019), subsurface transport
(Benson, Wheatcraft, and Meerschaert 2000; Schumer et al. 2003), image
processing (D’Elia, De los Reyes, and Trujillo 2019; Gilboa and Osher 2007),
multiscale and multiphysics systems (Alali and Lipton 2012; Askari et al.
2008; You, Yu, and Kamensky 2020), finance (Scalas, Gorenflo, and Mainardi
2000), and stochastic processes (D’Elia et al. 2017; Meerschaert and Sikorskii
2012).

However, it is often the case that nonlocal kernels defining nonlocal
operators are justified a posteriori and it is not clear how to define such
kernels to faithfully describe a physical system. The problem of learning an
appropriate kernel for a specific application is one of the most challenging
open problems in nonlocal modeling. The literature on techniques for learning
kernel parameters for a given functional form is vast, see, e.g., (Burkovska,
Glusa, and D’Elia 2020; D’Elia and Gunzburger 2016) for control-based
approaches, and (Pang et al. 2020; Pang, Lu, and Karniadakis 2019) for
machine-learning approaches. However, the use of machine learning to learn the
functional form of the kernel is still in its infancy, (Xu and Foster 2020;
Xu, D’Elia, and Foster 2020; You et al. 2020) being the only relevant works
that we are aware of.

In this work we use an approach similar to the one developed in (You et al.
2020) to learn nonlocal kernels whose associated nonlocal wave equation is
well posed by construction and can be used as an accurate surrogate for more
detailed, high-fidelity wave propagation models. In particular, we present an
application to wave propagation at the microscale in a heterogeneous solid. In
this context, the machine-learned nonlocal kernel embeds the material
constitutive behavior so that the material interfaces do not have to be
treated explicitly and, more importantly, the material microstructure can be
unknown. Furthermore, the corresponding nonlocal models allow for accurate
simulations at scales that are much larger than the microstructure.

Our main contributions are:
$\bullet$ The design of an optimization technique that bridges micro and
continuum scales by providing accurate and stable model surrogates for the
simulation of wave propagation in heterogeneous materials.
$\bullet$ The illustration of this method via one-dimensional experiments that
confirm the applicability of our technique and the improved accuracy compared
with state-of-the-art results.
$\bullet$ The demonstration of generalization properties of our algorithm
whose associated model surrogates are effective even on problem settings that
are substantially different from the ones used for training in terms of
loading and time scales.

## Nonlocal kernel learning

We introduce the high-fidelity (HF) model that faithfully represents the
system: for ${\Omega}\in{\mathbb{R}^{d}}$, the scalar function $u(x,t)$
solves, for $(x,t)\in\Omega\times[0,T]$

$\dfrac{\partial^{2}u}{\partial
t^{2}}(x,t)-\mathcal{L}_{\text{HF}}[u](x,t)=f(x,t),$ (1)

provided some boundary conditions on $\partial\Omega$ for $u(x,t)$ and initial
conditions at $t=0$ for $u$ and $\partial u/\partial t$ are satisfied. Here,
$\mathcal{L}_{\text{HF}}$ is the HF operator, which can either be a
differential or integral operator, and $f$ represents a forcing term.

We assume that solutions to this HF problem may be approximated by solutions
to a nonlocal problem of the form

$\dfrac{\partial^{2}u}{\partial t^{2}}(x,t)-\mathcal{L}_{K}[u](x,t)=f(x,t),$
(2)

for $(x,t)\in\Omega\times[0,T]$, augmented with nonlocal boundary conditions
on $\Omega_{\delta}$ (a layer of thickness $\delta$ that surrounds the domain)
and the same initial conditions on the variable $u$ and its derivative as in
(1). The forcing $f$ may coincide with the forcing term in (1) or it could be
an appropriate representation of the same.

We seek $\mathcal{L}_{K}$ as a nonlocal operator of the form

$\mathcal{L}_{K}[u](x,t)=\int_{\overline{\Omega}}K(|x-y|)\left(u(y,t)-u(x,t)\right)\,dy,$
(3)

where $K$ is a radial, sign-changing, kernel function, compactly supported on
the ball of radius $\delta$ centered at $x$, i.e., $B_{\delta}(x)$ and
$\overline{\Omega}=\Omega\cup\Omega_{\delta}$.

### The algorithm

To learn the kernel $K$, we assume that we are given $N$ pairs of forcing
terms and corresponding solutions to (1), normalized with respect to the
$L^{2}$ norm of each solution over $\Omega\times[0,T_{\rm tr}]$. These are
denoted by

$\mathcal{D}_{\rm tr}=\left\\{(u_{k}(x,t),f_{k}(x,t))\right\\}_{k=1}^{N},$ (4)

for $x\in\Omega$ and $t\in(0,T_{\rm tr}]$. Similarly to (You et al. 2020), we
represent $K$ as a linear combination of Bernstein basis polynomials:

$\displaystyle K\left(\frac{|y|}{\delta}\right)$
$\displaystyle=\sum_{m=0}^{M}\frac{C_{m}}{\delta^{d+2}}B_{m,M}\bigg{(}\bigg{|}\frac{y}{\delta}\bigg{|}\bigg{)},$
(5)

where the Bernstein basis functions are defined as

$B_{m,M}(x)=\begin{pmatrix}M\\\ m\\\ \end{pmatrix}x^{m}(1-x)^{M-m}\;\;$ for
$0\leq x\leq 1$
and where $C_{m}\in\mathbb{R}$. Note that, by construction, this kernel
guarantees that (2) is well-posed (Du, Tao, and Tian 2018).

We machine-learn the nonlocal model by finding optimal parameters
$\\{C_{m}\\}$ such that solutions $\hat{u}_{k}$ to (2), for $f=f_{k}$ and the
kernel function $K$ associated to $\\{C_{m}\\}$, are as close as possible to
the training variable $u_{k}$.

In this work we numerically approximate $\hat{u}_{k}$ by $\bar{u}_{k}$ using a
central-differencing scheme in time with time step d$t$, i.e.

$\displaystyle\bar{u}_{k}^{n+1}(x_{i})$
$\displaystyle=2\bar{u}^{n}_{k}(x_{i})-\bar{u}^{n-1}_{k}(x_{i})$ (6)
$\displaystyle+{\rm
d}t^{2}\left(\mathcal{L}_{K,h}[\bar{u}^{n}_{k}](x_{i})+f_{k}(x_{i},t^{n})\right),$

where $\bar{u}_{k}^{n+1}(x_{i})$ represents the $k$-th approximate solution at
time step $t^{n+1}$ and at discretization point $x_{i}$, and
$\mathcal{L}_{K,h}$ is an approximation of $\mathcal{L}_{K}$ by Riemann sum
with uniform grid spacing $h$. The optimal parameters are obtained by solving
the following optimization problem.

$\displaystyle\min_{C_{m}}$ $\displaystyle\frac{T_{\rm tr}}{{\rm
d}t^{3}\,N}\sum_{k=1}^{N}\sum_{n=1}^{T_{\rm tr}/{\rm
d}t}\big{\|}\bar{u}_{k}^{n+1}-u_{k}(t^{n+1})\big{\|}_{\ell^{2}}^{2}+\mathcal{R}(\\{C_{m}\\}),$
(7) s.t. $\;\;\bar{u}_{k}$ satisfies (6) and (8) $K$ satisfies physics-based
constraints. (9)

Here, the $\ell^{2}$ norm is taken over the space-discretization points
$x_{i}$, $\mathcal{R}(\cdot)$ is a regularization term on the coefficients
that improves the conditioning of the optimization problem, and (9) depends on
the physics of the problem (as an example, it may correspond to enforcing that
the surrogate model reproduces exactly a certain class of solutions).

## Dispersion in heterogeneous materials

We apply the learning algorithm described above to the propagation of waves in
a one-dimensional heterogeneous bar, like the one reported in Figure 1, with
an ordered microstructure, i.e. two materials with the same length alternate
periodically. Our goal is to learn a nonlocal model able to reproduce wave
propagation on distances that are much larger than the size of the
microstructure without resolving the microscales. The high-fidelity model we
rely on is the classical wave equation; the corresponding high-fidelity data
used for training and validation are obtained with the solver described below.

Figure 1: One-dimensional bar with ordered microstructure of period $2L$.
Material 1 and 2 have the same density and Young modulus $E_{1}$ and $E_{2}$.
The horizon $\delta$, the wave length $\lambda$, and the discretization size,
$h$, are reported for comparison.

### High-fidelity data

For both training and validation purposes we generate data using high-fidelity
simulations for the propagation of stress waves within the microstructure of
the heterogeneous, linear elastic bar. This method, which will be referred to
as Direct Numerical Solution (DNS), constructs an arbitrarily complex wave
diagram (also called an $x$-$t$ diagram), that treats the mutual interaction
and superposition of many wavefronts moving in either direction. The bar is
discretized into nodes such that it takes a constant amount of time $\Delta t$
for a wave to travel between nodes $\gamma$ and $\gamma+1$, regardless of the
elastic wave speed in the material between these two nodes. Therefore, in a
heterogeneous medium, the spacing between nodes is not constant. Each node
$\gamma$, at each time step $n$, has velocity $v_{\gamma}^{n}$ (note that, in
this case, the subscript refers to position, as opposed to the previous
section where it corresponds to a specific sample $k$). To compute the
velocities in the next time step, it is assumed that two waves moving in
opposite directions converge on the node $\gamma$ at time step $n$ (see Figure
2). The waves shown in the figure can have unequal slopes on the $x$-$t$
diagram because the materials on either side of node $\gamma$ can have
different waves speeds $c$. The jump conditions for the waves are applied that
relate the stress change $[\sigma]$ across a wave to the velocity change
$[v]$. These jump conditions have the following form:

$[\sigma]=\pm\rho c[v],$

where $\rho$ is the mass density, and where the $+$ and $-$ signs apply to
right-running waves and left-running waves respectively. From these
conditions, the velocity $v_{\gamma}^{n+1}$ can be computed explicitly from
the values at the adjacent nodes in time step $n-1$. Externally applied forces
can also be included in a straightforward way. After $v_{\gamma}^{n+1}$ is
computed, the updated displacement is approximated by

$u_{\gamma}^{n+1}=u_{\gamma}^{n}+\Delta tv_{\gamma}^{n+1}.$

Details of the method can be found in (Silling 2020).

This DNS solver has the important advantage of not using an approximate
representation of derivatives in space or time for the computation of the
velocity, which is, therefore, free from truncation error and other sources of
discretization error that are usually encountered with PDE solvers. This
allows us to model the propagation of waves through many thousands of
microstructural interfaces without the need to worry about what features of
the velocity are real and what are numerical artifacts.

Figure 2: Interaction of two waves in the DNS method. Each node may or may not
be located at a material interface.

We consider four types of data and use the first two for training and the last
two for validation of our algorithm. In all our experiments we set
$L\\!=\\!0.2$, $E_{1}\\!=\\!1$, $E_{2}\\!=\\!0.25$, $\rho\\!=\\!1$, and the
symmetric domain $\Omega\\!=\\!(-b,b)$. Discretization parameters for the DNS
solver are set to $\Delta t\\!=\\!0.01$ and $\max\\{\Delta x\\}\\!=\\!0.01$.

1) Oscillating source. We set $b=50$, $v(x,0)=u(x,0)=0$,

$f(x,t)\\!=e^{-\left(\frac{2x}{5kL}\right)^{2}}\\!e^{-\left(\frac{t-t_{0}}{t_{p}}\right)^{2}}\\!\cos^{2}\left(\frac{2\pi
x}{kL}\right),\,k=1,2,\ldots,20,$ $t_{0}=t_{p}=0.8$.

2) Plane wave. For $b=50$, $f(x,t)=0$ and $u(x,0)=0$, we prescribe
$v(-b,t)=\sin(\omega t)$ for $\omega=0.35,0.7,\cdots,3.85$.

3) Wave packet. For $b=133.3$, $f(x,t)=0$ and $u(x,0)=0$, we prescribe
$v(-b,t)=\sin(\omega t)\exp(-(t/5-3)^{2})$ for $\omega=2,3.9,5$.

4) Impact. For $b=$266.6, $f(x,t)=0$ and $u(x,0)=0$, we prescribe $v(x,0)=1$
for all $x\in[-b,-b+1.6]$ and $v=0$ outside of this interval. This initial
condition represents an impactor hitting the bar at time zero, generating a
velocity pulse of width roughly 3.2 that propagates into the interior of the
bar. The pulse attenuates and changes shape as it encounters the many
microstructural interfaces.

### Training procedure

For the optimization problem (7) we choose a Tikhonov regularization of the
form $\mathcal{R}(\\{C_{m}\\})=\frac{\epsilon}{M+1}\sum_{m=0}^{M}C_{m}^{2}$,
where the regularization weight $\epsilon$ is chosen empirically to guarantee
accurate predictions, as we explain later on. The physics-based constraints in
(9) are defined as follows and also discretized by Riemann sum; they are used
to explicitly prescribe values of $C_{M-1}$ and $C_{M}$:

$\displaystyle\sum_{m=0}^{M}C_{m}$
$\displaystyle\int_{0}^{\delta}\frac{y^{2}}{\delta^{3}}B_{m,M}\bigg{(}\frac{|y|}{\delta}\bigg{)}dy=\rho
c^{2}_{0},$ (10) $\displaystyle\sum_{m=0}^{M}C_{m}$
$\displaystyle\int_{0}^{\delta}\frac{y^{4}}{\delta^{3}}B_{m,M}\bigg{(}\frac{|y|}{\delta}\bigg{)}dy=-4\rho
c_{0}^{3}R,$

where $\rho$ is the density and $c_{0}$ is the effective wave speed for
infinitely long wavelengths. For $\rho=1$, it is given by
$c_{0}=(2/(1/E_{1}+1/E_{2})^{\frac{1}{2}}$. $R$ is the second derivative of
the wave group velocity with respect to the frequency $\omega$ evaluated at
$\omega=0$. Both parameters are obtained by simulating a very low frequency
plane wave propagating through the microstructure over a long distance using
DNS (Silling 2020). These parameters primarily affect simulations at large
times, $t>10$. However, due to practical limitations on computer resources,
our training simulations are restricted to $t\leq 2$. Therefore, we
incorporate these parameters as constraints obtained from DNS as indicated in
(10), rather than attempting to learn these through our algorithm. The first
constraint in (10) is also used for similar purposes in (Xu, D’Elia, and
Foster 2020) and prescribed in a weak sense by penalization.

Training is performed with DNS data of type 1) and 2). Parameters for the
nonlocal solver and the optimization algorithm are set to $h=0.05$, ${\rm
d}t=0.02$, $T_{\rm tr}=2$, $\delta=$1.2, $M=24$ and $\epsilon=0.01$. The
optimization problem (7) is solved with L-BFGS. Note that we empirically
choose $\delta$ and $\epsilon$ in such a way that the group velocity, defined
below, corresponding to the optimal kernel is as close as possible to the one
computed with DNS.

The optimal kernel, $K_{\rm opt}$, is reported in Figure 3; as expected from
the literature (Xu and Foster 2020; Xu, D’Elia, and Foster 2020; You et al.
2020; Weckner and Silling 2011), we observe a sign-changing behavior. We also
compute the corresponding dispersion $\omega(k)$ and group velocity
$v_{g}(\omega)=d\omega/dk$. For a given kernel $K$ and different frequencies
$k_{i}=0,\frac{2\pi}{200h},\cdots,\frac{2\pi}{h}$, the corresponding angular
frequency $\omega(k_{i})$ and group velocity $v_{g}(\omega(k_{i}))$ are
approximated by

$\displaystyle\omega(k_{i})^{2}$
$\displaystyle\approx\frac{1}{\rho}\sum_{q}K(|y_{q}|)(1-\cos(k_{i}y_{q}))h,$
$\displaystyle v_{g}(\omega(k_{i}))$
$\displaystyle\approx\dfrac{\omega(k_{i+1})-\omega(k_{i-1})}{k_{i+1}-k_{i-1}},$

where $y_{q}$ belong to a uniform grid of size $h$ in $(-\delta,\delta)$.

The dispersion curve is reported in Figure 4, its positivity indicates that
$K_{\rm opt}$ corresponds to a physically stable material model. The group
velocity is reported in the upper plot of Figure 5 in comparison with the
curve computed with DNS by observing the speed of a wave packet of a given
frequency as it moves through the microstructure. We also display the group
velocity associated with an alternative kernel obtained for the same material
by a completely different method (Silling 2020). This alternative kernel is a
constant, specifically, we have that $C_{m}\\!=\\!K_{\rm const}\\!=\\!0.7714$,
for $M\\!=\\!3$ and $\delta\\!=\\!0.15$.

Figure 3: Optimal kernel $K_{\rm opt}$ as a function of distance. Figure 4:
Dispersion curve associated with $K_{\rm opt}$.

Figure 5: Comparison of the group velocity. Upper: group velocity
corresponding to $K_{\rm opt}$, $K_{\rm const}$, and DNS data. Bottom: group
velocity corresponding to $K_{\rm opt}$ for different pairs
$(\delta,\epsilon)$.

It is well known that layered, periodic elastic media have a band structure
for wave propagation, see (Bedford and Drumheller 1994, pages 121-122). In the
present study, because it is not possible to reproduce the higher-frequency
pass bands with the coarse discretization that is used in fitting our nonlocal
kernel, we address only the first, low-frequency pass band, i.e.
$\omega\in(0,\omega_{\rm bs})$, where “bs” stands for band stop. Hence, the
optimal kernel is suitable only for wavelengths that are bigger than the
microstructure; this is enough to reproduce the physically most important
features of wave propagation in layered media for typical applications.

The profile of the group velocity shows the improved accuracy of our optimal
kernel that not only matches the behavior for low values of $\omega$, but also
catches the behavior at $\omega=\omega_{\rm bs}\approx 4$. This fact has
important consequences on the ability to reproduce wave propagation for values
of $\omega$ bigger than $\omega_{\rm bs}$. In order to justify the statement
above on the optimality of the parameters $\delta$ and $\epsilon$, we report
the group velocity profile in correspondence of different pairs; it is clear
from the profiles in Figure 5 that $(\delta,\epsilon)=(1.2,0.01)$ provides the
best match both in terms of curvature at $\omega=0$ and identification of the
band stop.

Figure 6: Velocity computed with $K_{\rm opt}$. Plots from top to bottom
correspond to: 1. $\omega_{1}=2$ at $t=100$; 2. $\omega_{2}=3.9$ at $t=320$;
3. $\omega_{3}=5$ at $t=100$.

Figure 7: Velocity computed with $K_{\rm const}$. Plots from top to bottom
correspond to: 1. $\omega_{1}=2$ at $t=100$; 2. $\omega_{2}=3.9$ at $t=320$;
3. $\omega_{3}=5$ at $t=100$. The optimal kernel $K_{\rm opt}$ obtained by
machine learning (Figure 6) clearly performs better than $K_{\rm const}$ for
the second and third cases ($\omega_{2}$ and $\omega_{3}$).

### Numerical validation

We test the performance of the optimal kernel $K_{\rm opt}$ on data sets of
type 3) and 4), i.e. the problem setting considered for validation has
different model parameters, including the domain, than the one used for
training and, hence, these tests serve as an indicator of the generalization
properties of our algorithm.

Wave packet. For data type 3) we numerically compute solutions to (2) using
$K_{\rm opt}$ and DNS data as nonlocal boundary conditions. We consider
solutions corresponding to three values of $\omega$:
$\omega_{1}\\!=\\!2<\\!\omega_{\rm bs}$,
$\omega_{2}\\!=\\!3.9\approx\omega_{\rm bs}$ and
$\omega_{3}\\!=\\!5>\\!\omega_{\rm bs}$. Note that the latter value is beyond
the band stop and, as such, corresponds to a zero group velocity, i.e. the
wave does not travel in time. In Figure 6 we report the velocity corresponding
to the computed displacement $\bar{u}$ at time $t=100$, $t=320$, and $t=100$
for $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$ respectively; as a reference,
we also report the exact DNS velocity. Our results indicate that our kernel
can accurately reproduce solutions of type 3) at times larger than $T_{\rm
tr}$ and for all values of $\omega$, even larger than $\omega_{\rm bs}$. This
is possible because the group velocity corresponding to $K_{\rm opt}$
reproduces the true group velocity very accurately, see Figure 5. In
particular, detecting the presence of a band stop allows us to accurately
predict the wave propagation for values of $\omega\\!>\\!\omega_{\rm bs}$. Due
to the poor accuracy of the group velocity associated to $K_{\rm const}$,
corresponding solutions are not as accurate for $\omega$ in the proximity of
$\omega_{\rm bs}$ and beyond. To illustrate this phenomenon, we report in
Figure 7 the behavior of the velocity corresponding to $K_{\rm const}$ at time
$t=100$, $t=320$ and $100$, respectively for $\omega_{1}$, $\omega_{2}$ and
$\omega_{3}$. Comparison with DNS data shows that, for $\omega_{2}$ the wave
associated with $K_{\rm const}$ is traveling faster than the exact one and,
for $\omega_{3}$, it keeps traveling while the exact wave is not propagating.

Impact. We use the optimal kernel to compute solutions corresponding to data
type 4). In Figure 8 we report the velocity profile at different time steps in
correspondence of $K_{\rm opt}$ and DNS data, displayed for comparison. Figure
9 displays the same results in correspondence of $K_{\rm const}$. These
results indicate that our optimal kernel can accurately predict the short- and
long-time wave propagation, as opposed to the constant kernel that
successfully predicts the long-time behavior only. We also point out that for
values of $(\delta,\epsilon)$ for which the group velocity is not accurate,
the predicted velocity and displacement exhibit non-physical oscillations,
which disappear in correspondence of pairs that guarantee an accurate group
velocity profile.

Figure 8: Velocity profile for the Impact problem at $T=20$ and $T=600$ with
$K_{\rm opt}$.

Figure 9: Velocity profile for the Impact problem at $T=20$ and $T=600$ with
$K_{\rm const}$. The optimal kernel $K_{\rm opt}$ obtained by machine learning
(Figure 8) provides better agreement with the DNS data than $K_{\rm const}$ at
the earlier time ($T=20$) by reducing the size of the oscillations that trail
the main pulse. For large $T$, the solution is dominated by the low frequency
components of the pulse, for which the two kernels behave similarly.

## Conclusion

We introduced a new data-driven, optimization-based algorithm for the
identification of nonlocal kernels in the context of wave propagation through
material featuring heterogeneities at the microscale. The corresponding
nonlocal model is well-posed by construction and allows for accurate
simulations at a larger scale than the microstructure. We stress the fact that
our algorithm does not require a priori knowledge of the microstructure (often
unknown and/or hard to model), but only requires high-fidelity measurements of
the displacements or the velocity. We also point out that our algorithm has
excellent generalization properties as the optimal kernel performs well at
much larger times than the time instants used for training and on problem
settings that are substantially different from the training data set.

One of the most important findings in this work is the key role of the group
velocity in the accuracy of the predictions; in fact, our criterion for the
choice of the horizon $\delta$ and the regularization weight $\epsilon$ is the
accurate prediction of the group velocity profile. Given the critical role of
such quantity, our future work includes the identification of the optimal
horizon by, possibly, embedding constraints on the group velocity to the
training procedure. Another natural follow-up work is the illustration of the
efficiency of our algorithm on two- and three-dimensional test cases.

## Acknowledgments

MD and SS are supported by the Sandia National Laboratories (SNL) Laboratory-
directed Research and Development program and by the U.S. Department of
Energy, Office of Advanced Scientific Computing Research under the
Collaboratory on Mathematics and Physics-Informed Learning Machines for
Multiscale and Multiphysics Problems (PhILMs) project. SNL is a multimission
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. Department of Energy’s National Nuclear
Security Administration under contract DE-NA-0003525. This paper,
SAND2020-13633, describes objective technical results and analysis. Any
subjective views or opinions that might be expressed in this paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government. HY and YY are supported by the National Science Foundation
under award DMS 1753031.

## References

  * Alali and Lipton (2012) Alali, B.; and Lipton, R. 2012.  Multiscale dynamics of heterogeneous media in the peridynamic formulation.  _Journal of Elasticity_ 106(1): 71–103.
  * Askari et al. (2008) Askari, E.; Bobaru, F.; Lehoucq, R.; Parks, M. L.; Silling, S. A.; and Weckner, O. 2008.  Peridynamics for multiscale materials modeling.  In _Journal of Physics: Conference Series_ , volume 125, 012078. IOP Publishing.
  * Bedford and Drumheller (1994) Bedford, A.; and Drumheller, D. S. 1994.  _Introduction to Elastic Wave Propagation_.  Wiley.
  * Benson, Wheatcraft, and Meerschaert (2000) Benson, D.; Wheatcraft, S.; and Meerschaert, M. 2000.  Application of a fractional advection-dispersion equation.  _Water Resources Research_ 36(6): 1403–1412.
  * Burkovska, Glusa, and D’Elia (2020) Burkovska, O.; Glusa, C.; and D’Elia, M. 2020.  An optimization-based approach to parameter learning for fractional type nonlocal models.  ArXiv:2010.03666.
  * D’Elia, De los Reyes, and Trujillo (2019) D’Elia, M.; De los Reyes, J.; and Trujillo, A. M. 2019.  Bilevel parameter optimization for nonlocal image denoising models.  ArXiv1912.02347.
  * D’Elia et al. (2017) D’Elia, M.; Du, Q.; Gunzburger, M.; and Lehoucq, R. 2017.  Nonlocal convection-diffusion problems on bounded domains and finite-range jump processes.  _Computational Methods in Applied Mathematics_ 29: 71–103.
  * Du, Tao, and Tian (2018) Du, Q.; Tao, Y.; and Tian, X. 2018.  A peridynamic model of fracture mechanics with bond-breaking.  _Journal of Elasticity_ 132(2): 197–218.
  * D’Elia and Gunzburger (2016) D’Elia, M.; and Gunzburger, M. 2016.  Identification of the diffusion parameter in nonlocal steady diffusion problems.  _Applied Mathematics & Optimization_ 73(2): 227–249.
  * Gilboa and Osher (2007) Gilboa, G.; and Osher, S. 2007.  Nonlocal linear image regularization and supervised segmentation.  _Multiscale Model. Simul._ 6: 595–630.
  * Ha and Bobaru (2011) Ha, Y. D.; and Bobaru, F. 2011.  Characteristics of dynamic brittle fracture captured with peridynamics.  _Engineering Fracture Mechanics_ 78(6): 1156–1168.
  * Meerschaert and Sikorskii (2012) Meerschaert, M.; and Sikorskii, A. 2012.  _Stochastic models for fractional calculus_.  Studies in mathematics, Gruyter.
  * Pang et al. (2020) Pang, G.; D’Elia, M.; Parks, M.; and Karniadakis, G. E. 2020.  nPINNs: nonlocal Physics-Informed Neural Networks for a parametrized nonlocal universal Laplacian operator. Algorithms and Applications.  To appear in Journal of Computational Physics.
  * Pang, Lu, and Karniadakis (2019) Pang, G.; Lu, L.; and Karniadakis, G. E. 2019.  fPINNs: Fractional Physics-Informed Neural Networks.  _SIAM Journal on Scientific Computing_ 41: A2603–A2626.
  * Scalas, Gorenflo, and Mainardi (2000) Scalas, E.; Gorenflo, R.; and Mainardi, F. 2000.  Fractional calculus and continuous time finance.  _Physica A_ 284: 376–384.
  * Schumer et al. (2003) Schumer, R.; Benson, D.; Meerschaert, M.; and Baeumer, B. 2003.  Multiscaling fractional advection-dispersion equations and their solutions.  _Water Resources Research_ 39(1): 1022–1032.
  * Silling (2020) Silling, S. 2020.  Propagation of a Stress Pulse in a Heterogeneous Elastic Bar.  Sandia Report SAND2020-8197, Sandia National Laboratories.
  * Silling (2000) Silling, S. A. 2000.  Reformulation of elasticity theory for discontinuities and long-range forces.  _Journal of the Mechanics and Physics of Solids_ 48(1): 175–209.
  * Trask et al. (2019) Trask, N.; You, H.; Yu, Y.; and Parks, M. L. 2019.  An asymptotically compatible meshfree quadrature rule for nonlocal problems with applications to peridynamics.  _Computer Methods in Applied Mechanics and Engineering_ 343: 151–165.
  * Weckner and Silling (2011) Weckner, O.; and Silling, S. A. 2011.  Determination of nonlocal constitutive equations from phonon dispersion relations.  _International Journal for Multiscale Computational Engineering_ 9(6).
  * Xu, D’Elia, and Foster (2020) Xu, X.; D’Elia, M.; and Foster, J. 2020.  Bond-Based Peridynamic Kernel Learning with Energy Constraint.  In preparation.
  * Xu and Foster (2020) Xu, X.; and Foster, J. 2020.  Deriving peridynamic influence functions for one-dimensional elastic materials with periodic microstructure.  ArXiv:2003.05520.
  * You, Yu, and Kamensky (2020) You, H.; Yu, Y.; and Kamensky, D. 2020.  An asymptotically compatible formulation for local-to-nonlocal coupling problems without overlapping regions.  _Computer Methods in Applied Mechanics and Engineering_ 366: 113038\.
  * You et al. (2020) You, H.; Yu, Y.; Trask, N.; Gulian, M.; and D’Elia, M. 2020.  Data-driven learning of robust nonlocal physics from high-fidelity synthetic data.  ArXiv:2005.10076.