{'fraction_non_alphanumeric': 0.04233072179689595, 'fraction_numerical': 0.029700020852572315, 'mean_word_length': 4.277472095582455, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Detecting weak beryllium lines with CUBES

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.

Beryllium has atomic number Z = 4 and a single stable isotope, 9 Be. 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 metallicity 1 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(Duncan et al. , 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 , , 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] ∼ −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Å Evans et al. 2016). The current CUBES design envisions a high-efficiency spectrograph able to operate in low-(R ∼ 6000) and mediumresolution (R ∼ 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.

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 Suzuki and Yoshii 2001). To test this, Pasquini et al. (2004Pasquini et al. ( , 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 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] −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 9 Be. Standard primordial (or Big Bang) nucleosynthesis predicts the formation of tiny amounts of Be, at the level of log(Be/H) ∼ −17.6 (Coc et al. 2014). Several alternative scenarios could result in an enhanced primordial Be abundance of up to log(Be/H) ∼ −15 (Boyd and Kajino 1989;Orito et al. 1997;Pospelov et al. 2008;Kusakabe et al. 2014Kusakabe et al. , 2017. Nevertheless, the Be abundances currently detected in extremely metal-poor stars are at the level of log(Be/H) ∼ −13.5 with a few upper limits down to log(Be/H) ∼ −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

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.

Hyperfine structure of the Be lines
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 Table 2 Constants needed to compute the energy shifts of the hyperfine sublevels (Yerokhin 2008 affect the determination of stellar chemical abundances. The adopted values of log(g f ) 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 −α . The cross section (σ ) and the α 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 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 9 Be is I = 3/2. The selection rules allow transitions between levels with ∆ F = ±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(g f ) of the hyperfine transitions are given in Table 3.

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 mainsequence stars with (T 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 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 (ξ ) of 1.50 km s −1 and using α-enhanced models with [α/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-toend 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 ∼ 23 000 with a sampling of ∼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 ∼ 200 and 100, for the bright and faint cases, respectively. For the subgiant, the values are SNR ∼ 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 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 IRAF 2 . 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.

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: 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 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.

The subgiant star case
As an example, we start looking at the four simulated one hour observations of a bright (V = 12.  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(Be/H) = −13.2. The fourth observed spectrum is fit with a slightly lower value, log(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 ±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 ±0.10 dex. In the fitting of all four simulated observations, the Kolmogorov-Smirnov test indicates that the line was detected. 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 ±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.
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 ±0.05, the R 2 values indicate a clear Fig. 4 Same as Fig. 2, fitting the 3130Å line for an input abundance of log(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. 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.
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 ±0.05 dex of the input, the R 2 plots indicate that abundances within ±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 ±0.15. Only for abundances below  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. 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 ±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 ±0.15 dex.
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 ±0.15 dex. The derived detection limits for each level of SNR are given in Table 4.

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 (  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 ±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 ±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 ±0.15 dex of the input value. The R 2 plots indicates that synthetic spectra with abundances within ±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.

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] ≤ −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 ∼ 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) ≥ −13.4 dex and of log(Be/H) ≥ −13.7 dex, for the main-sequence and subgiant cases, respectively. The typical uncertainty of the fitting only can be of the order of ±0.10-0.15 dex. The SNR ∼ 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 ∼ 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).
A search in the PASTEL catalog (Soubiran et al. 2016) and the SAGA database (Suda et al. 2008) returned 13 known stars with V ≤ 14 mag, log g ≥ 3.4, and [Fe/H] ≤ −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 ∼ 21 known objects, if we can push the magnitude limit down to V ∼ 14.3 mag. Populating the region with [Fe/H] ≤ −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) ∼ −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) ≥ −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.   (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) ≥ −13.7 for stars with −2.5 ≥ [Fe/H] ≥ −3.5 where CUBES can potentially be used to detect Be abundances in extremely metal-poor stars and improve the current statistics.