Benchmarking the fundamental parameters of Ap stars with optical long-baseline interferometric measurements

Context. The variety of physical processes at play in chemically peculiar stars makes it difficult to determine their fundamental parameters. In particular, for the magnetic ones, called Ap stars, the strong magnetic fields and the induced spotted stellar surfaces may lead to biased effective temperatures when these values are derived through spectro-photometry. Aims. We propose to benefit from the exquisite angular resolution provided by long-baseline interferometry in the visible to determine the accurate angular diameters of a number of Ap stars, and thus estimate their radii by a method that is as independent as possible of atmospheric models. Methods. We used the visible spectrograph VEGA at the CHARA interferometric array to complete the sample of Ap stars currently observable with this technique. We estimated the angular diameter and radius of six new targets. We estimated their bolometric flux based solely on observational spectroscopic and photometric data to derive nearly model-independent luminosities and effective temperatures. Results. We extend to 14 the number of Ap stars for which interferometric angular diameters have been measured. The fundamental parameters we derived for the complete Ap sample are compared with those obtained through a self-consistent spectroscopic analysis. Based on a model fitting approach of high-resolution spectra and spectro-photometric observations over a wide wavelength range, this method takes into account the anomalous chemical composition of the atmospheres and the inhomogeneous vertical distribution for different chemical elements. Regarding both the radii and the effective temperatures, the derived values from our interferometric observations and from self-consistent modelling are consistent within better than 2σ for nine targets out of ten. We thus benchmark nine Ap stars for effective temperatures ranging from 7200 and 9100 K, and luminosities ranging between 7 L and 86 L . Conclusions. These results will be key for the future derivation of accurate radii and other fundamental parameters of fainter peculiar stars for which both the sensitivity and the angular resolution of the current interferometers are not sufficient. Within the context of the observations of Ap stars with the Transiting Exoplanet Survey Satellite (TESS), these interferometric measurements are crucial for testing the mechanism of pulsation excitation at work in these peculiar stars. In particular, our interferometric measurements provide accurate locations in the Hertzsprung-Russell diagram for hot Ap stars among which pulsations may be searched for with TESS, putting to test the blue edge of the theoretical instability strip. These accurate locations could be used to derive masses and ages of these stars through a specific grid of models, and to test correlations between the properties of these peculiar stars and their evolutionary state.


Introduction
Magnetic fields are detected in a large variety of stars throughout the spectral type range and are known to play a crucial role during star formation and evolution. They have an impact on numerous physical phenomena at play within the stars and their environments such as accretion, mass loss, turbulence, and pulsations, and on some fundamental stellar properties such CHARA/VEGA observations. as rotational speed and stellar chemical abundances (Donati & Landstreet 2009). One of the best examples of how stellar magnetic fields impact stellar surface properties is found in the chemically peculiar Ap stars. These are stars of spectral types A and B, characterized by strong and large-scale organized magnetic fields typically of several kG, and abundance inhomogeneities of a few orders of magnitude with respect to solar values, leading to spotted surfaces (Preston 1974). They also exhibit low rotational speeds (≤100 km s −1 ) supposed to be due to magnetic braking (Mathys 2017 Notes. The objects previously observed are in italics, and the corresponding references are given in the Notes. References. (a) Perraut et al. (2016); (b) Perraut et al. (2015); (c) Bruntt et al. (2008); (d) Bruntt et al. (2010); (e) Perraut et al. (2013); ( f ) Romanovskaya et al. (2019b); (g) Perraut et al. (2011). stellar surfaces, the determination of their fundamental parameters by (spectro)-photometry can be biased. Using long-baseline interferometry to provide accurate angular diameters of these particular stars appears to be a promising approach to overcoming the difficulties in deriving accurate global parameters. However, this approach is also very challenging because of the stars' small angular size (less than 1 mas). The first interferometric observations of the largest Ap star α Cir by Bruntt et al. (2008) were performed with the VLTI (Glindemann et al. 2003). Then several campaigns on the brightest and the largest Ap stars were conducted with the CHARA array ) since this interferometric array with its baselines B as long as 330 m and its visible instrumentation allows targets as small as λ/2B = 0.2 mas to be partially resolved. To date, accurate angular diameters have been determined by optical long-baseline interferometry for eight Ap stars (given in italics in Table 1). In this paper we report on the angular diameter measurements of six additional Ap stars in the northern hemisphere, which extends the Ap sample for which we derive stellar radii by interferometry to 14 (hereafter referred to as the complete Ap sample). Our objective is to cover a wide range of luminosities and effective temperatures, and thus to probe the theoretical instability strip where the sub-category of Ap stars exhibiting pulsations, the rapidly oscillating Ap (roAp) stars, are expected (Balmforth et al. 2001;Cunha 2002;Cunha et al. 2019). The effective temperatures of our complete Ap sample range from about 7200 to 10 550 K and the luminosities from about 7 to 100 L (Table 1). Ultra-precise photometric light curves will be available for most known Ap stars once the NASA Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2014) completes its first two years of science observations. Having accurate stellar parameters and the necessary data to search for the presence of pulsations in the sample of stars considered in this work is an important step towards testing theoretical predictions for the driving of pulsations in roAp stars (Cunha et al. 2020).
Our aim is to derive stellar radii and effective temperatures that are as independent of stellar models as possible. Once validated these models could be used to derive radii of fainter and/or farther stars for which the sensitivity and/or the angular resolution of interferometry are currently not sufficient. When looking at the Renson catalogue of Ap, HgMn, and Am stars (Renson & Manfroid 2009), most of the targets are located in the southern hemisphere and the histograms of the stars for both hemispheres is centred on Vmag = 9 (Fig. 1, left), which is out of reach of the sensitivity of the current visible interferometric instrumentations whose best performance corresponds to a limiting magnitude of about 8-8.5 in R band (Ireland et al. 2008;Ligi et al. 2013). In addition, when looking at the expected angular sizes of the targets that are brighter than this limit, we find that only a few of them are above 0.2 mas (Fig. 1, right), a value that is already very challenging to measure.
We present the interferometric observations and the angular diameters and radii derived from our last CHARA/VEGA observations in Sect. 2. We then use (spectro-)photometric data to compute the bolometric fluxes based solely on observational data in Sect. 3. In Sect. 4 we compute the stellar effective temperatures and luminosities and discuss our results in Sect. 5.

Last interferometric observations of Ap stars
with VEGA/CHARA

Observations
During the 2016-2018 period, we observed six Ap stars with different telescope triplets of the CHARA array ) and the VEGA visible instrument (Mourard et al. 2009). We operated VEGA in parallel with the K-band CLIMB beam combiner acting as a fringe tracker (Sturmann et al. 2010). We used the medium spectral resolution of VEGA (R ∼ 6000) and recorded fringes on spectral bands centred around ∼700 or ∼730 nm, depending on the targets. We interleaved our target observations with calibrator observations, following the sequence calibrator-target-calibrator, with 40 or 60 blocks of 2500 short exposures (of 10 ms) per star. We swapped from one star to another every 30 min to correctly calibrate the instrumental transfer function. We used the SearchCal software (Bonneau et al. 2006) proposed by the JMMC 1 to find suitable calibrators. The observation log is given in Appendix A along with the spatial frequency coverage for all targets.

Data processing
We used the standard VEGA reduction pipeline (Mourard et al. 2009) to compute the squared visibilities over spectral bands of 15 or 20 nm according to the width of water absorption lines in the spectrum. Such large bandwidths lead to a mean spectral resolution of about 40, and our visibility measurements are thus mainly sensitive to the stellar photosphere. We followed the same processing as described in our previous papers (Perraut et al. 2011(Perraut et al. , 2013(Perraut et al. , 2015(Perraut et al. , 2016 to compute raw squared visibility for each block of 2500 individual frames. Even though we recorded data with telescope triplets, we cannot retrieve closure phases on our targets since the signal-to-noise ratio on this observable is strongly limited by the saturation limit of the VEGA photon-counting detectors (see Fig. 4 in Mourard et al. 2012).

Angular diameters and radii
We followed the same approach as in our previous papers to calibrate the target squared visibilities V 2 by using the angular   (Bonneau et al. 2006;Bourges et al. 2017) to compute the transfer function. We adopted a conservative uncertainty of 7% on the uniform angular diameter (Table 2). Our calibrators are stars of A-B spectral types. For these early-type stars, since the JSDC diameter estimation exhibits discrepancies with the surface brightness-colour relations (SBCRs) for V-K < 0, we adopted the uniform diameters provided by the JSDC2 catalogue which shows a better consistency with SBCRs (see Fig. 3d in Challouf et al. 2014). The calibrated squared visibility curves versus spatial frequencies are given in Fig. 2.
We determined the uniform-disc angular diameter θ UD of our targets in the R band using the model fitting tool LITpro 2 (Tallon-Bosc et al. 2008), and the limb-darkened angular diameters θ LD considering the linear limb-darkening coefficients in the R band from the Claret tables (Claret & Bloemen 2011) when using an effective temperature range in agreement with the values in Table 1 and a surface gravity log g between 4 and 4.5.
Stellar radii R (in solar radius, R ) can be derived from our angular diameter measurements θ LD (in mas), and the parallaxes π P (in seconds of arc) through the formula (1) 2 www.jmmc.fr/litpro_page.htm We used the parallaxes π P from the second Gaia DR2 release (Gaia Collaboration 2018) and computed the error bars on π P using the formula provided by Lindegren et al. (2018) at the IAU colloquium: Here σ i is the internal error provided in the Gaia DR2 catalogue, k = 1.08, and σ s = 0.021 mas is the systematic error for G < 13. All the angular diameters, the parallaxes we consider, and the radii are given in Table 4.

Bolometric fluxes
We follow the same approach as in our previous papers to derive the bolometric fluxes from observational spectral energy distributions (SEDs; e.g. Perraut et al. 2011Perraut et al. , 2013Perraut et al. , 2015Perraut et al. , 2016. So, in addition to the six Ap targets presented in the last section, we also computed the SED of HD 188041 and HD 204411, which was not done in Romanovskaya et al. (2019b).

Building the spectral energy distributions
To compute the bolometric fluxes from the SED, we gathered the following spectro-photometric data: -For wavelengths in the range 1265 Å < λ < 3180 Å, we used rebinned spectra from the Sky Survey Telescope obtained at the IUE Newly Extracted Spectra (INES) data archive 3 . We removed all bad pixels and measurements with negative flux using the quality flag (Garhart et al. 1997); -We collected the photometric data from TD1 (Thompson et al. 1978(Thompson et al. ) at 1565(Thompson et al. , 1965(Thompson et al. , 2365(Thompson et al. , and 2740 -For wavelengths in the visible range we used fluxcalibrated ELODIE (Prugniel & Soubiran 2001;Prugniel et al. 2007) and/or MILES (Sánchez-Blázquez et al. 2006) spectra, along with spectroscopic observations with the Boller & Chivens long-slit spectrograph mounted at the 2.1 m telescope at the Observatorio Astronómico Nacional (OAN) at San Pedro Mártir in Mexico (see Romanovskaya et al. 2019b for the data description and processing). -For wavelengths in the range 3300 Å < λ < 7850 Å we also collected the photometric data from the Adelman catalogue (Adelman et al. 1989) and converted the flux as described in our Appendix B. We used these data to calibrate the visible spectrum of HD 108662 in flux; -At longer wavelengths (in the J, H, and K bands) we collected the photometric data from the 2MASS All Sky Catalog of point sources 4 (Cutri et al. 2003) and corrected for the filter responses using the zeropoints from Cohen et al. (2003).
At the lower wavelength end of the spectral distribution, we performed a linear interpolation on logarithmic scale between by more than 4% over the 3400-3700 Å wavelength range (hereafter called constant flux). This corresponds to the Balmer jump, which is very pronounced for the A-type stars. We noticed that the IUE flux for HD 108662 is too low compared to this Adelman constant flux (Fig. 3, right, first row), which might come from flux losses of the IUE satellite for bright sources. To match IUE flux with Adelman UV flux, the IUE flux needs to be scaled up by a factor of 1.32 (Romanovskaya et al. 2020b).

Estimating the bolometric fluxes
Three targets have a complete SED (HD 108662, HD 108945, and HD 120198) and their bolometric flux values can be determined by computing the area under the curve of this distribution. When spectra and photometric data are both available in a spectral range, we always used the spectra to compute the area under the SED. For the objects for which the (spectro-)photometric data do not cover the full spectral range, we used the spectro-photometry dataset coming from the complete target sample to devise a way to estimate, for each star, the flux in wavelength ranges where data are missing. This was done by establishing the correlations between the different flux data, i.e. the TD1 flux integrated over different wavelength intervals, the Adelman constant flux, and the B and V magnitudes, as described in Appendix B. In practice we proceeded as follows: -for HD 4778, for which no Adelman data are available, we used the linear relation we found between the constant flux and the B magnitude (Fig. B.1) to estimate a constant flux of 0.99 × 10 −11 erg cm −2 s −1 Å −1 (white stars in Fig. 3, top left). Using the linear relation derived from the V magnitude led to the same constant flux within an accuracy better than 4%. We then computed the bolometric flux from flux integration.
-for HD 153882 and HD 204411 the INES spectra are missing, but TD1 data are available (blue triangles in Fig. 3, first column, third and fourth rows, respectively). For all the targets for which we have TD1 and INES spectra, we computed the integrated flux from 1975 to 3180 Å and showed that it strongly linearly correlates with the TD1 flux at 2365 Å (regression coefficient of 0.97; blue symbols and line on Fig. B.2). We used the derived linear law to determine an integrated flux of 7.9 × 10 −9 erg cm −2 s −1 for HD 153882 and of 1.46 × 10 −8 erg cm −2 s −1 for HD 204411 over the [1975 Å; 3180 Å] range. We also estimated this integrated flux by considering the correlation with the Adelman constant flux for the targets for which we have both INES and Adelman datasets (regression coefficient of 0.7; orange symbols and line in Fig. B.2). With this approach, we estimated an integrated flux over the [1975 Å; 3180 Å] range of 9.7 × 10 −9 erg cm −2 s −1 for HD 153882 and of 1.58 × 10 −8 erg cm −2 s −1 for HD 204411. The agreement between the two approaches is better than 8% for HD 204411 and of the order of 20% for HD 153882. In the following we used the flux values obtained with the first method since the linear regression coefficient in that case is very close to 1.
-for HD 220825 the situation is different since we do not have TD1 data and we could not apply the same methods as before. We estimated the integrated flux for the [1975 Å; 3180 Å] range by using its linear correlation with the Adelman constant flux (Fig. B.2). We obtained an integrated flux of 3.1 × 10 −8 erg cm −2 s −1 for HD 220825. We studied for all our targets the correlation between the integrated fluxes over the [1265 Å; 1975 Å]  As in our previous papers, we estimated the uncertainty associated with the bolometric flux by considering the following conservative uncertainties for each contribution: 3% uncertainty on the flux computed from the ELODIE, MILES, or OAN spectrum; 10% on the flux computed from the combined INES spectra; and 15% on the flux derived from interpolations and/or estimations from correlations. These largest uncertainties are on the fluxes over the lowest wavelength range, which generally correspond to about 12-15% of the total flux, and always to less than 25%.

Comparison with atmospheric models
The bolometric fluxes we obtain for our new targets are listed in Table 3. We derived uncertainties on the bolometric fluxes A101, page 5 of 13 A&A 642, A101 (2020) between 7 and 12%. We compared our values with the bolometric fluxes derived from a self-consistent spectroscopic analysis, when available. The latter approach is based on model fitting of spectro-photometric observations over a wide wavelength range, including the same INES, TD1, Adelman, and 2MASS datasets as ours. This method includes the modelling of the anomalous chemical composition of the atmospheres of these peculiar stars and takes into account the inhomogeneous vertical distribution for different chemical elements. Model atmosphere calculations are conducted with the LLMODELS code (Shulyak et al. 2004). Through an iterative process, a theoretical SED is produced and compared to the observed one; the radius is a free parameter Notes. The last three columns give the radii and the effective temperatures as derived from self-consistent modelling and when considering the Gaia DR2 parallaxes, and the magnetic field strengths. Four targets have not been modelled yet. ( * ) For these two targets we keep the parallaxes from HIPPARCOS because of the binarity and/or the excessive brightness. (+) We report the average values of the effective temperatures and luminosities provided in Perraut et al. (2013).
References. together with the surface gravity and the effective temperature; the stellar parallax is used to scale the model flux to the distance of the star to be compared with the observed flux at the Earth's orbit; and the iterations continue until the fundamental parameters used in the SED fitting converge. Thus, the abundance pattern is consistent with the physical parameters of the final stellar atmosphere. This approach is described in detail in several papers (e.g. Shulyak et al. 2009;Romanovskaya et al. 2019b). As shown in Table 3, the two methods provide bolometric fluxes that are in excellent agreement. We use our values, derived solely from data, combined with our angular diameters to compute the fundamental parameters of our Ap targets.

Fundamental parameters
We used our bolometric fluxes f bol (Table 3) and the parallaxes π P from Gaia DR2 (Gaia Collaboration 2018) to determine the luminosities of our targets as where C stands for the conversion factor from parsecs to metres. We used the measured angular diameters and bolometric fluxes to estimate the effective temperatures of our new targets with the equation where σ stands for the Stefan-Boltzmann constant (5.67 × 10 −5 erg cm −2 s −1 K −4 ).
All our results are given in Table 4. We also report the fundamental parameters of the Ap stars that were observed previously via interferometry. We recomputed the values when considering the Gaia DR2 parallaxes, except for the two brightest targets HD 128898 and HD 137909, the latter being a binary star. For both we kept the HIPPARCOS parallaxes (van Leeuwen 2007) that are more accurate than the Gaia DR2 values.

Positions in the Hertzsprung-Russell (HR) diagram
We used our fundamental parameters to put the stars of the complete Ap sample in the HR diagram along with evolutionary tracks from Girardi et al. (2000), when considering a solar metallicity for all the stars (i.e. an initial metallicity Z = 0.019 and an initial helium abundance Y = 0.273), and models with overshooting (Fig. 4). We also overplotted the blue and red edges of the theoretical instability strip (blue dash-dotted lines), as computed by Cunha et al. (2019). Our sample nicely covers this instability strip. It can be clearly seen that none of the roAp stars in our sample (denoted by stars in Fig. 4) is located close to the blue edge of this theoretical instability strip since all of them have effective temperatures log T eff < 3.90. This is in line with the position of the overall sample of roAp stars known to date, as discussed by Cunha et al. (2019), and potentially indicates that the excitation mechanism at work in roAp stars is not fully understood.

Discussion and concluding remarks
Thanks to our interferometric measurements of angular diameters of Ap stars, we can compare the derived interferometric radii with those computed from self-consistent modelling, which are available for ten stars from our complete sample A101, page 7 of 13 A&A 642, A101 (2020) Fig. 4. Complete sample of Ap stars in the Hertzsprung-Russell diagram. Stars denote the known rapidly oscillating Ap stars, while circles denote the Ap stars for which pulsations have not been detected yet. The dashed blue lines denote the blue and red edges of the theoretical instability strip. The thin dash-dotted lines correspond to evolutionary tracks for masses ranging from 1.4 M to 2.5 M and the thick black line denotes the ZAMS. The colour-coding gives the stellar magnetic field strength Bs. The numbers close to the symbols refer to Table 4. ( Table 4). We recomputed these modelled radii using the Gaia DR2 parallaxes for all stars except HD 128898 and HD 137909, as explained in the notes of Table 4. Radii and effective temperatures derived from interferometric observations and selfconsistent modelling are consistent within better than 2σ for nine of the ten targets (Fig. 5). The two radii determinations are in tension at a 4σ level for HD 108662, which is resolved by our observations and for which the scatter of the visibility points is the largest (Fig. 2, top right), i.e. the difference between the measured visibilities and the uniform-disc modelled visibilities can reach up to ±2.6 times the error on our visibility measurements. For this object such a scatter in visibilities remains when considering the 2017 dataset alone or the 2018 dataset alone. When looking at the weather conditions, for all observations, the Fried parameter is quoted to be larger than 9 cm, which corresponds to good conditions for such a source as bright as 5.2 in V and K bands. HD 108662 has the same calibrator as HD 108945, whose visibility curve does not exhibit such a large scatter. All these findings suggest that the observed scatter in the HD 108662 data might have a stellar origin. This star has a known companion at 140 , as revealed in the Washington Double Star catalogue (WDS; Mason et al. 2001). This companion (HD 108651) is a spectroscopic binary. Even if hierarchical systems are commonly observed for these spectral types (HD 108662 is an A0 star), HD 108662A is not known as a spectroscopic binary (no note in WDS; no entry in Pourbaix et al. 2004). To go a step further, we used the CANDID tool (Gallenne et al. 2015) to estimate the 3σ detection limit of a companion within 50 mas around HD 108662A. We defined two sets of observations, one for the 2017 observations (recorded in April and May) and one for the observations recorded in April 2018. Despite the lack of closure phases, we can rule out a companion two magnitudes fainter than the star (i.e. V ∼ 7.2). Based on previous studies, even if HD 108662 is quoted in Simbad as a variable star of α 2 CVn type and, as such, might exhibit stellar spots, it seems difficult to explain such a scatter in visibility as being due to a stellar spot; Shulyak et al. (2014) present several simulations of spotted surfaces of chemically peculiar stars and mention that the impact of a spot becomes noticeable once the squared visibility is below 0.3, provided that the visibility accuracy is better than 1%. In our case, even if the squared visibilities are as small as 0.25 along the longest baselines mainly oriented along a position angle of 100-130 • (see the spatial frequency coverage in Fig. A.1), the visibility differences we observed are much larger than 1%. The too-large radius of HD 108662 we derived from the interferometric angular diameter leads to a too-small effective temperature (see object 3 in Fig. 5) which does not allow the visible spectrum, and especially the H lines, to be correctly fitted (Romanovskaya et al. 2020b).
Benchmarking the Ap stars that have been already selfconsistently modelled is of strong interest for deriving accurate radii, and then other fundamental parameters, of fainter Ap stars for which both the sensitivity and the angular resolution of the current interferometers are not sufficient. As clearly illustrated in this work with HD 153882, one of the faintest and smallest targets in our sample, deriving an accurate angular diameter for targets as small as 0.2 mas is very challenging. We have to cope with calibration issues, and thus with large error bars on the radius and the effective temperatures. Thus, for most of the Ap stars, indirect spectroscopic analysis would be the relevant approach to derive accurate fundamental parameters.
In order to further test the excitation mechanism at work in roAp stars, which is not yet fully understood, it is important to look for pulsations in hotter Ap stars whose location in the HR diagram is known accurately. That is the case of the hotter stars in our sample (indicated by circles in Fig. 4) whose effective temperatures range between 8000 and 10 000 K (i.e. 3.90 < log T eff < 4.0). Seven of these stars will be observed by the TESS satellite during the currently ongoing second year of science operations, and it will be possible to establish if any of these A101, page 8 of 13 K. Perraut et al.: Benchmarking the fundamental parameters of Ap stars with optical long-baseline interferometric measurements Fig. 5. Comparison between the radii (left) and the effective temperatures (right) derived from interferometric measurements and those derived from the self-consistent spectroscopic modelling. The dash lines give the 1:1 relations. The numbers close to the symbols refer to Table 4. hotter stars exhibits pulsations with amplitudes within reach of the TESS observations. A detailed modelling of the stars, aimed at deriving their masses and ages based on a comprehensive grid of models, will be presented in a future work.
With a precise location of the Ap stars in the HR diagram, one can also consider searching for evidence of a correlation between the magnetic field strength and stellar age. This question is still a matter of debate. Looking at the distribution of the Ap magnetic stars in the HR diagram, Hubrig et al. (2000) claimed, for stars with masses smaller than 3 M , that magnetic fields only appear in stars that have already completed at least approximately 30% of their main sequence lifetime. However, this is contradicted by Bagnulo et al. (2003) who have observed the Ap star HD 66318. They report the largest magnetic field detected in an Ap star, but this target has only completed about 10-20% of its main sequence life. This is also in contradiction with the findings of the study of Kochukhov & Bagnulo (2006) on the evolutionary state of these magnetic stars, and with the distribution of roAp stars of a TESS sample presented in Fig. 1 of Cunha et al. (2019), where about 20 targets are very close to the zero age main sequence (ZAMS). Our accurate determination of the fundamental parameters of HD 220825 (object 14) is in direct line with such a distribution and clearly sets this 2 M magnetic Ap star very close to the ZAMS. More recently, based on a study of stars in open clusters, Landstreet (2011) argued that the typical magnetic field declines rapidly as the Ap stars evolve on the main sequence. With this in mind, we collected the surface magnetic field strength of our targets in the literature and reported them in Table 4. Looking at the position on the main sequence of each target with respect to the ZAMS represented by the thick black line in Fig. 4, no apparent correlation emerges in our results; for instance, the youngest star of our sample, HD 220825 (object 14), exhibits a stellar magnetic field smaller than 3 kG, while HD 137909 (object 8) located at the middle of the main sequence has a magnetic field more than two times higher. More generally, a decrease in the magnetic field strength with age would translate into green symbols close to and along the ZAMS, and black symbols farther away at the end of the main sequence. Perhaps we could not see a clear correlation because the number of stars in our study is limited and the star-to-star variations do not allow us to identify a possible statistical correlation between the two properties. It could also be that the correlation found by Landstreet (2011) is not present for stars with masses lower than 2.5 M . The stars in their study were typically of higher mass than the stars in our study, with their lowest mass bin including stars in the range of 2-3 M . In the future, the use of self-consistent atmospheric modelling will allow us to increase the number of stars in the mass range studied here, and thus to verify whether the evolutionary dependence found by the authors is also present for the lower mass Ap stars.