© ESO, 2014

1. Introduction

The bright, late-B star CU Vir (HR 5313, HD 124224, HIP 69389) is one of the best known intermediate-mass magnetic, chemically peculiar (CP) stars. This object shows abnormally strong lines of Si, weak lines of He and a notable spectroscopic (Deutsch 1952) and photometric (Hardie 1958) variability with a period of ~0.5 d. These periodic changes arise due to a combination of stellar rotation and an inhomogeneous distribution of chemical elements over the stellar surface. It is generally believed that high-contrast chemical spots are produced by an anisotropic atomic diffusion in the presence of a strong, globally-organised magnetic field (Michaud et al. 1981; Alecian & Vauclair 1981; Alecian & Stift 2010), although a detailed diffusion theory capable of explaining the surface structure of individual magnetic CP stars is still lacking.

Being one of the most rapidly rotating magnetic CP stars, CU Vir represents an optimal target for the reconstruction of the surface chemical spot distributions with the Doppler imaging (DI) method. Using this technique, maps of He, Mg, Si, Cr, and Fe were obtained for CU Vir by different authors (Goncharskii et al. 1983; Hiesberger et al. 1995; Hatzes 1997; Kuschnig et al. 1999). The study by Krtička et al. (2012) showed that these uneven chemical distributions provide a successful quantitative explanation for the observed spectrophotometric behaviour of CU Vir in the ultraviolet and optical wavelength regions.

The fast rotation and stability of the photometric light curves of CU Vir enable a very precise period determination and allows one to study subtle long-term variations of stellar rotation. Consequently, CU Vir is one of a few stars for which rotational period changes have been securely detected. Remarkably, this star appears to show period glitches (Pyper et al. 1998, 2013) or possibly a cyclic period variation (Mikulášek et al. 2011) rather than a constantly increasing period, as expected from magnetic spin down models (Ud-Doula et al. 2009). The origin of these period changes is currently unknown.

In addition to the photometric and spectroscopic variability typical of magnetic CP stars, CU Vir exhibits a unique radio emission signature, which is not seen in any other star. Trigilio et al. (2000) reported two sharp 100 percent circularly polarised radio pulses occurring at particular rotational phases. This emission was subsequently studied in different radio bands (Trigilio et al. 2008, 2011; Ravi et al. 2010; Lo et al. 2012) and interpreted as an electron cyclotron maser emission that originates in the stellar magnetosphere above one of the magnetic poles (see Lo et al. 2012, and references therein). The sharp radio emission pulses can also be used for a very precise measurement of the stellar rotational period. However, it is not clear if the phase shifts observed in the radio are related to the intrinsic variation of the stellar rotation or to instabilities in the region where the radio emission is produced (Ravi et al. 2010; Pyper et al. 2013).

Evidently, the magnetic field is a key parameter for understanding different variability phenomena observed in CU Vir. The strength and geometry of the surface magnetic field influences the photospheric distribution of chemical elements, while the structure of the magnetosphere at a distance of 2–3 stellar radii is the main ingredient of any possible non-thermal radio emission process. The presence of a kG-strength magnetic field in CU Vir was established by Landstreet & Borra (1977) and Borra & Landstreet (1980) using a Balmer line magnetometer. Their mean longitudinal magnetic field (⟨ B_z ⟩) measurements revealed a smooth, roughly sinusoidal variation with an amplitude of ≈600 G and reversing sign. Several oblique dipolar models were derived using these ⟨ B_z ⟩ data (Borra & Landstreet 1980; Trigilio et al. 2000), resulting in estimates of the polar field strength B_d = 3.0−4.5 kG and magnetic obliquity β = 70−90°, depending on the assumed inclination angle i. On the other hand, Hatzes (1997) and Glagolevskij & Gerth (2002) argued that the magnetic field geometry of CU Vir is better described by an offset dipole or by a combination of dipolar and quadrupolar components of comparable strengths. Observations of the mean longitudinal magnetic field available in the literature are insufficient to definitively distinguish between these alternative magnetic field models.

The goal of our investigation is to obtain a much more precise and detailed picture of the magnetic field topology of CU Vir by using new high-resolution circular polarisation observations. To interpret these data, we developed a new version of the Zeeman Doppler imaging (ZDI) methodology, which for the first time combines detailed polarised spectrum synthesis calculations based on realistic model atmospheres with a multi-line approach that is widely used in stellar spectropolarimetry.

The rest of this paper is organised as follows. Section 2 presents spectropolarimetric observations and outlines the corresponding analysis methodology, including calculation of the mean polarisation profiles and their interpretation in terms of the magnetic and chemical stellar surface maps. In Sect. 3, we present revised stellar parameters, discuss new longitudinal field measurements, and derive magnetic field topology and chemical abundance distributions with the help of ZDI inversions. The results of our investigation are summarised and discussed in the context of other recent studies in Sect. 4.

2. Methods

2.1. Spectropolarimetric observations

We collected 24 spectropolarimetric observations of CU Vir in the context of the Magnetism in Massive Stars (MiMeS) large program. The Stokes I and V spectra were obtained over a time span of 26 months in 2009–2011 with the Narval spectropolarimeter that is attached to the 2-m Télescope Bernard Lyot at the Pic du Midi Observatory. Narval is a fibre-fed, cross-dispersed echelle spectropolarimeter, covering the wavelength range from 3694 to 10 483 Å at a resolution of λ/ Δλ = 65 000.

The spectra of CU Vir were processed with the Libre-Esprit reduction package (Donati et al. 1997). The extracted echelle orders were subsequently renormalised to the continuum and merged using a set of dedicated IDL routines. The signal-to-noise ratio (S/N) of the resulting spectra varied between 700 and 1500 per extracted pixel, as measured at λ = 5400 Å.

Each Narval observation of CU Vir consisted of a sequence of four sub-exposures obtained at two different orientations of the quarter wave retarder plate relative to the beam splitter. Combining the resulting left- and right-hand polarised spectra, as discussed by Donati et al. (1997) and Bagnulo et al. (2009), allows one to minimise spurious polarisation and to derive a diagnostic null spectrum, which is useful for identifying residual instrumental artefacts.

The rotational phases of all CU Vir observations analysed in our study were evaluated with the help of the variable-period ephemeris suggested by Mikulášek et al. (2011). This allowed us to phase the Narval data with the earlier magnetic field measurements by Borra & Landstreet (1980). We note that the period variability discussed by Mikulášek et al. (2011) and others in recent publications has a negligible impact on the relative phases of the Narval spectra, since our observations were collected during a relatively short period. The mean value of the rotational period during Narval observations, 0.5207130 d, is also identical within error bars to the constant period of 0.5207137 ± 0.0000010 d, which is proposed by Pyper et al. (2013) for the time period 1993–2012.

Information on the observing dates, heliocentric Julian dates, rotational phases, and S/N of individual Narval observations of CU Vir is provided in the first four columns of Table 1.

2.2. Multi-line polarisation analysis

Owing to a large rotational Doppler broadening of the lines in the spectrum of CU Vir and to a moderate strength of its magnetic field, the circular polarisation signatures can be detected at a significance level of 3–5σ in only a few very strong metal lines and for a limited range of rotational phases even with our very high quality spectropolarimetric observations. To enhance the analysis of the Stokes V spectra, we applied a multi-line technique of least-squares deconvolution (LSD, Donati et al. 1997; Kochukhov et al. 2010), which combines hundreds of spectral lines into high S/N mean profiles.

The line mask necessary for application of LSD was obtained from the VALD data base (Piskunov et al. 1995; Kupka et al. 1999), using the stellar parameters discussed below (see Sect. 3.1). In total, we used 359 spectral lines with a central depth (before instrumental and rotational broadening) of ≥10 per cent of the continuum. The lines were chosen in the 4000–6865 Å wavelength interval, avoiding the regions affected by the broad hydrogen Balmer line wings. Three sets of the Stokes I and V profiles were calculated: the one based on all metal lines, the mean Si profile (based on 63 lines), and the mean Fe profile (based on 202 lines). The LSD profiles derived from the full line mask were used for the longitudinal field measurements, while the mean Si and Fe profiles were employed for the magnetic inversions. All LSD profiles were calculated within ± 350 km s^-1 of the laboratory line positions, using 3 km s^-1 velocity bins.

The LSD profiles for specific elements were recovered using the multi-profile version of the LSD method discussed by Kochukhov et al. (2010). This means that, in practice, we calculated LSD profiles for a given element taking the blending of its lines by all other elements into account. In all cases, we also calculated mean profiles from the diagnostic null spectra, applying the same line masks and using the same wavelength regions as in the Stokes V LSD analysis. No evidence of spurious polarisation signatures was found.

The application of LSD led to more than a tenfold increase of the S/N of the Stokes V spectrum, which reaches a typical value of 2 × 10⁴ per velocity bin (see 5th column in Table 1). This enabled a secure detection of the circular polarisation signatures for all rotational phases. An example of the LSD Stokes V profile obtained using the full metal line mask is presented in Fig. 1. This particular observation, corresponding to the rotational phase 0.424, yields the highest amplitude circular polarisation signature detected in our data. Interestingly, the Stokes V amplitude reaches maximum for the rotational phase interval 0.3–0.5, when the mean longitudinal magnetic field changes sign. During this part of the rotation cycle, the Stokes V profile has a predominantly symmetric shape, indicating the presence of at least two regions of different field polarity on the stellar surface. Such symmetric Stokes V profiles cannot be assessed by the spectropolarimetric techniques, which are limited to ⟨ B_z ⟩ measurements.

2.3. Zeeman Doppler inversions with realistic interpretation of LSD profiles

Numerous previous studies demonstrated the LSD to be a remarkably powerful method for detecting weak stellar magnetic fields, for precise measurements of the mean longitudinal field, and, potentially, other field moments. However, the usage of LSD profiles for recovering the physical information about the stellar atmosphere or geometrical information on the surface distributions of temperature, chemical abundance, and magnetic field is more controversial. All types of modelling of LSD profiles attempted so far have relied on a number of simplifying, often severe, assumptions about the nature of the mean profiles and their relation to the mean line parameters that are inferred from the LSD line mask.

In practice, the interpretation of LSD profiles in the context of DI or ZDI studies proceeds as follows (see Fig. 2a). First, one adopts appropriate stellar atmospheric parameters and calculates the intensities of spectral features contributing to the LSD line mask. Then LSD profiles are derived for this mask and for given observational data. This step also provides weighted mean values for different atomic parameters that correspond to the LSD mask. For the next step, one assumes in one way or another that LSD profiles behave similarly to a single spectral feature with the previously found mean parameters. In particular, the response to the magnetic field is assumed to be that of a Zeeman triplet with the mean effective Landé factor. Using this approximation, local line profiles are computed. Disk integration is performed, and a comparison between the observed and computed LSD profiles is carried out.

Previous studies differ in the details of theoretical calculations of the local LSD profiles. Most ZDI investigations of cool stars (e.g., Petit et al. 2004; Marsden et al. 2011) assume that the local Stokes I profile is represented by a temperature-independent Gaussian function, and the local Stokes V profile is given by the derivative of the Gaussian, according to the weak field approximation. It is not entirely clear how the width and intensity of this Gaussian profile are chosen, especially when the observed LSD profiles are strongly distorted by cool spot signatures. A handful of studies (Donati et al. 2006; Morin et al. 2008) used a more sophisticated Unno-Rachkovsky analytical solution of the polarised radiative transfer (PRT) equation for calculating local LSD profiles. This method does not rely on magnetic field being weak but, again, requires a number of input line parameters that cannot be determined uniquely. Finally, Folsom et al. (2008) and Kochukhov et al. (2013) attempted to approximate the local LSD profiles with a calculation based on a realistic stellar model atmosphere and numerical PRT modelling. In this case, it is assumed that the local LSD profile behaves as a single line of the dominant ion (typically Fe i or Fe ii) with some parameters (wavelength, excitation potential, and Landé factor) determined from the LSD mask, and other parameters (oscillator strength and damping) adjusted to match the LSD profile shape that is obtained from a full theoretical spectrum synthesis calculation.

Notwithstanding different degrees of sophistication for various LSD profile modelling approaches used so far, the single-line approximation represents the most significant limitation of the LSD methodology in the context of DI and ZDI. Given a wide range of spectral line responses to temperature, abundance, and magnetic field, it is unclear if any single line can match an average profile constructed from hundreds or thousands absorption features. To this end, Kochukhov et al. (2010) showed that the validity range of the single-line LSD methodology is very limited even with the most detailed PRT interpretation. The LSD Stokes I and V profiles can be reasonably matched by an adjusted synthetic line only up to the field strength of ≈2 kG. The Stokes Q and U profiles cannot be matched at any field strength. Furthermore, the response of the Stokes I LSD profile to chemical abundance variation differs drastically from that of any single line, making an accurate chemical analysis or a realistic abundance DI mapping impossible with LSD profiles. A similar problem affects reconstruction of temperature maps from the observed LSD profiles of cool active stars.

2.3.1. Theoretical local LSD profiles

In this study, we have developed a new method of interpreting LSD profiles that overcomes the main limitations described above. We abandon the single-line LSD profile interpretation in favour of deriving the local LSD profiles from the full spectrum synthesis calculations. This approach, recently applied by Kochukhov & Sudnik (2013) for forward Stokes V profile calculations, is schematically illustrated in Fig. 2b. We start by generating a library of theoretical four Stokes parameter spectra, including all relevant absorption lines and covering the entire observed wavelength range, for a grid of abundances/temperatures, magnetic field strengths, orientations of the field vector with respect to the line of sight, and the limb angles. These spectra are convolved with a Gaussian instrumental profile. Then, LSD profiles are calculated for each of these theoretical spectra using the same line mask as applied to the observed spectra. The relative wavelength dependent weights are taken from a representative observation. Simultaneously, we also determine the mean continuum intensity, which retains only the abundance/temperature and limb angle dependence for the case of unpolarised continuum. The local LSD profiles and continuum intensities for arbitrary parameters are then determined by a linear interpolation within the resulting theoretical LSD profile data base. After interpolation, the Stokes Q and U profiles are transformed (rotated) according to the local orientation of the field vector in the plane perpendicular to the line of sight. In the final step, we perform the disk integration and compare with observations.

The new LSD profile interpretation methodology attributes no particular meaning to the mean profiles. In particular, we do not assume that they behave as a single line, require them to be smooth and have certain symmetry properties, nor impose any limits on the field strength. In our approach LSD is merely a method to compress information and numerically compare observations and theoretical calculations. For reasons of numerical tractability, the interpretation of LSD Stokes profiles suggested here assumes that it is possible to interchange the operations of disk integration and LSD. This approximation is exact given the linear nature of the disk integration. Its validity has been verified numerically with the help of forward calculations for different star spot scenarios (magnetic, temperature, and abundance inhomogeneities).

In the present study of CU Vir, we used the LLmodels code (Shulyak et al. 2004) to generate model atmospheres for two independent grids of Si and Fe abundances. For both elements, calculations covered a range from −2.0 to + 2.5 relative to the solar abundance (log  N_Si/N_tot = −4.53 and log N_Fe/N_tot = −4.54, Asplund et al. 2009) with a step of 0.25 dex. Abundances of other elements were either set to the mean values of the DI maps derived by Kuschnig et al. (1999) or assumed to be solar. Theoretical spectra were calculated with the help of the SynMast PRT code (Kochukhov 2007; Kochukhov et al. 2010) using a line list containing 2185 lines. This includes all spectral features deeper than 1% retrieved from the VALD data base for the region of interest. Stokes profiles were evaluated for 17 values of the magnetic field strength between 0 and 7 kG using a step of 250 G for the 0–1 kG interval and a step that is twice as large for B> 1 kG. At each field strength value and for each atmospheric model, calculations were performed for a 15 × 15 grid of limb angles and magnetic vector orientations with respect to the line of sight. Theoretical LSD profiles were evaluated for 72 675 resulting model Stokes spectra using the same multi-component LSD approach as applied to the observational data.

An example of the synthetic LSD Stokes I profiles for three different Si abundances is shown in Fig. 3. We illustrate an attempt to reproduce these profiles with the single-line spectrum synthesis calculations that are adjusted to match the intermediate (solar-abundance) Si LSD profile. Even after such an adjustment, the single-line spectrum synthesis does not reproduce all details of the reference theoretical LSD profile. Moreover, when the Si abundance is varied by ± 1 dex, the single-line calculation fails to correctly predict the LSD profile intensity owing to its much steeper abundance dependence than shown by the theoretical LSD profiles.

2.3.2. Magnetic and abundance mapping

The procedure of reconstructing the magnetic field topology and chemical spot maps adopted in our study generally follows that of the Invers13 code (Kochukhov et al. 2012, 2013) with the exception that the local Stokes spectra are obtained by a linear interpolation within the LSD profile table, which is computed as described above rather than calculated on the fly. The stellar surface is divided into 1876 zones with roughly equal areas (see Piskunov & Kochukhov 2002). A chemical map is specified by the element abundance values in each of these surface zones. On the other hand, the maps of radial, meridional, and azimuthal magnetic field components are specified with the help of a general spherical harmonic expansion (Donati et al. 2006), where we use the real spherical harmonic functions to describe the mode with degree ℓ and order mwith (7)and (8)Here θ and φ are the colatitude and longitude at the stellar surface, and P_ℓ,m(θ) is the associated Legendre polynomial. Its derivative with respect to θ, required to evaluate Z_ℓ,m(θ,φ), can be computed with (9)The free parameters of this harmonic expansion are coefficients α_ℓ,m, β_ℓ,m, and γ_ℓ,m. They characterise contributions of the radial poloidal, horizontal poloidal, and horizontal toroidal magnetic field components, respectively. The expansion is carried out up to a sufficiently large ℓ_max, as to fit the details of the observed Stokes profiles. In the present study of CU Vir, we used ℓ_max = 10, yielding 360 independent magnetic variables.

An iterative reconstruction of the abundance and magnetic maps, aimed to fit the observed Stokes I and V profiles, is carried out simultaneously and self-consistently, using the modified Levenberg-Marquardt optimisation algorithm (Piskunov & Kochukhov 2002). For the regularisation, which is a necessary ingredient of any Doppler imaging code, we use the Tikhonov functional (e.g. Tikhonov & Arsenin 1977) for the abundance map and a penalty function for the magnetic field with corresponding individual regularisation parameters. This means that we constrain the inversion by requiring the final abundance map to be as smooth as allowed by observations and the magnetic map to contain the least possible contribution of the higher-order harmonic modes. Alternatively, the magnetic penalty function can be substituted with a conservative field parameterisation. For example, the classical oblique dipolar field is realised by setting α_ℓ,m = β_ℓ,m and γ_ℓ,m = 0 for ℓ_max = 1.

In comparison to the direct reconstruction of magnetic field practiced in our previous ZDI studies of Ap stars (Kochukhov et al. 2004; Kochukhov & Wade 2010), the usage of spherical harmonic expansion has the advantage of yielding the null magnetic flux through the stellar surface and hence automatically satisfying the Maxwell’s equations. In addition, since individual terms of the expansion are orthogonal, one can conveniently quantify contributions of different spatial scales to the stellar magnetic field topology by comparing the relative energies (B² integrated over the stellar surface) of harmonic components. The disadvantage of the harmonic expansion and the associated penalty function is that this method is difficult to apply to stars showing fine-structured bipolar field geometries. However, this problem concerns only certain types of active cool stars rather than early-type magnetic stars, such as the one studied here.

3. Results

3.1. Stellar parameters

Determination of the atmospheric parameters of magnetic CP stars is notoriously difficult due to non-solar chemical abundance patterns and variability. Using standard photometric T_eff calibrations trained on normal stars may lead to significant systematic errors, especially for the Si-type CP stars, such as CU Vir (Netopil et al. 2008). In this situation, it is preferable to rely on more detailed model atmosphere analyses, which attempt to account for chemical peculiarity, at least to some extent. To this end, the most relevant previous studies of CU Vir are those by Kuschnig et al. (1999), Shulyak et al. (2004), and Lipski & Stȩpień (2008). The first of these publications suggested T_eff = 13 000 K and log  g = 4.0 based on spectroscopy. This set of parameters was also used by Krtička et al. (2012) in the modelling of the optical and ultraviolet spectral energy distribution (SED) of CU Vir. On the other hand, Shulyak et al. (2004) and Lipski & Stȩpień (2008) directly fitted the observed SED, finding T_eff = 12 750 K and log  g = 4.0–4.5. Taking these results into account, we adopted T_eff = 12 750 ± 250 K in the present paper.

Combining T_eff with the observed V magnitude, the trigonometric parallax 12.63 ± 0.21 mas (van Leeuwen 2007), and an empirical bolometric correction BC = − 0.79 ± 0.1 (Lipski & Stȩpień 2008), we found the stellar luminosity to be L = 100 ± 11 L_⊙ and the radius R = 2.06 ± 0.14 R_⊙. Then, a comparison with the Padova theoretical stellar evolutionary tracks (Bertelli et al. 2009) allows us to place CU Vir in the H-R diagram. As illustrated in Fig. 4, the parameters of our target correspond to the stellar mass M = 3.06 ± 0.06 M_⊙. The surface gravity is log  g = 4.30 ± 0.06 and the stellar age, is constrained to be ≤100 Myr. The stellar radius obtained here is consistent with the values R = 2.3 ± 0.1 R_⊙ (Kochukhov & Bagnulo 2006) and R = 2.2 ± 0.2 R_⊙ (North 1998), which are derived from photometry and evolutionary tracks. A similar radius of R = 1.9 ± 0.1 R_⊙ was estimated by Krtička et al. (2012) from the observed UV flux. Lipski & Stȩpień (2008) also reported R = 2.2 R_⊙ from the optical SED.

Taking the projected rotational velocity v_esini = 145 ± 3 km s^-1 obtained from the DI inversions in this paper into account, we can refine the inclination angle to i = 46.5 ± 4.1° using the oblique rotator relation. This is close to i = 43 ± 7°, which was estimated by Trigilio et al. (2000), and is compatible with v_esini = 147 ± 2 km s^-1, i = 60° (Hatzes 1997) and with v_esini = 150 km s^-1, i = 60° (Hiesberger et al. 1995), given a typical uncertainty of 10–20° for the inclination derived by DI fits. On the other hand, Kuschnig et al. (1999) adopted v_esini = 160 km s^-1 and i = 30°, yielding a stellar radius of R = 3.3 R_⊙. This radius is substantially larger than other estimates and is unexpected for an unevolved main sequence star with the parameters of CU Vir.

Different sets of v_esini, i, and R discussed here correspond to different equatorial rotational velocities and, hence, a different impact of rapid rotation on the observed stellar properties. Using standard relations describing the shape of rapidly rotating stars and the corresponding gravity darkening (e.g., Collins 1963), we assessed the implication of different v_e. Figure 5 shows that the set of parameters used by Kuschnig et al. (1999) actually yields a noticeable rotational distortion of the stellar surface shape. In this case, the equatorial velocity amounts to 76 per cent of the critical value, and the pole-to-equator temperature difference reaches 3700 K. However, these effects were ignored by Kuschnig et al. (1999), which makes their DI inversions inconsistent with the basic stellar parameters adopted in their study. On the other hand, the smaller R and lower v_esini adopted by us and by Hatzes (1997) result in a nearly spherical stellar surface, v_e/v_crit ≤ 0.37, and a much lower temperature surface variation (500–800 K), which can be neglected in DI of CU Vir.

Table 2 summarises all parameters of CU Vir assumed or derived in our study.

3.2. Longitudinal magnetic field variation

The mean line of sight magnetic field, ⟨ B_z ⟩, was obtained from the normalised first moment of the Stokes V LSD profiles, as described by Kochukhov et al. (2010). For this analysis, we used the LSD profiles computed using all metal lines. These measurements, provided in Table 1, show a smooth, single-wave rotational variation with ⟨ B_z ⟩ changing between about 900 G and −700 G. The median error of our mean longitudinal field estimates is 60 G.

In Fig. 6, we compare our ⟨ B_z ⟩ measurements derived from the high-resolution Narval spectra with the photopolarimetric magnetic observations by Borra & Landstreet (1980), which were obtained more than 30 years ago from the hydrogen Balmer line wings. Despite using different diagnostic techniques and different spectral features, the two data sets show a reasonable agreement when phased together using the ephemeris with variable rotational period (Mikulášek et al. 2011). At the same time, our metal line ⟨ B_z ⟩ measurements are on average 2.5 times more precise and indicate a slightly higher amplitude of the longitudinal field variation.

Using the times of radio emission maxima at 1.4 GHz published by Trigilio et al. (2000, 2011) and Leto et al. (2006), we can establish that the peaks of radio emission occur at the rotational phases 0.34 ± 0.01 and 0.76 ± 0.01 for the variable-period ephemeris adopted in our study. These phases are not distributed symmetrically with respect to any of the two ⟨ B_z ⟩ extrema and are offset by −0.07 (first peak) and −0.04 (second peak) from the phases of the zero longitudinal magnetic field.

The longitudinal field curve of CU Vir exhibits an asymmetry between the two magnetic poles: the negative field extremum is somewhat sharper than the positive field maximum. This indicates a departure from a purely dipolar field configuration: for example, the presence of an offset dipolar field or a superposition of the dipolar and quadrupolar magnetic components. Adopting the latter field parameterisation and assuming that the quadrupole is axisymmetric and aligned with the dipole (e.g. Landstreet & Mathys 2000), we obtained B_d = 4.0 ± 0.2 kG, B_q = −7.6 ± 2.2 kG, and β = 79 ± 2° for a fixed inclination angle i = 46.5°. The resulting fit to the longitudinal field curve (solid curve in Fig. 6) is generally satisfactory, although not perfect, and yields . The remaining discrepancies are probably related to the effects of chemical spots and brightness variation across the stellar surface. A pure oblique dipolar model gives B_d = 3.8 ± 0.2 kG, β = 76 ± 2° and corresponds to a significantly worse fit (, dashed curve in Fig. 6) to the new ⟨ B_z ⟩ measurements. In any case, the value of β agrees well with the results of previous determinations (Hatzes 1997; Trigilio et al. 2000) using the ⟨ B_z ⟩ data of Borra & Landstreet (1980).

Notwithstanding the improvement in the description of the ⟨ B_z ⟩ variation brought by considering quadrupolar terms, the magnetic field geometry inferred from the longitudinal field curve appears to be misleading. As we show below, this field configuration does not reproduce the observed Stokes V profile behaviour and, therefore, cannot be accepted as a realistic representation of the magnetic field structure at the surface of CU Vir.

3.3. Magnetic field topology and chemical abundance distributions

The ZDI reconstruction of the magnetic field topology and abundance distributions was performed separately for the Si and Fe LSD profiles. In Figs. 7 and 8, we compare the observed Stokes I and V spectra with the final fit achieved by the magnetic inversion code. Results of magnetic inversion using the Si LSD profiles are illustrated in Fig. 9. This plot shows spherical projections of the field modulus (), as well as the horizontal () and radial (B_r) field components. The bottom row shows the vector magnetic field map, highlighting the regions of different field polarity.

We find that the magnetic field topology of CU Vir is dipolar-like in the sense that the distribution of the radial field component is dominated by two regions of opposite field polarity. However, other details of the field topology deviate significantly from what is expected for a simple oblique dipolar magnetic field. The maximum local field strength of 4.0 kG is found for the negative radial field region visible at the rotational phase ≈0.6. On the other hand, the positive longitudinal field extremum corresponds to an extended surface region with a relatively low field strength of 1–2 kG. In other words, CU Vir exhibits a large asymmetry between the negative and positive magnetic poles. The field topology is also clearly non-axisymmetric with the strongest positive and negative radial field regions separated by about 60° in longitude. According to our inversion results, it is the presence of both regions on the stellar disk at the rotational phase interval 0.3–0.5 that is responsible for a strong crossover signature in the observed LSD Stokes V profiles of CU Vir.

The outcome of the magnetic inversions using the Fe Stokes I and V LSD profiles is presented in Fig. 10 using the same format as the Si magnetic maps in Fig. 9. The Fe magnetic maps agree closely with the results of Si ZDI inversion. Magnetic maps derived from the mean Fe line exhibit the same field strength distribution asymmetry, the same non-axisymmetric character, and similar ranges in the maps of different magnetic field components. We emphasise that analysis of the two sets of LSD profiles was carried out completely independently, starting from zero magnetic field in both cases. The Fe profile inversion results support all aspects of the magnetic field topology found with the ZDI using Si LSD profiles. This confirms robustness of the new ZDI inversion procedure and rules out a possibility of significant random artefacts in the resulting magnetic maps.

As discussed above, ZDI with spherical harmonic parameterisation of the field topology allows one to quantify the contribution of different harmonic modes to the stellar magnetic field topology. In particular, we consider the relative energies obtained by integrating B² of individual harmonic components over the whole stellar surface. Figure 11 compares the energies of the poloidal and toroidal harmonic modes with ℓ = 1–5. The poloidal components dominate over the toroidal ones. The dipolar field comprises 52–60% of the total magnetic field energy. The quadrupolar (ℓ = 2) and octupolar (ℓ = 3) modes are non-negligible, yielding 23–30% of the total magnetic energy together. The contribution of modes with ℓ ≥ 4 is insignificant.

The penalty function applied to minimise the unwarranted contribution of the high-order harmonic modes ensures that our ZDI inversion code recovers the simplest possible magnetic field geometry compatible with observations. Nevertheless, to verify the necessity of including complex field components, we carried out additional inversions in which we severely limited the maximum angular degree ℓ of the harmonic expansion. In particular, we attempted to fit the observed Si LSD profiles with an arbitrary abundance map and a pure dipolar poloidal field topology. As demonstrated by Fig. 12, this model provides a good description of the Stokes I spectra but fails to reproduce the circular polarisation profiles of CU Vir. Theoretical Stokes V profiles cannot match observations in the rotational phase interval 0.3–0.5 when the Stokes V amplitude reaches a maximum. The optimal dipolar field parameters (B_d = 3.0 kG and β = 57°) suggested by this LSD profile fit are different from the corresponding parameters (B_d = 3.8 kG and β = 76°) inferred from the observed ⟨ B_z ⟩ curve.

Figure 12 also shows that the axisymmetric dipole plus quadrupole field model obtained by fitting the longitudinal field curve is inadequate for describing the high-resolution polarisation observations of CU Vir. The corresponding theoretical Stokes V profile variation roughly reproduces our observations in the phase interval 0.4–0.5 but also predicts strong and sharp polarisation signatures at phases 0.2–0.4 and 0.6–0.8. No such features are seen in the observed Stokes V spectra in either individual lines or LSD profiles.

Finally, Fig. 13 presents the surface distributions of Si and Fe abundances reconstructed from the LSD profiles of these elements. The Si abundance map spans the range from + 2.2 to −1.2 dex relative to the solar abundance. The Fe map shows a narrower range from + 1.4 to −0.2 dex. Distribution of both elements is dominated by a large spot of relative abundance depletion coinciding with the positive field region, which is characterised by a weaker predominantly radial field. On the other hand, an overabundance of both elements is correlated with the stronger negative field region and the neighbouring strong positive field region. As one can see from the horizontal field maps in Figs. 9 and 10, this is also the only part of the stellar surface where significant horizontal magnetic field is found. The main Fe overabundance zone is extended in longitude, tracing an arc close to the negative magnetic pole. A similar feature, but spread over a larger area, is visible in the Si abundance map.

4. Discussion

4.1. Comparison with previous DI studies

Previously, conventional abundance DI mapping was carried out for CU Vir by Hatzes (1997) and Kuschnig et al. (1999). The first study presented a Si equivalent width map that was derived from a single line; the second study obtained abundance maps for five different elements, including Si and Fe for which multiple lines were used. Both analyses gave a high weight to the strongest Si ii lines. Consequently, their results may be distorted by an interplay of horizontal and vertical abundance inhomogeneities and, therefore, are less representative than our chemical spot maps derived from the average line profiles. In addition, previous DI studies completely neglected magnetic field and ignored local modifications of the atmospheric structure.

Despite these differences, the main features in our Si and Fe abundance maps agree reasonably well with the distributions obtained by Kuschnig et al. (1999). Taking the 0.5 phase offset between the ephemeris adopted in the two studies into account, we both find an underabundance area around phase 0.0 with an Fe underabundance zone that trails slightly behind the one for Si. Similar to our study, the maps by Kuschnig et al. (1999) also exhibit arch-like features around phase 0.5–0.6 when the Si and Fe abundances reach their maximum. Kuschnig et al. reported a surface abundance variation from +2.5 to −1.5 for Si and from +1.3 to −0.7 for Fe. These ranges are somewhat larger compared to our study, which may be related to a lower inclination angle adopted by these authors and/or their neglect of the continuum intensity variation due to spots.

4.2. Verification of the Si abundance map

In this study, we attempted to derive chemical abundance distributions from LSD profiles, which represent an average of many spectral features with different characteristics. Since this approach is used here for the first time, it is of interest to check how the resulting chemical spot maps reproduce the Stokes I and V profiles of individual spectral lines. In the particular case of CU Vir, such an assessment is best carried out for Si, since there are many reasonably unblended lines of this element in the stellar spectrum. Using the results of magnetic and chemical abundance mapping with the Si LSD profiles (see Figs. 9 and 13), we computed about a dozen of Si ii spectral lines and compared their predicted profiles with observations. Figure 14 illustrates this comparison for five representative spectral features.

We find that the LSD Si abundance map reproduces weak and intermediate strength Si ii lines with the excitation potential of the lower level E_i ≈ 10 eV (e.g. Si ii 5041, 5056 Å) fairly well. However, a number of discrepancies is found for other lines. For example, the higher excitation Si ii lines, such as 4276 and 5268 Å (E_i = 13–15 eV), appear substantially stronger in observations than in our spectrum synthesis. We had to increase their oscillator strengths by 0.5–1.5 dex relative to the values given in VALD to match the observed profiles. On the other hand, very strong silicon lines such as Si ii 6347 and 6371 Å (the latter not shown in Fig. 14) exhibit a somewhat different profile shape, especially in the phase interval of 0.23–0.42, when a distinct flattening of the line cores is observed.

This line-to-line scatter is almost certainly caused by an inhomogeneous vertical distribution of Si in the atmosphere of CU Vir. If the Si abundance decreases with height, as found for magnetic CP stars of similar T_eff (Glagolevskii et al. 2005; Ryabchikova et al. 2006) and predicted by theoretical radiative diffusion models (LeBlanc et al. 2009; Alecian & Stift 2010), one might expect anomalously strong high-excitation lines (formed in deep atmospheric layers) and profile distortions in strong lower excitation features (sensitive to a large range of heights, including high atmospheric layers). This qualitative scenario describes reasonably well our observations of CU Vir. This star appears to be an interesting target for studying how the vertical stratification of Si changes across the stellar surface.

In conclusion, our LSD-based Si abundance map provides a satisfactory description of the intensity and circular polarisation spectra of typical individual Si ii. It does not, however, explain the profiles of certain features likely to be affected by chemical stratification.

4.3. Magnetosphere and radio emission of CU Vir

Our target shows a complex and interesting behaviour in the radio domain. Similar to some other magnetic chemically peculiar stars, CU Vir shows a smoothly varying, rotationally modulated gyrosynchrotron emission (Leone et al. 1996). In addition, unlike any other known radio-loud stellar source, it exhibits two 100 percent circularly polarised radio pulses for every rotation period (Trigilio et al. 2000). These radio pulses are attributed to the electron cyclotron maser (ECM) mechanism, which produces emission directed at a certain narrow range of angles to the magnetic field lines. This emission is believed to be produced in the stellar magnetosphere at a distance of about two stellar radii in the vicinity of one of the magnetic poles. Different modifications of the ECM model have been examined for CU Vir (see Lo et al. 2012), but it is still not entirely clear which one is the most appropriate.

Our work has provided several important constraints for the nature of the magnetosphere of CU Vir and its relation to the unique coherent radio emission of this star. First, our magnetic observations are essentially concurrent with the recent radio studies. This allows us to establish the phase relation between the longitudinal magnetic field curve and the times of arrival for the radio pulses (see Sect. 3.2) with a much higher certainty than in the previous studies, which had to rely on the ⟨ B_z ⟩ measurements by Borra & Landstreet (1980) that were acquired several decades ago. Along these lines, we report new precise phase offsets between the radio pulses at 1.4 GHz and the phases of the zero longitudinal magnetic field.

Furthermore, we found non-negligible deviations from a dipolar magnetic field topology for CU Vir. All previous radio emission models considered purely dipolar field configurations, occasionally speculating over the presence of non-dipolar field components as a means to explain the asymmetry of the radio pulses (e.g., Trigilio et al. 2000). Our ZDI maps characterised these non-dipolar field components and provided detailed information about the surface magnetic field topology that can be used for realistic modelling of the stellar magnetosphere and the radio emission processes. Although such analysis is beyond the scope of our study, it is already clear that our magnetic maps are generally compatible with the asymmetry of the beamed radio emission. Evidently, the negative field region is more compact, is geometrically much better defined, and has a stronger field compared to the positive field zone. Therefore, association of the radio pulses with only the negative surface field extremum comes as no surprise given our ZDI results.

Magnetospheric models that have been invoked invoked to explain the radio emission of CU Vir estimated that the magnetic field channels the mass loss out to the Alfvén radius of ~15 R (Trigilio et al. 2004). This is much larger than the Keplerian radius, thus providing conditions for the accumulation of material in co-rotating clouds or in warped disk (Townsend & Owocki 2005). For CU Vir, several studies suggested the presence of a torus of relatively cold, dense material surrounding the star at the magnetic equatorial plane (Trigilio et al. 2004; Leto et al. 2006). Lo et al. (2012) argued that refraction from this cold torus represents the most promising version of the ECM mechanism to explain the frequency dependence and shapes of the radio pulse profiles. The circumstellar material in the torus may produce characteristic spectroscopic variation in the hydrogen Balmer lines, which is similar to the behaviour commonly observed in hotter early-B magnetic stars (Townsend et al. 2005; Petit et al. 2013). We searched for such variations using the Narval spectra of CU Vir. The Hβ–δ lines exhibit smooth changes in their Stark wings due to a non-uniform He distribution over the stellar surface (Kuschnig et al. 1999; Shulyak et al. 2004). The maximum line width is found at phase ≈0.10; the minimum occurs at phase ≈0.65. In addition, Hα shows a subtle variation, primarily on the red side of the line core, which is similar to the spectroscopic changes described earlier by Shore et al. (2004). Relative to the mean profile, the line core alternates between a lack of absorption at phase ≈0.24 and an excess at phase ≈0.75. This variability in Hα is, again, a likely result of He abundance spots rather than cold, dense gas trapped in the stellar magnetosphere. Thus, our observations provide no direct evidence for the existence of circumstellar material in the vicinity of CU Vir.

4.4. Conclusions

In this study, we developed a new Zeeman Doppler imaging procedure, which combines the LSD multi-line technique with the detailed magnetic spectrum synthesis using realistic model atmospheres for the first time. This new magnetic inversion methodology was successfully applied to a set of high-quality circular polarisation observations of the well-known magnetic chemically peculiar star CU Vir. The main conclusion of our investigation can be summarised as follows:

• The new ZDI methodology enables a simultaneous reconstruction of the stellar magnetic field topology and a distribution of chemical abundance (or temperature for cool active stars) from the LSD Stokes profiles without making simplifying approximations that are typical of previous ZDI studies. In particular, the new ZDI method does not assume that a response of the LSD profiles to the variation of magnetic field, chemical abundance, or temperature is equivalent to that of a single spectral line.

• The revised fundamental parameters of CU Vir suggest that this is a young magnetic star observed from an intermediate inclination. It is a relatively rapid rotator, presumably due to its youth, but its equatorial rotational velocity is not high enough to significantly distort the spherical shape of the stellar surface.

• The magnetic field topology of CU Vir is generally dipolar-like but is characterised by a strong field modulus asymmetry between the two poles and by a major deviation from an axisymmetric configuration. This field geometry is mainly poloidal, and nearly all field energy is contained in the harmonic modes with the angular degree ℓ = 1–3.

• Our calculations prove that neither the purely dipolar magnetic field model nor a more complex dipole plus quadrupole field topology derived from the longitudinal field curve of CU Vir can describe the observed circular polarisation signatures in spectral lines.

• The chemical abundance maps of Si and Fe are dominated by a large underabundance region, where the field is positive, predominantly radial, and weak. In contrast, the spots of element overabundance occur at the stronger negative and positive field regions where a significant horizontal field component is also present.

• The magnetic and chemical abundance maps obtained from the Si LSD profiles match the variability of the Stokes I and V spectra in individual Si ii lines. However, observations show some line-to-line scatter in the Si profile shapes and intensities, which we tentatively attribute to a vertical stratification of silicon in the atmosphere of CU Vir.

• The magnetic field maps derived in our study explain the asymmetry of the pulsed radio emission of CU Vir and provide detailed constraints for future magnetospheric modelling.

Acknowledgments

O.K. is a Royal Swedish Academy of Sciences Research Fellow, supported by grants from the Knut and Alice Wallenberg Foundation, Swedish Research Council, and the Göran Gustafsson Foundation. T.L. acknowledges support from the FWF NFN project S11601-N16 “Pathways to Habitability: from Disks to Active Stars, Planets and Life”, and the related FWF NFN subproject S116 604-N16, as well as the financial contributions of the Austrian Agency for International Cooperation in Education and Research (CZ-10/2012). E.A., C.N., and the MiMeS collaboration acknowledge financial support from the Programme National de Physique Stellaire (PNPS) of INSU/CNRS. The computations presented in this paper were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under project snic2013-11-24.

References

All Tables

All Figures