How magnetic activity alters what we learn from stellar spectra

The study by Galarza et al. (2019) on the solar twin HD 59967 revealed that the EWs of iron lines can vary along the stellar activity cycle of quantities that depend on the mean line-centre optical depth $\tau_{\lambda}$. Namely, lines with log $\tau_{\lambda}\geq$-1 do not show any significant variation along the cycle, but those that form in the external layers of the photosphere have EW values that change with chromospheric activity along the stellar cycle. This effect is also shown in Fig. 1-top where we plot the EW measurements of the Ba II line 5853$\rm\AA$ measured from the HD 59967 spectra at different phases of its activity cycle: at higher activity levels, the EW is clearly higher.

We perform a linear fit of the EWs-log R${}^{\prime}_{\rm HK}$ relation for this and other atomic lines in the HD 59967 spectra. The resulting slopes are plotted in Fig. 1-middle as a function of $\langle$EW$\rangle$, the median EW of the line over all exposures. The $\langle$EW$\rangle$ value can be used as a proxy of the line optical depth $\tau_{\lambda}$, as stronger lines form higher above the optical surface of the stellar atmosphere (Gray, 1992). The red squares in Fig. 1-middle panel and the relative error bars represent the binned averaged slope values and standard deviations calculated at different $\langle$EW$\rangle$ intervals. From the plot we observe that lines with $\langle$EW$\rangle\gtrsim$50 m$\rm\AA$ have slopes that are typically positive, indicating that their EWs increase with the activity index. All the other lines have slopes that are typically consistent with zero, meaning that their EWs have not varied significantly during the stellar cycle.

These trends are likely caused by stellar magnetic fields and cool spots affecting the formation of absorption lines observed in stellar spectra. Magnetic fields can affect spectral lines both directly, through the Zeeman effect, and indirectly, due to magnetically induced changes the temperature and density of the atmospheric plasma in the line formation region. In the spectrum of a star permeated by magnetic fields of some 10mT or more, the Zeeman pattern of many absorption lines is considerably wider than the thermal Doppler profile resulting in artificially larger EWs (Babcock, 1949). The magnitude of magnetic intensification of the line depends on the strength of the magnetic field and on the $\tau_{\lambda}$ of the single line: absorption features that form near the top of the stellar photosphere, where magnetic fields are stronger, undergo a stronger magnetic intensification than lines that form in lower layers. Following the study of Babcock (1949) on magnetic intensification of stellar absorption lines, many theorists have suggested the presence of magnetic fields to explain the observed phenomenologies of pre-main-sequence stars, such as chemical anomalies or the lithium spread in young clusters (e.g., Uchida & Shibata 1984; Leone & Catanzaro 2004; Leone 2007; Oksala et al. 2018).

On the other hand, the indirect effects of magnetic fields arise from the fact that stellar photospheres become more transparent near magnetic concentrations due to the lower density of the plasma. This allows one to probe into deeper and hotter layers of the stellar atmosphere. The hotter temperatures that the radiation “feels” in these regions weakens the absorption lines with higher potential energy (Fabbian et al., 2012). Since we do not observe weakening of lines as a function of the activity index, we consider that the magnetic intensification predicted by Babcock (1949) is a more likely explanation for the EW modulation than other indirect effects caused by magnetic fields. However, it is also possible that indirect effects due to the presence of strong magnetic fields in the vicinity of stellar spots and plage regions, have also affected the absorption lines creating a certain degree of scatter in the EW measurements, probably further modulated by the rotation of the star.

Finally, the same EW modulation can also be explained by the variation of the stellar surface covered by cold spots along the activity cycle. In fact, cold stellar spots can make the stellar photosphere appear cooler, increasing the EW of lines with low energy potential (Gray, 1992). Unfortunately, lines with low energy potentials tend to form at smaller $\tau_{\lambda}$, which means it is impossible to clearly determine if the main cause of the EW variation is the Zeeman broadening, cool stellar spots or a combination of the two. It is also possible that the relative importance of the direct and indirect effects of magnetic fields and stellar spots changes as the star ages, due to the drastic variation of spot filling factors, number of faculae, and strength of magnetic fields that stars undergo across the pre-main-sequence and early stages of the main sequence phases.

Our analysis confirms the conclusions given by Flores et al. (2016) and Galarza et al. (2019) for the solar twins HD59967 and HD45184, respectively. However, the goal of this paper is to pass from the analysis of single objects to the study of a larger number of stars of different ages and typical activity levels. To do so, we first perform a linear fit of the relation between the EW-log R${}^{\prime}_{\rm HK}$ slopes and $\langle$EW$\rangle$ for all the absorption features with $\langle$EW$\rangle\geq$50 m$\rm\AA$. The resulting linear function is represented in Fig. 1-middle as a red solid line. The slope of this linear function is another important parameter that hereafter we call $\Sigma_{\rm EW}$. Similarly, we calculate $\Sigma_{\rm EW}$ for all the other stars in our sample. In Fig. 1-bottom we plot $\Sigma_{\rm EW}$ values along with their uncertainties as a function of $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$ for all stars in our sample, where $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$ is the median of the activity indexes measured at all epochs. The red squares represent the binned-averaged $\Sigma_{\rm EW}$ values and their standard deviations. The plot shows that the variation on EW during the activity cycle becomes significant for $\langle$log R${}^{\prime}_{\rm HK}$$\rangle\gtrsim-$5.0 dex. Accordingly to the age-log R${}^{\prime}_{\rm HK}$ relation calibrated by Lorenzo-Oliveira et al. (2018) on solar twin stars, stars with log R${}^{\prime}_{\rm HK}>-$5.0 dex are typically younger than $\sim$5 Gyr. The sensitivity of the EWs to the variation of chromospheric activity during the stellar cycle increases with the median stellar activity and, consequently, it is more significant for younger stars.

Another interesting outcome of this analysis is that the indirect effects of magnetic fields on stellar spectra are negligible in relation to the Zeeman broadening of atomic lines or the effect of cool stellar spots. This is a general behaviour of all the most active stars in our sample, regardless of other possible key factors, such as the morphology of magnetic fields, that can vary from star to star and that could determine the magnitude by which magnetic fields influence stellar spectra (Moore et al., 2015; Shchukina & Trujillo Bueno, 2015).

Finally, we test whether Voigt profiles measure changes in EWs along the stellar cycle that are different than those obtained with Gaussian profiles. In Fig. 2 we show the EWs of the Ba II line at 5853 $\rm\AA$ measured with Voigt profiles in HD59967 spectra as a function of log R${}^{\prime}_{\rm HK}$. These EWs are typically $\sim$8 m$\rm\AA$ larger than those obtained with Gaussians (i.e., see Fig 1-top), because a Voigt profile can better capture the damping wings of stronger lines. However, the EW variation as a function of log R${}^{\prime}_{\rm HK}$ traced by Voigt profile is marginally consistent with that observed with Gaussians: while the first gave a EW-log R${}^{\prime}_{\rm HK}$ slope equal to 33$\pm$5 m$\rm\AA$, the use of Gaussians produced a slope of 43$\pm$4 m$\rm\AA$.

In Section 3.1 we have shown that the modulation in chromospheric activity during the stellar cycle can modify the EW of absorption features of a multiplicative factor that depends on $\langle$EW$\rangle$. Stars with typical log R${}^{\prime}_{\rm HK}>-$5.0 dex are affected by this phenomenon.

As a consequence, the stellar parameters inferred from the simultaneous search for three spectroscopic equilibria of iron lines (i.e., excitation equilibrium, ionization balance, and the relation between log N_FeI and the reduced equivalent width EW/$\lambda$) can also be indirectly influenced by the stellar chromospheric activity. This is clearly visible in the four panels of Fig. 3. Each panel shows the sensitivity of the four atmospheric parameters (T_eff, log g, [Fe/H], and $\xi$) to the variation in chromospheric activity. For example, $\partial$T_eff/$\partial$log R${}^{\prime}_{\rm HK}$ as a function of $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$ for all stars in our sample. A negative $\partial$T_eff/$\partial$log R${}^{\prime}_{\rm HK}$ means that the T_eff value determined though our analysis for a particular star decreases as a function of chromospheric activity during the stellar cycle, while a $\partial$T_eff/$\partial$log R${}^{\prime}_{\rm HK}$ consistent with zero indicates that the T_eff value has not changed during the stellar cycle.

As shown before, chromospheric activity can induce an increment of lines with typically large EWs. Therefore, higher $\xi$ values are required to balance the relation between abundances and EW/$\lambda$. The effect becomes more prominent at larger $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$. This explains why $\partial\xi$/$\partial$log R${}^{\prime}_{\rm HK}$ increases with $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$, as observed in the lower-right panel of Fig. 3. An increase of $\xi$ delays the saturation of the curve of growth of each absorption line and, as a consequence, decreases the inferred abundance of the corresponding element (Gray, 1992). This effect is clearly visible in the lower-left panel of Fig. 3, where $\delta$[Fe/H]/$\delta$log R${}^{\prime}_{\rm HK}$ decreases with the increase of the stellar $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$. The effect described in Section 3.1 also impacts the spectroscopic determination of T_eff (see upper-left panel of Fig. 3), but it does not significantly affect log g (see upper-right panel of Fig. 3).

A detailed inspection of Fig. 3 can provide insights on how the chromospheric activity of stars affects our ability to infer T_eff, [Fe/H], and $\xi$ from stellar spectra. With this aim, we simultaneously model the $\partial$P_i/$\partial$log R${}^{\prime}_{\rm HK}$ - $\langle$log R${}^{\prime}_{\rm HK}\rangle$ relations for the i^th atmospheric parameters P_i through Markov-chain Monte Carlo simulations. For the procedure we adopt a model that switches between a null dependence from stellar activity at low log R${}^{\prime}_{\rm HK}$ values to a linear dependence $\partial$P_i/$\partial$log R${}^{\prime}_{\rm HK}$ from stellar activity at high log R${}^{\prime}_{\rm HK}$ values. The model is described as it follows:

where x_j is the $\langle$log R${}^{\prime}_{\rm HK}$$\rangle$ value of the j^th star and $\tau_{i}$ is the switchpoint of the i^th parameter. The model assumes priors for a_i and $\tau$_i that are Normal distributions $\mathcal{N}$($\mu$,$\sigma$), where $\mu$ is the mean and $\sigma$ the standard deviation. Namely, the priors for a${}_{\rm T_{eff}}$, a_[Fe/H], and a_ξ are $\mathcal{N}$(5$\times$10³ K, 10$\times$10³ K), $\mathcal{N}$(0.0 dex, 3 dex), and $\mathcal{N}$(1 km s^-1, 5 km s^-1), respectively. The prior for $\tau_{i}$ is $\mathcal{N}$($-$5 dex, 1 dex). We also assume that data points have Gaussian uncertainties of variance s${}_{ij}^{2}$ + $\zeta_{i}^{2}$, where s_ij is the uncertainty in $\partial$P_i/$\partial$log R${}^{\prime}_{\rm HK}$ for the j-star and $\zeta_{i}$ is a parameter that accounts for the possibility that s_ij are underestimated and that HARPS observations have not homogeneously sampled the entire stellar cycles. The priors for $\zeta_{i}$ is a half Cauchy function with $\gamma$ parameter equal to 10 K, 1 dex, and 10 km s^-1 for T_eff, [Fe/H] and $\xi$ respectively. We ran the simulation with 10,000 sample, half of which are used for burn-in, and employing the No-U-Turn Sampler (Hoffman & Gelman, 2011). The script was written in Python using the pymc3 package (Salvatier et al., 2016).

The convergence of the simulations have been checked by inspecting the traces for each parameter and their autocorrelation plots. The 90$\%$ confidence intervals of the posteriors are listed in Table 6: they are well within the ranges allowed by the priors. The 90$\%$ confidence intervals of the models resulting from the inference are represented in Fig. 3 as blue areas.

In our model, the spectroscopic determination of the P_i parameter is influenced by activity cycle only if the star has a $\langle$log R${}^{\prime}_{\rm HK}\rangle$ $\geq$ $\tau_{i}$. The $\tau$ switchpoint inferred for T_eff occurs with the highest probability at $\sim$$-$4.5 dex, which corresponds to solar twins of 1 Gyr (Lorenzo-Oliveira et al., 2018). While the switchpoint inferred for [Fe/H] and $\xi$ happens at $\sim$$-$5.0, meaning that the determination of these parameters of solar twins younger than 4-5 Gyr is possibly affected by magnetic fields or stellar spots. However, it must be noted that these results are heavily dependent on the master list of absorption features employed in the spectroscopic analysis. In fact, a reduced use of absorption lines with $\langle$EW$\rangle$$\leq$50 m$\rm\AA$ in the spectroscopic analysis would mitigate the impact of magnetic fields or spots on the determination of stellar parameters. However, in this case it would be more difficult to probe microturbulence due to a lack of strong lines.

Finally, we can predict how a star identical to the Sun (i.e., T_eff=5770 K, [Fe/H]=0.0 dex, $\xi$=1.0 km s^-1) would look as a function of its chromospheric activity if analysed with the same linelist employed in this study. To do so, we assume that the variation of the i^th parameter for the j^th star $\Delta$P_ij is equal to zero if the stellar activity index is smaller than the switchpoint (i.e., x_j$<$$\tau_{i}$), while for larger activity indexes it is equal to the integral of Eq. 1. Our model also includes a second term that depends on $\zeta_{i}$ and that accounts for the possibility that the measured uncertainty $\partial$P_i/$\partial$log R${}^{\prime}_{\rm HK}$ is underestimated. Therefore:

The $\Delta$P_i probability distribution at each log R${}^{\prime}_{\rm HK}$ has been inferred from the a_i, $\tau_{i}$, and $\zeta_{i}$ posteriors resulting from the MCMC simulation. In Fig. 4 we show the 90$\%$ confidence intervals of these distributions. The blue areas highlight the range in log R${}^{\prime}_{\rm HK}$ values covered by our dataset, while the red areas are an extrapolation at larger activity indexes of the model fitted to the observations. The maximum variations in stellar parameters traced by blue areas are around 100 K, 0.06 dex and 0.35 km s^-1 for T_eff, [Fe/H], and $\xi$, respectively. These values are similar to the typical uncertainties in atmospheric parameters provided by large spectroscopic surveys observing at optical wavelengths, such as the Gaia-ESO and GALAH surveys (Smiljanic et al., 2014; Buder et al., 2018). However, these surveys targeted pre-main-sequence clusters and stars that are significantly more active than those analysed here. Thus, even if spectroscopic surveys use different methods of analysis and linelists than those employed in our study, it is possible that the magnetic fields or star spots have affected to a certain extent some scientific outcomes of these collaborations (see Section 4 for a discussion on the scientific implications of our result).

Since our dataset samples the artificial variation of stellar parameters up to activity indexes of log R${}^{\prime}_{\rm HK}$$\sim$$-$4.3 dex, the red areas in Fig. 4 represent an extrapolation of this phenomenon at higher indexes, up to log R${}^{\prime}_{\rm HK}$$\sim$$-$4.0 dex, values that are typically measured in star forming regions (Mamajek, 2008). Interestingly, $\xi$ values measured for members of these young associations are typically around 2 km s^-1 (e.g., James et al. 2006; Santos et al. 2008), in agreement with the prediction in Fig. 4 - right panel. This indicates that magnetic fields and stellar spots play an important role in shaping atomic lines in stellar spectra also at the high activity regimes typical of pre-main-sequence stars. Also the prediction of a T_eff variation of $\sim$600 K at log R${}^{\prime}_{\rm HK}$$\sim$$-$4.0 dex is partially consistent with the results of other studies that compared spectroscopic and photometric T_eff values for stars in pre-main-sequence stars (Morel et al., 2003; Baratella et al., 2020) finding differences up to 400 K. However, these studies are not able to probe the effect of stellar spots, which have a roughly similar quantitative effect on the temperatures derived from spectroscopic and photometric data (see Fekel et al. 1986; Morel et al. 2003). Thus, the question around the full effect of stellar activity on spectroscopic T_eff of T Tauri stars is still open. For instance, it is possible that in extremely active stars (e.g., log R${}^{\prime}_{\rm HK}\sim$-4.0 dex) the Zeeman effect and the consequences of cold stellar spots are counteracted but not fully compensated by the indirect effects that magnetic fields have on the stellar photosphere, which allow stellar spectra to probe into deeper layers where temperatures are higher. This could be possible due to the significantly higher strength of magnetic fields on T Tauri stars compared to other objects older than a few 10 Myr.

The selective intensification of the strongest absorption lines in stellar spectra affects the determination of chemical abundances in two opposite ways. On one side, the EW intensification produces an increase of chemical abundances. On the other hand, the growth of $\xi$ has the effect of lowering the chemical abundances. Therefore, the net effect on abundances depends on the relative importance of these two reactions to stellar activity, which in turn depends on the linelist that is used for the analysis.

In Fig. 5 we show the dependence of elemental abundances from stellar activity $\partial$[X/H]/$\partial$log R${}^{\prime}_{\rm HK}$ as a function of $\langle$log R${}^{\prime}_{\rm HK}\rangle$. While each black dot represents one star of the sample, the red squares and their error bars correspond to the binned averaged $\partial$[X/H]/$\partial$log R${}^{\prime}_{\rm HK}$ value and its standard deviation at different $\langle$log R${}^{\prime}_{\rm HK}\rangle$ intervals. These values are weighted by the $\partial$[X/H]/$\partial$log R${}^{\prime}_{\rm HK}$ uncertainties. The general behaviour traced by the red symbols shows that the growth in $\xi$ compensates the abundance increase for all the elements. Similarly to what we have observed for [Fe/H] (see Fig. 3), the variation of $\xi$ along the cycle of the most active stars is large enough to lower the abundances of Si, Sc, Ti, Mn, Ni, and Cu. This is not surprising, in fact, according to our linelist, these elements are mostly determined through the measurement of absorption lines that form deep into the stellar photosphere (i.e., $\tau_{\lambda}$$>$1) and that are not significantly intensified by the stellar magnetic fields or cold spots. Furthermore, most of these lines are medium-strong (e.g., EW$\sim$40-70 m$\rm\AA$) and sensitive to $\xi$. Therefore, for these elements the increase of line EWs is not enough to compensate for the growth in $\xi$.

There is a second class of elements formed by C, Na, Al, S, V, Co, Y, and Zr for which the lines used are also formed deep in the atmosphere, but these lines are so weak that they are not sensitive to $\xi$. Therefore, the change in the abundances of these elements as a function of $\langle$log R${}^{\prime}_{\rm HK}\rangle$ is nearly zero.

Finally, we identify a third group of elements (i.e., Mg, Ca, Cr, and Ba) that are detected mostly through strong lines formed in the upper atmosphere. These lines are both sensitive to $\xi$ and intensified by the stellar magnetic fields or cold spots, in a way that one effect counterbalance the other. Therefore, their change in abundances is also nearly zero.

The chemical analysis of young ($\lesssim$100 Myr) stars is extremely important in the context of Galactic chemical evolution as it provides strong constraints to the models of stellar nucleosynthesis. In addition, young stars have not had time to disperse along the Galactic disk, therefore their chemical content is representative of the ISM’s composition at the location where they are observed today. In contrast to models of Galactic chemical evolution (e.g., Minchev et al. 2013; Sanders & Binney 2015; Frankel et al. 2018), independent studies have consistently found that the youngest stars in our Galaxy have a metal content lower than the Sun. For instance, the metal content measured in different star forming regions located within 500 pc from the Sun (i.e., Chamaeleon, Corona Australis, Lupus, Orion Nebula Cluster, Rho Ophiuchi, Taurus) is on average equal to [Fe/H]=$-$0.07$\pm$0.03 dex (Cunha et al., 1998; Santos et al., 2008; D’Orazi et al., 2011; Biazzo et al., 2011a, 2012; Spina et al., 2014, 2017), which is $\sim$15$\%$ lower than the Solar metallicity. Interestingly, members of star forming regions are found to have $\xi$ values significantly higher than that of the Sun and typically within the range 2.0-2.5 km s^-1 (e.g., James et al. 2006; Santos et al. 2008; Baratella et al. 2020). The typical log R${}^{\prime}_{\rm HK}$ index of these stars is $\sim$$-$4.0 dex (Mamajek & Hillenbrand, 2008). The comparison of these values to the prediction in Fig. 4 suggests that the youngest stars of the solar vicinity appear to be artificially metal poor and with high $\xi$ as a consequence of their high activity levels (e.g., strong magnetic fields, high coverage of cold stellar spots) that have selectively intensified the strongest atomic lines. According to our predictions, these star forming regions should actually contain the same amount of metals of the Sun.

The [Fe/H]-log R${}^{\prime}_{\rm HK}$ prediction in Fig. 4 also provides a convincing explanation for the metal content measured in the sub-clusters of the Orion association. Orion is one of the nearest regions (d$\sim$350-450 pc) of ongoing star formation where both low and high mass stars are formed. It is a complex composed of different sub-clusters with different ages, as the stellar formation burst has spread across the association during the last 20 Myr, triggered by supernovae explosions (Bally, 2008). Since type II supernovae are sites of major nucleosynthesis, these explosions may also chemically enrich parts of the surrounding interstellar gas, and hence the newly formed next generation of stars (e.g., Reeves 1972; Cunha & Lambert 1992, 1994). Therefore, one would expect that the sequential star formation occurring in Orion should result in a peculiar chemical enrichment with the youngest regions being enhanced in metals relative to older ones. Instead, the Orion Nebula Cluster, which is the youngest region of the Orion association, has been found to be the metal-poorest (D’Orazi et al., 2009; Biazzo et al., 2011a, b). Our analysis suggests a revised analysis of Orion’s metal content that should properly take into account the effects of magnetic fields and cold spots on stellar spectra. Providing the first evidence of self-enrichment in a young stellar association would give fundamental insights into stellar nucleosynthesis and, most importantly, on the role that supernovae explosions have in the sequential collapse of molecular clouds, hence on the origin of stars and stellar clusters.

The so-called $barium$ $puzzle$ is still one of the most debated open questions around the production of $s$-process elements in the Milky Way. It originated when D’Orazi & Randich (2009) measured [Ba/Fe] ratios in young ($<$50 Myr) open clusters in the solar vicinity ($<$500 pc) which showed to have $\sim$0.3 dex higher Ba than the value predicted by models of stellar nucleosynthesis (Travaglio et al., 1999; Busso et al., 2001). Further, in contrast to the anomalous Ba overabundance, the abundances of other $s$-process elements such as Y, Zr, La and Ce relative to Fe was found to be Solar (D’Orazi et al., 2012). Interestingly, further studies in young open clusters have shown that additional channels of nucleosynthesis, such as the intermediate neutron-capture process (Cowan & Rose, 1977), cannot explain the Ba overabundance compared to other neutron-capture elements (Mishenina et al., 2015).

A new piece of this puzzle was provided by Reddy & Lambert (2015), who analysed five young (5-200 Myr) local associations and found that they cover an abnormally large range [Ba/Fe] ratios, from +0.07 to +0.32 dex. A further analysis of solar twin stars by Reddy & Lambert (2017) finally provided some clues to the solution of the $Ba$ $puzzle$. Namely, they showed a trend of increasing abundances from the Ba II 5853 $\rm\AA$ line with stellar activity among coeval stars. Therefore, they speculated that the high Ba abundance measured in young associations is not nucleosynthetic in origin but associated with the level of stellar activity. Specifically, they argued that a $\xi$ value derived from Fe lines that form at much larger $\tau_{\lambda}$ in the photosphere is not sufficient to represent the true broadening imposed by the turbulence of the upper photospheric layers where the Ba II lines form.

Our analysis proceeds in this framework, demonstrating that a relation exists between the EW of lines formed at small $\tau_{\lambda}$ into the stellar atmosphere, such as the Ba II lines, and the stellar activity (see Fig. 1, top and middle panels). This dependence is visible in our data only for stars with log R${}^{\prime}_{\rm HK}\gtrsim$-5.0 dex (Fig. 1, bottom panel).

However, in apparent contradiction to Fig. 1, Fig. 5 does not show any clear evidence of a systematic positive dependence of Ba abundances from stellar activity. This is not surprising because, as we pointed out in Section 3.3, the increase in $\xi$ as a function of stellar activity has the important consequence of lowering chemical abundances that are derived at high activity levels. Therefore, in our analysis, the increase in $\xi$ is large enough to counterbalance the effect that the Zeeman broadening and stellar spots would have on Ba abundances.

Since the increment in $\xi$ is highly dependent on the line list employed in the analysis and in particular to the number of lines with high EW/$\lambda$ that can trace the turbulence in the upper stellar layers, the use of different line lists can result in a different sensitivity of the $\xi$ parameter to stellar activity. Therefore, a different line list can also produce very different abundances of Ba, whose lines are extremely sensitive to $\xi$. This explains the large spread in Ba abundances found in young nearby associations by different teams (D’Orazi & Randich, 2009; Reddy & Lambert, 2015, 2017). For instance, the $\xi$ applied by Reddy & Lambert (2017) for the calculation of Ba abundances was taken from Nissen (2015) and based on a list of weak Fe (EWs$\leq$70 m$\rm\AA$) formed quite deep in the atmosphere. Therefore, they were using $\xi$ values that do not reflect the extra broadening of absorption lines in the upper layers of active stars, where Ba lines are formed. On the other hand, our line list includes Fe lines with formation depths similar to those of the Ba lines. This explains the apparent contradiction of the large Ba abundances obtained by Reddy & Lambert (2017) with the lack of Ba variation in Fig. 5. In a similar way, while our masterlist contains 78 Fe I lines, the one used by D’Orazi & Randich (2009) was probably too small (only 33 Fe I lines) to adequately probe the turbulence in the upper stellar layers. In fact, while KG-type stars younger than 50 Myr have $\xi$ that are typically greater than 1.5 km s^-1 (see also D’Orazi et al. 2012), the $\xi$ values estimated in D’Orazi & Randich (2009) are within 0.7 and 1.2 km s^-1 and very close to their first guess values. In conclusion, the Zeeman effect or stellar spots have intensified the lines used to determined the Ba abundance by D’Orazi & Randich (2009), but in their analysis this effect was not counteracted by any $\xi$ increase as in our analysis, leaving the Ba abundances anomalously high.

Do stars swallow their own planets? The major consequence of planet engulfment would be a chemical enhancement of the host star due to the pollution of rocky material. If the accreting star has a sufficiently thin convective zone, the planetary material is not too diluted and can produce a significant increase of the atmospheric metallicity, which can be reliably detected (Spina et al., 2015; Church et al., 2019). In fact, such dilution will not yield an indiscriminate abundance rise of all the metals, but likely will produce a characteristic chemical pattern that mirrors the composition observed in rocky objects with mostly refractory elements (i.e., those with higher condensation temperature) being overabundant relatively to volatiles (Chambers, 2010).

The chemical signatures of planet engulfment events have been found among members of binary systems (e.g., Ramírez et al. 2015; Teske et al. 2016; Oh et al. 2018; Tucci Maia et al. 2019a; Nagar et al. 2019) and open clusters (Spina et al., 2015, 2018b; D’Orazi et al., 2019). Members of the same stellar association are born at the same time and from the same gas, therefore they should be chemically identical. This indicates that the chemical anomalies found among members of the same association cannot be explained by processes of nucleosynthesis. However, this work poses the suspicion that these chemical anomalies are not actually due to planet engulfment events, but instead caused by different activity levels of the members of the same stellar association.

The red and blue symbols in Fig. 6 represent the average of the $\partial$[X/H]/$\partial$log R${}^{\prime}_{\rm HK}$ values for stars with log R${}^{\prime}_{\rm HK}\in$[-4.7,-4.5] and log R${}^{\prime}_{\rm HK}\in$[-4.5,-4.3], respectively, as a function of the condensation temperature T_cond of the X-element listed in Lodders (2003). From this plot it is evident that some elements are sensitive to the variation in chromospheric activity during the stellar cycle, while others are not. We also observe that there is no clear relation between the sensitivity of an element to the stellar activity and its condensation temperature. Therefore, from these data there is no evidence in support of the possibility that chromospheric activity could mimic signatures of planet engulfment events in the chemical composition of stars. Instead, there are elements suggesting that the chemical anomalies found so far are genuine and not related to stellar activity. For example, a number of chemically anomalous Sun-like stars are older than 5 Gyr (Ramírez et al., 2015; Ramirez et al., 2019; Tucci Maia et al., 2019b), therefore - according to our results - their activity levels are not high enough to produce the observed differences in elemental abundances, i,e., $\Delta$[Fe/H]$\geq$0.05 dex. Furthermore, even when the anomaly has been found among members of young stellar associations, there is no relation between the Fe abundance of the stars and their $\xi$ value (Spina et al., 2015, 2018b; D’Orazi et al., 2019), as one would expect if the chemical anomaly was driven by an especially high (or low) activity level of the anomalous star compared to the other siblings.

Similar to a DNA profile, one could use the individual chemical patterns of stars that are not in clusters today to trace them back to a common site of origin (Freeman & Bland-Hawthorn, 2002). This approach - commonly called “chemical tagging” - is a powerful tool for Galactic archaeology, which aims at recovering the remnants of the ancient building blocks of the Milky Way (e.g., clusters, super-clusters, moving groups) that are now dispersed, reconstructing their star formation history and the migration rate of stars within the Milky Way (e.g., Bland-Hawthorn et al. 2010). In fact, the main motivation behind the large-scale spectroscopic surveys of the current decade (e.g., APOGEE, Gaia-ESO, GALAH; Holtzman et al. 2015; Gilmore et al. 2012; De Silva et al. 2015) has been the acquisition of large and homogeneous sets of spectroscopic data from different environments within the Galaxy to trace its history in space and time.

Regardless of the precision achievable in elemental abundances, the success of chemical tagging relies on the significance of critical factors, including the level of chemical homogeneity within members of open clusters and the chemical diversity between open clusters. Even if it is now established that processes of atomic diffusion (Dotter et al., 2017) and planet engulfment events (Laughlin & Adams, 1997) can imprint chemical inhomogeneities among members of the same stellar association, a growing number of studies based on high-precision analysis of solar twin stars are showing that most cluster members on the same evolutionary phase are chemically identical at the typical precision levels reached by large spectroscopic surveys (Liu et al., 2016b, a; Spina et al., 2018a; Nagar et al., 2019).

On the other hand, our analysis shows that magnetic fields or stellar spots can reduce the chemical diversity between stars of different ages, posing a serious challenge for chemical tagging. In fact, while younger stars in our Galaxy should be chemically richer, they also tend to appear poorer in metals than older stars due to the increasing levels of stellar activity (see Fig. 4 - middle panel). In fact, recent studies have shown how challenging is to reconstruct and reassemble the dissolved stellar associations in the Solar vicinity or identify the dispersed family of open clusters solely based on the chemical composition of stars (e.g., Blanco-Cuaresma & Fraix-Burnet 2018; Ness et al. 2018; Simpson et al. 2019; Casey et al. 2019). Therefore, new methods of spectroscopic analysis that could consider the effect of chromospheric activity in the stellar spectra (e.g., Baratella et al. 2020) would pave the way for chemical tagging in the Milky Way disk.

On the other hand, in the context of chemical tagging, it is also possible to use the effect line intensification on stellar spectra to our advantage. In fact, the evidence that equivalent widths of lines can vary during the stellar cycle of quantities which depend on $\tau_{\lambda}$, allows us to use abundance ratios from lines of the same element formed at different depths in the stellar atmosphere to identify young and active stars in the field or among candidate members of young associations.

It is well known that a star, as it ages, evolves along a determined track in the Hertzsprung-Russell diagram that - to a first approximation - depends on the stellar mass and metallicity. Therefore, if the atmospheric parameters and the absolute magnitude M_V of the star are known with enough precision, it will be possible to determine reasonable estimates of its age and mass (Vandenberg & Bell, 1985; Lachaume et al., 1999). In Fig. 4 we show that Sun-like stars with higher chromospheric activity tend to appear cooler and metal poorer than what they really are. As a consequence, their ages and masses determined though isochrones are systematically overestimated (see also Galarza et al. 2019). This has fundamental implications for all the studies that aim to trace the nucleosynthetic history of elements in the Galaxy through stellar [X/Fe]-age relations and chemical clocks (e.g., Nissen 2015, 2016; Spina et al. 2016b, a, 2018b; Feltzing et al. 2017; Tucci Maia et al. 2016; Bedell et al. 2018).

In the recent years, spectroscopic and photometric Galactic surveys have enabled the computation of three-dimensional interstellar extinction maps thanks to accurate stellar atmospheric parameters and line-of-sight distances. This can be achieved by comparing the observed colours to those computed through the stellar parameters and a set of isochrones (e.g., Schultheis et al. 2015; Schlafly et al. 2017). This technique is extensively used also to infer interstellar extinction, which is particular important for young stars in star forming regions or pre-main-sequence clusters that are still partially embedded in the parental cloud. In fact, knowledge of the bolometric luminosities of stars in these young associations is vital for a broad variety of topics in stellar astrophysics. These include the studies of mass segregation, dynamical and structural properties of stellar associations before their dissolution, age spread in star forming regions, rate of star formation in giant molecular clouds and its dependence on time and stellar mass, the analysis of photoevaporation of circumstellar disks and its impact on planet formation (e.g., Luhman 2008; Sacco et al. 2017; Prisinzano et al. 2019).

Young stars are very active and their spectra are heavily affected by the strong magnetic fields. According to our analysis, the T_eff of these stars is systematically underestimated due to the effect of chromospheric activity (Fig. 4 - left panel), which results in an underestimation of the stellar extinction, hence to an overestimation of the bolometric luminosity. Therefore, depending on the technique of spectroscopic analysis and the employed line list, stellar activity may have affected the results of observational studies of pre-main-sequence populations.

The possibility to disentangle the effect of stellar activity and jitter from the radial velocity modulation of stars is of paramount importance for spectroscopic surveys that aim at planet detection. This can be achieved through different indicators, such as the bisectors of the spectral cross-correlation function (Queloz et al., 2001), H$\alpha$ (Bonfils et al., 2007; Robertson et al., 2014) and log R${}^{\prime}_{\rm HK}$ (Noyes et al., 1984; Delisle et al., 2018). More recently, new approaches using these activity indicators and statistical techniques such as Gaussian Processes (e.g., Haywood et al. 2014; Rajpaul et al. 2015; Jones et al. 2017; Delisle et al. 2018) or Moving Average (e.g., Tuomi et al. 2013) have significantly improved our ability to mitigate the impact of stellar activity on the planetary signal. Recent works have searched for new activity indicators, showing how stellar activity affects spectral lines in different ways (Thompson et al., 2017; Davis et al., 2017; Wise et al., 2018; Zhao & Tinney, 2020), which opens up the possibility of using the wealth of information contained in high-resolution spectra to verify the authenticity of a planetary signal (e.g., Dumusque 2018).

Our analysis has shown that not all lines have the same sensitivity to stellar activity, but that only the strongest ones (i.e., those with EW$>$50 m$\rm\AA$ or with log $\tau_{\lambda}$$<$-1 dex; see also Galarza et al. 2019) vary their EW during the activity cycle (see Fig. 1). Therefore, by masking all the lines for which no variation is expected along the stellar cycle and by detecting all the other lines, it is possible to further reduce the noise from stellar activity. This is a very promising possibility that could significantly increase the potential of high-resolution spectrographs in the hunt for planets around Sun-like stars. This is especially true for the spectrographs that do not observe the Ca lines, such as Veloce-Rosso (Gilbert et al., 2018) and Minerva-Australis (Addison et al., 2019).