Early science with SPIRou: near-infrared radial velocity and spectropolarimetry of the planet-hosting star HD 189733

SPIRou is the newest spectropolarimeter and high-precision velocimeter that has recently been installed at the Canada-France-Hawaii Telescope on Maunakea, Hawaii. It operates in the near-infrared and simultaneously covers the 0.98-2.35 {\mu}m domain at high spectral resolution. SPIRou is optimized for exoplanet search and characterization with the radial-velocity technique, and for polarization measurements in stellar lines and subsequent magnetic field studies. The host of the transiting hot Jupiter HD 189733 b has been observed during early science runs. We present the first near-infrared spectropolarimetric observations of the planet-hosting star as well as the stellar radial velocities as measured by SPIRou throughout the planetary orbit and two transit sequences. The planetary orbit and Rossiter-McLaughlin anomaly are both investigated and modeled. The orbital parameters and obliquity are all compatible with the values found in the optical. The obtained radial-velocity precision is compatible with about twice the photon-noise estimates for a K2 star under these conditions. The additional scatter around the orbit, of about 8 m/s, agrees with previous results that showed that the activity-induced scatter is the dominant factor. We analyzed the polarimetric signal, Zeeman broadening, and chromospheric activity tracers such as the 1083nm HeI and the 1282nm Pa\b{eta} lines to investigate stellar activity. First estimates of the average unsigned magnetic flux from the Zeeman broadening of the FeI lines give a magnetic flux of 290+-58 G, and the large-scale longitudinal field shows typical values of a few Gauss. These observations illustrate the potential of SPIRou for exoplanet characterization and magnetic and stellar activity studies.


Introduction
The study and characterization of exoplanet systems has greatly benefited from the instrumental developments in high-resolution optical spectroscopy since the pioneer work by Walker et al. (1995), Mayor & Queloz (1995), and Marcy & Butler (1992) and in particular, from the improvement of simultaneous wavelength calibrations (Baranne et al. 1996) in the advent of precise radialvelocity (RV) measurements. While recent years have seen very interesting developments in the optical with instruments suh as ESPRESSO (Pepe et al. 2014) and EXPRES (Petersburg et al. 2020), the strongest emphasis on precision-RV instrumentation Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. Based on observations obtained with SPIRou, an international project led by Institut de Recherche en Astrophysique et Planétologie, Toulouse, France. has been in the near-infrared (NIR) domain, with CARMENES-NIR (Becerril et al. 2017), GIARPS (Claudi et al. 2017), IRD (Kotani et al. 2018), HPF (Mahadevan et al. 2014), CRIRES+ (Dorn et al. 2016), and SPIRou (Donati et al. 2020). This new branch of RV NIR spectrographs has started a new era, in which the characterization of exoplanet systems widens toward a spectral domain full of riches: not only can different planet hosts be reached (young stars, and the coolest red dwarfs), but different atomic and molecular species are accessible to extend their characterization. In the extrasolar systems that can be observed both at optical and NIR wavelengths, the relative amplitude of stellar activity and planet signatures can, for instance, be used to distinguish between the two phenomena (Hébrard et al. 2014;Reiners et al. 2010;Zechmeister et al. 2018). Moreover, Zeeman broadening increases with wavelength, which offers new ways to explore the surface activity of host stars that plagues the radial velocity technique . In this era of new NIR spectrographs with precise RV measurements, SPIRou offers the Article number, page 1 of 18 arXiv:2008.05411v1 [astro-ph.EP] 12 Aug 2020 A&A proofs: manuscript no. paper advantage of including the full Y JHK coverage (0.98-2.35 µm) and a polarimeter to monitor the evolution of circular polarization in lines. This is an additional tool for stellar characterization (Hébrard et al. 2016).
We study early SPIRou science data of the well-known hot-Jupiter system HD 189733. Although SPIRou is optimized for stars cooler than K types, whiose wealth of information contained in the spectral domain is greater, it seemed interesting to revisit this system with a new instrument and study the spectroscopic content and behavior in the NIR domain. HD 189733 is a moderately active K2 dwarf star that is transited by a hot-Jupiter planet in a 2.218-day circular orbit (Bouchy et al. 2005). It has a distant M-dwarf companion (Bakos et al. 2006b). Since the planet was discovered and the first spectroscopic transit observations were reported, it has been extensively studied with highresolution optical spectrographs (e.g., Winn et al. 2006;Bakos et al. 2006a;Redfield et al. 2008;Triaud et al. 2009).
The stellar activity of the planet host HD 189733 has been observed with space-based photometry with MOST (Croll et al. 2007), unveiling a 3% flux modulation due to the rotation of contrasted features in the photosphere. Radial velocities obtained with SOPHIE in turn exhibit a stellar activity jitter amplitude of about 10 m/s rms when the planetary orbit is removed (Boisse et al. 2009).
In spectropolarimetry, the star was also extensively followed up with CFHT/ESPaDOnS and TBL/NARVAL for its strong magnetic activity. Features induced by star-planet interactions were likewise searched for Fares et al. 2010Fares et al. , 2017. The magnetic field of HD 189733 is characterized by a complex topology. The mean value of the total field is about 40 G, and there is a strong contribution of the toroidal component (Fares et al. 2017).
Recently, HD 189733 transits have been observed using CRIRES (Brogi et al. 2016;Brogi & Line 2019), GIANO (Brogi et al. 2018) and the NIR arm of CARMENES. These observations and previous ones in transmission spectroscopy led to the discovery of water, carbon monoxide, and helium in the extended atmosphere of the planet (Birkby et al. 2013;Salz et al. 2018;Alonso-Floriano et al. 2019).
In this paper, we present data on the orbital signal, on the Rossiter-McLaughling radial-velocity anomaly, on the circularly polarized Zeeman signatures in photospheric atomic lines, on the average small-scale surface magnetic flux, and on Paβ and HeI chromospheric lines of HD 189733 as seen by SPIRou. The transmission spectroscopy analysis will be presented in a forthcoming paper. In section 2 we provide a short description of the instrument, the observing strategy, the data collection, and the methods used in the data analysis. In section 3 we present our radial-velocity results. In section 4 we present our analyses of the stellar activity, of polarization features, and Zeeman broadening before we conclude in the last section.

SPIRou
SPIRou is the new NIR spectropolarimeter and high-precision velocimeter installed in 2018 at the Cassegrain focus of the 3.6 m Canada-France-Hawaii Telescope atop Maunakea (Donati et al. 2020). Its wide spectral range (0.98-2.35 µm), high resolving power (about 70 000), its temperature-controlled environment (Challita et al. 2018), and precise light injection devices (Parès et al. 2012) are key for obtaining time series of stellar spectra with extreme stability . The CFH telescope is guided with a system that includes a tip-tilt plate operating at a frequency of 50 Hz for the observations of HD 189733 and a NIR imaging camera sensitive to Y JH bands (Barrick et al. 2018). This system then allows the operator to position the star on the 1.29 circular aperture. The achromatic polarimeter located in the Cassegrain unit is composed of two ZnSe quarterwave rhombs and a Wollaston prism, splitting the beam into the two orthogonal states of the selected polarization (circular polarization was used in this study). The spectrograph is fed with three fluoride fibers: two science fibers collecting the light out of the polarimeter, and one calibration fiber for the spectral reference lamp (Micheau et al. 2018). A short octogonal fiber in the spectrograph contributes to a better scrambling to increase the RV precision. The SPIRou cryogenic spectrograph is operated in vacuum, at 73 K, and is thermally regulated at a submillikelvin level (Reshetov et al. 2012). The design of the echelle spectrograph is described in Thibault et al. (2012) and contains an R2 diffraction grating and a cross-disperser train consisting of three prisms. Finally, the detector used in SPIRou is a 15-micron science-grade H4RG from Teledyne Systems with up-the-ramp readout mode; the detector characteristics are presented in Artigau et al. (2018). The orders are tilted relative to the detector rows, providing a good numerical sampling of the resolution element.
SPIRou can be used in spectropolarimetric and spectroscopic mode. In both modes, the star light passes through the polarimeter that is located in the front end. In the spectroscopic mode, each exposure is acquired in a single rhomb position, while in the polarimetric mode, four sub-exposures with different rhomb configurations are necessary to retrieve each Stokes parameter, as usual in optical spectropolarimeters (Donati et al. 1997). Radial velocities of the star can be measured in both modes, and the same set of calibrations is used for the two SPIRou modes.
SPIRou is operated with an extensive calibration plan, including daytime and nighttime observations. Each afternoon a series of internal calibrations is run, including detector dark frames and instrumental background (fixed time of 600s), order localization and flat-field frames using a tungsten lamp, accurate wavelength solution using exposures from a stabilized Fabry-Perot (Perruchot et al. 2018) and uranium-neon hollow-cathod lamps. The Fabry-Perot is also used simultaneously to stellar observations in the reference channel, allowing a monitoring of the residual spectrograph shifts over a timescale of the night. Extensive tests have shown that the calibrated instrument internally delivers radial velocities at a precision level better than 0.5 m/s rms over days. Nighttime observations for SPIRou calibration plan include a sky background of 10 min, at least one nightly radial-velocity standard and a few bright hot stars to monitor the telluric absorption spectrum in various conditions (for a total 10 min of observing time, nightly). SPIRou has demonstrated its sensitivity, resolving power, and radial-velocity precision on cool nonactive stars and was accepted as a guest instrument at CFHT in January 2019. The on-sky performances of SPIRou are outlined in a forthcoming paper (Donati et al. 2020). SPIRou has been open to the CFHT community since February 2019.

Observations
HD 189733 was observed with SPIRou between July 2018 and June 2019 as part of the science verification program and of the SPIRou Legacy Survey (SLS) 1 . Fourty-eight and eighty-six spectra were obtained in polarimetric and in spectroscopic mode, respectively. The polarimetric mode has been used in circular polarisation at a temporal sampling of about one Stokes V sequence per night when SPIRou was operated at CFHT between July 28 and October 25, 2018, for a total of 12 visits. In addition to recovering the planetary orbit signature, this sparse data set allows us to probe the magnetic field of the stellar host, which is already known from optical measurements with ESPaDOnS and NARVAL at much denser time sampling Fares et al. 2010Fares et al. , 2017. It also gives the option of investigating how stellar activity in the NIR is modulated by the 12d stellar rotation period. The spectroscopic mode was used for the observation of two planetary transit sequences, on 21-22 September, 2018, and on 14-15 June, 2019. The first partial transit sequence was obtained during a commissioning night and the baseline before ingress is missing. The second transit sequence has a longer baseline with equal baselines before ingress and after egress as the in-transit duration. Table 1 lists the nights with observations of HD 189733 considered in this paper, and the average properties of the N spectra taken in each of these nights. For polarimetric measurements, 4 consecutive exposures were taken, while for each of the two transits 36 and 50 exposures were taken in spectroscopic mode. The signal-to-noise ratio (S/N) per 2.28-km/s pixel bin in the middle of the H band ranges from 110 to 280 for 245 s exposures, with a median value of 230. Conditions were photometric during both transit sequences, with an average seeing estimated on the guiding images of 0.77 and 0.80 .

Data processing
The spectra were reduced with the SPIRou data reduction system, APERO, using version 0.5.0. The nonlinearity of the detector was accounted for at the readout level, and the slope of individual pixels was calculated along the readout ramp within the linearity regime. The pipeline uses the optimal extraction method (Horne 1986) and the daily-updated calibration database to produce 2D spectra on which the radial-velocity measurement and (whenever relevant) the polarimetric calculation are performed. For RV measurements, the two science channels are extracted and treated together; for polarization, they are extracted individually and combined, as we describe in section 4.3. The Blaze function is estimated from the flat-field in each channel and removed from extracted spectra. The wavelength calibration is obtained from the combination of exposures using the UNe hollow-cathod lamp and the thermally controlled Fabry-Pérot etalon (FP). The measurement of the FP cavity length is anchored on the absolute line positions of the UNe lamp, as originally proposed by Bauer et al. (2015) and described in detail for SPIRou by Hobson et al. (2020). The offset of the wavelength solution is then estimated for each exposure from the calibration of the simultaneous FP observation and the accurate calculation of the Barycentric Earth Radial Velocity (BERV).
Prior to cross-correlating the extracted spectra with a numerical mask to obtain the radial velocity, it is critical to correct for the telluric contamination that plagues the NIR domain. In the SPIRou pipeline, we use the PCA-based correction technique that was first developed and tested on HARPS data, as described in Artigau et al. (2014), which was then implemented for SPIRou. Nightly observed telluric standards are collected into a library covering a widely sampled parameter space in atmospheric conditions, and this library is used to reconstruct with PCA the best-match telluric contamination spectrum of a given science observation. We used five principal components for telluric correction. Figure 1 shows one example spectrum of HD 189733 that was extracted with the SPIRou pipeline. In the top plot, the black dots show the full extracted spectrum, and the red spectrum was obtained after correction for telluric contamination, excluding the most opaque atmospheric bands. The blaze function was not removed at this stage. In the bottom plot, we show a small part of the H band, after the blaze function was removed. In black and red we show the extracted stellar spectrum before and after correction of the most appropriate PCA-derived telluric spectrum, which is shown in green.
After the extracted spectra were corrected from the telluric contribution, the nominal RV calculation offered by the AP-ERO pipeline was the cross-correlation function with a numerical mask, that was then fit by a Gaussian. The mean velocity of the Gaussian is the radial velocity. This is well adapted for slow rotators of K type such as HD 189733 (Pepe et al. 2002;Bouchy et al. 2005), but has been shown to prevent the recovery of RV information content for fast rotators (e.g., Galland et al. 2005;Yu et al. 2019) or M dwarfs in the optical (e.g., Astudillo-Defru et al. 2015).

Masks for cross-correlation
The telluric-corrected spectra can be used to create a high-S/N template spectrum free of telluric lines, which is the result of combining all the spectra, excluding areas that are affected by a telluric line at the time of observation. We used the template of HD 189733 to produce a weighted mask that contains the central wavelength of the lines and a weight corresponding to their depth. The lines were first identified as a local minimum, more than 4σ lower than the variation in the local continuum computed at a distance greater than the typical stellar line full-width at half-maximum (FWHM). This criterion rejected blended lines of similar depth. A Gaussian fit was then made around these local minima to determine the central position of the line and to reject features whose appearance was too different from that of a stellar line. See Pepe et al. (2002) for a full description of the method of correlation with a weighted binary mask.
One difference between optical and NIR data is the role of residual telluric lines in the construction of the line list. Telluric lines show relative and variable motion with respect to the stellar lines, therefore the mask needs to be free of regions where the residual telluric correction still contaminates the spectrum and adds spurious RV variations. We found that telluric lines deeper than 40% leave significant residuals. Stellar lines that overlap telluric features with a depth larger than 40% were then masked out for a BERV corresponding to the maximum range associated with HD 189733. For shallower telluric lines, only the excess photon noise due to telluric residuals was included in the line weight (ranging from 1.5 times the photon noise for a line of 10% depth to 7 times the photon noise for a line of 30% depth). This mask-construction method is similar to method 3 described in Figueira et al. (2016). The produced K2 mask finally contained 1924 atomic and molecular lines in the whole SPIRou domain. As an experiment, other partial masks were made from the original ones: a first set of three masks included only lines in the J, H, and K bands from the K2 mask.
We then used the out-of-transit observations of the nights of September 21-22, 2018, and June 14-15, 2019, which were obtained in spectroscopic mode, to produce an alternative stellar mask. The line-list filtering makes use of a template spectrum obtained from the mean of all out-of-transit spectra, where we fit every line in the mask using a multiple Gaussian model to account for blends. After performing a standard least-squares fit, we compared the fitted central wavelength of each line with the Table 1. SPIRou observations of HD 189733. The exposure time corresponds to one exposure. The S/N is per pixel and per exposure at 1.7 µm. The barycentric Julian date is shifted by -2450000. SPIRou has been used in both polarimetric and spectroscopic modes ('Pol' and 'Sp'). Zeropoint (ZP) indicates the mean CFHT/SKYPROBE extinction during the visit, when available, in magnitude (Cuillandre et al. 2016). The last two columns "RV phn" and "RV rms" show the average photon noise of the sequence and the scatter observed at a given epoch over the N data points when the planet signal is removed, respectively. mask values and removed the lines in which the difference was greater than half of the FWHM, using FWHM∼12 km/s as the reference value. We also removed lines whose fitted amplitudes were lower than twice the mean flux error within the spectral region around the line. Finally, we removed lines in which the fitted FWHM lay outside an acceptable range of widths (a conservative range of 0.12 to 24 km/s was used here). The final K2filtered mask contains 1405 lines.
A last set of masks was created from ATLAS9 stellar models (Castelli & Kurucz 2003) and the VALD database, specifically, the model corresponding to an effective temperature of 5000 K and logg of 5.0. Originally, this mask contained 6102 atomic lines in the SPIRou domain; when most opaque atmospheric bands at  and  were removed, 4121 lines remained. These constitute the Castelli5000 mask. From this line list, we then selected 1289 lines with the highest Landé factors (high Landé with a factor larger than 1.4 and a median factor of 1.57) and 1323 lines with the lowest Landé factors (low Landé with a factor smaller than 1.1 and a median factor of 0.91). The Landé factor of each line affects the line broadening because of the unsigned magnetic field in active regions of the stellar surface. By selecting the lines with the highest Landé factor, one amplifies the effect of the magnetic field on the line broadening, and potentially, on the radial-velocity jitter. One could thus expect the cross-correlation function (CCF) of the high Landé mask to have a larger width and smaller depth than the CCF obtained with the low Landé mask in the presence of a strong magnetic field. The largest number of lines, those with average Landé factors of about the median 1.25, remains in this selection to enhance any contrasted behavior.  Table 2 shows the number of lines, the average FWHM and spectral shifts obtained on the cross-correlation functions with these masks. The number of lines is given in two ways: the smaller number is the number of lines in the line list, and the larger number in parenthesis includes all duplicate occurrences of lines when consecutive orders overlap (hence the larger difference in the J-band mask, where the overlaps are wider).
It might be expected that the K-band lines are either of similar width or slightly wider than the J-band lines because the Zeeman-splitting behavior increases with wavelength squared. We observe a 10% wider H-band CCF compared to the J-band CCF. However, the K-band CCF in turn is 17% narrower than the H CCF. The K band includes a larger fraction of molecular lines, which have a narrower profile than atomic lines, and the average profiles might reflect this as the dominant factor. Because these CCFs are not built from similar line densities, this discrepancy may also be due to a lower precision of the K-band CCF because the FWHM has a scatter four times larger for the K band than the H band. In addition, the comparison between the low-Landé and high-Landé masks shows no systematic broadening for the second one, within a dispersion of 0.2 km/s. This requires further investigations using other stars with various activity levels. The comparison of CCF width between the theoretical mask (Castelli) and the empirical mask (K2) shows that the second is 1 km/s narrower than the first. This may indicate that some reference wavelengths in the theoretical line list are slightly misplaced.
Finally, the K2-J line list shows a spectral shift compared to the reference K2 mask. If this shift were further investigated, it might show that some lines are contaminated by telluric residuals or residual blended stellar lines.
The RV precision of each wavelength bin is proportional to a RV = √ N× < I >, where < I > is the average intensity of the lines and N is the number of lines in the bin. The bottom plot of Figure 2 shows that the distribution of a RV over the spectral domain is quite different between the theoretical and empirical lists: while it is similar in the two sets of masks, it is twice lower in the H and K bands in the empirical mask than in the theoretical mask. The low-and high-Landé masks have a similar wavelength distribution.

Dispersion and conditions
In this section, we compare photon noise RV estimates and observed RV scatter in this data set featuring the K star HD 189733.
The pipeline estimates the RV uncertainty of each individual exposure using the optimal weight method on the spectrum, as introduced by Connes (1985) and Bouchy et al. (2001). Be-cause this is performed on the extracted spectrum rather than the CCF, this RV uncertainty indicator does not depend on the mask and underestimates the real value because telluric lines are still included. Moreover, while the photon-noise estimator in the pipeline gives relevant values for most spectra, some spectra had erroneous values that did not scale with S/N. For these, we estimated the photon-noise uncertainty by using the scaling  observed in the majority of spectra. The scaling function is where S/N is the S/N per 2.28-km/s pixel at 1700 nm (note that this relation can only be used for SPIRou observations of a star whose spectral type and rotational velocity are similar to those of HD 189733). The chosen threshold for correction was ±10% of the expected value. After this correction, the average RV precision of the data set is 1.8 m/s. Individual nightly averages of the estimated photon-noise RV uncertainty are given in Table 1 in the penultimate column, RV phn.
As each circular-polarisation sequence consists of four individual spectra, it is possible to extract and derive the spectrum position of all subexposures as opposed to the position of the combined Stokes I spectrum, and derive an estimate of the RV scatter inside a sequence, although small-number statistics applies here. The last column, RV rms, in Table 1 shows the nightly RV dispersion values. The lowest dispersion values obtained in this data set and with version 0.5.0 of the data reduction system are 3-4 m/s, or about twice the average photon noise. This reaches higher values when the extinction varies or the S/N is low. Worse seeing conditions in turn do not lead to significantly higher dispersion values in this limited data set.

Orbit
All RVs as measured with the K2-filtered mask are listed in Table  6. The radial velocity measurements obtained with SPIRou were then phased according to the following ephemeris (Baluev et al. 2019): As shown in Figure 3, the recovery of the planetary orbit as known from the literature is secured from this sparse data set, although the time series is spread over three commissioning runs from July to October 2018, and various changes were made in the instrument between runs (both hardware and software). It shows that daily calibrations operate as expected, by deriving the correct velocity shift from spectral lamp observations. The hardware changes during the six-month commissioning period were not considered significant enough to justify adding a free offset between runs in this study, because there were only a few visits.
The rms of these O-C residuals is 10.8 m/s rms with the K2 cross-correlation mask, when combined RVs per polarimetric C. Moutou et al: HD 189733 seen by SPIRou , and with the 33.5 m/s peak-to-peak variations reported by Triaud et al. (2009) between several transit sequences (HARPS). These were obtained with optical spectrographs.
We obtain a dispersion of O-C residuals of: 41.1 m/s with the K2 mask limited J, 12.2 m/s in H, and 24.7 in K. The dispersion is much smaller in H within a polarimetric sequence and over the full time series.

Rossiter-McLaughlin effect
The two transit sequences obtained in spectroscopic mode can be used to measure the Rossiter-McLaughlin effect. We can compare this with previous results from optical measurements. The K2-filtered mask was used for this analysis. The improvement in RV precision with this mask is about 1 m/s compared to the K2 mask, that is, an improvement of about 20%, but more importantly it removes a residual systematic slope of about 0.07 km/s/day, which is probably due to residuals from telluric lines that slowly move over the stellar spectrum with Earth rotation during the sequences.

Classical Rossiter-McLaughlin fit
The Rossiter-McLaughlin effect was fit using the model developed by Ohta et al. (2005). This model derives an accurate analytic formula for the radial velocity anomaly and takes the stellar limb darkening into account. In order to determine the limbdarkening coefficient (u) for HD 189733 b in the infrared range, we used 3D stellar model atmospheres from Hayek et al. (2012). A simple linear limb-darkening law was used to derive an approximate value of u = 0.433 for the J band.
Five free parameters were used to simultaneously fit the two transits of HD 189733 b: the sky-projected obliquity λ, the projected rotation velocity v sin i , the transit epoch τ, and two systematic velocities γ 1 and γ 2 . All other transit and Keplerian parameters were fixed at the values reported in the literature (see Table 3).
We sampled the posterior distributions of λ, v sin i , and τ using the Markov chain Monte Carlo (MCMC) software emcee (Foreman-Mackey et al. 2013), assuming uniform priors. We used 50 walkers and 2000 steps of the MCMC and discarded the first 500 steps. Best-fit values are the medians of the distributions, with 1 σ uncertainties that are derived by taking limits at 34% on either side of the median, as listed in Table 3. Figure 4 shows the RV measurements of HD 189733 taken during the two transits of planet b on September 21, 2018, and June 14, 2019, along with the best-fit model. We found a sky-projected obliquity λ = −3.6 +1.5 −1.4 degrees and vsini = 3.29 +0.09 −0.09 kms −1 . When we fit this separately, the second sequence has a smaller offset on the λ value than the reference optical value, −1.9 ± 1.8 degrees, compared to when the two sequences were combined, probably because the first sequence is not complete. It starts slightly after ingress and has no baseline before the transit. These partial sequences understandably lead to increased systematics. These values obtained with one or two SPIRou transits agree within 2σ with the values obtained from data in the optical range: λ = −0.85 The dispersion of residuals after the subtraction of the model is 5.31 m/s for both transits, or 4.08 m/s and 6.04 m/s for the September 2018 and June 2019 transit, respectively. Out-oftransit RV data have a dispersion of 3.82 and 5.25 m/s for the first and second transit sequences (15 and 27 baseline data points), respectively, after the orbital slope is removed. This number likely represents the current instrumental scatter that was measured on early-science SPIRou data of a bright K star with this version of the pipeline. For comparison, when the full K2 mask is used or when the measurement is limited to the J, H, and K bands, the dispersion of the residuals for both transits are 5.9, 26, 9, and 17 m/s, respectively, when the two transits are fit together, with the residual slope mentioned earlier. The final fit parameters did not change significantly for different values of the limbdarkening coefficient u corresponding to each band (u Y = 0.467, u H = 0.343 and u K = 0.286). In Figure 5, we show the correlation diagrams for all five free parameters. The best-fit value of the transit epoch, τ, is consistent with the value of (τ 0 ) in Baluev et al. (2019) by propagating over different orbits using the orbital period from Table 2

Doppler tomography
We applied another technique, Doppler tomography, to the same data to determine the obliquity measurement for HD 189733 b and checked its consistency with the classical RM method. The 50 CCFs from the second transit were used for this analysis. We considered a model of stellar CCF that includes a limb-darkened rotation profile convolved with a Gaussian corresponding to the intrinsic photospheric line profile and instrumental broadening, following the approach of Collier Cameron et al. (2010). The Gaussian bump due to the planet occultation was also modeled in the stellar line profile, whose spectral location depends on the planet position in front of the stellar disk during the transit. We made a tomographic model that depends on the same parameters as the classical RM fit (Sect. 3.2.1), and added the local line profile width, s (nonrotating local CCF width) expressed in units of the projected stellar rotational velocity (Collier Cameron et al. 2010;Dalal et al. 2019). The free parameters we used to fit the Gaussian bump were λ, vsin i , γ and s. The τ was fixed to τ 0 Article number, page 7 of 18 A&A proofs: manuscript no. paper  This is the first time that two different techniques, a classical RM fit and Doppler tomography, were used to measure obliquity in the NIR range. They are both compatible with the results obtained from the optical data. The data from the second sequence modeled alone gave a result closer to the optical data, although with slightly larger error bars; this is explained by the incomplete transit sequence of the first epoch, although it was observed in better conditions than the second epoch.

Activity jitter
HD 189733 is known to have an average rotation period close to 12 days (Henry & Winn 2008) while a differential rotation of 0.11 rad/d has been measured by Zeeman Doppler imaging (Fares et al. 2010(Fares et al. , 2017. We used an average period of 12 days with the following ephemeris:  Figure 7 shows the behavior of the residual RVs from the planetary orbit subtraction as a function of the rotational phase of the parent star. As earlier, the error bars reflect the relative scatter of the four consecutive RV measurements that were made within a polarimetric series (last column of Table 1). The residual RVs are quite dispersed in this data set, which extends over three months (eight rotational cycles). The peak-to-peak variations are about 30 m/s. Considering the eight first visits, which were obtained within a dozen days, that is, during a single rotational cycle, a more organized pattern appears with a maximum near phase 0.1 and a minimum near phase 0.4 before another rise (points with a red square in Figure 7). This suggests that cool spots cluster around rotation phase 0.3 at the time of our first eight observations. The rapid evolution of the stellar photosphere is well known (Boisse et al. 2009;Lanza et al. 2011), and it is not surprising that the five last RVs do not seem in phase with the first RVs, as they are separated by 47 days from the first batch, and spread over another 34 days. The data set is not sufficient for a relevant jitter analysis or for attempting to correct for it, but the possible dependence of the residuals on rotational phase may indicate that at least part of the RV residuals is due to stellar activity.
The FWHM of the CCF averaged over each polarimetric sequences does not correlate with the RV residuals of the orbit (Pearson coefficient -0.08). It may differ from observations made in the optical where the FWHM seems to be a good indicator of the stellar activity when it is due to a spot or a plage . This difference might also indicate that in SPIRou data, the CCF broadening is dominated by Zeeman broadening rather than spot and plage temperature contrast, and hence is less prone to rotational modulation. Alternatively, residual systematic noise in the intensity profile might affect its shape. Even in the more frequently studied cases of optical CCFs, the origin of variations in the FWHM is not always clear and may be more sensitive to instrumental systematics, as shown earlier for SOPHIE data of HD 189733 (Dumusque 2014) and other cases observed with HARPS-N Mayo et al. 2019).

Paβ and HeI activity indices
In order to compare the RV jitter to an independent activity criterion, we searched the spectra for chromospheric lines that would be sensitive to stellar modulated phenomena. We selected the Paschen β line at 1282 nm and the HeI triplet at 1083 nm. The Paβ spectral region is shown in Figure 8 (top) together with the reconstructed spectrum of telluric contamination (a byproduct of the telluric correction), which has very little effect in this domain. Two methods were compared. In the first, a Gaussian fit was adjusted to the core of the line in all spectra as well as to the mean spectrum. An average equivalent width per polarimetric sequence was then calculated after the mean value was removed. In the second method, the mean spectrum was first removed from each spectrum and the residual in the stellar line was integrated over the line width. The error of each visit was calculated to be the scatter within the polarimetric sequence. The two methods give similar trends, and the second method has slightly less noise than the first and is considered in the analysis. The telluricsubtracted spectra are used.
As a sanity check, this method was used on other lines close to the HeI and Paβ regions. Three lines at 1084.680, 1282.520, and 1283.495 nm were then measured in the same way as HeI and Paβ. The first two are contaminated by telluric residuals at some values of the BERV, and the last line is not contaminated. This comparison allowed us to set a threshold on any significant variability, accounting for the variable noise in the spectra. The main conclusion is that there is no significant difference between the variable behavior of the HeI and Paβ lines with respect to the check lines.

Spectropolarimetry
The Stokes V Zeeman signatures for our SPIRou observations were calculated using the ratio method as described in Donati et al. (1997) and Bagnulo et al. (2009), where we used a sequence of four exposures consisting of two pairs of observations that each correspond to different retarder position. For each exposure in the sequence we extracted the two orthogonalpolarization science channels individually and then applied the blaze flux correction. Then a wavelength correction was applied to each exposure to compensate for the BERV. Finally, the flux values of each exposure were interpolated using a B-spline to match the wavelength sampling of the first exposure in the sequence. Then we calculated the degree of polarization for each spectral point using the ratio method. We measured the pseudocontinuum of the polarization spectrum using a moving median with a window width of about 1000 points. Then we fit a loworder polynomial to the measured continuum points to model the continuum in the whole spectrum, which was further subtracted from the data.
In order to obtain the circular polarization (Stokes V) in the stellar lines, we used the least-square deconvolution technique (LSD) (Donati et al. 1997) with the Castelli mask that contains the Landé factor of each line, that is, the sensitivity of each line to the stellar magnetic field. These factors are currently only available through physical experiment or theory and are lacking for most molecular lines in the NIR spectrum. The exclusion parameter we used to remove spectral domains that are contaminated by telluric residuals corresponds to telluric lines that are deeper than 5%. Between the Landé factor restriction and the cut for telluric residuals, only 1366 lines were used in the polarization profile. An example of Stokes V signatures is shown in Figure 9. It also shows the LSD profile of the null spectrum (N), obtained from a different combination of individual spectra and aimed at confirming the potential instrumental artifacts (Donati et al. 1997). Table 4 lists all polarimetric sequences and their characteristics: stellar rotational phase, maximum signal and standard deviation in the Stokes V profile, and longitudinal field B . The latter was obtained from the following equation: B = −2.14 × 10 11 vV(v)dv where the mean wavelength λ 0 is 1512 nm and the average Landé factor g eff is 1.25. The values for B span a few Gauss for the few profiles in which a signature is detected in the Stokes V profile, in accordance with the values reported in the optical with ESPaDOnS Fares et al. 2010Fares et al. , 2017. The median error on B is 1.7 Gauss, which is to be compared with a median error of 1 Gauss in optical literature values, while exposure times that were used with ESPaDOnS are two to three times longer. It is expected that with similar exposure times, the error on the longitudinal field would have been similar between ES-PaDOnS and SPIRou. Moreover, the number of available lines with a known Landé factor is larger in the wavelength domain of NARVAL (0.37-1.0 µm) than in SPIRou (0.98-2.35 µm) for a K dwarf star. The Zeeman sensitivity, however, is higher in the NIR than the optical, which more or less compensates the lower number of lines. This explains the similar error bars on B for comparable exposure times. Correcting for telluric contamination before calculating the LSD profiles does not significantly change the Stokes V profile (at a level less than 1σ), which is expected because the ratio method acts as a filter against any signatures that do not change in a short time or with the polarization state, including telluric features. As seen previously, the theoretical mask is not optimal for deriving the radial velocity at high precision because it does not include molecular lines, in contrast to the empirical mask. Further iterations between the theoretical line lists and observed spectra are required to improve the Stokes V profiles, for instance, by including strong molecular features with known Landé factors. Figure 9 (bottom) shows the variation in the longitudinal magnetic field with respect to the stellar rotation phase. As for previous indicators, the variation is smooth during the first rotational cycle, while additional visits acquired several cycles later show no clear behavior. Although the longitudinal field and the RV residuals show indications of rotationally modulated evolution, these two quantities do not appear to be directly correlated with each other in this data set, similarly as has been observed for other stars (Hébrard et al. 2016;Hussain et al. 2016). Furthermore, the data collection is not sufficient to try modeling the magnetic topology of the stellar surface with the few detected signatures spread over several months.

Zeeman broadening
Zeeman broadening provides an alternative diagnostic of the stellar magnetic field, because it infers a magnetic field strength from line widths in the Stokes I spectrum. This provides a measure of magnetic field strengths on small scales, which are invisible in Stokes V because regions with opposite sign cancel each other out (Morin et al. 2013). Because Zeeman broadening scales with wavelength squared while most other broadening processes scale with wavelength, the ability to detect this effect is much higher with SPIRou than with optical spectrographs.
The approach taken here was to directly model carefully chosen spectral lines with a wide range of effective Landé factors and minimum blending. We used the Fe i lines at 1534.3788, 1538.1960, 1561.1145, and 1564.8510 nm. The 1564.8510 nm line was used by Valenti et al. (1995) and Lavail et al. (2017). Table 5 lists the line characteristics used in this study. The choice for these specific lines is the result of a trade-off analysis using the following criteria: large (and, respectively, very small) Landé factors, large line depths, no blending, and long wavelengths. Atomic and molecular lines in the K band that are used   in M stars (e.g., Flores et al. 2019), for instance, are too weak in the spectrum of HD 189733.
Synthetic spectra including Zeeman splitting were calculated using the Zeeman spectrum synthesis code (Landstreet 1988;Table 5. Atomic data of the Fe i lines used in the spectral synthesis including Zeeman broadening. g eff is the effective Landé factor, log g f the oscillator strength, and log C 6 is the van der Waals damping coefficient. log g f for Fe i 1564.85 is from Valenti et al. (1995), other log g f values are from fits to the solar spectrum by Petit et al. (2020, in prep.), and other data are from Kurucz (2014 Wade et al. 2001;Folsom et al. 2016). This produces localthermal-equilibrium spectra in all four Stokes parameters, although here we focused on Stokes I. A grid of MARCS model atmospheres (Gustafsson et al. 2008) was used as input. While Zeeman accurately reproduces atomic lines in LTE for K stars (e.g., Folsom et al. 2016), it does not calculate molecular line spectra. We therefore avoided molecular lines. Atomic data for four Fe i lines were taken from the VALD database (Kupka et al. 1999;Ryabchikova et al. 2015), but the oscillator strengths were not sufficiently accurate for this analysis. For the Fe i 1564.8510 nm line we adopted the empirical log g f of Valenti et al. (1995). For the other three lines we fit the oscillator strengths using an observation of the solar spectrum. The observation was obtained with SPIRou in sunlight reflected off the moon. Synthetic spectra were calculated for a solar model atmosphere and oscillator strength log g f values were fit through a χ 2 minimization routine.
Spectral synthesis was performed using the stellar effective temperature (4875 K) and gravity (4.56) from interformetric observations (Boyajian et al. 2015). A vsin i value of 3.47 km/s was inferred from the equatorial rotation period and inclination of the rotation axis from Fares et al. (2017) and the stellar radius of Boyajian et al. (2015). We assumed a microturbulence of 0.9 km/s, which is typical for this spectral type (e.g., Doyle et al. 2013). These parameters were considered as fixed. Then the spectral fit derived the following parameters: macroturbulence, metallicity, a magnetic field strength, and a filling factor for this magnetic field. A radial-tangential macroturbulence profile was used. The magnetic field was assumed to have a uniform radial geometry, covering a fraction of the surface given by the filling factor. This is clearly simpler than the real magnetic geometry, but Zeeman broadening is relatively insensitive to magnetic geometry, and this is the simplest geometry consistent with the observations, and is consistent with models used in previous Zeeman broadening studies (e.g., Valenti et al. 1995).
The four telluric-subtracted spectra were coadded for the visit with the highest S/N (September 23) and were then generalized to all epochs. A plot of the fitted lines for September 23 is shown in Figure 10 together with a comparison of the best fit of the four lines with a null magnetic field. This modeling shows that a nonzero magnetic field is required to fit the wings of the higher Landé factor lines, while the model also fits the width of the zero Landé factor line. The red line shows the line profiles when spectroscopic parameters are fixed and without Zeeman broadening: all three lines with significant g appear to be deeper. Finally, the best-fit model with additional magnetic flux parameters is shown as the blue line and matches all line profiles.
The best-fit macroturbulence is 3.5±0.3 km/s and metallicity of 0.09±0.02 (in 1.3σ agreement with the literature value of -0.03±0.08 (Torres et al. 2008)). On average for all epochs, HD 189733 has a magnetic field of 1.9±0.2 kG that covers 15±2% of the star; the average magnetic flux B f is 290±58 G. Individual uncertainties are obtained from the covariance matrix, scaled by the square root of the reduced χ 2 . A few sources of errors are not included (on parameters that were considered as fixed: Teff, gravity, and microturbulence), therefore these uncertainties may be slightly underestimated. The last column in Table 4 lists the individual B f values corresponding to each polarimetric sequence. While there is a mild correlation between the longitudinal field and the field flux, as shown in Figure 11, the latter cannot be related to RV residuals.
This value of the magnetic field strength and flux for HD 189733 is intermediate between nonactive, slow-rotators M dwarfs such as Gl 411 or Gl 380, which have typical fields lower than 1 kG (Moutou et al. 2017;Flores et al. 2019), and fast-rotator M dwarfs whose field usually ranges from 2 to 5 kG (Shulyak et al. 2019). With a rotational period of 12 days, HD 189733 is similar to Gl 49, DS Leo, Gl 208 or CE Boo, whose magnetic field fluxes lie between 0.5 and 2 kG according to far-red measurements of FeH and Ti (see Figure 5 of Shulyak et al. (2019)). Comparisons of the filling factor of 15% we determined with the filling factor of other solar-type dwarfs agree well with data in the compilations by Cranmer & Saar (2011) and See et al. (2019), as this parameter is highly variable. For a Rossby number of 0.403 as estimated for HD 189733 (Vidotto et al. 2014), Figure 4 in See et al. (2019) shows a range of magnetic filling factors of 5 to 50%.

Discussion and conclusion
We observed the planet-hosting star HD 189733 with SPIRou during its science verification phase and early operations. The SPIRou data allowed the recovery of the orbital signature and spectropolarimetric features in agreement with previous work done at optical wavelengths.
From the flux in the spectra of HD 189733, we can estimate the photon-noise RV uncertainty in the total NIR spectrum to be about 1.8 m/s for an S/N of about 250, using the method described in Bouchy et al. (2001). A more comprehensive study of the photon-noise distribution in SPIRou spectra was not included in the objectives of this article, and will be the object of future studies. In particular, it would be interesting to study whether the photon noise or the instrumental or processing scatter (or which portion of each) is responsible for the observed difference between the relative J, H, and K observed RV dispersion values. The scatter in the K band RV data is currently twice larger than in the H band (Table 2), in line with the relative weight of these bands shown in Figure 2. However, the situation is less clear for the Y J bands: while their contribution to the RV precision is expected to be as large as (or larger than) the H band, the observed scatter is about three times larger. Moreover, the CCF calculated with the Y J lines alone shows a redshift compared to the CCF obtained with both the H and K bands. Whether this is due to a differential effect of stellar activity or to a stronger effect of uncorrected absorption or emission lines in the Y J bands remains to be investigated. Furthermore, a comparison with similar analyses based on RV content of spectra of various spectral domains and resolving power (Figueira et al. 2016;Cloutier et al. 2018;Reiners & Zechmeister 2020) would help in understanding the specific limitations of SPIRou data for this type of stars.
The observed RV scatter over a timescale of an hour has a median of 6 m/s and it is ∼10 m/s over a timescale of a few months. If the short-term scatter represents the combination of photon-noise contribution and instrumental or processing jitter, then the activity jitter of HD 189733 in the period July to October 2018 is about 8 m/s rms. There is an indication of rotation modulation in this additional scatter. Because HD 189733 is known to be a star with an active surface, a comparison with the RV jitter that is usually attributed to stellar activity from optical data cannot be definite at this stage: it would require contemporaneous observations. A residual activity jitter of about 8 m/s, however, is plausible, given the range of values from 9 to 15 m/s of the literature values obtained in the optical domain (see Section 3.1 and references therein). While observing the stellar surface in the NIR range might damp the jitter because the spot contrast is lower than in the optical (Reiners et al. 2010), it might also enhance it with larger Zeeman broadening ). A precise comparison would therefore clearly require simultaneous campaigns with NIR and optical spectrographs that would ideally be equivalently equipped with spectropolarimetry and high-precision velocimetry for a precise diagnostics.
Further ways to improve the instrumental or processing RV precision budget are being investigated. One way is to increase the data set with a larger BERV span, in order to improve the stellar template. Overall telluric correction will also become more robust with the use of additional principal components, as allowed by a larger data set. Another way is to use template matching (as in Astudillo-Defru et al. (2015), e.g.) when this betterquality template is available. With more adjusted telluric corrections, the line list can also be enlarged with a net improvement expected on the RV precision as well as on the retrieval of activity indices such as bisector span and FWHM of the CCF. A more accurate wavelength calibration is also being implemented in the SPIRou pipeline. Secondary effects, such as contamination from the calibration channel or persistence from previous observations, are finally being investigated. A main difference between NIR and optical RV data sets is indeed that spectra are not independent of each other. Not only is the telluric subtraction dependent on the relative velocities, but flux contamination may occur between channels or in time (detector persistence, and smearing or electronic crosstalk from the calibration channel). While these effects are not expected to dominate the error budget for HD 189733, they should nevertheless be estimated when Fig. 10. Observed (black line and dots) and modeled spectra of HD 189733 in four lines of Fe I with various Landé factors (see the g value in the line label). Three models are compared: a model spectrum with fixed spectroscopic parameters and a magnetic field fixed to zero (red line); a model with fitted spectroscopic parameters in which all lines contributing equally to the fit, even the broadened ones and no magnetic field (orange line), and a model in which all parameters are fit (blue line) (see text). The dashed line represents the spectrum that is not corrected for the telluric lines. these corrections will be implemented and validated by testing. The data set presented in this paper is particularly adapted to test further improvements of the pipeline, when available, because the achromatic signatures of the planet are well known.
For the observation of two transits of HD 189733 it was possible to model the Rossiter-McLaughlin signal in two different ways (direct modeling and Doppler tomography), which both gave results that agree well with the optical analysis from the literature. The Doppler tomography result from the second, complete, transit, gives a projected obliquity of −0.5 ± 1.3 degrees, which agrees very well with the modeling of four transit sequences observed with HARPS, which allowed measuring −0.4 ± 0.2 degrees (Cegla et al. 2016). The spectroscopic analysis of these transits in search for planet atmospheric signatures is ongoing and will be published in a forthcoming paper.
The fast-evolving activity of the stellar surface of HD 189733 has been studied and modeled multiple times in the past (e.g., Boisse et al. 2009;Lanza et al. 2011;Fares et al. 2017) because in many ways, it affects the measurement of the system parameters, including atmospheric composition and planetary environment conditions (Pont et al. 2007;Sing et al. 2011;Bourrier et al. 2020). The stellar magnetic field of HD 189733 was also analyzed in the SPIRou data set. From the intensity spectrum, we were able to fit the unsigned magnetic field contained in active regions, for the Zeeman broadening, they produce in NIR atomic lines of various Landé factors. This analysis gave a magnetic field of 1.9 kG distributed over 15% of the stellar surface, corresponding to B f value of 290±58 G. These val-