© The Authors 2026

1 Introduction

Signals of magnetic activity on stellar surfaces are known to vary over different timescales. This variability can be measured in several ways, such as starspots (e.g. Messina & Guinan 2002), emission in chromospheric lines (e.g. Baliunas et al. 1995), and even from direct measurements of stellar magnetic fields. When measuring the magnetic field directly, we rely on the Zeeman effect on spectral lines, either from intensity spectra or spectropolarimetry (see e.g. Donati & Landstreet 2009, for a review). These different methods can characterise magnetic fields on different spatial scales on cool stars, due to the signal cancellation in the spectropolarimetric observations caused by nearby surface regions with opposite magnetic polarities.

The variation of large-scale magnetic fields has been studied extensively using spectropolarimetric observations. Either by using longitudinal magnetic field measurements or complete magnetic field geometries from Zeeman Doppler imaging (ZDI, see e.g. Kochukhov 2016). Using these methods, both shortterm variations, such as rotational modulation, and long-term evolution, such as polarity reversals of the magnetic field geometry caused by the magnetic cycle (e.g. Boro Saikia et al. 2018; Bellotti et al. 2025b), have been observed.

In addition to the large-scale magnetic field, there are also plenty of unresolved magnetic features on the stellar surface. This small-scale field would also include starspots and faculae, as their size is insufficient to be resolved by ZDI, especially when accounting for a random distribution of magnetic field polarities that will cancel out any polarimetric signals. This means that whenever measuring magnetic fields using polarimetry, contributions from many of the magnetic field structures will not be recovered. In fact, when comparing the total and large-scale magnetic field strengths of Sun-like stars, Kochukhov et al. (2020) found that less than ~10% of the magnetic field strength is recovered from ZDI, showing that the small-scale magnetic field is the dominant component for these stars. It is possible to obtain the total unsigned magnetic field strength¹ on stellar surfaces by applying the Zeeman broadening technique to unpolarised spectra (e.g. Robinson 1980). This is done with a combination of broadening and intensification of spectral lines caused by the magnetic field present on the stellar surface. While not sensitive to the magnetic field direction, this diagnostic provides information about the tangled and small-scale surface magnetic field structures, as the effect is not subject to the same signal cancellation as spectropolarimetric measurements. The temporal variation of the total unsigned magnetic field has, until recently, been studied in less detail compared to the large-scale field as multiple spectra are commonly averaged before performing the Zeeman broadening analysis (e.g. Shulyak et al. 2019; Reiners et al. 2022). Snapshot observations have hinted at statistically significant variations over timescales of weeks to years (e.g. Hahlin et al. 2023; Pouilly et al. 2024). Time-series observations covering both rotational and longer timescales have shown indications of both rotational modulations (e.g. Kochukhov & Lavail 2017; Donati et al. 2023) and long-term evolution (e.g. Bellotti et al. 2023) of the total unsigned magnetic field. These clear detections of small-scale variability have, up until recent work on a set of M dwarfs by Cristofari et al. (2025), been limited to very active M dwarfs. Furthermore, Lavail et al. (2019) found no strong signs of rotational modulation on T Tauri stars. For more moderately active Sun-like stars, Kochukhov et al. (2020) found no significant variations from optical spectra between different observational epochs, and Hahlin et al. (2023) reported variations from snapshot observations over the time span of a few months using H-band observations. Furthermore, Petit et al. (2021) found no rotational modulation of ε Eridani while Lehmann et al. (2015) found a statistically significant evolution of the small-scale magnetic field with an approximately 3-year cycle.

The surface-averaged small-scale magnetic field on the Sun, as calculated by disc integrating the line-of-sight magnetic field from SDO images (see Haywood et al. 2016, for details), shows that it constantly changes due to the rotation and evolution of active regions. Furthermore, the disc-integrated total unsigned magnetic field on the Sun appears to trace effects such as radial velocity variations, making the small-scale field a promising tool to improve our precision in radial velocity surveys. However, only the Sun has a surface sufficiently resolved for the magnetic field strength of individual active regions to be measured. To apply this approach in other stars, the small-scale magnetic field would have to be extracted from disc-integrated stellar spectra.

To better understand how magnetic fields influence stellar properties, it would be interesting to characterise the variation of the small-scale magnetic field in more detail. This could help in better understanding the role of the small-scale field on stars other than the Sun, as well as how it interacts with other phenomena used to study stellar activity. Another point of interest is to look at how different spatial scales of the stellar magnetic field are connected to each other. By quantifying how the large- and small-scale magnetic fields are correlated, we would be able to better understand how the underlying small-scale field structures influence the observed magnetic geometry measured from polarised observations.

To investigate the variability of small-scale magnetic fields, we selected a collection of stars with already well-characterised large-scale magnetic fields from Folsom et al. (2016, 2018). We focus on the stars TYC 6349-0200-1 (HD 358623), HIP 76768 (HD 139751), HH Leo (HD 96064, HIP 54155), and Mel 2-55 (HIP 16908). The selection was based on two criteria. The first was that they showed clear and relatively simple rotational modulation of the longitudinal magnetic field. Secondly, the stars were also selected to cover the complete age range (~10 − 600 Myr) studied in Folsom et al. (2016, 2018). This allowed us to test the possibility of recovering total unsigned magnetic field variations on a range of active stars. As the time series used for this study are generally limited to a few weeks to about a month, this study primarily examines variations over a few rotational phases. While longer variations of the total unsigned magnetic field are obviously of interest to monitor how the small-scale field evolves over a stellar activity cycle (Bellotti et al. 2023; Lienhard et al. 2023, e.g.), the shorter time window should allow us to study the total magnetic field on the stellar surface without significant evolution of the structures present on the surface.

The paper is organised as follows: The spectroscopic observations are described in Sect. 2, including additional data reduction steps. In Sect. 3, the process for determining the total unsigned magnetic field on the stellar surface is outlined and the results for each star are presented. In Sect. 4, we discuss the variability, or lack thereof, in the total unsigned magnetic field measurements. Then the field measurements are examined in the context of other activity indicators and in the context of the variability of other activity indicators. The results of this work are summarised in Sect. 5.

2 Observations

2.1 ESPaDOnS and NARVAL spectroscopy

The observations analysed in this study were obtained with ESPaDOnS² at the Canada France Hawaii Telescope and NARVAL at the Bernard Lyot Telescope. Both spectropolarimeters have a resolving power of ~68 000, cover the wavelength range between 3700 and 10 000 Å, and are capable of observing any of the four Stokes parameters for spectropolarimetric analysis. The data were automatically reduced by the LIBRE-ESPRIT package (Donati et al. 1997) and has been used previously (with the exception of the 2017 data of HH Leo) to measure activity indices (the S index, Hα index, and Ca IRT index), longitudinal magnetic fields, and surface magnetic field distributions using ZDI of the stars (see Folsom et al. 2016, 2018). The information about the observations of each star can be found in Table 1.

The ESPaDOnS observations were obtained during a time window of ~2 weeks, which means that any variation in magnetic field strength or activity signals detected in these time series is unlikely to originate from any long-term evolution of the magnetic field on the stellar surface due to cycles or other forms of activity. The data from NARVAL of HH Leo was obtained over a longer time frame, almost three months in 2015 and about one month in 2017. For this reason, these observations could exhibit a larger evolution of the surface field within the time series. With the exception of HH Leo, the cadence of observations ranges from about two times per night to every other night.

2.2 Telluric correction and continuum normalisation

For the total unsigned magnetic field analysis of most stars, as described in Sect. 3, we utilise a Ti I multiplet with lines located between 9600–9800 Å. Unfortunately, this region is significantly affected by telluric absorption, making the lines challenging to use without some form of telluric removal. We use MOLECFIT (Smette et al. 2015) to remove the telluric contamination. Following the procedure from Kochukhov & Lavail (2017), we modify the fits headers for compatibility with MOLECFIT and perform the telluric removal. For this step, we only input the wavelength region between 9630–9900 Å of the observations.

We also performed an additional continuum normalisation of all spectra using a third-order polynomial. At this point, we also manually investigate the spectra for any issues in the telluric removal. We find that some spectra still have significant telluric residuals in some lines, which results in a rejection of a few spectra for some observation sequences (between 0 and 3 for each star). As the method used for the total magnetic field measurements used in Sect. 3 is dependent on high-quality spectra, we also rejected a handful of observations from the variability analysis with median S/N < 70 in the wavelength region³ around the lines used for magnetic field investigations. This is motivated by the fact that at these S/N levels, the magnetic field results became increasingly unstable depending on the input parameters. In some cases, this could create outliers that bias the investigation of variability. With these rejection steps, all stars and epochs still have a sufficient number of observations to provide a complete phase coverage according to the reported periods. In addition, including the rejected spectra in the analysis appears to have little influence on the overall periodic modulation reported on each star.

3 Total unsigned magnetic field inference

We aim to measure how the properties of the surface magnetic field from disc-integrated intensity spectra vary over time. This property is obtained from the distortion of spectral lines caused by the Zeeman effect. In contrast to polarised observations, the signal in the intensity spectra is insensitive to the orientation of the magnetic field (see e.g. Kochukhov 2021). While this means that we do not gain information about the geometry of the local field, it also means that we do not suffer from any signal cancellation of nearby surface elements. This effectively means we can recover the complete magnetic field strength on the stellar surface. This includes field structures that are not resolvable with ZDI, such as starspots.

We use the magnetic inference code from Hahlin et al. (2023) to determine the total magnetic field. It uses Markov chain Monte Carlo (MCMC) sampling from the Solar Bayesian Analysis Toolkit⁴ (SoBAT; Anfinogentov et al. 2021) to find both magnetic and non-magnetic parameters simultaneously. The non-magnetic free parameters used are the chemical abundance (given as log(N/N_total)), non-magnetic broadening (projected rotational velocity v sin i or macroturbulence v_mac), and radial velocity. For the magnetic field, we use a multi-component model to describe the surface distribution of the magnetic field, this means that the surface magnetic field consists of a fixed set of magnetic regions, evenly distributed over the stellar disc. The total unsigned magnetic field (B_I) at the stellar surface is then determined from, $$B_I=\sum _i{f}_i{B}_i.$$(1)

With f_i representing the filling factor of each magnetic field strength used in the model. We set the strengths of these regions to be steps of 2 kG. Another aspect of this code is that the error can also be used as a free parameter, this is justified by the fact that in high S/N studies using high resolution spectroscopy the residuals are often dominated by systematic effects rather than noise. By including the error as a free parameter, we can give the inference more freedom to find a suitable fit. The consequence of implementing this is primarily more conservative error estimates while the median parameters tend to be unaffected.

As the code has previously been used on individual observations, we have made some modifications in order to treat a time series with more consistency. This involves optionally fixing non-magnetic parameters such as abundance and non-magnetic broadening to not let them vary across observations in the same time-series. As these parameters should be stable, any variation in non-magnetic parameters could hide the magnetic variation. In order to constrain the non-magnetic parameters, we use the time averaged spectra, find the best fit parameters (both magnetic and non-magnetic), and then fix the abundance and non-magnetic broadening parameters for each individual spectrum. In order to limit the number of filling factors used, we use the Bayesian information criterion (BIC; Sharma 2017), to weight the improvement of fit with the increased complexity of adding additional filling factors. The BIC has the following form, $$\mathrm{B}\mathrm{I}\mathrm{C}=-2\ln p(Y\mid \hat{\theta })+d\ln n,$$(2)

where $p(Y\mid \hat{\theta })$ is the best fit likelihood, d is the number of parameters and n is the number of data points. We select the model producing the lowest BIC value with the motivation that any increased complexity beyond the lowest BIC will not have a significant improvement to the fit.

As most of our stars have T_eff in the 4000–5000 K range, we elect to primarily use the magnetically sensitive multiplet of Ti I lines in our study. These lines have several advantages as outlined in Kochukhov & Lavail (2017). First, the lines being from the same multiplet means that the relative oscillator strengths are well-known. Second, there is a large range of sensitivities to magnetic fields, the multiplet has a magnetic null line (a spectral line not affected by the magnetic field) making the lines ideal to constrain non-magnetic parameters simultaneously with the magnetic field.

The synthetic spectral grid is generated with line lists from VALD (Ryabchikova et al. 2015), MARCS model atmospheres (Gustafsson et al. 2008) and the polarised radiative transfer code SYNMAST (Kochukhov 2007). We use stellar parameters T_eff, log g, v_mic for each star given by Folsom et al. (2016) or Folsom et al. (2018) as shown in Table 2. The grid is calculated before the MCMC sampling begins, where we fix T_eff and log g using bi-linear interpolation within the MARCS grid. We generate the spectral grid using abundances in a range of ±0.15 from an initial test fit with a step size of 0.05 and magnetic field strengths in steps of 2 kG. Each generated spectrum is then disc-integrated for different strengths of non-magnetic broadening. Here, the non-magnetic broadening parameter has a step size of 0.5 km s⁻¹ in the grid. The macroturbulent velocity is determined by the empirical relationship from Doyle et al. (2014) (unless it is used as a free parameter) and v sin i-values are taken from Folsom et al. (2016, 2018). While both HIP 76768 and TYC 6349-0200-1 are outside of the validity range of the empirical relationship for macroturbulence, their high v sin i-values mean that v sin i is the dominant process for the line shape. In either case, changing macroturbulence only influences the optimal v sin i-value without any significant effect on the average magnetic field. The resulting median parameter values from the obtained posterior distributions for each star can be seen in Table 3. The fit and corner plot for one case (TCY 6349-0200-1) are shown in Appendix D

3.1 HIP 76768

3.1.1 Time-averaged spectra

We produce a time-averaged spectrum from the ESPaDOnS observations by combining the time series observations described in Sect 2. This was done by cross-correlating each spectra with a reference spectra to identify any radial velocity offsets between different observations. Then, we co-added the spectra, using a weighted mean, into the same wavelength grid by shifting each spectra with its radial velocity. As LIBREESPRIT shifts the spectra into a barycentric rest frame, the cross-correlation has a marginal effect on the output spectrum. It does, however, ensure that we correct for other radial velocity variations in the observations.

Using the stellar parameters from Table 2, we generate a synthetic grid and use MCMC sampling to find the magnetic filling factors, abundance, and rotational velocity of HIP 76768. Comparing different numbers of filling factors, in the case of HIP 76768 we explored all models up to 12 kG. We find that the BIC (Eq. (2)) favours a rather complex model (see Table C.2 for values) consisting of six components with field strengths of 0, 2, 4, 6, 8, and 10 kG. The magnetic field obtained for HIP 76768 from the time-averaged spectrum, labelled as ⟨B_I⟩, was found to be 2.7 kG.

3.1.2 Verification of stellar parameters

In order to verify that our choice of stellar parameters from Folsom et al. (2016, 2018) does not significantly influence our results, we re-run the MCMC inference for the mean spectrum of HIP 76768 while also including effective temperature and surface gravity as additional free parameters. Our initial guess is the median parameters of the original analysis in Sect. 3.1.1. We use uniform priors for both T_eff and log g. The procedure is the same, except that we expand the T_eff and log g grid to the neighbouring points in the MARCS grid. The median values, as well as the difference between the original analysis, can be seen in Table 4.

From the result, we can see that there is a close agreement between the obtained and adopted stellar parameters. Both T_eff, log g, and v sin i are all within the uncertainties reported for HIP 76768 by Folsom et al. (2016). We note that our own uncertainties for T_eff, log g, and ε_Ti are quite low. This is a common occurrence for the abundance parameter in other magnetic field studies (e.g. Hahlin et al. 2023) and likely comes from the fact that all lines investigated here are from the same multiplet and therefore share a similar response to changing stellar parameters. The good agreement between standard spectroscopic analysis methods and methods including a proper magnetic field description has also been found by Cristofari et al. (2023), showing that magnetically sensitive lines can be used for stellar characterisation provided the magnetic field effects are accounted for. The magnetic field value changes by slightly above 0.1 kG, this is a variation of about 1σ compared to the case when T_eff and log g were held fixed. This reflects the results of other investigations of magnetic sensitivity to stellar parameter choice (e.g. Kochukhov et al. 2020; Hahlin & Kochukhov 2022; Hahlin et al. 2023) that find similar variations.

We find that when we fit models with fewer filling factors, obtained parameters can be quite different. For example, when using a model with magnetic field strengths of 0, 2, and 4 kG, the obtained temperature is 4660 K. However, the χ² is a factor of 3.6 larger compared to the 6-component model. This shows that the model is significantly worse at fitting the data. In addition, if the magnetic field is not included in the analysis of magnetically sensitive lines, there is a significant risk of introducing systematic errors into fundamental stellar parameters.

While we are able to obtain non-magnetic parameters self-consistently with our MCMC inference, we elect not to use it for the entire sample. The reason is the increased computational cost, in part due to adding two extra dimensions to the spectral grid. In addition, we found that the convergence time of the MCMC sampling for adding T_eff and log g significantly increased compared to when they were held fixed. Given that our results for HIP 76768 are consistent with literature values, we therefore elected to use the established values for these parameters for our analysis.

3.1.3 Individual observations

With the non-magnetic median parameters from Table 3, we then determine the magnetic parameters of each individual observation with the non-magnetic parameters fixed see Table C.1. For consistency, we use the same number of filling factors as in the mean spectrum case. The field strengths obtained for each observation lies in the range 2.6–3.3 kG, this variation is significant compared to the ~0.1 kG uncertainties which indicates that there are substantial variations in the disc averaged magnetic field strength over a time period of a few weeks.

3.2 TYC 6349-0200-1

The analysis of TYC 6349-0200-1 is carried out with the same procedure as described for HIP 76768 in Sect. 3.1. For this star, we find a field strength from the time averaged spectra of around 1.2 kG with individual measurements in the range between 1.1–1.3 kG, this variation does not appear to be significant. Besides the weaker field, we also find that a simpler model is more favoured for the filling factors, specifically a model with 4 components with strengths of 0, 2, 4, and 6 kG.

The intermediate filling factor is consistent with 0 (see Fig. D.1). While this could indicate a surface with highly variable magnetic field strengths, the correlation between f₄ and the other filling factors also indicates that there is cross-talk between the filling factors. In fact, the rotational broadening of TYC 6349-0200-1 is significantly larger (~15 km s⁻¹) compared the broadening of spectral lines caused by the Zeeman effect (~5 km s⁻¹ for a 2 kg field). This means that we are more sensitive to intensification effects related to the strength of the lines, rather than the detailed line shape, making individual filling factors more difficult to infer. While not favoured by the BIC (see Table C.2), we find a similar behaviour from the three-component model where the f₂−component is consistent with zero and most of the contribution comes from the f₄−component. This further highlights the difficulty to disentangle the filling factors at larger rotation rates as different choices of filling factor can produce very different results without substantially changing the magnetic field.

3.3 Mel 25-5

For Mel 25-5, we elect to use v_mac as the free parameter for non-magnetic broadening. The reason for this is that Mel 25-5 has the lowest v sin i in the sample at 3.28 km s⁻¹ (Folsom et al. 2018), this means that the rotational broadening is likely comparable to the broadening caused by macroturbulence. As v sin i is more closely tied to the geometry and period of the star, we elect to keep it fixed at 3.28 km s⁻¹ from Folsom et al. (2018). Instead, we let v_mac be the free parameter. Other than that, the analysis is performed in the same way as described in Sect. 3.1. For this star, we find a magnetic field strength of ~0.7 kG using a multi-component model with 0, 2, 4, and 6 kG. The individual measurements have a variation between 0.6–0.8 kG, which is a bit above two σ between peak to peak.

3.4 HH Leo

HH Leo is the hottest star in the sample by about 500 K compared to the other stars. This is a problem when using the Ti I lines for magnetic inference. While these lines have been used extensively on stars with spectral types M and K, they are not commonly used for hotter stars. In the case of HH Leo the lines are much weaker compared to the other stars in this study. For this reason, we instead choose to use the magnetically sensitive Fe I lines identified by Kochukhov et al. (2020). As these lines are located at around 5500 Å within the optical, there is no significant telluric contamination, which means that we do not need to process the spectra using MOLECFIT. Similarly to the Ti I lines, these lines originate from the same multiplet, have a range of magnetic sensitivities, and include a line that is insensitive to the magnetic field.

Besides the choice of magnetically sensitive lines, the star was also observed twice. One of the datasets is from 2015 and was used in Folsom et al. (2018) to measure the large-scale magnetic field. Another, unpublished data set from 2017 was also available. To explore variations on larger timescales, we analyse these two sets separately by making a time-average spectrum of each epoch individually and then perform two independent analyses with these two averaged spectra. This yields field strengths from the time-averaged spectra of 0.43 and 0.50 kG for the 2015 and 2017 epochs, respectively. In both cases the favoured model contained three components with strengths of 0, 2, and 4 kG. Given the uncertainty of ~0.05 kG, this epoch variation is slightly significant. For the individual measurements, both epochs have peak-to-peak variations of ~0.2 kG around the corresponding average result. While this indicates potential variation over the two years, the scatter of individual measurements is sufficiently large to create a substantial overlap. We also note that, even though the two epochs are fitted independently, the non-magnetic parameters obtained from the two epochs are in agreement with each other.

One notable aspect with the HH Leo results is that there is little variation in the BIC-values (see Table C.2). When testing other model comparison parameters (AIC and WAIC, see Sharma 2017), we find that there is not a consistent agreement of which model is preferred between them. This indicates that neither the two- or three-component is strongly preferred. We elect to persist with the BIC criterion posed for the other stars for consistency, and note that the obtained ⟨B⟩-values for the two-component model is within 0.01 kG of the three-component model. The model choice would therefore have a negligible impact on the result. One benefit with the three-component model is that the degeneracy between the magnetic and non-magnetic parameters are reduced, without changing the median values.

Fig. 1 Magnetic field periodicity of HIP 76768. Left: total unsigned magnetic field variation of HIP 76768 phase-folded using the rotational period from Messina et al. (2010). The dashed red line represents the best-fit sinusoidal function to the data with the fixed rotational period from Table 2. The dashed grey line represents the average of the individual measurements, while the dotted line represents the total unsigned magnetic field obtained from the time-averaged spectrum. Right: lomb–Scargle periodogram obtained from the magnetic field measurements of HIP 76768. The dashed lines indicate false alarm probabilities. Also shown is the expected period from Table 2. Similar plots for all stars can be seen in Fig. A.1.

4 Discussion

4.1 Total unsigned magnetic field variations

As some of the stars exhibited statistically significant variations in the total unsigned magnetic field measurements, it is worthwhile to explore this in more detail. To check for rotational modulation, we fit a sinusoidal function to our observed data where the frequency is fixed by the rotational period of the star from the references in Table 2. To see if the period is recoverable on its own, we also calculate a Lomb–Scargle periodogram (Zechmeister & Kürster 2009) on the same time series.

4.1.1 HIP 76768

The phase-folded result and the sinusoidal fit to the HIP 76768 magnetic field strength can be seen in Fig. 1. From this, we can clearly see that the measured magnetic field strength follows the rotational phase of the star. When comparing the variation with the field from the average spectra, we see that the average spectra produces a slightly weaker field strength compared to most of the individual measurements, although still within the observed variation. When looking at the fit quality, we find that the reduced-χ² is around 0.7, this could indicate that the sinusoidal function is over-fitting the trend, or that we have slightly overestimated the uncertainties by making the observational errors a free parameter in our inference.

This periodicity is further supported by the Lomb–Scargle periodogram shown in Fig. 1. Here we see a peak in the periodogram at a frequency that matches the period from Messina et al. (2010) adopted by Folsom et al. (2016). The other strong peaks appear to be multiples of the frequency corresponding to the strongest peak (i.e. harmonics). This shows that, while the sinusoidal fit might be slightly over-fitted, there is a clear preference for the magnetic field to vary with a frequency corresponding to the rotational period of the star. As these observations were obtained within about 12 days, or three rotational cycles, it is unlikely that any other significant magnetic field evolution would have happened during the observations. This not only points to the fact that the total magnetic field can be significantly different in different regions on the stellar surface but also that the rotational period can be recovered self-consistently from the magnetic field measurements.

4.1.2 TYC 6349-0200-1

The phase-folded results of TYC 6349-0200-1 can be seen in Fig. A.1. As mentioned in Sect. 3.2, the variation of the individual measurements appeared not to be significant. This results in a phase-folded variation that does not show any strong rotational modulation. Upon closer investigation, it appears that the sinusoidal variation that we do see is primarily driven by the outlier point close to the 0.5 phase in Fig. A.1. When that point is manually removed, the other data shows no sinusoidal signal as a function of phase. This is supported by the periodogram in Fig. A.1, that shows no strong peak near the rotational period reported by Folsom et al. (2016). Its lack of rotational modulation is similar to the results from Lavail et al. (2019) that found no rotational modulation on T Tauri stars. As TYC 6349-0200-1 is the youngest star in our sample, while showing the weakest magnetic field variation, our results suggest that young stars have a much more even distribution of magnetic fields compared to their more evolved counterparts. This is also supported by Yamashita et al. (2025), which finds lower than expected longitudinal magnetic field variations in young stars compared to the main-sequence trend.

Fig. 2 Total unsigned magnetic field measurements of HH Leo during 2017. The dashed red line marks the field strength obtained from the time-averaged spectra. The horizontal bar represents the rotational period.

4.1.3 MEL 25-5

The phase-folded data of Mel 25-5, along with its sinusoidal fit, can be seen in Fig. A.1. While the observational sampling is biased to phases between 0–0.5, it does appear to be a rotationally modulated variation in the measured magnetic field strengths. However, this signal is much weaker than in HIP 76768 and while there is a peak corresponding to the rotational period in the periodogram in Fig. A.1, it has a large false alarm probability (>0.1). If the tentative rotational modulation is real, a big limitation for this star is likely the relatively large uncertainties of the magnetic field strengths, although our uncertainties might be slightly overestimated due to setting the error as a free parameter, as well as a phase-biased sampling. This indicates that for accurate constraints on the rotational modulation of magnetic fields on Mel 25-5, we would either need observations with higher S/N or a method able to constrain the magnetic field with higher precision. Alternatively, better sampling of the rotational period might be able to mitigate the noise limitations.

4.1.4 HH Leo

We show the phase-folded data of the two observation epochs separately in Fig. A.1. Although variations of a few σ are detected, neither epoch shows any variation that would hint at a rotationally modulated magnetic field. One aspect that could obscure the rotational modulation of these targets is the longer time period over which the observations were obtained. Over a month or longer, cool stars could undergo sufficient changes to magnetically active regions of the surface to change the phase of any rotational modulation. Specifically, Folsom et al. (2018) reported a worse reduced χ² for the 2015 ZDI map that was attributed to some evolution of the large-scale field even when accounting for differential rotation. In addition, the total magnetic field variation observed in the 2017 dataset appears to follow some trend as the strongest values are concentrated during the beginning and end of the observational epoch, as shown in Fig. 2. This could represent the evolution of surface structures, rather than rotational modulation, on the timescales of ~1 month. However, no similar trends could be seen in the 2015 data.

The mean spectra of the two epochs appear to have a slight difference in the measured magnetic field strength at the 1σ-level, as shown in Fig. 3. This variation is well within the scatter of each epoch but could hint at an evolution of the average activity level of HH Leo. Such trends have been detected on similar stars, such as ξ Boo A where a tentative total unsigned magnetic field variation has been reported over a few years (see Kochukhov et al. 2020; Hahlin et al. 2025).

Fig. 3 Posterior distribution of the total unsigned magnetic field strengths obtained from the time-averaged spectra of HH Leo from the two epochs. The distribution from the 2015 data is shown in red, and the shaded region represents the 2017 data. The black and blue lines represent median and 68% credence regions of the 2015 and 2017 datasets, respectively.

4.2 Comparison with other activity measurements

Many studies on the variability of stellar magnetic fields published over the past decades focused on indirect activity indices. While works such as Reiners et al. (2022) have demonstrated a relationship between activity indicators and the total unsigned magnetic field for a large number of stars, the variation for an activity indicator at any given field strength can vary by an order of magnitude or more. Furthermore, their time-averaged spectra do not indicate if variations in an activity indicator correspond to a variation in the underlying magnetic field strength.

For the longitudinal magnetic field (B_ℓ) and results from ZDI, similar results are reported by Marsden et al. (2014), Vidotto et al. (2014), Brown et al. (2022) and Bellotti et al. (2025a), who find correlations between activity indicators, age, and magnetic field strengths with similar scatter as the total magnetic field measurements. See et al. (2016) have also demonstrated that magnetic cycles determined from magnetic geometry can differ from activity cycles measured with other indicators. From monitoring of individual stars, studies (e.g. Amazo-Gómez et al. 2023; Rescigno et al. 2024) have found a lack of correlation between the longitudinal field and activity indicators. However, Rescigno et al. (2024) found that the RMS of B_ℓ over different time windows correlate well with the solar RV variations, indicating that large variations in B_ℓ on rotational timescales are connected to enhanced stellar activity. While activity and magnetic fields are clearly connected, whether these proxies show coherent variability with the magnetic flux over different timescales remains unclear. Evaluating this connection could help unveil the underlying mechanisms that drive stellar activity.

Figure 4 illustrates the temporal evolution of the discintegrated total unsigned magnetic fields derived in this study, alongside simultaneous activity indices reported by Folsom et al. (2016, 2018), namely the S index from Ca II H&K, the Hα index, and the Ca IRT index from the Ca II infrared triplet. We observe a coherent variability, indicating that the temporal evolution of these activity indicators is closely linked to changes in the small-scale photospheric magnetic field. As highlighted in Fig. B.1, the total unsigned field is positively correlated with the activity traces for all stars but TYC 6349-0200-1. For this target, the Spearman correlation coefficient indicates a positive association only between B_I and the S index.

Figure 5 presents the scatter diagrams of the relationship between magnetic field strength and activity indices, with symbol colours representing the measurements for the four targets analysed in this study. What stands out in this figure is that there is no linear relationship between B_I and the activity indices when considering the full sample. This suggests that, while activity indices can effectively track temporal variations in the magnetic field for a given star, they cannot reliably determine its absolute magnetic field strength. A possible interpretation is that the way magnetic flux is distributed among plages, networks, and spots may significantly impact the S index, as suggested by the solar analysis of Cretignier et al. (2024).

These findings have important implications for studies that compare stars based on their activity indices. For example, Hα and Ca IRT indices would indicate that TYC 6349-0200-1 and HIP 76768 have similar activity levels, whereas HIP 76768 has a significantly stronger magnetic field in reality. Similarly, the S index fails to recognise that HIP 76768 has a stronger magnetic field, erroneously suggesting that TYC 6349-0200-1 should be more magnetic. A similar discrepancy arises between HH Leo and Mel 25-5, where the Ca IRT index incorrectly identifies which star has a stronger magnetic field. It is worthwhile to note that we would see a similar behaviour if using the large-scale magnetic field strength from Folsom et al. (2016, 2018) as the average B_ZDI of our targets (see Table 2) corresponds to between 5 and 3% of our obtained B_I.

Our results clearly indicate that activity indices are powerful tools for tracing magnetic variability over time and are, therefore, well suited to detect magnetic cycles. On the other hand, we also showed that relying solely on activity proxies to compare different stars can introduce significant biases, as they do not provide a comparable measure of the absolute stellar magnetic field strength. This limitation may contribute to the apparent spread observed in empirical relationships such as the magnetic activity–rotation relation (Newton et al. 2017; Astudillo-Defru et al. 2017; Gomes da Silva et al. 2021; Fritzewski et al. 2021; Anthony et al. 2022; Boudreaux et al. 2022), where part of the scatter could be spurious, arising not from intrinsic stellar properties (e.g., stellar mass) but from how magnetic flux manifests differently on stellar surfaces and translates into chromospheric emission. As shown in Fig. B.2, similar conclusions can be drawn when considering how the unsigned magnetic flux, defined as ${\mathrm{\Phi }}_I=4\pi R_{\star }^2{B}_I$, scales with the activity indices, given that the stellar radius does not vary significantly in the sample of stars examined (see Table 2).

Fig. 4 Comparison between the variation of total magnetic fields (left y-axis) and other activity indicators (right y-axis). From top to bottom rows display the results for the four targets studied with increasing age. From left to right columns show the results for the S index, Hα index, and Ca IRT index.

Fig. 5 Correlations between the total unsigned magnetic field and other activity indicators. Top-left: S index. Top-right: Hα index. Bottom-left: Ca IRT index. Symbol colours identify the four targets analysed in this study (see legend).

Fig. 6 Similar to Fig. 4, but for the comparison between the rotational variation of total unsigned magnetic fields and other large-scale polarimetric field measurements from Folsom et al. (2016, 2018). Left-column: Magnetic field strength obtained with Zeeman–Doppler Imaging technique. Right-column: disc-integrated longitudinal magnetic field modulus.

4.3 Magnetic variations on different spatial scales

Figure 6 presents a comparison between the variability of the total unsigned magnetic field and two diagnostics of the large-scale field reported by Folsom et al. (2016, 2018): longitudinal magnetic fields and average large-scale field strengths derived from the visible disc of the ZDI maps, $B_{\text{ZDI}}=\sqrt{B_r^2+B_{\theta }^2+B_{\varphi }^2}$.

We do not find a consistent relationship when comparing the total magnetic field variability with disc-integrated large-scale field strengths derived from ZDI maps at the corresponding rotational phases. Among our targets, one shows no correlation, two show a positive correlation, and one displays an anti-correlation (see Fig. B.1 and the correlation coefficients therein). This lack of a systematic trend suggests that magnetic field strengths obtained from ZDI may not reliably trace variations in stellar activity. This interpretation is supported by Lehmann et al. (2021), who demonstrated that surface-averaged magnetic field strengths obtained from ZDI are not expected to reflect the evolution of the solar activity cycle, particularly when small-scale components dominate the total field.

When comparing the total magnetic field strength with the absolute value of B_ℓ, we find a moderate anti-correlation for HIP 76768 and Mel 25-5, with Spearman coefficients of ρ = −0.6 and −0.58, respectively. One interpretation of this result is that regions with the strongest magnetic field (higher B_I) also possess more complex and tangled field structures, leading to an enhanced cancellation of the polarisation signal and a decrease in |B_ℓ|. Another possibility, as signal cancellation can also be caused by well separated surface features, could be that as more magnetically active regions come into view, increasing the small-scale field strength on the visible disc, the probability of these regions cancelling each other out increases.

5 Conclusion

From this work, it is clear that we can detect rotational modulation of the disc-integrated total unsigned magnetic field on active stars. While this has been reported in the past (e.g. Kochukhov & Lavail 2017; Donati et al. 2023), we have also shown, similar to Cristofari et al. (2025), that the rotational period can be self-consistently recovered from the total unsigned magnetic field variation. Recovery of rotational periods has also been seen using other spectroscopic methods related to activity, such as starspot variability (Tang et al. 2024) and equivalent widths of activity sensitive lines (Gomes da Silva et al. 2025). This apparent variability has some interesting consequences; it is already known that activity variation can influence the determination of accurate stellar parameters (e.g. Spina et al. 2020). The shortterm variations detected here could further complicate things, as accurate determination of stellar parameters would need to be performed simultaneously with an activity analysis to avoid any systematic effects from variability.

When comparing the magnetic field with other activity indicators, we find coherence between the two. Given that there is no clear correlation between the activity and the large-scale field on rotational timescales (see Fig. B.1), it shows that the chromospheric emission is connected to unresolved magnetic regions on the stellar surface. This could be similar to reconnection due to sunspots that drive the chromospheric emission of the Sun. While activity indices appear to follow the magnetic field variation well on individual stars, there is a significant discrepancy between different stars. This means that activity indices appear to be well suited for monitoring activity evolution, but they are less reliable for studying activity relationships of different stars.

One clear limitation found in this work is the precision of magnetic field measurements in the context of studying variability. Even for stars with an age of a few hundred million years, the activity level is not sufficient to clearly detect rotational modulation. While studies in the near-infrared (e.g. Hahlin et al. 2025; Cristofari et al. 2025) have obtained higher magnetic field precision thanks to the lines available in that domain, we would still be unable to investigate stars with activity levels comparable to the Sun. Such work would require the development of more precise methods for the measurement of stellar magnetism, such as an application of the multiline technique described by Lienhard et al. (2023) on stars other than the Sun.

Acknowledgements

We thank the referee for their comments to help improve the paper. A.H. acknowledges support by the Swedish Research Council (projects 2019-03548 and 2023-03667). B.Z. acknowledges financial support from CAPES-PrInt (program #88887.939462/2024-00). C.P.F. acknowledges funding from the European Union’s Horizon Europe research and innovation programme under grant agreement No. 101079231 (EXOHOST), and from the United Kingdom Research and Innovation (UKRI) Horizon Europe Guarantee Scheme (grant number 10051045). K.A. acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP), which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s Excellence Strategy – EXC-2094 – 390783311. Based on observations obtained at the Canada–France–Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada–France–Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site.

References

Appendix A Magnetic variability

Fig. A.1 Magnetic variability for the stars studied in this work. Left column: Phase-folded magnetic field data along with the best fit sinusoidal function with periods from Table 2. The average magnetic field, both from the averaged spectra and the sinusoidal fit are shown as dotted and dashed lines, respectively. Right column: Periodograms of the magnetic sample, red lines marks the literature periods used for the sinusoidal fitting. False alarm probabilities are marked with dashed lines.

Appendix B Activity correlations

Fig. B.1 Correlations between the total unsigned magnetic field, activity traces and other large-scale polarimetric field measurements. The Spearman correlation coefficient is shown in each panel.

Fig. B.2 Similar to Fig. 5, but for correlations between the unsigned magnetic flux and other activity indicators.

Appendix C Individual measurements

Appendix D Representative fit

Fig. D.1 Corner plot from the MCMC inference of the mean spectrum of TYC 6349-0200-1. The contours correspond to 0.5, 1, 1.5, 2σ for a 2D Gaussian distribution. Histogram in the second column shows 95% upper limits of the f4 parameter.

Fig. D.2 Best fit obtained from the MCMC inference of TYC 6349-0200-1.

All Tables

All Figures

Fig. 1 Magnetic field periodicity of HIP 76768. Left: total unsigned magnetic field variation of HIP 76768 phase-folded using the rotational period from Messina et al. (2010). The dashed red line represents the best-fit sinusoidal function to the data with the fixed rotational period from Table 2. The dashed grey line represents the average of the individual measurements, while the dotted line represents the total unsigned magnetic field obtained from the time-averaged spectrum. Right: lomb–Scargle periodogram obtained from the magnetic field measurements of HIP 76768. The dashed lines indicate false alarm probabilities. Also shown is the expected period from Table 2. Similar plots for all stars can be seen in Fig. A.1.

In the text

	Fig. 2 Total unsigned magnetic field measurements of HH Leo during 2017. The dashed red line marks the field strength obtained from the time-averaged spectra. The horizontal bar represents the rotational period.
In the text

	Fig. 3 Posterior distribution of the total unsigned magnetic field strengths obtained from the time-averaged spectra of HH Leo from the two epochs. The distribution from the 2015 data is shown in red, and the shaded region represents the 2017 data. The black and blue lines represent median and 68% credence regions of the 2015 and 2017 datasets, respectively.
In the text

	Fig. 4 Comparison between the variation of total magnetic fields (left y-axis) and other activity indicators (right y-axis). From top to bottom rows display the results for the four targets studied with increasing age. From left to right columns show the results for the S index, Hα index, and Ca IRT index.
In the text

	Fig. 5 Correlations between the total unsigned magnetic field and other activity indicators. Top-left: S index. Top-right: Hα index. Bottom-left: Ca IRT index. Symbol colours identify the four targets analysed in this study (see legend).
In the text

	Fig. 6 Similar to Fig. 4, but for the comparison between the rotational variation of total unsigned magnetic fields and other large-scale polarimetric field measurements from Folsom et al. (2016, 2018). Left-column: Magnetic field strength obtained with Zeeman–Doppler Imaging technique. Right-column: disc-integrated longitudinal magnetic field modulus.
In the text

Fig. A.1 Magnetic variability for the stars studied in this work. Left column: Phase-folded magnetic field data along with the best fit sinusoidal function with periods from Table 2. The average magnetic field, both from the averaged spectra and the sinusoidal fit are shown as dotted and dashed lines, respectively. Right column: Periodograms of the magnetic sample, red lines marks the literature periods used for the sinusoidal fitting. False alarm probabilities are marked with dashed lines.

In the text

	Fig. B.1 Correlations between the total unsigned magnetic field, activity traces and other large-scale polarimetric field measurements. The Spearman correlation coefficient is shown in each panel.
In the text

	Fig. B.2 Similar to Fig. 5, but for correlations between the unsigned magnetic flux and other activity indicators.
In the text

	Fig. D.1 Corner plot from the MCMC inference of the mean spectrum of TYC 6349-0200-1. The contours correspond to 0.5, 1, 1.5, 2σ for a 2D Gaussian distribution. Histogram in the second column shows 95% upper limits of the f4 parameter.
In the text

	Fig. D.2 Best fit obtained from the MCMC inference of TYC 6349-0200-1.
In the text