© The Authors 2025

1 Introduction

Binary stars constitute a significant fraction of the stellar population in the Galaxy. Observations indicate that over half of massive stars are born in binary or multiple systems (Duchêne & Kraus 2013; Sana et al. 2012), with possible interactions between components shaping their subsequent evolution. One of the key processes in such systems is mass transfer, which occurs when one star fills its Roche lobe and begins transferring material to its companion. This mass exchange profoundly alters the evolution of both stars, often leading to the rejuvenation, expansion, and spin-up of one component or stripping of the other (Podsiadlowski et al. 1992; de Mink et al. 2013).

Mass transfer in close binaries leads to a range of structural and orbital consequences. As mass is accreted by the initially less massive star, the binary orbit responds according to angular momentum conservation. The increasing mass of the accretor and corresponding mass ratio reversal typically result in orbital expansion during the later stages of the interaction (Soberman et al. 1997; Renzo et al. 2023). However, the transfer of mass from the more massive donor initially causes the orbit to shrink (e.g. Renzo et al. 2019). At the same time, the internal structure of the accreting star undergoes significant changes (Renzo & Götberg 2021). As material is deposited onto the surface, it is gradually incorporated into the stellar interior through thermohaline, convective, and rotational mixing processes, particularly near the core-envelope boundary (Neo et al. 1977). This mixing can supply the convective core with fresh hydrogen, effectively rejuvenating the star. The phenomenon of core rejuvenation leads to an increase in convective core mass and alters the subsequent nuclear burning history of the accretor. Moreover, the associated mixing processes can also alter the convective core size, internal chemical gradients, and rotational profile (Moe & Di Stefano 2017; Renzo & Götberg 2021). These changes in internal stratification have direct consequences for the star's pulsational properties (Guo 2021; Miszuda et al. 2022; Wagg et al. 2024; Henneco et al. 2024, 2025).

Asteroseismology offers a unique window into stellar interiors by analysing oscillation modes sensitive to different regions within the star (e.g. Aerts et al. 2010). In particular, gravity modes (g modes) probe stable, stratified regions and, specifically, the chemical gradient and convective core boundaries; whereas the pressure modes (p modes) offer constraints on the sound-speed profile in outer layers (Miglio et al. 2008). In massaccreting stars, modifications to the internal structure due to accretion are expected to leave characteristic signatures in both g (Wagg et al. 2024; Henneco et al. 2024, 2025) and p mode frequencies (Miszuda 2025b). In particular, the rejuvenation from mass transfer shifts the p mode to higher frequencies (Guo et al. 2017) and the g mode period spacing to larger values (Guo & Li 2019). Arras et al. (2006) found that mass accretion can shift the edge of g mode instability strip of the white dwarf star in cataclysmic variables.

Recent observational and theoretical works have identified such asteroseismic signatures in post-mass transfer binaries. Studies of slowly pulsating B-type stars revealed period spacing deviations that are consistent with altered chemical profiles from mass accretion (Wagg et al. 2024). Similarly, models of premain sequence stars indicate that accretion histories shape their pulsational behaviour and internal structure (Zwintz & Steindl 2022).

In this work, we investigate the impact of mass accretion on the internal mixing, structure, and pulsational properties of massive, 10 M_☉ accretor stars using detailed binary evolution modelling. By comparing accretor models with single-star counterparts, we identified key structural changes and their observational consequences, advancing the understanding of the role of mass transfer in stellar evolution.

2 Methodology

2.1 Evolutionary calculations

We used the Modules for Experiments in Stellar Astrophysics code (mesa, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023, version 23.05.1) to construct a toy model under the non-rotating approximation. The model was created both as a product of binary interactions (through mass gain) and as a single star, allowing us to compare the binary model to its isolated counterpart. The details on the variety of physical processes, which MESA relies on, can be found in Appendix A.

In our evolutionary computations, we adopted the AGSS09 (Asplund et al. 2009) initial chemical composition and used the OPAL opacity tables, supplemented by data from Ferguson et al. (2005) for lower temperatures. For high-temperature regimes, as well as hydrogen-poor or metal-rich conditions, we used CO enhanced tables. We assumed a metallicity of Z=0.014 and an initial hydrogen abundance of X₀=0.723.

Convective instability in the models was treated using the Ledoux criterion, combined with the turbulent convection theory based on the Kuhfuss (1986) model, employing a mixing-length parameter of α_MLT=0.5. In regions that were stable according to the Ledoux criterion, but unstable by the Schwarzschild criterion, we applied the semiconvective mixing following the formalism of Langer et al. (1985) with a scaling factor α_sc=0.1. To address regions exhibiting an inversion in the mean molecular weight, such as those formed during mass accretion, we applied the thermohaline mixing using the formalism of Kippenhahn et al. (1980), with an α_th=1 coefficient. As our models did not include rotational mixing, we introduced a minimum diffusive mixing coefficient of D=100 cm² s⁻¹ to smooth out numerical noise or discontinuities in internal profiles, such as the BruntVäisälä frequency or composition gradients. Additionally, we accounted for overshooting beyond the formal convective boundaries, using an overshooting parameter, f_ov, to capture turbulent motions extending into the radiative zone. We applied an exponential scheme of Herwig (2000) on the top of the H-burning core with the value of f_ov=0.02.

To ensure numerical convergence, we explored the numerical stability of our results against changes in both spatial and time step resolution, testing various configurations and ultimately adopting mesh_delta_coeff = 0.2 and time_delta_coeff = 0.5. We have also tested the impact of the adopted value of the minimum diffusive mixing coefficient, ultimately adopting D=100 cm² s⁻¹. We selected this value after testing various options, specifically D=10,20,50, 100, and 200, identifying it as the lowest value that effectively smoothed out numerical noise in the aforementioned profiles without erasing any significant physical structures.

Our initial binary model consists of a 10 M_☉ donor in a 3-day circular orbit with a 7 M_☉ accretor. The initial orbital period was selected so that the donor fills its Roche lobe and initiates mass transfer while still on the main sequence. For simplicity, we assumed fully conservative mass transfer (MT), following the description of Kolb & Ritter (1990). In this scenario, the initially less massive star accretes 3 M_☉, leading to a full mass ratio reversal. Once the accretor reached 10 M_☉, all binary interactions were stopped, and we assumed that the subsequent evolution proceeded as a single star, following the standard path until central hydrogen depletion. To simulate this behaviour in MESA, we used the routine DETACH_BINARY by Mathieu Renzo¹.

For the purposes of comparison, we also calculated single-star evolution for stars with masses from 7 M_☉ to 10 M_☉ with a 0.5 M_☉ step with the physics input as described above.

2.2 Pulsational calculations

To calculate pulsations using the previously described MESA models, we employed the GYRE code (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020; Sun et al. 2023, version 7.1), in the adiabatic approximation. For radial pulsations (ℓ=0), we scanned the frequency range from 1 to 150 d⁻¹ using a linear grid. For higher mode numbers (non-radial pulsations, ℓ=1,2,3,4), we applied an inverse grid sampling between 0.1 and 3.0 d⁻¹ and a linear grid from 3.0 to 150 d⁻¹. The GYRE calculations were performed for all models in our MESA grid. These computations were applied to both single-star models and accretor models from binary evolution alike.

3 Results

3.1 Binary evolution

The early evolution of the system, initially containing 10 M_☉ and 7 M_☉ components, follows the single-star evolutionary paths, namely, where both stars detach from the zero age main sequence (ZAMS) and start synthesizing hydrogen in their cores. Their subsequent evolution is illustrated in the Hertzsprung-Russell (HR) diagram, in Figure 1. As hydrogen burning proceeds in their convective cores, both components start to increase their radius and move through the HR diagram at different paces. By the time of 10 M_☉ star synthesized almost all central hydrogen into helium (approximately 10% hydrogen left in the core), it has grown enough to fill its Roche lobe (at a radius of approximately 9.3 R_☉) and starts to transfer mass onto its companion (I_don in Figure 1). At this point the system is approximately 23.55 Myr old and the components have begun to diverge from their single track equivalents.

As the donor star begins to transfer mass through the inner Lagrange point L₁ onto its companion during a case A mass transfer episode (Kippenhahn & Weigert 1967), it remains in thermal equilibrium. The mass transfer proceeds stably, yet rapidly, with the rate reaching up to ·M ∼10^−3.5 M_☉ yr⁻¹. The MT occurs on the nuclear timescale, and is primarily driven by the gradual expansion of the donor's radius resulting from internal evolutionary changes. As mass is transferred onto the initially less massive companion (starting from point I_acc in Figure 1), the accretor adjusts to the increasing mass by gradually expanding in size and evolving off its original evolutionary track. During MT, the accretor remains close to thermal equilibrium, while accommodating the incoming material. As a result, it experiences a steady increase in mass, luminosity, effective temperature, and radius. We allowed the accretor to grow in terms of mass until it reached 10 M_☉, namely, until a full mass ratio reversal (II acc in Figure 1). After this point, we artificially detached the system, which allowed the accretor to regain thermal equilibrium in the Kelvin-Helmholtz time scale. We then followed the subsequent evolution of the accretor as a single 10 M_☉ star (starting from point III _acc in Figure 1) until the end of the main sequence phase (I V_acc in Figure 1).

Fig. 1 The HR diagram with the evolutionary tracks of the binary system components. The donor star is shown with the purple line and the accretor with the gray line, both plotted alongside single-star reference models for different masses (dashed lines). The mass transfer/accretion rate, log |·M|, is overplotted along the tracks using colour coding, as indicated by the colour bar. The key evolutionary changes are marked: the donor fills its Roche lobe and begins mass transfer (Idon), the accretor starts to gain mass and departs from its original track (Iacc), full mass ratio reversal when the accretor reaches 10 M☉ and the following detachment of the system (IIacc and IIdon), thermal readjustment of the accretor (IIIacc), and the accretor's single-star evolution until the end of the main sequence (I Vacc).

3.2 General properties of the accretor model

In mass-accreting models, the incoming mass has a severe effect on the internal structure. As the transferred material binds to the star, the induced higher central temperature boosts the nuclear energy generation rate. As the core luminosity goes up, the radiative gradient, ∇_rad, increases, leading to a larger convective core region according to the convective instability criterion, as discussed in more detail later in this section.

The evolution of the convective core during this process is shown in Figure 2, where the core mass, M_cc, is plotted against the central hydrogen abundance, X_c. For the single-star model, the convective core steadily recedes from ZAMS to terminal age main sequence (TAMS). In contrast, the initially smaller convective core of the accretor follows a similar trend only until X_c ≍ 0.5 (I_acc, see Figure 1), after which it expands due to mass accretion. Within just 0.02 Myr, M_cc increases from 1.62 M_☉ to 3.3 M_☉(II_acc), reflecting the sudden influx of fresh fuel. Once mass transfer ceases, the star regains thermal equilibrium and the core mass stabilizes at M_cc=2.87 M_☉ with X_c ≍ 0.6 (III _acc). From that point onward, the evolution of the accretor's convective core closely tracks that of the single-star counterpart, although it remains slightly more massive by about 1%. This offset likely results from enhanced mixing across the core-envelope boundary (CEB, Neo et al. 1977; Renzo et al. 2023), which altered the chemical profile and increased the core mass.

The structural changes brought about by accretion are clearly visible in the evolution of internal mixing and hydrogen abundance profiles, shown in Figure 3, based on the convection, overshoot, thermohaline, and user-defined artificial mixing processes. Mass accretion is relatively short-lived and rapid and occurs between the left, second and third top panels. This figure allows us to clearly trace the evolution of the convective core boundary in mass coordinates, both during the standard core recession and the core rejuvenation phase.

The effect of rejuvenation is twofold. As the freshly accreted material settles onto the stellar surface, an inversion in the mean molecular weight profile forms and gives rise to thermohaline mixing, which contributes to the internal redistribution of chemically stratified layers (Kippenhahn et al. 1980; Stancliffe & Glebbeek 2008; Wagg et al. 2024). The expansion of the convective core during accretion leads to the engulfment of layers with a hydrogen gradient, left after the core recession, mixing hydrogen-rich material from the envelope into the core (Neo et al. 1977). At the same time, this process compresses the surrounding chemically stratified region, steepening the hydrogen gradient in the adjacent layers. As a result, the mean molecular weight, μ, drops sharply in this narrow zone, and an inflection in the density profile emerges. Specifically, a slightly increased density in the location of the previously engulfed regions is followed by a steeper decline due to drop in μ, compared to the single models as shown in Figure 4. We refer to this feature as a local bump in the density profile. This bump reflects a deviation from the otherwise monotonic density gradient and emerges just above the growing convective core, in layers enriched in hydrogen as a signature of structural reconfiguration due to accretion. As the core expands during the mass accretion, it shifts this bump into higher layers, near the end of the hydrogen gradient region. This feature is stable, meaning that it persists at least until the end of the main sequence phase (see Figure 4). The presence of this bump is also consistent with earlier findings by Renzo et al. (2023) and Wagg et al. (2024), though the models of the just mentioned authors describe a local ‘dip’ in the density profile.

A key consequence of the density bump is the development of off-centre convective zones (oCZs). These zones appear in the chemically stratified region above the convective core during the late and post-RLOF phases, as seen in the right panels of Figure 3. Their occurrence has been reported by Renzo & Götberg (2021) and their physical origin was analysed by Miszuda (2025a). The formation of oCZs is governed by the Ledoux criterion for convective instability. According to this condition, a region becomes unstable to convection if the actual temperature gradient exceeds the Ledoux gradient, $$\nabla >{\nabla }_{\mathrm{L}}={\nabla }_{\mathrm{ad}}+\frac{\phi }{\delta }{\nabla }_{\mu },$$(1) where ∇ is the temperature gradient, ∇_ad is the adiabatic temperature gradient, and ∇_μ is the gradient of the mean molecular weight, $$\nabla =\frac{\mathrm{d}lnT}{\mathrm{d}lnP},{\nabla }_{\mathrm{ad}}={\left(\frac{\partial lnT}{\partial lnP}\right)}_{\mathrm{ad}}\text{and}{\nabla }_{\mu }=\frac{\mathrm{d}ln\mu }{\mathrm{d}lnP}.$$(2) The thermodynamic derivatives, $$\delta =-{\left(\frac{\partial ln\rho }{\partial lnT}\right)}_{P,\mu }\text{and}\phi ={\left(\frac{\partial ln\rho }{\partial ln\mu }\right)}_{P,T},$$(3) are obtained from the equation of state (e.g. Kippenhahn & Weigert 1990). In radiative layers, where energy is transported by radiation, the actual temperature gradient, ∇, is equal to the radiative temperature gradient, ∇_rad. Therefore, in these layers, the Ledoux instability criterion takes the form of $${\nabla }_{\text{rad}}>{\nabla }_{\mathrm{L}}.$$(4) In the post-mass transfer models, the appearance of a local density bump leads to an increase in opacity, which enhances ∇_rad. If this exceeds the stabilising contribution from ∇_μ, the stability Ledoux criterion is violated and a convective zone forms away from the centre, giving rise to the oCZ.

The oCZs play a critical role in reshaping the hydrogen abundance profile. As they develop within chemically stratified zones, they homogenise the local composition by erasing sharp gradients. The mixing within the oCZ reduces the mean molecular weight gradient leading to nearly uniform hydrogen abundance profiles followed by new discontinuities at the oCZ boundaries. These discontinuities increase ∇_L and can suppress convection locally, effectively splitting or terminating the oCZs. At the same time, the chemically stratified region continues to be compressed by the outward growth of the convective core during the rejuvenation. This compression can again steepen ∇_rad, allowing it to exceed ∇_L and temporarily revive convection. During the mass transfer episode, multiple density bumps form at the boundaries of successive oCZs (Miszuda 2025a) altering the hydrogen abundance profiles. However, only the final feature, formed just before the disappearance of the last oCZ at the end of the mass-accretion phase, persists into the post-mass transfer phases. Because it forms sufficiently far from the convective core, it is not smoothed out by overshooting and thus remains as a lasting structural relic of the accretion history. The impact of this feature on the chemical structure is illustrated in Figure 5, where we compare hydrogen abundance profiles at various evolutionary phases between the single-star (dashed lines) and accretor (solid lines) models. This abundance kink remains clearly visible for X_c ≤ 0.5, for 3.5<m/M_☉<3.8, highlighting its persistence throughout the post-mass transfer main sequence evolution.

The structural imprint of the density bump is also reflected in the sound-speed profile. Since the adiabatic sound speed is defined as $$c_s^2={\Gamma }_{1}P/\rho ,$$(5) where Γ₁ is the first adiabatic exponent, P is the pressure, and ρ is the density; hence, any local modulation in the density at fixed pressure alters the sound speed. In the region near the CEB, where pressure gradients are smooth due to hydrostatic equilibrium, the density bump causes a clear deviation in the c_s profile, as shown in Figure 6. This modulation has observable consequences, as it affects the propagation of p modes and may alter mode trapping behaviour (e.g. Dziembowski et al. 1993; Miglio et al. 2008). We investigate these effects in more detail in the following sections.

Fig. 2 Evolution of mass of the convective core, Mcc, in a function of central hydrogen abundance, Xc, for the single-star (orange lines) and the accretor (blue lines) models. We mark key evolutionary changes of the accretor according to Figure 1.

Fig. 3 Mixing profiles, D, during the evolution of the accretor. Left: evolution from Xc=0.7 to 0.0 in steps of 0.1. Right: profiles at selected evolutionary stages, before, during and after the RLOF. The shaded areas indicate the regions over which the mass varies, with the current mass indicated by dashed lines. The black lines, as denoted on the right-hand ordinates of the right panels show the corresponding profiles of hydrogen abundances.

Fig. 4 Density, log ρ, profiles at selected central hydrogen abundances, Xc, for the single-star (orange lines) and the accretor (blue lines) models. The local bump discussed in the text is clearly visible in the binary models.

Fig. 5 Hydrogen abundances, X, in a function of a mass for different evolutionary stages of 10 M☉ single (orange lines) and accretor (blue lines) models. The vertical lines mark features corresponding to Figure 7. For clarity, we only show a fraction of the inner mass coordinate.

Fig. 6 Sound speed, cs, profiles at selected central hydrogen abundances, Xc, for the single-star (orange lines) and the accretor (blue lines) models.

3.3 Asteroseismic properties of the accretor model

There are two characteristic frequencies that determine the behaviour of stellar oscillations. These are the Lamb frequency (Lamb 1917), $$S_{\ell }^2=\frac{\ell (\ell +1)c_s^2}{r^2},$$(6) and the Brunt-Väisälä frequency (Väisälä 1925; Brunt 1927), $$N^2=g(\frac{1}{{\Gamma }_{1}P}\frac{\mathrm{d}P}{\mathrm{d}r}-\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}r}).$$(7) Here, ℓ is the spherical harmonic degree, c_s is the adiabatic speed of sound, g is the local gravitational acceleration, and r is the radial coordinate.

The propagation of modes is determined by the relation between the oscillation frequency, ω, and the characteristic frequencies, N and S_ℓ². Modes for which ω²<N² and ω²<S_ℓ² propagate as gravity modes (g modes), for which buoyancy acts as the restoring force. They are confined to radiative zones, where N²>0, while convective zones act as evanescent barriers (e.g. Cox 1980; Aerts et al. 2010). In chemically homogeneous, non-rotating and non-magnetic stars, the g modes with consecutive radial orders, n, and the same spherical degree, ℓ, exhibit equally spaced periods, Δ P. According to the asymptotic properties of g modes, valid for ℓ ≪ n, this spacing remains constant and is given by (Tassoul 1980): $$\Delta P_{\ell ,\text{asym}}=\frac{2{\pi }^2}{\sqrt{\ell (\ell +1)}}{\left({\int }_{C}N{r}^{-1}\mathrm{d}r\right)}^{-1},$$(8) where the integration is performed over the g mode propagation cavity, C, with respect to the radial coordinate. This expression shows that the period spacing is determined by the integral of the Brunt-Väisälä frequency and is therefore highly sensitive to the presence of chemical gradients and sharp structural features in the deep interior of the star. In case of the main sequence massive stars, with convective core and radiative envelope, the observed value of Δ P is directly related to the extent of the core (e.g. Moravveji et al. 2015, 2016; Mombarg et al. 2019; Pedersen et al. 2018, 2021).

In contrast, modes with ω²>N² and ω²S_ℓ² propagate as pressure modes (p modes), for which pressure acts as the restoring force. These propagate in the outer layers, where the sound speed and density is lower and the acoustic cavity is located. In massive main sequence stars such as β Cephei variables, only low- to intermediate-order p modes are observed. These modes of moderate radial order can penetrate into the deep stellar interior, making them valuable probes of the region near the convective core. Because of this dependence, p mode frequencies can provide a complementary asteroseismic diagnostic to g mode periods, especially in stars that have undergone structural changes due to processes such as mass accretion. Changes in the internal c_s profile, for instance, due to convective core expansion, envelope compression, or compositional readjustment, can subtly shift the p mode frequencies and thus leave observable signatures in the mode spectrum.

3.3.1 The Brunt-Väisälä frequency

The study of the Brunt-Väisälä frequency N profile is essential to understand how mass accretion modifies the internal structure of massive stars. In particular, it allows us to trace the chemical and thermal stratification that govern the propagation of buoyancy modes in a star. Following the ideal gas plus radiation pressure for a fully ionised gas approximation, N can be expressed (e.g. Unno et al. 1989) as: $$N^2\approx \frac{g^2\rho }{P}[\frac{4\beta -3}{\beta }({\nabla }_{\mathrm{ad}}-\nabla )+{\nabla }_{\mu }],$$(9) where β=P_gas/P is the ratio of gas pressure to total pressure.

The chemical gradient, ∇_μ, is extremely important in the asteroseismology of massive stars with convective cores. As the star evolves on the main sequence, the fully mixed core shrinks leaving a gradient of hydrogen (see Figure 5) and helium abundance. This causes a positive contribution to ∇_μ and, hence, an emergent rise in N (Kippenhahn & Weigert 1990). As the star evolves and the convective core recedes, the peak extends to the deeper parts of the star in a transition region adjacent to the core. The chemical gradient acts stabilising on the convection (Ledoux 1947), and thus this region defines a boundary for the g mode cavity as the trapping occurs in the μ gradient zone (Dziembowski et al. 1993).

Via expansion of the convective core, the presence of the oCZs produce a deviation from the otherwise monotonic increase in hydrogen profile (see Figure 5, solid lines) in the form of a plateau followed by an additional steep chemical gradient close to the outer boundary of the main chemical stratification. That gradient causes additional contribution to ∇_μ (see top panels of Figure 7) and thus to N (see lower panels of Figure 7) producing the bimodal distribution. As the star progresses through the main sequence, that feature persists, at least until TAMS, after which the star restructures its interior and develops a deep convective zone, which potentially can wash out this effect. The stability of the bimodal structure in N across various evolutionary stages is visible in Figure 7.

Fig. 7 Variation of radiative ∇rad (dotted lines) and chemical ∇μ gradients (dashed lines) on the top panels along with Lamb Sℓ=1 (dashed lines) and Brunt-Väisälä N frequencies (solid lines) on the bottom panels in a function of mass for 10 M☉ single (orange lines) and accretor (blue lines) models. The vertical lines correspond to features shown in Figure 5. We show the results at four different evolutionary stages, Xc=0.7,0.5,0.3 and 0.1.

3.3.2 g modes

The presence of two distinct peaks in the N profile leads to the trapping of specific modes within specific stellar regions, thereby modifying their periods and inducing deviations from the asymptotic period spacing (e.g. Dziembowski et al. 1993; Miglio et al. 2008; Wu et al. 2018). These effects were recently explored by Wagg et al. (2024) and Henneco et al. (2025) for intermediate-mass slowly pulsating B star (SPB) pulsator models and massive merger models, respectively. We revisit them here, as we see a similar behaviour for the g modes.

As the accretor evolves on the early main sequence, the period spacing closely follows the asymptotic value. However, once the convective core begins to recede and a chemical gradient forms at its boundary, the Δ P(P) pattern develops regular oscillations around the asymptotic value, typically within a P ∼1 day range. During the mass transfer episode, the spacing pattern becomes increasingly irregular, with no clear periodicity or coherence. This behaviour reflects rapid structural changes, the emergence of the off-centre convective zones, and the development of a highly structured Brunt-Väisälä profile shaped by composition gradients and complex mixing. In the final stages of the main sequence, the period spacing pattern stabilises and regains regularity with a modulated oscillatory behaviour, while the amplitude alternates between increasing and decreasing in a quasi-periodic manner. The amplitude varies in a cyclic fashion, resembling a beating pattern. This behaviour likely results from mode trapping between multiple internal layers, shaped by persistent chemical gradients near the convective core boundary, and reflects the increasingly stratified structure of the stellar interior as the star approaches the terminal age main sequence. Such a behaviour can be seen in Figure 8, in which we plot the dipole mode period spacing Δ P for the accretor and single star models, corresponding to four separate evolutionary cases: pre-MT case at X_c=0.7 and post-MT at X_c=0.5,0.3, and 0.1. With horizontal lines, we mark the asymptotic Δ P values (see Equation (8)), using solid lines for the accretor model and dashed lines for the equivalent single model. At X_c=0.7, both models exhibit nearly constant period spacings, with almost all modes lying close to the asymptotic value. This indicates a very small chemical composition gradient. Additionally, due to different initial masses, the asymptotic values themselves differ between the models, reflecting the dependence of the Brunt-Väisälä integral on the stellar size and internal structure. As evolution proceeds, after the mass transfer and mass equalization, the period spacing patterns begin to overlap, yet subtle differences persist between the models. Characteristic dips in the period spacing diagram start to emerge, associated with the development of a μ gradient region just outside the convective core. These dips are caused by mode trapping and result from the local increase in the BruntVäisälä frequency. Their presence and depth provide valuable asteroseismic diagnostics of the internal chemical stratification and the evolutionary state.

At the bottom panel of Figure 8, we also show the evolution of the ℓ=1 asymptotic Δ P_{ℓ, asym} values as a function of central hydrogen abundance X_c. For the accretor model, Δ P_{ℓ, asym} remains nearly constant at Δ P_{ℓ, asym} ∼0.115 d from the ZAMS up to X_c ≍ 0.5. Following the onset of mass accretion (I_acc, see Figure 1), Δ P_{ℓ, asym} increases rapidly, reaching a peak value of Δ P_{ℓ, asym} ∼0.16 d at X_c ≍ 0.6 (II_acc) and later stabilizes at Δ P_{ℓ, asym} ∼0.155 d (III acc). After mass transfer ends, the accretor evolves similarly to the single-star model, but with asymptotic Δ P_{ℓ, asym} values consistently higher by about 2%. This offset can be linked directly to changes in the Brunt-Väisälä frequency profile induced by mass accretion and structural readjustment.

We confirm the findings of Wagg et al. (2024), namely, that the period spacing pattern has a larger amplitude for the accretor compared to the single-star models. Additionally, for specific periods, the accretor's pattern can be out of phase with the single star's. This may be the most evident for X_c=0.1, in the second-to-last panel of Figure 8. The observed changes in period spacings between the accreting and single-star models indicate that mass transfer leaves detectable asteroseismic signatures. The magnitude of these differences, particularly in the g mode regime, suggests that they could be within the observational capabilities of missions such as Kepler (Koch et al. 2010), TESS (Ricker et al. 2015) or their successors, such as PLATO1 (Rauer et al. 2014), provided that a sufficient number of consecutive modes can be identified.

Fig. 8 Period spacing diagrams showing dipole g modes, for both the accretor model (orange line) and the single-star model (blue line). The horizontal lines mark the asymptotic value of the period spacing. We show results for only four different evolutionary stages, Xc= 0.7, 0.5, 0.3, and 0.1. The bottom panel shows the comparison of the ℓ=1 asymptotic period spacing evolution between the accretor and the single model. We mark points corresponding to key evolutionary changes of the accretor according to Figure 1.

3.3.3 p modes

In the case of p modes, we are primarily interested in low-degree, low to medium-order modes, since they are most likely to probe the deeper parts of the stellar interiors. For such modes, we calculated the accretor-to-single-star-frequency ratio, $v_{\ell ,n_{pg}}^{\mathrm{acc}}/v_{\ell ,n_{pg}}^{\sin }$, for an equal degree and order. We present the result in Figure 9 for four separate evolutionary cases, X_c=0.7,0.5,0.3, and 0.1. Similarly as in Figure 8, the top panel shows the pre-MT case, where we have a mismatch between the stellar masses. However, as shown in the following panels, distinct deviations between accretor and single star frequencies appear in the low- to intermediate-order regime. These are superimposed on a systematic offset observed in the asymptotic region, where the frequency ratio $v_{\ell ,n_{\mathrm{pg}}}^{\text{acc}}/v_{\ell ,n_{\mathrm{pg}}}^{\text{sin}}$ reaches approximately 1.0012 for X_c=0.5,0.996 for X_c=0.3, and 0.9925 for X_c=0.1. To illustrate this in more detail, we selected the representative model at X_c=0.1 and list the corresponding frequencies and their exact ratios in Table 1. As can be seen, the most affected frequencies are in the 5−15 d⁻¹ range, which is the typical observed range for β Cep stars (e.g. Stankov & Handler 2005). These non-asymptotic p modes are more sensitive to conditions in the near-core region, including the sound-speed gradient and the local stratification shaped by previous mass accretion. As a result, their frequency deviations provide direct seismic evidence of internal restructuring. To quantify this sensitivity and link the frequency shifts to structural differences, we computed the weight functions. This allowed us to identify the specific layers where the eigenfunctions contribute most significantly to the observed frequency deviations.

It was shown by Kawaler et al. (1985) and later revised by Townsend & Kawaler (2023) that the variational expression for the eigenfrequency of a mode can be written as $${\sigma }^2\left(\mathbf{y}\right)=\frac{{\int }_0^R[C(\mathbf{y},r)+\mathcal{N}(\mathbf{y},r)+\mathcal{G}(\mathbf{y},r)]\rho r^2\mathrm{d}r}{{\int }_0^R\mathcal{T}(\mathbf{y},r)\rho r^2\mathrm{d}r},$$(10) where $$C(\mathbf{y},r)=g^2\ell (\ell +1)S_{\ell }^{-2}{(y_2-y_3)}^2,$$(11) $$\mathcal{N}(\mathbf{y},r)=r^2{N}^2{y}_1^2,$$(12) $$\mathcal{G}(\mathbf{y},r)=-\frac{gr}{U}{[y_4+(\ell +1)y_3]}^2$$(13) are weight functions and y(r)={y₁(r), y₂(r), y₃(r), y₄(r)}= {ξ_r/r,(σ²r/g)(ξ_h/r), Φ^′/(gr),(1/g) d Φ^′/dr} is a vector of eigenfunctions introduced by Dziembowski (1971), corresponding to the radial displacement, the horizontal displacement, the gravitational potential perturbation and its derivative. Weight functions highlight the internal regions of the star that contribute to the mode's frequency, revealing its sensitivity to specific layers and hence establishing its main formation regions. These have been used extensively to examine the mode trapping properties in subdwarf B stars (Charpinet et al. 2000; Guyot et al. 2025).

The function in the denominator, $$\mathcal{T}(\mathbf{y},r)=r^2[y_1^2+\ell (\ell +1){\left(\frac{g}{r{\sigma }^2}\right)}^2{y}_2^2],$$(14) is proportional to the kinetic energy density and the U quantity, appearing in Equation (13), is the homology invariant (see e.g. Kippenhahn & Weigert 1990).

The frequency of a given mode, v_{ℓ, n}, can be expressed as $$v_{\ell ,n}^2=V(x=1),$$(15) where V is a cumulative integral of the weight functions W over a normalised radius x=r/R, $$V\left(x\right)={\int }_0^{x}W\left(x^{\text{'}}\right)\mathrm{d}{x}^{\text{'}}$$(16) and $$W\left(x\right)=\left(\frac{GM}{R^3}\right){\left(\frac{86400}{2\pi }\right)}^2·\frac{C\left(x\right)+\mathcal{N}\left(x\right)+\mathcal{G}\left(x\right)}{{\int }_0^R\mathcal{T}\left(r\right)\rho r^{2}dr}$$(17) The above is connected with Equation (10) via $$v=V{\left(1\right)}^{1/2}={\left(\frac{GM}{R^3}\right)}^{1/2}{\left(\frac{86400}{2\pi }\right)}^{-1}\sigma .$$(18)

The weight functions help to identify the regions in the star to which a given mode is most sensitive and their comparison allows us to understand how structural differences translate into changes in mode frequencies. To investigate how the internal structure affects the oscillation frequencies, we examined the behaviour of the weight functions W for selected modes. We focus on those with the largest frequency deviations between the two models at X_c=0.1, namely, the modes with (ℓ, n)=(2,1) and (ℓ, n)=(4,0) (see Figure 9 and Table 1). Additionally, to illustrate the response of gravity modes to the internal structure, we show the weight functions for selected g modes, namely, (ℓ, n)=(2,−2) and (ℓ, n)=(2,−3). We deliberately selected the low-order modes that can be trapped in the external peak region of the Brunt-Väisälä frequency (see Figure 7), and are therefore particularly sensitive to the stratification region left by the oCZs. As a result, their character may differ significantly from that of the corresponding modes in single-star models, reflecting the structural imprint of the mass transfer history. These weight functions are shown in Figure 10, together with their cumulative integrals V. The insets zoom into the region 0.07<r/R_☉<0.15 (corresponding to 1<m/M_☉<4.5).

The analysis of p mode frequency ratios further confirms that mass accretion subtly, but measurably alters the sound-speed profile of the star. These variations, although smaller than those observed in the g mode period spacings, offer complementary diagnostics sensitive to the layers close to the core and can thus provide additional constraints in asteroseismic modelling of mass-accreting systems. In particular, the differences in weight functions Δ W=W_accretor−W_single and the cumulative integrals Δ V=V_accretor−V_single indicate that the dominant contribution to the frequency shifts originates from the region 2<m/M_☉ <4, which coincides with the zone of strongest divergence in the sound-speed profiles c_s between the accretor and the single-star model (see Figure 6). The lower limit is, however, set by the size of the chemical gradient and thus dependent on the evolutionary phase we selected for study. This region encompasses the outer edge of the convective core and adjacent stratified layers shaped by the accretion history. The frequency deviations can thus be interpreted as a seismic signature of the structural imprint left by the mass transfer event. The magnitude of these p mode frequency shifts, typically of the order of 0.5%, correspond to absolute differences of ∼0.03−0.07 d⁻¹, in the 5−15 d⁻¹ range. These are well above the typical observational uncertainties for β Cep stars (∼0.003%, e.g., Jerzykiewicz et al. 2005; Burssens et al. 2023), and are therefore very likely to be detectable. However, as pointed out by Wagg et al. (2024), we are still most likely to confuse the observed frequencies originating from binary evolution with those of a much younger single star. This degeneracy can lead to a misinterpretation of the star's evolutionary state and age, potentially obscuring the effects of past mass transfer in asteroseismic modelling.

Fig. 9 Accretor and single-star frequency ratios for a given modes at different evolutionary stages, Xc=0.7,0.5,0.3, and 0.1.

Fig. 10 Weight functions, W, and their cumulative integrals, V, for the specific p modes (top panels) listed in Table 1 at Xc=0.1 and for low-order g modes, showing most distinct deviations between accretor (blue lines) and single-star (orange lines) frequencies.

3.3.4 High-order p modes

Even though high-order p modes are not observed in β Cephei pulsators, we explore the asymptotic properties of the p modes, as their asymptotic behaviour might offer a valuable theoretical framework for comparison.

It is well known that a local glitch such as a discontinuity in sound speed or density can induce periodic variations in the p mode frequencies (Gough 1990). For solar-like oscillators, the bottom of the convective envelope and the partial ionization zone of helium, have been inferred from high-frequency p modes (Roxburgh & Vorontsov 1994; Basu et al. 1994). The observables include the small frequency separation δ₀₂=v_{n, ℓ}−v_{n−1, ℓ+2} or the second-order frequency differences v_{n−1, l}−2 v_{n, l}+v_{n+1, l} (White et al. 2011; Verma et al. 2014).

From the asymptotic analysis, the small separation over large separation ratio, δ₀₂/Δ v, can be approximated as $${\delta }_{02}/\Delta v\approx \frac{(2\ell +3)}{4{\pi }^{2}v}(-\int \frac{\mathrm{d}c}{\mathrm{d}r}\frac{\mathrm{d}r}{c}+\frac{c\left(R\right)}{R}).$$(19) Thus, the imprints of mass transfer on the sound speed gradient near the core directly affect the magnitude of small separation. The above approximation has been refined and applied to red giants recently by Ong et al. (2025) and Reyes et al. (2025).

Figure 11 presents the ratio of the small to large frequency separations, computed using GYRE and from the leading-order asymptotic approximation for both the accretor and single-star models. As expected from theory, the asymptotic estimate of δ v₀₂/Δ v decreases monotonically with increasing radial order of n_pg, a trend that aligns well with the GYRE calculations in the high-order regime. Notably, the accretor model exhibits a slightly small to larger separation ratio compared to the single-star model. This difference can be attributed to a larger integral of the sound-speed gradient in the accretor model (see Equation (19)). Superimposed on the overall decline is a periodic modulation, which arises from the density discontinuity of the convective core boundary.

Following Roxburgh & Vorontsov (1994), the amplitude of this variation depends on the magnitude of the density discontinuity, and the period of the variational component can be approximated by $$\tilde{P}=(T/t_{\text{core}})\Delta v,$$(20) where t is the acoustic radius $$t\left(r\right)={\int }_0^r\frac{\mathrm{d}r}{c},$$(21) and T=t(r=R), t_core=t(r=r_core). Since the large frequency separation is given by Δ v ≍ 1 /(2 T), we can estimate the frequency of the variation as $$\tilde{v}=1/\tilde{P}=2t_{\text{core}}.$$(22)

We measured the periodic component in δ v₀₂/Δ v by recording the local minima and calculating the averaged period from high radial order p modes from GYRE (n_pg>30) and find a nice agreement with the expected frequency ṽ. We note that δ v₀₂/Δ v of the accretor model (solid lines) has a slightly larger period than that of the single model (dashed lines), indicating a larger T/t_core ratio.

These results highlight the sensitivity of high-order p modes to structural discontinuities introduced by mass transfer and demonstrate how asymptotic diagnostics, even if not directly observable in β Cephei stars, can provide valuable insights into the internal structure of post-interaction models.

Fig. 11 Ratio of the small frequency separation δ v02 derived from ℓ=0 and ℓ=2 p modes, to the large frequency separation, Δ v, is shown as a function of radial order. The asymptotic values are represented by dashed lines, while the solid lines correspond to the results obtained from GYRE. Both the single-star (orange lines) and the accretor (blue lines) models are included for comparison.

4 Summary

In this study, we investigate the structural and asteroseismic consequences of mass accretion onto a working example of a main sequence star with a final mass of 10 M_☉, representative of a typical β Cephei pulsator. We compare a single-star model and an accretor model that has undergone conservative mass transfer, focusing on key internal diagnostics such as density, chemical composition, sound speed, Brunt-Väisälä frequency, and g mode and low-order p mode frequencies.

Mass accretion significantly modifies the internal structure and asteroseismic properties of the β Cep accretor. The most prominent structural change is the formation of a local bump in the density profile just above the convective core. This feature arises during rejuvenation when the core expands and engulfs hydrogen-rich layers, creating a sharp μ gradient and a nonmonotonic density bump. The bump persists throughout the main sequence evolution and corresponds to a lasting chemical discontinuity visible at the outer-boundary of the μ gradient region.

These structural modifications are reflected in the BruntVäisälä frequency profiles. While the single-star model shows a smooth Brunt-Väisälä profile shaped solely by core retreat and composition gradients, the accretor exhibits a bi-modal profile that results in additional mode trapping.

The g mode period spacing diagram offers a sensitive asteroseismic probe of these structural changes. Before accretion, both models show near-constant Δ P values, closely matching the asymptotic expectation. However, during and after mass transfer, the period spacing of the accretor becomes irregular due to rapid internal restructuring. Once the star reaches a post-accretion state, the Δ P patterns regain coherence but develop characteristic modulations. These features result from mode trapping in the layered structure around the convective core and trace the evolving composition gradients. The quasi-periodic modulations in period spacing, including beating-like behaviour, serve as direct signatures of internal stratification shaped by mass accretion.

Low- and intermediate-order p modes are likewise affected. Frequency ratios v^acc/v^sin reveal systematic shifts up to ∼0.5%, corresponding to absolute differences of 0.03−0.07 d⁻¹ in the observable 5−15 d⁻¹ range for β Cephei stars. These shifts arise primarily from differences in the sound-speed profile near the core due to the local density bump. Weight function analysis confirms that these layers dominate the contribution to the observed frequency deviations.

Finally, we detect changes in the δ v₀₂/Δ v patterns. The accretor shows a slightly longer modulation period, indicating a larger T/t_core ratio. This further confirms the seismic signature of a structurally more extended core boundary region caused by prior mass accretion.

In summary, we demonstrate that mass accretion leaves measurable imprints on the internal structure of mass-accreting stars, especially in density, chemical gradients, and sound-speed profiles. These translate into observable differences in both g and p mode pulsations. Period spacing modulations, mode trapping features, and systematic frequency shifts provide robust asteroseismic diagnostics of prior mass transfer events. Such signatures are within the sensitivity of current high-precision photometry, making asteroseismology a powerful tool to probe the evolutionary history of massive stars in binary systems. This study, although focused strictly on a high-mass stellar model, provides a general methodology and a clear prediction of the seismic signatures expected from mass accretion. The approach and results presented here are supported by independent modelling efforts, such as those of Wagg et al. (2024) and Henneco et al. (2024, 2025), further reinforcing the diagnostic potential of pulsations for unveiling past binary interaction.

Future work ought to concentrate on the effect of the adopted physics (in particular, rotation and the associated mixing) on the asteroseismic properties of post-mass-transfer stars. Differential rotation and enhanced mixing can substantially alter the internal chemical stratification and thus modify the driving and damping of pulsation modes, leading to different instability domains compared to single-star counterparts (see the discussion on mode driving in Miszuda et al. 2021 and Wagg et al. 2024). This could directly impact which modes are preferentially excited and how their frequencies are distributed. In addition, employing machine-learning techniques to systematically identify discrepancies between dynamically inferred stellar parameters and those derived from seismic modelling could provide a powerful diagnostic tool to distinguish genuine binary interaction products from apparent single stars.

Data availability

We make all files needed to recreate our MESAbinary and GYRE results publicly available at Zenodo: 10.5281/zenodo.15829740

Acknowledgements

This work was supported by the (AM) Polish National Science Centre (NCN), grant number 2021/43/B/ST9/02972, (RHDT) NASA grants 80NSSC24K0895 and 80NSSC23K1517 and (ZG) European Research Council (ERC) under the Horizon Europe program (Synergy Grant agreement 101071505: 4D-STAR). While partially funded by the European Union, views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. AM would like to express his gratitude to G. Handler for his support and valuable comments on this manuscript. Calculations have been carried out using resources provided by Wrockaw Centre for Networking and Supercomputing (http://wcss.pl), grant no. 265. Used software: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), GYRE (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020; Sun et al. 2023), pyMESAREADER (https://billwolf.space/py_mesa_reader/index.html).

Appendix A MESA input physics

The MESA code builds upon the efforts of many researchers who advanced our understanding of physics and relies on a variety of input microphysics data. The MESA EOS is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010) EOSes. Radiative opacities are primarily from the OPAL project (Iglesias & Rogers 1993, 1996), with data for lower temperatures from Ferguson et al. (2005) and data for high temperatures, dominated by Compton-scattering from Buchler & Yueh (1976). Electron conduction opacities are from Cassisi et al. (2007). Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010) plus additional tabulated weak reaction rates from Fuller et al. (1985), Oda et al. (1994) and Langanke & Martínez-Pinedo (2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996). The MESA-BINARY module allows for the construction of a binary model and the simultaneous evolution of its components, taking into account several important interactions between them. In particular, this module incorporates angular momentum evolution due to mass transfer. Roche lobe radii in binary systems are computed using the fit of Eggleton (1983). mass transfer rates in Roche lobe overflowing binary systems are determined following the prescriptions of Ritter (1988) and Kolb & Ritter (1990).

References

All Tables

All Figures

Fig. 1 The HR diagram with the evolutionary tracks of the binary system components. The donor star is shown with the purple line and the accretor with the gray line, both plotted alongside single-star reference models for different masses (dashed lines). The mass transfer/accretion rate, log |·M|, is overplotted along the tracks using colour coding, as indicated by the colour bar. The key evolutionary changes are marked: the donor fills its Roche lobe and begins mass transfer (Idon), the accretor starts to gain mass and departs from its original track (Iacc), full mass ratio reversal when the accretor reaches 10 M☉ and the following detachment of the system (IIacc and IIdon), thermal readjustment of the accretor (IIIacc), and the accretor's single-star evolution until the end of the main sequence (I Vacc).

In the text

	Fig. 2 Evolution of mass of the convective core, Mcc, in a function of central hydrogen abundance, Xc, for the single-star (orange lines) and the accretor (blue lines) models. We mark key evolutionary changes of the accretor according to Figure 1.
In the text

Fig. 3 Mixing profiles, D, during the evolution of the accretor. Left: evolution from Xc=0.7 to 0.0 in steps of 0.1. Right: profiles at selected evolutionary stages, before, during and after the RLOF. The shaded areas indicate the regions over which the mass varies, with the current mass indicated by dashed lines. The black lines, as denoted on the right-hand ordinates of the right panels show the corresponding profiles of hydrogen abundances.

In the text

	Fig. 4 Density, log ρ, profiles at selected central hydrogen abundances, Xc, for the single-star (orange lines) and the accretor (blue lines) models. The local bump discussed in the text is clearly visible in the binary models.
In the text

	Fig. 5 Hydrogen abundances, X, in a function of a mass for different evolutionary stages of 10 M☉ single (orange lines) and accretor (blue lines) models. The vertical lines mark features corresponding to Figure 7. For clarity, we only show a fraction of the inner mass coordinate.
In the text

	Fig. 6 Sound speed, cs, profiles at selected central hydrogen abundances, Xc, for the single-star (orange lines) and the accretor (blue lines) models.
In the text

Fig. 7 Variation of radiative ∇rad (dotted lines) and chemical ∇μ gradients (dashed lines) on the top panels along with Lamb Sℓ=1 (dashed lines) and Brunt-Väisälä N frequencies (solid lines) on the bottom panels in a function of mass for 10 M☉ single (orange lines) and accretor (blue lines) models. The vertical lines correspond to features shown in Figure 5. We show the results at four different evolutionary stages, Xc=0.7,0.5,0.3 and 0.1.

In the text

Fig. 8 Period spacing diagrams showing dipole g modes, for both the accretor model (orange line) and the single-star model (blue line). The horizontal lines mark the asymptotic value of the period spacing. We show results for only four different evolutionary stages, Xc= 0.7, 0.5, 0.3, and 0.1. The bottom panel shows the comparison of the ℓ=1 asymptotic period spacing evolution between the accretor and the single model. We mark points corresponding to key evolutionary changes of the accretor according to Figure 1.

In the text

	Fig. 9 Accretor and single-star frequency ratios for a given modes at different evolutionary stages, Xc=0.7,0.5,0.3, and 0.1.
In the text

	Fig. 10 Weight functions, W, and their cumulative integrals, V, for the specific p modes (top panels) listed in Table 1 at Xc=0.1 and for low-order g modes, showing most distinct deviations between accretor (blue lines) and single-star (orange lines) frequencies.
In the text

	Fig. 11 Ratio of the small frequency separation δ v02 derived from ℓ=0 and ℓ=2 p modes, to the large frequency separation, Δ v, is shown as a function of radial order. The asymptotic values are represented by dashed lines, while the solid lines correspond to the results obtained from GYRE. Both the single-star (orange lines) and the accretor (blue lines) models are included for comparison.
In the text