The variable magnetic field of V889 Her and the challenge of detecting exoplanets around young Suns using Gaussian process regression

Discovering exoplanets orbiting young Suns can provide insight into the formation and early evolution of our own solar system, but the extreme magnetic activity of young stars obfuscates exoplanet detection. Here we monitor the long-term magnetic field and chromospheric activity variability of the young solar analogue V889 Her, model the activity-induced radial velocity variations and evaluate the impacts of extreme magnetism on exoplanet detection thresholds. We map the magnetic field and surface brightness for 14 epochs between 2004 and 2019. Our results show potential 3-4 yr variations of the magnetic field which evolves from weak and simple during chromospheric activity minima to strong and complex during activity maxima but without any polarity reversals. A persistent, temporally-varying polar spot coexists with weaker, short-lived lower-latitude spots. Due to their different decay time-scales, significant differential rotation and the limited temporal coverage of our legacy data, we were unable to reliably model the activity-induced radial velocity using Gaussian Process regression. Doppler Imaging can be a useful method for modelling the magnetic activity jitter of extremely active stars using data with large phase gaps. Given our data and using Doppler Imaging to filter activity jitter, we estimate that we could detect Jupiter-mass planets with orbital periods of $\sim$3 d. A longer baseline of continuous observations is the best observing strategy for the detection of exoplanets orbiting highly active stars.


INTRODUCTION
Young solar analogues show extreme magnetic activity that contrasts against the solar case, yet our Sun is thought to have behaved similarly during its infancy.By studying the magnetism of young solar analogues we can look back at the Sun's possible early main sequence life and see how the conditions would have been within our own space environment.The detection of exoplanets around these stars could also give insight into the formation and early evolution of planetary systems around stars like our own, making young solar analogues important targets for both magnetic activity studies and in the search for exoplanets.
Today's Sun undergoes ∼ 22 yr magnetic cycles, wherein the largescale magnetic field polarity reverses twice (every ∼ 11 yr) and between polarity reversals the magnetic field is converted from a dipolar, poloidal structure at magnetic minima, to a slightly more toroidal and complex geometry at magnetic maxima (DeRosa et al. 2012).The emergence of magnetic flux at the solar surface manifests in the growth and decay of dark spots, bright faculae and plages throughout an 11 yr 'spot cycle', and also gives rise to a chromospheric activity cycle that is measurable using the emissions in the cores of the ★ E-mail: emma.brown@usq.edu.auchromospheric Ca ii H and K spectral lines (Wilson 1978).Solar-like magnetic cycles have been identified in three other mature Sun-like stars based on periodic reversals of the large-scale magnetic field polarity that coincide with chromospheric activity cycles; 61 Cyg A (1.3 -6 Gyr, K5V, ∼ 14 yr magnetic cycle, Boro Saikia et al. 2016Saikia et al. , 2018)),  Boo (0.9 Gyr, F8V, ∼ 240 d cycle, Fares et al. 2009;Mengel et al. 2016;Jeffers et al. 2018) and HD 75332 (1.88 Gyr, F7V, ∼ 1.06 yr cycle, Brown et al. 2021).
Young solar-type stars have shown comparatively chaotic magnetic behaviour, with magnetic fields that are an order of magnitude or more stronger than those seen in mature stars, often multiple coexisting activity cycles (Oláh et al. 2016;Boro Saikia et al. 2022;Jeffers et al. 2022) and polarity reversals that do not necessarily correspond to chromospheric activity cycles (Boro Saikia et al. 2022).One of the closest young solar-type stars, the ∼ 0.5 Gyr old (Barnes 2007) K2 dwarf  Eridani, potentially alternates between chaotic and cyclic variability (Metcalfe et al. 2013;Petit et al. 2021), and the results from Jeffers et al. (2022) are consistent with a dynamo of two cycles where the shorter 3 yr cycle is a modulation of the longer 13 yr cycle.Young solar-type stars typically show strongly toroidal and axisymmetric magnetic fields (See et al. 2016;Folsom et al. 2018) compared to the Sun's dominantly poloidal magnetic field, with the most rapidly rotating main-sequence stars often showing strong rings of azimuthal field around the rotational axis (e.g.Marsden et al. 2006;Jeffers & Donati 2008;Jeffers et al. 2011).Tomographic imaging of surface temperature contrasts, known as Doppler Imaging (DI), reveals large polar spots, or collections of spots, on young stars that are sustained for many years, in addition to comparatively short-lived low-to mid-latitude spots (e.g.Marsden et al. 2006;Jeffers & Donati 2008;Jeffers et al. 2011;Lehtinen et al. 2022).This is in stark contrast to the Sun, where spots are restricted to low latitudes, cover only 0.3 to 0.4 percent of the solar surface at any time (Solanki & Unruh 2004) and last only days to months.
The extreme magnetic activity of young stars impacts on the ability to detect exoplanets using the radial velocity (RV) method, which is an important adjunct to other exoplanet detection techniques because it places unique constraints on the planetary mass.Exoplanet detection using the RV method relies on identifying Doppler shifts in stellar spectral lines to measure small changes in the RV of the star that are induced by orbiting objects.Magnetic activity complicates the detection of true Doppler shifts because spots and bright regions cause asymmetries in spectral line profiles that shift the line centre as they rotate in and out of view.These line shifts result in apparent, quasi-periodic RV variations with amplitudes up to hundreds of m s −1 , which can obscure or even mimic the RV signatures of exoplanets (Nava et al. 2020;Barnes et al. 2011;Jeffers et al. 2014).Other magnetic activity signatures, including a net convective blueshift across the stellar disk (Meunier et al. 2010) and unresolved magnetic spots (Lisogorskyi et al. 2020), can each impart RV amplitudes as large as 10 m s −1 , with this alone dwarfing the stellar RV signatures of planets up to 20 M ⊕ .
Young solar analogues make particularly difficult targets in the search for exoplanets, being near-solar mass which makes true exoplanet RV signatures small, and highly active, making activityinduced RV amplitudes large.However, Jeffers et al. (2014) predicts that hot Jupiters (HJs) can be detected for young G and K dwarfs if they are observed over a sufficient time-base of 100 epochs.So far, for stars with masses within 10 per cent of  ⊙ and ages up to 500 Myr, shown in Figure 1, only one has had an exoplanet detected using the RV method alone.
Methods to model and filter activity-induced RVs, such as the popular method of Gaussian Process Regression (GPR) and the 'DI' approach (see Section 3), have mostly been applied to near-solar mass stars with low-to-moderate levels of magnetic activity (e.g.Haywood et al. 2014;Heitzmann et al. 2021).GPR has also been tested against simple simulated spot patterns (Perger et al. 2021), but the reality is that young solar analogues have more complex spot topologies, with both long-lived polar spots and rapidly varying lower-latitude features (e.g.Marsden et al. 2006;Şenavcı et al. 2021;Willamo et al. 2022b).Heitzmann et al. (2021) tested the techniques of DI and GPR on the moderately active, 17-32 Myr, 1.3  ⊙ G2 dwarf HD 141943, using legacy observations from the non-stabilized University College London Echelle Spectograph (UCLES3 ) at the Anglo Australian Telescope (AAT).They found GPR to be the more successful technique, able to recover synthetic planetary signals with semi-amplitudes ≥ 100m s −1 (i.e. a 1 M Jup planet with orbital period of 5 d, for HD 141943).However, the spot pattern in HD 141943 that consists of low-and high-latitude features with similar strengths and area coverage, still contrasts against the intense polar spots that are ubiquitous among the most active young Suns and their significantly weaker low-latitude features.For active, low-mass stars, such Only the star TAP26 (red) has had an exoplanet detected using only stellar RVs.For the remaining stars (black) exoplanets were detected using transits, direct imaging or astrometry.Data taken from The Extrasolar Planets Encyclopaedia2 .as the 0.344  ⊙ , M3.5 dwarf EV Lac (Jeffers et al. 2022), the use of low-resolution DI activity models to remove the spot-induced RV component have been shown to improve the planet-mass sensitivity by a factor of three.The effectiveness of the DI and GPR methods for modelling the activity-induced RVs of highly active solar-mass stars has not yet been evaluated.
The young solar analogue V889 Her (HD 171488, G2V, Montes et al. 2001, 30-50 Myr old, ∼1.06 M ⊙ , Strassmeier et al. 2003) is one of the brightest solar analogues that is rotating fast enough for active regions on its surface to be resolved using the technique of high-resolution DI.This makes it an important target to gain insight into the magnetic dynamos operating inside the most active Sunlike stars, and an ideal candidate to model activity-induced RVs using both the DI and GPR methods.V889 Her has been monitored spectropolarimetrically with the NARVAL4 spectropolarimeter at Télescope Bernard Lyot (TBL, at the Observatoire Pic du Midi, France, Aurière 2003), HARPSpol (High Accuracy Radial velocity Planet Searcher polarimeter) at the 3.6-m ESO telescope (La Silla Observatory, Chile, Snik et al. 2011) and the SEMEL polarimeter (SEMpol) at the AAT (Siding Springs, Australia, Semel et al. 1993).This large data set provides an opportunity to study the long term variability of magnetism in V889 Her, and evaluate the impact of an extended baseline of observations on exoplanet detection capabilities for young solar analogues.
The large-scale magnetic field structure and brightness topology of V889 Her have been reconstructed in a number of previous studies (Strassmeier et al. 2003;Marsden et al. 2006;Jeffers & Donati 2008;Järvinen et al. 2008;Huber et al. 2009;Frasca et al. 2010;Jeffers et al. 2011;Willamo et al. 2022a), which have shown a strongly toroidal field geometry and a large spot (or collection of small, unresolved spots) that covers the visible pole of the star.The polar spot seems to have persisted from at least 1998 to 2011 (Strassmeier et al. 2003;Willamo et al. 2022a).Long term photometric monitoring by Lehtinen et al. (2016) and Willamo et al. (2019b) has shown V889 Her to have a brightness cycle with a period of ∼ 7 − 9 yr, but long-term S-index observations have not been available to detect a chromospheric activity cycle.Willamo et al. (2019b) did not find any consistent anti-correlation between the long-term brightness variations and spot filling factor (percent of the stellar surface covered in spots), as would be expected for a rapidly rotating star (Radick et al. 1990), which they attributed to unreliable filling factor estimates.Willamo et al. (2019b) also found that V889 Her could have experienced a grand maximum in its magnetic activity between 2007 and 2017, based on a long-term minimum in brightness (if taken as a maximum in spot coverage/activity).Additionally, V889 Her has been monitored with NASA/MIT's Transiting Exoplanet Survey Satellite (TESS 5 ).Potential transit signals have been flagged in the TESS data validation report summaries, with periods ranging between ∼ 10 and 382 d, but none of these candidates have been elevated to a TESS Object of Interest (TOI).
In this work we study the magnetic and chromospheric activity variability of V889 Her, and model the activity-induced stellar RV shifts in an attempt to uncover the true stellar RV and search for RV shifts induced by exoplanets.The remainder of this paper is organised as follows; in Section 2 we describe the data used in this study.In Section 3 we outline the techniques used to reconstruct the large-scale magnetic field and brightness topology of V889 Her, extract precise RVs, chromospheric Ca ii H&K activity and model the activity jitter.In Section 4 we present the results of brightness and magnetic field modelling, and analyse the long-term evolution of the large-scale magnetic field and brightness topology of V889 Her.The activity-induced RV variability is discussed in Section 5, as well as the challenges encountered in modelling the activity-induced RVs using GPR and DI to search for exoplanets.Finally, our conclusions are summarized in Section 6.

Spectropolarimetry with NARVAL, HARPSpol and SEMpol
Our data include both previously published and new observations of V889 Her covering the period from 2004 May to 2019 July.We used observations from the NARVAL spectropolarimeter, HARP-Spol polarimeter in conjunction with the HARPS spectrograph, and SEMEL polarimeter coupled with UCLES.HARPS spectra have a resolution of ∼110000 and spectral coverage from 3800 to 6900 Å compared to the data from NARVAL with a resolution of ∼65000 and spectral coverage from 3700 to 11000 Å, and UCLES spectra with a resolution of ∼70000 and wavelength coverage 4400 to 6800 Å.
The observations are summarized in Table 1 and further details are provided in the supplementary materials, with previously published data indicated where applicable.Each of the instruments is configured for spectropolarimetry with the polarimeter mounted at the Cassegrain focus, to split incoming circularly polarized (Stokes V) light into two orthogonally-polarized beams, both of which are fibre-fed to the spectrograph and recorded on the detector.A single Stokes V observation consists of a series of 2 or 4 sub-exposures, where the positions of the orthogonallypolarized beams are switched between sub-exposures.The Stokes V observation is determined by dividing sub-exposures with orthogonal polarization states, which removes spurious polarization signals to a first-order degree.Adding sub-exposures provides the unpolarized, Stokes I observation.For series' of 4 sub-exposures, a 'Null' spectrum can also obtained by dividing spectra with identical polarization states, and this provides a measure of noise and a diagnosis of the reliability of the polarimetric measurement.For further information on this process see Donati et al. (1997).

Data reduction and calibration
The observations were reduced and calibrated using the automated software package libre-ESpRIT (based on ESpRIT, Echelle Spectra Reduction: an Interactive Tool, Donati et al. 1997).libre-ESpRIT reduces the observations using bias and flat-field exposures averaged from series' of exposures taken on the observing night.A first-order wavelength calibration is carried out using a Thorium-Argon arc lamp exposure.For NARVAL and UCLES spectra, the wavelength calibration is then further refined using telluric lines as RV references, as was first done by Donati et al. (2003).Telluric corrections are carried out for each observation individually, and this reduces the relative RV shifts to 100 m s −1 for UCLES spectra (Donati et al. 2003) and 30 m s −1 for NARVAL spectra (Donati et al. 2008;Moutou et al. 2007).For HARPS spectra we adopt an instrumental precision of 5 m s −1 based on Trifonov et al. (2020).

LSD spectral line profiles
We used the technique of Least-squares deconvolution (LSD, Donati et al. 1997) to derive high signal-to-noise spectral line profiles from the observed spectra.LSD assumes that the spectral line asymmetries caused by magnetic activity are roughly equal across all spectral lines.Therefore, the entire observed spectrum of an active star can be thought of as a convolution of a 'mean' spectral line profile (i.e. the LSD profile) and the spectrum of the star when it is magnetically inactive.To model the spectrum of V889 Her we used a line mask for a  eff = 5750 , log  = 4.50 and log (M/H) = 0.00 star from Marsden et al. (2014).The line mask excluded lines with depth less than 10 per cent of the continuum, and strong lines such as the chromospheric Ca ii H&K lines.LSD was performed using the automated routine included in libre-ESpRIT (Donati et al. 1997), and using the LSD normalization parameters from Marsden et al. (2014).Uncertainties for each of the LSD spectral pixels are scaled up from the spectral uncertainties as described in Wade et al. (2000a,b).V-band photometry obtained with the Tennessee State University T3 0.4 m Automatic Photoelectric Telescope (APT) located at Fairborn Observatory, Arizona, and spanning the years 1994 to 2018, was also previously published in Lehtinen et al. (2016) and Willamo et al. (2019b), and provided on VizieR (Willamo et al. 2019a).The observations were made at a cadence of ∼ 1 − 3 d.They were made differentially with respect to the comparison star HD 171286 and the check star HD 170829, which is a constant star that is observed at the same time as the variable target to confirm that the comparison star has constant brightness.For further information on the acquisition and reduction of the differential photometry see Lehtinen et al. (2016) and Willamo et al. (2019b).

Image reconstruction with ZDIpy
The tomographic imaging technique of DI and its derivative, Zeeman Doppler Imaging (ZDI), reconstruct the brightness contrasts and 3dimensional large-scale magnetic field across the stellar surface from series' of Stokes I and V LSD profiles respectively.In this work we used the ZDIpy code by Folsom et al. (2018), which is a python implementation of the extensively used C code described in Donati & Brown (1997) and uses a principle of maximum entropy image reconstruction, wherein the feature content of the reconstructed map is minimized.The resulting spot or magnetic map provides a lower limit for the amount of spot coverage or magnetic field strength, and reliably recovers the fractions of different magnetic field components (Lehmann et al. 2019).
ZDIpy divides the stellar disk into equal-area tiles.For each tile the Stokes I profile is modelled with a psuedo-Voigt profile (convolution of a Gaussian and Lorentzian, approximated according to Humlícek 1982), which is a more realistic approximation of the Stokes I profile compared to a Gaussian (Folsom et al. 2018).The Stokes V profile is derived using the weak-field approximation (Donati et al. 1997).The radial, azimuthal and meridional magnetic field components are each represented as a series of spherical harmonics, with the order of the harmonic expansion limited by a factor   that is related to the complexity of the magnetic field (e.g. = 1 for a dipolar field,  = 2 for quadrupolar field etc.).The full, disk-integrated line models are derived by summing the contributions of all tiles across the stellar surface.The contribution of each tile is scaled by its projected area and relative Doppler shift, taking into account the impact of limb darkening.The brightness or magnetic field (depending on which is being fitted) in each tile is iteratively varied, and the disk-integrated model computed.A model surface brightness or magnetic field map is simultaneously reconstructed from the model line profiles for each iteration, using a principle of maximum entropy image reconstruction, wherein the image feature content is minimized (Donati et al. 2006).The final maximum entropy map is that which achieves a target reduced- 2 fit to the observations.We used a target reduced- 2 value of 0.1 for Stokes I fitting, and 0.9 for Stokes V fitting (which is just below the noise level).It is usual to over-fit the Stokes I line profiles because the signal-to-noise ratio for the LSD Stokes I profiles is underestimated (Petit et al. 2004a;Wade et al. 2000a,b).
ZDIpy is capable of modelling the stellar surface brightness using a 'filling factor' or spot-only approach, whereby the star is modelled as a quiet surface (relative brightness = 1) with dark regions (relative brightness < 1), or the surface can be modelled as a combination of bright (relative brightness > 1) and dark (relative brightness < 1) regions.The spot-only approach has been used in all previous spot studies of V889 Her (Strassmeier et al. 2003;Marsden et al. 2006;Jeffers & Donati 2008;Järvinen et al. 2008;Huber et al. 2009;Frasca et al. 2010;Willamo et al. 2019bWillamo et al. , 2022a) so here we use both methods to compare the resulting spot maps and determine if bright features like faculae are negligible in the DI analysis.Brightness mapping is a somewhat simplified approach compared to temperature mapping; in brightness mapping each surface pixel has a relative brightness, while in temperature mapping each pixel has a spot and a quiescent photosphere temperature, as well as a filling factor for the portion of the pixel covered by spots (e.g.Collier Cameron 2001).In a single line analysis, this spot temperature and filling factor are effectively degenerate, an issue that brightness mapping sidesteps by using a simplified model with only one free parameter per pixel.In brightness mapping spots are assumed to only produce changes in line shape, not equivalent width.Thus the reconstructed brightness maps are essentially only sensitive to variations in the LSD profile shapes, not equivalent widths.
ZDIpy can model the magnetic field with or without taking into account a brightness map.When a brightness map is used as input for the magnetic field inversion, the code adjusts the inversion to account for the suppression of magnetic flux in dark regions.We tested our magnetic inversions when using spot-only and spot + bright maps as input, and found no significant difference in the recovered magnetic field.The magnetic maps presented in our results take the spot-only maps as input.
ZDIpy can account for differential rotation of surface features using a solar-like differential rotation law: where Ω() is the angular rotation rate at a latitude , Ω eq is the angular rotation rate at the equator, and dΩ the difference in rotation rate between the equator and the poles (i.e. the differential rotation rate), all measured in rad d −1 .
Table 1 provides a summary of the spectropolarimetric observations organised into 14 epochs, each with a sufficient number of observations (at least 5) and adequate rotational phase coverage to map the stellar surface using DI and ZDI.The time-spans covered by the data subsets range from 1.64 rotation cycles (∼ 3 d) in 2011 May 14-16 and 2011 May 17-19, through to 8.20 rotational cycles (∼ 11 d) in 2018 July.Even when adopting significant differential rotation in our models, we found that we needed to split the 2011 May data into small sub-sets to achieve a good fit of the ZDI model to the observations, consistent with Willamo et al. (2022a) and the rapid surface evolution seen in the TESS light curve for Sector 40 in Figure A1.This suggests that field evolution occurs rapidly, on a time-scale as short as ∼3 d.For other epochs, such as 2018 July through to 2019 July, we were able to easily fit to the noise level with much longer series' of observations.The difference could be due to the higher resolution of the HARPSpol data, which may better recover the evolution of small scale magnetic features, or possibly evolution of the stellar surface is not always as rapid as in 2011 May.We tested shorter series' of observations for 2018 June through to 2019 July, but did not find that this significantly affected the results we report in Section 4.3.The magnetic and bright/dark features in our reconstructions for these epochs (Figures B1 to B3) do not show any obvious smearing that would suggest more rapid evolution of the magnetic field was occurring.

Stellar and model parameters
We modelled the Stokes I LSD line profiles using a solar line model with the properties shown in Table 2.The line model parameters were fitted to a solar LSD profile, derived from a NARVAL spectrum of sunlight reflected from the moon.The line model is considered to be a suitable approximation for V889 Her which has a similar spectral type to the Sun.
ZDIpy also requires input of the centre-of-mass RV of the star, its line-of-sight projected rotational velocity ( sin ), inclination angle (), equatorial rotation period (P eq ) and rate of differential rotation (dΩ).To select appropriate parameters we used a process of  2 minimization, whereby we iteratively varied input parameters and recomputed the stellar model to identify the parameters that minimize the reduced- 2 fit to the observations.
For the RV,  sin  and inclination angle, we varied the parameters individually and fitted to the Stokes I LSD line profiles.The bestfit RV fluctuated slightly between epochs, as shown in Table B1.
It is not clear if these long-term RV shifts are true centre-of-mass fluctuations, but to optimize the DI and ZDI models we have adopted the best-fit RV on an epoch-to-epoch basis.The best-fit  sini and inclination angle are shown in Table 2, and are reasonably consistent with previous work (Strassmeier et al. 2003;Marsden et al. 2006;Jeffers & Donati 2008;Jeffers et al. 2011;Willamo et al. 2019b).
The inclination angle is typically the most difficult parameter to constrain using  2 minimization.The 60 ± 5 • we derived using  2 minimization is consistent within the error bars with an inclination angle of 64.5 • derived using our measured equatorial rotation period and  sin , adopting a stellar radius of 1.09 ⊙ (Strassmeier et al. 2003) and by assuming solid-body rotation.As discussed by Willamo et al. (2022a), the possible detection of transiting exoplanets by TESS suggests a higher inclination angle if we assume that the planet's orbital axis is aligned with the star's axis of rotation.We tested our magnetic field and surface brightness reconstructions using a greater inclination angle of 85 • and the models converged at much higher reduced- 2 values, indicating a poor fit of the model to the data.
Given that none of the potential transits detected by TESS have been elevated to a TOI, we have adopted the best-fit inclination angle of 60 • for ZDI.
To determine the optimal equatorial rotation period and rate of differential rotation we varied the parameters together according to Equation 1.We fit this using the Stokes V profiles, since V889 Her shows more latitudinal distribution of magnetic features compared to spots, which are mostly observed over the pole both in our study and in previous work (Strassmeier et al. 2003;Marsden et al. 2006;Jeffers & Donati 2008;Jeffers et al. 2011;Willamo et al. 2019bWillamo et al. , 2022a)).The HJD corresponding to rotational cycle zero (i.e. the date where the magnetic field is assumed to be unsheared by differential rotation) is given in Table 1 for each epoch.
In Figure 2 (left panel) we show the reduced- 2 landscape for our 2007 May data set, which has excellent phase coverage for tracing the differential rotation of magnetic features (29 observations across 3.87 rotational cycles).The model can be fitted to the noise level, which suggests that the stellar surface is not significantly evolv-ing during the observations.In Figure 2 the colour represents the reduced- 2 fit (with a fixed entropy) for each combination of P eq and dΩ, and we also show the 1- and 2- uncertainty contours around the best-fitting values.In the right panel of Figure 2, we show the 1- uncertainty regions for several of our epochs that have adequate phase coverage to measure differential rotation.The best-fit rotational parameters derived for each epoch vary, which is consistent with the varying parameters derived across previous studies of V889 Her.Given that the uncertainty regions for these epochs overlap, roughly forming a larger uncertainty ellipse, we take these apparent epoch-to-epoch variations as an indication that the uncertainty of our measurement is larger than suggested by the individual 1- uncertainty regions, similar to Marsden et al. (2007).For all our models we have adopted the median equatorial rotation period of P eq =1.328 d and median rotational shear of dΩ = 0.319 rad d −1 .
Finally, ZDIpy uses a maximum degree of harmonic expansion,   , to control the complexity of the magnetic field reconstruction.We tested increasing   , and found there to be no improvement in the reduced- 2 fit nor any obvious changes in the reconstructed magnetic field geometry for   > 18, therefore we have adopted   = 18.

Uncertainty analysis
We tested the sensitivity of our magnetic field and brightness models to changes in the P eq , dΩ, RV and  sin  within their uncertainties.The variation bars shown in Table B1 and Figure 4 for the DI and ZDI results relate to the changes in the model outputs while we simultaneously varied the input parameters within their uncertainties.The largest impact comes from variations in the  sin  which relates to the rotational broadening of the line profile.Changes in the  sin  are reflected as significant increases in spot coverage.

Chromospheric activity
We measured chromospheric activity for each RV-corrected observation (Section 3.4) using the S-index, which is a measure of the  excess flux in the cores of the Ca ii H and K absorption lines that is related to non-thermal heating of the chromosphere.S-index observations carried out with NARVAL and HARPS were converted to the common Mount Wilson scale (Wilson 1978) using the formula where   and   are the fluxes in two triangular bandpasses centred on the cores of the Ca ii H and K lines (3968.469Å and 3933.663Å) with widths of 2.18 Å (at the base), and    and    are the fluxes in two rectangular 20 Å bandpasses centered on the continuum either side of the H&K lines at 3901.07 Å and 4001.07Å.The calibration coefficients a, b, c, d and e for NARVAL data are available in Table 4 from Marsden et al. (2014), and the coefficients for HARPS data were taken from Gomes da Silva et al. ( 2020).The wavelength coverage of UCLES spectra do not allow for the S-index to be calculated.

Longitudinal magnetic field
The longitudinal magnetic field,   (G), is the line-of-sight magnetic field strength averaged over the visible stellar surface.For each RVcorrected observation we measured   directly from the LSD profiles using the equation where V(v) and I(v) are the Stokes V and I LSD profiles respectively,  (580 nm) is the central wavelength of the LSD profile,  is the mean Landé factor (1.22), c is the speed of light in km s −1 and  is velocity in km s −1 (Donati et al. 1997).We integrated over a velocity domain of -52 to 52 km s −1 with respect to the line centre to include the entire Stokes V polarization signal while minimizing noise.  measurements for each observation are provided as supplementary materials, with error bars determined by propagating the uncertainties computed during the reduction process for each spectral bin of the normalized spectrum through Equation 3.

Precise RVs from LSD profiles
RVs measured from LSD spectral line profiles can be used to detect exoplanets, as demonstrated by Lienhard et al. (2022).We measured precise RVs for each individual observation by computing the firstorder moment (FOM) of the continuum-normalized Stokes I LSD profile, according to where   is the continuum intensity level,  () is the intensity at a velocity  (km s −1 ), and the integral is computed over a velocity range of -76 to 30 km s −1 in the heliocentric corrected frame.RVs for each individual observation are included in Table 1 of the Supplementary Materials.To determine the uncertainties in each RV measurement, we derived the uncertainty in the FOM by propagation of the uncertainties in the LSD line profile, and added this in quadrature to the instrumental uncertainty from Section 2.2.

Photometric rotation period
The TESS light curves in Figure A1 show clear rotational modulation related to brightness contrasts on the stellar surface.The light curve from sector 40 is particularly interesting, with the shape of the modulation changing significantly over ∼ 15 days.The beating and sub-period structure in sector 40 is likely due to rapid evolution of ≥ 2 spots at different latitudes and the presence of significant differential rotation.We measured periods of 1.329 d through to 1.432 d for the modulations using the Lomb-Scargle periodogram.These are in reasonable agreement with our adopted equatorial rotation period and differential rotation rate in Table 2, and would reasonably correspond to the differentially rotating equatorial and polar spots.

Photometry smoothing
The APT differential photometry shows both rotational modulation and significant long-term variability, as has been previously discussed by Lehtinen et al. (2016) and Willamo et al. (2019b).We used a moving average to smooth out short-term variations and observe longer-term periodicities.Willamo et al. (2019b) previously took the mean photometric magnitude from a sliding 30 d window of observations, with a minimum of 14 observations required to be in the window.For the same data presented by Willamo et al. (2019b), we have used a more vigorous smoothing approach to search for brightness fluctuations on a timescale of months to years, taking the mean magnitude from a 180 d sliding window which is similar in length to an observing season.The raw and smoothed data are shown in Figure 5, and we discuss this further in Section 4.3.

Activity filtering and an RV search for exoplanets
The stellar RVs we measured from spectral line profiles are a superposition of true centre-of-mass RVs, and apparent spectral linecentre shifts caused by surface temperature contrasts (i.e.dark spots and bright regions) that change the line shape.We modelled these activity-induced RV fluctuations using two methods, the DI method and Gaussian Process Regression (GPR), and searched for underlying periodic RV fluctuations that could indicate the presence of a planetary companion.

The 'DI method'
The DI method (e.g.Donati et al. 2014;Barnes et al. 2017) makes use of the synthetic Stokes I profiles that were generated to model the surface brightness topology of V889 Her (Section 3.1.1),and which give an estimate of apparent RV shifts related to surface brightness features.We measured RVs from both the raw and synthetic LSD profiles using the method described in Section 3.4, and then subtracted the synthetic RVs directly from the raw RVs.We used the Generalized Lomb-Scargle periodogram (Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009) to search for periodicities in the residual RVs, and assessed the significance of peaks using False Alarm Probabilities (FAPs) calculated using a bootstrap re-sampling approach.Following Zechmeister et al. ( 2009), we generated 10000 re-sampled data sets, for each set retaining the timestamps of the RV observations and assigning randomized RV values (with replacement) to each time stamp.We then used the GLS periodogram to find the highest-power peak for each re-sampled data set, and used the maximum peak heights from all 10000 data sets as an estimate of the probability distribution of power maxima.The power values that are above 99.9 and 99 percent of all power maxima are taken as the 0.1 and 1 percent FAPs respectively.

Gaussian process regression
The second method was to model the activity-induced RV shifts and their temporal evolution as a Gaussian Process (GP).GPR treats the activity-related RVs as correlated Gaussian noise, where the covariance matrix, C, specifies the correlation between each pair of RV observations.We used the GPR implementation of Heitzmann et al.
(2021), which follows Haywood et al. (2014), with each entry in the covariance matrix computed from the quasi-periodic kernel where the hyper-parameter  1 (kms −1 ) is the semi-amplitude of the RV signal related to activity,  2 (d) is the decay parameter, taken to represent the typical lifetime of stellar surface features,  3 (d) is the recurrence time-scale of features (∼ P rot ) and  4 is a smoothing parameter that indicates how complex the signal is within one period.We fixed  4 at 0.31 according to Perger et al. ( 2021), but we found that having a fixed or variable  4 made little difference to our findings in Section 5.3.The uncertainty in each data point,   , and an extra uncorrelated noise parameter,   , are added in quadrature and applied to the diagonal of the covariance matrix (i.e. the variance of the data).We also include an RV offset (RV  ) for each observing instrument.We fitted the GP to the raw RVs using the nested sampling Monte Carlo routine pymultinest (Buchner et al. 2014), a python implementation of multinest (Feroz et al. 2009).Priors were taken from Table 3.In GPR the rotational spot signal is allowed to change its amplitude and phase over time, which is physically motivated by the short-term changes in mid-to low-latitude spots on V889 Her.Given that the large polar spot on V889 Her is very long-lived, we set the priors for the decay parameter to trace the formation and decay of comparatively short-lived lower latitude features.GPR can be used to search for exoplanets by also fitting an activity + planet model to the data, and using Bayesian model comparison to assess the significance of the activity + planet model over the activity-only model.The RV shifts due to an exoplanet in a circular orbit follow the equation where RV planet (t) is the RV shift at a time  (km s −1 ), K is the RV semi-amplitude (km s −1 ),    is the orbital period (d) and Φ the orbital phase [0, 1].The Bayesian evidence for each model was computed with the nested sampling MC algorithm during parameter fitting, by integrating over all model parameter spaces according to Equation 7 where L() is the likelihood and () the prior probability for the set of model parameters , and N is the number of model parameters.
According to Trotta (2008) the activity + exoplanet model can be Table 3. Priors adopted for each of the hyperparameters in equations 5 and 6.U and G represent Uniform and Gaussian distributions.Bracketed values are the minimum and maximum for uniform priors, and the mean and standard deviation for Gaussian priors.RV max is the maximum unsigned RV,  RV is the standard deviation of the RV and T is the time-base of the data set.

Hyper-parameter Prior
1 (km s −1 ) U (0,RV max ) accepted over the activity-only model when the difference in the log of the Bayesian evidence is Δ ln  > 5 .

MAGNETIC FIELD AND BRIGHTNESS TOPOLOGY OF V889 HER
Table B1 provides a summary of results from the surface brightness and magnetic field mapping.The full set of surface brightness and magnetic field maps are provided in Figures B1 through to B3.The fits of the model Stokes I and V line profiles to the observations, and the measured RV, S-index and longitudinal magnetic field strength for each individual observation are also provided as supplementary materials.

Brightness maps
The bright+spot, and spot-only maps are all dominated by the presence of a large, slightly off-centre and long-lived polar spot (or multiple, unresolved small spots) that extends down to 40 − 50 • latitude.This is consistent with all previous surface brightness models of V889 Her (Strassmeier et al. 2003;Marsden et al. 2006;Jeffers & Donati 2008;Jeffers et al. 2011;Willamo et al. 2019bWillamo et al. , 2022a)).The results also confirm the presence of lower-latitude spots in the visible hemisphere of V889 Her.Previous work by Marsden et al. (2006), Jeffers et al. (2011) and Willamo et al. (2019b) also showed smallscale dark features at the mid-to low-latitudes whereas Jeffers & Donati (2008) did not find evidence for such spots in 2005 May to June using data from the MuSiCoS spectrograph (Baudrand & Bohm 1992;Donati et al. 1999) coupled with the TBL, which were unavailable for reanalysis in this work.This difference could be related to the slightly lower resolution of MuSiCoS spectra (∼ 35000), or it could be that the lower-latitude spots are not always present, although all of our maps show at least weak low-latitude features.It is not clear from our data if there is a relationship between the strength of the polar spot and the existence of lower latitude features, but the mean spot coverage does show some correlation with the strength of the magnetic field (Figure 4), which we discuss further in Section 4.3.
Between the bright+spot and spot-only models, the most striking difference is that the bright+spot maps show weak but significant spot coverage in the 'southern' hemisphere of the star.Although there are likely to be spots on both hemispheres, it is not clear if the spots shown in the southern hemisphere are artefacts that compensate for the extensive bright regions in the spot+bright models.The spots in the southern hemisphere are reflected in Table B1 as a significantly larger spot coverage for the bright+spot models compared to the spot-only models.Although we found that the bright+spot models typically converged to lower  2 values compared to the spot-only models, this is to be expected for any model with a greater number of variables, and it is not clear which of the inversion techniques results in the most reliable spot model.
The percentage of the stellar surface covered in spots varies between 3.6 percent in 2008 May and 7.3 percent in 2011 May 14-19, for the spot-only models.Our spot map for 2004 Sept. is similar to that produced by Marsden et al. (2006) using the same data.They obtained a spot coverage of 5.5 percent, which is in reasonable agreement with our 6.9 percent in the spot-only model, considering that Marsden et al. (2006) used a different ZDI implementation from Brown et al. (1991) and Donati & Brown (1997).

Large-scale magnetic field maps
The magnetic field maps in Figures B1 through to B3 indicate a persistent wreath of positive azimuthal field that surrounds the visible pole of V889 Her at all epochs except 2008 May when the magnetic field is at its weakest.Although the weaker field shown in 2008 May may be partly related to the fairly sparse phase coverage of the data, there is similar phase coverage in 2010 Aug., whereas the recovered magnetic field is comparatively strong.Therefore we are confident that the changes in the field strength are not completely related to differences in the phase coverage of data across epochs.
At latitudes below ∼ 60 • the azimuthal field shows mixed polarity regions that have low-to-moderate field strengths (compared to the positive polar wreath).From 2018 June toward the end of our data, the mid-latitudes become dominated by negative regions of azimuthal field, which form a second wreath surrounding the polar wreath.Similar double wreaths persisting over multiple years have been observed for other active stars, like the rapidly rotating K1 subgiant HR 1099 (Petit et al. 2004b), which has a rotation rate approximately 10 times the solar rotation rate (compared to ∼ 20 times the solar rotation rate for V889 Her).Large-scale azimuthal wreaths have also been reproduced in magneto-hydrodynamoic (MHD) simulations of rapidly rotating Suns by Brown et al. (2011).The presence of such largescale azimuthal magnetic structures at the stellar surface has been taken as an indication of a fundamentally different dynamo process operating in active stars, compared to that operating within the Sun (e.g.Donati et al. 2003).
In contrast to the azimuthal field, the radial field shows mixed polarities at high latitudes, consistent with previous observations of V889 Her by Marsden et al. (2006), Jeffers & Donati (2008) and Jeffers et al. (2011).The radial field over the visible pole and at high latitudes is dominantly negative for most epochs, including in 2011 May 14-15 and 2011 May 17-19, which is consistent with Willamo et al. (2022a).However, the radial field over the pole becomes dominantly positive in 2007 Nov., 2010 Aug., 2013 Sept. and 2019 June, as indicated by Figure 3, which shows the fractional magnetic field strength as a function of latitude for each of the radial, azimuthal and meridional field components.The fractional magnetic field was calculated using the formula where  () and () are the fractional and average magnetic field strengths of each component at latitude , and  is the thickness of individual latitude bands.The switch to a dominantly positive radial field over the visible pole does not indicate a full polarity reversal, with the opposite hemisphere of the star remaining dominantly positive throughout all of our observations.Periods of singular large-scale polarity, where the large-scale magnetic field has the same polarity in both stellar hemispheres, have also been reproduced in the simulations of Brown et al. (2011) for rapidly rotating Suns.Interestingly, the dominant field polarity of the meridional component over the pole also changes for 2007 Nov., 2010 Aug. and 2013 Sept. The meridional field polarity even appears to reverse in the hidden hemisphere in 2019 June, but reverts back before 2019 July 10-19.The meridional field is the weakest field component, possibly because ZDI becomes less sensitive to low-latitude meridional fields with increasing stellar inclination (Donati & Brown 1997).
Our magnetic field reconstructions show a stronger magnetic field compared to previously published results from Marsden et al. (2006), who measured a mean field strength of ∼ 30 G in 2004 Sept. compared to our ∼ 87 G.One reason for this is that we have used different modelling approaches.ZDIpy is capable of fitting to the Stokes V profiles while taking into account the brightness image, so it can account for regions of reduced flux and increase the magnetic field accordingly.As shown in Figure C1, magnetic models reconstructed in consideration of the brightness topology show considerably stronger magnetic structures, in particular the azimuthal polar wreath, when compared to a magnetic model that does not take into account the brightness topology at all.Conversely, the field strengths measured here for 2011 May 14-16 and 2011 May 17-19 are an order of magnitude smaller compared to those measured by Willamo et al. (2022a), potentially due to the use of different ZDI software.

Long-term activity and magnetic field evolution
Figure 4 compares key results from the surface brightness and magnetic field mapping to the chromospheric activity, surface-averaged longitudinal magnetic field strength and radial velocity of V889 Her for the period from 2004 to 2020.The variation bars on the magnetic field parameters in Figure 4 were determined by varying the ZDI model input parameters within their uncertainties, and recomputing the brightness and magnetic maps.All of the trends in magnetic field properties that we discuss in this section are more significant compared to the uncertainty bars.We also show the results of our photometry smoothing in Figure 5, where we used a moving 180 d average to smooth out short-term variations in surface brightness.We have drawn green vertical lines in Figure 5 to indicate the epochs where the large-scale magnetic dipole latitude reaches local minima in Figure 4.The synthesis of these data suggests that there are shortterm, quasi-periodic magnetic field fluctuations in V889 Her, the first to be observed for such a young solar analogue.
Between 2007 May and 2011 May 14-16, V889 Her appears to undergo one full oscillation in terms of its magnetic field strength measured from ZDI.The chromospheric S-index and surface-averaged longitudinal magnetic field strength,   , correlate loosely with the field strength measured from ZDI, with both reaching local minima in 2008 May, and local maxima in 2010 Aug.We determined Pearson's correlation coefficients of 0.52 and 0.39 respectively for S-index and   data binned by epoch, with p-values (probability of finding a correlation where there is none) of 0.15 and 0.30.The magnetic maps also indicate that the field is simple during the magnetic and activity minima in 2008 May, with the fraction of magnetic energy stored in high-order modes ( ≥ 3) at its lowest, while the poloidal field fraction and dipolar field strength are at maxima (Table B1).As the field strength and chromospheric activity grow toward maxima around 2010/2011, the dipolar component of the field decreases and higherorder components increase as the field structure becomes complex and dominantly toroidal.Our observations show the negative pole of the large-scale radial dipole dipping toward the equator on two occasions during this period, in 2007 Nov. and 2010 Aug.These epochs are also when the radial field becomes dominantly positive over the visible stellar pole in Figure 3.
The transformation of the magnetic field during this 3-4 yr period has some similarities to the solar cycle, during which the Sun's magnetic field gains strength and is converted from simple and weak to complex and strong.A key aspect missing for V889 Her, at least in our observations, is the cyclic reversal of the large-scale magnetic dipole, which occurs at around magnetic maxima in the Sun and which regenerates the poloidal field but with opposite polarity.In V889 Her, the large-scale magnetic dipole instead approaches the equator around magnetic maxima, and then the poloidal field is regenerated with the same polarity as before.MHD simulations by Brown et al. (2011) have showed similar polarity reversal 'attempts' for a young Sun rotating at 5 times the solar rotation rate (compared to ∼20× the solar rotation rate for V889 Her).The simulations showed that full reversals can occur in the 5Ω-rotating young Sun model, but not every time magnetic energies undergo a full oscillation.There   were several failed polarity reversal attempts for each successful one.Further observations will be required to determine whether or not full reversals of the magnetic dipole of V889 Her are possible.
Although we cannot confirm the presence of any chromospheric activity cycles from our S-index data, at least two potential cycles appear to be present in the smoothed photometry in Figure 5, which we would expect to relate to the formation and decay of spots in a rapidly rotating star such as V889 Her (Radick et al. 1990).There is a long cycle of 7-9 yrs that is clearly evident in the photometry between 1994 and 2007, and is consistent with the long-term cycle previously reported by Strassmeier et al. (2003), Lehtinen et al. (2016) and Willamo et al. (2019b).There are also short-term fluctuations in the photometry which occur as rapidly as every ∼ 1 − 1.5 yr, around half the period of our observed magnetic field and chromospehric activity fluctuation shown in Figure 4. Interestingly, the magnetic field reversal attempts indicated by our magnetic maps each coincide roughly with local minima in the short-term photometric fluctuations.The fact that the polarity reversal attempts seem to follow photometric cycles supports the presence of repetitive magnetic fluctuations, even though our magnetic field maps have only captured a single fluctuation.Jeffers et al. (2022) recently showed for  Eri, which has coexisting ∼ 3 yr and ∼ 13 yr chromospheric activity cycles, that the magnetic flux generation and emergence follow the longer chromospheric activity cycle, but are modulated by the shorter cycle.They found that large-scale polarity reversals occur with the longer cycle only, and if V889 Her behaves in the same way this might explain why we have observed no polarity reversals here for V889 Her in phase with its 3-4 yr fluctuations.The short-term magnetic modulations may be a consequence of a near-surface dynamo.MHD simulations by Brun et al. (2022) with a near-surface dynamo for rapidly rotating stars produced short cycles (∼ 1 yr), although these are marked by magnetic polarity reversals, with possible quasi-biennial oscillations.
It is clear from Figure 4 that the magnetic field of V889 Her can also evolve much more rapidly compared to the 3-4 yr fluctuations, and that the apparent solar-like variability does not strictly hold across all of our observations.In 2019 June, the large-scale dipole of V889 Her crosses over the equator, but still does not fully reverse.The 2019 July 10-19 and 2019 July 20-29 observations seemingly capture the field evolution directly after an attempted reversal, and the magnetic field evolves rapidly during this time.Between 2019 June and 2019 July 10-19, the dipolar axis moves from the equator to 74 • latitude, and the field strength drops as its structure becomes more simple.∼10 d later, in 2019 July 20-29, the magnetic field has a similar complex structure to 2019 June and has an even stronger field strength, whereas the dipole latitude remains relatively high.This rapid, non-solar like evolution of the magnetic field could be similar to the chaotic field evolution observed for  Eri during its short activity cycle.Metcalfe et al. (2013) and Petit et al. (2021) have also noted that  Eridani potentially alternates between chaotic and cyclic variability.We speculate that the short term evolution of the magnetic field in these young stars could potentially reflect the phase of the observations in the longer activity cycle.In the work of Jeffers et al. (2022),  Eri is approaching a possible grand minimum in activity when its magnetic field is observed to vary chaotically.In our observations of V889 Her, the short-term variability resembles the solar cycle as the star approaches a grand maximum in activity around 2017 (Willamo et al. 2019b), inferred from the photometry in Figure 5, and the field evolution is seemingly more chaotic afterwards.
The relationship between spot coverage and the large-scale magnetic field evolution is not entirely clear from our surface brightness and magnetic field maps.Both the spot-only and brightness + spot models suggest that the spot fraction reaches a local minimum in 2008 May, coinciding with minima in the magnetic field strength both from ZDI and the   , as well as a minimum in the chromospheric S-index.After 2008 May the spot evolution indicated by the spot-only and brightness + spot models are not consistent, and we cannot be sure which model is the most reliable.Based on the bright+spot models, a single spot cycle in V889 Her would seem to span two large-scale dipole reversal attempts in 2007 Nov. and 2010 Aug., contrary to the solar case.Given that the photometric brightness varies much more rapidly than this, it is possible that the long spot cycle inferred by our DI results is an effect of infrequent sampling.The presence of both a long-lived polar spot and short-lived low-latitude features in V889 Her also makes deciphering patterns in mean spot coverage and photometric brightness difficult.Willamo et al. (2019b) compared spot occupancy to the same (unsmoothed) photometric data we show here.Their spot occupancy data showed rapid variability (∼2-4 yrs between local maxima) that would seem to correspond well with the magnetic field and photometric fluctuations shown by our results.For example, their spot filling factors peak between 2009 Sept and 2012 Aug., which is loosely consistent with peaks in the magnetic field strength, complexity and chromospheric activity shown by our observations.The peak amplitude of long-term photometric variations in Figure 5 is ∼0.09 mag across the epochs of our DI observations, which corresponds to brightness variations of ∼ 9 per cent.The spot coverage from DI is, as expected, lower than the variability from photometry but is of a similar order.

Differential rotation
We measured a high level of differential rotation from several of our magnetic field maps, as shown in Figure 2.For epochs with results not shown in Figure 2, which have poorer phase coverage and sampling frequency, the  2 minimization tests either didn't converge to a clear minimum, or we measured zero rotational shear.We consider these to be less reliable measurements.Modelling by Brown et al. (2011) suggested that differential rotation would change throughout a magnetic cycle due to dynamo waves, and observations of some cool stars support this (Petit et al. 2004a), but it is not clear from our data if the evolution of the differential rotation parameters for V889 Her follow any trend with respect to the magnetic fluctuations.Like Marsden et al. (2007) did for IM Peg, we take the apparent epochto-epoch variations in the rotational shear of V889 Her, and the fact that the 1- uncertainty regions from different epochs overlap, as an indication that the uncertainty of our measurement is larger than suggested by individual epochs.

RV variations on both short and long time-scales
The RV as measured from the Stokes I LSD observations shows both short-term fluctuations that relate to the rotation of magnetic surface features in and out of view, and also longer-term RV variability, the nature of which is unclear.
In Figure 6 we show the mean and standard deviation of the observed RVs from Figure 4 for each of our DI epochs (black markers).The standard deviation of the short-term RV variations is the greatest in 2013 Sept.This is also when the fractional spottedness at the equator is the greatest, according to Figure 7 which shows the variation in fractional spottedness with stellar latitude for each of our DI models.The fractional spottedness was calculated in the same way as the fractional magnetic field strength, using Equation 8, but with the average field strength at a latitude  (()) replaced by the average spot occupancy.Figure 7 also shows high spot fractions at the mid-to-low latitudes in 2008 May and 2009 May, and the RVs in Figure 6 show large standard deviations for these epochs.These data suggest that the amplitude of rotational RV modulations in V889 Her is most impacted by low-latitude surface features.
The short term RV oscillations are not about a consistent RV, but rather a slightly varying RV across our observational epochs.This is consistent with the long-term RV variations observed over several months to years by Huber et al. (2009).For comparison against the observed RV in Figure 6, we also show in red the line-centre RV that optimized the DI model fit for each of the epochs.The apparent line-centre RV variations in Figure 6 are not sufficiently larger than their uncertainties for us to believe they are true RV fluctuations, nor are the data adequate to carry out a period search, but we have adopted the fluctuating RVs to optimize the DI and ZDI model fits.In 2011 May 17-19 the offset between the observed and the best-fit line-centre RV is the most significant, and this is also when the spot occupancy is the greatest in Figure 4. Apart from this, there is no clear relationship between the two data sets in Figure 6, and it is not clear if the long-term RV shifts can be explained by the magnetic activity evolution of V889 Her, are indicative of a long-period companion, or a combination of the two.The HARPS data all show mean measured RVs that are higher compared to the AAT and NARVAL data, but we do not expect that the long-term RV fluctuations can be explained by an instrumental offset, because we corrected for the heliocentric velocity of the observatory toward the star and used telluric lines to reduce night-to-night shifts in the NARVAL and UCLES RVs, as discussed in Section 2.2.which are the residual values after subtracting the synthetic RVs from the observed RVs.DI filtering results in a 136 m s −1 improvement in the RV RMS (raw RVs have RMS ∼ 250 m s −1 compared to ∼114 m s −1 for filtered RVs).Power spectra for the observed and residual RVs are also shown in Figure 8, where significant frequencies in the data may correspond to a higher periodogram power.The power spectrum for the observed RVs is dominated by forests of peaks around the 1 d observing cadence and its aliases (0.5 d, 0.25 d etc.).There are also forests of strong peaks around the rotation period (marked with a red line), an alias of the rotational signal at a frequency of 1 − 1 P rot d −1 (marked with a green line) which is related to the observing cadence, as well as aliases of both of these signals that recur at frequency intervals of ∼1 d −1 (dashed red and green lines).The power spectrum for the filtered RVs shows residual peaks at periods of 1 d and 0.5 d, related to the observing cadence.The results are similar for all other data sets we tested.The results for other selected data sets are provided as supplementary materials.

DI activity modelling
Although the power spectra in Figure 8 prove we were able to remove the activity-related rotational RV signal and its aliases by subtracting the synthetic RVs, we were still unable to recover any significant periodicities in the residuals outside those related to the observing cadence, and this is an impact of the limited temporal coverage of our data.For V889 Her, longer periods of dense observations could improve exoplanet detection capabilities, but care will need to be taken when using the DI activity filtering technique because of the underlying long-term RV variability.We are not sure of the nature of the long-term RV shifts, but the varying line-centre RVs adopted for each DI epoch undoubtedly bias the DI activity filtering and thus the exoplanet search.The impact would likely be the greatest when searching for long-period companions, such as planets orbiting within the habitable zone.

GPR activity modelling
We also tested GPR to model the magnetic activity of V889 Her, using a quasi-periodic GP kernel which has been shown to be highly effective at reproducing rotationally modulated activity-induced RVs (e.g.Haywood et al. 2014;Angus et al. 2018;Heitzmann et al. 2021;Perger et al. 2021).However, our GPR models did not converge to a reasonable solution for any of our activity-only or activity + exoplanet tests, so here we explore the challenges of applying GPR to V889 Her as a test case for young solar analogues.

Issue 1: Different decay time-scales for polar and lower-latitude spots
A key issue was constraining the hyper-parameter  2 which relates to the time-scale for active region growth and decay.When we applied GPR to single epochs of data  2 did not converge to a realistic solution (Figure D1).Potentially the RV signal related to the longlived polar spot obfuscates the signal from lower-latitude features, which we would expect to decay rapidly over ∼ 3 to 15 d based on our brightness and magnetic maps, and the TESS light curve in Figure A1.The polar-spot evolves on a timescale much longer than a single epoch of data and Nicholson & Aigrain (2022) found that  2 breaks down when this occurs.When combining data from multiple epochs the GP could fit to either short (several days) or long (several thousand days) spot decay timescales depending on the prior used for  2 (see Figure D2).There was not a significant difference in the log of the Bayesian evidence between these models, so the preferred model is unclear.

Issue 2: Strong differential rotation
When fitting to multiple epochs of data the posterior distribution of the rotation period,  3 , seemed to be correlated with the decay timescale.A slightly longer rotation period was associated with a longer spot decay timescale (likely related to the polar spot), and vice versa.When we fitted the GP to the entire ∼ 15 yr of data,  3 even showed a bimodal distribution (Figure D3).This is probably related to the significant differential rotation exhibited by V889 Her, which means that the RV modulations caused by spots at different latitudes will have different periods.

Issue 3: Insufficient temporal coverage of data
Given the above 2 issues, the temporal coverage of our legacy data is insufficient to model the activity-induced RVs for V889 Her using GPR.Our data is not overall very different from that used by Heitzmann et al. ( 2021 Prot .Aliases of both of these signals repeat at frequency intervals of 1 d −1 , and are indicated by the dashed red and green lines.Inset, we also include the power spectrum of the observed RVs in the frequency domain, where the aliasing of the rotational signal is more clear.The horizontal dashed grey lines indicate the periodogram power levels that correspond to a 0.1 percent and 1 percent FAP. of V889 Her, at 2.198±0.002d, and lower differential rotation rate at Ω = 0.1331 +0.0095 −0.0094 rad d −1 compared to 0.319 rad d −1 for V889 Her (although they measured DR from the Stokes I profiles while we used Stokes V profiles).Their reconstructed DI images indicate a polar spot that appears to be somewhat smaller compared to the polar spot on V889 Her, and they recovered low-latitude features with a similar strength to the polar spot and much larger area coverage compared to in V889 Her, for which the low-latitude features are weak and small relative to the polar spot.The strength of the lower latitude features relative to the polar spot, and the low differential rotation rate in HD 141943, may explain why the GP converges easily for the moderately active star in Heitzmann et al. (2021).
Another issue in modelling magnetic activity for V889 Her is the underlying long-term RV variations and the fact that our observations with different instruments have no temporal overlap.The mean RV from HARPS is ∼ 300 m s −1 higher compared to the AAT and NAR-VAL RVs, differences that the GP can easily account for with large instrumental offsets, but this dilutes the contributions of activity and any exoplanets.Given that we corrected for the heliocentric velocity of the observatory toward the star and used telluric lines to reduce night-to-night shifts in the NARVAL and UCLES RVs, we do not expect such large RV differences between the data sets to be entirely instrumental.
It is clear from our results that DI is more suitable compared to GPR for modelling the magnetic activity of V889 Her using our current data set.DI has more flexibility to model complex magnetic activity, compared to our quasi-periodic GP model.However, DI cannot resolve small scale magnetic structures, which Jeffers et al. (2014) and Lisogorskyi et al. (2020) showed can prevent the detection of up to 20  ⊕ planets in temperate zones around Sun-like stars.Therefore it is important that we improve data-driven techniques such as GPR to model the magnetic activity of young Suns.

Recommendations for modelling extreme magnetic activity with GPR
Although Perger et al. (2021) found that a quasi-periodic GP can fit well individually to solar-like low-latitude spot configurations, or polar spots, our results suggest that the quasi-periodic GP may not be able to constrain the activity evolution for a combination of the two when significant differential rotation is present, at least when using legacy data with large gaps.Future observing plans should focus on obtaining long time-series of continuous high-cadence observations, ideally with nightly observations spanning a few months as is carried out with the RedDots exoplanet search program (Jeffers et al. 2020).Exploration of alternative GP covariance kernels and simultaneous modelling of RVs and activity indicators are also recommended.We performed some preliminary tests combining a quasi-periodic and stable periodic GP component (Kossakowski et al. 2022) to model contributions of the low-latitude and polar features respectively to the activity-induced RVs, but the GP did not converge for our tests.

Exoplanet detection limits with DI
We used V889 Her to estimate the detection limits for HJs orbiting extremely active stars when using existing legacy data and filtering activity jitter using DI.We injected synthetic planetary RV shifts into our data, re-computed our DI models and then searched the DIsubtracted RVs for the injected planetary signals.The planet-induced RV shifts were calculated using equation 6, assuming a circular orbit, and the shifts were applied to the observed LSD line profiles.We limited our tests to planets with orbital periods of 3 to 7 d (orbiting at ∼0.04 to 0.07 AU), and masses of 1 to 2    .
Given our ∼ 15 yr data set we estimate that we could, in principle, detect Jupiter-mass planets with orbital periods of 3 d. Figure 9 (top) shows the power spectrum for the entire series of RVs with the injected 3 d, 1   planetary signal, and the power spectrum of the DI-subtracted RVs.The DI-subtracted RVs show a 3 d peak which is significant with respect to the 0.1 percent FAP level.
At orbital periods of 4-5 d the planetary signal is buried under a forest of peaks with periods around 1/(1 +    ), which are aliases of the stellar rotational signal.The rotational signal and its aliases are removed by DI such that 2    planets at these orbital periods would be detectable with a FAP < 1 percent, as shown in the second panel of Figure 9. 2    planets with orbital periods of 6-7 d are detectable in the DI-subtracted RVs with a FAP < 1 percent.However, for these longer orbital periods there are several strong signals shown in the DI-subtracted RVs and selecting the correct peak is non-trivial.This issue arises as the orbital period becomes similar to the length of most our individual epochs of observations.

CONCLUSIONS
In this study we used previously published and new spectropolarimetric data spanning the years 2004 to 2019, with archival photometry taken between 1994 and 2018, to investigate the variability of the magnetic field, surface activity and RV of the young solar analogue V889 Her.We also tested the techniques of Doppler Imaging (DI) and Gaussian Process regression (GPR) to filter activity-induced RV variations and search the residual RVs for exoplanets.
Our magnetic maps indicate possible quasi-periodic magnetic field fluctuations with a potential period of 3-4 yr.The magnetic field evolution during this period loosely resembles the solar magnetic cycle; the field is simple and dominantly poloidal during magnetic and chromospheric activity minima, and is complex and dominantly toroidal during a magnetic and chromospheric activity maxima.Smoothed photometric data show a 7-9 yr brightness cycle and short-term brightness oscillations with a period of 1-1.5 yr, roughly half the length of the magnetic field oscillation we observed.Contrary to the solar case, no large-scale polarity reversals were observed for V889 Her, but several attempts at polarity reversals were observed.Long-term monitoring will be required to determine if the star undergoes large-scale magnetic polarity reversals and if it truly switches between solar-like and non-solar like behaviour.
Our surface brightness maps confirmed the presence of both a large, long-lived polar spot and multiple smaller short-lived spots at low-to mid-latitudes.The observed stellar RV showed significant rotational modulation, its amplitude tied to the spot coverage at low latitudes.The RVs also showed underlying long-term fluctuations of unconfirmed origin.
We were unable to reliably fit a quasi-periodic GP model to the stellar activity-induced RV variations due to three main challenges: (i) There are competing RV signals from the long-lived and temporally varying polar spot and rapidly evolving low-latitude features; (ii) There is significant differential rotation; and (iii) The temporal coverage of our legacy data is insufficient.
To help overcome these challenges we recommend that alternative GP kernels be tested and high cadence observations of V889 Her (or similar stars) be carried out for as long as possible, ideally with nightly observations spanning a few months as is carried out with the RedDots exoplanet search program (Jeffers et al. 2020).When using existing legacy data, the DI technique can be a useful tool for modelling the activity-induced RV variability of active young Suns, with the flexibility to model complex surface spot distributions.Using the DI method we were able to remove the activityrelated rotational RV signal and its aliases.Given our ∼ 15 yr set of RV observations we estimate that we could, in principle, detect the presence of a 1    planet with an orbital period of ∼ 3 d and, with less certainty, 2    planets with orbital periods between 3 and 7 d.
Midi-Pyrénées and Université Toulouse III -Paul Sabatier.We acknowledge use of the simbad and VizieR data bases operated at CDS, Strasbourg, France.This work has also made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna.ELB is supported by an Australian Postgraduate Award Scholarship.SVJ acknowledges the support of the German Science Foundation (DFG) priority program SPP 1992 'Exploring the Diversity of Extrasolar Planets' (JE 701/5-1).In addition to the python packages referenced in Section 3, this research has made use of the following packages: astropy (Astropy Collaboration et al. 2013), corner (Foreman-Mackey 2016), loguniform (MIT licence; Joao Faria), matplotlib (Hunter 2007), NumPy (Harris et al. 2020), and SciPy (Virtanen et al. 2020).

APPENDIX A: TESS PHOTOMETRY
Figure A1 shows TESS photometry from sectors 26, 40 and 53, where V889 Her was observed in 2-minute integrations.Power spectra derived using the GLS periodogram (Zechmeister & Kürster 2009) are also shown for each data set.

APPENDIX C: IMPACT OF BRIGHTNESS TOPOLOGY AS AN INPUT TO MAGNETIC FIELD MODELLING
Figure C1 shows the impacts on the final reconstructed magnetic field of using pre-computed DI maps for the magnetic field inversion.

APPENDIX D: RESULTS OF RV ACTIVITY FILTERING -GPR METHOD
Figures D1 to D3 show corner plots indicating the posterior distributions and best-fitting parameters for a quasi-periodic GP fitted to various sets of RV observations.Figure D1 (left) shows the corner plot for 2007 May when using a logarithmic prior for  2 .The logarithmic prior allows for greater exploration of short length scales.In practice this resulted in the decay parameter converging toward the lower edge of the prior (i.e.toward zero, see Figure D1) such that the final numerical model did not represent a realistic physical solution.With a uniform prior (right),  2 did not converge to a unique solution, regardless of the upper limit of the prior.Without any strong evidence for the mean or standard deviation of the spot decay timescale, we avoided using an informative Gaussian prior for  2 , which would bias the model and artificially boost the posterior likelihood.
Figure D2 shows the posterior distributions of the model parameters for the combined HARPS data using two different priors for  2 .On the left the GP converges toward a short decay timescale of ∼ 9 d, which may be related to the lengths of individual observing runs.On the right the value for  2 is larger than the total length of the data (∼ 1130 d), and is possibly tied to the long-lived polar spot which is present for the entire data set.B1 when the spot-only (filling factor) model is considered during magnetic field reconstruction.Bottom: The magnetic field when no brightness model is considered in the reconstruction of the magnetic model.The azimuthal and meridional fields are clearly weaker when the presence of the dark polar spot is not considered in the magnetic field reconstruction.There is minimal difference in the radial magnetic field reconstruction across the three models.There is also minimal difference in the three reconstructed field components when using the brightness + spot or the spot-only models as input for the field reconstruction.

Figure 1 .
Figure 1.Age versus mass for planet-hosting stars with masses within 10 percent of M ⊙ and ages below ∼ 500 Myr.Only the star TAP26 (red) has had an exoplanet detected using only stellar RVs.For the remaining stars (black) exoplanets were detected using transits, direct imaging or astrometry.Data taken from The Extrasolar Planets Encyclopaedia 2 .

V889
Her was observed by TESS in 2-minute short-cadence integrations in sectors 26, 40 and 53.The Pre-search Data Conditioned Standard Aperture Photometry (PDCSAP) light curves are shown in Figure A1, which show relative flux variations of up to ∼ 0.15 × 10 5 e s −1 which are assumed to be astrophysical.

Figure 2 .
Figure 2. Left: Example plot of the reduced- 2 landscape when fitting to the Stokes V profiles from 2007 May.The darkest red region represents the best-fit combination of the equatorial rotation period (d) and rotational shear (rad d −1 ).The 1 and 2 contours are also shown.Right: The 1 contours from multiple epochs are plotted on top of each other.The mean and uncertainty of the adopted equatorial period and rotational shear are also shown with the crosshair.

Figure 3 .
Figure 3. Fractional magnetic field strength per latitude band for the radial (top), azimuthal (middle) and meridional (bottom) field components.For epochs addressed in Section 4.2 the data are coloured, and the rest of the epochs are shown as dashed grey lines.The horizontal dashed line at -60 • latitude indicates the lower extents of our magnetic maps.

Figure 4 .
Figure 4. Time-series plots of the surface brightness properties, large-scale magnetic field parameters and magnetic activity metrics for V889 Her between 2004 and 2020.From top to bottom the axes show (a) mean spot coverage for spot-only models, (b) mean spot coverage for bright+spot models, (c) mean magnetic field strength from ZDI, (d) fraction of the poloidal field stored in octopolar and higher-order modes, (e) latitudinal location of the large-scale magnetic dipole (negative pole), with dashed vertical lines indicating possible reversal attempts (f) chromospheric S-index, (g) surface-averaged longitudinal magnetic field strength, and (h) the radial velocity measured from the first-order moment of the observed LSD profile.The shaded regions for plots (a) through (e) indicate the sensitivity of our ZDI results to variations in the P eq ( +0.012 −0.011 d), dΩ (±0.114 rad d −1 ),  sin  (±0.5 km s −1 ) and RV (see Table B1).In axes (f) to (h) solid lines join the mean values for each observational epoch.The grey and blue vertical shaded bands indicate data from UCLES and HARPS respectively.Un-shaded data is from NARVAL.Additionally, in (c), (d), (e) and (h) the data shown with open circles were provided by the authors of Jeffers & Donati (2008) (note that these data have no 'variation bars').In (h) the grey filled circles are RVs from Frasca et al. (2010).
Figure 4. Time-series plots of the surface brightness properties, large-scale magnetic field parameters and magnetic activity metrics for V889 Her between 2004 and 2020.From top to bottom the axes show (a) mean spot coverage for spot-only models, (b) mean spot coverage for bright+spot models, (c) mean magnetic field strength from ZDI, (d) fraction of the poloidal field stored in octopolar and higher-order modes, (e) latitudinal location of the large-scale magnetic dipole (negative pole), with dashed vertical lines indicating possible reversal attempts (f) chromospheric S-index, (g) surface-averaged longitudinal magnetic field strength, and (h) the radial velocity measured from the first-order moment of the observed LSD profile.The shaded regions for plots (a) through (e) indicate the sensitivity of our ZDI results to variations in the P eq ( +0.012 −0.011 d), dΩ (±0.114 rad d −1 ),  sin  (±0.5 km s −1 ) and RV (see Table B1).In axes (f) to (h) solid lines join the mean values for each observational epoch.The grey and blue vertical shaded bands indicate data from UCLES and HARPS respectively.Un-shaded data is from NARVAL.Additionally, in (c), (d), (e) and (h) the data shown with open circles were provided by the authors of Jeffers & Donati (2008) (note that these data have no 'variation bars').In (h) the grey filled circles are RVs from Frasca et al. (2010).

Figure 5 .
Figure 5. V889 Her brightness variations as shown by differential V-band photometric magnitudes taken between 1994 and 2018.The photometric data, originally published by Lehtinen et al. (2016) and Willamo et al. (2019b), are shown as light grey data points, and our moving 180 d averages are shown as black points.Green dashed vertical lines indicate the epochs of possible 'attempted reversals' in Figure 4, and red vertical lines indicate all other ZDI epochs.

Figure 8 (
Figure 8 (top panel) shows the observed RVs and our synthetic, activity-induced RVs from DI, for a combination of data from all DI epochs.The 'filtered' RVs are shown in the bottom panel of Figure 8,

Figure 6 .
Figure 6.Comparison of the observed mean RV (black) for each DI epoch to the line-centre RV adopted to optimize the DI model fit (red, taken from Table B1).Squares, circles and triangles indicate data from UCLES, NARVAL and HARPS respectively.The bars for the measured RVs show the standard deviation of the individual observations, and for the modelled line-centre RVs the bars represent the uncertainties in the best-fit line-centre RVs.The data in red have been shifted to the right for clarity.

Figure 7 .
Figure 7. Fractional spot coverage versus stellar latitude for all of our DI maps.For epochs addressed in Section 5.1 the data are coloured, and the rest of the epochs are shown as dashed grey lines.

Figure 8 .
Figure 8. Left: Observed (black) and synthetic (red) RVs of V889 Her for combined data from all DI epochs (top) and residual RVs calculated by subtracting the synthetic from the observed RVs (bottom).Right: Power spectra for the observed (top) and residual (bottom) RVs, computed using the GLS periodogram.The solid red vertical line indicates the stellar rotational signal, and the solid green line indicates its alias at a frequency of 1 − 1Prot .Aliases of both of these signals repeat at frequency intervals of 1 d −1 , and are indicated by the dashed red and green lines.Inset, we also include the power spectrum of the observed RVs in the frequency domain, where the aliasing of the rotational signal is more clear.The horizontal dashed grey lines indicate the periodogram power levels that correspond to a 0.1 percent and 1 percent FAP.
B3 show the brightness and magnetic field maps for each of our 14 observational epochs.Model fits corresponding to each set of maps are shown in Figures 1 to 4 of the supplementary materials.

Figure B1 .
Figure B1.Left to right: reconstructed brightness + spot, spot-only, radial field, azimuthal field and meridional field maps for V889 Her, from 2004 Sept. through to 2010 Aug.The maps are flattened polar projections, showing the visible pole at the centre and with latitudes extending down to -60 • .Radial ticks indicate the phases of observations.

Figure C1 .
Figure C1.Top: The magnetic field in 2007 May when the brightness + spot model is considered during magnetic field reconstruction.Middle: The magnetic field as shown in FigureB1when the spot-only (filling factor) model is considered during magnetic field reconstruction.Bottom: The magnetic field when no brightness model is considered in the reconstruction of the magnetic model.The azimuthal and meridional fields are clearly weaker when the presence of the dark polar spot is not considered in the magnetic field reconstruction.There is minimal difference in the radial magnetic field reconstruction across the three models.There is also minimal difference in the three reconstructed field components when using the brightness + spot or the spot-only models as input for the field reconstruction.

Table 1 .
Jeffers & Donati (2008)imetric observations of V889 Her used for DI and ZDI.For each epoch, columns list the number of observations, the number of rotational cycles covered by the observations, the HJD of the first (HJD start) and last (HJD end) Stokes V observations, and the HJD corresponding to rotational cycle 0 for the epoch (close to the mid-point of the data, but an integer multiple of rotations from HJD=2453200).The final columns show the instrument that the observations were sourced with and indicate whether the data has been previously published.Full details for individual observations, including observations not used for tomographic imaging, are provided as supplementary materials.*The2005Maydata originally published inJeffers & Donati (2008)were not reanalysed here, but have been included in Table1for completeness.