© The Authors 2025

1. Introduction

With the discovery of thousands of exoplanets in the last few decades (5678 as of June 2024; Nasa Exoplanet Archive¹), the need to precisely characterise their host stars has become crucial in order to accurately derive the properties of the systems (e.g. Adibekyan et al. 2018) and understand their formation and evolution. From their birth, stars shape the evolution of exoplanetary systems. At the same time, host stars’ properties may be significantly influenced by the presence of interacting close-by planetary companions, potentially resulting in anomalous rotational periods (e.g. Privitera et al. 2016; Ilić et al. 2024) and magnetic activity, with brighter X-ray luminosity emissions (e.g. Shkolnik et al. 2003; Poppenhaeger & Wolk 2014; Pillitteri et al. 2022; Ilić et al. 2023) (among others). Understanding how the evolution of stars is impacted by the interactions with their planets is thus fundamental to disentangling their intrinsic evolution from SPI-induced phenomena on the one hand and to helping determine properties that could not be derived otherwise (e.g. tracing back their rotational history) on the other.

In this context, GJ 504 (a.k.a. HD 115383, TIC 397587084) is a promising candidate for a star that underwent strong star-planet interactions. GJ 504 is an isolated G0 spectral type star that is slightly more massive than the Sun (M_⋆ ∼ 1.22 M_⊙) and hosts a directly imaged substellar companion at a projected distance of ∼43.5 AU (Kuzuhara et al. 2013). Because of the significant uncertainties in determining GJ 504’s evolutionary state, with age estimations varying from hundreds of megayears to several gigayears (Valenti & Fischer 2005; Takeda et al. 2007; Holmberg et al. 2009; da Silva et al. 2012; Kuzuhara et al. 2013; Fuhrmann & Chini 2015; D’Orazi et al. 2017; Di Mauro et al. 2022), the nature of the companion is still hotly debated, and very different values for its mass have been proposed in the literature (1 ≲ M_c/M_J ≲ 25) (Kuzuhara et al. 2013; Fuhrmann & Chini 2015).

Establishing the age of isolated stellar objects is particularly hard. Different methods could be used (Soderblom et al. 2014), including gyrochronology and activity indicators, for the age-rotation-activity relation (Skumanich 1972; Barnes 2007) or comparison of classical spectroscopic parameters with model isochrones (e.g. Pont & Eyer 2004). Kuzuhara et al. (2013) estimated the age of GJ 504 via several methods, but they considered the one derived from gyrochronology and activity indicators, based on P_rot = 3.329 days (Donahue et al. 1996; Messina et al. 2003) and chromospheric activity indices (Ca II H and K lines, log(R_HK′) = −4.45), as the most likely one, giving the value $\mathrm{t}={160}_{-60}^{+350}\mathrm{Myr}$. A few years later, Fuhrmann & Chini (2015) carried out a detailed analysis of high-resolution, high-quality spectra, which led to a major revision of the stellar gravity parameter log(g) = 4.23 ± 0.10, which is about 8% smaller than the value in Kuzuhara et al. (2013). From the isochrone-fitting method, they found an age much closer to the one of the Sun ($\mathrm{t}=4.5_{-1}^{+2}\mathrm{Gyr}$) (Fuhrmann & Chini 2015). To explain the short rotational period, intense chromospheric activity, and bright X-ray luminosity (Log(L_X/L_Bol) = − 4.42, Voges et al. 1999; Wright et al. 2011) and to reconcile these indicators with the isochronal age, Fuhrmann & Chini (2015) invoked the engulfment² of a close-by planetary companion that would have spun up the stellar surface and enhanced its activity levels (Oetjens et al. 2020). D’Orazi et al. (2017) reassessed the fundamental properties of GJ 504, finding that a comparison of their spectroscopic parameters with isochrones provided an age between 1.8 and 3.5 Gyr, with a most probable value of t ≈ 2.5 Gyr. They also envisaged a possible engulfment scenario to reconcile the different age indicators and tested this hypothesis by means of a tidal-evolution code. They found that the engulfment of a hot Jupiter, with an initial mass no larger than ≈3 M_J and initial orbital distance of ≈0.03 AU, could be a very likely scenario. Bonnefoy et al. (2018) revisited the system by means of high-contrast imaging and interferometric and radial-velocity observations. From their analysis, they retrieved an interferometric radius for the host star of R = (1.35 ± 0.04) R_⊙, which is compatible with two isochronal ages: (21 ± 2) Myr and (4.0 ± 1.8) Gyr. The mass of the substellar companion would correspond to $1.3_{-0.3}^{+0.6}{\mathrm{M}}_{\mathrm{J}}$ or ${23}_9^{+10}{\mathrm{M}}_{\mathrm{J}}$, respectively. The authors also revised the almost pole-on line-of-sight stellar-rotation-axis inclination, which is $18.6_{-3.8}^{+4.3}$ degrees or $162.4_{-4.3}^{+3.8}$ degrees. They also excluded the presence of additional objects (with 90% probability) more massive than 2.5 and 30 M_J with semi-major axes in the range of 0.01 − 80 AU for the young and old isochronal ages, respectively. More recently, Di Mauro et al. (2022) attempted to employ asteroseismic techniques on the observational data collected by the Transiting Exoplanet Survey Satellite space mission (TESS, Ricker et al. 2014) to accurately characterise GJ 504. Unfortunately, the non-detection of solar-like oscillations, which Di Mauro et al. (2022) ascribed to the high level of the stellar magnetic activity, hindered their analysis. Nevertheless, the results deduced by TESS photometric data – supported by the Mount Wilson Observatory long-term campaign spanning nearly 30 years – and modelling procedures allowed the authors to refine the fundamental parameters of GJ 504. Among these, they derived an age ≤2.6 Gyr, which is in agreement with previous findings by Kuzuhara et al. (2013) and D’Orazi et al. (2017) within the quoted uncertainties. From the TESS light curves, they also identified a clear modulation corresponding to a rotational period of P_rot = 3.4 days, which confirms the average value used in Kuzuhara et al. (2013). Finally, from Mount Wilson Observatory data they obtained the detection of a main magnetic cycle of 11.97 years, which, together with the relatively short rotational period, located GJ 504 before the transition to the weakened magnetic braking regime, as theorised by van Saders et al. (2016). During this transition the large-scale magnetic field would fail to efficiently brake the stellar surface.

The aim of our work is to study the peculiar properties of GJ 504 by testing the impact of a putative, close-by planet’s inward migration on the evolution of the stellar rotational period and X-ray luminosity. For this reason, we first looked for the optimal stellar parameters representative of GJ 504 by means of a minimisation procedure based on classical spectroscopic and interferometric parameters and correspondingly computed best-fit stellar models (see Sect. 2.1). Subsequently, we coupled the stellar models to our SPI code (see Sect. 2.2), in which the evolution of the stellar surface-rotation rate and X-ray luminosity is computed simultaneously to the evolution of the orbit of a close-by planetary companion; this is driven by the dissipation of tides within the host star (see Sect. 4). We envisage two evolutionary scenarios:

• The star-without-close-by-planet scenario (Sect. 4.1), in which no massive or close-by companions affect the evolution of the host star

• The star-with-close-by-planet scenario (Sect. 4.2), in which a putative, close-by planet strongly impacts the evolution of its host star³

We thus compare the results of the simulations with observational data for the rotational period and X-ray luminosity. For the latter one, in addition to the ROSAT (Truemper 1982) data, we give a detailed analysis of the X-ray data taken with the ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021) on board the Russian Spektrum-Roentgen-Gamma mission (SRG; Sunyaev et al. 2021) (see Sect. 3) and compare it with previous results found in Foster et al. (2022). We searched for an X-ray counterpart of GJ 504 among the sources listed in the five all-sky surveys of eROSITA (eRASS; Merloni et al. 2024). The first four surveys are completed and count six months of observation, while the last survey was suspended after almost two months because the science operations of the instrument were paused. Finally, in Sect. 5 we draw our conclusions.

2. Method and physics

To carry out our study on GJ 504, we proceeded according to the following steps. Firstly, we derived the optimal parameters for the star (initial mass, radius, chemical composition, etc.) by means of a minimisation technique based on classical spectroscopic and interferometric parameters (unfortunately asteroseismic indicators are not available for this star); secondly, best-fit stellar models were computed based on the optimal parameters found before; finally, the stellar models are coupled to our SPI code in which the evolution of the star is computed by envisaging two potential scenarios. These are the star-without-close-by-planet and the star-with-close-by-planet scenarios.

2.1. Optimal parameter search and stellar model

To compute stellar models for GJ 504, we started by deriving optimal stellar parameters by means of a two-step, global, and local-minimisation procedure. For the global minimisation, we employed the SPInS software (Lebreton & Reese 2020), which is based on a Markov chain Monte Carlo (MCMC) approach and Bayesian statistics, to derive probability distribution functions for stellar parameters. For the local minimisation, we used the Levenberg-Marquardt algorithm (e.g. Miglio & Montalbán 2005; Farnir et al. 2020).

In both minimisation procedures, we used models computed with the Liège stellar evolution code (CLES) (e.g. Scuflaire et al. 2008). For the setup of input physics, we considered the FreeEOS (Irwin 2012) equation of state, AGSS09 (Asplund et al. 2009) abundances, and OPAL (Iglesias & Rogers 1996) opacity tables for solar mixture. The classical mixing-length theory was applied for convection, with a solar calibrated value α_MLT = 2.01. For the outer boundary conditions, we used Vernazza et al. (1981).

In this context, we aimed to determine the values of four free parameters (M_⋆, age, X_0, and Z₀) using four observational constraints taken from Di Mauro et al. (2022) (see references therein): Teff = (6205 ± 20) K, R = (1.35 ± 0.04) R_⊙, [Fe/H] = 0.22 ± 0.04, and log(g) = 4.29 ± 0.07.

The results obtained from the first step of modelling with SPInS are M_⋆ = (1.30 ± 0.05) M_⊙, X₀ = 0.70 ± 0.02, Z₀ = 0.025 ± 0.002, and age =(1.81 ± 0.46) Gyr. While this solution is limited to the fixed input physics used for the computation of the model grid, in the second step the local minimisation procedure allows us to further explore the parameter space thanks to the use of the Levenberg-Marquardt algorithm. With the computation of models on the fly, it is possible to investigate the impact of changing the input physics on the stellar parameters. In particular, we tested the effect of including a moderate amount of overshooting (α_Ov = 0.1 H_p, with H_p being the pressure scale height). We also changed the outer boundary conditions, using the Eddington T(τ) relations, and the corresponding solar calibrated value for the mixing length α_MLT = 1.8. The optimal parameters derived by means of this procedure are M_⋆ = (1.29 ± 0.21) M_⊙, X₀ = 0.70 ± 0.09, Z₀ = 0.025 ± 0.002, and age =(2.11 ± 1.75) Gyr.

Both solutions agree very closely on the mean values, but from the local minimisation procedure we derived larger uncertainties, especially for the age. This is due not only to the broader input physics used in the Levenberg-Marquardt algorithm, but also to the intrinsic differences in comparison to SPInS for the computation of the uncertainties.

In general, these solutions appear to be in good agreement with the ones found in Di Mauro et al. (2022), although a smaller uncertainty on the age estimation was found in our case. According to our analysis, an age as young as ≈200 Myr seems to be disfavoured for this star. Nevertheless, it is worth stressing that only seismic constraints would likely allow us to derive more precise age estimates (Soderblom 2010).

We thus proceeded by computing the evolutionary sequence of GJ 504 using CLES, for which the derived optimal parameters (M_⋆ = 1.29 M_⊙, X₀ = 0.70, Z₀ = 0.025) were provided as input and constraints. According to this model, GJ 504 is a main-sequence (MS) star, with a central abundance of hydrogen of X_c = 0.37, mass of the convective envelope of M_env = 2.5 × 10⁻³ M_⋆, radius of the convective envelope of R_env = 0.18 R_⋆, and an inner convective core with M_cc = 1.3 × 10⁻² M_⋆ and M_cc = 4.3 × 10⁻² R_⋆. In Fig. 1, we show the evolutionary track corresponding to the best-fit model. For a detailed description of the physics included in CLES, we refer the interested reader to Scuflaire et al. (2008).

Fig. 1. Hertzsprung-Russell diagram of GJ 504 as obtained from best-fit stellar model (Teff = 6202 K, log(g) = 4.29, R = 1.35 R⊙, L/L⊙ = 2.43) computed with CLES.

As mentioned above, the computed evolutionary sequence of GJ 504 is provided to our SPI code to study its past evolution. Given the significant uncertainties derived on GJ 504’s parameters, in the following we discuss our results considering the broadest age interval ((0.36 − 3.86) Gyr) and show both average values and their corresponding uncertainties in the figures for comparison. While different age values would correspond to different stellar masses and chemical compositions, considering the physics included in our SPI code we can reasonably assume that our conclusion remains valid within the 4D space (M_⋆, age, X_0, and Z₀) of variation for our best-fit stellar model.

2.2. Star-planet-interaction code

The stellar models computed for GJ 504 are provided as input to our SPI code. In this code it is possible to compute the rotational evolution of the star while accounting for two main types of star-planet interaction: gravitational-tidal and radiative (e.g. Cuntz et al. 2000; Vidotto 2020; Strugarek 2024). While a more detailed description of the physics implemented in the code is provided in Privitera et al. (2016), Rao et al. (2018), Pezzotti et al. (2021), Fellay et al. (2023), in the following we just recall the equations of interest for the evolution of a putative planet orbiting closely around GJ 504.

2.2.1. Host-star surface-rotation evolution

The evolution of the stellar surface-rotation rate is computed assuming that the host-star rotates as a solid body on the pre-main-sequence (PMS) and MS phases (Rao et al. 2021). Given the relative shallowness of the convective envelope characterising GJ 504 and the fact that the differential rotation in this region is less significant compared to the almost flat profile of internal rotation in the core, this is a reasonable assumption, which has been supported by seismic analyses conducted on MS, solar-like stars (García et al. 2007; Nielsen et al. 2015; Benomar et al. 2015, 2018; Bétrisey et al. 2023) and γDor stars (Saio et al. 2021). The equation describing the variation of the stellar angular momentum reads as

$$\begin{array}{c}\hfill {\dot{L}}_{\star }={\dot{L}}_{\mathrm{wb}}+{\dot{L}}_{\mathrm{tides}},\end{array}$$(1)

where ${\dot{\mathrm{L}}}_{\mathrm{wb}}$ accounts for the rate of change of angular momentum due to magnetised winds (Matt et al. 2015, 2019), with equation

$$\begin{array}{c}\hfill {\dot{L}}_{\mathrm{wb}}=\{\begin{array}{cc}-{\mathrm{T}}_{\odot }{\left(\frac{{\mathrm{R}}_{\star }}{{\mathrm{R}}_{\odot }}\right)}^{3.1}{\left(\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{0.5}{\left(\frac{{\tau }_{\mathrm{conv}}}{{\tau }_{\mathrm{conv}\odot }}\right)}^{\mathrm{p}}{\left(\frac{{\mathrm{\Omega }}_{\star }}{{\mathrm{\Omega }}_{\odot }}\right)}^{\mathrm{p}+1},\hfill & \mathrm{Ro}>\frac{{\mathrm{Ro}}_{\odot }}{\chi },\hfill \\ -{\mathrm{T}}_{\odot }{\left(\frac{{\mathrm{R}}_{\star }}{{\mathrm{R}}_{\odot }}\right)}^{3.1}{\left(\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{0.5}{\chi }^{\mathrm{p}}\left(\frac{{\mathrm{\Omega }}_{\star }}{{\mathrm{\Omega }}_{\odot }}\right),\hfill & \mathrm{Ro}\le \frac{{\mathrm{Ro}}_{\odot }}{\chi },\hfill \end{array}\end{array}$$(2)

where R_⊙ and M_⊙ are the radius and mass of the Sun; R_⋆ and M_⋆ are the radius and mass of the considered star; τ_conv is the convective turnover timescale (see Eq. (6)); Ro is the stellar Rossby number, defined as the ratio between the star’s rotational period (P_rot) and τ_conv; and Ro_⊙ is the solar Rossby number. The braking constant T_⊙ = 8 × 10³⁰ erg is calibrated to reproduce the surface-rotation rate of the Sun, and the coefficient p is taken as equal to 2.1. The quantity χ ≡ Ro_⊙/Ro_sat indicates the critical rotation rate for stars with a given τ_conv/τ_conv⊙, defining the transition from the saturated to unsaturated regime. In this work, χ is considered equal to 13 (Pezzotti et al. 2021). For what concerns the braking of the stellar surface by magnetised stellar winds, it is possible to account for the onset of a ‘weakened-magnetic-braking’ regime (van Saders et al. 2016) when the stellar Rossby number (Ro) is equal to a critical value of Ro_cri = 0.92 Ro_⊙ (Metcalfe et al. 2024): once a star enters this regime, the efficiency of the braking at its surface is hampered by a shift in the morphology of the magnetic field, from larger to smaller spatial scales (Réville et al. 2015; Garraffo et al. 2016), and (or) by an abrupt change in the mass-loss rate (Ó Fionnagáin & Vidotto 2018). In our code, we account for this effect by simply turning off the magnetic braking component in the computation of the surface angular momentum evolution for Rossby numbers larger than Ro_cri. The component ${\dot{L}}_{\mathrm{tides}}$ in Eq. 1 accounts for the exchange of angular momentum between the star and the planetary orbit due to tidal dissipation (Rao et al. 2021), and its equations read as

$$\begin{array}{c}\hfill {\dot{L}}_{\mathrm{tides}}=-\left[\frac{1}{2}{M}_{\mathit{pl}}{\left(\frac{\dot{a}}{a}\right)}_{\mathrm{tides}}\right]\times \sqrt{\mathrm{G}({\mathrm{M}}_{\star }+{\mathrm{M}}_{\mathrm{pl}})\mathrm{a}},\end{array}$$(3)

where the term ${\left(\frac{\dot{\mathrm{a}}}{\mathrm{a}}\right)}_{\mathrm{tides}}$ indicates the net variation of the planetary orbital distance (a) due to tidal dissipation (see Eq. (5), (7)), M_pl is the planetary mass, and G is the universal gravitational constant.

2.2.2. Tidal interaction

In our code, star-planet interactions related to the dissipation of equilibrium and dynamical tides within the stellar convective envelope are considered, which can lead to planetary migration through the transfer of angular momentum from the planetary orbit to the stellar surface. In the hypothesis that the planet is on a circular, coplanar orbit around the host star, the net change of its orbital distance can be written as

$$\begin{array}{c}\hfill \left(\frac{\dot{\mathrm{a}}}{\mathrm{a}}\right)=-\frac{{\dot{\mathrm{m}}}_{\mathrm{pl}}}{{\mathrm{m}}_{\mathrm{pl}}+{\mathrm{M}}_{\star }}+{\left(\frac{\dot{\mathrm{a}}}{\mathrm{a}}\right)}_{\mathrm{eq}}+{\left(\frac{\dot{\mathrm{a}}}{\mathrm{a}}\right)}_{\mathrm{dyn}},\end{array}$$(4)

where the first term is related to planetary mass loss due to photoevaporation, assuming that the mass is lost in space and does not contribute to the exchange of angular momentum between the star and the planetary orbit. The last two terms refer to the contribution of equilibrium and dynamical tides, respectively. Equilibrium tides are most efficiently dissipated in low-mass stars’ convective envelopes due to turbulent friction induced by convection, while their damping is less efficient in other regions (Zahn 1966, 1977). In our approach, we therefore only account for the dissipation of equilibrium tides in convective envelopes, and the equation accounting for their contribution reads as

$$\begin{array}{c}\hfill {\left(\frac{\dot{\mathrm{a}}}{\mathrm{a}}\right)}_{\mathrm{eq}}=\frac{\mathrm{f}}{{\tau }_{\mathrm{conv}}}\frac{{\mathrm{M}}_{\mathrm{env}}}{{\mathrm{M}}_{\star }}\mathrm{q}(1+\mathrm{q}){\left(\frac{{\mathrm{R}}_{\star }}{\mathrm{a}}\right)}^8(\frac{{\mathrm{\Omega }}_{\star }}{{\omega }_{\mathrm{pl}}}-1),\end{array}$$(5)

where a is the orbital distance, ȧ is its derivative with respect to time, M_env is the mass contained in the stellar convective envelope, q is the ratio between planetary and stellar mass, Ω_⋆ is the host star’s surface-rotation rate, ω_pl is the orbital frequency, and R_⋆ is the stellar radius. The term τ_conv is the convective turnover timescale, computed as in Rasio et al. (1996):

$$\begin{array}{c}\hfill {\tau }_{\mathrm{conv}}={\left(\frac{{\mathrm{M}}_{\mathrm{env}}{({\mathrm{R}}_{\star }-{\mathrm{R}}_{\mathrm{b}})}^2}{3{\mathrm{L}}_{\star }}\right)}^{1/3}.\end{array}$$(6)

where R_b is the radius at the bottom of the convective envelope and L_⋆ is the bolometric luminosity of the host star. The term f in Eq. (5) is a numerical factor obtained from integrating the viscous dissipation of the tidal energy across the convective zone (Villaver & Livio 2009), which is f = (P_orb/2τ_conv)² (Goldreich & Nicholson 1977) when τ_conv < P_orb/2, otherwise f = 1.

Contrarily to equilibrium tides, dynamical tides might be efficiently dissipated in both convective and radiative regions, depending on the stellar and planetary companion properties. In convective regions, they are excited in the form of inertial waves driven by the Coriolis force whenever the tidal frequency (ω_t) ranges between [−2Ω_⋆,2Ω_⋆] (Ogilvie & Lin 2007). We account for the impact of dynamical tides in convective envelopes in the form of a frequency-averaged tidal dissipation of inertial waves (Ogilvie 2013; Mathis 2015; Bolmont & Mathis 2016). In particular, we assume a two-layer model (core-envelope) in which each of them is characterised by a uniform density (Mathis 2015). If the planet is on a circular-coplanar orbit, dynamical tides are active whenever ω_pl < 2Ω_⋆ (Ogilvie & Lin 2007), where ω_pl is the orbital frequency of the planet. The equation describing the contribution of this type of tides on the evolution of the planetary orbital distance reads as

$$\begin{array}{c}\hfill {(\dot{\mathrm{a}}/\mathrm{a})}_{\mathrm{dyn}}=\left(\frac{9}{2\mathrm{Q}{\prime }_{\mathrm{d}}}\right)\mathrm{q}{\omega }_{\mathrm{pl}}{\left(\frac{{\mathrm{R}}_{\star }}{\mathrm{a}}\right)}^5\frac{({\mathrm{\Omega }}_{\star }-{\omega }_{\mathrm{pl}})}{\mid {\mathrm{\Omega }}_{\star }-{\omega }_{\mathrm{pl}}\mid },\end{array}$$(7)

with Q′_d = 3/(2D_ω) (modified tidal dissipation factor) and ${\mathrm{D}}_{\omega }={\mathrm{D}}_{0\omega }{\mathrm{D}}_{1\omega }{\mathrm{D}}_{2\omega }^{-2}$. The ‘D’ terms are defined as follows:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathrm{D}}_{0\omega }=\frac{100\pi }{63}{ϵ}^2\frac{{\alpha }^5}{1-{\alpha }^5}{(1-\gamma )}^2,\hfill \\ {\mathrm{D}}_{1\omega }={(1-\alpha )}^4{(1+2\alpha +3{\alpha }^2+\frac{3}{2}{\alpha }^3)}^2,\hfill \\ {\mathrm{D}}_{2\omega }=1+\frac{3}{2}\gamma +\frac{5}{2\gamma }(1+\frac{\gamma }{2}-\frac{3{\gamma }^2}{2}){\alpha }^3-\frac{9}{4}(1-\gamma ){\alpha }^5,\hfill \end{array}\end{array}$$(8)

where α = R_b/R_⋆, β = M_b/M_⋆, $\gamma =\frac{{\alpha }^3(1-\beta )}{\beta (1-{\alpha }^3)}$, and $ϵ={\mathrm{\Omega }}_{\star }/\sqrt{{\mathrm{GM}}_{\star }/{\mathrm{R}}_{\star }^3}$. The terms M_b and R_b are the mass and the radius of the radiative core, considered as the region of the star below the base of the convective envelope. The term D_ω is the frequency-averaged tidal dissipation (Ogilvie 2013). We notice that assuming this schematic stratification for the stellar structure, together with the use of a frequency-averaged dissipation rate for the dynamical tide, adds a certain degree of uncertainty to our results. Nevertheless, while this approach suffers from some schematic simplifications, it has the advantage of providing us with relevant orders of magnitude of tidal dissipation rates accounting for the evolution of the structural and rotational parameters of the host star (Mathis 2015; Bolmont & Mathis 2016; Rao et al. 2018; Barker 2020).

In general, for the formalism treated above, tides widen (respectively shrink) the planetary orbit when it is beyond (respectively inside) the corotation radius, which is defined as the distance at which the orbital and host star rotational periods are equal, namely:

$$\begin{array}{c}\hfill {\mathrm{a}}_{\mathrm{cor}}={[\mathrm{G}({\mathrm{M}}_{\star }+{\mathrm{M}}_{\mathrm{pl}})/{\mathrm{\Omega }}_{\star }^2]}^{\frac{1}{3}}.\end{array}$$(9)

Here, G is the universal gravitational constant, and M_⋆ is the mass of the host star. Concerning the formalism used in this work for dynamical tides in the stellar convective envelope, these are efficiently excited in the form of inertial waves when the planetary orbital distance is larger than a minimum critical value, which is defined as ${\mathrm{a}}_{\min }=4^{-\frac{1}{3}}\times {\mathrm{a}}_{\mathrm{cor}}$ (Ogilvie & Lin 2007).

In the case of early G-type and late F-type stars, a convective core is present in addition to the envelope during the MS and is separated by a radiative region. Nevertheless, dynamical tides excited in convective cores are weakly dissipated. In this case indeed, inertial waves would propagate in a fully convective sphere (the inner core), without the possibility to reflect and lead to the formation of shared wave attractors, which may induce strong dissipation (Ogilvie & Lin 2007, 2004). The dissipation of inertial waves in convective cores is thus not considered in our SPI code.

In general, dynamical tides might also be excited in stellar radiative regions in the form of gravity waves and cause the migration of planets by means of thermal, (weakly) non-linear, or resonance-locking dissipative effects (Goodman & Dickson 1998; Barker & Ogilvie 2010; Weinberg et al. 2012; Ivanov et al. 2013; Essick & Weinberg 2016; Fuller 2017; Ma & Fuller 2021). For solar-type stars with radiative cores, the most likely dominant mechanisms for the dissipation of tidally excited gravity waves are non-linear dissipation, triggered by planetary companions with mass M_p ≳ 0.3 M_J (Essick & Weinberg 2016); or wave breaking, whose triggering requires the planetary mass to be larger than a critical value, which strongly varies as a function of the stellar age from ∼ 10² − 10³ M_J at t_⋆ ∼ 10⁸ yr to ∼ 10⁻¹ − 10⁻² M_J at t_⋆ ∼ 10¹⁰ yr for a 1 M_⊙ star (Barker 2020; Lazovik 2021).

Stars slightly more massive than the Sun (1.1 ≲ M_⋆/M_⊙ ≲ 1.6) are generally characterised by a rather complex structure on the MS, which is composed of an inner convective core surrounded by a radiative region topped by a shallow convective envelope. In this configuration, the tidally excited gravity waves that travel from the radiative-convective interface towards the stellar centre are reflected at the convective core interface, typically resulting in inefficient dissipation (Barker 2020). An alternative scenario has been proposed for the efficient dissipation of gravity waves propagating in such mixed cores, which consist of the conversion of gravity waves into magnetic waves (Lecoanet et al. 2017, 2022; Rui & Fuller 2023). In this context, tidally excited gravity waves travelling from the radiative-convective envelope boundary towards the stellar centre – upon encountering a sufficiently strong magnetic field generated by the dynamo at work in the convective core – would be fully converted into outwardly propagating magnetic waves; finally, they would be dissipated by radiative or ohmic diffusion (see Duguid et al. 2024, and references therein). By using stellar models with initial masses between 1.2 and 1.6 M_⊙, Duguid et al. (2024) showed that the wave conversion mechanism may operate for a significant fraction (of the order of gigayears) of stars’ MS lifetimes. This happens whenever the local radial magnetic field in proximity of the convective core is larger than a critical value at which the radial wave number of tidally excited gravity waves and magnetic waves match (Fuller et al. 2015; Lecoanet et al. 2017, 2022; Rui & Fuller 2023). Differently from the wave-breaking mechanism, wave conversion does not require a plantet’s minimum-mass threshold to be triggered. Another promising dissipation mechanism of tidal gravity waves in this type of star is tidal resonance-locking, in which a planet locks into resonance with a tidally excited stellar gravity mode. This process, similarly to the wave conversion one, can operate for planets of any mass in stars with convective cores (Ma & Fuller 2021).

It is worth stressing that all the tidal excitation and dissipation processes mentioned above strongly depend on the structural and rotational properties of the host star and on its age. The significant uncertainties derived in Sect. 2.1 on the mass, age, and initial chemical composition of GJ 504 – essentially due to a lack of asteroseismic constraints – make it difficult to understand which tidal excitation mechanism might have dominated its past evolution in the hypothetical presence of a close-by planetary companion. While the last two mechanisms mentioned above could have played a major role in this sense, given that GJ 504 could harbour an inner convective core according to our best-fit model, in this work we solely focused on the impact of equilibrium and dynamical tides dissipated within the stellar convective envelope, and we leave the investigation of the impact of tides dissipated in the radiative regions to a future work.

2.2.3. Radiative interaction

In the SPI code, it is possible to compute the erosion of the planetary atmosphere due to the impact of stellar X-ray and extreme ultraviolet (EUV) flux (e.g. Watson et al. 1981; Lammer et al. 2003; Jin et al. 2014). The evolution of the X-ray luminosity is computed by carrying out a recalibration of the Rx-Ro relationship of Johnstone et al. (2021) as in Pezzotti et al. (2021), where Rx = L_X/L_Bol is the ratio between the X-ray luminosity and the bolometric one, while Ro is the stellar Rossby number. The L_X and P_rot evolution are tightly linked by means of the dynamo-activity-rotation feedback loop (e.g. Parker 1955; Wilson 1966; Kraft 1967); therefore, it is important to simultaneously compute these two quantities. For the computation of the EUV flux, we referred to Johnstone et al. (2021).

Depending on the planetary and system properties, the mass-loss rate from the planet may be computed either by using an energy-limited formula (Watson et al. 1981; Erkaev et al. 2007), in which the factor accounting for the evaporation efficiency is estimated by following Salz et al. (2016) or Caldiroli et al. (2022) or by interpolating (eventually extrapolating) the grid of upper-atmosphere models of Kubyshkina et al. (2018), Kubyshkina & Fossati (2021). Regarding the computation of the planetary radius, depending on the initial mass, fraction of atmosphere, and composition, we used fitting formulae from Lopez & Fortney (2014), Chen & Rogers (2016) or mass-radius relations (Otegi et al. 2020; Bashi et al. 2017). For the range of planetary masses considered in this work, we used Bashi et al. (2017).

3. Analysis of the X-ray data from eROSITA

3.1. Method and tools

We looked for X-ray detections of GJ 504 in the five all-sky surveys of eROSITA (eRASS; Merloni et al. 2024) using the tool EROSE provided by the German consortium. First, we retrieved the sky map where the source is located during the five surveys. Then, we performed the source detection using the eSASSusers_240410 software release (Brunner et al. 2022) within the 0.2 − 5.0 keV energy band, which is usually adopted for M dwarfs (Magaudda et al. 2022). We compiled the list of the detected sources for each of the five surveys using the dedicated eSASS pipeline, ERMLDET, for which we used a threshold detection maximum likelihood of 6.0. Finally, we propagated the Gaia-DR3 coordinates of GJ 504 for their proper motion and at the epoch of the five eROSITA surveys and cross-matched within 10″ these coordinates with those in the result of our source detection. The results are shown in Table 1, where we present the eROSITA survey to which the observation refers (col. 1); the actual observational period of GJ 504 (col. 2); the separation between the proper-motion-corrected coordinates and the X-ray source position from eROSITA (col. 3); the number of source counts with the detection likelihood (cols. 4 and 5); and the source count rate (col. 6).

3.2. eROSITA spectral analysis

We created the spectrum (shown in Fig. 2), response matrix, and ancillary file for each eROSITA detection using the SRCTOOL pipeline with the AUTO option for the choice of the source and background region sizes. We grouped each spectrum with ten photons per channel for all eRASS, except for eRASS 1 & 5, for which we adopted 15 photons per channel. We fitted the four spectra with XSPEC version 12.13 (Arnaud 1996), adopting the XSPEC library for the stellar abundances compiled by Asplund et al. (2009). We used a one-temperature thermal APEC model with global coronal abundance (Z) fixed at 0.3 Z_⊙, as it is typically considered for stellar coronae (Maggio et al. 2007). We calculated the model-dependent X-ray fluxes with the dedicated XSPEC pipeline FLUX initially adopting the eROSITA energy band. For a consistent comparison with the X-ray luminosity from ROSAT (Voges et al. 1999), we also performed a calculation in the ROSAT energy band (0.1 − 2.4 keV), finding that the difference in emitted flux within the two bands is not significant, as previously shown by Magaudda et al. (2022). We used the model-dependent flux with a Gaia-DR3 distance of 17.58 pc to calculate the X-ray luminosity (L_x). The results of our spectral analysis are summarised in Table 2, where we show the eROSITA survey (col. 1), the coronal temperature, and emission measure with their 1σ uncertainties (cols. 2 and 3), the reduced chi-squared value (${\chi }_{\mathrm{red}}^2$) with the degrees of freedom (cols. 4 and 5), and the X-ray luminosity in the 0.1 − 2.4 keV energy band. We notice that the X-ray luminosity derived in this work with an average value of 10^29.42 erg/s is in line with that of Foster et al. (2022) (10^29.4 erg/s).

Fig. 2. eROSITA spectra of GJ 504 plotted together with best-fitting apec model and the residuals (both shown in green); see text in Sect. 3.2 and Table 2.

4. Compared evolutionary scenarios: Ω_surf and L_X

Once the best-fit model of GJ 504 is computed as described in Sect. 2.1, it is then coupled with our SPI code in order to compute the evolution of the host-star’s surface-rotation rate (Ω_surf) and X-ray luminosity (L_x) in two different evolutionary scenarios: on one hand, we simulate the evolution of the host star by assuming that there is no close-by, massive planet (star without close-by planet; see Sect. 4.1); on the other hand, we assume that a massive planet formed in a short orbit (a_in ≤ 0.1 AU) (or migrated close to its host star before the dissipation of the protoplanetary disc), and we studied the impact of its eventual inward migration driven by tidal dissipation and engulfment by the host star on the stellar properties (star with close-by planet; see Sect. 4.2).

With the rotational history of GJ 504 being unknown, we compute models representative of different rotators, from the super-slow to the fast one, using different values for the initial surface-rotation rate (Ω_in = 1, 2, 3.2, 5, 6, 7, 8, 9, 10, 18 Ω_⊙). Regarding the lifetime of the protoplanetary disc, we considered τ_dl = 2 Myr for Ω_in = 18 Ω_⊙ and τ_dl = 6 Myr for the other values. We recall that the choice of the Ω_in and τ_dl is determined from the distribution of surface-rotation rates observed for stars in star-forming regions and young open clusters at various ages (Eggenberger et al. 2019)⁴.

4.1. Star without close-by planet (SwoP)

In this evolutionary scenario, we assume that there is no close-by, massive planet orbiting GJ 504. Even if the detection of a sub-stellar companion through direct imaging was announced by Kuzuhara et al. (2013), with a mass of 1 ≲ M_c/M_J ≲ 25, this object is too far away from its host star (≈43.5 AU) to have any significant impact on the evolution of its surface-rotation rate and X-ray luminosity, at least according to the physics included in our SPI code. In this context, in Fig. 3 we collect the results obtained in terms of evolution of the surface-rotation rate, stellar Rossby number, and X-ray luminosity.

Fig. 3. Host star’s evolutionary tracks in the SwoP scenario. Top panel: Surface-rotation-rate evolution versus age for our optimal stellar model, with Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙. The dashed red line shows the evolution of the critical rotation velocity (Ωcrit), defined as the velocity at which the centrifugal acceleration at the equator equals the gravity. The magenta and black markers show GJ 504’s surface rotation rate, with age uncertainties derived from the global and local minimisation modelling, respectively. Middle panel: Evolution of stellar Rossby number (RO). RO for GJ 504 is indicated by the red, black, and blue circles corresponding to the lower, mean, and upper values of the largest age uncertainty, respectively. The grey and black stars represent the values obtained from τconv as in Wright et al. (2011) and Wright et al. (2018). Analogously, the magenta markers represent the same quantities for the smallest age uncertainty. Bottom panel: Evolution of X-ray luminosity for each of the considered rotators. The black and blue markers show the X-ray luminosities from Voges et al. (1999), Wright et al. (2011) (ROSAT) and in this work (eROSITA), respectively. Analogously, the magenta markers represent the same quantities for the smallest age uncertainty. In all panels, the starting point of the tracks corresponds to the dissipation of the protoplanetary disc (2 Myr for Ωin = 18 Ω⊙; 6 Myr for the other ones).

4.1.1. SwoP: Surface-rotation rate

In the top panel of Fig. 3, we show the evolution of the surface rotation rate normalised to the value of the Sun (Ω_⊙ = 2.9 × 10⁻⁶ s⁻¹) for each of the considered rotational histories as a function of the stellar age. As expected in the case of stars without close-by planets, the evolution of the surface-rotation rate is initially governed by the contraction of the stellar radius during the PMS phase until the star reaches the zero-age main sequence (ZAMS) at about 25 Myr. Subsequently, the braking of the stellar surface due to magnetised winds takes over, determining the overlap of the majority of the tracks at ≈500 Myr, except for the 1 Ω_⊙ model, which globally evolves at lower values. By comparing the theoretical tracks with the rotation rate derived for GJ 504 in Di Mauro et al. (2022) (P_rot = 3.4 days, Ω_{GJ 504} = 2.13 × 10⁻⁵ s⁻¹), we notice that there is only compatibility if we consider the youngest possible value for the stellar age (0.36 Gyr), and this is only for Ω_in ≥ 10 Ω_⊙; otherwise, Ω_{GJ 504} is larger than what was expected from our models.

We analysed whether considering the onset of weakened magnetic braking at Ro = Ro_cri could result in a better overlapping of the theoretical tracks and the observational data. Because of the impact of this mechanism indeed, the braking of the stellar surface would be stalled after a certain age, leading to larger surface-rotation rates at the age of GJ 504. The trend of the evolutionary tracks in this circumstance, computed only for Ω_in = 3.2, 5, 18 Ω_⊙, is shown by the dashed grey lines in the top panel of Fig. 3. It is possible to noticed that at the onset of weakened magnetic braking (∼2 Gyr), the star has already been significantly braked, and there is no compatibility with the observational value.

4.1.2. SwoP: Stellar Rossby number

In the middle panel of Fig. 3, the stellar Rossby number evolutionary tracks are presented. To make a comparison with the Rossby number obtained by using P_rot = 3.4 d, and given the significant uncertainty on GJ 504 age, we estimated three different R_O at the lower, mean, and upper limits of the error bar. Using Eq. (6), we derive R_Olow = 0.65, R_Om = 0.62, R_Oup = 0.43. Similarly to what we found for the surface rotation rate, an agreement with the evolutionary tracks was only obtained for the smallest value of the host star ages, namely for R_Olow and for Ω_in ≥ 10 Ω_⊙. It is worth noticing that we used a model-dependent approach to compute the Rossby number by means of Eq. (6) for τ_conv in order to be consistent with the physics included in the SPI code. However, it is possible to compare our estimation of τ_conv with the values derived by prescriptions based on purely observational properties. This is the case of the formulae in Wright et al. (2011, 2018) based on the (V − K_S) colour, which were calibrated on a sample of solar and late-type MS stars. In order to make a comparison with our model-dependent convective turnover timescale computation, we used Eq. (10) from Wright et al. (2011) and Eq. (5) in Wright et al. (2018), with (V − K_S) = 1.29 (Wright et al. 2011). We finally obtain τ_convW11 = 10.3228 d and τ_convW18 = 9.1727 d. In Fig. 4, we compare these two values with the evolutionary track. At ≈2.11 Gyr, there is a difference of about 8 d between the purely observational and theoretical estimations of τ_conv. Consequently, the corresponding Rossby numbers also show significantly diverse values in the middle panel of Fig. 3, where R_Om = 0.62, R_OW11 = 0.33 and R_OW18 = 0.37.

Fig. 4. Comparison between convective turnover timescale obtained from CLES best-fit model using Eq. (6), and the values computed as in Wright et al. (2011) (grey star) and Wright et al. (2018) (black star).

In the same panel, we also indicate the value of the critical Rossby number (R_Ocri, dotted black line, Metcalfe et al. 2024), above which the star may transition to the weakened magnetic braking regime (van Saders et al. 2016). According to our theoretical models, a star such as GJ 504 would cross R_Ocri at an age of ∼1 − 2 Gyr. As mentioned in the previous sections, at this point the star has been significantly braked by the magnetised winds to values well below Ω_GJ504. Thus, magnetic braking seems not to be responsible for the faster rotation of the star.

4.1.3. SwoP: X-ray luminosity

In the bottom panel of Fig. 3, the X-ray luminosity tracks are shown. Each track corresponds to a different rotational history. We compare the tracks with observational data from Voges et al. (1999), Wright et al. (2011) for ROSAT – with L_X = 10^29.47 erg/s (black dot) – and from our analysis of eROSITA, from which we derive an average value of L_X = 10^29.42 erg/s. Differently than what we observed for Ω and Ro, in this case an agreement between the theoretical tracks and the observed value is not observed, even for the lowest age value allowed for GJ 504. It is worth recalling that low and solar-like stars may experience significant variations in the emission of X-ray luminosity (over one order of magnitude), which are usually linked with the periodicity of their magnetic cycle.

4.2. Star with close-by planet (SwP)

In the star-with-close-by-planet scenario, we test the impact of the inward migration of a giant planet on the surface-rotation rate and X-ray luminosity of the host star, GJ 504. On the basis of analytical considerations, Fuhrmann & Chini (2015) proposed that the minimum mass needed for the planet to reproduce the observed stellar spin up is M_pl ≳ 2.7 M_J. In their work, the authors assumed that the angular momentum of the star is entirely related to the transfer of angular momentum from the planetary orbit, which occurred once the planet had filled its Roche lobe. In this work, we account for the intrinsic angular momentum of the host star together with that transferred by the planetary orbit. We thus studied the impact of equilibrium and dynamical tides dissipated within the host star’s convective envelope on the migration of the planet and its feedback on the stellar surface-rotation rate and X-ray luminosity evolution.

We first explored a range of initial orbital distances (0.02 ≤ a (AU) ≤ 0.1) and planetary masses (1 ≤ M_pl/M_J ≤ 10), and we selected stellar models representative of slow, moderate, and fast rotators – namely Ω_in = 3.2, 5, 18 Ω_⊙ (Eggenberger et al. 2019) – to find for which combination the computed stellar rotational period and X-ray luminosity are compatible with the observed values. In Fig. 5, we show the corresponding results. From these computations, we find that for the slow and moderate rotators it is possible to determine a region of compatibility in the a − M_pl plane, which is indicated by the black hatched area in the top and middle panels of Fig. 5. Contrarily, for the fast rotator case this compatibility is not retrieved. In this case, the planetary engulfment generally occurs at too early evolutionary stages (∼4 Myr), when the star is spinning up due to its contraction in the PMS; thus, the contribution to the stellar angular momentum is not only less important, but also quickly erased by subsequent magnetic braking. A further difference between the slow-moderate and fast rotators is that while in the first two cases the orbit of the planet shrinks for 0.03 ≲ a(AU) ≲ 0.06, in the last one the orbit expands for 0.025 ≲ a(AU) ≲ 0.1. This is due to the different geometry of the corotation radius (a_cor) and, consequently, a_min, which at the beginning of the evolution is smaller for larger values of Ω_in (see Eq. (9)). To clarify this point, in Fig. 6 we show an example of orbital evolution for a 1 M_J planet with a_in = 0.02, 0.025, 0.03, 0.035, 0.04, 0.06, 0.08, 0.1 AU, and Ω_in = 3.2, 18 Ω_⊙. In the slow rotator case, most of the tracks (a_in ≤ 0.06 AU) begin the evolution below a_min; thus, dynamical tides are not at work, and the orbit of the planet remains constant. Given that the initial surface-rotation rate of the star is quite low and the efficiency of tidal dissipation is tightly linked with this quantity, even after crossing a_min at about 10 − 20 Myr the orbit is not significantly affected by dynamical tides. Equilibrium tides finally deflect the orbit for a_in = 0.02, 0.025 AU at ∼2 Gyr, making the planet migrate towards its host star, since a_in < a_cor. It is worth noting that we followed the evolution of the planet until it reached its Roche limit⁵ (Zhang & Penev 2014). Following Benbakoura et al. (2019), at this point we assume that it is instantaneously depleted by tidal interactions and transfers its orbital angular momentum to the host star. As indicated in Metzger et al. (2012), when the ratio between the planetary and stellar density is in the 1 < ρ_pl/ρ_⋆ < 5 range and the planet overflows its Roche lobe, an unstable mass transfer takes place, tearing apart the planet within a timescale of several hours. For the range of planetary masses considered here, for the slow and moderate rotators, the ratio ρ_pl/ρ_⋆ at the engulfment (∼2 Gyr) is about 2–3. Thus, the assumption of instantaneous planetary disruption seems to be reasonable. Concerning the fast rotator case, as shown in Fig. 6, dynamical tides are much more efficiently dissipated, and planets starting the evolution below the corotation radius quickly migrate towards the host star, reaching the Roche limit. At this evolutionary stage, the ratio ρ_pl/ρ_⋆ > 5, and, according to Metzger et al. (2012), the planet would spiral and finally plunge into the stellar atmosphere. This event probably lasts more than the aforementioned instantaneous depletion. Even if we assume an instantaneous angular momentum transfer, to maximise its effect, given that the star is still spinning up at this stage it does not significantly impact the rotational evolution of its surface (at least according to the physics included in this work). Planets starting their evolution above the corotation radius (0.025 < a_in (AU) < 0.08) are instead pushed outwards.

Fig. 5. Different orbital evolutions for a planet with 1 ≤ Mplin/MJ ≤ 10 and 0.02 ≤ ain (AU) ≤ 0.1 computed for Ωin = 3.2 Ω⊙ (top panel), Ωin = 5 Ω⊙ (middle panel), and Ωin = 18 Ω⊙ (bottom panel). The salmon-coloured hatched area shows the parameter space for which, after the engulfment, it is possible to reproduce the rotational period and X-ray luminosity of GJ 504 when considering the largest age uncertainty. The lime-hatched area indicates a subset of the salmon-hatched area, for the smaller age uncertainty. The grey shaded area shows the region below the Roche limit.

Fig. 6. Orbital evolution for 1 MJ planet with 0.02 ≤ ain (AU) ≤ 0.1 and Ωin = 3.2, 18 Ω⊙ (top and bottom panels, respectively). The solid blue lines show the evolution of the orbital distance. The dotted red and dashed black lines indicate the evolution of amin and acor, respectively. The Roche limit is shown by the black dotted line. The cyan and magenta shaded areas represent the extension of the stellar core and envelope along the evolution of the star, respectively.

In the contextof the orbital evolution, we computed the mass loss from the planet due to the stellar XUV flux. Given the large planetary masses considered in these simulations, the mass loss due to the stellar high energy irradiation turns out to be negligible (e.g. Owen & Wu 2013).

After exploring the parameter space in M_pl, a_in, and Ω_in, in the following we focus on studying the orbital evolution of a planet with a mass of 3 M_J, which was arbitrarily chosen in the range of compatibility, and with a_in = 0.025 AU, for which we analyse a finer grid in Ω_in = 1, 2, 3.2, 5, 6, 7, 8, 9, 10, 18 Ω_⊙, as for the star-without-close-by-planet scenario. The goal is to investigate for which range of rotational histories the simulations reproduce the observed P_rot and L_X. The corresponding results are shown in Fig. 7⁶.

Fig. 7. Host star’s evolutionary tracks in the SwP scenario. Evolution of stellar surface-rotation rate (top panel), stellar Rossby number (middle panel), and X-ray luminosity (bottom panel) under the impact of planetary inward migration. For the stellar initial surface-rotation rate, we considered the following values: Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙. For the planetary mass, we considered Mpl = 3 MJ, and for the initial planetary orbital distance, we considered ain = 0.025 AU. The meaning of the other quantities is the same as in Fig. 3.

It is worth recalling that the results found in this section rely on the assumption that the planet formed (or migrated before the dispersal of the planetary disc) very closely to the host star, sitting on a coplanar-circular orbit. However, other mechanisms could have brought the planet to a very short orbit after its formation, with an initial planetary architecture characterised by a certain inclination and eccentricity. In particular, the presence of sub-stellar companion (at a projected distance of ∼43.5 AU) could have triggered the inner planet’s orbital shrinking by means of plant-planet scattering (e.g. Rasio & Ford 1996) and (or) Kozai-Lidov oscillations (Kozai 1962; Lidov 1962; Mazeh et al. 1997). The latter consists of a dynamical mechanism affecting the orbit of the inner planet with a periodically exchanging variation between its eccentricity and inclination. The excitation of the planetary eccentricity by means of these processes, together with tidal dissipation within the inner planet at every passage at the periastron, would potentially relocate the planet to a much shorter and circular orbit, which would eventually be short enough to induce tidal dissipation in the stellar host. In a more extreme scenario, for relatively large eccentricity values the planet could be directly engulfed at one of its periastron passages. Depending on the properties of the perturber, the initial orbital parameters of the system, and the structure of the inner planet, the timescales for these mechanisms to be triggered and result in an efficient planetary engulfment and spin-up of the host star might significantly vary. In the case of GJ 504, given the limited constrains available on the system properties, it is challenging to establish whether planet-planet scattering and (or) Lidov-Kozai oscillations could be triggered. Since stellar tides could still play a major role in a scenario when also considering these two processes, future simulations accounting for this mechanism could benefit from the results found in this work.

4.2.1. SwP: Surface-rotation rate

In the top panel of Fig. 7, we show the evolution of the stellar surface-rotation rate for each of the considered rotational histories as a function of the stellar age. As for the star-without-close-by-planet scenario, the PMS phase is dominated by the structural contraction until the star reaches the ZAMS, at which point the magnetic braking takes over. Since the planetary companion orbits very closely to the host star, because of the impact of equilibrium and dynamical tides, it migrates inwards. The majority of the orbital angular momentum is transferred to the star when the planet reaches the closest orbits to the Roche limit, and thus the star experiences a kick that results in a spike in the Ω vesus Age track. For Ω_in ≥ 3.2 Ω_⊙, the larger the value of Ω_in, the earlier the engulfment of the planet occurs. On the contrary, for Ω_in < 3.2 Ω_⊙ – since the orbital migration is only driven by equilibrium tides – the engulfment occurs at approximately the same age (∼2.5 Gyr). As mentioned above, for the fastest rotator considered here (18 Ω_⊙), the planet is engulfed during the PMS, and the angular momentum transferred from the planetary orbit does not contribute significantly to the global one of the star. This is shown from the yellow track in the top panel of Fig. 7, where at the moment of the engulfment (∼10 Myr) only a small increase of Ω is visible. For the other rotators, the planetary engulfment induces a stellar spin-up to about 11 − 12 Ω_⊙. The quasi-vertical increase of the surface-rotation rate is followed by a smoother decay, driven by the magnetic braking. Given that the timescale of the rotation-rate decay is larger than the spin-up one, if a planetary engulfment caused the acceleration of GJ 504, it is more probable that we are observing the star in the decelerating phase. In the less likely hypothesis that the planet has not been engulfed yet, but it is at the edge of its Roche limit, in transferring its angular momentum to the star, we can try to estimate the radial-velocity (RV) semi-amplitude signal induced on its host. Using Eq. (1) of Cumming et al. (1999), by considering the lowest initial mass required for the compatible companion to reproduce the observed P_rot and L_X for GJ 504, namely 2 M_J, the line-of-sight inclination from Bonnefoy et al. (2018) (i = 18.6 degrees), and an orbital distance of a_in = 0.01 AU (P_orb = 8 d), we obtain K = 52.18 m/s. D’Orazi et al. (2017) ruled out the presence of a massive companion close enough to tidally spin up the central star on the basis of RV monitoring performed at Lick Observatory (Fischer et al. 2014). The observed RV dispersion in this case is 25.7 m/s, which may be comparable to the activity level of the star. They thus excluded the presence of close-by companions more massive than 0.5 − 1 M_J. Nevertheless, if the planet is undergoing tidal destruction when crossing the Roche limit, the actual RV semi-amplitude could be much smaller than the ∼52 m/s estimated above, since it would only be caused by some planetary remnants at this stage. As reported in D’Orazi et al. (2017), this would require purposely designed observations.

In terms of compatibility between the observational value for the surface rotation rate and the evolutionary tracks, if we look at the decaying branches for each Ω_in, it is possible to narrow down the initial surface-rotation rate to values ≤10 Ω_⊙. It is worth noting that the overlap between the evolutionary tracks and the observational Ω_⋆ would favour ages between 360 Myr and 3 Gyr (approximately), while no overlap is obtained for older ages.

4.2.2. SwP: Stellar Rossby number

Once the star is spun up by the planetary engulfment, its Rossby number decreases. This is presented in the middle panel of Fig. 7, where for each of the considered rotators the Rossby number sharply decreases when the planet reaches its Roche limit. In this context, the tracks predict smaller values of R_O, which become compatible with both R_Om and R_Oup, even though the maximum age for compatibility is ∼3 Gyr, as for the surface-rotation rate. Better compatibility is also observed with the Rossby number computed by means of the Wright et al. (2011, 2018) formulae.

4.2.3. SwP: X-ray luminosity

The final comparison between the X-ray luminosity tracks and the observational data obtained from ROSAT (Voges et al. 1999; Wright et al. 2011) and eROSITA are shown in the bottom panel of Fig. 7. Very interestingly, with the engulfment of a massive enough planet (M_pl ≥ 2 M_J) and for Ω_in ≲ 7 Ω_⊙, an overlapping of the tracks and the observational data is retrieved. Conversely to the comparison with the surface-rotation rate, the interval for the possible values of Ω_in compatible with L_XGJ504 is narrower. Therefore, in order to simultaneously reproduce both the rotational period and X-ray luminosity of GJ 504, the initial surface-rotation rate of the star should be ≤7 Ω_⊙.

5. Conclusion

In this work, we studied the properties of GJ 504 and investigated whether with the engulfment of a planetary companion the observed surface rotation rate and X-ray luminosity can be reproduced. On the basis of the best-fit stellar model of GJ 504 derived by means of a two-step minimisation technique, and by using our SPI code, we computed tracks for the stellar surface-rotation rate and X-ray luminosity assuming two evolutionary scenarios: star without close-by planet, in which there is no close-by companion influencing the evolution of the star; and star with close-by planet, in which a massive planet orbits closely and migrates inwards due to tides dissipated within the host star. The rotational history of the star being unknown, we selected a broad range of initial surface-rotation rates, from super-slow (1 Ω_⊙) to fast-rotator (18 Ω_⊙). When comparing the tracks obtained in the star-without-close-by-planet evolutionary scenario, we find that there is only compatibility with the surface-rotation rate of GJ 504 for Ω_in ≥ 10 Ω_⊙, within the considered age range (0.36 − 3.86) Gyr. There is no compatibility with the observed X-ray luminosity, for which the tracks predict smaller values.

In the star-with-close-by-planet scenario, we initially explored the parameter space in planetary initial mass (1 ≤ M_pl(M_J)≤10), orbital distance (0.02 ≤ a_in(AU)≤0.1), and host-star surface-rotation rate (Ω_in = 3.2, 5, 18 Ω_⊙) to find the set of conditions needed for the planet to efficiently spin up the star after the engulfment. As a result, we find that for M_pl ≥ 2 M_J, 0.02 ≤ a_in(AU) ≤ 0.035 and host star Ω_in < 18 Ω_⊙, the engulfed planet is able to spin up the host star to values compatible with both the surface-rotation rate and X-ray luminosity observed for GJ 504.

Given the broad range of the considered host-star surface-rotation rates, within the aforementioned parameter space we arbitrarily focused on studying the orbital evolution of a 3 M_J planet, at a_in = 0.025 AU, and a finer grid in Ω_in = 1, 2, 3.2, 5, 6, 7, 8, 9, 10, 18 Ω_⊙, with the aim of breaking down the degeneracy on the GJ 504 initial surface-rotation rate. In this star-planet configuration, we find that the surface-rotation rate and X-ray luminosity of GJ 504 can be both reproduced for Ω_in ≤ 7 Ω_⊙. Moreover, for 3.2 ≤ Ω_in/Ω_⊙ ≤ 7, the larger the initial surface-rotation rate, the earlier the engulfment occurs. An interesting result is that the overlap between the evolutionary tracks and the observational Ω_⋆ seems to favour ages between 360 Myr and 3 Gyr for GJ 504, approximately, while no overlap is obtained for older ages. In Appendix B, we show the results obtained by doing a similar exercise, but this time we fixed Ω_in = 3.2 Ω_⊙ and let the initial orbital distance vary between 0.02 and 0.1 AU. Analogously to what was found above, the shorter the initial orbital distance, the earlier the engulfment occurs. Moreover, in this case the overlap between the tracks and the observational values also tends to favour ages for GJ 504 up to ∼2.5 Gyr. These results clearly show that despite the fact we managed to determine a region in M_pl-a_in-Ω_in within which the theoretical tracks are compatible with the observational constraints for GJ 504, the problem still remains strongly degenerate.

According to the two evolutionary scenarios analysed in this work, only the star-with-close-by-companion one is able to produce tracks compatible with both the surface-rotation rate and X-ray luminosity of GJ 504, which results to be a strong candidate for a star that underwent strong star-planet interactions. It is interesting to notice that the planetary inward migration and consequent destruction at the Roche limit should have occurred late enough in the MS, otherwise the braking at the stellar surface due to magnetised winds would easily have erased any trace of spin-up.

It is worth stressing that despite the fact that our work supports the engulfment scenario as a means to reconcile the age estimation discrepancies between activity indicators and classical parameters, thus favouring a more evolved evolutionary stage for GJ 504, the problem of deriving accurate stellar parameters remains. This will be crucial to establishing the nature of GJ 504 and characterising the impact of SPI. In this context, asteroseismology is the only tool able to address this need. Future observations with TESS in the announced ultra-short cadence mode in the extended campaign might open up the opportunity of detecting solar-like oscillations.

Acknowledgments

We thank the reviewer Dr. Yaroslav Lazovik for his thorough and helpful comments, that helped to improve the quality and clarity of the manuscript. CP thanks the Belgian Federal Science Policy Office (BELSPO) for the financial support in the framework of the PRODEX Program of the European Space Agency (ESA) under contract number 4000141194. GB is funded by the Fonds National de la Recherche Scientifique (FNRS). EM is supported by Deutsche Forschungsgemeinschaft under grant STE 1068/8-1. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS/NRTA software system developed by the German eROSITA consortium. VVG is an F.R.S.-FNRS Research Associate. SB acknowledges funding from the Dutch Research Council (NWO) with project number OCENW.M.22.215 of the research program ‘Open Competition Domain Science- M’.

References

Appendix A: Star with 1 M_J companion

The evolution of the host star surface-rotation rate and X-ray luminosity is significantly affected by a planetary companion at a_in = 0.02 AU, and with an initial mass 1 M_J. As we show in Fig. A.1, the star spins up when the planet migrates inwards and reaches its Roche limit, and simultaneously its X-ray emission increases. Nevertheless, a 1 M_J planet is not massive enough to transfer the required angular momentum from its orbit to the star to increase its X-ray luminosity at the point to be compatible with the observational data. The minimum mass required, according to our computations, is 2 M_J.

Fig. A.1. Evolution of the stellar surface-rotation rate (top panel), stellar Rossby number (middle panel) and X-ray luminosity (bottom panel), under the impact of planetary inward migration. For the stellar initial surface rotation rate we considered the following values Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙, for the planetary mass Mpl = 1 MJ, and initial planetary orbital distance ain = 0.02 AU. The meaning of the other quantities is the same as in Fig. 3.

Appendix B: Star with 3 M_J companion and Ω_in = 3.2 Ω_⊙

In Fig. 7 we show how the host star surface-rotation rate and X-ray luminosity are impacted by the inward migration of a planet with M_in = 3 M_J and a_in = 0.025 AU, while for the host-star a range of initial surface rotation rates is explored (Ω_in = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω_⊙). Similarly, in Fig. B.1 we show the results obtained by fixing Ω_in to 3.2 Ω_⊙ and varying the initial orbital distance in the range a_in = 0.02, 0.025, 0.03, 0.035, 0.04, 0.06, 0.08, 0.1 AU. In this context, the shorter the orbital distance and the earlier the engulfment of the planet occurs.

Fig. B.1. Evolution of the stellar surface-rotation rate (top panel), stellar Rossby number (middle panel) and X-ray luminosity (bottom panel), under the impact of planetary inward migration. In this case, we fix the initial mass of the planet (3 MJ), the initial surface rotation rate of the host star (Ωin = 3.2 Ω⊙), while we vary the initial orbital distance (ain = 0.02, 0.025, 0.03, 0.035, 0.04, 0.06, 0.08, 0.1 AU). The meaning of the other quantities is the same as in Fig. 3.

All Tables

All Figures

	Fig. 1. Hertzsprung-Russell diagram of GJ 504 as obtained from best-fit stellar model (Teff = 6202 K, log(g) = 4.29, R = 1.35 R⊙, L/L⊙ = 2.43) computed with CLES.
In the text

	Fig. 2. eROSITA spectra of GJ 504 plotted together with best-fitting apec model and the residuals (both shown in green); see text in Sect. 3.2 and Table 2.
In the text

Fig. 3. Host star’s evolutionary tracks in the SwoP scenario. Top panel: Surface-rotation-rate evolution versus age for our optimal stellar model, with Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙. The dashed red line shows the evolution of the critical rotation velocity (Ωcrit), defined as the velocity at which the centrifugal acceleration at the equator equals the gravity. The magenta and black markers show GJ 504’s surface rotation rate, with age uncertainties derived from the global and local minimisation modelling, respectively. Middle panel: Evolution of stellar Rossby number (RO). RO for GJ 504 is indicated by the red, black, and blue circles corresponding to the lower, mean, and upper values of the largest age uncertainty, respectively. The grey and black stars represent the values obtained from τconv as in Wright et al. (2011) and Wright et al. (2018). Analogously, the magenta markers represent the same quantities for the smallest age uncertainty. Bottom panel: Evolution of X-ray luminosity for each of the considered rotators. The black and blue markers show the X-ray luminosities from Voges et al. (1999), Wright et al. (2011) (ROSAT) and in this work (eROSITA), respectively. Analogously, the magenta markers represent the same quantities for the smallest age uncertainty. In all panels, the starting point of the tracks corresponds to the dissipation of the protoplanetary disc (2 Myr for Ωin = 18 Ω⊙; 6 Myr for the other ones).

In the text

	Fig. 4. Comparison between convective turnover timescale obtained from CLES best-fit model using Eq. (6), and the values computed as in Wright et al. (2011) (grey star) and Wright et al. (2018) (black star).
In the text

Fig. 5. Different orbital evolutions for a planet with 1 ≤ Mplin/MJ ≤ 10 and 0.02 ≤ ain (AU) ≤ 0.1 computed for Ωin = 3.2 Ω⊙ (top panel), Ωin = 5 Ω⊙ (middle panel), and Ωin = 18 Ω⊙ (bottom panel). The salmon-coloured hatched area shows the parameter space for which, after the engulfment, it is possible to reproduce the rotational period and X-ray luminosity of GJ 504 when considering the largest age uncertainty. The lime-hatched area indicates a subset of the salmon-hatched area, for the smaller age uncertainty. The grey shaded area shows the region below the Roche limit.

In the text

Fig. 6. Orbital evolution for 1 MJ planet with 0.02 ≤ ain (AU) ≤ 0.1 and Ωin = 3.2, 18 Ω⊙ (top and bottom panels, respectively). The solid blue lines show the evolution of the orbital distance. The dotted red and dashed black lines indicate the evolution of amin and acor, respectively. The Roche limit is shown by the black dotted line. The cyan and magenta shaded areas represent the extension of the stellar core and envelope along the evolution of the star, respectively.

In the text

Fig. 7. Host star’s evolutionary tracks in the SwP scenario. Evolution of stellar surface-rotation rate (top panel), stellar Rossby number (middle panel), and X-ray luminosity (bottom panel) under the impact of planetary inward migration. For the stellar initial surface-rotation rate, we considered the following values: Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙. For the planetary mass, we considered Mpl = 3 MJ, and for the initial planetary orbital distance, we considered ain = 0.025 AU. The meaning of the other quantities is the same as in Fig. 3.

In the text

Fig. A.1. Evolution of the stellar surface-rotation rate (top panel), stellar Rossby number (middle panel) and X-ray luminosity (bottom panel), under the impact of planetary inward migration. For the stellar initial surface rotation rate we considered the following values Ωin = 1, 2, 3.2, 4, 5, 6, 7, 8, 9, 10, 18 Ω⊙, for the planetary mass Mpl = 1 MJ, and initial planetary orbital distance ain = 0.02 AU. The meaning of the other quantities is the same as in Fig. 3.

In the text

Fig. B.1. Evolution of the stellar surface-rotation rate (top panel), stellar Rossby number (middle panel) and X-ray luminosity (bottom panel), under the impact of planetary inward migration. In this case, we fix the initial mass of the planet (3 MJ), the initial surface rotation rate of the host star (Ωin = 3.2 Ω⊙), while we vary the initial orbital distance (ain = 0.02, 0.025, 0.03, 0.035, 0.04, 0.06, 0.08, 0.1 AU). The meaning of the other quantities is the same as in Fig. 3.

In the text