© The Authors 2025

1 Introduction

Herbig Ae/Be stars are pre-main-sequence objects with pronounced emission line features and an infrared excess indicative of dust in their circumstellar disks. There is observational evidence that Herbig Ae and late-type Herbig Be stars are intermediate-mass analogues of T Tauri stars, but with convectively stable envelopes that do not support the dynamo action found in the fully convective T Tauri stars (e.g. Hubrig & Schöller 2021). A number of magnetic studies have been attempted in the last few years, indicating that about two dozen Herbig Ae/Be stars very probably possess globally organised magnetic fields reminiscent of fields of the classical Ap/Bp stars, of which the magnetic Herbig stars may be the precursors (e.g. Hubrig et al. 2004; Alecian et al. 2013b; Hubrig et al. 2015; Järvinen et al. 2015, 2018, 2019).

Magnetic fields in these stars might be fossils of the early star formation epoch, in which the magnetic field of the parental magnetised core was compressed into the innermost regions of the accretion disks (e.g. Banerjee & Pudritz 2006). Alternatively, it is plausible that the weak magnetic fields detected in a number of Herbig Ae/Be stars are just leftovers of the fields generated by pre-main-sequence dynamos during the convective phase (e.g. Hubrig & Schöller 2021, and references therein). If this scenario is valid, we should expect a significantly larger number of Herbig stars possessing weak magnetic fields. Indeed, fields greater than 200 G are very rare, with most stars possessing fields of about 100 G and less (Hubrig et al. 2015). The studies by Järvinen et al. (2018, 2019) also support this scenario: the authors report the presence of relatively weak longitudinal magnetic field strengths in a few Herbig Ae/Be stars in the range from 120 G to 210 G.

To understand the origin of magnetic fields in Herbig Ae/Be stars, knowledge of the field geometry is indispensable, but requires information on the stellar rotation period. Previous studies of upper main-sequence stars with large-scale organised magnetic fields demonstrated that their magnetic fields are dominated by dipolar fields tilted with respect to the rotation axis (following the oblique dipole rotator model; Stibbs 1950). Thus, their rotation periods were frequently determined by studying the periodicity in the available magnetic data.

Currently, the magnetic periods and corresponding magnetic phase curves presenting the dependence of mean longitudinal magnetic field strength on the rotation phase have only been reported for three Herbig Ae/Be stars, HD 101412 (Hubrig et al. 2011), V380 Ori (Alecian et al. 2009), and HD 190073 (=V1295 Aql; Alecian et al. 2013a). However, V380 Ori has not been confirmed to be a Herbig star and is probably a very young main-sequence star (Hubrig et al. 2019; Shultz et al. 2021). A rotation period of 42.076 d for HD 101412 was reported by Hubrig et al. (2011), who used a combination of photometric observations and measurements of the longitudinal magnetic field based on spectropolarimetric observations with the low-resolution FOcal Reducer low dispersion Spectrograph (FORS 1; Appenzeller et al. 1998) at the Very Large Telescope (VLT). This star has also been reported to possess a strong surface magnetic field of 3.5 kG (Hubrig et al. 2010). The measured mean longitudinal magnetic field exhibits a single-wave variation during the stellar rotation cycle. Such a variability is usually considered as evidence for a dominant dipolar contribution to the magnetic field geometry. Spectroscopic signatures of magnetospheric accretion in HD 101412 have been discussed by Schöller et al. (2016).

The presence of a magnetic field on the surface of the Herbig Ae/Be star HD 190073 has been known of for several years. The first measurement of a longitudinal magnetic field in HD 190073 was published by Hubrig et al. (2006), who indicate the presence of a longitudinal magnetic field ⟨B_z⟩ = 84 ± 30 G at a 2.8σ level, measured with FORS 1. Subsequently, this star was studied by Catala et al. (2007). The authors used observations with the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS; Donati et al. 2006a) installed at the Canada-France-Hawaii Telescope, at a spectral resolution of R ≈ 65 000 and, employing the least-squares deconvolution (LSD) introduced by Donati et al. (1997), confirmed the presence of a weak longitudinal magnetic field, ⟨B_z⟩ = 74 ± 10 G, at a higher significance level. The applied LSD method allows us to construct average photospheric profiles of both the Stokes I and Stokes V parameters, and the diagnostic null N profile, by deconvolving the observed spectra using a line mask with lines identified in these spectra.

Additional observations of HD 190073 using the Narval spectropolarimeter (Aurière 2003) at the 2 m Telescope Bernard Lyot at Pic du Midi and the High Accuracy Radial velocity Planet Searcher polarimeter (HARPSpol; Snik et al. 2008) at the 3.6 m European Southern Observatory (ESO) telescope on La Silla have subsequently been obtained by Alecian et al. (2013a) and Järvinen et al. (2015, 2019) between 2011 and 2019. Using the numerous observations of this star during 2011 May–November and 2012 July–October, and employing the LSD technique, Alecian et al. (2013a) revealed variations of the magnetic field strength on timescales of days to weeks. Furthermore, in contrast to the 2011 measurements, which showed a longitudinal magnetic field of positive polarity for all but one measurement, a sizeable number of measurements in 2012 indicated a field of negative polarity. To explain the observed annual change in the magnetic field configuration, the authors suggested that a phenomenon such as an interaction between the fossil field and the ignition of a dynamo field generated in the newly born convective core occurred, perturbing the magnetic field. In the study of Alecian et al. (2013a), the best fit to the 2011–2012 data set corresponding to a rotation period of 39.8 ± 0.5 d reproduced all but one of the 2012 observations, but a number of the 2011 data points were discordant with this period. No significant period could be identified in the 2011 data set, while a period of 40 ± 5 d was clearly detected in the 2012 data set when analysing these data sets separately.

Since HD 190073 is rather bright, with m_V = 7.7, it has been intensively studied in recent years in various wavelength domains. More accurate fundamental parameters have recently been reported by Guzmán-Díaz et al. (2021) from the analysis of the spectral energy distribution and Gaia Early Data Release 3 parallax and photometry (Riello et al. 2021). The following stellar parameters have been derived: spectral type B9, effective temperature T_eff = 9750 ± 125 K, mass M = 6.0 ± 0.2 M_⊙, radius R = 9.68 ± 0.44 R_⊙, and age 0.30 ± 0.02 Myr. In their study, Aarnio et al. (2017) used a radial velocity v_r = −1.2 ± 1.3 km s⁻¹ and v sin i = 3.19 ± 2.45 km s⁻¹.

The presence of a super-Jupiter-mass candidate around HD 190073 at a separation of about 1.1″ has been announced by Rich et al. (2022) using the Gemini Planet Imager (Macintosh et al. 2014). More recently, Ibrahim et al. (2023) studied HD 190073 using the Center for High Angular Resolution Astronomy (CHARA) and the Very Large Telescope Interferometer (VLTI) arrays. The authors report that their modelling is consistent with a near face-on disk with an inclination ≤ 20° and an average radius of a ringlike structure around the star of 1.4 ± 0.2 mas (1.14 au), which is interpreted as the dust destruction front. The observations indicate the disk around HD 190073 to have skewness on the finer spatial scale, possibly due to a low-mass companion. The sub-au structure that is detected in the image seems to move between the two epochs inconsistently, with Keplerian motion, pointing at dynamics effects from the outer disk. The authors speculate that the motion could be caused by interactions in the outer disk with an object embedded in the inner disk.

In addition, X-ray emission of log L_x = 29.86 ± 0.34 erg s⁻¹ has recently been detected by Anilkumar et al. (2024) from observations with the Chandra X-ray Observatory (Weisskopf et al. 2000). The most popular mechanisms proposed to explain X-ray emission from Herbig Ae/Be stars involve contribution from a low-mass companion in the close proximity of the Herbig Ae/Be stars (e.g. Stelzer et al. 2009), wind shocks, or magnetically confined winds (e.g. Günther & Schmitt 2009). These very recent findings make a reinvestigation of the magnetic properties of HD 190073 especially invaluable in view of a possible magnetospheric interaction between the magnetic host star, its low-mass companion, and the accretion disk. As for the giant planetary companion, additional observations are needed to confirm that it is co-moving with the host star.

After obtaining additional HARPSpol observations, we decided to reexamine the magnetic field configuration of HD 190073. This included a careful inspection of the available spectropolarimetric data and the determination of the rotation period based on our own magnetic field measurements. As the investigation of the rotation modulation of the hydrogen emission lines using dynamical spectra phased with a redetermined rotation period suggests the presence of a compact magnetosphere, we characterised this magnetosphere using the Nirvana MHD code developed in the Leibniz-Institut für Astrophysik Potsdam (Ziegler 2004, 2005). Further, we studied the rotational modulation of a few metal lines and of the He I 10 830 Å triplet, usually used as a diagnostic of magnetically driven accretion and outflows.

This paper is structured as follows: details about the observations of HD 190073 considered in this study are described in Sect. 2. In Sect. 3, we present our magnetic field measurements that we used to determine the rotation period of HD 190073. In Sect. 4, we study the line profile variability of several spectral lines over the rotation period, including the variability of hydrogen emission lines, which clearly shows the presence of a magnetospheric emitting region. The parameters of the detected magnetosphere are explored in Sect. 5 and the results obtained are discussed in Sect. 6.

2 Observational material

2.1 High-resolution HARPSpol observations

Seven high-resolution (R ≈ 115 000) polarimetric spectra of HD 190073 have been recorded with HARPSpol between 2011 and 2017 (Järvinen et al. 2015, 2019). Two more observations were acquired with this instrument in 2019, on June 16 and 18. The latter has been omitted from our analysis due to a very low signal-to-noise ratio. Polarimetric HARPSpol spectra cover wavelengths from 3780 to 6910 Å, with a small gap between 5259 and 5337 Å. The observations obtained were reduced using a dedicated HARPS data reduction pipeline available at the ESO La Silla observatory. Details on the normalisation of the spectra to the continuum level can be found in Hubrig et al. (2013).

2.2 Archival observations

The reduced Narval observations described in the work of Alecian et al. (2013a) were downloaded from the PolarBase (Petit et al. 2014) archive¹. Polarimetric spectra obtained with this instrument have a spectral resolution R ≈ 65 000 and cover wavelengths from 3750 to 10 500 Å. The data set of the available spectral material for 2011 contains 40 polarimetric spectra, whereas for 2012 the data set includes 25 spectra. The numbers of observations in 2011 and 2012 reported in the work of Alecian et al. (2013a) were 31 and 20, respectively. We also downloaded one available ESPaDOnS spectrum obtained in 2012. ESPaDOnS can be considered as a twin brother of Narval with the same wavelength coverage and a similar spectral resolution.

Furthermore, to investigate the rotational variability of the He I 10 830 Å triplet with the transition $2{\mathrm{s}}^3{\mathrm{S}}_1-2{\mathrm{p}}^3{\mathrm{P}}_{0,1,2}^0$ we collected archival X-shooter observations of HD 190073 acquired in the years 2010, 2015, and 2017. It has previously been demonstrated that in magnetic Herbig Ae/Be stars the line parameters of the He I 10 830 Å triplet show cyclic variability related to the stellar magnetic period (e.g. Schöller et al. 2016). The spectral data with this instrument are obtained simultaneously over the entire spectral range from the near-UV to the near-IR in three different arms. The spectral resolution in the red arm is R ≈ 11 000. X-shooter is operated by ESO on the VLT on Cerro Paranal, Chile.

3 Magnetic field measurements and the search for rotational modulation

As in our previous studies using HARPSpol data (see e.g. Järvinen et al. 2018, 2019), to increase the accuracy of the mean longitudinal magnetic field (⟨B_z⟩) determination, we employed the singular value decomposition (SVD) technique, following the description given by Carroll et al. (2012). The mean longitudinal magnetic field is usually determined by computing the first-order moment of the SVD Stokes V profile according to Mathys (1989). The parameters of the lines used to calculate the SVD profiles were taken from the Vienna Atomic Line Database (VALD3; Kupka et al. 2011). Only lines that appear to be unblended or minimally blended in the Stokes I spectra were included in the line mask. The resulting profiles were scaled according to line strength and sensitivity to the magnetic field. The final line mask contains 688 lines and includes Si, Sc, Ti, V, Cr, Mn, Fe, Ni, and Zr lines. The presence of a magnetic field in the SVD profile was evaluated according to Donati et al. (1997), who defined that a Zeeman profile with a false alarm probability (FAP) ≤ 10⁻⁵ is considered as a definite detection, 10⁻⁵ < FAP ≤ 10⁻³ as a marginal detection, and FAP > 10⁻³ as a nondetection.

An inspection of the 65 Narval SVD Stokes I, V, and N profiles acquired in 2011 and 2012 reveals that a few observations are either too noisy and therefore not reliable (2011 Oct. 14, 2012 Jul. 23, and 2012 Nov. 20) or show strong deviations in the profile shapes recorded on consecutive nights (observed in 2011 on Jul. 4, 7, and 8 as well as on Oct. 10, 11, and 12). The noisy profiles were immediately excluded from the analysis, whereas the deviating profiles were excluded at a later point.

Already, a cursory review of the shapes of the hydrogen line profiles in the Stokes I spectra reveals striking differences between the profiles recorded in 2011 and in 2012, probably pointing to a long-term variability caused by an interaction with an unknown component in the HD 190073 system. In Fig. 1, we present the observed overplotted and mean Hβ profiles. The absorption component of the P Cygni profiles of Hβ, related to wind outflow and shifted to the blue by about 200 km s⁻¹, is significantly deeper and the emission component of the profiles is lower in 2012 compared to the profile shapes recorded in 2011. Since HD 190073 possesses a globally organised magnetic field, we assume that the emission component of the profile is related to a magnetosphere around HD 190073. Following the results of the abovementioned studies of this star, we speculate that the observed annual change in the hydrogen line profiles from 2011 to 2012 is due to an interaction with a secondary body in the system. Such an interaction could cause an increase in the wind outflow, producing partial screening of the magnetosphere. With respect to the marked variability of the hydrogen line profiles on a daily timescale, we discuss in Sect. 4 that this variability can be caused by a change of the viewing angle between the observer and the magnetosphere during stellar rotation. Similar changes in hydrogen line profiles over a rotation period are frequently observed in studies of magnetic O and B stars and usually used to determine magnetic periods (e.g. Küker et al. 2024).

Our period search was performed using a nonlinear leastsquares fit to the multiple harmonics utilising the Levenberg–Marquardt method (Press et al. 1992). To detect the most probable period, we calculated the frequency spectrum with a specific number of trial frequencies within the region of interest. A weighted linear least-squares fit was used for each frequency to fit a sine curve and bias offset. A sine wave is expected if the magnetic field is similar to a dipole field inclined with respect to the rotation axis. Based on the result of the fit, we performed a statistical test to check the null hypothesis on the absence of periodicity, that is, to check the statistical significance of the amplitude of the fit (Seber 1977). The resulting F-statistics can be thought of as the total sum, including covariances of the ratio of the harmonic amplitudes to their standard deviations.

Using only the ⟨B_z⟩ measurements obtained from Narval spectra during the years 2011 and 2012, we obtain a period P = 39.6 d (FAP 5.8 × 10⁻⁶), which is in agreement with the period reported by Alecian et al. (2013a). The results of our period search in these data are presented in Fig. A.1. Similar to the study of Alecian et al. (2013a), we do not find any significant period in the data acquired in 2011 (Fig. A.2). For the Narval observations recorded in 2012, we obtain a definite detection of periodicity P = 51.6 d (FAP 9.5 × 10⁻⁵; Fig. A.3). Including in this data set additional ⟨B_z⟩ measurements using HARPSpol and ESPaDOnS observations, an even better FAP (1.3 × 10⁻⁸) refers to the period P = 51.7 d. The corresponding periodogram based on these observations and the phase curve are presented in Fig. 2. The summary of our search for periodicities in our ⟨B_z⟩ measurements using different data sets is presented in Table 1, which also includes the period search based on metal lines described in Sect. 4.1. The ephemeris based on the best visibility of the positive magnetic field extremum is described as $$T_{max(\mathrm{B}\mathrm{s}\mathrm{M}\mathrm{J}\mathrm{D})}=56108.06+51.70(6)\times E.$$(1)

A minimum dipole strength, B_d, of 102 G is estimated for HD 190073 employing the relation B_d ≥ 3|⟨B_z⟩_,all| (Babcock 1958), with a strongest measured ⟨B_z⟩ of 34 ± 3 G. Assuming a rotation period, P, of 51.7 ± 0.06 d and a stellar radius, R_*, of 9.68 ± 0.44 R_⊙ from Guzmán-Díaz et al. (2021), we obtain an equatorial velocity v_eq = 9.47 ± 0.43 km s⁻¹. Using v sin i = 3.19 ± 2.45 km s⁻¹ from the study of Aarnio et al. (2017), the inclination, i, that is, the angle between the rotation axis and the line of sight, is 19.7 ± 15.8°, which is not very accurate, but fully in agreement with interferometric observations (Ibrahim et al. 2023).

The general description for the strength of the observed longitudinal magnetic field for a simple centred dipole was presented by Preston (1967): $$⟨B_{\mathrm{z}}⟩=\frac{1}{20}\frac{15+u}{3-u}{B}_d(\cos \beta \cos i+\sin \beta \sin i\cos 2\pi t/P),$$(2)

where ⟨B_z⟩ is the longitudinal magnetic field, β is the angle between the rotation axis and the magnetic axis, i is the angle between the rotation axis and the line of sight, P is the rotation period of the star, u is the limb-darkening coefficient, and B_d the strength of the dipolar magnetic field. The relative amplitude of variation of the fitted longitudinal magnetic field phase curve is usually characterised by the parameter r, representing the ratio between ⟨B_z⟩_min and ⟨B_z⟩_max. In our case, with ⟨B_z⟩_max = 31 ± 3 G and ⟨B_z⟩_min = −15 ± 3 G, we have r = −0.48 ± 0.11. Using $$r=\frac{{⟨B_{\mathrm{z}}⟩}_{\text{min\hspace{0.17em}}}}{{⟨B_{\mathrm{z}}⟩}_{\text{max\hspace{0.17em}}}}=\frac{\cos \beta \cos i-\sin \beta \sin i}{\cos \beta \cos i+\sin \beta \sin i},$$(3)

we calculated an obliquity angle β = 82.9 ± 6.4°. Assuming a limb-darkening coefficient of 0.5 (Claret 2019), we obtain a polar magnetic field strength of B_d = 222 ± 66 G, typical for Herbig Ae/Be stars (e.g. Hubrig et al. 2015).

Fig. 1 Variability of Hβ profiles in Narval observations. Top panels: overplotted Narval Hβ profiles from 2011 (blue) and 2012 (yellow). Bottom panels: mean Hβ profiles for both years. The plots on the left panels include the full wings of the Hβ line whereas the plots on the right side show the central part of the profiles.

Fig. 2 Periodogram and the magnetic phase curve for the period P = 51.7 ± 0.06 d. Left: periodogram based on all Narval and one ESPaDOnS ⟨Bz⟩ measurements from 2012 and eight HARPSpol ⟨Bz⟩ measurements recorded until 2019. The red line is for observations and the blue line for the window function. Right: ⟨Bz⟩ measurements (magenta circles for Narval and ESPaDOnS and blue circles for HARPS) phased with P = 51.7 ± 0.06 d. The sinusoidal fit is shown with the red line.

4 Line profile variability

The strong variability of the hydrogen emission lines on timescales of days to years in the spectra of HD 190073 has already been reported in a number of studies (e.g. Aarnio et al. 2017; Kozlova et al. 2019; Adigozalzade et al. 2024), but the reason for such a behaviour has remained unclear. Rotational modulation of a magnetically confined wind as the most likely origin for the variability of emission Balmer lines observed in stars possessing magnetic fields was suggested by Donati et al. (2006b). Later on, Sundqvist et al. (2012) calculated the Hα line profile for more than 100 snapshots of a 2D MHD simulation for a magnetic massive star with a dynamical magnetosphere and demonstrated a good overall agreement between their model and observations. The presented simulations strongly supported the assumption that the dynamical magnetospheric model captures the key physics responsible for the Hα variability.

To try to understand the variability character of the hydrogen lines of HD 190073, in addition to the Narval and ESPaDOnS observations from 2011–2012, we downloaded from the Polar-Base archive all available Narval and ESPaDOnS observations distributed over the years 2005–2008, 2013, and 2016–2017. We analysed them together with the HARPSpol spectra acquired between 2011 and 2019. Unfortunately, since only single or very few spectra have been recorded per year, it is impossible to get a complete picture of the variability on the long timescale necessary to understand whether an interaction with a second body in the system indeed takes place and at which time intervals. Our discussion of the observed most extreme changes in the line profiles in the period from 2005 to 2019, and the corresponding figures, are presented in Fig. B.1. In addition, we show the variability observed in the line profiles based on medium-resolution spectra, which are described and discussed in Appendix B and are displayed in Fig. B.2.

Previous studies of highly variable emission line profiles observed in magnetic main-sequence stars demonstrate the great usefulness of dynamical plots to extract important information on the spectral variability over the rotation period (e.g. Küker et al. 2024). Therefore, we used dynamical spectra constructed for hydrogen emission lines observed in 2011 and 2012. To achieve the highest contrast in these spectra, we used differences between individual hydrogen line profiles and the profiles with the lowest core emissions. Even with our refined period, for the data acquired in 2012 the dynamical spectrum presented in Fig. 3 does not show much structure because of the low number of observations and rather large gaps in the phase coverage. Remarkably, as shown in Fig. 4 based on the data acquired in 2011, we detect a clear ringlike magnetospheric structure centred around the phase 0 corresponding to the best visibility of the positive magnetic pole. This is the first snapshot of a magnetosphere around a Herbig Ae/Be star. Taking into account the inclination, i, and the magnetic field geometry with the obliquity, β, of about 83° discussed in Sect. 3, we may also expect to see a magnetospheric ring of lower intensity around the negative magnetic pole at the phase 0.5, but no definite conclusion can be deduced due to the low number of spectra at this rotational phase.

The dynamical spectrum based on the 2011 Narval spectra phased with the period reported by Alecian et al. (2013a) is presented in Fig. D.1. While hinting at a ringlike structure in the variability pattern, it does not reveal this feature in a prominent way.

Fig. 3 Dynamical spectra for the data acquired in 2012. Observed Narval 2012 overplotted line profiles (bottom panels) of Hα (left), Hβ (middle), and Hγ (right), and dynamical spectra of differences between the individual hydrogen line profiles and the profile with the lowest core emission highlighted by the red colour. The phases are based on the period determined in Sect. 3, Prot = 51.7 d.

4.1 Metal lines

Because metal lines in the spectra of HD 190073 are strongly contaminated by emission, for our variability study we selected a few lines of special interest. One of them, the emission O I line at 8446 Å shows a significantly fainter profile in 2012 compared to the profile recorded in 2011. Moreover, in 2006 the emission part of the O I line profile showed three emission peaks of similar intensity (Aarnio et al. 2017), but in 2011 and 2012 the peak in the middle is much lower than the other two peaks. The dynamical spectra and the overplotted profiles for 2011 and 2012 are presented in Fig. C.1. The character of the variability is not very clear, but it is obvious that in 2011 the emission strength was lower at phase 0 at the best visibility of the positive magnetic pole and stronger closer to the negative pole.

It appears that the variability of the intensity of the metal absorption lines is similar to that of the O I 8446 Å line with equivalent widths stronger around the negative pole. As an example, we show in Fig. C.2 the overplotted line profiles, the corresponding periodogram, and the phase curve for the Fe II 6149 Å line observed in 2012 and not strongly affected by emission. This line is rather weak and the intensity changes are rather small. From the periodicity search we obtain P_rot = 51.1 ± 0.6 d, similar to the period determined from magnetic field measurements. No significant period has been detected in the observations acquired in 2011.

4.2 The He I 10 830 Å triplet

As mentioned in Sect. 2, the He I 10 830 Å triplet is an important diagnostic for tracing the structure and kinematics of accretion and outflow in young stars. Due to the nearly face-on disk orientation, we do not expect to detect significant accretion signatures in HD 190073, but should be able to trace the structure of the outflowing wind. On the left side of Fig. 5, we show the P Cygni profiles of this triplet observed in the years 2010, 2015, and 2017, covering the rotation phases from 0.1 to 0.76. The blueshifted absorption component extends in some spectra to more than −400 km s⁻¹, suggesting the presence of a strong stellar wind. It is intriguing that the blueshifted absorption component of the line profile shows a prominent feature in the red wing that persists over the rotation phase ranges 0.1–0.2 and 0.5–0.8 and decreases the contribution of the emission component, probably by strong wind screening. Between these phase ranges, we see an increase in a weaker feature in the blue wing of the absorption component. To better understand the character of the variability, we present on the right side of Fig. 5 a dynamical spectrum of differences between individual He profiles and the profile observed at the rotation phase 0.44, showing the least complex structure. We can see that the feature in the blue wing appears only in phases when the red feature disappears. Since we do not have a good rotation phase coverage, we can only speculate that the observed variability is related to the rotation modulation of the stellar magnetosphere. Without future near-infrared magnetic field measurements using this line, it remains unclear whether the emission component of the He I 10 830 Å triplet comes from the magnetosphere or from magnetically controlled accretion onto the star, with the accretion streams guided by the magnetic field lines of the magnetosphere, or from both. Inverse P Cygni profiles with redshifted absorption reaching several hundred kilometres per second are usually considered as evidence for magnetospheric accretion. These are not observed for HD 190073 at any rotational phase.

Fig. 4 As Fig. 3, but for 2011.

Fig. 5 Variability of the He I 10 830 Å triplet. Left: line profiles of the He I 10 830 Å triplet observed in the X-shooter spectra collected during 2010 (black), 2015 (red), and 2017 (blue). The phases indicated on the right are based on the rotation period of 51.7 d determined in Sect. 3. Right: dynamical spectrum for the He I 10 830 Å triplet. It is constructed using differences between individual He profiles and the profile observed at the rotation phase 0.44 showing the least complex structure.

5 Numerical simulations of the magnetosphere of HD 190073

Knowledge of the fundamental parameters and the dipole field strength of HD 190073 allows us to characterise its magnetosphere. We ran a series of 2D MHD simulations of a line-driven wind originating from HD 190073 using the NIRVANA MHD code (Ziegler 2004, 2005), with a setup similar to that of Küker (2017). The initial magnetic field is a dipole rooted in the star. We used spherical polar coordinates (r,θ) assuming symmetry with respect to the magnetic field axis; to be specific, the colatitude, θ, is measured against that axis. The simulation box is a spherical shell, with the inner boundary at the stellar surface and the outer boundary at 20 R_*. In the meridional direction, the box spans the full interval [0,π], that is, no symmetry with respect to the equatorial plane is assumed.

The main difference between the model of Küker (2017) and the one used here is that we solved the equation for the gas energy as well as the equations for the magnetic field and the gas motion; in other words, the assumption of an isothermal gas has been dropped. We used the cooling curve by Schure et al. (2009) to account for radiative cooling of the gas for temperatures up to 10⁸ K. The cooling rate is $$L=\int n_e{n}_H{\mathrm{\Lambda }}_N(T)dV,$$(4)

where n_e and n_H are the number densities of electrons and protons, respectively. The cooling curve, Λ_N(T), for solar metallicity is tabulated in Schure et al. (2009). For temperatures above 10⁸ K we used the cooling law from Rybicki & Lightman (1979) $$\mathrm{\Lambda }=L_0{T}^{1/2}(1+4.4\times {10}^{-10}T),$$(5)

where the constant, L₀, has been adjusted to produce a smooth transition from the Schure et al. (2009) cooling curve.

The physical parameters used in the simulations are listed in Table 2. The stellar parameters were adopted from Guzmán-Díaz et al. (2021). Catala et al. (2007) found a value of 1.4 × 10⁻⁸ M_⊙ yr⁻¹ for the mass loss rate but stated that this value is an upper limit. We therefore adopted a lower value of 10⁻⁹ M_⊙ yr⁻¹. The value for v_∞ results from an isothermal simulation run without a magnetic field. We note that values for the mass loss rate and terminal velocity refer to the start solution, not the evolved state.

We started the simulations with a spherically symmetric isothermal wind with constant mass flow as in the isothermal case, with a temperature of 10⁵ K rather than the effective temperature of the star. To avoid a discontinuity, we also used this value on the inner boundary. The gas velocity is given by $$v_r=v_{\mathrm{\infty }}{(1-\frac{R}{r})}^{\beta },$$(6)

where R is the stellar radius, v_∞ the terminal velocity from an isothermal model without a magnetic field, and β < 1 a free parameter. In the current case, we have used a value of 0.8.

While the gas flow was allowed to evolve freely, the magnetic field was kept fixed on the surface of the star. The system evolves rapidly in the initial phase of the run but reaches a quasisteady state after a simulation time of about 2 × 10⁵ s. Figure 6 shows the mass density and the magnetic field from a snapshot of the evolved system, taken at t = 5.8 × 10⁵ s. The dipole axis is in the x direction. The magnetic field is still close to the original dipole near the star at low (magnetic) latitudes. Field lines originating at high latitudes have opened up and a current sheet has formed in the equatorial plane. Outflowing gas from the star is trapped in the region of closed field lines and forms a magnetosphere. At the outer edge of the magnetosphere, where the field lines open up, gas flows originating from mid-latitudes converge and merge in the equatorial plane and a disk forms. The top panel of the figure shows an increased thickness of the disk near the outer boundary. This feature, which is also visible on the left side of Fig. 7, is transient, as it moves outwards and leaves the simulation box through the outer boundary. This ring of outwards-moving gas originates from close to the outer edge of the magnetosphere and is accompanied by reconnection in the equatorial plane. In the magnetosphere, the gas moves inwards along the field lines in a kind of zigzag motion, as can be seen in the bottom panel of Fig. 6.

The gas density in the magnetosphere is much higher than in the wind and disk regions of open field lines. The density distribution does not look fundamentally different from the isothermal case, but the temperature distribution does. The temperature of the outflowing gas is much higher than the effective temperature of the star, as shown in the left panel of Fig. 8. It is particularly high in the disk and in the low-density parts of the magnetosphere. The high-density parts of the latter are much cooler though, closer to the effective temperature of the star.

The adopted value for the dipolar magnetic field strength was slightly higher than the value derived from the observations but lies within the error margin. A lower value of 222 G will lead to a slight reduction in the confinement parameter if the values of the mass loss rate and terminal velocity remain the same. Given the uncertainty in the latter two quantities, a change of the dipole magnetic field strength to 222 G has only a minor impact on the results.

We find a smaller magnetosphere than predicted from isothermal models. The extent of the magnetosphere is given by the radius (ud-Doula et al. 2008) $$R_c\approx R_{\ast }+0.7(R_A-R_{\ast }),$$(7)

where the Alfvén radius R_A is given by $$\frac{R_A}{R_{\ast }}\approx 0.3+{\eta }_{\ast }^{1/4}$$(8)

and the confinement parameter, η_*, by $${\eta }_{\ast }=\frac{B_{\mathrm{e}\mathrm{q}}^2{R}_{\ast }^2}{\dot{M}{v}_{\mathrm{\infty }}}.$$(9)

With the values listed in Table 2, we find η_* = 800, R_A/R_* = 5.6, and R_C/R_* = 4.2 for the isothermal case. Using B₀ = 222 G, we would get η_* = 630, R_A/R_* = 5.3, and R_c/R_* = 4.0.

While causing the observed variability, the rotation is too slow to be dynamically important. The impact of rotation is measured by the parameter $$W=\frac{V_{\mathrm{r}\mathrm{o}\mathrm{t}}}{V_{\mathrm{o}\mathrm{r}\mathrm{b}}},$$(10)

where V_rot is the equatorial rotation velocity of the star at its (rotational) equator and V_orb the Keplerian velocity at the same point. We find W = 0.03, which means that the magnetosphere is a dynamical magnetosphere rather than centrifugal. We therefore limit the simulations to the nonrotating case.

The effect of including thermodynamics is demonstrated in Fig. 6, showing closed field lines up to a radius of 3.0 for the last closed field line. The disk also appears less pronounced and less fractured. Reconnection events leading to the ejection of gas from the outer edge of the magnetosphere happen episodically, but less often than in the isothermal case, with rings of gas moving outwards. The mass loss rate is reduced by an order of magnitude from the start solution, while the terminal velocity is increased by a factor of four. The right panel of Fig. 8 shows the spatial distribution of the cooling function. It reflects the density distribution, with values largest in the cool, high-density regions of the magnetosphere. The disk of hot outflowing gas is also quite pronounced.

Our simulations are in the r − θ plane. Figures 6 and 8 therefore show an edge-on view of the system. As the inclination angle, i, is 19.7 ± 15.8 degrees and the obliquity angle, β, has a value of 82.9 ± 6.2 degrees, this view is a good representation of what we observe. For an impression of a pole-on view, we have, however, created a three-dimensional data set using axisymmetry. In the left panel of Fig. 7 we show the gas density in the equatorial plane. The magnetosphere appears as a ring around the star, while the disk forms a fainter halo. The ring shows some substructure, with gaps at certain radii caused by the zigzag flow along the field lines in the magnetosphere. There is also a ring of enhanced gas density in the outer part of the box. This ring is not a static feature but has broken off the inner disk and is moving outwards away from the star.

In the right panel of Fig. 7, we show a zoom-in of the vertically (parallel to the dipole axis) integrated cooling function in the r−ϕ plane. The magnetosphere appears as a bright disk and shows some substructure in form of dark and bright rings and a rather sharp outer edge.

Overall, our simulations show a field geometry, gas distribution, and flow pattern similar to those from previous simulations that assumed the gas to be isothermal at the stellar effective temperature. The inclusion of heat transport and cooling, however, leads to a weaker impact of the magnetic field. The magnetosphere is smaller than in the isothermal case and the region outside the last closed field line shows much less structure. We attribute this to the much higher temperatures, which lead to higher gas pressure and thus larger values of plasma beta, the ratio between gas pressure and magnetic pressure.

Fig. 6 Density distribution and magnetic field geometry for our model of HD 190073. We used a mass loss rate of 10−9 M⊙ yr−1 and a polar magnetic field strength of 250 G. The white lines represent magnetic field lines, the underlying colour contour plot the mass density. Bottom: zoom-in of the magnetosphere.

Fig. 7 Mass density distribution in the equatorial plane (left) and vertically integrated cooling function (right).

Fig. 8 Temperature distribution (left) and cooling function (right) for the same snapshot as Fig. 6.

6 Discussion

Combining our own and archival spectropolarimetric observations obtained between 2012 and 2019 with three different spectropolarimeters, we determine for HD 190073 a rather long magnetic period P = 51.70 ± 0.06 d. No significant period was detected in the data acquired in 2011 between May and November. Monitoring of the variability of the Hβ profiles recorded from 2005 up to 2023 indicates the presence of a significant distortion of the profile shapes, causing striking variability in the depth and position of the blueshifted absorption line component and the intensity of the emission line component. Apart from the variability of the hydrogen line profiles over the rotation cycle, the more substantial changes appear on monthly or annual timescales. Since several studies suggest the presence of a binary or protoplanetary companion in the system HD 190073 (e.g. Baines et al. 2006; Rich et al. 2022; Anilkumar et al. 2024; Ibrahim et al. 2023), it is possible that the observed profile distortions can be explained by a magnetospheric interaction between the magnetic host star and the secondary component.

From our magnetic field measurements in combination with the fundamental parameters presented in the literature, we estimate an obliquity angle β = 82.9 ± 6.4° and a dipole strength B_d = 222 ± 66 G. It is remarkable that the two stars with known magnetic phase curves, HD 101412 and HD 190073, have rather long rotation periods, which are not typical for other Herbig Ae/Be stars, and that in both stars the magnetic axis is located close to the stellar equatorial plane. Based on the work of Romanova et al. (2003), Hubrig et al. (2014) suggested that for HD 101412 with a similar obliquity angle β = 84 ± 13° and a long rotation period of about 42 d, many polar field lines probably thread the inner region of the disk, while the closed lines cross the path of the disk matter, causing strong magnetic braking.

The strong variability of the hydrogen emission lines Hα, Hβ, and Hγ has been studied using dynamical spectra. Due to the rather small number of acquired Narval spectra in 2012 July–October, our dynamical spectra do not reveal any clear periodic features. However, in the data acquired in 2011 May–November, we detect for all three hydrogen emission line profiles clear ringlike magnetospheric structures appearing at the rotation phase 0 during the best visibility of the magnetic pole. These spectra present the first snapshot of a magnetosphere around a Herbig Ae/Be star. Our study of the variability of metal lines indicates that their strength varies in an opposite way to the strength of the hydrogen emission lines, with stronger lines observed closer to the negative magnetic pole. Intriguingly, variable absorption features changing their shape over the rotational period have been detected in the wings of the blueshifted absorption component of the He I 10 830 Å triplet. We speculate that these features are related to the rotational modulation of the observed magnetosphere.

For the first time, 2D MHD simulations have been carried out for a Herbig Ae/Be star. Using the NIRVANA MHD code and involving nonisothermal gas, our simulations show that the magnetosphere should be rather compact with a radius of about three stellar radii, but the wind flow extends over tens of stellar radii. With the radius of a disk of 1.14 au around HD 190073 (Ibrahim et al. 2023), the distance between the star and the disk is about 25 stellar radii. It is not clear whether an interaction between the stellar magnetosphere and the surrounding accretion disk would explain the sub-au structure inconsistent with Keplerian motion, as Ibrahim et al. (2023) report, and their speculation about the presence of an object embedded in the inner disk.

As mentioned in Sect. 1, several studies of HD 190073 indicate the possible presence of a lower-mass binary companion very close to the host star and a more distant giant planet component, making this system a valuable laboratory to study magnetic interaction between a host star and its companions. There are currently two ways to confirm the reported discoveries of companions. Studying the presence of a companion within inner disks of Herbig Ae/Be stars, a resolution below one milliarcsecond is required, which may be achievable with a new generation of optical interferometric instruments with longer baselines. Current near-infrared interferometric surveys of the inner regions in Herbig Ae/Be stars offer a spatial resolution of the order of one milliarcsecond. Follow-up polarized differential imaging observations with the Gemini Planet Imager are also needed to confirm that the detected giant planet and the host star are co-moving.

Acknowledgements

We thank the referee for their comments. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 187.D-0917, 089.D-0383, 097.C-0277, 099.C-0081, 0103.C-0240, 385.C-0131, 095.C-072, and 099.C-0342. Based on observations collected at the Bernard Lyot Telescope at the Pic du Midi de Bigorre, France, which is managed by the Observatoire Midi Pyrénées, and on observations collected at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii.

Appendix A Period analysis

Our period search based on 2011 and 2012 Narval spectra combined is presented in Fig. A.1 and separately in Figs. A.2 and A.3, respectively.

Fig. A.1 Periodogram based on all Narval measurements from the years 2011 and 2012 (left). The red line is used for the observations and the blue line for the window function. ⟨Bz⟩ measurements phased with P = 51.6 ± 1.1 d (right). The sinusoidal fit is shown with the red line.

Fig. A.2 As Fig. A.1, but for the Narval observations in 2011. The ⟨Bz⟩ measurements are phased with P = 56.4 ± 2.7 d and the sinusoidal fit is shown with the red line.

Fig. A.3 As Fig. A.1, but for the Narval observations in 2012. The ⟨Bz⟩ measurements are phased with P = 51.6 ± 1.1 d and the sinusoidal fit is shown with the red line.

Appendix B Variability of hydrogen lines

In Fig. B.1, we present the individual and mean (for details, see Table B.1) line profiles of Hβ obtained over different years from 2005 to 2019. In the following, we discuss the changes of these profiles over time.

The Hβ profiles between 2005 and 2008 show very little variability. The line emission component is the weakest in the 2007 Narval profile (light green). This profile shows also the strongest red wing absorption component among all profiles. The strongest absorption in the blue wing is detected in the 2005 ESPaDOnS profile (purple) and the 2008 Narval profile (black). After 2008, the profiles showed much more variability. Among all the profiles, the emission component is strongest in the 2017 HARPSpol spectrum (magenta), followed by the 2016 HARPSpol profile (brown), and the 2009 Narval profile (red). The blueshifted absorption component is strongest in the 2019 HARPSpol spectrum (navy blue). It also shows one of the weakest emission components, comparable to that of the Narval 2007 profile, albeit the red emission peak is stronger in 2019 than in 2007. Both the 2012 ESPaDOnS and Narval spectra (light and dark yellow, respectively) show similarly deep absorption line components, but clearly weaker than seen in 2019. Also the 2016 HARPSpol profile has almost equally strong absorption in the blue wing as was detected in 2012. The red absorption wing shows very little variability — the strongest absorption is seen in the 2007 Narval profile and the weakest in the 2019 HARPSpol profile.

Fig. B.1 Variability of Hβ profiles over different years. Top: the observations have been acquired from 2005 to 2019 using different instruments: E indicates ESPaDOnS, N Narval, and H HARPSpol. Left: larger velocity range to show the full extent of the line profiles. Right: zoom-in to the line core. The profile in magenta with the strongest emission contribution refers to the HARPSpol observation in 2017 whereas the profile in navy blue colour highlights the HARPSpol observation in 2019. The middle plot shows the profiles until 2008 and the bottom plot the profiles of the later years.

Similar to Fig. B.1, in Fig. B.2, we present individual and mean (for details, see Table B.1) line profiles of the Hβ line based on medium-resolution observations. The observations were obtained with the Cassegrain Echelle Spectrograph (CES) at the Special Astrophysical Observatory (SAO) of the Russian Academy of Sciences in 2015. After 2015, the observations were carried out using the ShAO Fiber-Echelle Spectrograph (ShAFES), developed jointly by specialists from ShAO and SAO. For more details, we refer to Adigozalzade et al. (2024). From 2015, there are 22 spectra with a resolution R ≈ 14 000. After 2015, the resolution is higher (R ≈ 28 000), but the number of spectra varies between one and seven per year. It has already been discussed by Adigozalzade et al. (2024) that the individual Hβ line profiles obtained during 2015 show variability from night to night. The mean profile of 2015 (light blue) has the strongest emission component and the weakest absorption wings on both sides. The 2016 profile (brown) shows the weakest emission component whereas the absorption wings show an average behaviour.

Fig. B.2 As Fig. B.1, but based on medium-resolution spectra. The plot shows the variability of the Hβ profiles over different years from 2015 to 2023. The mean profile based on 2015 spectra is plotted with a dashed line since the resolution in those spectra was lower. It also has the strongest emission contribution. The deepest blue wing absorption is seen in the mean profile from 2020 (dark green). However, also the 2022 profile (salmon) and the mean of the 2023 profiles (dark yellow) have a deep absorption in the blue wings.

The profile of 2020 (dark green) shows the strongest absorption in the blue wing. The single 2022 profile (salmon) and the average of 2023 (light yellow) also have both very deep absorptions in the blue wing. Also the emission line components are among the strongest, but the shapes are very different. In the 2022 profile the red emission peak is much higher than the blue one, whereas in the 2023 profile the difference in height is not that large, with the redder peak still being stronger.

Appendix C Variability of metal lines

The metal lines of HD 190073 are strongly contaminated by emission and show also variability over time. The dynamical spectra and overplotted Narval profiles for 2011 and 2012 of the O I line at 8447 Å are shown in Fig. C.1. The variability seen in O I is not as clear as in the hydrogen lines, but still the 2011 data set shows that the emission strength is lower at phase 0 at the best visibility of the positive magnetic pole.

Fig. C.1 As Fig. 3, but for the O I 8446 Å line profiles. Dynamical spectra for the years 2011 (left) and 2012 (right) are constructed using differences between individual O I 8446 Å line profiles and the line profile with the lowest core emission (shown in red).

In the case of the metal absorption lines, the scale of variations is even smaller, but the equivalent widths of Fe II 6149 Å vary similar to the OI line. The line and EW variability of the Fe II line in 2011 Narval spectra is presented in Fig. C.2. The EWs were also used for a period search, which led to P_rot = 51.1 ± 0.6 d. This period is slightly shorter than the one obtained from the mean longitudinal magnetic field measurements P_rot = 51.7 ± 0.06 d, but is comparable within the errors.

Fig. C.2 Variability of the Fe II 6149 Å lines. Top: 2012 Narval line profiles overplotted for Fe II 6149 Å. Middle: periodogram based on EWs measured from the profiles shown in the top panel. The red line is for the observations and the blue line for the window function. Bottom: magenta circles represent EWs measured from the 2012 Narval spectra using the Fe II 6149 Å line and are phased with P = 51.1 d. The obtained period is slightly shorter than the one obtained from the mean longitudinal magnetic field measurements, but within the errors still similar. The sinusoidal fit is shown with the red line.

Appendix D Dynamical spectrum using all 2011 Narval observations

In Fig. D.1 we show the dynamical spectrum for the Hβ line based on Narval observations obtained in 2011 and phased with the period reported by Alecian et al. (2013a).

Fig. D.1 As Fig. 3, but using Narval observations obtained in 2011 and the period reported by Alecian et al. (2013a).

References

All Tables

All Figures

	Fig. 1 Variability of Hβ profiles in Narval observations. Top panels: overplotted Narval Hβ profiles from 2011 (blue) and 2012 (yellow). Bottom panels: mean Hβ profiles for both years. The plots on the left panels include the full wings of the Hβ line whereas the plots on the right side show the central part of the profiles.
In the text

Fig. 2 Periodogram and the magnetic phase curve for the period P = 51.7 ± 0.06 d. Left: periodogram based on all Narval and one ESPaDOnS ⟨Bz⟩ measurements from 2012 and eight HARPSpol ⟨Bz⟩ measurements recorded until 2019. The red line is for observations and the blue line for the window function. Right: ⟨Bz⟩ measurements (magenta circles for Narval and ESPaDOnS and blue circles for HARPS) phased with P = 51.7 ± 0.06 d. The sinusoidal fit is shown with the red line.

In the text

	Fig. 3 Dynamical spectra for the data acquired in 2012. Observed Narval 2012 overplotted line profiles (bottom panels) of Hα (left), Hβ (middle), and Hγ (right), and dynamical spectra of differences between the individual hydrogen line profiles and the profile with the lowest core emission highlighted by the red colour. The phases are based on the period determined in Sect. 3, Prot = 51.7 d.
In the text

	Fig. 4 As Fig. 3, but for 2011.
In the text

Fig. 5 Variability of the He I 10 830 Å triplet. Left: line profiles of the He I 10 830 Å triplet observed in the X-shooter spectra collected during 2010 (black), 2015 (red), and 2017 (blue). The phases indicated on the right are based on the rotation period of 51.7 d determined in Sect. 3. Right: dynamical spectrum for the He I 10 830 Å triplet. It is constructed using differences between individual He profiles and the profile observed at the rotation phase 0.44 showing the least complex structure.

In the text

	Fig. 6 Density distribution and magnetic field geometry for our model of HD 190073. We used a mass loss rate of 10−9 M⊙ yr−1 and a polar magnetic field strength of 250 G. The white lines represent magnetic field lines, the underlying colour contour plot the mass density. Bottom: zoom-in of the magnetosphere.
In the text

	Fig. 7 Mass density distribution in the equatorial plane (left) and vertically integrated cooling function (right).
In the text

	Fig. 8 Temperature distribution (left) and cooling function (right) for the same snapshot as Fig. 6.
In the text

	Fig. A.1 Periodogram based on all Narval measurements from the years 2011 and 2012 (left). The red line is used for the observations and the blue line for the window function. ⟨Bz⟩ measurements phased with P = 51.6 ± 1.1 d (right). The sinusoidal fit is shown with the red line.
In the text

	Fig. A.2 As Fig. A.1, but for the Narval observations in 2011. The ⟨Bz⟩ measurements are phased with P = 56.4 ± 2.7 d and the sinusoidal fit is shown with the red line.
In the text

	Fig. A.3 As Fig. A.1, but for the Narval observations in 2012. The ⟨Bz⟩ measurements are phased with P = 51.6 ± 1.1 d and the sinusoidal fit is shown with the red line.
In the text

Fig. B.1 Variability of Hβ profiles over different years. Top: the observations have been acquired from 2005 to 2019 using different instruments: E indicates ESPaDOnS, N Narval, and H HARPSpol. Left: larger velocity range to show the full extent of the line profiles. Right: zoom-in to the line core. The profile in magenta with the strongest emission contribution refers to the HARPSpol observation in 2017 whereas the profile in navy blue colour highlights the HARPSpol observation in 2019. The middle plot shows the profiles until 2008 and the bottom plot the profiles of the later years.

In the text

Fig. B.2 As Fig. B.1, but based on medium-resolution spectra. The plot shows the variability of the Hβ profiles over different years from 2015 to 2023. The mean profile based on 2015 spectra is plotted with a dashed line since the resolution in those spectra was lower. It also has the strongest emission contribution. The deepest blue wing absorption is seen in the mean profile from 2020 (dark green). However, also the 2022 profile (salmon) and the mean of the 2023 profiles (dark yellow) have a deep absorption in the blue wings.

In the text

	Fig. C.1 As Fig. 3, but for the O I 8446 Å line profiles. Dynamical spectra for the years 2011 (left) and 2012 (right) are constructed using differences between individual O I 8446 Å line profiles and the line profile with the lowest core emission (shown in red).
In the text

Fig. C.2 Variability of the Fe II 6149 Å lines. Top: 2012 Narval line profiles overplotted for Fe II 6149 Å. Middle: periodogram based on EWs measured from the profiles shown in the top panel. The red line is for the observations and the blue line for the window function. Bottom: magenta circles represent EWs measured from the 2012 Narval spectra using the Fe II 6149 Å line and are phased with P = 51.1 d. The obtained period is slightly shorter than the one obtained from the mean longitudinal magnetic field measurements, but within the errors still similar. The sinusoidal fit is shown with the red line.

In the text

	Fig. D.1 As Fig. 3, but using Narval observations obtained in 2011 and the period reported by Alecian et al. (2013a).
In the text