# First Observation of the Complete Rotation Period of the Ultraslowly Rotating Magnetic O Star HD 54879

C. Erba, C. Folsom, A. David-Uraz, G. Wade, S. Seadrow, S. Bellotti, L.

Fossati, V. Petit, M. Shultz

► To cite this version:

C. Erba, C. Folsom, A. David-Uraz, G. Wade, S. Seadrow, et al. First Observation of the Complete Rotation Period of the Ultraslowly Rotating Magnetic O Star HD 54879. The Astrophysical Journal, 2024, 977 (1), pp.84. 10.3847/1538-4357/ad865a. hal-04919007

HAL Id: hal-04919007

https://hal.science/hal-04919007v1

Submitted on 30 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# First Observation of the Complete Rotation Period of the Ultraslowly Rotating Magnetic O Star HD 54879

C. Erba<sup>1,10</sup>①, C. P. Folsom<sup>2,10</sup>①, A. David-Uraz<sup>3,4</sup>, G. A. Wade<sup>5</sup>, S. Seadrow<sup>6</sup>, S. Bellotti<sup>7,8</sup>①, L. Fossati<sup>9</sup>①, V. Petit<sup>6</sup>, and M. E. Shultz<sup>6</sup>①

Received 2024 January 17; revised 2024 October 11; accepted 2024 October 11; published 2024 December 4

## Abstract

HD 54879 is the most recently discovered magnetic O-type star. Previous studies ruled out a rotation period shorter than 7 yr, implying that HD 54879 is the second most slowly rotating known magnetic O-type star. We report new high-resolution spectropolarimetric measurements of HD 54879, which confirm that a full stellar rotation cycle has been observed. We derive a stellar rotation period from the longitudinal magnetic field measurements of $P = 2562_{-58}^{+63}$ days (about 7.02 yr). The radial velocity of HD 54879 has been stable over the last decade of observations. We explore equivalent widths and longitudinal magnetic fields calculated from lines of different elements, and conclude the atmosphere of HD 54879 is likely chemically homogeneous, with no strong evidence for chemical stratification or lateral abundance nonuniformities. We present the first detailed magnetic map of the star, with an average surface-magnetic-field strength of 2954 G, and a strength for the dipole component of 3939 G. There is a significant amount of magnetic energy in the quadrupole components of the field (23%). Thus, we find HD 54879 has a strong magnetic field with a significantly complex topology.

*Unified Astronomy Thesaurus concepts: Early-type stars (430); Stellar magnetic fields (1610); Massive stars (732); Spectropolarimetry (1973); O stars (1137)*

## 1. Introduction

The large-scale spectropolarimetric surveys of OB stars conducted over the last two decades (e.g., BOB, MiMeS; T. Morel et al. 2015; G. A. Wade et al. 2016) have detected surface magnetic fields in around ~7% of massive OB stars. These fields typically have dipolar topologies, are globally organized, and are strong (~1 kG) and stable over timescales of decades (V. Petit et al. 2013; J. H. Grunhut et al. 2017; M. E. Shultz et al. 2018, 2019). Magnetic O-type stars are rare; thus, the detailed characterization of their field strengths and geometries is an important contribution to a developing understanding of this unique population.

The surface magnetic fields of OB stars significantly impact the star's long-term evolution (see, e.g., Z. Keszthelyi 2023 for a recent review). V. Petit et al. (2017) and C. Georgy et al. (2017) used the mass-loss quenching formalism of A. ud-Doula & S. P. Owocki (2002) to demonstrate that magnetic massive stars retain much more mass during their main-sequence evolution than inferred from nonmagnetic evolutionary models. Similarly, the effect of magnetic braking, following A. ud-Doula et al. (2009), has been incorporated into massive-star evolutionary models, showing a rapid decline of the surface rotation

(G. Meynet et al. 2011; Z. Keszthelyi et al. 2019, 2020, 2021, 2022). An extreme example of magnetic braking is found in the O8f?p star HD 108, which has the longest rotation period of all known O stars (54 yr; Y. Nazé et al. 2001; M. Shultz & G. A. Wade 2017; G. Rauw et al. 2023). HD 108 thus stands in stark contrast to the other 12 well-established magnetic O-type stars in our galaxy whose rotation periods are typically on the order of several days to months<sup>11</sup> (V. Petit et al. 2013; J. H. Grunhut et al. 2017).

The topic of this study, HD 54879, is a main-sequence O9.7V star (A. Sota et al. 2011) that is distinguished by its strong magnetic field and very slow rotation. Spectropolarimetric measurements of the star were first obtained as a part of the B-fields in OB stars (BOB) survey (e.g., T. Morel et al. 2015), with the definite detection of the star's magnetic field initially announced by N. Castro et al. (2015) and L. Fossati et al. (2015). Both papers reported the observed negative extremum of the star's longitudinal magnetic field ($\langle B_z \rangle$) curve to be ~-600 G, implying a strong surface magnetic field of at least 2 kG. Additional spectropolarimetric monitoring by G. A. Wade et al. (2020) and S. P. Järvinen et al. (2022) showed a slow increase of the longitudinal field to an observed positive extremum of $\langle B_z \rangle \sim 74$ G, before once again passing through magnetic null.

To date, estimates of the rotation period of HD 54879 have been based on fits to the incomplete monitoring of the

<figure><img src="image_2.png" /></figure>

<footnote n="1"><sup>1</sup> Department of Physics and Astronomy, East Tennessee State University, Johnson City, TN 37614, USA; christi.erba@gmail.com</footnote>

<footnote n="2"><sup>2</sup> Tartu Observatory, University of Tartu, Observatoorium 1, Tõravere, 61602, Estonia</footnote>

<footnote n="3"><sup>3</sup> Department of Physics and Astronomy, Howard University, Washington, DC 20059, USA</footnote>

<footnote n="4"><sup>4</sup> Center for Research and Exploration in Space Science and Technology, and X-ray Astrophysics Laboratory, NASA/GSFC, Greenbelt, MD 20771, USA</footnote>

<footnote n="5"><sup>5</sup> Department of Physics and Space Science, Royal Military College of Canada, PO Box 17000, Kingston, ON K7K 7B4, Canada</footnote>

<footnote n="6"><sup>6</sup> Department of Physics and Astronomy, and Bartol Research Institute, University of Delaware, Newark, DE 19716, USA</footnote>

<footnote n="7"><sup>7</sup> Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="8"><sup>8</sup> Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, IRAP/UMR 5277, 14 avenue Edouard Belin, F-31400, Toulouse, France
<sup>9</sup> Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, A-8042 Graz, Austria</footnote>

<footnote n="10"><sup>10</sup> co-first authors.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>

<footnote n="11"><sup>11</sup> The exception to this is HD 191612 (O6f?p-O8fp), with a rotation period of about 538 days (N. R. Walborn et al. 2010; G. A. Wade et al. 2011).</footnote>
<table><caption>**Table 1**
Archival and Newly Obtained (Indicated by * in the Instrument Column) Spectropolarimetric Observations of HD 54879, from HARPSpol (H), ESPaDOnS (E), and Narval (N)</caption><thead><tr><th>HJD<br>-2455000</th><th>⟨B<sub>z</sub>⟩<br>(G)</th><th>⟨N<sub>z</sub>⟩<br>(G)</th><th>RV<br>(km s<sup>-1</sup>)</th><th>Peak<br>S/N</th><th>φ</th><th>Inst</th></tr></thead><tbody><tr><td>1770.4993</td><td>-631 ± 12</td><td>-17 ± 11</td><td>27.77 ± 0.03</td><td>270</td><td>0.00</td><td>H</td></tr><tr><td>1971.0658</td><td>-561 ± 12</td><td>-1 ± 10</td><td>27.64 ± 0.06</td><td>795</td><td>0.08</td><td>E</td></tr><tr><td>1971.1066</td><td>-568 ± 12</td><td>14 ± 10</td><td>27.61 ± 0.06</td><td>810</td><td>0.08</td><td>E</td></tr><tr><td>2030.5195</td><td>-542 ± 19</td><td>-4 ± 18</td><td>27.74 ± 0.06</td><td>459</td><td>0.10</td><td>N</td></tr><tr><td>2032.5133</td><td>-544 ± 18</td><td>-1 ± 17</td><td>27.69 ± 0.06</td><td>479</td><td>0.10</td><td>N</td></tr><tr><td>2033.5159</td><td>-531 ± 19</td><td>1 ± 17</td><td>27.76 ± 0.06</td><td>480</td><td>0.10</td><td>N</td></tr><tr><td>2092.5485</td><td>-504 ± 12</td><td>-24 ± 12</td><td>27.76 ± 0.03</td><td>233</td><td>0.12</td><td>H</td></tr><tr><td>2095.5057</td><td>-509 ± 19</td><td>15 ± 19</td><td>27.76 ± 0.03</td><td>156</td><td>0.13</td><td>H</td></tr><tr><td>2310.6113</td><td>-491 ± 37</td><td>31 ± 36</td><td>27.61 ± 0.06</td><td>246</td><td>0.21</td><td>N</td></tr><tr><td>2338.6680</td><td>-445 ± 38</td><td>-10 ± 37</td><td>27.60 ± 0.06</td><td>229</td><td>0.22</td><td>N</td></tr><tr><td>2357.6202</td><td>-459 ± 37</td><td>2 ± 36</td><td>27.67 ± 0.06</td><td>228</td><td>0.23</td><td>N</td></tr><tr><td>2374.5860</td><td>-444 ± 33</td><td>16 ± 33</td><td>27.69 ± 0.06</td><td>255</td><td>0.23</td><td>N</td></tr><tr><td>2409.5524</td><td>-492 ± 36</td><td>-19 ± 35</td><td>27.68 ± 0.06</td><td>260</td><td>0.25</td><td>N</td></tr><tr><td>2439.4167</td><td>-538 ± 49</td><td>-72 ± 49</td><td>27.56 ± 0.06</td><td>191</td><td>0.26</td><td>N</td></tr><tr><td>2736.9835</td><td>-410 ± 16</td><td>-3 ± 15</td><td>27.47 ± 0.06</td><td>525</td><td>0.38</td><td>E</td></tr><tr><td>2758.8779</td><td>-409 ± 17</td><td>6 ± 16</td><td>27.46 ± 0.06</td><td>518</td><td>0.38</td><td>E</td></tr><tr><td>2775.9878</td><td>-409 ± 15</td><td>-17 ± 14</td><td>27.49 ± 0.06</td><td>587</td><td>0.39</td><td>E</td></tr><tr><td>2880.7390</td><td>-348 ± 17</td><td>12 ± 17</td><td>27.61 ± 0.06</td><td>553</td><td>0.43</td><td>E</td></tr><tr><td>3008.1272</td><td>-224 ± 12</td><td>-7 ± 12</td><td>27.68 ± 0.06</td><td>644</td><td>0.48</td><td>E</td></tr><tr><td>3066.0351</td><td>-182 ± 12</td><td>-5 ± 12</td><td>27.55 ± 0.06</td><td>649</td><td>0.50</td><td>E</td></tr><tr><td>3128.9253</td><td>-125 ± 14</td><td>8 ± 14</td><td>27.38 ± 0.06</td><td>514</td><td>0.53</td><td>E</td></tr><tr><td>3557.8753</td><td>75 ± 15</td><td>-2 ± 15</td><td>27.58 ± 0.06</td><td>658</td><td>0.70</td><td>E</td></tr><tr><td>3744.1304</td><td>-24 ± 13</td><td>1 ± 12</td><td>27.70 ± 0.06</td><td>678</td><td>0.77</td><td>E</td></tr><tr><td>4269.8770</td><td>-637 ± 14</td><td>2 ± 12</td><td>27.54 ± 0.06</td><td>729</td><td>0.97</td><td>E*</td></tr><tr><td>4630.8714</td><td>-500 ± 11</td><td>-3 ± 9</td><td>27.53 ± 0.06</td><td>830</td><td>1.12</td><td>E*</td></tr><tr><td>5239.0541</td><td>-427 ± 11</td><td>0 ± 10</td><td>27.57 ± 0.05</td><td>779</td><td>1.35</td><td>E*</td></tr></tbody></table>

**Note.** The longitudinal magnetic field (⟨B_z⟩) and null spectrum (⟨N_z⟩) are calculated using profiles computed with a least-squares deconvolution method, with an integration range of +5 to +50 km s⁻¹. The phase φ is calculated using the ephemeris and 2562 days period described in Equation (1). The peak S/N is around 550 nm for all spectra.

longitudinal magnetic field curve. These indicated the rotation period was longer than five years (G. A. Wade et al. 2020). The most recent sinusoidal fit reported by S. P. Järvinen et al. (2022) suggested a stellar rotation period of 7.2 yr, placing HD 54879 as the second slowest rotator among the known magnetic O-star population.

HD 54879's spectra show clear evidence of strong Hα emission, interpreted by N. Castro et al. (2015) as likely originating from circumstellar material. T. Shenar et al. (2017) reported variability in the star's Hα profiles over both short-term (days) and long-term (months to years) timescales, with the former being attributed to stochastic processes in the stellar wind (e.g., J. O. Sundqvist et al. 2012; A. ud-Doula et al. 2013), and the latter to variability associated with stellar rotation. These conclusions were supported with additional monitoring of the star reported by S. Hubrig et al. (2020) and S. P. Järvinen et al. (2022).

In this paper, we report the most recent spectropolarimetric measurements of HD 54879. The data, described in Section 2, indicate the ⟨B_z⟩ curve has once again turned over, providing evidence that observations now span one complete rotation cycle. We perform a new magnetic analysis in Section 3, and we provide a refined measurement of the stellar rotation period. We examine the radial velocity and Hα variability of the star, and investigate claims of differences in the strength of the longitudinal magnetic field of the star when the field is measured from lines of individual elements. We also present the first magnetic map of HD 54879. Our findings are discussed in Section 4.

## 2. Observations

### 2.1. ESPaDOnS Spectropolarimetry

Three new sets of high-resolution (R ~ 68,000) spectropolarimetric sequences of HD 54879 were obtained using the ESPaDOnS echelle spectropolarimeter at the Canada-France-Hawaii Telescope (CFHT) on 2021 February 24, 2022 February 20, and 2023 October 21.<fnref n="12" /> A log of these observations is provided in Table 1. The ESPaDOnS instrument covers an approximate wavelength range of 3700–10500 Å across 40 spectral orders. Each spectropolarimetric sequence consists of four subexposures, and yields four Stokes I spectra, one Stokes V spectrum, and two null (N) spectra, which were extracted using the methods described by J.-F. Donati et al. (1997). A detailed description of the reduction and analysis of ESPaDOnS data is provided in G. A. Wade et al. (2016). Each subexposure in the 2021 and 2022 measurements had a duration of 867 s; similarly, each of the 2023 subexposures had a duration of 880 s. The peak signal-to-noise ratio (S/N) per 1.8 km s⁻¹ pixel near 550 nm of the recently obtained Stokes V spectra ranges between 729 and 830. Thus, the new data are of similar quality to the available archival observations from ESPaDOnS.

### 2.2. Archival Data Sets

In Table 1 and in our corresponding analysis, we include all available archival spectropolarimetric observations of

<footnote n="12">¹² Program codes 21AC15, 22AC28, and 23BC19.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top to bottom: the Stokes V, null, and Stokes I profiles extracted from the three most recent observations of HD 54879, corresponding to Julian dates (−2455000) 4269.8770 (blue), 4630.8714 (red), and 5239.0541 (black) from Table 1. The integration range of +5 to +50 km s⁻¹ is indicated by the gray dashed lines in each panel. The V spectra deviate significantly from zero while the null spectrum does not, indicating the detection of the star’s magnetic field in these data.</figcaption></figure>

HD 54879 that were obtained with HARPSpol, mounted on the 3.6 m telescope at the European Southern Observatory’s (ESO) La Silla Observatory (N. Piskunov et al. 2011); with the Narval spectropolarimeter, mounted on the Bernard Lyot Telescope (TBL) at the Pic du Midi Observatory (M. Aurière 2003); and with ESPaDOnS at CFHT. The original reduction and analysis of these data were presented by N. Castro et al. (2015), G. A. Wade et al. (2020) and S. P. Järvinen et al. (2022). We have performed a complete reanalysis of the archival data, in conjunction with the new observations, to provide as consistent and homogeneous a set of results as possible. Our extracted longitudinal magnetic field and null profile (⟨N<sub>z</sub>⟩) values are consistent within uncertainties to the data reported by G. A. Wade et al. (2020), and are presented in Table 1.

# 3. Magnetic and Spectroscopic Analysis

## 3.1. Longitudinal Magnetic Field

The first step of the spectroscopic analysis was to normalize the spectra to their continuum. This was done by fitting low-degree polynomials through carefully selected continuum points, and spectral orders were normalized individually. We started from unnormalized extracted spectra for both the new and archival observations, to provide a consistent normalization for all data.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Stokes I (left column) and Stokes V (right column) profiles of the most recent spectropolarimetric measurements of HD 54879 (red) compared to archival data of the star at approximately the same rotation phase φ (black), but one cycle previous to the new data. The profiles are consistent with each other.</figcaption></figure>

We used least-squares deconvolution (LSD; J.-F. Donati et al. 1997) to calculate pseudoaverage line profiles for the intensity, Stokes V, and null spectra. This technique effectively boosts the S/N of the Zeeman signatures present in the Stokes V profiles, to enable the detection of stellar magnetic fields. Here we use a new open-source implementation of the LSD algorithm, LSDpy, contained within the SpecpolFlow software package (C. P. Folsom et al. 2024). This follows the description of O. Kochukhov et al. (2010), and provides results fully consistent with their iLSD code. We used the custom line mask of G. A. Wade et al. (2020), based on data from the VALD3 (N. E. Piskunov et al. 1995; T. A. Ryabchikova et al. 1997; F. Kupka et al. 1999; F. G. Kupka et al. 2000; T. Ryabchikova et al. 2015). This mask includes only metal lines, excludes He lines, and uses lines with stable Stokes I profiles to minimize any wind contribution. The observed He lines do not meet the self-similarity hypothesis of LSD, and display a variety of shapes due to varying amounts of Stark broadening, thus they were rejected. We verified that this mask is well optimized for the spectral lines present in observations from all three instruments, and used it consistently for all the spectra. We used a normalizing wavelength of 500 nm, effective Landé factor of 1.2, and depth of 0.1 for the LSD calculation. Every LSD Stokes V profile produced a definite detection,<fnref n="13" /> with the corresponding null profiles producing nondetections.

The Stokes I, V, and null LSD profiles from the most recently obtained ESPaDOnS data are shown in Figure 1. The Stokes V profiles show typical Zeeman signatures with a large amplitude (corresponding to the ⟨B<sub>z</sub>⟩ measurements in Table 1). In Figure 2, we compare these profiles to archival observations obtained at similar rotation phases (using the ephemeris found below). The profiles have nearly identical morphologies, amplitudes, and polarities, supporting the stability of the magnetic field and our derived period.

<footnote n="13">¹³ The false alarm probabilities (FAP) were calculated using the following criteria: definite detection, FAP < 10⁻⁵; nondetection, FAP > 10⁻³ (J.-F. Donati et al. 1997).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Upper panel: the phased time series of $\langle B_z \rangle$ measurements. The lower x-axis references the date of the observations, with the archival data sets in black dots, and the data from the 2021–2023 observing epochs in red. The timestamps are phased according to an ephemeris with $P = 2562$ days and $t_0 = 2456774$ (HJD), reported along the top axis. The dashed and solid black lines show the first- and second-order fits to the $\langle B_z \rangle$ data, respectively. Middle panel: the phased time series of $H\alpha$ equivalent width measurements, with colors as above. Overplotted is the best fit using a two-component harmonic function with the same period and ephemeris as the upper panel. Lower panel: radial velocity measurements as a function of HJD, with colors as above.</figcaption></figure>

The longitudinal magnetic field was measured from the LSD profiles using a first-moment technique (e.g., G. A. Wade et al. 2000) with an integration range of +5 to +50 km s⁻¹, and the normalizing wavelength and Landé factor of the LSD profiles. This is consistent with the procedure from G. A. Wade et al. (2020). We report the $\langle B_z \rangle$ and $\langle N_z \rangle$ (the $\langle B_z \rangle$ calculation performed on the null profile) for all observations in Table 1. Measurements from the three new polarimetric sequences indicate a dramatic decrease in the magnitude of the longitudinal magnetic field (from −24 G to −637 G), with a subsequent, smaller increase (to −500 G and −427 G). These values are close to the initial measurements of the longitudinal magnetic field strength from N. Castro et al. (2015) and G. A. Wade et al. (2020), strongly suggesting the magnetic field has now been observed over one complete rotational cycle.

Generally, the $\langle B_z \rangle$ and $\langle N_z \rangle$ values we find here are in good agreement with G. A. Wade et al. (2020). Comparing results for the same observations, they mostly differ by <1σ of the joint error and always by <2σ. Our results are in reasonable agreement with N. Castro et al. (2015), who report $\langle B_z \rangle$ for the first HARPSpol observation calculated with three different methods. Compared with their “Bonn” LSD analysis (the most similar to our procedure) our result is stronger by −39 G with a

joint error bar of 14 G, which is an acceptable agreement given differences in the line masks used. Our results consistently disagree by >3σ with the metal-line analysis of S. Hubrig et al. (2020), who reanalyzed the ESPaDOnS data of G. A. Wade et al. (2020) but found differing $\langle B_z \rangle$ values. S. Hubrig et al. (2020) present results for two observations obtained on the same night (JD 2456971) with $\langle B_z \rangle$ values differing by ~5σ (81 G), which suggests they substantially underestimated the real uncertainty in their results. Our results disagree with S. P. Järvinen et al. (2021) by >3σ for about half of the observations we have in common, based on their LSD analysis using all metal lines (their most similar results). S. P. Järvinen et al. (2021) similarly disagree with G. A. Wade et al. (2020). They partially agree with S. Hubrig et al. (2020), with six of nine of the measurements in common differing by <2σ. S. P. Järvinen et al. (2021) find large differences in $\langle B_z \rangle$ over short periods of time (up to 196 G in 39 days), which are not seen in other studies of the same data (and not seen in some of their alternative line masks). This suggests that there are significant systematic errors or underestimated uncertainties in some of their values.

The longitudinal magnetic field curve of HD 54879, presented in the upper panel of Figure 3, shows significant
<figure><img src="image_6.png" /><figcaption>Figure 4. Periodograms from fitting the ⟨Bz⟩ data using first- (red dashed line) or second-order (black line) Fourier series. In the inset frame, dotted lines represent χ²_ν contours for 1, 2, and 3σ confidence levels, about the best period for the second-order Fourier series fit.</figcaption></figure>

departures from what is expected from a dipolar topology (see also Section 3.5 below). To date, topologically complex magnetic fields with strong departures from a dipole have only been reported in one other O-type star (NGC 1624-2; A. David-Uraz et al. 2021; S. P. Järvinen et al. 2021), although complex geometries have also been found in a few early B-type stars (e.g., J.-F. Donati et al. 2006; O. Kochukhov et al. 2011). In this figure, we also show a first-order sinusoidal fit to the data (dashed line), which would be expected for a purely dipolar field. Also included is a second-order harmonic fit to the data (solid line), which is equivalent to including both dipolar and quadrupolar components.

Since the ⟨Bz⟩ measurements now cover a full rotation cycle, we can perform a periodogram analysis. We proceeded by fitting Fourier functions to the data for a grid of fixed periods, and recording the resulting reduced chi-squared (χ²_ν) for each period. For a first-order function this is equivalent to a Lomb-Scargle periodogram, except in units of χ²_ν rather than power, but this approach can be generalized for more complex magnetic fields by using a second- or third-order Fourier series. Periodograms for both the first- and second-order Fourier series fits are shown in Figure 4. The first-order fit produced a χ²_ν minimum near our adopted period, but with χ²_ν > 30, reflecting a strong departure from a dipolar magnetic field. The second-order fit reaches χ²_ν = 1.37, reflecting an adequate fit to the observations. Using the second-order fit, and placing the negative magnetic extremum at phase zero, we adopt the ephemeris

$$HJD = 2456774 + 2562_{-58}^{+63}E. \quad (1)$$

Thus the current best-fit rotation period is P = 2562<sub>-58</sub><sup>+63</sup> days, or about 7.02 yr.

The departure of the ⟨Bz⟩ curve from that expected for a pure dipole is highly statistically significant, which indicates a more complex magnetic topology. The success of the second-order fit to the ⟨Bz⟩ curve suggests that there is a strong quadrupolar component to the magnetic field. As noted above, such complex fields have been observed around cooler stars, such as τ Sco (B0.2V; J.-F. Donati et al. 2006), α² CVn (B9p;

J. Silvester et al. 2014), and HD 37776 (B2V; O. Kochukhov et al. 2011), but had not yet been so prominently observed in the ⟨Bz⟩ curve of an O-type star.

## 3.2. Radial Velocity Stability

With the LSD profiles calculated in Section 3.1, we calculated precise radial velocities for HD 54879. The LSD calculations used a mask with most available metal lines in the spectra. Radial velocities were determined by fitting a Gaussian to the Stokes I profiles using χ² minimization, and taking the centeroid. Uncertainties were taken from the square root of the diagonal of the covariance matrix, and thus only incorporate the pixel uncertainties in the LSD Stokes I profiles.

A plot of the radial velocity data as a function of time is in Figure 3 (lower panel), and this reveals no periodicity or significant trends with time. The radial velocities are highly stable within observations from one instrument, and are in good agreement between instruments. The full average is 27.62 km s⁻¹ with a standard deviation of 0.10 km s⁻¹. For just the 14 ESPaDOnS observations, the average is 27.56 km s⁻¹, with a mean uncertainty of 0.06 km s⁻¹, and standard deviation 0.08 km s⁻¹. For the nine Narval observations, the average is 27.67 km s⁻¹, with mean uncertainty and standard deviation of 0.06 km s⁻¹. For the HARPSpol data, the average is 27.76 km s⁻¹ with a mean uncertainty of 0.03 km s⁻¹.

We thus conclude that the radial velocity of HD 54879 has been stable to better than 0.1 km s⁻¹ for the last 9.5 yr. This definitively rules out proposals that this is an SB1 binary (e.g., S. Hubrig et al. 2019), unless the period is very long (>>10 yr). This stability also argues against most other kinds of metal-line variability, such as chemical spots, although we investigate that in more detail below.

## 3.3. Measuring the Magnetic Field from Lines of Individual Elements

In their recent paper, S. P. Järvinen et al. (2022) derive measurements of the longitudinal magnetic field strength of HD 54879 using single-element LSD line masks. They report that the ⟨Bz⟩ measurements obtained from LSD profiles
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Measurements of the longitudinal magnetic field from LSD profiles of He, O, and Si lines. The black curve is the best-fit second-order sine curve from Section 3.</figcaption></figure>

calculated using masks consisting individually of O, Si, and He lines are discrepant, while the LSD Stokes *I* profiles do not show variability. The authors thus conclude that these results are reflective of “different formation depths, with the He lines formed much higher in the stellar atmosphere compared to the silicon and the oxygen lines, and nonlocal thermodynamic equilibrium (NLTE) effects.” This sentence may be implying that O, Si, and He have vertically nonuniform abundance distributions, i.e., that they are stratified. However, later in their article the authors appear to claim a nonuniform *surface* distribution of O, Si, and He: “Thus, we conclude that the elements O, Si, and He are not uniformly distributed over the stellar surface of HD 54879...” This latter sentence appears to imply that O, Si and He have nonuniform lateral distributions, possibly in addition to being vertically stratified.

The different LSD masks associated with each of the elements they analyzed are comprised of different collections of lines of different depths, Landé factors, and widths, and certainly correspond to somewhat different measurement systems. Without calibration it is not surprising that the LSD profiles and longitudinal fields extracted using these different masks would be somewhat different. Nevertheless, we carried out an independent analysis of the observations using single-element line masks. We used the iLSD code of O. Kochukhov et al. (2010) in multimask mode. (It is not clear which LSD code S. P. Järvinen et al. 2022 employed.) Starting from the line mask of G. A. Wade et al. (2020), we created three new line masks containing only O lines, Si lines, and He lines, along with their associated antimasks (i.e., the original mask excluding the lines of each of those elements). We then extracted LSD profiles for each element, taking into account the blending by other elements using the antimasks. Finally, we measured the longitudinal magnetic field from each profile. Those measurements are shown in Figure 5.

As reported by S. P. Järvinen et al. (2022), there are some differences in the amplitudes of the ⟨*B*<sub>z</sub>⟩ variations for individual elements. However, the various ⟨*B*<sub>z</sub>⟩ curves have the same essential shapes, i.e., they are not strongly discrepant like those of chemically peculiar B-type stars such as HD 184927 (I. Yakunin et al. 2015). Moreover, we observe no significant variations of line profiles or equivalent widths (EWs) of these elements, either in LSD profiles or individual spectral lines. In Figure 6 we examine EWs of lines chosen to be free from blending, strong enough to be clearly seen, but not

<figure><img src="image_8.png" /><figcaption>**Figure 6.** EWs from individual lines of He, O, and Si. No significant variability is evident.</figcaption></figure>

so strong as to have large saturation effects, extreme NLTE effects, or a strong influence from the magnetosphere. These lines of evidence all point to small differences in the measurement systems corresponding to the various line masks, leading to small differences in the inferred longitudinal fields.<fnref n="14" />

In their papers, both N. Castro et al. (2015) and G. A. Wade et al. (2020) show that the spectrum of HD 54879 is remarkably stable, especially for metal lines. This stability is contrasted with the variability of the Hα profile (see Section 3.4). These observations seem to indicate that HD 54879’s dynamical magnetosphere does not demonstrably contribute in any significant way to the metal lines. Presumably, if there were a significant contribution from the dynamical magnetosphere, this would result in variability that would be visible in the EW measurements, which as Figure 6 illustrates, is not observed. Moreover, the high-resolution HARPS spectrum presented in N. Castro et al. (2015) is in good agreement with a FASTWIND synthetic spectrum corresponding to the inferred parameters of HD 54879, further supporting the assertion that most spectral lines are not significantly perturbed by the presence of the magnetosphere. We thus conclude that there is no compelling evidence for chemical stratification or lateral abundance nonuniformities in the atmosphere of HD 54879.

<footnote n="14">¹⁴ ⟨*B*<sub>z</sub>⟩ may be more sensitive to differences in the lines used for HD 54879 than for most stars. The star has a strong magnetic field, with substantial cancellation in the line-of-sight component, which causes detectable Zeeman broadening (see Section 3.5). Thus differences between lines in Landé factor and splitting pattern may be more important here, especially for masks with few lines.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Hα line from spectra taken with HARPS, NARVAL, and ESPaDOnS (see Table 1). The archival data are plotted as black solid lines, while the new data sets are indicated in red. The continuum region used for local normalization are highlighted with green vertical shading. The other vertical lines correspond to the integration bounds that were tested for the calculation of EWs, i.e., ±400 km s⁻¹ (solid black), ±500 km s⁻¹ (dashed black), ±600 km s⁻¹ (dashed–dotted black), ±700 km s⁻¹ (dotted black), and ±800 km s⁻¹ (dashed blue). We can see from this figure that HD 54879 exhibits significant Hα emission at all rotational phases that were observed.</figcaption></figure>

### 3.4. Hα Emission

HD 54879 exhibits atypically strong Hα emission for its spectral type (T. S. Boyajian et al. 2007; N. Castro et al. 2015), with variability on both short (days) and long (years) timescales (T. Shenar et al. 2017). Mass-loss rate estimates derived from this Hα emission are comparatively too large (by 1 order of magnitude; N. Castro et al. 2015), and so are considered unreliable.<fnref n="15" /> N. Castro et al. (2015) suggested that HD 54879's Hα emission is likely magnetospheric in origin. For a star with a circumstellar magnetosphere, outflowing line-driven stellar winds can become trapped by the magnetic field and are forced to corotate with the star. Typically, in slowly rotating stars, rotation is not fast enough to provide sufficient centrifugal support for the plasma to overcome gravity. The confined wind material will fall back to the stellar surface under the influence of gravity on a dynamical timescale, forming a “Dynamical Magnetosphere” (DM; J. O. Sundqvist et al. 2012). In the spectra of stars hosting DMs, Hα emission has been observed to have a single broad peak close to line center (e.g., V. Petit et al. 2013). The multipeaked Hα emission from HD 54879, illustrated in Figure 7, suggests its circumstellar environment is shaped by a complex magnetic field.

To explore the variation of the Hα profile and how it relates to the rotation period derived from $\langle B_z \rangle$ measurements, we compute EW measurements using the spectra listed in Table 1.

First, we refine the normalization that has already been applied to the spectra for the magnetometric analysis by defining continuum regions on either side of the Hα line (653.9–654.1 nm blueward, and 658.85–659.05 nm redward), and renormalizing the Hα line to a linear fit to those averaged regions. The resulting lines are shown overplotted in Figure 7.

We then conducted some tests to determine the best integration range over which to compute our EW measurements. We tried integration bounds at ±400, 500, 600, 700, and 800 km s⁻¹ with respect to the Hα line center (accounting for the systemic velocity of 27.62 km s⁻¹; see Section 3.2). As can be seen in Figure 7, all magnetospheric variability appears to be contained with the bounds at ±400 km s⁻¹, with larger bounds including undesirable photospheric lines overlapping with the wings of the Hα line profile. EWs computed with these various bounds followed a similar trend from one spectrum to another, with uncertainties increasing with larger bounds. This indicates that indeed, the smallest bounds contain the bulk of the magnetospheric variability, and therefore they were used for the rest of the analysis (even though they do not span the full Hα line profile).

Figure 3 (middle panel) shows the clearly coherent variation of the Hα EWs with time. The measurements are fit with a second-order harmonic curve (black line) with a period corresponding to that derived from the $\langle B_z \rangle$ analysis. Most of the EW variability is consistent with this sinusoidal fit, modulated at the rotation period. However, we also performed an empirical harmonic fit to the Hα data with a partly unconstrained period (it was allowed to vary between 2400 and 2900 days). We found that a two-component function is

<footnote n="15">¹⁵ Mass-loss rate estimates for magnetic stars derived with spherically symmetric (nonmagnetic) models will be inherently unreliable, since the magnetic channeling of the wind creates a fundamentally aspherical environment (see, e.g., J. O. Sundqvist et al. 2012; W. L. F. Marcolino et al. 2013; C. Erba et al. 2021).</footnote>
significantly preferred over a single sinusoid (with a reduced $\chi^2$ of 1.031 rather than 1.977). This method (with the two-component function) yielded a best-fit period of 2827 days (~7.75 yr), about 10% longer than the stellar rotation period estimated from the best fit to the longitudinal magnetic field data. Minor discrepancies in estimates of the stellar rotation period derived from fitting $H\alpha$ EWs and $\langle B_z \rangle$ measurements are not unexpected, since the $H\alpha$ emission can be modulated by processes other than the stellar rotation period (see, e.g., $\xi^1$ CMa, B0.5 IV; C. Erba et al. 2021). Future work will be needed to explore and interpret the modulation of the $H\alpha$ emission in HD 54879's circumstellar environment.

Shorter-term variability, as observed by T. Shenar et al. (2017) and also evidenced in our measurements (see Figure 3), could conceivably cause systematics in the determination of the period using EW measurements, although it appears unlikely to fully explain the discrepancy between the periods derived from both methods. Further investigation, both observationally (more coverage) and theoretically (e.g., modeling $H\alpha$ emission using analytical prescriptions; S. P. Owocki et al. 2016) will be needed to resolve this tension. Nevertheless, the $H\alpha$ EW variation supports our overall conclusion that HD 54879 rotates extremely slowly.

### 3.5. Magnetic Geometry

The strong, complex field of HD 54879, combined with its long rotation period and assumed young age,<fnref n="16" /> make it a particularly interesting candidate for Zeeman Doppler Imaging (ZDI; J. F. Donati & S. F. Brown 1997; N. Piskunov & O. Kochukhov 2002), a technique used to create a surface map of the star's magnetic field. To date, only one O-type star has been mapped in detail using ZDI (Plaskett's star; J. H. Grunhut et al. 2022). The typically large macroturbulent broadening of O-star spectral lines (and therefore Zeeman signatures) makes ZDI mapping challenging. HD 54879 has a relatively low macroturbulence ($v_{mac} = 8 \pm 3$ km s⁻¹, N. Castro et al. 2015; $v_{mac} < 4$ km s⁻¹, T. Shenar et al. 2017; $v_{mac} = 8$ km s⁻¹, G. Holgado et al. 2022), although its slow rotation implies $v \sin i$ is still much less than the macroturbulence. Nevertheless, with the very high S/N and relatively complete rotation-phase coverage provided by our observations, the most important large-scale components of the magnetic geometry can be reconstructed with ZDI.

For this analysis we adopted the rotation period (2562±58 days) and ephemeris from the $\langle B_z \rangle$ data. Using the rotation period and a radius one can derive the equatorial rotation velocity ($v_{eq}$). N. Castro et al. (2015) reported $R = 6.7_{-0.9}^{+1.0} R_\odot$ based on comparing spectroscopic $T_{eff}$ and log g with evolutionary tracks from BONNSAI. T. Shenar et al. (2017) reported $R = 6.1 \pm 1.5 R_\odot$, based on $T_{eff}$ and luminosity with the Stefan-Boltzmann law. Here we adopted the value of T. Shenar et al. (2017), although the exact value has very little impact on the ZDI results, giving us $v_{eq} = 0.12 \pm 0.03$ km s⁻¹.

### 3.5.1. Unno–Rachkovsky Line Model

An initial attempt to fit the Stokes V LSD profiles using ZDI and a line model relying on the weak field approximation (specifically C. P. Folsom et al. 2018a) failed. The width of the Stokes V profile varies and is narrower around phases 0.6–0.8, with the peaks closer to line center than most other phases. With $v \sin i \le 0.12$ km s⁻¹, a weak field model cannot reproduce this variation in Stokes V width. A model that incorporates Zeeman broadening due to a strong magnetic field is needed.

We used a modified version of the ZDIpy code (C. P. Folsom et al. 2018a), that incorporates an Unno-Rachkovsky model for the local line profile. This code is described by S. Bellotti et al. (2023), and follows the work of J.-F. Donati et al. (2008) and J. Morin et al. (2008). The code still uses a spherical harmonic description of the magnetic field (J.-F. Donati et al. 2006), the maximum entropy fitting routine of J. Skilling & R. K. Bryan (1984), and the entropy formulation of M. P. Hobson & A. N. Lasenby (1998) calculated from the spherical harmonic coefficients. The Unno-Rachkovsky calculation follows the description in E. Landi Degl'Innocenti & M. Landolfi (2004). This is essentially a Milne-Eddington line model extended to solve the polarized radiative-transfer equation. Since we are modeling LSD profiles rather than individual lines, we assume a simple triplet Zeeman splitting pattern. Local line opacity and anomalous dispersion profiles are taken to be Voigt and Faraday-Voigt profiles, respectively. J. Morin et al. (2008) and S. Bellotti et al. (2023) include filling factors in their line model, but these are not used here (or equivalently $f_i = f_v = 1$).

The model used here allows separate parameters for limb darkening (in continuum brightness), and the slope of the source function in the Milne-Eddington atmosphere ($\beta$ in equation [9.105] of E. Landi Degl'Innocenti & M. Landolfi 2004). This decouples the effects of brightness decreasing toward the limb and local line strength decreasing toward the limb. In the calculations, the local line profiles are normalized by the local continuum and then scaled by the local limb darkening. However, the linear approximation of a Milne-Eddington atmosphere implies a linear limb darkening. The local continuum flux is

$$I_c = B_0(1 + \beta \cos \theta) \quad (2)$$

(see E. Landi Degl'Innocenti & M. Landolfi 2004, Equations (9.105–9.112)). $B_0$ is the surface value of the source function, and $\theta$ is the angle between the line of sight and surface normal. The local continuum flux relative to the flux at disk center is

$$I_c/I_c^0 = (1 + \beta \cos \theta)/(1 + \beta). \quad (3)$$

A standard linear limb-darkening law (e.g., D. F. Gray 2005) is

$$I_c/I_c^0 = 1 - \eta + \eta \cos \theta \quad (4)$$

for a limb-darkening coefficient $\eta$. One can derive $\beta$ directly from $\eta$:

$$\beta = \eta/(1 - \eta). \quad (5)$$

We used this relation to set the slope of the source function to be consistent with the desired continuum limb darkening. This is convenient since tabulations of $\eta$ are available for a wide range of stars.

<footnote n="16">¹⁶ In a standard picture, a very long rotation period implies the star has had time to experience spindown, and is therefore old. However, recent work investigating the effect of strong magnetic fields in main-sequence stars (Z. Keszthelyi et al. 2019, 2020) and studying the outcomes of stellar merger events (F. R. N. Schneider et al. 2019, 2020) suggests there may be other pathways that produce young, slowly rotating magnetic massive stars. Complexity of the magnetic field could also be an indicator of stellar age, since higher-order multipolar components of the field decay first (P. Charbonneau 2013).</footnote>
### 3.5.2. Optimal Model Parameters

The ZDI model line requires a wavelength and effective Landé factor, which were set to the normalizing values for the LSD profile (500 nm and 1.2). A limb-darkening coefficient of 0.33 was used (A. Claret & S. Bloemen 2011), which implies a slope for the source function in the Unno–Rachkovsky model of 0.4925. The local line opacity profile Gaussian width (combination of thermal and turbulent broadening), Lorentzian width (pressure broadening), and strength were fit using the Stokes I and V profiles. Additional broadening by a Gaussian instrumental profile of R = 65,000 was assumed. The strong magnetic field produces a non-negligible amount of broadening in Stokes I, thus this fitting process was done iteratively with the magnetic mapping process described below. The Gaussian (thermal + turbulent) broadening dominates the total line broadening, but is ~1 km s⁻¹ smaller in our final model than suggested by a nonmagnetic model. The final Gaussian width was 7.2 km s⁻¹ (√2 σ of the Gaussian), as the midpoint between optimal values for Stokes I and Stokes V. This range leads to a potential systematic uncertainty of ±0.5 km s⁻¹, which spans the optimal values for Stokes I and V. The Lorentzian width was 1.1 km s⁻¹ (γ, half-width at half-maximum).

The very low v sin i of HD 54879 limits the information that can be reconstructed with ZDI. Smaller surface features, corresponding to higher-degree spherical harmonics, are not resolved in the Stokes V observations. However, even lacking Doppler resolution, the few lowest-degree spherical harmonics can be resolved by rotational modulation with dense phase coverage. We restricted the maximum degree *l* of the reconstructed spherical harmonics to *l*<sub>max</sub> = 5. In testing a larger *l*<sub>max</sub>, the higher-degree spherical harmonic coefficients tended to 0 due to the regularization. Even modes with *l* > 3 had relatively small coefficients, suggesting we have little sensitivity to these modes. The toroidal components of the magnetic field are also virtually undetectable in Stokes V at this v sin i. At higher v sin i the toroidal components are mostly detectable in V toward the limb of the star as regions of field with opposite orientation, but when v sin i = 0 these regions cancel out. Consequently we restricted the magnetic geometry to a potential poloidal field. In terms of the vector spherical harmonic coefficients<fnref n="17" /> defined by J.-F. Donati et al. (2006; equations [2–4]), this restriction is equivalent to setting α<sub>l,m</sub> = β<sub>l,m</sub> and γ<sub>l,m</sub> = 0. If we allow the toroidal field to be free, it tends toward 0 due to the regularization; but this does not represent an absence of toroidal field, only our lack of sensitivity to this field component.

The inclination of the rotation axis (*i*) is a significant uncertainty for this magnetic mapping. We fixed *i* by searching for the value that optimizes the ZDI result. We ran a grid of ZDI models varying *i* from 1° to 90° (in steps of 1°), with *v*<sub>eq</sub> = 0.12 km s<sup>-1</sup>, and chose the model with the highest entropy that converged the target χ<sup>2</sup><sub>ν</sub> (set to 1.5). To find a formal uncertainty on this inclination we repeated the grid search, but instead ran the code fitting to a target entropy rather than target χ<sup>2</sup><sub>ν</sub>, with the target entropy set to the previous best value. This produced a curve of χ<sup>2</sup><sub>ν</sub> as a function of *i*, for models with the same entropy, and the variation in χ<sup>2</sup><sub>ν</sub> about the minimum can be used to estimate an uncertainty. This follows

the procedure of P. Petit et al. (2002), discussed more by C. P. Folsom et al. (2018a) for *i*, and the approach was found to work well for cool stars with very low v sin *i* by C. P. Folsom et al. (2018b) and C. P. Folsom (2020). This produced a best inclination of *i* = 74° ± 3°, although this may underestimate the real uncertainty due to the very high S/N of the LSD profiles. Varying the Gaussian line width by ±0.5 km s⁻¹ varies *i* by at most 1°, so this is a small contribution to the *i* uncertainty. Varying *v*<sub>eq</sub> by its ±0.03 km s⁻¹ uncertainty, varying the rotation period by its uncertainty, or changing the limb-darkening coefficient by ±0.1, all change *i* by less than 1°.

We used a similar grid search approach to verify the rotation period, following P. Petit et al. (2002). Running a grid of ZDI models for a range of periods, we find P = 2580 ± 30 days, which is in good agreement with the value from ⟨B<sub>z</sub>⟩ of 2562<sup>+63</sup><sub>-58</sub> days. Since the period from the ZDI search is more model dependent, and subject to more systematic uncertainties that are harder to quantify, we adopt the period from ⟨B<sub>z</sub>⟩ in the ZDI analysis and the rest of this study.

### 3.5.3. The Magnetic Map

The final best fit to the Stokes V LSD profiles is shown in Figure 8, and the corresponding magnetic map is shown in Figure 9. This result fits the data to χ²<sub>ν</sub> = 1.5, which we consider acceptable given the small error bars in the LSD profiles. The surface average magnetic field strength is 2954 G, reaching a maximum of 6961 G at a point 20° from the rotation axis near phase 0. The dipole component has a polar strength of 3939 G and obliquity of the positive pole relative to the rotation axis of 168° at a longitude of 128° (phase 0.64; rotation phase runs in the opposite direction to longitude in this model). The magnetic geometry can be described in terms of the fraction of magnetic energy (proportional to ∮B² dΩ) in different spherical harmonic components. We find 74% of the magnetic energy in the dipole (*l* = 1) components, 23% in the quadrupole (*l* = 2) components, and 2% in the octupole (*l* = 3) components. The field is largely axisymmetric (about the rotation axis) with 89% of the energy in axisymmetric (*m* = 0) components. Parameters describing the magnetic field are presented in Table 2.

The very low v sin i of HD 54879 leads to some large uncertainties in the magnetic map. Most obviously, only the largest-scale components of the magnetic field are reconstructed. Components up to the octupole are resolved through the good rotation-phase coverage, but beyond that (for *l* > 3) the geometry is largely unresolved. The toroidal components of the magnetic field are similarly unresolved, due to the lack of Doppler resolution. The toroidal field is typically weak in O and B stars (e.g., J. H. Grunhut et al. 2022), but in a few cases it is important (J.-F. Donati et al. 2006).

The uncertainty on v sin i has a negligible impact on the magnetic field model, as do likely variations in the limb-darkening coefficient. Varying the rotation period by ±1σ changes the magnetic field strength and geometry (energy fractions) by less than 1%. Increasing the inclination by 1σ increases the magnetic field strength by ~500 G and shifts 10% of the energy into the dipole component. Decreasing the inclination by the same amount decreases the magnetic field strength by ~500 G but shifts 20% of the dipolar energy to higher multipoles. Increasing the Gaussian width by 0.5 km s⁻¹ increases the magnetic field strength by ~100 G and shifts ~3% of the magnetic energy away from the dipole, and also

<footnote n="17">¹⁷ In this formalism α<sub>l,m</sub> describes the radial field, β<sub>l,m</sub> describes the tangential poloidal field, and γ<sub>l,m</sub> describes the tangential toroidal field.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** ZDI model fits (red) to the observed Stokes V LSD profiles (black). Rotation cycles are labeled on the right. Profiles are shifted vertically for clarity. Error bars are plotted but are often smaller than the points.</figcaption></figure>

causes the ZDI code to converge very slowly to the target $\chi_{\nu}^2$ of 1.5 (implying possible overfitting, and that for this parameter combination a larger target $\chi_{\nu}^2$ likely should be chosen). Decreasing the Gaussian width by 0.5 km s$^{-1}$ decreases the magnetic field strength by ~30 G, and shifts the magnetic energy toward the dipole by ~1%.

The variation in $H\alpha$ EW correlates with the magnetic field map. The maxima in $H\alpha$ emission occur near phases 0.1 and 0.6, which correspond to places where the positive radial field crosses above the equator, and maxima in the meridional field occur. Minima in $H\alpha$ emission occur near phases 0.4 and 0.9, where the negative radial field spans the largest latitude range. This correlation between surface-magnetic-field structure and $H\alpha$ emission strongly suggests that the stellar wind is being influenced by the magnetic field. Material from the wind can become trapped in closed magnetic field loops above the surface of the star, leading to enhanced $H\alpha$ emission. HD 54879 has been noted to have an unusually strong $H\alpha$ emission for its spectral type (T. S. Boyajian et al. 2007; N. Castro et al. 2015). This could be explained by the very

strong axisymmetric dipolar component of the magnetic field, which should lead to large closed field loops being observable at all rotation phases. Modeling of the magnetosphere around the star is needed to fully explore this hypothesis, such as the analytic dynamical magnetosphere approach (S. P. Owocki et al. 2016; M. S. Munoz et al. 2020) or magnetohydrodynamic simulations.

**4. Discussion and Conclusions**

The most recent spectropolarimetric data of HD 54879 indicate that the longitudinal magnetic field has once again reached the negative magnetic extremum and begun increasing, implying the full stellar rotation period has now been observed. Our best empirical fit to the longitudinal field data yields a stellar rotation period of $P = 2562_{-58}^{+63}$ days, which is about 7.02 yr.

We revisited the longitudinal magnetic field variation measured from LSD profiles constructed from lines of individual elements, specifically He, O, and Si. While those longitudinal field curves show some small differences in
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Magnetic map from ZDI for the radial, azimuthal, and meridional components of the magnetic field. Presented in a Hammer projection with coordinates of latitude and rotation phase. The color scale gives magnetic field in G. Ticks at the bottom indicate phases of observations. The gray area is the unobservable portion of the star.</figcaption></figure>

amplitude, their similar shapes, compounded by the lack of any significant line profile or EW variability of these elements, leads us to attribute these discrepancies to small differences in the measurement systems associated with the different line masks. In contrast to S. P. Järvinen et al. (2022), we see no evidence for chemical stratification or lateral abundance nonuniformities.

Unlike the majority of other known magnetic O stars, HD 54879's longitudinal field curve shows significant evidence for a complex magnetic topology with a strong quadrupolar component. This is supported by the magnetic map for HD 54879 that we present—the first magnetic map of this star—which suggests a surface-averaged magnetic field strength of 2954 G. The dipole component has a polar strength of 3939 G and is nearly aligned with the rotation axis ($\beta$ = 168°). There is a relatively large quadrupolar component to the magnetic field, with 74% of the magnetic energy in the dipole, and 23% in quadrupolar components. HD 54879 is indeed one of the most strongly magnetic O-type stars identified to date.

During ZDI modeling we found that the total Gaussian line broadening is 7.2 km s⁻¹, including turbulent and thermal broadening. This is unusually low for an O-type star (e.g., S. Simón-Díaz et al. 2017). J. O. Sundqvist et al. (2013) investigated turbulent broadening in magnetic O stars, finding only NGC 1624-2 with a lower turbulence than we find for HD 54879, while other stars had turbulent broadening of

<table><caption>**Table 2**
Parameters Describing the Magnetic Map</caption><thead><tr><th>Stellar Parameter</th><th>Value</th></tr></thead><tbody><tr><td><i>P</i></td><td>2562<sup>+63</sup><sub>-58</sub> days</td></tr><tr><td><i>v sin i</i></td><td>0.116 ± 0.029 km s<sup>-1</sup></td></tr><tr><td><i>i</i></td><td>74° ± 3°</td></tr><tr><td><strong>Magnetic Parameter</strong></td><td><strong>Value</strong></td></tr><tr><td>mean |<i>B</i>|</td><td>2954 G</td></tr><tr><td>maximum |<i>B</i>|</td><td>6961 G</td></tr><tr><td>dipole <i>B</i><sub>p</sub></td><td>3939 G</td></tr><tr><td>dipole obliquity</td><td>168°</td></tr><tr><td>quadrupole <i>B</i><sub>max</sub></td><td>2299 G</td></tr><tr><td>quadrupole <i>B</i><sub>min</sub></td><td>-2331 G</td></tr><tr><td>quadrupole max. locations<br>(colatitude, longitude)</td><td>87°, 114°</td></tr><tr><td>quadrupole min. locations<br>(colatitude, longitude)</td><td>93°, 294°</td></tr><tr><td>dipolar energy</td><td>74%</td></tr><tr><td>quadrupolar energy</td><td>23%</td></tr><tr><td>octupolar energy</td><td>2%</td></tr><tr><td>total axisymmetry</td><td>89%</td></tr><tr><td>dipolar axisymmetry</td><td>96%</td></tr><tr><td>quadrupolar axisymmetry</td><td>72%</td></tr></tbody></table>

≥20 km s⁻¹. This study relied on known rotation periods from magnetic O stars to remove the ambiguity between rotational and turbulent broadening. They proposed that the very low turbulent velocity in NGC 1624-2 was a result of a strong magnetic field suppressing turbulence. Additionally, they provide a method for estimating the depth to which a magnetic field may stabilize the atmosphere against turbulence, by looking for the depth where the magnetic and gas pressures are equal. That depth is parameterized by the local temperature ($T_0$; J. O. Sundqvist et al. 2013, Equation (3)). Calculating $T_0$ using $T_{eff}$ and log g from T. Shenar et al. (2017) and our mean $|B|$ (2954 G) yields $T_0 = 54 \pm 7$ kK (with the uncertainty dominated by log g). In the strongest magnetic regions this may rise to ~82 kK. J. O. Sundqvist et al. (2013) propose that the turbulence may originate in the atmosphere between 100 and 200 kK, possibly from an iron opacity bump at ~160 kK. J. MacDonald & V. Petit (2019) investigated this question in more detail, and provided estimates of the magnetic field strength necessary to inhibit convection from this iron opacity bump. For the parameters of HD 54879, a vertical field of ~5 kG is needed, which is below our mean field but close to the maximum field value. It appears that the magnetic field of HD 54879 is not sufficient to stabilize the atmosphere down to the iron opacity bump. However, comparing with a TLUSTY model atmosphere for $T_{eff} = 30$ kK and log g = 4.5 (T. Lanz & I. Hubeny 2003), the magnetically stabilized region of the atmosphere ($T < T_0$) reaches down to an optical depth of ~15. Alternatively, HD 54879 is a very cool and relatively high log g O-type star, so if the generation of turbulence is particularly sensitive to $T_{eff}$ and log g (or mass and luminosity) in this parameter range, that may also explain the low turbulent velocity.

In parallel, we also investigated how Hα emission varies with rotation, finding that the EW of the line varies fairly coherently with $\langle B_z \rangle$. The correlation between the magnetic results and the EWs overall supports the interpretation that the variability in Hα is largely driven by a magnetically confined
wind. We present a fit to the Hα EWs using the rotation period derived from $\langle B_z \rangle$. However, assuming that these variations are rotationally modulated, an independent fit using a two-component harmonic function yields a different (longer) period, albeit within the same order of magnitude. While additional stellar monitoring will be needed to resolve this discrepancy, we conclude that the Hα EW variation supports the notion that HD 54879 is an extreme slow rotator: the magnetic O-type star with the second longest rotational period known to this day.

Some properties of the magnetosphere can be estimated without detailed modeling, although the results are rather approximate for a complex magnetic field and dynamical magnetosphere. We calculated the Alfvén radius ($R_A$) and wind confinement parameter ($\eta_*$) from A. ud-Doula et al. (2008). This requires the mass-loss rate in the absence of a magnetic field ($\dot{M}_{B=0}$), which we estimated from J. S. Vink et al. (2001) following their recommendation for the terminal wind speed ($v_\infty$) as 2.6 of the escape velocity (H. J. G. L. M. Lamers et al. 1995). Using $T_{\text{eff}} = 30,500$ K, $\log L/L_\odot = 4.45$, $R_*=6.1 R_\odot$, and $M_* = 14 M_\odot$ from T. Shenar et al. (2017) yields $\log \dot{M}_{B=0} = -7.9$ [$M_\odot \text{yr}^{-1}$] and $v_\infty = 2400 \text{ km s}^{-1}$. We used the dipole approximation of A. ud-Doula et al. (2008) with our dipolar component field strength (3939 G at the pole, 1970 G at the equator), although this neglects some magnetic field strength and complexity, which leads to $\eta_* = 3500$ and $R_A/R_* = 8.0$. The Kepler radius (A. ud-Doula et al. 2008), calculated with our period, is $R_K/R_* = 310$. Thus, the star's $R_A$ is large (larger than all but one other O star from V. Petit et al. 2013), but still much smaller than $R_K$, implying that HD 54879 has one of the largest dynamical magnetospheres among O-type stars, and within that rotation is nearly negligible.

The large Alfvén radius will substantially reduce the real mass-loss rate. Realistic calculations accounting for the 3D nature of the magnetic field are nontrivial, and here we limit ourselves to the most commonly used rough approximations.<fnref n="18" /> The dipole field approximation from A. ud-Doula et al. (2008) for the reduction in mass-loss rate yields $\dot{M}/\dot{M}_{B=0} = 0.09$ ($\log \dot{M} = -8.9$ [$M_\odot \text{yr}^{-1}$]). The large $R_A$ implies significant angular momentum loss and a relatively short spindown timescale of $\tau_J \sim 2.6$ Myr (A. ud-Doula et al. 2009; V. Petit et al. 2013). Combining this with the current very slow rotation of the star implies a spindown age of $t_J \sim 22$ Myr (V. Petit et al. 2013). This is discrepant with the age of 4 or 5 Myr (N. Castro et al. 2015; T. Shenar et al. 2017) from nonmagnetic evolutionary tracks. While the star may have simply reached the main sequence as a very slow rotator, this suggests more sophisticated calculations of the mass-loss and angular momentum loss rates may be needed, accounting for the 3D nature of the magnetic field and wind (e.g., A. ud-Doula et al. 2013; F. A. Driessen et al. 2019; M. S. Munoz et al. 2020; A. ud-Doula et al. 2023).

The evolutionary history of HD 54879 has yet to be fully explored. Previous work by N. Castro et al. (2015) and T. Shenar et al. (2017) used the Bayesian statistical tool and method BONNSAI (F. R. N. Schneider et al. 2014) to infer the age of HD 54879. This method infers the probability distribution of stellar parameters such as age via the comparison of observables to any set of stellar evolution models. Both studies

used a set of main-sequence, single, nonmagnetic massive-star evolution models with Milky Way compositions that were produced by I. Brott et al. (2011). There was statistical agreement between the inferred ages of $5 \pm 1$ Myr (T. Shenar et al. 2017) and $4.0_{-1.2}^{+0.8}$ (N. Castro et al. 2015). Additionally, T. Shenar et al. (2017) inferred an initial mass of $16 \pm 1 M_\odot$. However, these models do not account for the effects of magnetic fields such as the reduction of mass loss with respect to evolutionary timescales, so there are added uncertainties that have not been evaluated. Acknowledging the results of this work, which reveal the complex magnetic topology of the star, suggests that these previous estimates will need to be revised.

## Acknowledgments

This work is based on observations obtained at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The observations at the CFHT were performed with care and respect from the summit of Maunakea, which is a significant cultural and historic site.

We acknowledge the reverence and importance that Maunakea holds within the Hawaiian community. Hundreds of historic sites, archeological remains, shrines, and burials are on its slopes and summit. CFHT operates on the land of the Native Hawaiian people. We stand in solidarity with Native Hawaiians in their demands for shared governance of Maunakea and to preserve this sacred space for Native Hawaiians.

The authors also gratefully acknowledge support for this work through the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy—EXC-2094—390783311.

C.P.F. acknowledges funding from the European Union's Horizon Europe research and innovation program under grant agreement No. 101079231 (EXOHOST), and from the United Kingdom Research and Innovation Horizon Europe Guarantee Scheme (grant No. 10051045).

A.D.U. acknowledges support from NASA under award number 80GSFC21M0002. S.S.S. acknowledges support from Delaware Space Grant College and Fellowship Program (NASA grant 80NSSC20M0045).

The authors also thank Dr. Zsolt Keszthelyi for the insights he offered during discussions of this project in the early stages of this manuscript.

The authors also extend their thanks to the anonymous referee, who provided several thoughtful comments that led to the improvement of this work.

## Data Availability

The data used in this work can be accessed via the public archives maintained by the Canadian Astronomy Data Centre (ESPaDOnS; https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/cfht/; follow the link and then click “ESPaDOnS” under “Data” on the right-hand side of the web page) and by Polarbase (ESPaDOnS and Narval; http://polarbase.irap.omp.eu/).

The authors can also be contacted with any questions about the data included or analyses completed in this work.

<footnote n="18">¹⁸ An alternative approximate mass-loss rate estimate is provided by T. Shenar et al. (2017), who found $\log M_{B=0} \approx 9.0$ [$M_\odot \text{yr}^{-1}$] and $\log \dot{M} \lesssim 10.2$ [$M_\odot \text{yr}^{-1}$], illustrating some of the uncertainties involved.</footnote>
## Software Availability Statement

The spectral normalization performed for this work was computed using the normPlot software package. Least-squares deconvolution profiles were calculated using the LSDpy software package. Both normPlot and LSDpy are contained within the SpecpolFlow software package, which is available for installation from https://pypi.org/project/specpolFlow/. SpecpolFlow is also available via GitHub at https://github.com/folsomcp/specpolFlow. Detailed documentation and tutorials for SpecpolFlow, together with installation instructions, are available on the project's website at https://folsomcp.github.io/specpolFlow/. The packages are free and open access, and are licensed under the MIT (normPlot, LSDpy) or GPL-2.0 (SpecpolFlow) open-source software licenses.

normPlot is available for individual installation via the Python Package Index (PyPI) at https://pypi.org/project/normPlot/. More information on normPlot can be found at https://github.com/folsomcp/normPlot or on the SpecpolFlow website.

LSDpy is available for individual installation via the Python Package Index (PyPI) at https://pypi.org/project/LSDpy/. More information on LSDpy can be found at https://github.com/folsomcp/LSDpy or on the SpecpolFlow website.

## ORCID iDs

C. Erba https://orcid.org/0000-0003-1299-8878

C. P. Folsom https://orcid.org/0000-0002-9023-7890

S. Bellotti https://orcid.org/0000-0002-2558-6920

L. Fossati https://orcid.org/0000-0003-4426-9530

M. E. Shultz https://orcid.org/0000-0003-1387-5044

## References

Aurière, M. 2003, in EAS Publications Series, 9, Magnetism and Activity of the Sun and Stars, ed. J. Arnaud & N. Meunier, 105

Bellotti, S., Morin, J., Lehmann, L. T., et al. 2023, A&A, 676, A56

Boyajian, T. S., Gies, D. R., Baines, E. K., et al. 2007, PASP, 119, 742

Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115

Castro, N., Fossati, L., Hubrig, S., et al. 2015, A&A, 581, A81

Charbonneau, P. 2013, Solar and Stellar Dynamics, Saas-fee Advanced Course (Berlin: Springer)

Claret, A., & Bloemen, S. 2011, A&A, 529, A75

David-Uraz, A., Petit, V., Shultz, M. E., et al. 2021, MNRAS, 501, 2677

Donati, J. F., & Brown, S. F. 1997, A&A, 326, 1135

Donati, J.-F., Howarth, I. D., Jardine, M. M., et al. 2006, MNRAS, 370, 629

Donati, J.-F., Jardine, M. M., Gregory, S. G., et al. 2008, MNRAS, 386, 1234

Donati, J.-F., Semel, M., Carter, B. D., Rees, D. E., & Cameron, A. C. 1997, MNRAS, 291, 658

Driessen, F. A., Sundqvist, J. O., & Wade, G. A. 2019, in IAU Symp. 346, High-mass X-Ray Binaries: Illuminating the Passage from Massive Binaries to Merging Compact Objects, ed. L. M. Oskinova (Cambridge: Cambridge Univ. Press), 45

Erba, C., David-Uraz, A., Petit, V., et al. 2021, MNRAS, 506, 5373

Erba, C., Shultz, M. E., Petit, V., et al. 2021, MNRAS, 506, 2296

Folsom, C. P., Bouvier, J., Petit, P., et al. 2018a, MNRAS, 474, 4956

Folsom, C. P., Fossati, L., Wood, B. E., et al. 2018b, MNRAS, 481, 5286

Folsom, C. P., Ó. Fionnagáin, D., Fossati, L., et al. 2020, A&A, 633, A48

Fossati, L., Castro, N., Schöller, M., et al. 2015, A&A, 582, A45

Folsom, C. P., Erba, C., Petit, V., et al. 2024, JOSS, submitted

Georgy, C., Meynet, G., Ekström, S., et al. 2017, A&A, 599, L5

Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (3rd ed.; Cambridge: Cambridge Univ. Press)

Grunhut, J. H., Wade, G. A., Folsom, C. P., et al. 2022, MNRAS, 512, 1944

Grunhut, J. H., Wade, G. A., Neiner, C., et al. 2017, MNRAS, 465, 2432

Hobson, M. P., & Lasenby, A. N. 1998, MNRAS, 298, 905

Holgado, G., Simón-Díaz, S., Herrero, A., & Barbá, R. H. 2022, A&A, 665, A150

Hubrig, S., Jarvinen, S. P., Scholler, M., & Hummel, C. A. 2020, MNRAS, 491, 281

Hubrig, S., Küker, M., Järvinen, S. P., et al. 2019, MNRAS, 484, 4495

Järvinen, S. P., Hubrig, S., Schöller, M., et al. 2021, MNRAS, 501, 4534

Järvinen, S. P., Hubrig, S., Schöller, M., et al. 2022, MNRAS, 510, 4405

Keszthelyi, Z. 2023, Galax, 11, 40

Keszthelyi, Z., de Koter, A., Gotberg, Y., et al. 2022, MNRAS, 517, 2028

Keszthelyi, Z., Meynet, G., Georgy, C., et al. 2019, MNRAS, 485, 5843

Keszthelyi, Z., Meynet, G., Martins, F., de Koter, A., & David-Uraz, A. 2021, MNRAS, 504, 2474

Keszthelyi, Z., Meynet, G., Shultz, M. E., et al. 2020, MNRAS, 493, 518

Kochukhov, O., Lundin, A., Romanyuk, I., & Kudryavtsev, D. 2011, ApJ, 726, 24

Kochukhov, O., Makaganiuk, V., & Piskunov, N. 2010, A&A, 524, A5

Kupka, F., Piskunov, N., Ryabchikova, T. A., Stempels, H. C., & Weiss, W. W. 1999, A&AS, 138, 119

Kupka, F. G., Ryabchikova, T. A., Piskunov, N. E., Stempels, H. C., & Weiss, W. W. 2000, BaltA, 9, 590

Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269

Landi Degl'Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines, Vol. 307 (Dordrecht: Kluwer)

Lanz, T., & Hubeny, I. 2003, ApJS, 146, 417

MacDonald, J., & Petit, V. 2019, MNRAS, 487, 3904

Marcolino, W. L. F., Bouret, J.-C., Sundqvist, J. O., et al. 2013, MNRAS, 431, 2253

Meynet, G., Eggenberger, P., & Maeder, A. 2011, A&A, 525, L11

Morel, T., Castro, N., Fossati, L., et al. 2015, in IAU Symp. 307, New Windows on Massive Stars, ed. G. Meynet (Cambridge: Cambridge Univ. Press), 342

Morin, J., Donati, J.-F., Petit, P., et al. 2008, MNRAS, 390, 567

Munoz, M. S., Wade, G. A., Nazé, Y., et al. 2020, MNRAS, 492, 1199

Nazé, Y., Vreux, J.-M., & Rauw, G. 2001, A&A, 372, 195

Owocki, S. P., ud-Doula, A., Sundqvist, J. O., et al. 2016, MNRAS, 462, 3830

Petit, P., Donati, J.-F., & Collier Cameron, A. 2002, MNRAS, 334, 374

Petit, V., Keszthelyi, Z., MacInnis, R., et al. 2017, MNRAS, 466, 1052

Petit, V., Owocki, S. P., Wade, G. A., et al. 2013, MNRAS, 429, 398

Piskunov, N., & Kochukhov, O. 2002, A&A, 381, 736

Piskunov, N., Snik, F., Dolgopolov, A., et al. 2011, Msngr, 143, 7

Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525

Rauw, G., Naze, Y., ud-Doula, A., & Neiner, C. 2023, MNRAS, 521, 2874

Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, PhyS, 90, 054005

Ryabchikova, T. A., Piskunov, N. E., Kupka, F., & Weiss, W. W. 1997, BaltA, 6, 244

Schneider, F. R. N., Langer, N., de Koter, A., et al. 2014, A&A, 570, A66

Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Natur, 574, 211

Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2020, MNRAS, 495, 2796

Shenar, T., Oskinova, L. M., Järvinen, S. P., et al. 2017, A&A, 606, A91

Shultz, M., & Wade, G. A. 2017, MNRAS, 468, 3985

Shultz, M. E., Wade, G. A., Rivinius, T., et al. 2018, MNRAS, 475, 5144

Shultz, M. E., Wade, G. A., Rivinius, T., et al. 2019, MNRAS, 490, 274

Silvester, J., Kochukhov, O., & Wade, G. A. 2014, MNRAS, 440, 182

Simón-Díaz, S., Godart, M., Castro, N., et al. 2017, A&A, 597, A22

Skilling, J., & Bryan, R. K. 1984, MNRAS, 211, 111

Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24

Sundqvist, J. O., Petit, V., Owocki, S. P., et al. 2013, MNRAS, 433, 2497

Sundqvist, J. O., ud-Doula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21

ud-Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413

ud-Doula, A., Owocki, S. P., Russell, C., Gagne, M., & Daley-Yates, S. 2023, MNRAS, 520, 3947

ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2008, MNRAS, 385, 97

ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2009, MNRAS, 392, 1022

ud-Doula, A., Sundqvist, J. O., Owocki, S. P., Petit, V., & Townsend, R. H. D. 2013, MNRAS, 428, 2723

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Wade, G. A., Bagnulo, S., Keszthelyi, Z., et al. 2020, MNRAS, 492, L1

Wade, G. A., Donati, J.-F., Landstreet, J. D., & Shorlin, S. L. S. 2000, MNRAS, 313, 851

Wade, G. A., Howarth, I. D., Townsend, R. H. D., et al. 2011, MNRAS, 416, 3160

Wade, G. A., Neiner, C., Alecian, E., et al. 2016, MNRAS, 456, 2

Walborn, N. R., Sota, A., Maíz Apellániz, J., et al. 2010, ApJL, 711, L143

Yakunin, I., Wade, G., Bohlender, D., et al. 2015, MNRAS, 447, 1418