# A magnetic massive star has experienced a stellar merger

A. J. Frost, H. Sana, L. Mahy, G. Wade, J. Barron, J.-B. Le Bouquin, A.

Mérand, F. R. N. Schneider, T. Shenar, R. H. Barbá, et al.

► To cite this version:

A. J. Frost, H. Sana, L. Mahy, G. Wade, J. Barron, et al. A magnetic massive star has experienced a stellar merger. Science, 2024, 384, pp.214-217. 10.1126/science.adg7700. insu-04836883

HAL Id: insu-04836883

https://insu.hal.science/insu-04836883v1

Submitted on 18 Jul 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
A magnetic massive star has experienced a stellar
merger

A. J. Frost *¹,², H. Sana¹, L. Mahy³,¹, G. Wade⁴, J. Barron⁵,⁴, J.-B. Le Bouquin⁶, A. Mérand⁷, F. R. N. Schneider⁸,⁹, T. Shenar¹⁰,¹, R. H. Barbá†¹¹, D. M. Bowman¹²,¹, M. Fabry¹, A. Farhang¹³, P. Marchant¹, N. I. Morrell¹⁴ and J. V. Smoker²,¹⁵

¹Institute of Astronomy, KU Leuven, 3001 Leuven, Belgium

²European Southern Observatory, Santiago, Chile

³Royal Observatory of Belgium, B-1180 Brussels, Belgium

⁴Department of Physics & Space Science, Royal Military College of Canada, Kingston Ontario K7K 0C6, Canada

⁵Department of Physics, Engineering & Astronomy, Queen's University, Kingston Ontario K7L 3N6, Canada

⁶Université Grenoble Alpes, Centre national de la Recherche Scientifique, Institut de Planétologie et d'Astrophysique de Grenoble, F-38000 Grenoble, France

⁷European Southern Observatory Headquarters, 85748 Garching bei München, Germany

⁸Heidelberger Institut für Theoretische Studien, 69118 Heidelberg, Germany

⁹Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, 69120 Heidelberg, Germany

¹⁰The School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 69978, Israel

¹¹Departamento de Física y Astronomía, Universidad de la Serena, La Serena, Chile

¹²School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne NE1 7RU, UK

¹³School of Astronomy, Institute for Research in Fundamental Sciences, 19395-5531 Tehran, Iran

¹⁴Las Campanas Observatory, Carnegie Observatories, La Serena, Chile

¹⁵UK Astronomy Technology Centre, Royal Observatory, Edinburgh EH9 3HJ, UK

<footnote>*Corresponding author - e-mail: abigail.frost@eso.org</footnote>

<footnote>†Deceased</footnote>
Massive stars (those $\ge 8$ solar masses at formation) have radiative envelopes that cannot sustain a dynamo, the mechanism which produces magnetic fields in lower mass stars. Despite this, ~7% of massive stars have observed magnetic fields, the origin of which is debated. We use multi-epoch interferometric and spectroscopic observations to characterise HD 148937, a binary system of two massive stars, finding that only one star is magnetic and that it appears younger than its companion. The system properties and a surrounding bipolar nebula can be reproduced by a model where two stars merged (in a previous triple system) to produce the magnetic massive star. Our results provide observational evidence that magnetic fields form in at least some massive stars through stellar mergers.

Stars with initial masses larger than eight solar masses ($M_{\odot}$) release large amounts of energy into their immediate surroundings and their host galaxies (1). Such massive stars end their lives explosively as supernovae and gamma-ray-bursts, and produce neutron stars and black holes. In close binary systems, pairs of neutron stars or black holes can merge, producing a burst of gravitational waves (2). Massive stars can also experience mergers prior to their explosions, but it is unclear how frequently this occurs or what effect it has on their stellar evolution. Furthermore, over 90% of massive stars exist in binaries and higher-order multiple systems (3), raising the chances of a merger during their lifetime.

If a massive star has a magnetic field, the amount of mass lost through stellar winds is expected to be reduced compared to a non-magnetic massive star. This leaves more stellar mass available at the end of the star's life to form a compact object (4). However, the origin of magnetism in massive stars is not well understood. Lower mass stars like the Sun sustain magnetic fields when convective heating in their envelopes causes the circulation of charged stellar material which then acts like a dynamo. As stars become massive, however, their envelopes
change from convective to radiative, meaning they cannot sustain magnetic fields in this way. Nevertheless, approximately 7% (5, 6) of single O-type stars, a group of stars with birth masses $M \gtrsim 15$ M$_{\odot}$, display large-scale magnetic fields of hundreds to thousands of Gauss (7–9). Of?p stars are a subset of O-type stars that show evidence of magnetism in their optical spectrum and are characterised by their unusually strong C III $\lambda$4650 and N III $\lambda$4634 to $\lambda$4641 emission lines (10).

Several potential origins of magnetic fields in massive stars have been proposed. One possibility is that they could be remnants of the magnetic fields present in material from which the stars formed (11) that were later sustained through convection before the star reached the main-sequence and started burning hydrogen (12). However, it is unclear whether such fields would survive once the stars reach the main sequence (13). Alternatively, magnetic fields could be produced due to the mixing of stellar material during a stellar interaction or merger (14, 15).

## The HD 148937 system

HD 148937 (RA: 16:33:52.387, DEC: -48:06:40.476) is an Of?p star (10). The width of the H$\alpha$ emission line in its optical spectra displays short-period (7.03 d) variability (16–18). This has been considered an indirect indication of a dipolar magnetic field, a stellar wind confined by that magnetic field, and a co-rotating dynamical magnetosphere with an estimated field strength of 1020±300 G (19). Furthermore, HD 148937 is surrounded by a complex bipolar nebula enriched with carbon, nitrogen and oxygen (20).

Interferometric observations have shown that HD 148937 is a binary system containing two stars of almost equal near-infrared brightness (21). Spectroscopic study of the system has suggested two possible orbital periods for the binary: 18.1 yr or 26.2 yr, with corresponding orbital eccentricities of 0.58 and 0.75 respectively (22). This spectroscopic study also implied that the two stars in the binary are a mid-O-type and late-O-type, with only the latter showing signatures
of a magnetically confined stellar wind, in the form of strong Balmer emission (including the Hα emission) and Hε emission (22).

Interferometric observations of HD 148937

We monitored HD 148937 for nine years using the Very Large Telescope Interferometer, or VLTI, at Paranal Observatory in Chile (23). Observations were performed with two different near-infrared instruments; the Precision Integrated-Optics Near-infrared Imaging ExperIMENT (hereafter PIONIER) (24) and the GRAVITY instrument (25). Both instruments determine the interferometric visibilities (a measure of the target's spatial extent and how well resolved it is by the interferometer) and their phases (which indicate its symmetry). Additionally GRAVITY provides spectra in the K-band (from 1.98 to 2.40 µm) at a spectral resolving power R = 4000. We use least-squares minimisation geometrical modelling to measure the separation, orientation and the flux ratio between the primary (defined as the brightest object) and secondary stars for each of the ten observations made throughout the 9 yr observing campaign. We find the K-band brightness of the secondary star to be 94.6±0.41% that of the primary across the averaged GRAVITY observations, detailed in (26). A strong Balmer emission line (Br γ) is present at 2.16 µm in the GRAVITY spectrum, and given the fit to the interferometric data it can only be associated with the primary star, not the secondary. Strong Bry emission is an indirect sign of stellar winds confined by a magnetic field, e.g. (27), and has been used as an infrared indicator of magnetospheres in hot stars (28–30). Because the Bry line arises from only the primary star, we confirm the previous suggestion (22) that only one star in HD 148937 has a strong magnetic field.

From the multi-epoch interferometric data we determine the astrometric positions of both
stars at ten points across their orbit, as can be seen in Figure 1. Fitting an orbital model to
this astrometry combined with the radial velocity data from archival spectra (26) rules out the
<figure><img src="image_2.png" /><figcaption>Figure 1: **The best-fitting orbital model (black solid line) to the astrometric data only, projected onto the plane of the sky.** The relative offsets of the stars are shown as red dots, with the uncertainties of these positions shown as ellipses. The point of periastron passage is shown as a blue square.</figcaption></figure>

shorter of the two previously proposed orbital periods from the archival spectroscopic data
alone (22). Combining this astrometry with the geometric distance of 1155±28 parsecs (pc)
measured by the Gaia satellite (Data Release 3) (31) and with radial velocities from archival
optical spectroscopy (22), we constrain the orbit. We find an orbital period *P* = 25.76±0.82 yrs,
orbital eccentricity *e* = 0.7782±0.0051, orbital inclination *i* = 84.07±0.10° and a total mass of
the system of *M*<sub>total</sub> = *M*<sub>1</sub> + *M*<sub>2</sub> = 56.52 ± 0.75 M<sub>☉</sub> (Table S4).

We apply spectral disentangling to archival spectra to separate them into individual spectra
for both stars. The spectral disentangling technique that we use adopts a grid-based approach
(32–34), fixing orbital parameters constrained by the interferometric observations and adjusting
only the semi-amplitudes (*K*₁ and *K*₂) of each star’s radial velocity curve as free parameters (26).
We find *K*₁ = 28.4<sup>+3.2</sup><sub>-3.6</sub> km s<sup>-1</sup> and *K*₂ = 31.9<sup>-3.4</sup><sub>+3.7</sub> km s<sup>-1</sup> for the primary (more massive) and
secondary star, respectively. Combined with the above constraint on the total mass, these values
correspond to a dynamical mass of the primary star of *M*₁ = 29.9<sup>+3.4</sup><sub>-3.1</sub> M<sub>☉</sub> and a dynamical mass
of the secondary of $M_2 = 26.6_{-3.4}^{+3.0} M_\odot$.

To constrain the physical properties of the stars (26), we compare the disentangled spectra of each component to atmospheric models using a chi-square ($\chi^2$) metric. The atmosphere models are generated with the CMFGEN software (35), as suitable for O-type stars. CMFGEN does not consider spectral features due to the magnetism, and therefore the fit to some spectral lines impacted by the magnetic field is poor in (for example, the NII lines) (26). Such emission features were excluded from the $\chi^2$ fit procedure. We find that the model with the smallest $\chi^2$ has effective temperatures of $T_{\text{eff}} = 37.2_{-0.4}^{+0.9}$ and $35.0_{-0.9}^{+0.5}$ kK, and surface gravities of $\log g = 4.00_{-0.09}^{+0.09}$ and $3.61_{-0.09}^{+0.02}$ for the primary and secondary stars, respectively, where g is in units of cm s$^{-2}$.

The primary star thus appears to be hotter and less evolved than the secondary. The secondary star is enriched in nitrogen and depleted in carbon and oxygen (N/H ratio 8.74 ± 0.10) with respect to a baseline value of $7.78\pm0.10$ (36). While the primary appears to be N-rich, the presence of strong emission lines due to the magnetically confined winds prevents us from quantifying this.

Finally, we also find that the primary is the fastest rotator of the system, with a projected equatorial velocity of $v_{\text{eq}} \sin i = 165 \pm 20 \text{ km s}^{-1}$ compared to $67 \pm 15 \text{ km s}^{-1}$ for the secondary. It has been shown that magnetic fields in massive stars cause momentum loss and slow their rotation (37), a phenomena known as magnetic braking. Thus, this derived fast rotation rate implies that the currently observed magnetic field has not been able to slow down the magnetic star yet. Furthermore, light curves taken with the Transiting Exoplanet Survey Satellite (TESS) suggest a misaligned magnetic axis for the primary star (see Supplementary Text).
An age discrepancy within the binary

The primary is the more massive star in the HD 148937 system and is thus expected to have
evolved the fastest. However, a comparison of the bolometric luminosities and effective tem-
peratures of the two stars, as derived from our atmospheric analysis (26), with evolutionary
tracks (38) reveals that the primary star appears younger than the secondary star. This is dis-
played in Figure 2. We use single star stellar evolution models (38), a Bayesian comparison
method and the Bonn Stellar Astrophysics Interface (BONNSAI) (39) to further quantify this and
also include these results in Figure 2. We consider two cases for the secondary star: one in-
cluding the observed nitrogen enrichment and one without. The BONNSAI results indicate that the
magnetic primary has an estimated age of $2.68_{-0.36}^{+0.28}$ Myr, whilst the secondary has an estimated
age of $4.10_{-0.27}^{+0.29}$ Myr without using the nitrogen enrichment estimate, or $6.58_{-0.82}^{+0.26}$ Myr account-
ing for it. Therefore, the secondary star is older, regardless of whether the nitrogen is considered
or not. This age difference of at least ~1.4 Myr is significant, and allows us to reject the null
hypothesis of the two stars being coeval (i.e., that they have followed the same evolution since
their formation and essentially evolved as single stars) at the 99.5% confidence level.

HD 148937 has no nearby O-star neighbours, and this low local surface density of massive
stars makes a capture scenario involving stars of different ages improbable (40,41). We therefore
infer that both stars in the system formed together at the same time, but the primary star must
have undergone a rejuvenation event, producing the apparent age difference between the two
stars. One possibility is that a mass-transfer event between the two stars has rejuvenated the
primary magnetic star, as has been proposed for Plaskett's star (42,43). In this case, the initially
more massive star of the pair (which in this scenario would have to be the current secondary)
would have grown into a giant or supergiant star before its lower mass companion. Whilst doing
so, it would have exceeded in size the boundary at which its material remains gravitationally
bound (defined as its Roche lobe). Some of the material overflowing its Roche lobe would then
<figure><img src="image_3.png" /><figcaption>Figure 2: Hertzsprung-Russell diagram (HRD) illustrating the difference in luminosity and temperature between the primary ('P') and secondary ('S') stars of HD 148937. Effective temperature ($T_{\text{eff}}$) is shown on the horizontal axis, whilst the bolometric luminosity ($L$) of the stars is shown on the vertical-axis, where $L$ is in units of solar luminosities ($L_{\odot}$). The black, thick line indicates the zero-age main sequence (ZAMS, labelled). Coloured lines are isochrones and evolutionary tracks, with blue lines representing stars with an initial rotation of 165 km s$^{-1}$, whilst red lines represent stars with an initial rotation of 490 km s$^{-1}$. Dashed lines are isochrones for stellar populations with different ages. Solid lines are evolutionary tracks for various initial masses computed at Galactic metallicity (38). Best-fit atmospheric measurements (squares) are given with their $1\sigma$ confidence intervals while the shaded areas give the 2D BONNSAI posterior distributions obtained, in shade of blue, with low initial rotational velocity solutions (which applies to the primary and secondary in the case of no nitrogen enrichment) and, in shade or red, with the high initial rotational velocity solutions (which applies to the secondary accounting for the nitrogen enrichment). The different countour levels corresponds to the 1, 2 and 3$\sigma$-confidence intervals as computed from the highest peak in the posterior distribution.</figcaption></figure>

have been accreted by the companion star (the current magnetic primary). Such a Roche lobe
overflow (RLOF) event causes mass and angular momentum gain, mixing, and the rejuvenation
of the accretor. However, if this process occurred we would expect that the donor star (the
current secondary) would still almost fill its Roche lobe and be visibly much larger than the
primary, which disagrees with the radii that we determine (see Table S4). A former RLOF
event would also reduce the eccentricity of the orbit through tidal forces, producing a close to
circular orbit, which is inconsistent with the eccentricity we measure. Alternatively, today's
primary star could have been rejuvenated in a merger event such that it appears younger than its
companion.

Constraints from the bipolar nebula

The bipolar nebula around HD 148937 could also be formed by a merger. The nitrogen abun-
dance of the nebula (20) is far higher than what can be expected from the surface N enrichment
of the secondary star (26). The most enriched material is in the most distant regions of the neb-
ula from the binary (20) and this level of enrichment is only expected deep in stellar interiors.
This could be explained if a stellar interior was violently disrupted during the production of the
nebula. Removal of the outer hydrogen envelope of a massive star could expose nitrogen-rich
material, which would then be ejected through strong stellar winds, but the star responsible
would then appear to be a Wolf-Rayet star, not a main-sequence O-type star, so we reject this
possibility. Envelope stripping through RLOF can also produce nebulae, but we already rejected
this hypothesis given the sizes we derive for the stars. An alternative mechanism is mass loss
during a merger (44, 45).

In a merger scenario, one expects the lifetime of the nebula to be short. A kinematic age of ~3 kyr was first estimated for the nebula (46) and more recently high-resolution multi-object spectroscopic observations have allowed a minimum age of 7.5 kyr to be determined (47). This is much younger than the ~1.5 Myr it would take for magnetic braking to occur (37). Given the high rotation rate we find for the magnetic star, it is therefore possible that the nebula and
magnetic field were produced by the same event.

Other mechanisms that could form the bipolar nebula include a giant eruption, red supergiant mass-loss or strong winds in a pre-supernova evolutionary stage (48–50). However, each of these other pathways require one or both of the stars to have evolved off the main sequence which is inconsistent with our measured atmospheric parameters. We therefore consider only the merger scenario to be plausible.

<figure><img src="image_4.png" /><figcaption>Figure 3: Merger and rejuvenation models compared to the measured properties of the HD 148937 system. The lower horizontal-axis shows the age of the secondary star ($t_S$) as a reference whilst the age difference between the rejuvenated primary and secondary ($\Delta t$) is shown on the vertical-axis. Each coloured/style of line shows a different model representing a scenario where two stars could have merged to form the current magnetic star in HD 148937. Shaded regions with contours are the observational constraints on the age of the secondary in the case where nitrogen enrichment is considered (grey) and the case where it is not. The upper horizontal-axis shows the orbital period ($P_{\text{orb}}^{\text{merger}}$) the two stars in a 30+5M$_{\odot}$ binary would have to have had at the time of the merger.</figcaption></figure>
## Binary evolution and merger models

We further explore the merger scenario using hydrodynamical models (26). The models calculate the masses of the stars during the merger and the rejuvenation of the product star (51). We find that the merger of two stars can reproduce both the current ~30$M_{\odot}$ mass of the magnetic primary star and its apparent age discrepancy with the current secondary star. We identify several models that match (within $1-\sigma$) the measured age discrepancy and masses of both stars, both with and without considering the nitrogen abundance. This is illustrated in Figure 3. In each case, the total mass of the stars which merged to produce the current magnetic primary star is 35 $M_{\odot}$. Therefore, between 1.7 and 8.2 $M_{\odot}$ would have to have been lost during the merger event, which is consistent with a previously estimated mass-range of the nebula of 1.6 $M_{\odot}$ to 12.6 $M_{\odot}$ (20).

## Summary and implications

We conclude that HD 148937 was originally a higher-order multiple system, likely a triple system, with a close inner binary. This inner binary underwent a merger a few thousand years ago, which produced a magnetic field in the merged star and the nebula surrounding the system.

Our inferred history of the system provides observational support that mergers are a viable source of magnetism in massive stars, as already suggested theoretically (15). The fraction of O stars that are predicted to experience a merger is 8±3% (52), similar to the ~7% fraction which are observed to have magnetic fields (5, 6). Thus, the merger mechanism appears to be the dominant origin of magnetic fields in massive stars.
References

1. N. Langer, *Ann. R. Astron. Astroph.* **50**, 107 (2012).

2. B. P. Abbott, et al., *Physical Review Letters*, **116**, 061102 (2016).

3. H. Sana, et al., *Science* **337**, 444 (2012).

4. V. Petit, et al., *Mon. Not. R. Astron. Soc.* **466**, 1052 (2017).

5. G. A. Wade, et al., *Mon. Not. R. Astron. Soc.* **456**, 2 (2016).

6. J. H. Grunhut, et al., *Mon. Not. R. Astron. Soc.* **465**, 2432 (2017).

7. M. Schöller, et al., *Astron. Astroph.* **599**, A66 (2017).

8. M. Schöller, et al., *Astron. Astroph.* **599**, A66 (2017).

9. J. H. Grunhut, et al., *Mon. Not. R. Astron. Soc.* **465**, 2432 (2017).

10. N. R. Walborn, *Astron. J.* **77**, 312 (1972).

11. J. Braithwaite, H. C. Spruit, *Nature* **431**, 819 (2004).

12. R. Arlt, G. Rüdiger, *Mon. Not. R. Astron. Soc.* **412**, 107 (2011).

13. D. Moss, *Astron. Astroph.* **403**, 693 (2003).

14. L. Ferrario, J. E. Pringle, C. A. Tout, D. T. Wickramasinghe, *Mon. Not. R. Astron. Soc.* **400**, L71 (2009).

15. F. R. N. Schneider, et al., *Nature* **574**, 211 (2019).

16. Y. Nazé, et al., *The Astrophysical Journal*, **135**, 1946 (2008).
17. Y. Nazé, et al., Astron. Astroph. **520**, A59 (2010).

18. G. A. Wade, et al., Mon. Not. R. Astron. Soc. **419**, 2459 (2012).

19. S. Hubrig, et al., Astron. Astroph. **490**, 793 (2008).

20. L. Mahy, et al., Astron. Astroph. **599**, A61 (2017).

21. H. Sana, et al., Astrophys. J. Suppl. **215**, 15 (2014).

22. G. A. Wade, et al., Mon. Not. R. Astron. Soc. **483**, 2581 (2019).

23. X. Haubois, et al., Optical and Infrared Interferometry and Imaging VIII, A. Mérand, S. Sallum, J. Sanchez-Bermudez, eds. (2022), vol. 12183 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 1218306.

24. J. B. Le Bouquin, et al., Astron. Astroph. **535**, A67 (2011).

25. GRAVITY Collaboration, et al., Astron. Astroph. **602**, A94 (2017).

26. Materials and methods are available as supplementary materials.

27. A. Ud-Doula, S. P. Owocki, R. H. D. Townsend, Mon. Not. R. Astron. Soc. **385**, 97 (2008).

28. M. E. Oksala, et al., Astron. Astroph. **578**, A112 (2015).

29. J. P. Wisniewski, et al., Astrophys. J. Letters **811**, L26 (2015).

30. S. D. Chojnowski, et al., Mon. Not. R. Astron. Soc. **516**, 2812 (2022).

31. C. A. L. Bailer-Jones, J. Rybizki, M. Fouesneau, M. Demleitner, R. Andrae, Astron. J. **161**, 147 (2021).

32. M. Fabry, et al., Astron. Astroph. **651**, A119 (2021).
33. T. Shenar, et al., Astron. Astroph. **639**, L6 (2020).

34. L. Mahy, et al., Astron. Astroph. **664**, A159 (2022).

35. D. J. Hillier, D. L. Miller, Astrophys. J. **496**, 407 (1998).

36. N. Grevesse, M. Asplund, A. J. Sauval, P. Scott, Astroph. Sp. Sc. **328**, 179 (2010).

37. A. Ud-Doula, S. P. Owocki, R. H. D. Townsend, Mon. Not. R. Astron. Soc. **392**, 1022 (2009).

38. I. Brott, et al., Astron. Astroph. **530**, A115 (2011).

39. F. R. N. Schneider, et al., Astron. Astroph. **570**, A66 (2014).

40. S. F. Portegies Zwart, S. L. W. McMillan, M. Gieles, Ann. R. Astron. Astroph. **48**, 431 (2010).

41. E. Maraboli, F. Mantegazza, G. Lodato, European Physical Journal Plus **138**, 152 (2023).

42. N. Linder, et al., Astron. Astroph. **489**, 713 (2008).

43. J. H. Grunhut, et al., Mon. Not. R. Astron. Soc. **428**, 1686 (2013).

44. R. Hirai, P. Podsiadlowski, S. P. Owocki, F. R. N. Schneider, N. Smith, Mon. Not. R. Astron. Soc. **503**, 4276 (2021).

45. T. Morris, P. Podsiadlowski, Science **315**, 1103 (2007).

46. C. Leitherer, C. Chavarria-K., Astron. Astroph. **175**, 208 (1987).

47. B. Lim, Y. Nazé, S.-J. Chang, D. Hutsemékers, Astrophys. J. **961**, 72 (2024).

48. Y. H. Chu, R. R. Treffers, K. B. Kwitter, Astrophys. J. Suppl. **53**, 937 (1983).
49. J. Mackey, et al., *Nature* **512**, 282 (2014).

50. V. V. Gvaramadze, et al., *Mon. Not. R. Astron. Soc.* **454**, 219 (2015).

51. E. Glebbeek, E. Gaburov, S. Portegies Zwart, O. R. Pols, *Mon. Not. R. Astron. Soc.* **434**, 3497 (2013).

52. S. E. de Mink, H. Sana, N. Langer, R. G. Izzard, F. R. N. Schneider, *Astrophys. J.* **782**, 7 (2014).

53. A. Mérand, Flexible spectro-interferometric modelling of OIFITS data with PMOIRED (2024).

54. M. Fabry, spinOS: the SPectroscopic and INterferometric Orbital Solution finder (2024).

55. V. Lapeyrere, et al., *Proc. SPIE* (2014), vol. 9146 of *Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series*, p. 91462D.

56. A. Mérand, *Optical and Infrared Interferometry and Imaging VIII*, A. Mérand, S. Sallum, J. Sanchez-Bermudez, eds. (2022), vol. 12183 of *Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series*, p. 121831N.

57. J. B. Le Bouquin, et al., *Astron. Astroph.* **601**, A34 (2017).

58. J. Bodensteiner, et al., *Astron. Astroph.* **641**, A43 (2020).

59. T. Shenar, et al., *Astron. Astroph.* **639**, L6 (2020).

60. J. F. Donati, *Solar Polarization*, J. Trujillo-Bueno, J. Sanchez Almeida, eds. (2003), vol. 307 of *Astronomical Society of the Pacific Conference Series*, p. 41.

61. P. Petit, et al., *Pub. Astron. So. Pac.* **126**, 469 (2014).
62. P. Hadrava, *Astron. Astroph. Supplement*, **114**, 393 (1995).

63. Y. Nazé, N. R. Walborn, F. Martins, *Revista Mexicana de Astronomía y Astrofísica*, **44**, 331 (2008).

64. L. Mahy, et al., *Astron. Astroph. 634*, A118 (2020).

65. K. Pavlovski, H. Hensberge, *Binaries - Key to Comprehension of the Universe*, A. Prša, M. Zejda, eds. (2010), vol. 435 of *Astronomical Society of the Pacific Conference Series*, p. 207.

66. S. V. Marchenko, A. F. J. Moffat, P. R. J. Eenens, *Pub. Astron. So. Pac.* **110**, 1416 (1998).

67. J. F. González, H. Levato, *Astron. Astroph. 448*, 283 (2006).

68. L. Mahy, et al., *Astrophys. J.* **708**, 1537 (2010).

69. L. Mahy, F. Martins, C. Machado, J. F. Donati, J. C. Bouret, *Astron. Astroph. 533*, A9 (2011).

70. R. A. Remillard, J. E. McClintock, *Ann. R. Astron. Astroph. 44*, 49 (2006).

71. N. Zacharias, et al., *American Astronomical Society Meeting Abstracts* (2004), vol. 205 of *American Astronomical Society Meeting Abstracts*, p. 48.15.

72. J. A. Cardelli, G. C. Clayton, J. S. Mathis, *Astrophys. J.* **345**, 245 (1989).

73. F. Martins, B. Plez, *Astron. Astroph. 457*, 637 (2006).

74. T. Morris, P. Podsiadlowski, *Mon. Not. R. Astron. Soc.* **365**, 2 (2006).

75. M. MacLeod, E. C. Ostriker, J. M. Stone, *Astrophys. J.* **868**, 136 (2018).
76. F. R. N. Schneider, P. Podsiadlowski, N. Langer, N. Castro, L. Fossati, *Mon. Not. R. Astron. Soc.* **457**, 2355 (2016).

77. F. R. N. Schneider, et al., *Mon. Not. R. Astron. Soc.* **495**, 2796 (2020).

78. G. R. Ricker, et al., *Journal of Astronomical Telescopes, Instruments, and Systems 1*, 014003 (2015).

79. Lightkurve Collaboration, et al., *Lightkurve: Kepler and TESS time series analysis in Python*, Astrophysics Source Code Library, record ascl:1812.013 (2018).

80. N. R. Lomb, *Astroph. Sp. Sc.* **39**, 447 (1976).

81. J. D. Scargle, *Astrophys. J.* **263**, 835 (1982).

82. S. Burssens, et al., *Astron. Astroph.* **639**, A81 (2020).

83. M. H. Montgomery, D. O'Donoghue, *A derivation of the errors for least squares fitting to time series data*, vol. 13 (1999).

84. A. Schwarzenberg-Czerny, in, *Interplay of Periodic, Cyclic and Stochastic Variability in Selected Areas of the H-R Diagram*, C. Sterken, ed. (2003), vol. 292 of *Astronomical Society of the Pacific Conference Series*, p. 383.

85. D. M. Bowman, et al., *Nature Astronomy* **3**, 760 (2019).

86. D. M. Bowman, et al., *Astron. Astroph.* **616**, A77 (2018).

87. D. W. N. Stibbs, *Mon. Not. R. Astron. Soc.* **110**, 395 (1950).
**Acknowledgements**

Our co-author Rodolfo Barbá sadly passed away before this work was published. We offer our posthumous thanks for his enthusiasm for this work and his contributions to the project.

**Funding:** This project was conceived and executed during the European Union's Horizon 2020 research and innovation programme grant of HS (agreement number 772225: MULTIPLES). AJF was hired using this grant during 2020-2023 as was MF between 2019-2020 and TS between 2020-2022. FRNS also held a European Union's Horizon 2020 research and innovation programme grant (number 945806: TEL-STARS) and was additionally supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster). DMB acknowledges a senior postdoctoral fellowship from Research Foundation Flanders (FWO; grant number 1286521N), the Engineering and Physical Sciences Research Council (EPSRC) of UK Research and Innovation (UKRI) in the form of a Frontier Research grant under the UK government's ERC Horizon Europe funding guarantee (SYMPHONY; grant number [EP/Y031059/1]), and a Royal Society University Research Fellowship (grant number: URF\R1\231631).

**Author contributions:** AJF led the project, analysed the GRAVITY data, and wrote the manuscript. HS co-led the project, prepared the observations of the PIONIER and GRAVITY data, calculated the absolute mass and luminosities and did the evolutionary modelling. AJF and HS did the orbital fitting. LM performed the spectral disentangling, atmospheric analysis and an independent estimate of the luminosities. GW organised the optical spectroscopic data and provided crucial insight into the previous spectroscopic analyses of the system (22). JB and DMB analysed the TESS data. JBLB reduced and analysed the PIONIER data. AM wrote the
software used to analyse the GRAVITY data, reduced the GRAVITY data and assisted with its
analysis. FRNS performed the merger modelling. TS discussed the results and cross-checked
models. MF wrote the code for orbital fitting and assisted with the first attempts at the orbital
fitting. RHB, AF, NIM and JVS observed the raw optical spectroscopic data. PM contributed
to the theoretical interpretation. All co-authors contributed to the discussion and provided feed-
back during the manuscript preparation.

**Competing interests:** We declare no competing interests.

**Data and materials availability:** The GRAVITY data are available through the ESO raw science archive (https://archive.eso.org/eso/eso_archive_main.html) under the Program ID 60.A-9168 (PI: H Sana). The reduced PIONIER data archived at the OIData portal (http://oidb.jmmc.fr/index.html) and can be found by searching for target HD 148937 and using the dates listed in Table S1. They are also accessible through the ESO Archive Science portal http://archive.eso.org/scienceportal/home under program IDs 189.C-0644, 093.C-0503, 596.D-0495, 5100.D-0721, and 105.20FR, (PI: H Sana). The TESS data are available via the MAST data portal https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html (PI: George Ricker) under target name HD 148937. The previously unpublished ESPaDOnS data are available at the PolarBase database of high resolution spectropolarimetric stellar observations at http://polarbase.irap.omp.eu/, and can be retrieved by searching for target HD 148937 and using the dates listed in Table S5.

The code used to analyse the GRAVITY data is available at https://github.com/amerand/PMOIRED and on Zenodo (53). The orbital modelling code is available at https://github.com/matthiasfabry/spinOS and via Zenodo at (54). The BONNSAI tool and models are supplied by the BONNSAI web-service at www.astro.uni-bonn.de/stars/bonnsai.
<figure><img src="image_5.png" /></figure>

Supplementary Materials for:

**A magnetic massive star has experienced a stellar merger**

A. J. Frost, H. Sana, L. Mahy, G. Wade, J. Barron, J.-B. Le Bouquin, A. Mérand,

F. R. N. Schneider, T. Shenar, R. H. Barbá, D. M. Bowman, M. Fabry, A. Farhang, P. Marchant,

N. I. Morrell and J. V. Smoker

Corresponding author: abigail.frost@eso.org

**The pdf file includes:**

* Materials & Methods

* Supplementary Text

* Figures S1 to S15

* Tables S1 to S8

* References 55 to 87
# Methods and Materials

## Interferometric observations

Multi-epoch interferometric observations were performed using two different instruments of the Very Large Telescope Interferometer (VLTI) facility of the European Southern Observatory (ESO) at Cerro Paranal in Chile, which operates at optical and infrared (IR) wavelengths (23). Ten epochs were obtained over a time span of 9 years, from October 2012 to September 2021. Our observations used the four Auxiliary Telescopes (ATs) at the VLTI. Interferometers combine the light from multiple telescopes to observe astronomical sources and as a result probe scientific sources in Fourier space. Depending on the positions of the telescopes in an interferometer with respect to the observed source (baselines), different points in Fourier space (or ‘u-v’ points) are sampled. For AT observations with the VLTI, these points correspond to the different stations the AT telescopes can be placed at, which are named with a combination of letters and numbers (e.g. K0) and can be combined to describe the overall baseline configuration the interferometer is using (e.g. A0-G1-J2-K0). The key observables provided by optical and IR interferometry include the visibilities, closure phases and differential phases. Visibilities correspond to the amplitude of the waves of light received at each telescope and describe the spatial extent of the source (with an unresolved source having visibilities of 1) whilst the phases can tell the observer about the symmetry of the object (e.g. closure phases equal to 0 are associated with a perfectly symmetric source). Each science observation is bracketed by an observation of a calibrator in optical/IR interferometry, to allow the visibilities and closure phases to be calibrated according to the observing conditions at the time of the observation. Table S1 summarises all the interferometric observations, including the stations of the telescopes used at the VLTI.
<table><caption>Table S1: Journal of the interferometric observations. The first column lists the instrument; the second column lists the modified Julian date (MJD) of the observations while the last column gives the configuration of the interferometer. Each letter-number combination (e.g. A0) in the telescope configuration column corresponds to a different station, that is a different location at which one of the VLTI Auxiliary Telescopes can be placed (23).</caption><thead><tr><th>Instrument</th><th>MJD</th><th>Telescope<br>Configuration</th></tr></thead><tbody><tr><td>PIONIER</td><td>56088.066</td><td>A0-K0-G1-I1</td></tr><tr><td>PIONIER</td><td>56868.001</td><td>K0-A1-G1-J3</td></tr><tr><td>GRAVITY</td><td>57557.202</td><td>A0-G1-J2-K0</td></tr><tr><td>GRAVITY</td><td>57559.006</td><td>A0-G1-J2-K0</td></tr><tr><td>PIONIER</td><td>57623.991</td><td>A0-G1-J2-J3</td></tr><tr><td>GRAVITY</td><td>57646.999</td><td>A0-G1-J2-K0</td></tr><tr><td>PIONIER</td><td>57900.117</td><td>B2-K0-D0-J3</td></tr><tr><td>PIONIER</td><td>57995.031</td><td>A0-G1-J2-J3</td></tr><tr><td>PIONIER</td><td>58227.193</td><td>A0-G1-J2-J3</td></tr><tr><td>PIONIER</td><td>59477.016</td><td>A0-G1-J2-J3</td></tr></tbody></table>
K-band interferometry with the GRAVITY instrument

Interferometric data were obtained in June and September 2016 with the GRAVITY instru-
ment (25) (at the Very Large Telescope Interferometer, VLTI) as part of the instrument's sci-
ence verification programme. The data for HD 148937 were taken at spectral resolving power
R=4000 in single-field mode with the Auxiliary Telescopes (ATs). The observables retrieved
by GRAVITY include visibilities, closure phases, and differential phases in addition to the flux
of the source.

Our GRAVITY data were reduced and calibrated using the standard GRAVITY pipeline
(55). The GRAVITY data were analysed with the PMOIRED software (56), which we used to cre-
ate a geometrical model to represent the HD 148937 system from which synthetic observables
were derived to fit to the observed data. Specifically, we used a model of two uniform disks
to fit the data, as illustrated in Figure S1. The position of the primary was fixed at the origin.
The diameters of both uniform disks were fixed to 0.2 milliarcseconds (mas) so they would be

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

Figure S1: **Model images of the primary and secondary star for two epochs in 2016 based on the geometric fits to the GRAVITY data.** a) shows the model image corresponding to the parameters derived from the fit to the June data, whilst b) shows the same for the September data. The primary star (shown by a purple x and labelled ‘1’) in each case is fixed at (0,0) and the position of the secondary (shown by an orange x and labelled ‘2’) is described in relation to this to the East (E) and the North (N).
unresolved at VLTI baselines, as expected for main sequence O-type stars at kiloparsec (kpc) distances. The visibility amplitude, closure phase, differential phase and the normalised flux were all simulated during the fitting process. GRAVITY data were taken on two nights in June 2016 and one night in September 2016. Of the GRAVITY data sets, one data set on each night shows reduced quality across the G1A0 baseline, with reduced wavelength bins visible across the differential phases. We tested including and excluding these data and found negligibly different results, so ultimately all were included.

The total normalised flux of the system shows a strong Bry line (see Figure S2). The best-fitting model includes a Lorentzian line profile in the spectrum of the primary star, but not in

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

Figure S2: **GRAVITY spectra for each star in HD 148937 (taken on 2016-06-18) showing normalised flux against wavelength.** Also shown are the spectra normalised to the median (s1+s2, red solid line) and the ratio of the two spectra (green dotted line). In the spectra of the primary star (blue line) a Bry emission line is visible at ~2.166 µm. Additionally, in the spectra normalised to the median and the ratio of the lines, a He II emission line is visible at ~2.189 µm. These spectra change negligibly in the September data, so these are not displayed.
that of the secondary, to reproduce the strong Br-γ line in the normalised flux. We tested fitting
a Bry line profile for the secondary and including contributions to the emission from both stars.
In these cases it reduced the fit quality for all the data, resulting in a negative flux for both a
Lorentzian or Gaussian line profile. This could possibly indicate a weak absorption line in the
secondary spectrum but, given its flux value was not statistically significant, we removed the
emission line from our secondary star model. The goodness of fit (as determined through the
reduced chi-square, χ²) and the positions and the K-band flux of the secondary star (with respect
to the primary) varied negligibly depending on whether emission line profiles were included or
not.

The best-fitting parameters are listed in Table S2. We find both stars are of similar brightness, with the secondary ~93% to ~96% as bright as the primary in the K-band across the datasets. The change in astrometric position observed between the epochs of GRAVITY data is illustrated in Figure S1. The interferometric data and the best fitting models to them are shown in Figures S3 and S4.

To test the robustness of the derived parameters and of their uncertainties, we also perform
a bootstrapping procedure where we resample the single set to create a variety of simulated
samples. Although bootstrapping usually results in larger uncertainties (as it negates the effects
of correlated data), it can refine the estimation of the companion's parameters and help to show
whether a dataset is consistent. The bootstrap plots associated with the errors on the fits are
shown in Figures S5 and S6 and we use these as our final uncertainties.
<table><caption>Table S2: Parameters derived from the model fits to the GRAVITY observations. The position of the primary was fixed at the origin (0,0). The diameters of both stars were also fixed to 0.2 mas. The first two columns are the date of the observations and the reduced $\chi^2$ goodness-of-fit of the model. The remaining columns show the derived model parameters of the binary system, namely: $f_K$, the flux ratio in the K-band of the secondary to the primary; $\rho$, the angular separation of the two stars; PA, the position angle of the secondary with respect to the primary, measured East (+90°) from North (+0°) in the range 0-360°; $f_{\text{line}}$, the flux ratio of the fitted emission line to the normalised flux; $w_{\text{line}}$, the width of the line; and $\lambda_{\text{line}}$, its wavelength. The uncertainties on $\lambda_{\text{line}}$ are statistical only and all uncertainties are obtained from the bootstrapping procedure.</caption><thead><tr><th>Calendar date</th><th>χ<sup>2</sup><sub>red</sub></th><th>f<sub>K</sub></th><th>ρ<br>[mas]</th><th>PA<br>[°]</th><th>f<sub>line</sub></th><th>w<sub>line</sub><br>[nm]</th><th>λ<sub>line</sub><br>[μm]</th></tr></thead><tbody><tr><td>2016-6-18 to 2016-6-20</td><td>1.19</td><td>0.9311±0.0020</td><td>3.632±0.018</td><td>65.4±0.00268</td><td>0.566±0.0033</td><td>1.039±0.013</td><td>2.1661±0.0000065</td></tr><tr><td>2016-9-15</td><td>1.87</td><td>0.9615±0.0061</td><td>5.449±0.0403</td><td>76.1±0.00025</td><td>0.5073±0.0073</td><td>1.208±0.050</td><td>2.1663±0.000022</td></tr></tbody></table>
<figure><img src="image_8.png" /></figure>

Figure S3: Plot showing the fits of the final geometric model (made with PMORED) to the GRAVITY data taken in June 2016. a) shows the u-v coverage of the observations - the points in Fourier space probed by the telescopes over the different baselines. The baseline between each pair of telescopes possible in the interferometric array are shown in different colours, with the letters and numbers (e.g. K0) corresponding to the different positions of the telescopes (23). The June data were combined in the fit because the nature of the source is very unlikely to differ over a couple of days. This causes the duplicate points in a) and the multiple lines in the remaining subplots. In the remaining subplots the data are shown in black with uncertainties in grey and the best-fitting model in red. The text in each subplot is a different colour corresponding to the different baselines over which the measurement was taken (as in subplot a). Subplot b) shows the normalised flux with wavelength. Subplots c) to f) are the closure phase measurements ('T3PHI') across different baselines, whilst subplots g) to j) are the differential phases ('DPHI') and subplots l) to q) are the visibilities |V|. The cyan squares represent the continuum which is computed using a linear fit.
<figure><img src="image_9.png" /><figcaption>Figure S4: As for Figure S3, but for the data taken in September 2016. These data were taken on one night, hence the lack of duplicated u-v points in subplot a), for example.</figcaption></figure>
**Figure S5: Bootstrap error calculations for the June 2016 GRAVITY/VLTI data.** The final calculated values of each of the parameters (labelled as per Table S2) are at the top of each column with their uncertainties. “2,x” is the x-position of the secondary star and “2,y” is the y-position. These are both in mas and with respect to the origin (0,0) at which the primary star was fixed, as shown in Figure S1. These were used to calculate the separation and position angle of the companion in Table S2 (ρ and PA respectively). The diagonal shows the 1D distribution of values from the simulated bootstrap data compared to the value from the fit. The subplots below the diagonal show the 2D distribution of the simulated bootstrap data for each covariance for each pair of values. ‘c’ (inset) is the correlation factor from the co-variance of the two variables. Grey points are individual fits from the bootstrapping, blue and orange errorbars/ellipses are the $1\sigma$ 1D/2D confidence levels from the bootstrap and from all the data, respectively.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>Figure S6: As for Fig. S5, but for the September 2016 GRAVITY data.</figcaption></figure>
H-band interferometry with the PIONIER instrument

We also observed HD 148937 with the Precision Integrated Optics Near Infrared ExpeRiment (PIONIER) instrument at the VLTI (24). This operates at *H*-band wavelengths and the obser-
vations were also taken with the ATs with a variety of baseline configurations. As with the
GRAVITY data, observations of HD 148937 were alternated with those of calibrator sources
to determine the fringe visibilities and closure phases. The PIONIER data were reduced and
calibrated using the *PNDRS* package (24). Seven concatenations were taken with Modified Julian
Dates listed in Table S3. The earliest of these observations was analysed in a previous publica-
tion (21). PIONIER, not having the spectrointerferometric capabilities of GRAVITY, provides
the observer with visibilities and closure phases only.

The PIONIER data were analysed through the same geometrical modelling methods and
code of (57). As for the GRAVITY data, the interferometric data were fit using a binary model
composed of uniform disks 0.2 mas in diameter to represent two stars that are unresolved at
VLTI baselines. The free parameters were the angular separation ρ of the two stars and the flux
ratio *f*<sub>H</sub> between the secondary and the primary star (defined to be the brightest object in the
*H*-band). The flux ratio was fixed for the first and fourth observations due to low data quality
(starred in Table S3) and fit for the others. The position angles (PA) derived from PIONIER
measurements suffer from a ±180° degeneracy because the two components are almost equal
flux in the *H*-band. We used a comparison with the GRAVITY data to lift this degeneracy. The
resulting astrometry from the PIONIER data is listed in Table S3. The smallest separations
correspond to the epochs of largest radial velocity (RV) variation in a previous spectroscopic
study (22) and are compatible with a time of periastron passage (*T*<sub>0</sub>) around 2013 to 2014.
<table><caption>Table S3: Parameters derived from the model fits to the PIONIER observations. The first two columns are the date of the observations and the reduced $\chi^2$ of the model. The remaining columns show the derived model parameters of the binary system: the flux ratio in the H-band ($f_H$) of the secondary to the primary; the separation $\rho$ of the two stars; the position angle (PA) of the secondary with respect to the primary and; the parameters of the astrometric uncertainty ellipse, its semi-major ($e_{max}$) and semi-minor ($e_{min}$) axes of the astrometric uncertainty ellipse, as well as the position angle ($e_{PA}$) of its semi-major axis.</caption><thead><tr><th>Calendar date</th><th>χ<sup>2</sup></th><th>f<sub>H</sub></th><th>ρ<br>[mas]</th><th>PA<br>[°]</th><th>e<sub>max</sub><br>[mas]</th><th>e<sub>min</sub><br>[mas]</th><th>e<sub>PA</sub><br>[°]</th></tr></thead><tbody><tr><td>2012-6-10</td><td>0.79</td><td>0.96*</td><td>21.06</td><td>100.42</td><td>0.87</td><td>0.35</td><td>130</td></tr><tr><td>2014-7-30</td><td>0.41</td><td>0.96</td><td>3.60</td><td>266.54</td><td>0.18</td><td>0.13</td><td>139</td></tr><tr><td>2016-8-23</td><td>0.35</td><td>0.96</td><td>5.00</td><td>73.52</td><td>0.17</td><td>0.12</td><td>154</td></tr><tr><td>2017-5-27</td><td>0.46</td><td>0.96*</td><td>10.67</td><td>84.25</td><td>0.95</td><td>0.19</td><td>116</td></tr><tr><td>2017-8-30</td><td>0.33</td><td>0.96</td><td>12.46</td><td>85.89</td><td>0.24</td><td>0.18</td><td>120</td></tr><tr><td>2018-4-19</td><td>1.47</td><td>0.95</td><td>16.67</td><td>88.05</td><td>0.23</td><td>0.18</td><td>56</td></tr><tr><td>2021-9-20</td><td>2.11</td><td>0.96</td><td>34.05</td><td>92.78</td><td>0.50</td><td>0.50</td><td>0</td></tr></tbody></table>

<footnote>*fixed during fitting</footnote>
Determining the orbit

Two possible orbital solutions for HD 148937 have been previously determined (22): a longer-
period, larger eccentricity solution with *e* = 0.75 and *P* ~ 26 yr, and a less eccentric, shorter-
period one with *e* = 0.58 and *P* ~ 18 yr. Previous work favoured the longer period solution
(*P* ~ 26 yr) because it was more consistent with the variability observed in the HeI λ5876 line
in that work.

We use the additional astrometry from the GRAVITY and PIONIER observations to further constrain the orbital solution. We use the SPectroscopic and INterferometric Orbital Solution software (SPINOS) to constrain the 3D orbit (32), combining previous radial velocity (RV) data (22) with the new astrometry from VLTI. We tested leaving the distance as a free parameter and found that it converged on a value of d = 1135±5 pc, consistent with the Gaia DR3 value (31). We therefore fixed our distance to the Gaia distance to reduce the number of free parameters.

With the distance fixed, the free parameters are the orbital period (P), eccentricity (e), inclination (i), and time of periastron passage (T₀). In addition, we adjust a common systemic velocity γ for both RV curves, the argument of the periastron passage of the secondary star (ω₂) measured with respect to the position of the primary star in the relative orbit, the argument of the ascending node (Ω) as well as the total mass Mtotal = M₁ + M₂ of the system.

As a first step, we constrained the orbit using only the interferometric data, i.e. without using
previous radial velocity (RV) measurements (22). This fit excludes the 18-yr period allowed by
previous study, as it lies outside the range of potential fit values which can be achieved using
the astrometric data. As a second step we included both the astrometric and RV measurements,
adopting an equal weight for both data-types during the model fitting. In each step, a first
solution was determined through a Levenberg-Marquardt minimisation of the χ² which was
then refined using Markov Chain Monte Carlo (MCMC) which was also calculated with SPINOS.
Figure S7 shows the fitted RV curves from this final orbital fit, while the MCMC plot associated
<figure><img src="image_12.png" /></figure>

Figure S7: Orbital model fitted to the RV data (22) and the interferometric data corresponding to the solution listed in Table S4. As is customary for the presentation of spectroscopic binaries, we display more than one period for the system which means that the data points between phases 0.9-1.2 are the same as those between phases -0.9-0.2. Blue and red circles represent data for the primary and secondary star, respectively, whilst the dotted lines correspond to the model for each star following the same colour-code. The absolute RV amplitude of the secondary star in this figure may have been impacted by the respective companion star and is not to be blindly trusted (see text).

with the minimisation is shown in Figure S8.

We find a high eccentricity ($e \sim 0.8$) and a near edge-on orientation of the orbital plane with respect to the line of sight ($i \sim 85^\circ$). The best-fitting orbital parameters are provided in Table S4. The root-mean-square (rms) residuals of the final fit are $rms_{RV1}=2.9 \text{ km s}^{-1}$ and $rms_{RV2}=4.2 \text{ km s}^{-1}$ for the primary and secondary RV curves, respectively. The rms of the relative astrometric orbit is $rms_{AS}=0.06 \text{ mas}$.

The previous study used as the source of the RVs (22) did not fully disentangle the spectral contribution of both components and focused on lines that are dominated by one or the other companion, respectively. However, even a small cross-contamination of a diagnostic spectral line of one star by a weak line of the companion star may significantly bias the measured RVs (32, 58, 59). We expect that such contamination would mostly impact the derived semi-
amplitudes, $K_1$ and $K_2$, and, to a lesser extent, the measured eccentricity. As we show later, the diagnostics lines used for RVs in (22) seem to show such small contamination. We thus refrain to give $K_1$ and $K_2$ values at this stage of the analysis and we refer to later sections of our Materials & Methods which describe how we constrain more reliable $K_1$ and $K_2$ values.

<table><caption>Table S4: Best-fitting parameters of the orbital modelling, found through fitting the interferometric data jointly with RV data of (22). The meaning of the symbols is described in the text. $T_0$ is expressed in MJD.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3" style="text-align: center;">Orbital fit</td></tr><tr><td>P</td><td>days</td><td>9390±300</td></tr><tr><td>e</td><td></td><td>0.7782±0.0051</td></tr><tr><td>i</td><td>°</td><td>84.07±0.10</td></tr><tr><td>T<sub>0</sub></td><td>day</td><td>56958.2±2.8</td></tr><tr><td>γ</td><td>km s<sup>-1</sup></td><td>-24.15±0.26</td></tr><tr><td>ω<sub>2</sub></td><td>°</td><td>340.10±0.41</td></tr><tr><td>Ω</td><td>°</td><td>277.27±0.26</td></tr><tr><td>M<sub>total</sub></td><td>M<sub>☉</sub></td><td>56.52±0.75</td></tr><tr><td colspan="3" style="text-align: center;">Goodness of the fit</td></tr><tr><td>Degrees of freedom</td><td></td><td>94</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td></td><td>0.74</td></tr><tr><td>rms<sub>RV1</sub></td><td>km s<sup>-1</sup></td><td>2.9</td></tr><tr><td>rms<sub>RV2</sub></td><td>km s<sup>-1</sup></td><td>4.2</td></tr><tr><td>rms<sub>AS</sub></td><td>mas</td><td>0.056</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure S8: Density plots of the MCMC chains projected on a 1D axis (top panels) and 2D planes for each of the adjusted parameters of the orbital model. The blue lines indicate the model with the lowest $\chi^2$. The contours give the 39.3, 67.5, 86.4% confidence. Outside the outer counters, realisations are plotted individually as cloud of points.</figcaption></figure>
Additional ESPaDOnS spectra of the system

In addition to the published optical spectroscopic data (22), we also use additional archival spec-
tra from the ESPaDOnS instrument (60), accessed via the PolarBase database of high resolution
spectropolarimetric stellar observations (61). ESPaDOnS is the Echelle SpectroPolarimetric
Device for the Observation of Stars at the Canada France Hawaii Telescope (CFHT). It is a
bench-mounted, high-resolution echelle spectrograph and spectropolarimeter with an operating
wavelength regime of 370 to 1,050 nm and resolving power ranging from 68,000 to 81,000. The
supplementary observations from 2009 and 2010 covered a wavelength range of 370 to 900 nm
and employed the ‘star + sky’ instrumental mode with R=68,000, whilst the 2018 data were
taken at a resolution of R=65,000. The data reduction for these observations was performed
using UPENA, the original data reduction package for ESPaDOnS. A full list of the data used are
presented in Table S5.
<table><caption>Table S5: Journal of the additional archival ESPaDOnS observations used in this work to complement those from previous study (22). We list the dates of the observations in Heliocentric Julian Date (HJD) for consistency with this previous work.</caption><thead><tr><td>Calendar date</td><td>HJD -2400000</td></tr></thead><tbody><tr><td>2009-05-08</td><td>54959.946</td></tr><tr><td>2010-06-19</td><td>55366.855</td></tr><tr><td>2010-06-20</td><td>55367.862</td></tr><tr><td>2010-06-21</td><td>55368.834</td></tr><tr><td>2010-06-22</td><td>55369.861</td></tr><tr><td>2010-06-23</td><td>55370.821</td></tr><tr><td>2010-06-24</td><td>55371.818</td></tr><tr><td>2010-06-25</td><td>55372.824</td></tr><tr><td>2010-07-24</td><td>55401.754</td></tr><tr><td>2010-07-25</td><td>55402.750</td></tr><tr><td>2010-07-27</td><td>55404.763</td></tr><tr><td>2010-07-28</td><td>55405.739</td></tr><tr><td>2010-07-29</td><td>55406.729</td></tr><tr><td>2018-06-21</td><td>58290.879</td></tr><tr><td>2018-06-22</td><td>58291.862</td></tr><tr><td>2018-06-24</td><td>58293.879</td></tr><tr><td>2018-06-23</td><td>58292.865</td></tr><tr><td>2018-06-25</td><td>58294.858</td></tr><tr><td>2018-06-26</td><td>58295.884</td></tr><tr><td>2018-06-27</td><td>58296.871</td></tr><tr><td>2018-06-28</td><td>58297.842</td></tr><tr><td>2018-06-29</td><td>58298.828</td></tr><tr><td>2018-06-30</td><td>58299.868</td></tr><tr><td>2018-07-01</td><td>58300.853</td></tr><tr><td>2018-07-04</td><td>58303.839</td></tr><tr><td>2018-07-05</td><td>58304.859</td></tr></tbody></table>
## Spectral disentangling

In order to better understand the stars in HD 148937, one requires a better understanding of their individual atmospheric properties. In order to retrieve these, we separate the spectral signature of the two stars using spectral disentangling in order to avoid cross-contamination that may bias the atmospheric parameters.

The spectral disentangling approach that we use separates the spectral signatures of each star without relying on previously measured RVs. To reduce the number of degrees of freedom, we fixed most of the orbital parameters to the values listed in Table S4 but allowed the semi-amplitudes of the RV curves ($K_1$ and $K_2$) to vary over a small grid, between 0-60 km s⁻¹ (see Figure S9). This grid-disentangling approach has been tested using artificial datasets and applied to other long-period binary star systems in previous work (32, 58, 59).

The spectroscopic data contains spectra sampled over different epochs. Because HD 148937 exhibits spectral variability with a period of 7.03 days (17, 18), we built a master spectrum at each epoch to remove such non-orbital variabilities. These master spectra were then used for the spectral disentangling process.

To separate the spectroscopic features of both components, we apply the grid-based approach using a Fourier disentangling code (62) on spectral lines including He I+II λ4026 and He I λ4471. The resulting χ²-map is shown in Figure S9 and indicates that the values of $K_1$ and $K_2$ are not correlated. The best-fitting values are $K_1 = 28.4_{-3.6}^{+3.2}$ km s⁻¹ and $K_2 = 25.4_{-14.9}^{+15.5}$ km s⁻¹. The resulting disentangled spectra are displayed in Figures S10 and S11. The spectrum of the primary (magnetic) component in HD 148937 bears similarity to the disentangled spectra of another Of?p star, HD 108, as determined by (63). Both stars display the same strong N III, C III and He I emission (present between 4600 – 4700 Å) and a Hβ P-Cygni profile at ~ 4860 Å.

Fourier spectral disentangling has the disadvantage of losing the continuum when the sys-
<figure><img src="image_14.png" /><figcaption>Figure S9: **Reduced $\chi^2$ map from the grid method of disentangling.** The minimal value at $K_1 = 28.4 \text{ km s}^{-1}$ and $K_2 = 25.4 \text{ km s}^{-1}$, is denoted with a red dot. The solid white contour is the $1\sigma$ level. The background gradient corresponds to the value of the $\chi^2$ across the grid, with the colour-bar as reference. The side panels are 1-D cut-through views of the axes of the $\chi^2$ map. The orange lines delimit the 68% confidence interval on the sum of the semi-amplitudes of both RV curves ($K_1 + K_2$) as discussed in the 'Dynamical masses' subsection.</figcaption></figure>

The light curve does not present total eclipses. This introduces distortions in the disentangled continuum, which can affect the measurement of some stellar parameters such as the surface gravity (20, 64, 65). In the case of HD 148937, the shapes of the lines of the magnetic primary could prevent the detection of the continuum through the spectral lines. We therefore also tried another separation technique, the shift-and-add method (66, 67). This technique has
the disadvantage of deforming the wings of broad lines (68, 69), yet one can mitigate this by using atmospheric models as a template. Because the secondary spectrum is clear of emission features, we modelled the output spectrum produced from the Fourier disentangling technique using the grid of models and performing a $\chi^2$ analysis on hydrogen and helium lines (mainly sensitive to the surface gravity and effective temperature, respectively). We used the CMFGEN (35) best-fit atmosphere model (described in the 'Atmospheric analysis' section) as input template for the shift-and-add technique. The $K_1$ value of 28.4 km s$^{-1}$ from the Fourier-based disentangling was consistent with previous work (22), so we keep it fixed and only apply the shift-and-add disentangling over the $K_2$-grid in steps of 5 km s$^{-1}$. The shift-and-add method yields $K_2 = 30.5_{-11.6}^{+13.5}$ km s$^{-1}$. The 5 km s$^{-1}$ difference between the $K_2$ values from the two methods is not significant given the uncertainties, thus the results from the two methods are in agreement.

## Dynamical masses

As shown in Figure S9, $K_1$ is better constrained than $K_2$ by our grid spectral disentangling. However better constraining $K_2$ is of interest as the stellar masses of the individual stars in a binary can be related to the semi-amplitudes using Kepler's third law and the binary mass function (70). Additional information on the mass of the system can therefore provide a further constraint on the semi-amplitudes of the individual stars and allow a better constraint to be made on $K_2$. In this context, we make use of the total mass measurement of the system, determined independently from the values of $K_1$ and $K_2$ via the astrometric orbital solution ($M_{total} = M_1 + M_2 = 56.52 \pm 0.75$ M$_\odot$ as seen in Table S4). Using the definition of $K_1$ and $K_2$ and Monte Carlo simulations to propagate the uncertainties on P, e, i and $M_{total}$, we convert the constraint on the total mass derived from our orbital fitting into constraints on the sum of the semi-amplitudes of the individual RV curves. We find $K_1 + K_2 = 61.34 \pm 0.93$ km s$^{-1}$ and show this constraint on Figure S9. With these much tighter constraints on both $K_1$ and $K_2$, we adopt $K_1 = 28.4_{-3.6}^{+3.2}$ and
$K_2 = 31.9_{-3.7}^{+3.4}$ km s⁻¹ as our final values, where the reversed uncertainty notation indicates that upper limits on $K_1$ correspond to lower limits on $K_2$. Similarly, we obtain $M_1 = 29.9_{-3.1}^{+3.4}$ and $M_2 = 26.6_{-3.4}^{+3.0} M_\odot$. This is a much tighter constraint than was obtained from the disentangling alone, but remains in agreement with both the grid disentangling and the shift-and-add results.

Our final values are summarised in Table S6, alongside with the derived linear dimensions of the system: the semi-major axis of the relative orbit ($a = a_1 + a_2$), the semi-major axes of the primary and secondary barycentric orbits ($a_1$ and $a_2$, respectively); the radii ($R$) of the stars ($R$) relative to the size of their Roche lobe ($R_{RL}$).

<table><caption>Table S6: Summary of the dynamical and geometrical parameters HD 148937 from the combined steps of our analysis.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Spectral disentangling</td></tr><tr><td>K<sub>1</sub></td><td>km s<sup>-1</sup></td><td>28.4<sup>+3.2</sup><sub>-3.6</sub></td></tr><tr><td>K<sub>2</sub></td><td>km s<sup>-1</sup></td><td>31.9<sup>-3.4</sup><sub>+3.7</sub></td></tr><tr><td colspan="3">Dynamical masses</td></tr><tr><td>M<sub>1</sub></td><td>M<sub>☉</sub></td><td>29.9<sup>+3.4</sup><sub>-3.1</sub></td></tr><tr><td>M<sub>2</sub></td><td>M<sub>☉</sub></td><td>26.6<sup>+3.0</sup><sub>-3.4</sub></td></tr><tr><td colspan="3">Linear dimensions</td></tr><tr><td>a</td><td>au</td><td>33.45 ± 0.73</td></tr><tr><td>a<sub>1</sub></td><td>au</td><td>15.8 ± 1.4</td></tr><tr><td>a<sub>2</sub></td><td>au</td><td>17.7 ± 1.4</td></tr><tr><td>R<sub>1</sub>/R<sub>RL</sub></td><td></td><td>(3.46 ± 0.29) × 10<sup>-3</sup></td></tr><tr><td>R<sub>2</sub>/R<sub>RL</sub></td><td></td><td>(4.64 ± 0.33) × 10<sup>-3</sup></td></tr></tbody></table>
Atmospheric analysis

We derive the atmospheric parameters of both components of HD 148937 by modelling the
disentangled spectra using the CMFGEN stellar atmosphere code (35). For this purpose, we built a
grid of synthetic atmosphere models covering a range of effective temperature (T<sub>eff</sub>) of 27 kK ≤
T<sub>eff</sub> ≤ 45 kK with steps of 1 kK and a range of surface gravities from 3.0 ≤ log(g/cm s<sup>-2</sup>) ≤ 4.3
with steps of 0.1 dex, where g is expressed in cm s<sup>-2</sup>. Each model was subsequently convolved
with:

i. a rotational profile corresponding to the a series projected rotational velocity (v<sub>eq</sub> sin i)
 ranging from 0 to 200 km s<sup>-1</sup> in steps of 10 km s<sup>-1</sup>;

ii. a radial-tangential profile corresponding to macroturbulence velocities (v<sub>macro</sub>) ranging
from 0 to 200 km s<sup>-1</sup> in steps of 10 km s<sup>-1</sup> and;

iii. an instrumental broadening representing by a Gaussian kernel with a full-width-at-half
maximum corresponding to the spectra resolving power of the observations.

We focus our comparison of the model to the data on hydrogen and helium lines, as these
are sensitive to the surface gravity (g) and effective temperature (T<sub>eff</sub>) respectively. The best-
fitting CMFGEN models of each star are displayed in Figure S10 and S11. Poor fits are visible for
the primary for some emission lines, notably the Balmer lines and the N<sup>II</sup> lines. This is due to
the fact that the CMFGEN models do not include processes such as magnetic winds, which are the
origin of these features. This is reflected in the errorbars of the different stellar parameters and
the emission features were excluded from the χ<sup>2</sup> fit procedure. We compute χ<sup>2</sup> for each model
of the grid and the global χ<sup>2</sup> distributions for the primary and the secondary components are
given in Figures S12 and S13.
<figure><img src="image_15.png" /></figure>

Figure S10: **CMFGEN** best-fit model (red) of the disentangled spectrum of the magnetic primary star (black). Each of the subplots displays a region of the spectrum focused on a different spectral line.

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

Figure S11: As for Figure S10, but for the non-magnetic, secondary star.
<figure><img src="image_17.png" /><figcaption>Figure S12: χ² distribution for the parameters used in the atmospheric analysis of the primary star which result in the final fit shown in Figure S10. The diagonal black and white plots show the χ² distribution as a function of each parameter and the colour plots below show 2D χ² maps. The red line marks the 1σ confidence level and the orange line shows the interpolated value of the reduced χ² across the steps of the grid. The colour-bar indicates the value of the reduced χ².</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure S13: **As for Figure S12, but for the secondary star.** The resulting model is shown in Figure S11.</figcaption></figure>
Uncertainties are computed using a $\chi^2$ increase threshold with respect to the best $\chi^2$ value, computed to encapsulate the 68% confidence interval on each parameters. This selects a family of models with an acceptable quality of fit. The coarseness of the grid that we use however requires us to interpolate between the grid data points to obtain more precise uncertainties (illustrated in Figures S12 and S13). Whenever possible, all statistical uncertainties are given with two significant digits, to avoid significant round-off errors. However, we note that the true precision is limited by the physics in the atmospheric models. Because of this, we adopt minimum errors of 0.2 kK on $T_{\text{eff}}$ and 0.05 dex on log g

For the secondary star we use further additional models to derive the nitrogen, carbon and oxygen contents at the surface of the star. In order to do this we fix the projected rotational velocity ($v_{\text{eq}} \sin i$) of the stars and their macroturbulent velocities ($v_{\text{macro}}$) to the values determined through the $\chi^2$ fitting. From this we determine that the secondary component is enriched in nitrogen and depleted in carbon and oxygen on its surface. While the primary appears to be nitrogen-rich, the spectrum is contaminated by emission lines from the magnetically-confined wind so we cannot estimate its nitrogen abundance. We provide the stellar parameters derived for each star in Table S7.

<table><caption>Table S7: Atmospheric and physical parameters of the two stars in HD 148937, with their $1\sigma$ confidence intervals. No value is derived for the nitrogen enrichment ($\epsilon_N$) for the primary star due to contamination from the lines associated with magnetism in this star. We denote this with the ‘...’ symbol.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Primary</th><th>Secondary</th></tr></thead><tbody><tr><td>$T_{\text{eff}}$</td><td>kK</td><td>37.2<sup>+0.9</sup><sub>-0.4</sub></td><td>35.0<sup>+0.2</sup><sub>-0.9</sub></td></tr><tr><td>log g</td><td>cm s<sup>-2</sup></td><td>4.00<sup>+0.09</sup><sub>-0.09</sub></td><td>3.61<sup>+0.05</sup><sub>-0.09</sub></td></tr><tr><td>$v_{\text{eq}} \sin i$</td><td>km s<sup>-1</sup></td><td>165±20</td><td>67±15</td></tr><tr><td>$v_{\text{macro}}$</td><td>km s<sup>-1</sup></td><td>160±38</td><td>78±12</td></tr><tr><td>$f_i/f_{\text{tot}}$</td><td>V-band</td><td>0.55±0.02</td><td>0.45±0.02</td></tr><tr><td>$\epsilon_N$</td><td>[log+12]</td><td>...</td><td>8.74 ± 0.10</td></tr><tr><td>log L/L<sub>☉</sub></td><td></td><td>5.28±0.06</td><td>5.19±0.07</td></tr></tbody></table>
Given the uncertainty on $K_1$ and $K_2$, and the possible impact on the extracted spectra, we also extract the stellar parameters using the $K_1$ and $K_2$ values up to $3\sigma$ away from the best-fitting values. Differences in the resulting atmospheric parameters remain small and within errorbars.

To compute the bolometric luminosities ($L$) of the two stars, we first estimate interstellar extinction $A_ν$ by matching the spectral energy distribution (SED) of HD 148937 with atmospheric models. The magnitudes of the system are taken from the NOMAD catalogue (71), and the fluxes come from the CMFGEN models computed for each component and assuming a flux ratio of primary to secondary in the V-band (500-580 nm) of $1.13\pm0.05$. We obtain $A_ν = 1.89\pm0.02$, accounting for the IR excess from the nebula. The SED is shown in Figure S14. Using the average Galactic extinction curve for the diffuse interstellar medium of $R_ν = 3.1$ (72), the interstellar reddening is computed to be $E(B-V) = 0.61\pm0.02$. Using $A_ν = 1.89$, the Gaia DR3 geometric distance of $1155\pm28$ pc (31) and the average GRAVITY K-band brightness ratio (Table S2), we then compute the absolute K-band magnitude of each object. Using the

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

Figure S14: **Spectral energy distribution used to estimate the stellar reddening towards HD 148937.** The model is shown in black and the data in red. The IR excess is due to the nebula.
relation $A_K/A_V = 0.123$, and bolometric corrections (73), we compute luminosities of the two components of $\log(L/L_\odot) = 5.28 \pm 0.06$ for the primary and $\log(L/L_\odot) = 5.19 \pm 0.06$ for the secondary.

To summarise our atmospheric analysis, we find that the primary is intrinsically more luminous, hotter and rotates faster than the secondary. The secondary also has a lower surface gravity than the primary, suggesting that it is more evolved despite its lower mass.
Stellar evolution modelling

We compare the effective temperatures, luminosities, abundances and masses we derive to
single star stellar evolution models from the BONNSAI web service (38, 39), which performs a
Bayesian comparison with a pre-computed model grid at solar metallicity. We consider two
sets of inputs for the secondary, one including the nitrogen measurements derived from the at-
mospheric analysis and one without them. BONNSAI identifies single-star evolutionary models
that reproduce the observed properties of the primary and secondary star in both cases (Ta-
ble S8). The models indicate that the primary is more massive, hotter and more luminous than
the secondary. The posterior probability distributions suggest a high or low initial rotational ve-
locity, depending on whether the nitrogen abundance measurement of the secondary is included
or not, respectively. High rotational mixing, hence a high initial spin, is required to reproduce
the nitrogen enrichment observed for the secondary. The magnetic primary has a model age of
2.68<sup>+0.28</sup><sub>-0.36</sub> Myr, whilst the secondary is 4.10<sup>+0.29</sup><sub>-0.27</sub> Myr when nitrogen enrichment is not considered
and 6.58<sup>+0.26</sup><sub>-0.82</sub> Myr when it is. Independent of the nitrogen constraints, an age difference is ev-
ident between the two components. We thus reject the hypothesis that the stars are coeval at
the 99.5% level of significance. Either the stars are not the same age, or they did not evolve as
single isolated stars.

Our analysis assumes that the secondary non-magnetic star evolved as a single star and so
can be used as a reference clock. We cannot, however, exclude the possibility that the secondary
had a nearby companion, such that HD 148937 was initially a quadruple formed by two close
binaries or a higher-order system. With such initial conditions, it is conceivable that the sec-
ondary itself interacted with a companion and so also underwent rejuvenation. In this scenario,
the comparison of the secondary to single-star evolutionary tracks would underestimate its true
age. This would further increase the age discrepancy. Quantitatively, our models would then
favour a more equal mass merger (Figure 3) and might offer an alternative avenue to explain
the nitrogen abundance and low (current) rotation rate of the secondary. We note however that
there is no evidence for this more complex scenario.

As a final consistency check, we increase the uncertainties of the atmospheric parameters
to 1.0 and 0.5 kK in T<sub>eff</sub> for the primary and secondary, respectively, and to 0.1 dex in log g, in
order to make sure that the age discrepancy that we found is not the result of underestimated
uncertainties. While we obtain a slightly lower age for the secondary, the age discrepancy that
we find between the primary and secondary remains significant and we still reject the hypothesis
that the stars are coeval at the 99.5% confidence level.
<table><caption>Table S8: Evolutionary parameters obtained from the comparison with solar metallicity single-star evolutionary models. The listed values of replicated observables and predicted stellar parameters give the mode of the posterior probability distributions and the highest posterior density intervals (68%). The symbol ‘...’ denotes that no value was determined. ‘N’ is nitrogen.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Primary</th><th colspan="2">Secondary</th></tr><tr><td></td><td></td><td></td><td>without N</td><td>with N</td></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Input parameters</td></tr><tr><td>M<sub>current</sub></td><td>M<sub>☉</sub></td><td>29.9<sup>+3.4</sup><sub>-3.1</sub></td><td>26.6<sup>+3.0</sup><sub>-3.4</sub></td><td></td></tr><tr><td>T<sub>eff</sub></td><td>kK</td><td>37.2<sup>+0.9</sup><sub>-0.4</sub></td><td>35.0<sup>+0.2</sup><sub>-0.9</sub></td><td></td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td></td><td>5.28 ± 0.06</td><td>5.19 ± 0.06</td><td></td></tr><tr><td>log g</td><td>[cm s<sup>-2</sup>]</td><td>4.00<sup>+0.09</sup><sub>-0.09</sub></td><td>3.61<sup>+0.05</sup><sub>-0.09</sub></td><td></td></tr><tr><td>v<sub>eq</sub> sin i</td><td>km s<sup>-1</sup></td><td>165 ± 20</td><td>67 ± 15</td><td></td></tr><tr><td>ε<sub>N</sub></td><td>[log+12]</td><td>...</td><td>...</td><td>8.74 ± 0.10</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Replicated observables</td></tr><tr><td>M<sub>current</sub></td><td>M<sub>☉</sub></td><td>30.0<sup>+1.3</sup><sub>-1.4</sub></td><td>27.4<sup>+1.3</sup><sub>-1.5</sub></td><td>24.2<sup>+1.4</sup><sub>-0.9</sub></td></tr><tr><td>T<sub>eff</sub></td><td>kK</td><td>37.75<sup>+0.60</sup><sub>-0.73</sub></td><td>33.46<sup>+0.75</sup><sub>-0.60</sub></td><td>33.82<sup>+0.83</sup><sub>-0.56</sub></td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td></td><td>5.25<sup>+0.05</sup><sub>-0.05</sub></td><td>5.24<sup>+0.06</sup><sub>-0.04</sub></td><td>5.23<sup>+0.06</sup><sub>-0.04</sub></td></tr><tr><td>log g</td><td>[cm s<sup>-2</sup>]</td><td>3.93<sup>+0.05</sup><sub>-0.05</sub></td><td>3.68<sup>+0.04</sup><sub>-0.04</sub></td><td>3.68<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>v<sub>eq</sub> sin i</td><td>km s<sup>-1</sup></td><td>170<sup>+15</sup><sub>-26</sub></td><td>60<sup>+21</sup><sub>-11</sub></td><td>60<sup>+9</sup><sub>-8</sub></td></tr><tr><td>ε<sub>N</sub></td><td>[log+12]</td><td>...</td><td>...</td><td>8.77<sup>+0.03</sup><sub>-0.06</sub></td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Model stellar parameters</td></tr><tr><td>M<sub>initial</sub></td><td>M<sub>☉</sub></td><td>31.4<sup>+1.5</sup><sub>-1.5</sub></td><td>28.4<sup>+1.5</sup><sub>-1.6</sub></td><td>26.2<sup>+1.5</sup><sub>-1.3</sub></td></tr><tr><td>v<sub>eq,initial</sub></td><td>km s<sup>-1</sup></td><td>190<sup>+76</sup><sub>-44</sub></td><td>80<sup>+83</sup><sub>-32</sub></td><td>510<sup>+4</sup><sub>-48</sub></td></tr><tr><td>v<sub>eq,actual</sub></td><td>km s<sup>-1</sup></td><td>180<sup>+72</sup><sub>-42</sub></td><td>70<sup>+60</sup><sub>-28</sub></td><td>360<sup>+21</sup><sub>-14</sub></td></tr><tr><td>Age</td><td>Myr</td><td>2.68<sup>+0.28</sup><sub>-0.36</sub></td><td>4.10<sup>+0.29</sup><sub>-0.27</sub></td><td>6.58<sup>+0.26</sup><sub>-0.82</sub></td></tr><tr><td>R</td><td>R<sub>☉</sub></td><td>9.69<sup>+0.73</sup><sub>-0.58</sub></td><td>12.32<sup>+0.76</sup><sub>-0.64</sub></td><td>12.21<sup>+0.55</sup><sub>-0.85</sub></td></tr><tr><td>X<sub>He</sub></td><td></td><td>&lt;0.27</td><td>&lt;0.27</td><td>0.37<sup>+0.06</sup><sub>-0.01</sub></td></tr><tr><td>ε<sub>N</sub></td><td>[log+12]</td><td>7.79<sup>+0.23</sup><sub>-0.12</sub></td><td>7.66<sup>+0.06</sup><sub>-0.04</sub></td><td>8.77<sup>+0.03</sup><sub>-0.06</sub></td></tr></tbody></table>
Binary merger scenarios for the primary star in HD 148937

Our derived age of the primary and secondary stars in HD 148937 are inconsistent, with the primary, magnetic star appearing to be younger than the secondary star by at least ~1.4 Myr. Because the atmospheric properties derived for the secondary indicate it is not filling its Roche lobe, recent binary mass-transfer phase cannot explain this rejuvenation. We consider a merger as the most likely cause.

For a merger to have occurred, HD 148937 must have originally been a triple (or higher-order multiple) system in which two or more components merged to form the now magnetic primary star. The current secondary star then serves as a reference clock and we regard its age (6.58$_{-0.82}^{+0.26}$ Myr when its nitrogen enrichment is considered or 4.10$_{-0.27}^{+0.29}$ Myr when it is not) as the age of the entire HD 148937 system. This assumes that the entire star system formed together, and we consider capture of an additional object as unlikely due to the low stellar density of HD 148937's neighbourhood.

During a stellar merger, mass is ejected and can form a bipolar nebula (44, 74, 75) as is observed in HD 148937. The nebula has been investigated using high-resolution spectroscopy taken with the GIRAFFE spectrograph at the Very Large Telescope (VLT) (47). Using a combination of integrated intensity maps, position-velocity diagrams of the two bright lobes, geometric modelling and Monte-Carlo radiative transfer techniques, the morphology and kinematics of the nebula have been determined (47). It was found that the outermost south-eastern lobe is blueshifted and the outermost north-western lobe is redshifted. Using the Gaia DR3 distance a lower-limit of the kinematic age of nebula is determined to be 7.5 kyr. This implies that the nebula is very young, so we assume that the current mass of the primary star is equal to its post-merger mass.

Rejuvenation by a stellar merger can explain the apparently younger age of the primary star and put constraints on the progenitor system (76). Rejuvenation in a merger occurs through
the mixing of fresh hydrogen fuel from the outer envelope into the cores of the involved stars (15, 77). Following a merger, the apparent age discrepancy $\Delta t = t_{\text{true}} - t_{\text{app}}$ between the true age of the merged star, $t_{\text{true}}$, and the apparent age inferred from single star models, $t_{\text{app}}$, is a function of the masses of the merging binary stars and the age of the binary at merger (76).

We adopt predictions of rejuvenation in binary-star mergers from models (76). These models assume that the stellar mass decreases with time due to stellar winds (which are a crucial consideration for massive stars). The models also assume that a fraction $\phi$ of the total mass of an interacting binary is lost during the merger, that the composition of the lost material is the same as the composition of the system, and that the mass of the merger product is equal to a fraction $1 - \phi$ of the original mass of the two stars. The average hydrogen mass fraction decreases linearly with time from its ZAMS value to its final value pre-merger. We can thus determine the fractional main sequence age of the merger product by comparing the hydrogen mass fraction of the merger product to the fractional main sequence age of a single non-merged star of the same mass. Rejuvenation in the core of the merger product is included in the calculation of the apparent fractional main sequence age. The modelling approach also uses results from smoothed-particle hydrodynamics (SPH) simulations of head-on collisions of high-mass mergers in open star clusters (51). The structure of the SPH model merger products is fed into a 1D stellar evolutionary code. An extra mixing parameter $\alpha$ is included and is set to 1.14, as expected for high-mass mergers (51, 76).

These models show that the merger of a close binary system can reproduce both the current 30 M$_{\odot}$ mass of the magnetic primary star and the observed age discrepancy between today's components. They can reproduce the observed stellar properties within $1\sigma$, both with and without including the nitrogen abundance. Four example merger cases are shown in Figure 3. Without including the nitrogen constraint, the merger of a 30M$_{\odot}$ star and a 5M$_{\odot}$ star can explain the observed age discrepancy within $1\sigma$. In the case where the nitrogen enrichment is consid-
ered, the merger of either a 20$M_{\odot}$ star and a 15$M_{\odot}$ star or a 25$M_{\odot}$ star and a 10$M_{\odot}$ star can explain the difference within $1\sigma$. The predicted ejecta masses in all these merger models are similar to the inferred mass of the bipolar nebula of HD 148937 (20).
# Supplementary Text

## Rotation

The presence of a strong magnetic field is expected to increase the combined moment of inertia of the star and its stellar wind, allowing a magnetically confined wind to efficiently remove angular momentum from the star (37). This phenomenon, known as magnetic braking, is predicted to rapidly spin down stars with magnetically confined winds on time scales of a few Myr (37). Using the physical parameters that we derived in Table S7 and the formulae in (37) we estimate the fractional critical rotation rate ($f_c = 0.255 \pm 0.086$), where critical rotation rate refers to rate at the equator of a rotating body beyond which the centrifugal force will exceed Newtonian gravity. We also calculate the wind confinement ($\eta_\star \approx 10$) and spin-down timescale ($\tau_{\text{spin}} \approx 1.5$ Myr), with uncertainties of ~50%. The fractional critical rotation rate that we derive is a factor of two larger than previously reported (18) for the magnetic star ($f_c = 0.12^{+0.11}_{-0.05}$). While the multiple nature of HD 148937 was not known in previous work, the two estimates are consistent within their uncertainties. The values of $\eta_\star$ and $\tau_{\text{spin}}$ are also consistent, because the mass-to-radius ratio is similar in our study and in previous work (18).

These estimates indicate that the magnetic field in HD 148937 should spin down the primary star on a time scale of $\tau_{\text{spin}} \sim 1.5$ Myr. This implies that the merger occurred very recently. A recent merger could explain another property of the system - the bipolar nebula with an estimated life-time of ~7500 yrs (47). This timescale is a few percent of the spin-down time, so the magnetic field has not had time to slow the rotation of the primary star. Therefore, we conclude that the current rotation rate is likely to be similar to the rotation rate immediately after the merger.

The theoretically expected rotation rate of a merger product is uncertain. Simple angular momentum conservation indicates that a merger product should be critically rotating imme-
<figure><img src="image_20.png" /></figure>

Figure S15: Mean magnitude ($\Delta T_p$) subtracted Sector 12 and 39 TESS light curves of HD 148937. BJD is barycentric Julian date. a) Extracted sector 12 TESS light curves (black points) normalised to have a mean of zero. The red line denotes a single frequency model fitted to the combined light curve. b) Same as a) but for sector 39. c) Lomb-Scargle periodogram of the full light curve. The single identified frequency at $v = 0.14359 \pm 0.00002 \text{ d}^{-1}$ (corresponding to $P_{\text{TESS}} = 6.9645 \pm 0.0008 \text{ d}^{-1}$) is marked by the red triangle and line. d) Full light curve phased to the identified period. The red line denotes the single frequency model fit and the white circles denote ten binned data points that are evenly spaced in phase.

diately after merging, meaning that the star is rotating so fast its centrifugal force is close to exceeding its own Newtonian gravity. However, hydrodynamical simulations (15, 77) show that large amounts of angular momentum can be extracted from the merger remnant by a disk that forms during the merger process, such that the merger remnant can spin down efficiently. For example, during the thermal relaxation the restructuring of the stellar interior is predicted to increase the moment of inertia by a factor of about 20, rapidly spinning down the merger product to about 12% critical within a few thousand years (15, 77). The fraction of the critical rotation for HD 148937's magnetic component is slightly larger than – but within $1.4\sigma$ of – these predictions from hydrodynamic simulations. The theoretical values were derived for a
9M☉+8M☉ binary so might not apply to a merger with a large initial mass ratio, as we suggest for HD 148937 (Figure 3).

Previous work has determined that the Hα line profile in HD 148937 varies with a spectroscopically determined period P=7.03 d, which was interpreted as rotational modulations of the magnetic field signal (16, 18). To determine the rotation periods of both stars in the system, we consider archival data from the Transiting Exoplanet Survey Satellite (TESS) (78). TESS has observed HD 148937 in two sectors (survey periods), sector 12 and sector 39. We examine both the 2-min cadence processed light curve and 30-min full frame image extracted light curves from TESS. A periodogram analysis shows dominant periods of 1.66 d in the 2-minute data and 6.97 d in the 30-min cadence data (Figure S15). The light curves were extracted using an eight-by-eight-pixel mask centred on the system. As a single TESS pixel corresponds approximately to a ~21" × 21" region of the sky, the extracted light curves contain contributions from both O star components. We correct for systematics using the LIGHTKURVE software's regression routine (79) using six single-scale co-trending basis vectors provided by the TESS pipeline.

We perform our frequency analysis on the combined light curve using a Lomb-Scargle (80, 81) pre-whitening procedure (82). To avoid over-interpretation of stochastic signals we adopt a conservative significance criterion requiring the signal-to-noise (S/N) to be greater than 5, calculated using a window size of 1 d⁻¹ in the final residual periodogram. Uncertainties on the parameters are estimated using a well-established correlation correction factor (83, 84). We identify a single noteworthy frequency ν = 0.14359±0.00002 d⁻¹ with amplitude A = 1.20±0.03 millimagnitude (mmag) in the combined light curve (with a signal to noise ratio of 6.9). This frequency corresponds to a period P_TESS = 6.9645±0.0008 d which is close to, but not consistent with, the spectroscopic period P_spec = 7.032 ± 0.003 d (16, 18). A second peak can be seen at 0.091 d⁻¹ with a magnitude of ~0.9 mmag, but is not significant. Additionally, its amplitude decreases to ~0.6 mmag after the first frequency is extracted. A single period does
not capture the entire set of variabilities of the TESS light curve (Figure S15), since massive
stars commonly exhibit stochastic low-frequency variability in time-series photometry (85).

The 7.03 d spectroscopic period could be reconciled with the observed rotation of the mag-
netic component if the star has a misaligned magnetic axis. A harmonic of the rotation period is
sometimes the dominant peak in the frequency spectra of rotationally-modulated oblique mag-
netic stars (86). If the axis is misaligned, we would observe the emission from the dipoles of
the magnetic star twice per rotation cycle (86, 87), which is within the uncertainty of the esti-
mated rotation rate from Table S8. If the system had gone through a chaotic exchange before
the merger, similar to (44), this could have produced or increased any misalignment.