# The effects of surface fossil magnetic fields on massive star evolution: III. The case of τ Sco

Z. Keszthelyi, G. Meynet, Fabrice Martins, A. de Koter, A. David-Uraz

► To cite this version:

Z. Keszthelyi, G. Meynet, Fabrice Martins, A. de Koter, A. David-Uraz. The effects of surface fossil magnetic fields on massive star evolution: III. The case of τ Sco. Monthly Notices of the Royal Astronomical Society, 2021, 504 (2), pp.2474-2492. 10.1093/mnras/stab893. hal-03260798

HAL Id: hal-03260798

https://hal.science/hal-03260798

Submitted on 24 Nov 2021

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
<figure><img src="image_2.png" /></figure>

# The effects of surface fossil magnetic fields on massive star evolution – III. The case of τ Sco

Z. Keszthelyi<sup>1,★</sup> G. Meynet,<sup>2</sup> F. Martins,<sup>3</sup> A. de Koter<sup>1,4</sup> and A. David-Uraz<sup>5,6,7</sup>

Accepted 2021 March 19. Received 2021 March 17; in original form 2020 October 26

## ABSTRACT

τ Sco, a well-studied magnetic B-type star in the Upper Sco association, has a number of surprising characteristics. It rotates very slowly and shows nitrogen excess. Its surface magnetic field is much more complex than a purely dipolar configuration which is unusual for a magnetic massive star. We employ the CMFGEN radiative transfer code to determine the fundamental parameters and surface CNO and helium abundances. Then, we employ MESA and GENEC stellar evolution models accounting for the effects of surface magnetic fields. To reconcile τ Sco's properties with single-star models, an increase is necessary in the efficiency of rotational mixing by a factor of 3–10 and in the efficiency of magnetic braking by a factor of 10. The spin-down could be explained by assuming a magnetic field decay scenario. However, the simultaneous chemical enrichment challenges the single-star scenario. Previous works indeed suggested a stellar merger origin for τ Sco. However, the merger scenario also faces similar challenges as our magnetic single-star models to explain τ Sco's simultaneous slow rotation and nitrogen excess. In conclusion, the single-star channel seems less likely and versatile to explain these discrepancies, while the merger scenario and other potential binary-evolution channels still require further assessment as to whether they may self-consistently explain the observables of τ Sco.

**Key words:** stars: abundances – stars: evolution – stars: individual: τ Sco – stars: magnetic field – stars: massive – stars: rotation.

## 1 INTRODUCTION

τ Scorpii (HD 149438) is a magnetic massive B0.2 V star in the Upper Sco association, which has been well-studied for almost an entire century (e.g. Struve & Dunham 1933; Unsöld 1942; Traving 1955; Aller, Faulkner & Norton 1966; Lamers & Rogerson 1978; Wolff & Heasley 1985; Kilian 1992; Howk et al. 2000). In addition to a large number of observations in the optical band, there is also a wealth of multiwavelength observations, including X-ray (e.g. Macfarlane & Cassinelli 1989; Cohen, Cassinelli & Waldron 1997; de Messieres et al. 2001; Cohen et al. 2003; Mewe et al. 2003; Ignace et al. 2010; Nazé et al. 2014; Fletcher et al. 2018), ultraviolet (e.g. Walborn & Panek 1984; Peters & Polidan 1985; Rogerson & Ewell 1985; Cowley & Merritt 1987; Snow et al. 1994), and infrared (e.g. Waters et al. 1993; Zaal et al. 1999; Repolust et al. 2005).

Spectropolarimetric observations of τ Sco by Donati et al. (2006) led to the discovery of a surface magnetic field, which is unusually complex compared to other B-type stars whose field measurements can usually be reconciled with a dipolar configuration (e.g. Petit

et al. 2013; Shultz et al. 2018). Additional observations confirmed these findings, showing that the magnetic energy density indeed resides in higher order spherical harmonic components, which clearly implies that the field is complex (Donati & Landstreet 2009; Shultz et al. 2019a,b). Based on the circular (Stokes V) spectropolarimetric data set, Kochukhov & Wade (2016) produced possible magnetic field maps of τ Sco. However, the specific geometry needs to be known to produce a unique magnetic field map. This relies on further observations using linear polarimetry (Stokes Q, U) to constrain the magnetic field modulus.

Interestingly, despite the non-uniqueness problem, the magnetic field maps produced from the spectropolarimetric observations are quite reminiscent of the magnetohydrodynamic (MHD) simulations of Braithwaite (2008) who demonstrated that an initial seed field can relax into a stable non-axisymmetric equilibrium. Indeed, the decade-long spectropolarimetry suggests that the field is stable and of fossil origin.

Using the magnetic oblique rotator model (Stibbs 1950), the rotation period of the star is constrained to 41.033 ± 0.002 d (Donati et al. 2006). Depending on the assumed stellar radius, this yields a surface equatorial rotational velocity of ≈5 km s⁻<fnref n="1" />. Considering measurements of the projected rotational velocity (<10 km s⁻<fnref n="1" />, e.g.

<footnote n="1"><sup>1</sup>Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Geneva Observatory, University of Geneva, Maillettes 51, 1290 Sauverny, Switzerland</footnote>

<footnote n="3"><sup>3</sup>LUPM, Université de Montpellier, CNRS, Place Eugène Bataillon, F-34095 Montpellier, France</footnote>

<footnote n="4"><sup>4</sup>Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium</footnote>

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

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

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

<footnote>* E-mail: z.keszthelyi@uva.nl</footnote>
Hardorp & Scholz 1970; Mokiem et al. 2005; Simón-Díaz et al. 2006; Nieva & Przybilla 2014; Cazorla et al. 2017a, and see also Appendix B), it is clear that the present-day surface rotation of $\tau$ Sco is very slow.

Several authors, e.g. Kilian (1992), Morel, Hubrig & Briquet (2008), Przybilla et al. (2010), and Martins et al. (2012) measured the surface abundance of CNO elements. Their results are in excellent agreement, strongly indicating nitrogen excess (log [N/H] + 12 = 8.15 ± 0.15, 8.15 ± 0.20, 8.16 ± 0.12, and 8.15 ± 0.06 by number fraction, respectively). This means that $\tau$ Sco is highly enriched in core-processed material, by at least a factor of 3 compared to the solar baseline. Helium abundance measurements are more uncertain. Previously obtained values (e.g. He/H = 0.085 by number fraction from Wolff & Heasley 1985; He/H = 0.10 ± 0.025 from Kilian 1992; $Y$ = 0.28 ± 0.03 by mass fraction from Przybilla et al. 2010) are very close to the solar helium baseline (Grevesse, Noels & Sauval 1996; Asplund, Grevesse & Sauval 2005; Asplund et al. 2009), which would be compatible with expectations of a single star in its early main sequence evolution. Although still overlapping within uncertainty, Mokiem et al. 2005 obtained He/H = 0.12$^{+0.04}_{-0.02}$ by number fraction, which could point towards a slight excess in the surface helium abundance.

Nieva & Przybilla (2014) measured the stellar parameters of $\tau$ Sco and based on its position on the Hertzsprung-Russell diagram concluded that it is a blue straggler star that is much younger than the association it belongs to, suggesting that it could possibly originate from a stellar merger. Indeed, Ferrario et al. (2009) suggested that the origin of fossil magnetic fields may be stellar merger events. Further work by Schneider et al. (2016) explored this scenario and argued that $\tau$ Sco could be rejuvenated via a merger process.

Recently Schneider et al. (2019) presented 3D MHD simulations of stellar mergers and Schneider et al. (2020) implemented the results into 1D stellar evolution models to follow the long-term evolution of the post-merger object. These simulations showed that i) a seed field can be amplified to a strong magnetic field in the merger process and ii) the post-merger object evolves seemingly as a single star, although with an unusual internal chemical composition. Nonetheless, some important observables of $\tau$ Sco (such as the $\approx 5$ km s$^{-1}$ rotational velocity and the factor of 3 nitrogen enrichment) are not yet compatible with the current modelling efforts as they predict either modest rotation (50 km s$^{-1}$, which is still an order of magnitude larger than the observed value) with no nitrogen excess, or fast rotation (400 km s$^{-1}$) with nitrogen excess (see figs 3 and 5 of Schneider et al. 2020). The only solution thus far that predicts sufficiently slow rotation (but no nitrogen excess) is a model with a constant surface magnetic dipole moment of $\mu_B = 10^{40}$ G cm$^3$ corresponding to an approximately 270 kG field at the surface of a star with a radius of 5 R⊙. Such a strong field is reached in the 3D MHD merger simulation (Schneider et al. 2019), however, since the field strength remains quasi-constant on the main sequence, it is far too strong compared to any known OBA star and, in particular, to $\tau$ Sco's measured surface field strength of a few hundred G (Donati et al. 2006; Shultz et al. 2018, and further details in Appendix B).

These significant findings pose important questions, for example, i) are all fossil magnetic fields generated via a stellar merger event, and ii) can we confidently identify signs of a past stellar merger event and its evolutionary consequences when studying a star? This makes $\tau$ Sco a valuable laboratory to further our understanding of magnetic massive stars and for this reason the current discrepancies between models and observations need to be studied in more detail.

The prospect of $\tau$ Sco being a merger product is intriguing. Of course, the validity of this scenario would be strengthened if one

can firmly exclude the possibility that the star formed and evolved in isolation, having acquired its fossil field during the assembly process. The purpose of this work therefore is, first, to re-assess stellar and surface properties, including its age and nitrogen abundance, and to comprehensively explore whether these empirical characteristics can be reconciled with single-star evolutionary models that account for surface magnetic field effects; and, second, to discuss implications of these properties for merger models.

This paper is part of a series in which we aim to explore the effects of surface fossil magnetic fields on massive star evolution. In the first paper of the series (Keszthelyi et al. 2019, hereafter Paper I), we used the Geneva stellar evolution code (GENEC; Eggenberger et al. 2008; Ekström et al. 2012; Georgy et al. 2013; Groh et al. 2019; Murphy et al. 2021) to explore the cumulative impact of magnetic mass-loss quenching, magnetic braking, and field evolution. In the second paper (Keszthelyi et al. 2020, hereafter Paper II), we implemented and studied massive star magnetic braking in the Modules for Experiments in Stellar Astrophysics (MESA) software instrument (Paxton et al. 2011, 2013, 2015, 2018, 2019), detailing the magnetic and rotational evolution, and confronting the models with a sample of observed magnetic B-type stars from Shultz et al. (2018). A key finding of Paper II is that presently slowly rotating magnetic stars may originate from either slow or fast rotators at the zero-age main sequence (ZAMS).

The paper is organized as follows. In Section 2, we perform atmospheric modelling, particularly focusing on the CNO elements and helium abundance. In Section 3, we carry out numerical experiments with stellar evolution models. In Section 4, we discuss the impact of the initial assumptions and confront empirical evidence with our results. Finally, we summarize our findings in Section 5.

# 2 ATMOSPHERIC MODELLING

## 2.1 Determination of fundamental parameters

We have used the code CMFGEN (Hillier & Miller 1998) to determine the fundamental stellar parameters of $\tau$ Sco. CMFGEN computes atmosphere models in non-LTE, spherical geometry, and includes outflows and line-blanketing. A complete description of the code is presented in Hillier & Miller (1998) and we refer the reader to this publication for further information.

To the best of our knowledge, there is no atmosphere code for hot, massive stars (including non-LTE effects, sphericity, and winds), which includes polarization in the radiative transfer equation. For lower mass stars, Khan & Shulyak (2006) showed that the extra line-blanketing caused by Zeeman splitting had little impact on the emergent spectrum, especially in the field strength regime of $\tau$ Sco. The expected Zeeman splitting of helium and CNO lines is undetectable in individual lines unless the magnetic field is extremely strong (which is why field-detection techniques relying on many lines are used; Donati & Landstreet 2009), hence the use of non-polarized synthetic spectra is justified. Krtička (2018) showed that polarized radiative transfer has little effect on radiative driving so that mass-loss rates should not be affected. However, the presence of a magnetic field confines the expanding wind, leading to non-spherical effects that need to be taken into account when studying the details of stellar winds in magnetic OB stars (e.g. ud-Doula 2017). In the present study, we focus on photospheric parameters which should not be severely affected by the use of 'non-magnetic' models.

To determine the main fundamental parameters, we relied on a grid of CMFGEN models covering the effective temperatures and surface gravities of hot main sequence stars. We also used
the spectrum of τ Sco presented by Martins et al. (2012). It is an average spectrum built from observations conducted with the spectropolarimeter ESPaDOnS at CFHT. It covers the wavelength range 3800–6800 Å at a spectral resolution of about 65 000. In this spectrum, we isolated the classical indicators to determine the effective temperature and surface gravity, namely He I 4471, He II 4542, He II 5412, Hγ and Hβ. We subsequently looked for the best-fitting models using a χ<fnref n="2" /> minimization process. We found $T_{\text{eff}} = 31\,500 \pm 1000$ K and $\log g = 4.2 \pm 0.1$ cm s⁻<fnref n="2" />. For consistency, the process was repeated with additional helium and Balmer lines and consistent results were obtained.

In this process, we used a projected rotational velocity of 9 km s⁻<fnref n="1" /> and a macroturbulence of 10 km s⁻<fnref n="1" />. The former was obtained from the Fourier transform method (see Simón-Díaz & Herrero 2007) applied to O III 5592 and He I 4713. To determine the level of macroturbulence, we convolved spectra from our grid and with parameters close to the final values with a rotational profile, an instrumental one, and a radial-tangential profile mimicking macroturbulence. We constrained the macroturbulence velocity by fitting He I 4713.

To constrain the luminosity of τ Sco, we computed the bolometric correction from $T_{\text{eff}}$ and the relation of Martins & Plez (2006). The extinction was estimated from the colour excess $E(B-V)$ using the intrinsic $(B-V)_0$ colour of Martins & Plez (2006). With the new distance of $d = 195 \pm 42$ pc from Gaia DR2, we finally obtained $\log(L/L_{\odot}) = 4.56 \pm 0.20$. Using the radius obtained from $T_{\text{eff}}$ and $L/L_{\odot}$, as well as the projected rotational velocity, and assuming an inclination of about 90°, we find a rotation period of 36 d. This is in excellent agreement with the 41 d found by Donati et al. (2006) from spectropolarimetric measurements.

## 2.2 Surface abundances

We determined the surface abundances of carbon, nitrogen, oxygen, and helium. For each element, we computed atmosphere models and synthetic spectra with the set of fundamental parameters derived in Section 2.1 but different chemical compositions. We then used a set of selected lines of each element to perform a χ<fnref n="2" /> minimization and find the best-fitting model. These lines included in the analysis are:

(i) Carbon: C III 4056, C III 4068, C III 4070, C II 4075, C III 4665, C III 4667, C II 5133, C II 5144, C II 5151, C II 5250, C III 5254, C IV 5802, C IV 5812, C III 5827, C III 6205, C II 6578, C II 6583, C III 6744.

(ii) Nitrogen: N II 3995, N II 4004, N II 4035, N III 4044, N III 4196, N III 4216, N II 4447, N III 4511, N III 4515, N II 4602, N II 4607, N III 4615, N II 4621, N III 4634, N III 4640, N II 4788, N II 4803, N II 4995, N II 5001, N II 5005, N II 5011, N II 5026, N II 5045, N II 5678, N II 5680.

(iii) Oxygen: O II 3913, O II 3954, O II 3983, O II 4277, O II 4278, O II 4283, O II 4305, O II 4316, O II 4318, O II 4321, O II 4354, O II 4367, O II 4370, O II 4397, O II 4415, O II 4417, O II 4453, O II 4489, O II 4491, O II 4592, O II 4597, O II 4603, O II 4610, O II 4662, O II 4676, O II 4678, O II 4700, O II 4705, O III 5592.

We obtained the following values of surface abundances in number fractions: $C/H = 1.8^{+0.9}_{-0.7} \cdot 10^{-4}$, N/H= $1.8^{+0.8}_{-0.6} \cdot 10^{-4}$, and $O/H= 4.1^{+2.0}_{-1.3} \cdot 10^{-4}$. For these determinations, we used a macroturbulent velocity of 2 km s⁻<fnref n="1" />. A larger value would lead to excessive broadening, as, for instance, evidenced in the N II doublet at 5001 Å. The reduced microturbulence compared to the study of Martins et al. (2012) explains the slightly larger N/H ($1.8 \cdot 10^{-4}$ in the present study versus $1.4 \cdot 10^{-4}$ by Martins et al. 2012). The best-fitting synthetic spectrum is shown in Fig. 1.

We also investigated the possibility that τ Sco is enriched in helium. To this extent, we ran models with He/H between 0.1 (the reference value) and 0.2. We used the following lines to estimate the goodness of fit:

(i) Helium: He I 4026, He II 4200, He I 4388, He I 4471, He II 4542, He I 4713, He I 4921, He II 5412, He I 6680.

The best fits were obtained for models with He/H = 0.11 ± 0.01. This corresponds to Y = 0.30 ± 0.02 in mass fraction, implying a modest enrichment compared to the baseline value of Y ini = 0.266. The results are summarized in Table 1.

## 3 EVOLUTIONARY MODELLING

Evolutionary models were computed with two different codes, Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019) and the Geneva stellar evolution code (GENEC, Eggenberger et al. 2008; Ekström et al. 2012; Georgy et al. 2013; Paper I). The two codes are sufficiently different that they provide independent tests. In particular, the implementations of angular momentum transport and loss are different (see Paper II, Appendix A).

### 3.1 MESA modelling: setup

The general model parameters we adopt in MESA are as follows. A solar metallicity of Z = 0.014 is adopted with the Asplund et al. (2009) mixture of metals; isotopic ratios are from Lodders (2003). The mixing efficiency in the convective core is treated by adopting $\alpha_{\text{MLT}} = 2.0$. Exponential overshooting is used above the convective core with $f_{\text{ov}} = 0.034$ and $f_0 = 0.006$ (Herwig 2000; Paxton et al. 2013). This would approximately correspond to extending the convective core by 25 per cent of the local pressure scale height (thus $\alpha_{\text{ov}} = 0.25$) in a non-rotating model.<fnref n="1" /> The OPAL opacity tables are used (Rogers & Iglesias 1992). The mass-loss scheme is adopted from de Jager, Nieuwenhuijzen & van der Hucht (1988).<fnref n="2" /> This mass-loss rate is then systematically scaled by the magnetic mass-loss quenching parameter $f_B$ and rotational enhancement $f_{\text{rot}}$ (see Paper II).

To account for rotationally induced instabilities which lead to chemical mixing, we use MESA's default parametrization in a fully diffusive approach. The diffusion coefficients arising from dynamical and secular shear instabilities (Endal & Sofia 1978; Pinsonneault et al. 1989), Eddington-Sweet circulation (Eddington 1925; Sweet 1950), Solberg-Høiland (Solberg 1936; Høiland 1941), and Goldreich-Schubert-Fricke instabilities (Goldreich & Schubert 1967; Fricke 1968) are included. The Spruit-Tayler dynamo (Spruit

<footnote n="1">¹A range of $\alpha_{\text{ov}} = 0.1-0.5$ is presently considered plausible based on asteroseismic measurements as well as multidimensional hydrodynamic simulations (e.g. Herwig 2000; Meakin & Arnett 2007; Moravveji et al. 2015; Pápicz et al. 2017;). See recently Kaiser et al. (2020) and references therein. We test other values in Appendix A, none the less, the precise choice of the overshooting parameter does not play an important role here as without efficient envelope mixing, the surface composition will remain unaltered.</footnote>

<footnote n="2">²This specific choice is made because the de Jager et al. (1988) mass-loss rates in this mass and temperature range are slightly higher than the Vink, de Koter & Lamers (2001) rates, consequently they (very minimally) also aid the spin-down of the star (e.g. Keszthelyi, Puls & Wade 2017). However, we tested models with both mass-loss prescriptions and found that they lead to completely negligible differences in this study. This is expected since the mass-loss rates during the early evolution of the models are on the order of $\dot{M} \approx 10^{-9} - 10^{-8} M_{\odot} \text{yr}^{-1}$.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Best-fitting model (red) compared to the observed spectrum (black) of τ Sco.</figcaption></figure>

<table><caption>**Table 1.** Stellar parameters and surface abundances in number fractions obtained from atmospheric modelling.</caption><tr><td>log (L/L<sub>☉</sub>)</td><td>4.56 ± 0.20</td></tr><tr><td>T<sub>eff</sub> (kK)</td><td>31.5 ± 1</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.2 ± 0.1</td></tr><tr><td>v sin i (km s<sup>-1</sup>)</td><td>9</td></tr><tr><td>v<sub>mac</sub> (km s<sup>-1</sup>)</td><td>10</td></tr><tr><td>v<sub>mic</sub> (km s<sup>-1</sup>)</td><td>2</td></tr><tr><td>C/H</td><td>1.8<sup>+0.9</sup><sub>-0.7</sub> · 10<sup>-4</sup></td></tr><tr><td>N/H</td><td>1.8<sup>+0.8</sup><sub>-0.6</sub> · 10<sup>-4</sup></td></tr><tr><td>O/H</td><td>4.1<sup>+2.0</sup><sub>-1.3</sub> · 10<sup>-4</sup></td></tr><tr><td>He/H</td><td>0.11 ± 0.01</td></tr></table>

2002; Tayler 1973) is not adopted. To consider the inhibiting impact of composition gradients, $\nabla\mu$ is scaled by $f_\mu = 0.05$ when calculating stability criteria for chemical mixing (Yoon, Langer & Norman 2006; Brott et al. 2011; Paxton et al. 2013).

In MESA, angular momentum transport is modelled as a diffusive process. This means that angular momentum is only transported in one (outward) direction. The stellar core has a constant angular velocity profile. Close to the core-envelope boundary (at $q = 0.4$, i.e. at the layer encompassing 40 per cent of the total mass), we set the lower boundary to apply magnetic braking (see below). Therefore, we assume that the magnetic field spreads through the stellar envelope, leading to a very efficient angular momentum transport. An appropriate transport equation, however, relies on the properties of the internal magnetic field (e.g. strength, geometry, and obliquity), which are unknown. Instead, we adopt a uniformly high diffusion coefficient ($D = 10^{16}$ cm$^2$ s$^{-1}$) throughout the stellar envelope (from the photosphere down to $q = 0.4$) to account for
<figure><img src="image_4.png" /><figcaption>**Figure 2.** MESA models are shown on the Kiel (left) and Hunter-P (right) diagrams. Upper panels: varying the initial mass for $v_{ini} = 350$ km s$^{-1}$. Middle panels: varying the initial rotation rates for $M_{ini} = 18$ M$_{\odot}$. Lower panels: adopting slow rotation with $v_{ini} = 150$ km s$^{-1}$ while varying mixing efficiency for $M_{ini} = 18$ M$_{\odot}$. The colour-coding scales with stellar age. Yellow crosses indicate an age of 6 Myr. Observations from various authors are shown. The rotation period (here and hereafter, $P_{rot} = 41.033 \pm 0.002$ d) is adopted from Donati et al. (2006).</figcaption></figure>
### 3.2.3 Increasing the mixing efficiency in a slow rotator model

These experiments show that rotational mixing and magnetic braking tend to work against each other. In MESA's diffusive scheme, when magnetic braking is efficient, the (surface and thus internal) rotational velocity is lower, and thus chemical mixing is reduced.

To simulate a more efficient mixing for fixed stellar parameters, we increase the scaling factor $f_c$ which multiplies the sum of all diffusion coefficients considered for chemical mixing<fnref n="6" /> (see Heger, Langer & Woosley 2000; Brott et al. 2011; Paxton et al. 2013).

Fig. 2 (lower panels) shows this parameter test for an $18 M_\odot$ model with slow rotation, adopting $v_{ini} = 150$ km s$^{-1}$. The model with ‘standard’ mixing efficiency ($f_c = 0.033$, calibrated by Heger et al. 2000) does not lead to any nitrogen enrichment. When $f_c$ is increased by a factor of 3, there is some enrichment which is compatible with the lower limit of the observations. Increasing $f_c$ by an order of magnitude allows to reach the upper limit of the measured nitrogen abundance of $\tau$ Sco - in this particular model setup. This means that for even lower initial rotation rates than we considered here, an extremely (likely unphysically) high mixing efficiency would need to be invoked. Thus we see no feasibility to decrease the initial rotation rate below $v_{ini} \approx 150$ km s$^{-1}$.

### 3.2.4 Increasing the braking efficiency in a fast rotator model

A major shortcoming of all the previous parameter tests has been that the current long rotation period of $\tau$ Sco – together with the position in the Kiel diagram and the observed surface nitrogen abundance – could not be self-consistently (i.e. with the same age) reproduced by the models. Therefore, we continue our thought experiment with testing whether a more efficient braking mechanism may help to overcome this problem. In the upper panels of Fig. 3, we show models with initially $M_{ini} = 18 M_\odot$ and fast rotation<fnref n="7" /> with $v_{ini} > 500$ km s$^{-1}$, using a ‘usual’ ($f_{MB} = 1$, cf. equation 1), a factor of 3, and a factor of 10 more efficient magnetic braking (lower panels of Fig. 3). When magnetic braking is more efficient by a factor of 10, the observed spin rate of $\tau$ Sco can be recovered within 6 Myr. However, the efficient spin-down leads to very inefficient rotational mixing and thus no nitrogen is seen on the stellar surface. When $f_{MB} = 3$, chemical mixing can still remain efficient during the early evolution and thus the observed values are well-approximated by the model. However, the spin rate of this model after 6 Myr is still an order of magnitude larger than the observed one.

### 3.2.5 Varying both the mixing and braking efficiencies

The only experiment which succeeded to reproduce the long rotation period of $\tau$ Sco within 6 Myr (which is our criterion to obtain a self-consistent solution to match all key observables) is the one where the magnetic braking efficiency is increased by a factor of 10. To see if it is possible at all – within this thought experiment – to obtain nitrogen excess in this case, we compute models where additionally the mixing efficiency $f_c$ is increased by a factor of 10 for a model

with $v_{ini} = 250$ km s$^{-1}$, while it is increased by a factor of 3 for a model<fnref n="8" /> with $v_{ini} \approx 550$ km s$^{-1}$ (Fig. 3, middle panels).

In the lower panel of Fig. 3, we also show the measured $^{14}$N/$^{12}$C mass fraction as a function of $^{14}$N/$^{16}$O mass fraction. This confirms that the nitrogen abundance alone is not anomalous, instead the overall CNO abundances are consistent with the surfacing of stellar core material. That is, initially, carbon becomes depleted at the expense of nitrogen, whereas oxygen remains approximately constant (Przybilla et al. 2010; Maeder et al. 2014). After the CN cycle is in equilibrium, a slight oxygen depletion takes place leading to a further increase of the nitrogen abundance. It is thus clear that these observations cannot be reconciled with an initially slowly rotating single star model where the surface does not reflect the core abundance (see the Hunter-P diagram in Fig. 2). The evolutionary models begin with approximately N/C = 0.3 and N/O = 0.1, however, the tracks only represent the time evolution until mixing is efficient. When rotation becomes slow ($P_{rot} \gtrsim 10$ d in our experiment), the surface CNO abundance ratios remain unaltered. Thus main sequence models with inefficient mixing would appear as single points on this diagram, close to the ZAMS values.

These artificially engineered models, for the first time, can approximately reproduce not only the four strict observables but also the overall CNO abundance ratios with a consistent stellar age of less than 6 Myr. The question therefore is whether there might be a physically meaningful reason to justify large departures in the efficiency of commonly used prescriptions for rotational mixing and magnetic braking.

### 3.2.6 Tests with different magnetic field evolution

Thus far, we made the simplifying assumption that the surface magnetic field strength did not vary in time. However, magnetic field evolution is presently not well-constrained and therefore we now experiment with cases where the surface magnetic field strength can change in time.

The usual, first order estimate is based on magnetic flux conservation (following Alfvén's theorem, Alfvén 1942),

$$F \propto B_{\text{surf}}(t) R_{*}^{2}(t) = \text{const.} = B_{\text{surf}}(t=0) R_{*}^{2}(t=0), \quad (3)$$

where $B_{surf}(t=0)$ and $R_*(t=0)$ are the ZAMS magnetic field strength (which is assumed) and stellar radius (which is calculated). In this case, the magnetic field strength only varies as a function of the stellar radius.

Since an initially stronger magnetic field aids the spin-down of the star, considerations have been given to a magnetic field decay scenario, which we assume to have the form of

$$B_{\text{surf}}(t) = B_{\text{surf}}(t=0) \exp(-f_{\text{dec}} t/\tau), \quad (4)$$

where $f_{dec}$ is an arbitrary scaling factor to which we refer to as the decay efficiency, $t$ is the time, and $\tau$ is a characteristic time-scale. We set $\tau = 12$ Myr, which is the approximate main sequence lifetime of our models, so that the magnetic field strength would weaken to roughly 60 per cent of its initial value in 6 Myr with a decay efficiency of unity. With $f_{dec} = 2$ and 4, the field strength becomes about 30 per cent and 15 per cent, respectively, of its initial value in 6 Myr. We note that since the change in stellar radius is very modest in the first half of the main sequence, these arbitrary field decay

<footnote n="8">⁸In this fast-rotating model too, the rotational velocity undergoes some adjustment, therefore the initial value stated (at our ZAMS definition) should not be considered as an exact limit.</footnote>

<footnote n="6">⁶The Eddington-Sweet circulation is the dominant term throughout the stellar envelope.</footnote>

<footnote n="7">⁷The exact values of $v_{ini}$ are somewhat different at the ZAMS as a consequence of the model relaxation, applied magnetic braking, and the ZAMS definition itself. We use the criteria that 0.3 per cent of core hydrogen has already burnt to define the ZAMS. None the less, the exact initial rotational velocities do not significantly affect the trend that we show with this parameter test.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Upper panels: same as Fig. 2 but with $M_{ini} = 18$ M$_{\odot}$ and adopting fast rotation with $v_{ini} = 550$ km s$^{-1}$ while varying the braking efficiency. Middle panels: with $M_{ini} = 18$ M$_{\odot}$ and $v_{ini} = 250$ km s$^{-1}$, $f_c = 0.33$ and $f_{MB} = 10$ (dotted), and $v_{ini} \approx 550$ km s$^{-1}$, $f_c = 0.1$ and $f_{MB} = 10$ (dashed). Lower panel: $^{14}$N/$^{12}$C mass fraction as a function of $^{14}$N/$^{16}$O mass fraction. The derived CNO abundance ratios are shown with blue marker. The evolutionary models (same as in the middle panels) which can mix core material to the surface evolve towards the top right of the diagram. The end point of the model is where surface rotation is slow enough ($P_{rot} \approx 10$ d, see middle panel) that the surface abundance does not change anymore. Since this is reached in about the first 2–3 Myr, the further time evolution leads to no changes on this diagram.</figcaption></figure>
scenarios are practically identical to a (perhaps more commonly known) flux decay scenario. Contrasted with the flux conservation scenario (equation 3), we now have:

$$F \propto B_{\text{surf}}(t) R_{*}^{2}(t) = B_{\text{surf}}(t=0) \exp(-f_{\text{dec}} t / \tau) R_{*}^{2}(t). \quad (5)$$

In Fig. 4, models with different magnetic field evolution scenarios are shown. The top panels show models with magnetic flux conservation (equation 3), whereas the middle and lower panels show models with magnetic field decay (equations 4–5).

With magnetic field evolution, the models can have an initially stronger magnetic field and consequently initially more efficient magnetic braking. Although time dependent, this is somewhat similar to the previous tests where $f_{MB}$ was used to obtain a more efficient magnetic braking. This means that magnetic field evolution, in principle, allows for adopting a higher initial rotational velocity and thus a higher initial ('natural') efficiency of rotational mixing.

Assuming magnetic flux conservation, a field with an initial strength of 1 kG (assumed dipolar) well approximates the present-day field strength of $\tau$ Sco, however, leads to a small impact in its overall evolution compared to the reference model (600 G polar field strength constant in time). An initial 4 kG field produces efficient braking albeit no nitrogen enrichment. This model also leads to a magnetic field strength over the entire evolution (except close to the TAMS) which is far too strong (> 1 kG) to be compatible with spectropolarimetric measurements (see Appendix B).

When magnetic field decay is considered, the model with an initial 3 kG field strength (assumed dipolar) and decay efficiency $f_{dec}$ of unity produces sufficient nitrogen excess, however, does not spin-down fast enough (middle panels). The model with an initial 9 kG field strength does produce efficient braking, and while the long rotation period is recovered, no mixing is achieved. Increasing the mixing efficiency in this model by a factor of 3 leads to an acceptable, self-consistent match with the observables of $\tau$ Sco. However, the observable field strength would remain above 3 kG in the first 6 Myr, which is at odds with the observations.

Experiments with 30 kG initial field strengths allow to reproduce the current long rotation period (lower panels). Increasing $f_c$ by a factor of 10 is required to reach nitrogen excess. However, a decay efficiency of 2 still yields a magnetic field strength well above 1 kG in 6 Myr. $f_{dec} = 4$ results in an acceptable self-consistent solution.

From these tests, we conclude that while magnetic field evolution can, to some extent, alleviate the large magnetic braking efficiency ($f_{MB} = 10$) that was needed in the model with a constant 600 G field strength, it faces the challenge to reach a sufficiently low, sub-kG field strength in 6 Myr. This requires a high magnetic field decay efficiency. Importantly, in all tests, the need for a more efficient chemical mixing is still present.

### 3.3 GENEC modelling: setup

In GENEC, we adopt similar modelling assumptions as Ekström et al. (2012). A solar metallicity of $Z = 0.014$ is used with the Asplund et al. (2009) mixture of metals except for neon (Ekström et al. 2012), and isotopic ratios are from Lodders (2003). The adopted mixing efficiency in the convective core is $\alpha_{MLT} = 1.6$. A step overshooting method is applied with $\alpha_{ov} = 0.1$. The opacity tables are adopted from OPAL (Rogers & Iglesias 1992). Mass-loss rates are calculated following the prescription of Vink, de Koter & Lamers 2000 and Vink et al. 2001, multiplied by a factor of 0.85.

Chemical element transport is modelled as a diffusive process (Pinsonneault et al. 1989), adding up from three terms: vertical shear, horizontal turbulence, and meridional currents. Following

Zahn (1992), the latter two are combined into one effective diffusion coefficient, $D_{eff}$ (see Meynet et al. 2013 and references therein). In this work, we test two prescriptions. In one case (which will be referred to as D11 hereafter), the shear is adopted from Maeder (1997) as

$$D_{\text{shear}} = \frac{H_p}{g\delta} \frac{K}{[\frac{\varphi}{\delta}\nabla_{\mu} + (\nabla_{ad} - \nabla_{\text{rad}})]} \left(\frac{9\pi}{32}\Omega \frac{d\ln\Omega}{d\ln r}\right)^2, \quad (6)$$

where $H_p$ is the local pressure scale height, $g$ is the local gravitational acceleration, $\delta$ and $\varphi$ are derivatives from the equation of state, $K$ is the thermal diffusivity, $\Omega$ is the angular velocity, and $r$ is the distance from the centre. The horizontal turbulence is adopted from Zahn (1992) as

$$D_h = r|2V(r) - \alpha U(r)|, \quad (7)$$

where $\alpha = \frac{1}{2} \frac{d\ln(r^2\Omega)}{d\ln r}$ and $U(r)$ and $V(r)$ are the horizontal and vertical components of the meridional circulation.

In the second case (which we will refer to as D22 hereafter), the shear is adopted from Talon & Zahn (1997) as

$$D_{\text{shear}} = \frac{H_p}{g\delta} \frac{K + D_h}{[\frac{\varphi}{\delta}\nabla_{\mu}(1 + \frac{K}{D_h}) + (\nabla_{ad} - \nabla_{\text{rad}})]} \left(\frac{9\pi}{32}\Omega \frac{d\ln\Omega}{d\ln r}\right)^2, \quad (8)$$

and the horizontal turbulence is adopted from Maeder (2003) as

$$D_h = Ar(r\Omega(r)V(r)|2V(r) - \alpha U(r)|)^{1/3}, \quad (9)$$

where $A = (3/400n\pi)^{1/3}$ with $n$ being the number of axial rotations (Maeder 2003).

Angular momentum transport is modelled by using an advecto-diffusive equation which accounts for the radial component of meridional currents (the advective term) and shears (the diffusive term). The meridional currents are an advective process by nature. The shear term is modified when using the D11 ($D_{shear}$ adopted from Maeder 1997) or D22 ($D_{shear}$ adopted from Talon & Zahn 1997) schemes in the models. Both of these cases, without internal magnetic fields, allow for radial differential rotation and lead to a weaker core-envelope coupling than in solid-body rotating models. We do not test solid-body rotating models here because it was shown in previous works that it leads to less surface enrichment than differentially rotating models (Meynet, Eggenberger & Maeder 2011, Paper I).

The effects of the surface magnetic field are modelled via magnetic braking, which is implemented as a boundary condition for internal angular momentum transport (Meynet et al. 2011; Georgy et al. 2017, Paper I, Paper II). We refer the reader to Paper II Appendix B, where the MESA and GENEC implementations of magnetic braking are detailed and contrasted. Since in GENEC only the outermost layers are ascribed to lose specific angular momentum, the use of angular momentum transport without a strong coupling means that significant shears can develop in the outer part of the stellar envelope, while meridional currents remain efficient to transport chemical elements close to the stellar core. (Note that in MESA, we model the opposite scenario: shears remain efficient close to the core but (weak) meridional circulation dominates the transport in the outer envelope.) The equatorial magnetic field strength is set to $B_{eq} = 300$ G (corresponding to 600 G polar) and is kept constant over time.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** MESA models with $M_{ini} = 18$ M$_\odot$ and $v_{ini} > 500$ km s$^{-1}$ are shown. Upper panels: Models with magnetic flux conservation compared to the reference model ($B_p = 600$ G). Standard efficiencies are used for $f_c$ and $f_{MB}$. Middle and lower panels: Models with magnetic field decay. The magnetic braking efficiency is set with $f_{MB} = 1$ but the decay efficiency $f_{dec}$ and rotational mixing efficiency $f_c$ are varied (cf. equations 1 and 4). The two models with initially 8 kG fields in the middle right-hand panel overlap.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Shown are models computed with the Geneva code on the Kiel (left) and Hunter-P diagrams (right). The colour-coding scales with stellar age. Yellow crosses indicate an age of 6 Myr.</figcaption></figure>

## 3.4 GENEC modelling: results and analysis

Some of the major differences between GENEC and MESA are the treatment of angular momentum transport and the way magnetic braking is applied. We seek to probe here whether the same four strict observables of τ Sco could be reconciled when using the modelling assumptions as described above. Here, we do not introduce efficiency parameters, instead we test whether a change in mixing prescription alone could remedy the discrepancies.

### 3.4.1 Impact of initial mass, rotation, and mixing prescription

Fig. 5 shows a reference model (dashed line) with $M_{ini} = 18 M_{\odot}$ and $v_{ini} = 300 \text{ km s}^{-1}$, using the D22 scheme (see above) and assuming a 600 G polar field strength constant in time. The colour-coding scales with the stellar age. When decreasing the initial mass, only modest differences are seen (primarily, a shift in the effective temperature to lower values). Likewise, when the initial rotational velocity is increased (by 100 km s$^{-1}$ here) compared to the reference model, the impact remains small. Changing the mixing prescription, however, leads to a large difference. The model with the D11 scheme allows for a far more efficient chemical mixing, and due to the increased mean molecular weight on average inside the star, a more vertical evolution on the Kiel diagram.

From these model computations, we can conclude that within observational and modelling uncertainty the Kiel diagram predicts an age of at most 6 Myr.<fnref n="9" /> This is at odds with the model ages predicted from the Hunter-P diagram where we encounter the same problem as before. Namely, sufficient nitrogen enrichment and low spin velocity are not obtained with a self-consistent solution when using standard formulas for chemical mixing and magnetic braking. This is consistent with findings of Meynet et al. 2011 and Paper I, where

the model predictions typically yield notable surface enrichment in the second half of the main sequence only.

## 4 DISCUSSION

We have shown that *standard* single-star evolutionary models cannot simultaneously reproduce the observed log g, $T_{eff}$, N/H, and $P_{rot}$ of τ Sco.

Single-star models matching the observations on the Kiel diagram (fitting log g and $T_{eff}$) can either reproduce the N/H ratio but not the observed rotation, or, inversely can reproduce the observed rotation but not the surface enrichment. This is precisely the same problem encountered by the post-merger models of Schneider et al. (2020).

Possible resolutions, in the frame of single-star models, is to consider a more efficient chemical mixing and/or a more efficient magnetic braking. The latter one hypothetically could be replaced by a magnetic field decay scenario. Overall, we find that in order to reconcile all observables of τ Sco, rather extreme assumptions are necessary, which would not be compatible with observations of other stars.

### 4.1 Chemical mixing

Rotational mixing in massive stars is presently one of the most uncertain processes influencing stellar evolution. A large number of observations show discrepancies with current model predictions (e.g. Hunter et al. 2009; Martins et al. 2017; Cazorla et al. 2017b; Markova, Puls & Langer 2018).

Potential sources of these discrepancies are likely many-fold. Importantly, stellar evolution model computations remain restricted to a 1D treatment. Therefore, necessarily, scaling factors are introduced to model chemical mixing (Pinsonneault et al. 1989) and some mixing prescriptions are derived based on order-of-magnitude estimates (Brott et al. 2011). Even more, the relevant physical processes and their interactions are not yet fully understood (e.g. Maeder 2009) and some processes are not yet ubiquitously modelled and included in the computations, such as mixing by internal gravity waves (Decressin

<footnote n="9">⁹The values obtained from Simón-Díaz et al. (2006) may allow for a somewhat higher age, although the values we derive in this study would point to a lower age and perhaps very slightly lower initial mass than 17 $M_{\odot}$, based on the GENEC models.</footnote>
et al. 2009; Mathis et al. 2013; Rogers et al. 2013; Aerts & Rogers 2015; Rogers & McElwaine 2017; Bowman et al. 2019a,b, 2020).

Despite the notable uncertainties, a drastic increase in the efficiency of the present-day prescriptions of rotational mixing by a factor of 3 (in the initially fast-rotator MESA model) or a factor of 10 (in the initially moderate-rotator MESA model) seem excessive in a single-star model. We estimate that (along with an increase in magnetic braking efficiency) a similar increase in chemical mixing would also be required in the ‘D11’ GENEC model to obtain a self-consistent solution.

Let us recall here that the chemical mixing efficiency in rotating stellar evolution models is calibrated such that with an initial mean rotation rate, they reproduce the mean surface nitrogen abundances of B-type stars at the end of the main sequence phase (e.g. Brott et al. 2011). If we change this mixing efficiency for the purpose of fitting the data of τ Sco, then an appropriate physical reason should be given. At the moment, we are unaware of a physical cause which could be invoked in the single-star channel. Nevertheless, adopting a factor of 10 increase is not unprecedented in evolutionary modelling, for example, Aguilera-Dena et al. (2020) use such an increased mixing efficiency in their approach.

Interestingly, τ Sco stands out from the sample of magnetic massive stars as being a putative blue straggler. Blue stragglers are often associated with stellar mergers or quasi-chemically homogeneously evolving stars. Although quasi-chemically homogeneous evolution via long-term rapid rotation may help explain the surface enrichment of τ Sco, it is unclear if at all the star could suddenly become a very slow rotator as observed, therefore this evolutionary channel does not seem very favourable. Magnetic OB stars are typically found to be consistent with a usual ‘redward’ evolution after their initial spin-down (Paper I and Paper II). Furthermore, being a magnetic blue straggler poses the interesting question whether all blue stragglers would have a detectable magnetic field. Schneider et al. (2016) suggested that since the merger rate may be higher in blue stragglers, consequently the incidence rate of magnetism may also be higher. However, observational efforts dedicated to hot stars have not found such hints yet (Mathys 1988; Grunhut et al. 2017).

Morel et al. (2008), Aerts et al. (2014), and Martins et al. (2012), Martins et al. (2015) analysed the nitrogen enrichment of magnetic OB stars and identified a number of cases with nitrogen excess (including τ Sco). Presently, it remains elusive why some stars show excess and others do not. In Paper I, we proposed that the observable nitrogen enrichment of magnetic stars largely depends on presently unconstrained mixing processes (and their efficiencies), therefore the outcome is a multivariate function of a number of parameters (Aerts et al. 2014; Maeder et al. 2014), and not only a function of the magnetic field strength.

## 4.2 Magnetic field evolution and spin-down

Regarding its magnetic field characteristics, τ Sco clearly stands out of the known non-chemically peculiar magnetic B-type stars (this sample is discussed by Shultz et al. 2019b and confronted with evolutionary models in Paper II). The complexity of τ Sco’s magnetic field is unique. In general, a dipole-dominated model allows for reproducing magnetic field measurements (even if other higher order harmonics are present), although a number of cases indeed point to a quadruple-dominated geometry (Shultz et al. 2018, 2019a). Recently, David-Uraz et al. (2021) found that NGC 1624 - 2’s magnetic field is more complex than the typically assumed pure dipole. This evidence suggests that deviations from the pure dipole geometry may not be uncommon but require extensive monitoring to identify it.

The evolution of surface fossil magnetic fields remains largely uncertain. Several observational studies are consistent with magnetic flux conservation, whereas other empirical evidence points to a more rapid decline in field strength over time, suggesting a flux decay scenario (see Paper I and Paper II, and references therein).

From a theoretical standpoint, fossil fields are expected to slowly dissipate on an Ohmic time-scale, given by the induction equation of non-ideal MHD:

$$ \eta \nabla^2 \mathbf{B} = \frac{\partial \mathbf{B}}{\partial t}, \qquad (10) $$

where η is the magnetic diffusivity (assumed to be constant in space), B is the magnetic field vector, t is the time, and we assume that the fluid velocity is zero, accounting for an equilibrium fossil magnetic field (Braithwaite & Spruit 2017).

For a fully ionized plasma, the diffusion time-scale can be approximated as the ratio of the square of a characteristic length scale and the magnetic diffusivity:

$$ t_{\text{diff}} \sim \frac{R^2}{\eta}. \qquad (11) $$

For the Sun and low-mass stars, this formula typically leads to estimates of diffusion time-scales longer than the main sequence lifetime (e.g. Cowling 1945; Feiden & Chaboyer 2014). Here, we refrain from providing an estimate for more massive stars (and τ Sco in particular) since the uncertainty in an appropriate magnetic diffusivity is several orders of magnitudes (e.g. Charbonneau & MacGregor 2001). Since the outward diffusion of the magnetic flux is a possible mechanism to explain why flux decay may need to be invoked as a field evolution scenario, in this case, we would require the diffusion time-scale to be shorter than the ~12 Myr main sequence lifetime of our models. In conclusion, the scenario of magnetic flux decay still requires rigorous theoretical considerations (Braithwaite & Spruit 2017).

The magnetic braking formula by ud-Doula et al. (2009) is robust but relies on important assumptions such as the field geometry and alignment with the rotation axis. In principle, the complexity of τ Sco’s magnetic field geometry is expected to decrease, and not increase, the spin-down efficiency since the lowest order harmonic (i.e. the dipole component) has the most relevant contribution to spin-down. Nevertheless, the impact of a strong azimuthal field remains uncharacterized. Recent simulations (ud Doula et al., in preparation) show that oblique rotation does not significantly affect the spin-down scaling either.<fnref n="10" />

## 4.3 Alternative evolutionary scenarios

Although the single-star and the merger scenario have been suggested to explain the observables of τ Sco, nature may potentially offer different ways to explain them.

Two main channels might be considered. i) Higher mixing efficiency may be induced by a close companion star. Recent works have indeed shown that tidally induced mixing can supersede the efficiency compared to a single star (e.g. Song et al. 2013). This would require the presence of a main sequence companion star (see below). ii) The surface enrichment of τ Sco might have been caused by an earlier mass-transfer event from a more massive companion

<footnote n="10">¹⁰Deviations in the rate of angular momentum loss by approximately up to 30 per cent can be reached for oblique rotators (ud Doula, private communication). However, these deviations are expected to decrease the rate of angular momentum loss.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The Si III 4553 line in the spectra of τ Sco taken over 8 d with the CHFT/ESPaDOnS instrument.</figcaption></figure>

(but without a merger event), leaving behind a present-day helium star or a compact object. Since the proper motion and radial velocity of τ Sco are low, it is unlikely that it could have been ejected from such a system as a run-away star during the supernova explosion creating the compact remnant. However, in this hypothetical scenario, the system's orbit may have widened sufficiently that it is not a close binary system anymore.

In order to check the possibility that τ Sco is in a close binary system (with an orbital period of less than 20 d), we retrieved CFHT/ESPaDOnS spectra<fnref n="11" /> taken in between 2005 June 19 and 26. Fig. 6 displays these spectra around the Si III 4553 line. The radial velocity is remarkably stable over a week, with variations estimated to be less than 1 km s⁻¹. Other lines show the same behaviour. We also checked that no variability is present in spectra taken within the same night, a few hours apart. This hints at the absence of a massive companion on a short period orbit. Hence, effects of close binarity on the present-day properties and evolution of τ Sco are not found. The probability that we would view such a binary system at an inclination angle low enough to accommodate a < 1 km s⁻¹ projected orbital velocity is extremely low.

There is no clear consensus whether a long-period (> 20 d) companion is present or not. The study of Rizzuto et al. (2013) identifies a possible companion with a mass of roughly 1/3 that of the primary at an orbital separation of 2.84 AU; however, the literature is not conclusive about such a possible companion (Eggleton & Tokovinin 2008; Pecaut, Mamajek & Bubar 2012; Lafrenière et al. 2014; Grellmann et al. 2015).

While it is difficult to observationally identify a main sequence companion on a long-period orbit, it remains to be worked out if it could produce a sufficient boost in chemical mixing, providing a viable alternative explanation for τ Sco's observables.

## 4.4 Testing the single star versus the merger scenario

The major difficulty to evaluate the evolutionary history of τ Sco is that the presently considered scenarios (single star or merger) both lead to similar evolutionary tracks as the post-merger object is expected to evolve as a look-alike single star (Schneider et al. 2020). However, a few important points can be addressed. i) It is expected that the merger event leads to an abrupt mass removal from the system, and Glebbek et al. 2013 calculate that this mass-shedding could amount to about 10 per cent of the total mass of the system in the merger of equal mass stars. Such mass-shedding may be linked to the transient phenomenon known as luminous red nova (Blagorodnova et al. 2017), ultimately leading to a young nebula around the star (Langer 2012; Schneider et al. 2019, 2020). A young nebula is not seen around τ Sco, however, an asymmetric bow-shock appears present in the 22 µm WISE observations (Gvaramadze et al. 2019). The origin of this bow-shock remains to be characterized. ii) Asteroseismology could provide strong evidence regarding the near core chemical composition and rotation profile of the star (Aerts, Mathis & Rogers 2019; Bowman 2020) and these predicted characteristics are distinguishable between the post-merger model of Schneider et al. (2020) and our single-star models.<fnref n="12" /> iii) The surface helium content is a strong indicator of evolutionary status and we now discuss this point further.

### 4.4.1 Helium abundance

Fig. 7 shows those of our evolutionary models that have some helium enrichment over the main sequence. The models with increased magnetic braking efficiency and the models with field decay are not shown as they rapidly brake the rotation and do not lead to any notable helium enrichment. We find that within 6 Myr, none of the models produce a surface helium abundance compatible with our measured helium abundance of $Y = 0.30 \pm 0.02$ (Section 2). Some models, particularly the ‘D11’ GENEC model, may reach sufficiently high values but only close to the TAMS.

In fact, in single stars a large amount of helium is only expected to be mixed to the surface at the end of the main sequence evolution, while, for example, the build-up of nitrogen is much more rapid.

The reason is that for nitrogen a strong gradient is built very early between the core (rich in nitrogen) and the envelope (with only the initial nitrogen), and the velocity of diffusion increases with the gradient of the considered element (equation 3 of Meynet, Maeder & Mowlavi 2004).

During the early evolution, the contrast between the helium abundance in the core and envelope is much weaker thus the diffusion is much less rapid. Eventually, the gradient in helium composition will increase and, consequently, will allow for some surface enrichment.

Although in MESA we adopted a commonly used formalism with a scaling factor $f_μ$ to mitigate the impact of composition gradients (Yoon et al. 2006; Brott et al. 2011), the rapid decrease of the rotational velocity means that these models cannot mix a sufficient amount of helium to the stellar surface. Of course, an important provision is that the effect of the magnetic field on chemical mixing still needs to be explored. Especially for τ Sco's complex magnetic field, this remains an open question which warrants further investigation.

<footnote n="12">¹²τ Sco was observed with the *Transiting Exoplanet Survey Satellite* (TESS). Unfortunately, from the *TESS* data, no clear non-radial oscillations are detectable which prevents seismic modelling (Buzasi, Aerts, Bowman, private communication).</footnote>

<footnote n="11">¹¹The spectra were collected from the Polarbase database of the CFHT/ESPaDOnS and TBL/NARVAL spectropolarimeters, accessible at: http://polarbase.irap.omp.eu/</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Surface helium abundance versus time. The same MESA models are shown as in Fig 2, and the top panels of Figs 3 and 4, from top left to lower middle. The lower right panel shows the same GENEC models as in Fig. 5, including our measured helium abundance with blue (however, see Section 1 and Appendix B for a range of measurements by other authors). Note the different scales on the vertical axes.</figcaption></figure>

In summary, when adopting a baseline of $Y_{ini} = 0.266$ in our models, the measured helium abundance is at odds with single-star model predictions at a $1\sigma$ level. This clue suggests that invoking binarity may be necessary.

## 4.5 Age and abundance of Upper Sco and τ Sco's membership

The proper motion and radial velocity of τ Sco are broadly consistent with other stars identified as members of the Upper Sco association. The Gaia DR2 parallax places τ Sco near the edge of the association with a distance of 195 pc, compared to the association's mean distance of 145 pc (Gaia Collaboration 2018).

The age of the Upper Sco association is constrained by main-sequence turn-off fitting at about 11 Myr by several authors (Sartori, Lépine & Dias 2003; Pecaut et al. 2012; Feiden 2016; David et al. 2019). Some authors find a much lower value of about 5 Myr (de Geus, de Zeeuw & Lub 1989; Blaauw 1991; Preibisch et al. 2002; Dahm & Carpenter 2009). The differences in reported ages are attributed to modelling uncertainties and observational bias but not to an actual age spread in the association (Donaldson et al. 2017; Fang, Herczeg & Rizzuto 2017).

However, it is important to realize that Upper Sco's age determination plays a vital role in determining the star's evolutionary history. If Upper Sco's age is indeed about 11 Myr, then a rejuvenating mechanism (disfavouring single-star evolution) must be invoked since the apparent age of the star is less than 6 Myr. A detailed

kinematic study of τ Sco is much needed and tracing back its birthplace could help to better constrain its evolution.

Previous studies have found that the elemental abundance of Upper Sco is close to solar (Mamajek et al. 2013). This lends support to choosing the solar metallicity (and nitrogen abundance) as the baseline to model the chemical enrichment. Therefore, we do not expect that the observed nitrogen excess of τ Sco could somehow be due to an initially higher nitrogen abundance. We further test the impact of choosing various initial abundances in Appendix A.

In summary, there is no evidence suggesting that τ Sco does not belong to Upper Sco. Likewise, there is no evidence that an unusual star formation could have led to a late formation of τ Sco. However, the age of Upper Sco plays a decisive role in establishing τ Sco's evolutionary history. Therefore further work on constraining Upper Sco's age and chemical composition with similar modelling assumptions as taken here would be a logical next step.

## 5 CONCLUSIONS

τ Scorpii is a very special and unusual magnetic star. In this work, we address the question if it is possible at all to reconcile the observable properties of τ Sco with single-star evolution.

We perform atmospheric modelling and determine the CNO and helium abundances. The measured nitrogen enrichment (along with carbon depletion) implies that the stellar surface is mixed with core-processed material. The derived helium abundance is close to the baseline value but might hint a slight excess. This rules out that τ Sco
could be an evolved helium star. Equal-mass mergers are expected to lead to a substantial increase in the surface helium abundance (GlebbEEK et al. 2013; Schneider et al. 2016), whereas in single stars the surface helium abundance should only increase towards the end of the main sequence.

We present experiments with stellar evolution models from which we draw the following conclusions:

(i) Considering observational and modelling uncertainties, we estimate a maximum age of 6 Myr for τ Sco in the framework of the single-star scenario.

(ii) In order to produce surface enrichment by rotational mixing, an initial spin velocity higher than 150 km s⁻¹ is required. This disfavours that τ Sco could have been a very slow rotator at the ZAMS if it was a single star.

(iii) We can only reconcile the spin-down of the stellar surface with the observed 41 day rotational period (Donati et al. 2006) if we increase the magnetic braking efficiency by a factor of 10 compared to a standard prescription.

(iv) Alternatively, an increase in magnetic braking efficiency may not need to be invoked if one assumes an initially much stronger field (30 kG) which has rapidly decayed over time. There exists tentative evidence for magnetic flux decay from observations (e.g. Fossati et al. 2016; Shultz et al. 2019b), however, the decay rate remains to be quantified. Interestingly, the flux decay could potentially be caused by a companion star (Vidal et al. 2019).

(v) In a typical rapidly rotating massive star model, the observed nitrogen excess may be reached within 6 Myr. However, due to the necessary rapid spin-down by magnetic braking, in both of the above scenarios one needs to invoke a subsequent increase in the efficiency of rotational mixing by at least a factor of 3.

Presently, neither the single-star nor the post-merger models can satisfactorily explain the simultaneous slow rotation and nitrogen excess of τ Sco. Neither can the rotation be explained with a present-day sub-kG magnetic field strength, unless a seemingly extreme magnetic field evolution is invoked. Although there are ways to reconcile these discrepancies in both channels, binary evolution may provide a more versatile scenario to achieve this. We find no evidence for a present-day close companion star and the literature remains inconclusive regarding a possible companion on a long-period orbit.

A number of follow-up observations and modelling efforts can help shed more light on the nature of τ Sco.

(i) Extended series of Stokes Q, U profiles over the rotation period of the star would allow to measure the magnetic field modulus, aiding to constrain the surface field geometry.

(ii) A detailed observational investigation of a possible companion object (a long-period main sequence star, or a compact remnant) would be of great value to assess binary-evolution scenarios which may not lead to a stellar merger.

(iii) A theoretical explanation of tidally induced mixing from a highly eccentric, long-period companion or whether previous binary mass transfer could produce self-consistent solutions should be scrutinized in the context of τ Sco.

Finally, we emphasize that τ Sco does not seem to be a representative magnetic B-type star and its magnetic field is much more complex than those inferred for other B-type stars. Consequently, whichever evolutionary channel produced this special star, may not be able to provide a general explanation for the origin of fossil magnetic fields in massive stars. However, it is an intriguing prospect if unusually complex magnetic field geometries of OBA stars may be linked to stellar multiplicity.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for a constructive report which helped improve the manuscript. We thank John Hillier for making CMFGEN available to the community and for constant help with it. We thank the MESA developers for making their code publicly available. We appreciate discussions with Fabian Schneider, Henk Spruit, Conny Aerts, Derek Buzasi, Dominic Bowman, Cole Johnston, Jo Puls, Paco Najarro, Asif ud-Doula, Difeng Guo, and Sam Geen. GM acknowledges support from the Swiss National Science Foundation (project number 200020-172505). GM has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 833925, project STAREX). ADU gratefully acknowledges the support of the Natural Science and Engineering Research Council of Canada (NSERC). This work is supported by NASA under award number 80GSFC17M0002. This work made use of the Polarbase database (developed and maintained by CNRS/INSU, Observatoire Midi-Pyrénées and Université Toulouse III). This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.

## DATA AVAILABILITY

A full reproduction package is shared on zenodo in accordance with the Research Data Management plan of the Anton Panekoek Institute for Astronomy at the University of Amsterdam: 10.5281/zenodo.4633408

## REFERENCES

Aerts C., Rogers T. M., 2015, ApJ, 806, L33

Aerts C., Molenberghs G., Kenward M. G., Neiner C., 2014, ApJ, 781, 88

Aerts C., Mathis S., Rogers T. M., 2019, ARA&A, 57, 35

Aguilera-Dena D. R., Langer N., Antoniadis J., Müller B., 2020, ApJ, 901, 114

Alfvén H., 1942, Nature, 150, 405

Aller L. H., Faulkner D. J., Norton R. H., 1966, ApJ, 144, 1073

Anders E., Grevesse N., 1989, Geochim. Cosmochim. Acta, 53, 197

Asplund M., Grevesse N., Sauval A. J., 2005, in Barnes T. G. III, Bash F. N., eds, ASP Conf. Ser. Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, Astron. Soc. Pac., San Francisco, p. 25

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Blaauw A., 1991, in Lada C. J., Kylafis N. D., eds, NATO Advanced Study Institute (ASI) Series C Vol. 342, The Physics of Star Formation and Early Stellar Evolution. p. 125

Blagorodnova N. et al., 2017, ApJ, 834, 107

Bowman D. M. et al., 2019a, Nat. Astron., 3, 760

Bowman D. M. et al., 2019b, A&A, 621, A135

Bowman D. M., 2020, Frontiers in Astron. and Space Sci., 7, 70

Bowman D. M., Burssens S., Simón-Díaz S., Edelmann P. V. F., Rogers T. M., Horst L., Röpke F. K., Aerts C., 2020, A&A, 640, A36

Bragança G. A., Daflon S., Cunha K., Bensby T., Oey M. S., Walth G., 2012, AJ, 144, 130

Braithwaite J., 2008, MNRAS, 386, 1947

Braithwaite J., Spruit H. C., 2017, R. Soc. Open Sci., 4, 160271

Brott I. et al., 2011, A&A, 530, A115

Cazorla C., Morel T., Nazé Y., Rauw G., Semaan T., Daflon S., Oey M. S., 2017a, A&A, 603, A56

Cazorla C., Nazé Y., Morel T., Georgy C., Godart M., Langer N., 2017b, A&A, 604, A123

Charbonneau P., MacGregor K. B., 2001, ApJ, 559, 1094

Cohen D. H., Cassinelli J. P., Waldron W. L., 1997, ApJ, 488, 397

Cohen D. H., de Messières G. E., MacFarlane J. J., Miller N. A., Cassinelli J. P., Owocki S. P., Liedahl D. A., 2003, ApJ, 586, 495
Cowley C. R., Merritt D. R., 1987, ApJ, 321, 553

Cowling T. G., 1945, MNRAS, 105, 166

Dahm S. E., Carpenter J. M., 2009, AJ, 137, 4024

David T. J., Hillenbrand L. A., Gillen E., Cody A. M., Howell S. B., Isaacson H. T., Livingston J. H., 2019, ApJ, 872, 161

David-Uraz A. et al., 2019, MNRAS, 483, 2814

David-Uraz A., Petit V., Shultz M. E., Fullerton A. W., Erba C., Keszthelyi Z., Seadrow S., Wade G. A., 2021, MNRAS, 501, 2677

de Geus E. J., de Zeeuw P. T., Lub J., 1989, A&A, 216, 44

de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259

de Messieres G. E., Cardamone C., Cohen D. H., MacFarlane J. J., Owocki S. P., ud-Doula A., 2001, in American Astronomical Society Meeting Abstracts. p. 135.12

Decressin T., Mathis S., Palacios A., Siess L., Talon S., Charbonnel C., Zahn J. P., 2009, A&A, 495, 271

Donaldson J. K., Weinberger A. J., Gagné J., Boss A. P., Keiser S. A., 2017, ApJ, 850, 11

Donati J. F. et al., 2006, MNRAS, 370, 629

Donati J.-F., Landstreet J. D., 2009, ARA&A, 47, 333

Eddington A. S., 1925, Observatory, 48, 73

Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Ap&SS, 316, 43

Eggleton P. P., Tokovinin A. A., 2008, MNRAS, 389, 869

Ekström S. et al., 2012, A&A, 537, A146

Endal A. S., Sofia S., 1978, ApJ, 220, 279

Fang Q., Herczeg G. J., Rizzuto A., 2017, ApJ, 842, 123

Feiden G. A., 2016, A&A, 593, A99

Feiden G. A., Chaboyer B., 2014, ApJ, 789, 53

Ferrario L., Pringle J. E., Tout C. A., Wickramasinghe D. T., 2009, MNRAS, 400, L71

Fletcher C. L., 2018, PhD thesis, Department of Physics and Space Sciences. Florida Institute of Technology, Melbourne, Florida

Fletcher C. L., Petit V., Cohen D. H., Townsend R. H., Wade G. A., 2018, Contrib. Astron. Obs. Skaln. Pleso, 48, 144

Fossati L. et al., 2016, A&A, 592, A84

Fricke K., 1968, Z. Astrophys., 68, 317

Gaia Collaboration, 2018, A&A, 616, A1

Georgy C. et al., 2013, A&A, 558, A103

Georgy C., Meynet G., Ekström S., Wade G. A., Petit V., Keszthelyi Z., Hirschi R., 2017, A&A, 599, L5

Glebbeek E., Gaburov E., Portegies Zwart S., Pols O. R., 2013, MNRAS, 434, 3497

Goldreich P., Schubert G., 1967, ApJ, 150, 571

Grellmann R., Ratzka T., Köhler R., Preibisch T., Mucciarelli P., 2015, A&A, 578, A84

Grevesse N., Noels A., Sauval A. J., 1996, in Holt S. S., Sonneborn G., eds, ASP Conf. Ser. Vol. 99, Cosmic Abundances, Astron. Soc. Pac., San Francisco, p. 117

Groh J. H. et al., 2019, A&A, 627, A24

Grunhut J. H. et al., 2017, MNRAS, 465, 2432

Gvaramadze V. V., Maryeva O. V., Kniazev A. Y., Alexashov D. B., Castro N., Langer N., Katkov I. Y., 2019, MNRAS, 482, 4408

Hardorp J., Scholz M., 1970, ApJS, 19, 193

Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368

Herwig F., 2000, A&A, 360, 952

Hillier D. J., Miller D. L., 1998, ApJ, 496, 407

Howk J. C., Cassinelli J. P., Bjorkman J. E., Lamers H. J. G. L. M., 2000, ApJ, 534, 348

Hunter I. et al., 2009, A&A, 496, 841

Høiland K., 1941, Avhandliger Norske Videnskaps-Akademi i Oslo, I math.-naturv. Klasse, 11, 1

Ignace R., Oskinova L. M., Jardine M., Cassinelli J. P., Cohen D. H., Donati J. F., Townsend R. H. D., ud-Doula A., 2010, ApJ, 721, 1412

Kaiser E. A., Hirschi R., Arnett W. D., Georgy C., Scott L. J. A., Cristini A., 2020, MNRAS, 496, 1967

Keszthelyi Z. et al., 2020, MNRAS, 493, 518 (Paper II)

Keszthelyi Z., Puls J., Wade G. A., 2017, A&A, 598, A4

Keszthelyi Z., Meynet G., Georgy C., Wade G. A., Petit V., David-Uraz A., 2019, MNRAS, 485, 5843 (Paper I)

Khan S. A., Shulyak D. V., 2006, A&A, 448, 1153

Kilian J., 1992, A&A, 262, 171

Kilian J., Becker S. R., Gehren T., Nissen P. E., 1991, A&A, 244, 419

Kochukhov O., Wade G. A., 2016, A&A, 586, A30

Krtička J., 2018, A&A, 620, A176

Lafrenière D., Jayawardhana R., van Kerkwijk M. H., Brandeker A., Janson M., 2014, ApJ, 785, 47

Lamers H. J. G. L. M., Rogerson J. B. J., 1978, A&A, 66, 417

Langer N., 2012, ARA&A, 50, 107

Lodders K., 2003, ApJ, 591, 1220

Macfarlane J. J., Cassinelli J. P., 1989, ApJ, 347, 1090

Maeder A., 1997, A&A, 321, 134

Maeder A., 2003, A&A, 399, 263

Maeder A., 2009, Physics, Formation and Evolution of Rotating Stars, Astronomy and Astrophysics Library. Springer, Berlin Heidelberg

Maeder A., Przybilla N., Nieva M.-F., Georgy C., Meynet G., Ekström S., Eggenberger P., 2014, A&A, 565, A39

Mamajek E. E., Pecaut M. J., Nguyen D. C., Bubar E. J., 2013, in Protostars and Planets VI Posters

Markova N., Puls J., Langer N., 2018, A&A, 613, A12

Martins F. et al., 2015, A&A, 575, A34

Martins F., Plez B., 2006, A&A, 457, 637

Martins F., Escolano C., Wade G. A., Donati J. F., Bouret J. C., Mimes Collaboration, 2012, A&A, 538, A29

Martins F., Simón-Díaz S., Barbá R. H., Gamen R. C., Ekström S., A&A, 599, A30

Mathis S., Decressin T., Eggenberger P., Charbonnel C., 2013, A&A, 558, A11

Mathys G., 1988, The Messenger, 53, 39

Meakin C. A., Arnett D., 2007, Astrophys. J., 667, 448

Mewe R., Raassen A. J. J., Cassinelli J. P., van der Hucht K. A., Miller N. A., Güdel M., 2003, A&A, 398, 203

Meynet G., Maeder A., Mowlavi N., 2004, A&A, 416, 1023

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

Meynet G., Ekstrom S., Maeder A., Eggenberger P., Saio H., Chomienne V., Haemmerlé L., 2013, in Goupil M., Belkacem K., Neiner C., Lignières F., Green J. J., eds, Lecture Notes in Physics, Vol. 865. Springer Verlag, Berlin, p. 3

Mokiem M. R., de Koter A., Puls J., Herrero A., Najarro F., Villamariz M. R., 2005, A&A, 441, 711

Moravveji E., Aerts C., Pápics P. I., Triana S. A., Vandoren B., 2015, A&A, 580, A27

Morel T., Hubrig S., Briquet M., 2008, A&A, 481, 453

Murphy L. J. et al., 2021, MNRAS, 501, 2745

Nazé Y., Petit V., Rinbrand M., Cohen D., Owocki S., ud-Doula A., Wade G. A., 2014, ApJS, 215, 10

Nieva M. F., Przybilla N., 2012a, A&A, 539, A143

Nieva M.-F., Przybilla N., 2012b, A&A, 539, A143

Nieva M.-F., Przybilla N., 2014, A&A, 566, A7

Owocki S. P., ud-Doula A., Sundqvist J. O., Petit V., Cohen D. H., Townsend R. H. D., 2016, MNRAS, 462, 3830

Pápics P. I. et al., 2017, A&A, 598, A74

Paxton B. et al., 2013, ApJS, 208, 4

Paxton B. et al., 2015, ApJS, 220, 15

Paxton B. et al., 2018, ApJS, 234, 34

Paxton B. et al., 2019, ApJS, 243, 10

Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3

Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ, 746, 154

Peters G. J., Polidan R. S., 1985, in Hayes D. S., Pasinetti L. E., Philip A. G. D., eds, IAU Symp. Vol. 111, Calibration of Fundamental Stellar Quantities. p. 417

Petit V. et al., 2013, MNRAS, 429, 398

Petit V. et al., 2017, MNRAS, 466, 1052

Pinsonneault M. H., Kawaler S. D., Sofia S., Demarque P., 1989, ApJ, 338, 424
Preibisch T., Brown A. G. A., Bridges T., Guenther E., Zinnecker H., 2002, AJ, 124, 404

Przybilla N., Butler K., 2004, ApJ, 609, 1181

Przybilla N., Nieva M.-F., Butler K., 2008, ApJ, 688, L103

Przybilla N., Firnstein M., Nieva M. F., Meynet G., Maeder A., 2010, A&A, 517, A38

Puls J., Urbaneja M. A., Venero R., Repolust T., Springmann U., Jokuthy A., Mokiem M. R., 2005, A&A, 435, 669

Repolust T., Puls J., Hanson M. M., Kudritzki R. P., Mokiem M. R., 2005, A&A, 440, 261

Rizzuto A. C. et al., 2013, MNRAS, 436, 1694

Rogers F. J., Iglesias C. A., 1992, ApJS, 79, 507

Rogers T. M., McElwaine J. N., 2017, ApJ, 848, L1

Rogers T. M., Lin D. N. C., McElwaine J. N., Lau H. H. B., 2013, ApJ, 772, 21

Rogerson J. B. J., Ewell M. W. J., 1985, ApJS, 58, 265

Sartori M. J., Lépine J. R. D., Dias W. S., 2003, A&A, 404, 913

Schneider F. R. N., Podsiadlowski P., Langer N., Castro N., Fossati L., 2016, MNRAS, 457, 2355

Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., Springel V., 2019, Nature, 574, 211

Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., 2020, MNRAS, 495, 2796

Shultz M. E. et al., 2018, MNRAS, 475, 5144

Shultz M. E. et al., 2019b, MNRAS, 490, 274

Shultz M. E., Wade G. A., Rivinius T., Alecian E., Neiner C., Petit V., Wisniewski J. P., 2019a, MNRAS, 485, 1508

Simón-Díaz S., Herrero A., 2007, A&A, 468, 1063

Simón-Díaz S., Herrero A., Esteban C., Najarro F., 2006, A&A, 448, 351

Simón-Díaz S., Godart M., Castro N., Herrero A., Aerts C., Puls J., Telting J., Grassitelli L., 2017, A&A, 597, A22

Snow T. P., Lamers H. J. G. L. M., Lindholm D. M., Odell A. P., 1994, ApJS, 95, 163

Solberg H., 1936, Astrophys. Nor., 1, 237

Song H. F., Maeder A., Meynet G., Huang R. Q., Ekström S., Granada A., 2013, A&A, 556, A100

Spruit H. C., 2002, A&A, 381, 923

Stibbs D. W. N., 1950, MNRAS, 110, 395

Struve O., Dunham T. J., 1933, ApJ, 77, 321

Sweet P. A., 1950, MNRAS, 110, 548

Talon S., Zahn J.-P., 1997, A&A, 317, 749

Tayler R. J., 1973, MNRAS, 161, 365

Traving G., 1955, Z. Astrophys., 36, 1

ud-Doula A., 2017, Astron. Nachr., 338, 944

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

Unsöld A., 1942, Z. Astrophys., 21, 1

Vidal J., Cébron D., ud-Doula A., Alecian E., 2019, A&A, 629, A142

Vink J. S., de Koter A., Lamers H. J. G. L. M., 2000, A&A, 362, 295

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

Walborn N. R., Panek R. J., 1984, ApJ, 286, 718

Waters L. B. F. M., Marlborough J. M., Geballe T. R., Oosterbroek T., Zaal P., 1993, A&A, 272, L9

Weber E. J., Davis L. Jr, 1967, ApJ, 148, 217

Wolff S. C., Heasley J. N., 1985, ApJ, 292, 589

Yoon S. C., Langer N., Norman C., 2006, A&A, 460, 199

Zaal P. A., de Koter A., Waters L. B. F. M., Marlborough J. M., Geballe T. R., Oliveira J. M., Foing B. H., 1999, A&A, 349, 573

Zahn J.-P., 1992, A&A, 265, 115

Zorec J., Cidale L., Arias M. L., Frémat Y., Muratore M. F., Torres A. F., Martayan C., 2009, A&A, 501, 297

## APPENDIX A: FURTHER TESTS

There are a number of modelling assumptions which may affect the quantitative results. On the top panels of Fig. A1, we first test how adopting different mixtures of metals, that is, altering the initial baseline value of nitrogen, modifies our findings. To this extent, we use a reference MESA model ($M_{ini} = 18 M_{\odot}$, $v_{ini} = 350$ km s$^{-1}$, $B_{p} = 600$ G, $f_{c} = 0.033$, $f_{MB} = 1$) and change the initial abundances which are now adopted from Anders & Grevesse (1989), Asplund et al. (2005), Asplund et al. (2009) and the latter one updated by Przybilla et al. (2010) and Nieva & Przybilla (2012a). (In MESA, this correspond to choosing `initial_zfracs=1, 5, 6, and 8, respectively.) We find that the initial baseline value does not play a significant role.

The lower panels of Fig. A1 show the results of adopting different overshooting values. The exponential overshooting values tested here are $f_{ov} = 0.001, 0.020, 0.034, 0.090$ together with $f_{0} = 0.006$.

The lowest value formally leads to decreasing the convective core size since $f_{ov}$ is applied from the nominal convective core boundary minus a distance $f_{0}$ of the local pressure scale height. The latter values would correspond to extending the convective core by roughly $\alpha_{ov} = 0.1, 0.25$, and 0.8 of the local pressure scale height in non-rotating models.

An extreme value of overshooting as shown with $f_{ov} = 0.090$ could, in principle, mitigate the need for a much more efficient envelope mixing (which is controlled by $f_{c}$ in our approach). However, overshooting only mixes core material to the base of the envelope. Therefore, if envelope mixing is inefficient, even an extreme overshooting value will fail to produce a model which accounts for surface enrichment.
<figure><img src="image_10.png" /><figcaption>**Figure A1.** Upper panels: test with different initial metal fractions adopted from various authors. Lower panels: test with different exponential overshooting parameters. The reference MESA model in both cases has $M_{\text{ini}} = 18 \text{ M}_\odot$, $v_{\text{ini}} = 350 \text{ km s}^{-1}$, $B_p = 600 \text{ G}$, $f_c = 0.033$, $f_{\text{MB}} = 1$, $f_{\text{ov}} = 0.034$, and uses the Asplund et al. (2009) mixture of metals.</figcaption></figure>

## APPENDIX B: PREVIOUS DETERMINATIONS OF ATMOSPHERIC, ROTATIONAL, WIND, AND MAGNETIC PROPERTIES

In this section, we report a collection of previous determinations of stellar properties that are found in the recent (in the last ~20 yr) literature. This modern list (see Table B1) is not meant to be exhaustive, but rather serves to illustrate typical ranges and uncertainties existing in the determination of these parameters, and that can affect the interpretation of our evolutionary modelling results. The listed works use a range of photometric and spectroscopic methods to determine parameters, and therefore provide a good sense of the variance across methods.

We note that some authors contrast different methods, which yield, for example, $T_{\text{eff}} = 29.6 \pm 0.15 \text{ kK}$ in one case and $T_{\text{eff}} = 31.4 \pm 2 \text{ kK}$ in the other by Zorec et al. (2009).

The wind and magnetic properties have a significant impact on the modelling presented here, as their interplay determines how quickly the stellar surface spins down. However, certain complications arise specifically as a result of this interaction. Any empirical determina-

tion of the wind mass-loss rate $\dot{M}$ based on spherically symmetric models will not properly account for the structure and anisotropy of the outflow induced by the field, so theoretical determinations of the *wind-feeding rate* (see e.g. discussion by David-Uraz et al. 2019 on this topic), such as that adopted in this study, might be more straightforward. Previous determinations in the literature provide a range of values. On the lower end, Repolst et al. (2005) derive either $\dot{M} = 2 \times 10^{-8} \text{ M}_\odot \text{ yr}^{-1}$ using FASTWIND (Puls et al. 2005) to fit H- and K-band spectra, or values varying between $9 \times 10^{-9}$ and $2 \times 10^{-8} \text{ M}_\odot \text{ yr}^{-1}$ using optical spectra, for parameters derived respectively by Kilian et al. (1991) and Przybilla & Butler (2004). Similarly, on the higher end, Mokiem et al. (2005) derive $\dot{M} = 6.14 \times 10^{-8} \text{ M}_\odot \text{ yr}^{-1}$ by coupling optical spectra with their fitting algorithm, also using FASTWIND models. These values are roughly consistent with the adopted mass-loss rates in the evolutionary models.

Furthermore, $\tau$ Sco possesses a complex magnetic field, which cannot be simply parametrized using the usual conventions for dipolar fields. The most comprehensive description of $\tau$ Sco's field is provided by Donati et al. (2006). Their longitudinal field
<table><caption>**Table B1.** Stellar parameters determined in previous studies: spectroscopic distance, effective temperature, surface gravity, luminosity, projected rotational velocity, microturbulent velocity, macroturbulent velocity, and He and CNO abundances in number fractions.</caption><thead><tr><th>Sources</th><th>M05</th><th>P08</th><th>M08</th><th>P12</th><th>B12</th><th>NP12</th><th>C17</th><th>SD17</th></tr></thead><tbody><tr><td>Distance (pc)</td><td>–</td><td>152 ± 20</td><td>–</td><td>–</td><td>–</td><td>143 ± 9</td><td>–</td><td>–</td></tr><tr><td>T<sub>eff</sub> (kK)</td><td>31.9<sup>+0.5</sup><sub>-0.8</sub></td><td>32 ± 0.3</td><td>31.5</td><td>29.9 ± 5</td><td>~29.9</td><td>32 ± 0.3</td><td>31.2 ± 1.5</td><td>~32.4 ± 0.16</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.15<sup>+0.09</sup><sub>-0.14</sub></td><td>4.30 ± 0.05</td><td>4.05 ± 0.15</td><td>–</td><td>–</td><td>4.30 ± 0.05</td><td>4.3 ± 0.15</td><td>–</td></tr><tr><td>log (L/L<sub>☉</sub>)</td><td>4.39 ± 0.09</td><td>–</td><td>–</td><td>4.31 ± 0.16</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>v sin i (km s<sup>-1</sup>)</td><td>~ 5</td><td>4 ± 4</td><td>8 ± 2</td><td>10 ± 2</td><td>3 ± 2</td><td>4 ± 1</td><td>8</td><td>7 and 9</td></tr><tr><td>v<sub>mic</sub> (km s<sup>-1</sup>)</td><td>–</td><td>5 ± 1</td><td>2 ± 2</td><td>–</td><td>–</td><td>5 ± 1</td><td>–</td><td>–</td></tr><tr><td>v<sub>mac</sub> (km s<sup>-1</sup>)</td><td>–</td><td>4 ± 4</td><td>–</td><td>–</td><td>–</td><td>4 ± 1</td><td>–</td><td>10</td></tr><tr><td>He/H</td><td>0.12<sup>+0.04</sup><sub>-0.02</sub></td><td>0.098 ± 0.01</td><td>–</td><td>–</td><td>–</td><td>–</td><td>0.083 ± 0.025</td><td>–</td></tr><tr><td>log(C/H) + 12</td><td>–</td><td>8.30 ± 0.12</td><td>8.19 ± 0.14</td><td>–</td><td>–</td><td>8.30 ± 0.12</td><td>8.18 ± 0.12</td><td>–</td></tr><tr><td>log(N/H) + 12</td><td>–</td><td>8.16 ± 0.12</td><td>8.15 ± 0.20</td><td>–</td><td>–</td><td>8.16 ± 0.12</td><td>7.90 ± 0.13</td><td>–</td></tr><tr><td>log(O/H) + 12</td><td>–</td><td>8.77 ± 0.08</td><td>8.62 ± 0.20</td><td>–</td><td>–</td><td>8.77 ± 0.08</td><td>8.39 ± 0.21</td><td>–</td></tr></tbody></table>

M05 = Mokiem et al. (2005); P08 = Przybilla, Nieva & Butler (2008); M08 = Morel et al. (2008); P12 = Pecaut et al. (2012); B12 = Bragança et al. (2012); NP12 = Nieva & Przybilla (2012b); C17 = Cazorla et al. (2017a); SD17 = Simón-Díaz et al. (2017)

measurements vary between −51.9 and + 87.8 G, and they also derive an average field modulus over the entire surface of the star of roughly 300 G. Nevertheless, the surface mapping of the field (performed using Zeeman Doppler Imaging) offers the best description of its geometry. Other measurements of this star's field are expressed in a variety of ways, including reporting the coefficients of a multi-harmonic fit to the longitudinal field curve (Shultz et al. 2018) obtain the following coefficients: $B_0 = 0.0094 \pm 0.0009$ kG, $B_1 = 0.049 \pm 0.001$ kG, $B_2 = 0.032 \pm 0.001$ kG, $B_3 = 0.019 \pm 0.001$ kG), and reporting the strength of the dipolar component of the field (Shultz et al. 2019b obtain $B_p = 0.31 \pm 0.02$ kG).

Finally, to extrapolate the magnetic field structure into the wind, Donati et al. (2006) derive a source radius<fnref n="13" /> of 2 $R_*$, based on

their adopted value of $\dot{M} = 2 \times 10^{-8} M_{\odot} \text{yr}^{-1}$. However, X-ray modelling using an arbitrary (non-dipolar) implementation of the Analytic Dynamical Magnetosphere model (Owocki et al. 2016) finds that a factor of 2 larger source radius (and hence a lower wind-feeding rate), leading to higher shock temperatures, is required to properly reproduce the features of $\tau$ Sco's X-ray spectrum (Fletcher 2018, Fletcher et al. 2018).

This paper has been typeset from a TeX/Latex file prepared by the author.

<footnote n="13">¹³Similar to the Alfvén radius for a dipolar field but applied to a more complex geometry, this is related to the radius of closure of the largest closed magnetic field loop.</footnote>