# Unstable accretion in TW Hya: 3D simulations and comparisons with observations

M. M. Romanova, C. C. Espaillat, J. Wendeborn, J. -F. Donati, P. P. Petrov, R. V. E. Lovelace

► To cite this version:

M. M. Romanova, C. C. Espaillat, J. Wendeborn, J. -F. Donati, P. P. Petrov, et al. Unstable accretion in TW Hya: 3D simulations and comparisons with observations. Monthly Notices of the Royal Astronomical Society, 2025, 538, pp.480-502. 10.1093/mnras/staf148. insu-05100393

HAL Id: insu-05100393

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

Submitted on 6 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Unstable accretion in TW Hya: 3D simulations and comparisons with observations

M. M. Romanova,<fnref n="1" />,<fnref n="2" />* C. C. Espaillat,<fnref n="3" />,<fnref n="4" /> J. Wendeborn,<fnref n="3" />,<fnref n="4" /> J.-F. Donati,<fnref n="5" /> P. P. Petrov<fnref n="6" />
and R. V. E. Lovelace<fnref n="1" />,<fnref n="2" />

Accepted 2025 January 20. Received 2025 January 20; in original form 2024 November 25

## ABSTRACT

We investigate the origin of photometric variability in the classical T Tauri star TW Hya by comparing light curves obtained by *Transiting Exoplanet Survey Satellite (TESS)* and ground-based telescopes with light curves created using three-dimensional (3D) magnetohydrodynamic (MHD) simulations. TW Hya is modelled as a rotating star with a dipole magnetic moment, which is slightly tilted about the rotational axis. We observed that for various model parameters, matter accretes in the unstable regime and produces multiple hotspots on the star’s surface, which leads to stochastic-looking light curves similar to the observed ones. Wavelet and Fourier spectra of observed and modelled light curves show multiple quasi-periodic oscillations (QPOs) with quasi-periods from less than 0.1 to 9 d. Models show that variation in the strength and tilt of the dipole magnetosphere leads to different periodograms, where the period of the star may dominate or be hidden. The amplitude of QPOs associated with the stellar period can be smaller than that of other QPOs if the tilt of the dipole magnetosphere is small and when the unstable regime is stronger. In models with small magnetospheres, the short-period QPOs associated with rotation of the inner disc dominate and can be mistaken for a stellar period. We show that longer period (5–9 d) QPOs can be caused by waves forming beyond the corotation radius.

**Key words:** accretion discs – hydrodynamics – planet–disc interactions – protoplanetary discs.

## 1 INTRODUCTION

Classical T Tauri stars (CTTSs) show photometric variability on different time-scales, from seconds to decades (e.g. Herbst et al. 1994, 2002; Hartmann et al. 2016; Fischer et al. 2023). Observations of multiple CTTSs in clusters revealed that only a few CTTSs show periodic light curves. Most show quasi-periodic or stochastic-looking light curves (e.g. Alencar et al. 2010; Cody et al. 2014, 2022; Cody & Hillenbrand 2018; Robinson et al. 2022). Many light curves show bursts, which may indicate stochastic accretion (e.g. Stauffer et al. 2014, 2016).

To better understand the mass accretion process in CTTSs, a multi-epoch, multiwavelength spectroscopic and photometric monitoring campaign of four CTTSs (TW Hya, RU Lup, BP Tau, and GM Aur) was carried out in 2020/2021 (Epoch 1) and 2022/2023 (Epoch 2) as part of the *UV Legacy Library of Young Stars as Essential Standards (ULLYSES)* HST Director’s Discretionary Time Program (Roman-Duval et al. 2020), and the *Outflows and Disks Around Young Stars: Synergies for the Exploration of ULLYSES Spectra (ODYSSEUS)* programme (Espaillat et al. 2022). Light curves

were obtained using several ground-based telescopes at multiple wavelengths (see table 1 in Wendeborn et al. 2024b). All stars show significant variability in their light curves at different time-scales; some time-scales are associated with the period of the star (like in GM Aur), while in RU Lupi, the stochastic component dominates (Wendeborn et al. 2024a, b). Those authors conclude that there is a strong connection between mass accretion rate and photometric variability. Therefore, it is important to understand which type of accretion processes can produce such light curves.

CTTSs typically have a strong, ~kG-scale, magnetic field (e.g. Johns-Krull 2007; Donati & Landstreet 2009), which stops the accretion disc and the matter is expected to accrete on to the star in two funnel streams forming two ordered spots near magnetic poles (e.g. Königl 1991; Bouvier et al. 2007a; Hartmann et al. 2016). Early global three-dimensional (3D) magnetohydrodynamic (MHD) simulations of accretion on to a star with a tilted dipole magnetosphere confirmed the theory and have shown the magnetospheric accretion in two funnel streams (Romanova et al. 2003, 2004) which form two banana-shaped spots at the star’s surface, and the light curve is almost sinusoidal. Some observations suggest that accretion occurs through multiple streams (e.g. Ingleby et al. 2013; Johnstone et al. 2014; Robinson & Espaillat 2019).

<footnote n="1">¹Department of Astronomy, Cornell University, Ithaca, NY 14853-6801, USA</footnote>

<footnote n="2">²Carl Sagan Institute, Cornell University, Ithaca, NY 14853-6801, USA</footnote>

<footnote n="3">³Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="4">⁴Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="5">⁵CNRS/Université de Toulouse IRAP – Observatoire Midi-Pyrénées, 14 Av E Belin, F-31400 Toulouse Cedex, France</footnote>

<footnote n="6">⁶Crimean Astrophysical Observatory, Bakhchisaray, p. Nauchny, 298409, Republic of Crimea</footnote>

<footnote>*E-mail: romanova@astro.cornell.edu</footnote>
<table><caption>**Table 1.** Observational parameters of TW Hya. a1: Venuti et al. (2019), a2: Manara et al. (2014), a3: Baraffe et al. (2015), b1: (Hughes et al. 2007), b2: Robinson & Espaillat (2019), b3: Rhee et al. (2007), b4: Baraffe et al. (2015), Donati et al. (2024), b5: Venuti et al. (2019), b6: GRAVITY Collaboration (2020), c1: Setiawan et al. (2008), c2: Huélamo et al. (2008), c3: Donati et al. (2024), d1: Donati et al. (2011), d2: Donati et al. (2024), e1: Siwak et al. (2018), e2: Qi et al. (2004), e3: Teague et al. (2019).</caption><thead><tr><th>Parameter of TW Hya</th><th>Value</th></tr></thead><tbody><tr><td>Mass</td><td>0.6 M<sub>☉</sub><sup>a1</sup>0.79 M<sub>☉</sub><sup>a2</sup>0.8 M<sub>☉</sub><sup>a3</sup></td></tr><tr><td>Radius</td><td>R<sub>*</sub> = 0.85 ± 0.25<sup>b1</sup> 0.93 R<sub>☉</sub><sup>b2</sup>, R<sub>*</sub> = 1.1 R<sub>☉</sub><sup>b3</sup>, 1.16 ± 0.13 R<sub>☉</sub><sup>b4</sup>, 1.22 R<sub>☉</sub><sup>b5</sup>, 1.29 R<sub>☉</sub><sup>b6</sup></td></tr><tr><td>Period</td><td>3.56 d <sup>c1,c2</sup>3.606 ± 0.015 d <sup>c3</sup></td></tr><tr><td>Mag. field (dipole)</td><td>a few hundred Gauss <sup>d1</sup></td></tr><tr><td>Mag. field (octupole, almost aligned)</td><td>2.5–2.8 kG</td></tr><tr><td>Mag. field (dipole)</td><td>990–1190 G <sup>d2</sup></td></tr><tr><td>Mag. field (small-scale fields)</td><td>3–4 kG</td></tr><tr><td>Tilt of the dipole magnetosphere θ</td><td>(17–23)° ± (6–10)°</td></tr><tr><td>Inclination angle i</td><td>5–15°<sup>e1</sup>, 7°<sup>e2</sup> 5.6° (for the disc) <sup>e3</sup></td></tr></tbody></table>

Global 3D MHD simulations performed at a broader range of parameters and finer grid resolution have shown that in models with a small tilt of the dipole magnetosphere, the matter often accretes in the unstable regime, where it penetrates through the magnetosphere due to the magnetic Rayleigh–Taylor (interchange) instability (Romanova & Lovelace 2006). Earlier, this instability was proposed by Arons & Lea (1976) for mixing matter and the field in the external parts of the magnetosphere. However, the penetration of filaments (or tongues) into the deep layers of the magnetosphere was not expected. This unstable regime has been studied in detail by Romanova, Kulkarni, & Lovelace (2008) and Romanova & Kulkarni (2009) in 3D MHD simulations. Simulations show that matter accretes in several unstable tongues, which form at the inner edge of the disc, creating frequent and randomly located hotspots on the stellar surface. Fourier and wavelet analysis show that in a strongly unstable regime, quasi-periodic oscillation (QPO) associated with rotating tongues may dominate over the period of the star producing false ‘periods’.

In models with smaller magnetospheres and slower rotating stars, an ‘ordered unstable regime’ has been observed where matter accretes in one or two unstable tongues that rotate with the angular velocity of the inner disc. The period of their rotation may be much smaller than the period of the star (e.g. Romanova & Kulkarni 2009). Blinova et al. (2016) studied the unstable regime at multiple parameters and derived a boundary between stable and unstable regimes. An unstable accretion regime has also been observed in more recent 3D MHD simulations of other groups (e.g. Takasao et al. 2022; Parfrey & Tchekhovskoy 2024; Zhu et al. 2024).

Several groups used the paradigm of unstable regime to explain the light curves of CTTSs. For example Siwak et al. (2011) analysed the light curves of TW Hya obtained with the *Microvariability and Oscillations of STars (MOST)* satellite and concluded that irregular light variations could be caused by the stochastic variability of hotspots induced by unstable accretion. Takami et al. (2016) analysed spectra of the active star RW Aur A and suggested that it may switch between stable and unstable regimes during its evolution when comparing the complex variability of different spectral lines to that simulated by Kurosawa & Romanova (2013).

Overall, the unstable regime is frequently used to explain the quasi-periodic or stochastic light curves of CTTS. However, the observed light curves have never been compared with those obtained by dedicated numerical models of stars accreting in the unstable regime. Here, we choose one of the CTTS stars (TW Hya) targeted by the ODYSSEUS campaign and develop a 3D MHD model of a

star with the parameters of TW Hya. We analyse and compare the observed and modelled light curves and QPOs.

In Section 2, we present observational data of TW Hya and describe our numerical model in Section 3. In Section 4, we show the results of simulations and compare them with observations. In Section 5, we discuss different points and conclude in Section 6.

## 2 OBSERVATIONS OF TW HYA

TW Hya is a bright, nearby CTTS of the spectral type K6/K7 (when measured at blue/optical wavelengths; Manara et al. 2014). It retains a large massive disc at the age of 8–10 Myr (Sokal et al. 2018) and also contains a number of gaps and rings (e.g. Calvet et al. 2002; Andrews et al. 2016). The disc and the magnetosphere are both seen nearly pole-on (e.g. Qi et al. 2004; Donati et al. 2011). Table 1 summarizes the observational parameters of TW Hya.

### 2.1 Searching for a period of TW Hya

Many attempts have been made to determine the rotation period of TW Hya using photometric data obtained with ground-based telescopes. Rucinski (1988) did not find a significant period from observations in 1986, but reported a possible 2-d period from earlier data. Herbst & Koret (1988) found a 1.28-d period. Mekkaden (1998) obtained a 2.2-d period and pointed out the presence of sudden brightness, which could be due to short-lived hotspots that can occur at any photometric phase and could mask the periodicities. Kastner et al. (1999) examined the optical V-band photometry from *Hipparcos*, and showed that neither of the above periods produced a regular folded light curve. They suggested that quasi-random flaring, instead of a hotspot, could be causing the appearance and disappearance of periods in different data sets. Alencar & Batalha (2002) analysed spectra of TW Hya and performed periodogram analysis of the veiling and veiling-corrected line intensity variations using the Scargle (1982) periodogram estimator. They obtained three prominent periodicities at 1.4±0.1, 2.85±0.25, and 3.75±0.45 d. A similar study by Batalha et al. (2002) pointed to 4.4±0.4 d as the stellar rotation period. Lawson & Crause (2005) found a period of 2.80 d. So, ground-based observations suggested a range of periods for TW Hya.

Setiawan et al. (2008) detected periodic (3.56 d) radial velocity variability in their spectroscopic observations, interpreted as the rotation of a 10 Jupiter mass planet. However, Huélamo et al. (2008) has shown that this period can be explained by the rotation of a cool stellar spot, and they measured a stellar rotation period of 3.57 d.
The search for QPOs associated with magnetospheric accretion and hotspots became even more intense and interesting with data from space telescopes, which provide frequent and continuous observations. Rucinski et al. (2008) analysed photometric variability in TW Hya using data from the MOST space telescope on timescales from a fraction of a day to 7.5 weeks. A 3.7-d period was observed in the continuous 11-d observations performed with 0.07 d time resolution in 2007. However, this periodicity was absent in the second, four times longer MOST run in 2008. Instead, a spectrum of quasi-periods within the 2–9 d range was observed.

Siwak et al. (2011) studied the variability of TW Hya using the light curves obtained by the MOST satellite and the All Sky Automated Survey (ASAS) project over 40 d in 2009 with a temporal resolution of 0.2 d. A wavelet analysis of the combined MOST–ASAS data provided a rich picture of QPOs with periods of 1.3-10 d, similar to those discovered in the 2008 data. The authors concluded that the observed shortest oscillation period may indicate a stellar rotation period of 1.3 or 2.6 d, synchronized with the disc at 4.5 or 7.1$R_*$, respectively.

Siwak et al. (2014) presented an analysis of the 2011 photometric observations of TW Hya by the MOST satellite. The light variations were dominated by a strong, quasi-periodic 4.18-d oscillation with superimposed chaotic-looking flares. They conclude that the QPOs were probably produced by stellar rotation, with one large hotspot created by a stable accretion funnel, while the flaring component could be produced by small hotspots created at moderate latitudes by unstable accretion tongues.

Siwak et al. (2018) report on photometric variability observed by MOST during 2014, 2015, and 2017. In 2014 and 2017, the light curves showed stable 3.75- and 3.69-d QPOs, respectively. Both values appear closely related to the stellar rotation period, as they might be created by changing the visibility of a hotspot formed near the magnetic pole. These major light variations were superimposed on a chaotic, flaring-type activity caused by hotspots resulting from unstable accretion – a situation reminiscent of that in 2011 when TW Hya showed signs of a moderately stable accretion state. In 2015, only drifting QPOs were observed, similar to those present in 2008–2009 data and typical for magnetized stars accreting in a strongly unstable regime.

Observations of TW Hya with the Transiting Exoplanet Survey Satellite (TESS) were performed in Sectors 9, 36, and 63. Scargle–Lomb periodograms of TESS light curves show QPOs with a period close to the star’s period. However, QPOs with other periods are also observed and often have higher amplitudes (Wendeborn et al. 2024b). The ODYSSEUS programme’s ground-based telescopes observed TW Hya in 2021 (Epoch 1) and 2022 (Epoch 2). Scargle–Lomb periodograms show a QPO with a stellar period in Epoch 2. However, no QPO with a stellar period was observed in Epoch 1 (Wendeborn et al. 2024b). Therefore, observations of TW Hya show the star’s period during some observing runs but not during others.

## 2.2 Unstable regime. Fastness parameter

The unstable regime has been studied in multiple 3D MHD simulations by Blinova et al. (2016). They conclude that the ‘strength’ of instability and the boundary between stable and unstable regimes depends on the fastness parameter $\omega_s = \Omega_s/\Omega_K = (r_m/r_{cor})^{3/2}$, where $r_m$ is the magnetospheric (or truncation radius), where matter is stopped by the magnetosphere, and $r_{cor}$ is the corotation radius, where the angular velocity of the star matches the Keplerian angular velocity of the disc; $\Omega_K$ is the Keplerian angular velocity of the disc at $r = r_m$. At a small tilt of the dipole $\theta = 5^\circ$, they found a

boundary at $\omega_s \approx 0.6$ which corresponds to $r_m/r_{cor} \approx 0.71$, while at the large tilt $\theta = 20^\circ$, the boundary is at $\omega_s \approx 0.54$ which corresponds to $r_m/r_{cor} \approx 0.66$. They also observed that in smaller magnetospheres $r_m \lesssim 4.2R_*$, matter accretes in one or two ordered tongues if $\omega_s \approx 0.45$ ($r_m/r_{cor} \lesssim 0.59$). They rotate with the period of the inner disc and may provide a false period of the star, which can be much shorter than the actual period.

The corotation radius is determined by the period of the star and its mass. The period of TW Hya is known as $P_* \approx 3.56$ d (Huélamo et al. 2008; Setiawan et al. 2008). We present the corotation radius in the following form:

$$r_{\text{cor}} = \left[ \frac{GM_* P_*^2}{(2\pi)^2} \right]^{1/3} \approx 9.10 R_{\odot} \left( \frac{M_*}{0.8 M_{\odot}} \right)^{1/3} \left( \frac{P_*}{3.56 \text{ d}} \right)^{2/3}, \quad (1)$$

where $M_*$ is the mass of the star.

There are different estimates of TW Hya radius: $R_* = 0.85 \pm 0.25 R_\odot$ (Hughes et al. 2007); $R_* = 0.93 R_\odot$ (Robinson & Espaillat 2019); $R_* = 1.1 R_\odot$ (Rhee et al. 2007); $R_* = 1.16 \pm 0.13 R_\odot$ (Donati et al. 2024); $R_* = 1.29 \pm 0.19 R_\odot$ (GRAVITY Collaboration 2020). We adopt a radius of $R_* = 0.93 R_\odot$, mass $M_* = 0.79 M_\odot$, and stellar period $P = 3.5$ d and obtain a corotation radius $r_{cor} \approx 9.6 R_*$, which we use in most of our simulations.

There are different theoretical approaches for finding the magnetospheric radius. In one of them, it is suggested that the matter of the disc stops when magnetic pressure balances the total matter pressure: $B_d^2/8\pi = p + \rho v^2$, where $B_d$ is the magnetic field at the disc–magnetosphere boundary, $p-$ and $\rho v^2-$ are the thermal and ram pressure of matter in the disc. In accretion disc $p \ll \rho v^2$, $v \approx v_\phi$, so that $B^2/8\pi \approx \rho v_\phi^2$, where $v_\phi = (GM_*/r)^2$ is the Keplerian velocity of the disc. Substituting the magnetic moment $\mu_* \equiv B_* R_*^3$, we obtain a widely used formula (e.g. Lamb et al. 1973; Elsner & Lamb 1977):

$$r_m = kr_A = k \frac{\mu_*^{4/7}}{(2GM_*)^{1/7} M^{2/7}} = k \times 4.09 R_\odot \left( \frac{M_*}{0.8 M_\odot} \right)^{-1/7} \times \left( \frac{\dot{M}}{2.51 \times 10^{-9} M_\odot \text{yr}^{-1}} \right)^{-2/7} \left( \frac{B_{*\text{eq}}}{500G} \right)^{4/7} \left( \frac{R_*}{R_\odot} \right)^{12/7}, \quad (2)$$

where coefficient $k \sim 1$. Bessolaz et al. (2008) derived $k \approx 0.77 - 1$ by comparing the theoretical value with the results of 2D (axisymmetric) simulations. Kulkarni & Romanova (2013) derived $k = 0.55 - 0.72$ from similar comparisons but in 3D simulations. Blinova et al. (2016) derived $k \approx 0.5 - 0.9$ in stable regime and $k \approx 0.8 - 0.93$ in the unstable regime from multiple 3D simulations (see details in section A of their paper).<fnref n="1" />

In our research, we derive the magnetospheric radius from numerical models but use theoretical formulae for comparisons. The magnetospheric radius mainly depends on the accretion rate and magnetic moment (or, magnetic field and radius) of the star. Below, we briefly discuss the results of observational measurements of these values of TW Hya.

## 2.3 The magnetic field of TW Hya

Donati et al. (2011) used the spectropolarimetry method to measure the magnetic field of TW Hya from optical spectra secured with ESPaDOnS at the Canada–France–Hawaii Telescope (CFHT; Donati 2003). They concluded that the field is predominantly octupolar of

<footnote n="1">¹The formulae for $r_m$ is also slightly different from equation (2). The coefficient $k$ helps to adjust parameters to standard dependencies shown in equation (2).</footnote>
2.5–2.8 kG with a smaller dipole component of a few hundred Gauss. More recently, Donati et al. (2024) repeated observations of TW Hya in 2019, 2020, 2021, and 2022 using the Zeeman–Doppler Imaging<fnref n="2" /> and obtain that the large-scale field mainly consists of a 1 kG dipole tilted at about 20° to the rotation axis, whereas the small-scale field reaches strengths of up to 3–4 kG. They show that the strength of the dipole component varies from 990 to 1190 G, the tilt of the dipole moment varies from 17° up to 23°, and the phase varies too (see table 3 of Donati et al. 2011).

## 2.4 Accretion rate

The accretion rate of TW Hya obtained from continuum and line analyses was determined as $\dot{M} \approx 4 \times 10^{-10} - 6 \times 10^{-9} M_{\odot}yr^{-1}$ (e.g. Muzerolle et al. 2000; Alencar & Batalha 2002; Herczeg et al. 2002; Donati et al. 2011; Espaillat et al. 2019; Robinson & Espaillat 2019; Wendeborn et al. 2024a, c). An accretion rate $\dot{M} \approx 2.51 \times 10^{-9} M_{\odot}yr^{-1}$ has been derived from optical spectra taken mostly with the ESPaDOnS instrument (Herczeg et al. 2023). Those authors noted that this accretion rate might be underestimated by a factor of up to 1.5 because of uncertainty in the bolometric correction and another factor of 1.7 because of excluding the fraction of accretion energy that escapes in lines, especially $Ly\alpha$.<fnref n="3" /> If these factors are taken into account, then $\dot{M} \approx 6.40 \times 10^{-9} M_{\odot}yr^{-1}$.

Kastner et al. (2002) analysed the high-resolution X-ray spectrum of TW Hya and concluded that if accretion powers the X-ray emission, then the accretion rate should be ~ $10^{-8} M_{\odot}yr^{-1}$. Nayak et al. (2024) observed TW Hya with the Ultra-Violet Imaging Telescope (UVIT). Based on C iv line luminosity, they estimated an accretion luminosity ($0.1 L_{\odot}$) and a the mass accretion rate of $2.2 \times 10^{-8} M_{\odot}yr^{-1}$. Due to the wide range of accretion rates measured for TW Hya, we consider different possible accretion rates in this paper.

## 2.5 TESS observations of TW Hya

TW Hya was observed by TESS (Ricker et al. 2015) in 2019, 2021, and 2023 (Sectors 9, 36, 63, respectively). In Sector 9, observations were performed during $T = 25.20$ d with a temporal resolution of $\Delta T = 0.0195$ d (0.468 h = 28.08 min) in most cases. In Sector 36, $T = 25.08$ d and $\Delta T \approx 0.0078$ d (0.1872 h = 11.23 min). In Sector 63 $T = 26.53$ and $\Delta T = 0.0039$ (0.0936 h = 5.61 min).

The light curves show variability on different time-scales (see Fig. 1). The Lomb–Scargle periodograms are shown below the light curves for each sector. Red lines show periodograms for non-flagged data. In Sector 9, there are two prominent peaks with periods of 3.5–3.6 d (associated with the star’s rotation) and a higher peak at 4.7 d. There is also a smaller amplitude peak at 2.6 d. In Sector 36, one peak is at 3.8 d, and another peak is at 9 d. In Sector 63, there is one broad peak at 4.5–5.5 d and a peak of similar amplitude at 8.7 d. There are also several smaller amplitude peaks at shorter periods. Periodograms for the light curves, including flagged data (grey colour curves), look somewhat different and often show different peaks.

We also used the Morlet-6 wavelet to analyse TESS light curves. First, we carefully checked light curves for stellar flares using earlier developed approaches (e.g. Siwak et al. 2010, 2018; Günther et al.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top two panels: Top: TESS light curve in Sector 9 for TW Hya. Grey points have been flagged by either the default TESS pipeline or by tglc (typically due to high background flux), while red points are data that have not been flagged. Bottom: Lomb–Scargle periodograms for the associated light curves. The red curve is calculated using the non-flagged data, while the grey curve is calculated using all data, including flagged points. Shaded regions are $1\sigma$ uncertainties calculated using a bootstrap approach. Middle two panels: same but for Sector 36. Bottom two panels: same but for Sector 63.</figcaption></figure>

2020). We did not find any notable flares.<fnref n="4" /> Next, we interpolated the light curves to a strictly uniform grid (the data in TESS light curves are recorded at slightly non-equal time intervals). The new grid is much denser than the original one: it has temporal resolutions of 1.87, 0.625, and 0.216 min in Sectors 9, 36, and 63, respectively. These

<footnote n="2">²They observed TW Hya in the near-infrared with the SPIRou (Donati et al. 2020) high-resolution spectropolarimeter and velocimeter at the 3.6-m CFHT.</footnote>

<footnote n="3">³According to Arulanantham et al. (2023), the major part of the accretion energy escapes in $L_{y\alpha}$ radiation, and this factor may be significant (see also France et al. 2014).</footnote>

<footnote n="4">⁴In earlier studies of TW Hya, stellar flares were found in some sets of MOST observations, and not in others (Siwak et al. 2018).</footnote>
grids are 15–26 times denser than the original *TESS* grid, which
helped to resolve large and small bursts seen in the original light
curve. Next, we calculated the Morlet-6 wavelet using this uniform
grid. The top panels of Fig. 2 show the *TESS* light curves, which
include flagged data. Linear interpolation was used in parts of the
curve with no data points. The second panel from the top (for each
sector) shows the wavelet of the light curve. Wavelets typically show
quasi-periods of *P* ≈ 3.5–4.8 d which last *T* ≈ 5–15 d. This
quasi-period is close to the stellar period but varies in time. There
are also oscillations with periods of *P* ≈ 1–2 d. They usually last
for *T* ≈ 3–5 d. We stretched the wavelet in the vertical direction
(*P* < 2 d) and resolved oscillations with even lower periods, *P* ≈
0.2–0.9 d which last *T* ≈ 1–3 d (see third panels from the top
in Fig. 2). We observe that oscillations with shorter periods last a
shorter time. We also took a part of the wavelet with *P* < 0.7 d with
a duration of 5 d (see bottom panels for each sector). Wavelets show
QPOs with 0.2–0.3 d, which last approximately 1 d. There are also
QPOs with periods 0.07–0.1 d, which last for a part of the day.

There are also QPOs with long quasi-periods of *P* ≈ 7.5 – 8.8 d.<fnref n="5" /> They have a higher amplitude during *T* ≈ 10 – 15 d but can be seen at lower amplitudes throughout the observations. We discuss the possible origins of these quasi-periods later in the paper.

## 2.6 Multiwavelength observations of TW Hya

TW Hya was observed during a several weeks observational campaign by a set of ground-based telescopes<fnref n="6" /> in different wavebands (*u*, *B*, *V*, *g*, *r*, *i*, *z*) during two Epochs: Epoch 1 (2020–2021) and Epoch 2 (2022–2023; see details in Wendeborn et al. 2024b). Fig. 3 shows the light curves in the different photometric bands.

The Lomb–Scargle periodogram (see Fig. 4; taken from Wendeborn et al. 2024b) shows that in Epoch 1 (left panels) multiple periods were observed. However, the peak of QPO associated with the stellar period is smaller than the other peaks. In Epoch 2, the period of the star typically dominates. However, the peaks at approximately 6.5 and 8 d often have comparable amplitudes. A number of peaks with periods shorter than the period of the star are also observed.

Quasi-periods are expected to vary with time, and we performed
a Morlet-6 wavelet analysis of the light curves. For the wavelet
analysis, we took an interval of time from Epoch 2, where data were
recorded more frequently and regularly (mJD 59641–59767). We
subtracted *MJD* = 59641 d in each spectral band and obtained a
total of 109 d in bands *B*, *g*, *i*, and fewer *d* in other bands due to
restricted intervals of observations. We processed the light curves in
the following way. First, we removed outliers (the points where the
flux suddenly dropped, typically during one moment of time). Next,
we interpolated each light curve into a high-resolution uniform time
grid. For example the initial number of data points in *B*-band was 311.
In the new grid, the number of points is 4960 with a grid resolution
of 46 min. This grid reproduced high- and small-amplitude bursts
observed in the original light curve with high precision. Next, we
calculated the Morlet wavelets in each spectral band.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Morlet wavelet analysis of the light curves obtained by *TESS* during Sectors 9, 36, and 63 observations. Top panels: the light curves, obtained with homogeneous grid which include flagged points and interpolation. Second panels from top: Morlet wavelet obtained for light curves. Third panels from top: A part of wavelets for periods of *P* < 2 d. Bottom panels: A part of wavelets for periods of *P* < 0.7 d taken during the 5-d interval (marked with arrows in top panels).</figcaption></figure>

<footnote n="5">⁵The long-period QPOs may be connected with a short time-interval of observations. However, we observe this type of QPOs in wavelets obtained from much longer sets of ground-based observations (see Section 2.6).</footnote>

<footnote n="6">⁶Light curves were obtained with Las Cumbres Observatory Global Telescope (LCOGT), American Association of Variable Star Observers (AAVSO), and All-Sky Automated Survey for Supernovae (ASAS-SN). See table 1 of Wendeborn et al. (2024b) for details of observations.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top panel: optical light curves for TW Hya obtained during Epochs 1 (2021, left) and Epoch 2 (2022, right) of observations. Symbols corresponding to *u*, *B*, *g*, *V*, *r*, *i*, *z* data are labelled in the key. Data obtained by *TESS* are shown in Epoch 1. More details on the sources of the photometry are given in table 1 of Wendeborn et al. (2024b). Bottom panel: apart of Epoch 2 light curves (MJD 59660–59715) where the data were recorded more frequently. The vertical lines correspond to the times of the *HST* ULLYSES observations.</figcaption></figure>

Fig. 5 shows the wavelets and light curves in different wavebands. The wavelets look similar to those obtained from the *TESS* light curves. One can see a quasi-period of 3.2–4.0 d associated with the star’s rotation, which lasts 10–20 d, and also QPOs with smaller periods, which last shorter intervals of time. Oscillations with smaller periods last for shorter intervals of time. In addition, the wavelets show the presence of QPOs with a more extended period of 5.3–8.5 d, with typical periods of *P* ≈ 6.4 and *P* ≈ 8.5 d. These quasi-periods change with time and last 10–20 d. Similar long-period QPOs were observed in *TESS* wavelets.

Light curves in Fig. 5 show that the stellar magnitude decreases systematically with wavelength and the radiated energy increases systematically with wavelength. Fig. 6 shows the dependence of the time-averaged stellar magnitude on the wavelength of the wavebands. Vertical bars show the scatter in each waveband. One can see that the time-averaged stellar magnitude increases from 12.6 to 10.0 mag in Epoch 1 and from 12.1 to 9.9 mag in Epoch 2. The scatter gradually decreases from the *u*-band to bands with longer wavelengths.

## 3 NUMERICAL MODEL OF TW HYA

To model TW Hya, we took our earlier developed 3D MHD ‘cubed sphere’ code (Koldoba et al. 2002). The code has been used for modelling accretion on to magnetized stars with a tilted dipole magnetosphere (e.g. Romanova et al. 2003; Romanova & Owocki 2015), tilted dipole and rotational axes Romanova et al. (2021), and for modelling unstable regime of accretion (e.g. Romanova et al. 2008; Blinova et al. 2016). The numerical model is almost identical to that used in Blinova et al. (2016; see section 3 of their paper). Here, we briefly describe our numerical model and provide more details in Appendix A.

### 3.1 Set-up of the model and dimensionalization

We place a star of mass $M_*$ and radius $R_*$ to the centre of the simulation region. We place the cold, dense disc in the equatorial plane and a hot, rarefied corona above and below the disc and in
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left panels: Lomb–Scargle periodograms for TW Hya during Epoch 1. Right panels: the same but for Epoch 2. From Wendeborn et al. (2024b).</figcaption></figure>

the rest of the simulation region. The disc has an initial aspect ratio $h/r = 0.1$ determined at the inner edge of the disc, which is placed at a distance $R_d$ from the star. The density of the disc at $R_d$ is 100 times larger than that in the low-density region around the star. Initially, the disc and corona are in the rotational hydrodynamic equilibrium, where the corona above the disc rotates with the Keplerian velocity of matter in the disc (see Romanova et al. 2002 for details). This helps to decrease initial magnetic braking at the disc–corona boundary (which is present otherwise). A star has a dipole magnetic field with equatorial strength of $B_*$.

The magnetic moment of the dipole is tilted about the rotational axis of the star at an angle $\theta$. The rotational axis of the star coincides with that of the disc. The simulation region spreads from the inner boundary $R_{in}$ up to the outer boundary $R_{out} \approx 30R_{in}$. We used a grid with $N_r = 112$ grids in radial direction and $N \times N = 51 \times 51$ grids in each of six sides of the inflated cube. A set of 3D MHD equations has been solved using the Godunov-type method described in earlier works (Koldoba et al. 2002; Romanova et al. 2003). We describe initial and boundary conditions and other details of the numerical set-up in the Appendix A.

Equations are solved in dimensionless form using dimensionless parameters $\tilde{A}$. We introduce reference parameters $A_0$ and, after simulations, convert dimensionless values to dimensional ones. We take the reference scale, $R_0 = R_*/0.35$, where $R_*$ is the radius of the star; the reference mass, $M_0 = M_*$; Reference velocity is Keplerian velocity at $R_0$, $v_0 = (GM_0/R_0)^{1/2}$; period of rotation at $r = R_0$: $P_0 = 2\pi R_0/v_0$. To obtain the physical dimensional values A, the dimensionless values $\tilde{A}$ should be multiplied by the corresponding reference values $A_0$ as $A = \tilde{A}A_0$. For TW Hya, we take the following parameters: mass $M_0 = M_* = 0.79 M_*$, radius $R_* = 0.93 R_\odot$ and $R_0 = R_*/0.35 = 2.66 R_\odot$. Then reference velocity $v_0 = 238 \text{ km s}^{-1}$,

and reference period $P_0 = 0.565$ d. See other reference values in Appendix A.

In the model, we used the period of the star $P_* = 3.5$ d (which is close to period 3.56 d obtained by Setiawan et al. 2008 and Huélamo et al. 2008), and obtain the corotation radius $r_{cor} \approx 9.6 R_*$. We use it as a base for all simulations. In our model, we rotate a star with such a period that the corotation radius is located at 9.6 radii of the inner boundary. Another important parameter is the dimensionless magnetospheric parameter $\mu$ which we vary in the code and which determines the final magnetospheric radius (see more in Section A).

We calculated the kinetic energy flux at the stellar surface, $F = v_n \rho v^2/2$ (where $v$ is the total velocity and $v_n$ is its normal component), and assumed that all energy of the falling matter is converted into isotropic radiation. Then, we calculated the radiation flux towards the observer at two inclination angles: $i = 5^\circ$ and $i = 15^\circ$.

# 4 RESULTS

Simulations were performed at several values of the magnetospheric parameters $\mu$, which determined the final size of the magnetosphere. We calculated models at $\mu = 0.5, 1, 1.5, 2$. The size of the magnetosphere increases systematically with $\mu$. We obtain the smallest magnetospheres at $\mu = 0.5$ and the largest at $\mu = 2$. We obtain the size of the magnetosphere $r_m$ from simulations and show it in stellar radii. In all models, we obtain an unstable regime of accretion. In models with smaller $\mu$, we obtain a strongly unstable regime, where unstable tongues penetrate closer to the stellar surface, while in models with $\mu = 1.5$ and $2$, we obtain a mildly unstable regime, where unstable tongues are stopped by the inner magnetosphere at larger distances from the star, and matter accretes in funnel streams, going around the inner magnetosphere. We typically perform simulations at two tilts of the dipole moment, $\theta = 5^\circ$ and $20^\circ$, but we also calculate models with intermediate angles $\theta = 10^\circ$ and $15^\circ$ for comparisons. We show distances in stellar radii and time and periods in days. Other variables are shown in dimensionless form but can be converted to dimensional values using reference values from Table A1.

## 4.1 Small-sized magnetosphere, strongly unstable regime, model

First, we consider a model with magnetospheric parameter $\mu = 0.5$ and the tilt of the dipole magnetosphere $\theta = 20^\circ$. Fig. 7 shows the simulation results. The top panels show equatorial slices of density distribution in the inner part of the simulation region. Matter accretes in an unstable regime where equatorial ‘tongues’ carry matter from the inner disc towards the star. Two tongues dominate, which corresponds to an ordered unstable regime of accretion. The second from the top row of panels shows the density distribution and sample magnetic field lines in the $\mu-\Omega$ plane, where $\mu$ and $\Omega$ are unit vectors showing the directions of the magnetic moment and the star’s angular momentum, respectively. The third row of panels shows the flux distribution in hotspots viewed pole-on. Two spots are observed, and typically, one spot dominates because one tongue carries more matter than the other. Spots are shown in the pole-on projection.

In the unstable regime, the inner boundary varies its shape, so we determine the approximate value of the magnetospheric radius by placing a circle which corresponds to the sharp density drop at the disc–magnetosphere boundary (see also fig. A1 from Blinova et al. 2016). The dashed circles in the top panels show the location of the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Light curves and wavelets obtained from multiwavelength observations by multiple ground-based telescopes during a part of Epoch 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The time-averaged stellar magnitude versus wavelength of the band for different bands observed in Epochs 1 (circles and the lower, short-dashed curve) and 2 (diamonds and upper, long-dashed curve). Short- and long-dashed lines show the best fit for Epochs 1 and 2, respectively.</figcaption></figure>

magnetospheric boundary, which is $r_m \approx 3.2R_*$. At our corotation radius of $r_{cor} = 9.6R_*$, the ratio $r_m/r_{cor} \approx 0.33$ corresponds to a strongly unstable regime.

We calculated the radiation flux from spots towards the observer at angles of $i = 5^\circ$ and $i = 15^\circ$ (see left and right panels in the fourth row). The fifth row shows the Morlet wavelet spectrum of the light curves. At the inclination angle of $i = 5^\circ$ (left panel), variability with quasi-periods of 1.5–1.9 d dominates. In addition, quasi-periods of

0.8–1.3, 3.5, 4.4, 5.4, 7.2–7.5 d are present in the wavelets. At the observer's angle $i = 15^\circ$, QPOs with 1.5–1.9 d and other QPOs are also present. However, the QPOs with 3.3–3.4 d (close to the stellar rotational period) have also a high amplitude. The bottom panels show corresponding Fourier analysis of the light curves. The left panel shows that the period of the star is not present, but several shorter periods with maxima at 1.9, 1.6, 1.4, and 1.0 d dominate. In addition, there is a broader peak at 5.7 d and an even wider one at 7.7 d. The right panel shows that at $i = 15^\circ$ the main peak associated with the period of the star dominates, though peaks at the shorter and longer periods are also present.

This model shows that at a relatively small magnetosphere, an ordered unstable regime dominates (e.g. Romanova & Kulkarni 2009; Blinova et al. 2016) where matter accretes predominantly in one or two tongues that rotate with the angular velocity of the inner disc and produces short-period QPOs in light curves, wavelets, and Fourier spectrum.

We took a part of the light curve from the model $\mu 0.5\theta 20$ and considered three episodes where quasi-periods are observed (see top row of Fig. 8). The second row shows the corresponding wavelet where the QPO periods are marked. We see a good correlation between episodes of short-period QPOs observed in the light curves and wavelet. The bottom panels show that two tongues dominate during episodes 1 and 2, and one tongue dominates during episode 3. They originate due to the rotation of one or two unstable tongues in an ordered unstable regime (e.g. Romanova & Kulkarni 2009; Blinova et al. 2016).
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Results of simulations in the model $\mu=0.5\theta=20^\circ$ where the magnetospheric radius (measured in radii of the inner boundary) is the smallest. Top row: slices of density distribution in the equatorial plane at sample moments of time $t$. Second row from the top: the same, but in the $xz$-plane, the plane where the vectors of the angular momentum $\Omega$ and magnetic moment $\mu$ are located. Lines are sample magnetic field lines. *Third row:* the flux of energy distribution on the star’s surface as seen from the pole. *Fourth row:* fluxes towards the observer at inclination angles $i = 5^\circ$ (left) and $i = 15^\circ$ (right). *Fifth row:* Morlet wavelet of the observed fluxes. Bottom row: Fourier spectrum of the observed fluxes.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Figure explains the origin of the short-period QPOs in the model $\mu=0.5\theta=20^\circ$. Top row: a part of the light curve is shown where episodes of short-period QPOs are selected and marked with numbers 1, 2, and 3. The middle panel shows part of the wavelet up to a period of 3 d, where episodes of short-period QPOs are shown and marked. Numbers show an approximate value of quasi-period. The bottom row shows that two or one tongues dominate during these episodes of short-period QPOs.</figcaption></figure>
We calculated the radii in the disc corresponding to Keplerian rotation with QPO periods of $P = 1.8$ d (Sector 1 in the plot) and $P = 1.5$ d (Sector 3). These radii are $4.23R_*$ and $3.74R_*$ and are shown as dashed black circles in the bottom panel of Fig. 8. The dashed line shows the region of the disc where tongues form. These QPOs correspond to the rotation period of the strongest filament. Note that in Sector 3, the accretion rate is a bit higher, the inner disc comes closer to the star, and the period of QPOs is lower compared with Sector 1. The QPO in Sector 2 has twice as a lower period, $P = 0.8$ d, and can be explained by the rotation of two equal tongues forming at the inner disc at the radius $3.91R_*$ (corresponding to $P = 1.6$ d). The duration of QPOs is 5–7 d which corresponds to 3–4 Keplerian rotations of the inner disc. This model demonstrates the direct connection between QPO frequencies and Keplerian velocity and the location of the inner disc.

In TW Hya, this regime does not dominate in the currently presented light curves. However, in some earlier observations, the strong peaks at short-period QPOs and possible ‘stellar period’ of 1.3–2.8 d were suggested by other groups (e.g. Alencar & Batalha 2002; Lawson & Crause 2005; Siwak et al. 2011). Their observations may correspond to times when TW Hya switched to an ordered unstable regime.

In another example, Armeni et al. (2024) studied the photometric variability of another star – RU Lupi, which also shows the stochastic- looking light curves, and concluded that this star may accrete in the ordered unstable (in another terminology – the magnetic boundary layer regime – MBL; Romanova & Kulkarni 2009). In this star, the magnetospheric radius is expected to be small, $r_m \approx 2R_*$ (Armeni et al. 2024), and an ordered unstable/MBL regime is expected.

## 4.2 Mid-sized magnetospheres. The dependence on the dipole's tilt $\theta$

At the magnetospheric parameter $\mu = 1$, we obtain slightly larger magnetospheres with approximate magnetospheric radii of $r_m \approx 4.3$. We use this model to test the dependence of the variability on the tilt of the magnetic dipole. We calculated models with several tilts of the dipole: $\theta = 5^\circ, 10^\circ, 15^\circ$, and $20^\circ$.

### 4.2.1 Model $\mu 1\theta 5$ (tilt of the dipole: $\theta = 5^\circ$)

Fig. 9 shows the results of simulations in the model where the dipole is only slightly tilted about the rotational axis, $\theta = 5^\circ$. We observed that part of the time matter accretes in an unstable regime (with multiple tongues), and part of the time is in a strongly unstable regime (with two tongues). The light curves show some ordered oscillations with a quasi-period of 6.5 d and multiple short-period bursts with different quasi-periods of 2.8–3.1, 1.5–1.6, and shorter periods. In addition, the longer quasi-period of 8.2 d is present in the Fourier spectrum, but it is not seen in the wavelet, so it is insignificant. A quasi-period of 2.8–4.2 d is present in wavelet spectra but is not visible in Fourier spectra because the quasi-period changes with time. Note that spectra look similar at both inclination angles. We conclude that at a small tilt of the dipole, it is more tricky to find the period of the star from light curves.

### 4.2.2 Model $\mu 1\theta 20$ (tilt of the dipole: $\theta = 20^\circ$)

Fig. 10 shows the results of simulations in the model with an even larger tilt $\theta = 20^\circ$. The top three rows show a typical unstable regime of accretion. The light curve looks stochastic at an inclination of

$i = 5^\circ$ and more periodic at $i = 15^\circ$. Different quasi-periods are observed in both wavelets, and a quasi-period of 3.1–3.7 d associated with the star’s rotation is persistent at $i = 15^\circ$. Fourier analysis shows periods of 6.5, 3.5, and 4.1 d if $i = 5^\circ$ and a clear peak associated with the star’s rotation if $i = 15^\circ$.

Fig. 11 summarizes plots for Fourier spectra for models calculated for medium magnetospheres ($\mu = 1$). At $i = 5^\circ$ (left panels), the peak associated with the star’s rotation is not prominent, and it is typically smaller than other peaks. At $i = 15^\circ$ this peak is smaller than other peaks at tilts of the dipole of $\theta = 5^\circ$ and $\theta = 10^\circ$. However, it is seen and has the highest amplitude when $\theta = 15^\circ$ and $\theta = 20^\circ$. We suggest that different variability patterns in the light curves of TW Hya can be due to variation of the dipole tilt: at small tilts, the period of the star becomes less significant than other quasi-periods, and vice versa.

## 4.3 Larger magnetospheres

We also calculated models with a parameter $\mu = 1.5$ that provides larger magnetospheres. Here, we briefly discuss these models.

### 4.3.1 Model $\mu 1.5\theta 5$

Fig. A1 shows the results of simulations. The top row of panels shows that matter accretes in an unstable regime through multiple tongues. The magnetospheric radius is $r_m \approx 4.9R_*$. The second row shows that tongues are stopped at radii $(3-4)R_*$ by a strong dipole field, and matter falls from these radii to the stellar surface. The third row shows that multiple hotspots form at the star’s surface. Wavelets show different periods. Fourier analysis shows the period of the star at an inclination angle $i = 15^\circ$. However, a peak with a similar amplitude is observed at periods 5.8–7 d.

### 4.3.2 Model $\mu 1.5\theta 20$

Fig. A2 shows that accretion is still in the unstable regime. However, more matter flows above the main magnetosphere, like in a stable regime (see the second row of panels). Hotspots look irregular, like in other models. At $i = 15^\circ$, wavelet shows persistent quasi-period associated with rotation of the star. Fourier analysis shows a peak associated with the period of the star at $i = 15^\circ$ and no evident period at $i = 5^\circ$.

### 4.3.3 Model $\mu 2\theta 5$

At an even larger magnetosphere (obtained at $\mu = 2$), the magnetospheric radius is $r_m \approx 5.6R_*$. Fig. A3 of the Appendix shows similar unstable accretion but at a larger magnetosphere. Wavelet shows quasi-period close to the period of the star and other periods. Fourier shows several peaks at different periods at both inclination angles.

### 4.3.4 Model $\mu 2\theta 20$

At a larger tilt of the magnetosphere, $\theta = 20^\circ$ accretion occurs through both stable funnels and unstable tongues (see top two rows of Fig. A4). After $t > 30$ d, matter mainly accretes above the magnetosphere, and the flux becomes more stable. Fourier shows a clear peak associated with rotation of the star at $i = 15^\circ$, and several major peaks are observed at $i = 5^\circ$.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Same as in Fig. 7 but for the model $\mu=1\theta=5$.</figcaption></figure>

## 4.4 Variation of the magnetic field and photometric variability

Fig. 11 compares Fourier spectra for models with the same magnetospheric parameter $\mu = 1$ (corresponding to $r_m \approx 4.3R_*$) but at different tilts of the dipole magnetosphere $\theta$. One can see that at a very small tilt, $\theta = 5^\circ$ – there is no peak associated with the period of the star at both inclination angles of the observer. When $\theta = 10^\circ$ – the peak associated with the period of the star is present, but it is smaller than other peaks. At $\theta = 15^\circ$ and $20^\circ$, and inclination angle of $i = 15^\circ$, it dominates, but comparable with other peaks at $i = 5^\circ$. This result is important: if the tilt angle of the dipole moment slightly varies with time, then observers will detect or not detect the period of the star from their light curves.

Fig. 12 shows Fourier spectra for models with different magnetospheric radii (parameters $\mu$) and tilts of the dipole $\theta = 5^\circ$, and $\theta = 20^\circ$. The right panels of the figure show that at a larger tilt, $\theta = 20^\circ$, and $i = 15^\circ$ period of the star dominates over other QPO peaks. At $\theta = 5^\circ$ and $i = 15^\circ$, and larger magnetospheres ($\mu = 1.5$ and 2) the peak associated with the period of the star is larger than in model with smaller magnetosphere ($\mu = 1$), but still comparable with peaks of other QPOs. At a small inclination angle of $i = 5^\circ$ (two left columns), the peak associated with the period of the star has either a very small amplitude or is comparable with other peaks. Table 2 summarizes candidate ``periods'' and QPOs observed in Fourier spectra in different models. This analysis shows that in models with different sizes of the magnetosphere, the period of the star is less visible at a smaller tilt of the dipole magnetosphere.

Observations of TW Hya show that the wavelet and Fourier spectra vary from one observing season to the next and show different ‘periods,’ which are not real periods of the star. We suggest that the dipole component of the field changes from year to year: both the strength of the field and its tilt may vary with time due to convection processes inside the star. Our models show that the tilt of the dipole is the major factor: at small tilts $\theta \lesssim 10^\circ$, amplitudes of different QPOs are typically larger than the amplitude of the stellar period. At a larger tilt angle of the dipole, $\theta \gtrsim 10^\circ$, QPOs associated with the period of the star often dominate.

## 4.5 Best matching models

The above analysis shows that models with mid-size magnetospheres $4.3R_*$ and $4.9R_*$ show wavelets similar to those obtained with TESS and ground-based telescopes. These models were calculated using parameters $\mu = 1$ and $1.5$. However, at other times, models with smaller or larger magnetospheres may be applied. For example in years where the ‘period’ of $\sim 2$ d dominated, the magnetosphere may be smaller, like in our model $\mu = 0.5$ (where $r_m \approx 3.2R_*$). If the period of the star strongly dominates, and there are signs of high-velocity accretion on to the stellar surface and to higher latitudes (closer to the magnetic pole), then we expect a larger magnetosphere, which may correspond
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as in Fig. 7 but for the model $\mu 1\theta 20$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Fourier spectrum from light curves obtained in models with $\mu = 1$ but different tilt angles of the dipole $\theta = 5^\circ, \theta = 10^\circ, \theta = 15^\circ$, and $\theta = 20^\circ$. The observer's inclination angles are $i = 5^\circ$ (left panel) and $i = 15^\circ$ (right panel).</figcaption></figure>

to our model where $\mu = 2$ with $r_m \approx 5.6R_*$ or even larger.<fnref n="7" /> In models with $r_m/R_* = 3.2, 4.3, 4.9, 5.6$ and our corotation radius $r_{cor}/R_* = 9.6$, we obtain ratios $r_m/r_{cor} \approx 0.33, 0.45, 0.51, 0.58$, respectively. The strength of instability decreases when this ratio increases.

The magnetospheric radius depends on the strength of the magnetic field and the accretion rate $r_m \sim B_*^{4/7}/\dot{M}^{2/7}$. Different combinations of $B_*$ and $\dot{M}$ are possible. Table 3 shows combinations of these parameters for different sizes of the magnetosphere $r_m$. In the table, we calculate the accretion rate from the model using the time-averaged value of the dimensionless accretion rate and multiply it by $\dot{M}_0$ taken from Table A1. One can see that if the polar field is $B_* = 1,000$ G then the accretion rate should be $\dot{M} \approx 8.15 \times 10^{-9} \text{M}_\odot \text{yr}^{-1}$ or $\dot{M} \approx 2.0 \times 10^{-9} \text{M}_\odot \text{yr}^{-1}$ in models with $\mu = 1$ and $\mu = 1.5$, respectively. In the latter case, the accretion rate approximately corresponds to that found from optical observations of Herczeg et al. (2023) ($\dot{M} \approx 2.51 \times 10^{-9} \text{M}_\odot \text{yr}^{-1}$) and $r_m \approx 4.9R_*$. The former case ($\mu = 1$) may be relevant if the accretion rate is larger (see Section 2.4) and $r_m \approx 4.3$. Table 3 shows that if the magnetic field is smaller, then at the same accretion rates, $r_m$ is expected to be smaller, and the unstable regime stronger. For example at $B_* = 600$ G and $\dot{M} \approx 2.4 \times 10^{-9} \text{M}_\odot \text{yr}^{-1}$, the

<footnote n="7">⁷We note that it takes progressively longer time to simulate models with larger magnetospheres.</footnote>
<table><caption>**Table 2.** Main periods observed in Fourier spectra of different models at inclination angles of the observer *i* = 5° and 15°. Numbers in bold show the period of the star and periods close to this number.</caption><thead><tr><th>Model</th><th>Magn. moment</th><th>Tilt of the dipole</th><th><i>i</i> = 5°</th><th><i>i</i> = 15°</th></tr></thead><tbody><tr><td>μ0.5θ20</td><td>μ = 0.5</td><td>θ = 20°</td><td><i>P</i> = 5.7, 1.9, 1.6, 1.4, 1.0 d</td><td><i>P</i> = <b>3.5</b>, 5.7, 1.9, 1.6, 1.4, 1.0 d</td></tr><tr><td>μ1θ5</td><td>μ = 1</td><td>θ = 5°</td><td><i>P</i> = 6.5, 8.2, 2.8, 2.3 d</td><td><i>P</i> = 6.5, 8.2, 2.8, 2.3 d</td></tr><tr><td>μ1θ10</td><td>μ = 1</td><td>θ = 10°</td><td><i>P</i> = 6.0, 5.0, 3.0, 2.4, <b>3.5</b> d</td><td><i>P</i> = 6.0, 5.0, <b>3.5</b>, 2.4, 3.0 d</td></tr><tr><td>μ1θ15</td><td>μ = 1</td><td>θ = 15°</td><td><i>P</i> = 4.1, <b>3.5</b>, 2.7, 5.0, 6.1 d</td><td><i>P</i> = <b>3.5</b>, 4.2, 5.1, 6.1 d</td></tr><tr><td>μ1θ20</td><td>μ = 1</td><td>θ = 20°</td><td><i>P</i> = 6.5, <b>3.5</b>, 4.1 d</td><td><i>P</i> = <b>3.5</b>, 6.5, 4.1 d</td></tr><tr><td>μ1.5θ5</td><td>μ = 1.5</td><td>θ = 5°</td><td><i>P</i> = 5.8, 6.6, <b>3.5</b>, 3.2, 4.4 d</td><td><i>P</i> = <b>3.5</b>, 5.8, 6.6, 3.2, 4.4 d</td></tr><tr><td>μ1.5θ20</td><td>μ = 1.5</td><td>θ = 20°</td><td><i>P</i> = 7.2, 4.7, <b>3.7</b>, 3.2, 2.7, 1.7 d</td><td><i>P</i> = <b>3.5</b>, 7.2, 4.7, 2.7, 1.7 d</td></tr><tr><td>μ2θ5</td><td>μ = 2</td><td>θ = 5°</td><td><i>P</i> = 7.5, 5.7, 4.0, 2.9, 2.4, 1.5 d</td><td><i>P</i> = 7.5, 5.7, <b>3.5</b>, 2.9, 2.4, 1.5 d</td></tr><tr><td>μ2θ20</td><td>μ = 2</td><td>θ = 20°</td><td><i>P</i> = 7.9, <b>3.7</b>, 4.3, 2.8, 1.8 d</td><td><i>P</i> = <b>3.5</b>, 7.9, 4.3, 2.8, 1.8 d</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Fourier spectrum from light curves obtained in different models and shown at the inclination angle *i* = 5° (left two rows of panels) and *i* = 15° (right two rows of panels).</figcaption></figure>

<table><caption>**Table 3.** The accretion rate $\dot{M}(M_{\odot} \text{yr}^{-1})$ obtained in models with different values of the magnetospheric parameter $\mu$ and equatorial magnetic field of the star $B_{\text{eq}}$.</caption><thead><tr><th>B-field (Gauss)<br>polar</th><th>B-field (Gauss)<br>equatorial</th><th>μ = 0.5<br>r<sub>m</sub> ≈ 3.2R<sub>*</sub></th><th>μ = 1.0<br>r<sub>m</sub> ≈ 4.3R<sub>*</sub></th><th>μ = 1.5<br>r<sub>m</sub> ≈ 4.9R<sub>*</sub></th><th>μ = 2.0<br>r<sub>m</sub> ≈ 5.6R<sub>*</sub></th></tr></thead><tbody><tr><td>B<sub>*</sub> = 1,000</td><td>B<sub>eq</sub> = 500</td><td>6.64 × 10<sup>-9</sup></td><td>8.15 × 10<sup>-9</sup></td><td>2.00 × 10<sup>-9</sup></td><td>1.95 × 10<sup>-9</sup></td></tr><tr><td>B<sub>*</sub> = 800</td><td>B<sub>eq</sub> = 400</td><td>4.26 × 10<sup>-9</sup></td><td>2.33 × 10<sup>-9</sup></td><td>1.28 × 10<sup>-9</sup></td><td>1.24 × 10<sup>-9</sup></td></tr><tr><td>B<sub>*</sub> = 600</td><td>B<sub>eq</sub> = 300</td><td>2.40 × 10<sup>-9</sup></td><td>1.31 × 10<sup>-9</sup></td><td>7.19 × 10<sup>-10</sup></td><td>7.02 × 10<sup>-10</sup></td></tr></tbody></table>

magnetospheric radius $r_m \approx 3.2R_*$, and stronger unstable regime is expected.

Note that these results depend on the radius of the star. In our models, we took $R_* = 0.93 R_\odot$, and fixed the value $r_{\text{cor}}/R_* = 9.6$. If the stellar radius is larger (see Table 1) then the results are different. For example if we take $R_* = 1.16 R_\odot$, and $M_* = 0.8 M_\odot$ (Baraffe et al. 2015), then the corotation radius is smaller $R_{\text{cor}} \approx 7.9 R_*$. Test simulations at this radius and $\mu = 1$ and 1.5 show, that the magnetospheric radius is approximately the same as in models with larger $R_{\text{cor}}$. However, the ratios are larger: $r_m/r_{\text{cor}} \approx 0.40, 0.54, 0.62, 0.71$ in models with $\mu = 0.5, 1, 1.5, 2$, respectively, and unstable regime is not as strong as in our main models.

## 4.6 Long-period QPOs: waves in the disc

The light curves also show the longer period QPOs with quasi-periods of 6.4–6.6, 7.5, and 8.2 d in models $\mu1\theta5$, $\mu1\theta10$, and $\mu1\theta15$, respectively. These QPOs can be explained by waves in the accretion

disc. 3D MHD simulations show that a star with a tilted dipole magnetosphere excites different types of density and bending waves in the accretion disc (Romanova et al. 2013). Simulations show that if the magnetospheric radius is comparable with the corotation radius, then the bending wave forms near the closed rotating magnetosphere and rotates with an angular velocity of the star (see figs 1–3 from Romanova et al. 2013).<fnref n="8" /> However, if the magnetospheric radius is much smaller than the corotation radius (like in TW Hya), then slowly rotating waves form at the outer Lindblad resonance (see fig. 12 of the same paper for 3D and 2D plots and figs 13–15 for the location of these waves at different corotation radii).

We demonstrate the presence of such a wave using a model $\mu1\theta5$ where the slowly rotating wave is clearly observed. The light curve in Fig. 9 shows an episode of long-period QPO with a quasi-period

<footnote n="8">⁸Formation of a bending waves has been predicted theoretically (e.g. Terquem & Papaloizou 2000) and proposed to explain variability in AA Tau (e.g. Bouvier et al. 1999, 2007b).</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Top panel: a part of the light curve shown in Fig. 9 for the model µ1θ5, where the long-period QPO is observed. Middle panel: wavelet corresponding to this part of the light curve where we concentrate on the long period, 6.3–6.5 d QPO. The bottom row of panels: density distribution in the equatorial plane at different moments in the interval of 28.2–34.5 d, which approximately correspond to one rotation of the large-scale density wave about the star. The wave is marked with a small circle.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Accretion rate on to the star $\dot{M}$ (top curve) and the observed flux from spots $\dot{E}_{\text{obs}}$ (bottom curve). The latter was multiplied by the factor 7 to show them side by side. Vertical lines show approximate positions of maxima.</figcaption></figure>

of 6.4–6.6 d. We analysed this time interval in greater detail (see Fig. 13). We took one of the quasi-periods and plotted the density distribution in the disc frequently (see bottom panels of Fig. 13). We detected the presence of a density wave, which had approximately one complete rotation during this interval of time. The light curve originates from accretion on to the star. That is why we checked the variability in the accretion rate. Fig. 14 shows that the accretion rate also has oscillations with approximately the same quasi-period. It means that the remote density wave modulates the accretion rate, which leads to variation in the light curve. This analysis shows that the long-period QPOs observed in light curves of TW Hya can be caused by large-scale density wave that forms beyond the corotation radius. It is possible that many long-period QPO features are connected with large-scale density or bending waves in the disc.

Donati et al. (2024) reports the detection of a radial velocity signal of semi-amplitude 11.1 m s⁻¹ at a period of 8.3 d in the spectrum of TW Hya. They suggested that the origin may be attributed to either a non-axisymmetric density structure in the inner accretion disc or a 0.55 M☉ candidate close-in planet at an orbital distance of 0.075 ± 0.001 au. We suggest that this period can be associated with the wave in the disc. However, the period associated with density wave may vary in time (e.g. due to inhomogeneities in the disc). If this period has a high precision, then the hypothesis of the planet is much more probable.

## 4.7 Properties of hotspots

Unstable tongues propagate some distance in the equatorial plane (pushing magnetic field lines aside), then encounter the stronger field of the inner magnetosphere (at some radius $r_{\text{m,in}}$) and are lifted above and below the magnetosphere, forming short-lived funnel-looking streams. They deposit matter at some distance from magnetic poles. In models with smaller magnetospheres, tongs are lifted at $r_{\text{m,in}} \approx (1-2)R_*$, and spots are located at $\sim (35^{\circ} - 45^{\circ})$ from magnetic poles (see Fig. 15 for a model with $\mu = 0.5$). At larger magnetospheres, they are lifted at larger distances from the star, $r_{\text{m,in}} \approx (2-4)R_*$, and spots are located at $\phi \lesssim 30^{\circ}$ from magnetic poles (see Fig. 16 for model with $\mu = 1$).

Matter accreting from temporary funnel streams is accelerated by gravity and falls to a star with a free-fall velocity. If the matter falls from the distance $r_{\text{m,in}}$ then its velocity at the surface of the star is

$$V_{\text{ff}} = \left[ \frac{2GM_*}{R_*} \left( 1 - \frac{R_*}{r_{\text{m,in}}} \right) \right]^{1/2} \\ \approx 569 \frac{km}{s} \left[ \left( \frac{M_*}{0.79 M_{\odot}} \right) \left( \frac{R_*}{0.93 R_{\odot}} \right) \left( 1 - \frac{R_*}{r_{\text{m,in}}} \right) \right]^{1/2} . \quad (3)$$

If matter falls from the distance of $r_{\text{m,in}} = 2R_*$ then $V_{\text{ff}} \approx 402$ km s⁻¹. It is 464 km s⁻¹ in case of $r_{\text{m,in}} = 3R_*$, and 493 km s⁻¹ in case of $r_{\text{m,in}} = 4R_*$.

The spots have an inhomogeneous structure: the energy flux per unit area is larger in the central regions of the spots and smaller on the outskirts. The inhomogeneity of the energy distribution in hotspots has been predicted in 3D MHD simulations of stars accreting in stable regime (Romanova et al. 2004) and confirmed in observations by Espaillat et al. (2021) and Singh et al. (2024). In the unstable regime, spots are also inhomogeneous (see bottom panels in Figs 15 and 16). The area covered by the highest energy flux is much smaller than that covered by the lower energy fluxes. This could be reflected in fluxes of light curves observed in different wavebands. Fig. 6 shows that the stellar magnitude is the largest (and the flux is the

<footnote n="9">⁹Zhu et al. (2024) obtained the similar result in independent 3D MHD simulations.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Top panels: 3D plots of matter flow in the model $\mu0.5\theta20$ at sample times. The background shows one of the density levels. Lines are sample field lines. Bottom panels: distribution of the energy flux at the star's surface observed at an inclination angle $i = 90^\circ$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Same, but in the model $\mu1\theta10$.</figcaption></figure>

smallest) in the *u* band. The flux systematically increases with the wavelength. It is the largest in *i* and *z* wavebands.

The stochasticity in light curves arises from the fact that tongues and funnel streams are short-lived features: they deposit matter on to the stellar surface rapidly and form temporary hotspots on the surface. Typically, more powerful tongues live longer and provide longer living spots on the surface of the star, while less powerful tongues form spots that live for a shorter period of time. This may explain QPOs observed in wavelets: QPOs with larger power last longer, often a few days, while QPOs with smaller power last shorter time.

# 5 DISCUSSION

## 5.1 Comparisons with MOST-ASAS observations

Rucinski et al. (2008) and Siwak et al. (2011, 2014, 2018) analysed results of TW Hya observations obtained with MOST and ASAS telescopes. Most of their results are similar to results obtained in our analysis of light curves obtained with TESS, ground-based telescopes, and in numerical simulations. The common feature of all wavelets is the presence of QPOs with different quasi-periods. The duration of QPOs decreases when the quasi-period decreases. Rucinski et al. (2008) and Siwak et al. (2014, 2018) noted that in a

few instances period of QPOs decreases with time. Our observational and simulation data show that the period of QPOs may decrease, or increase, or be approximately the same. Additional work is required to track QPO periods with similar amplitudes, which can be done in the future. Siwak et al. (2011) noted that a shorter period oscillations are often observed during higher amplitude bursts (see e.g. fig. 1 from their paper). We observed similar features in Section 9 and 36 of TESS light curves but less so in Sector 63. In numerical simulations, this phenomenon is observed in models with larger sized magnetospheres, e.g. in models $\mu2\theta5$ and $\mu2\theta20$. Siwak et al. (2018) noted that the time duration of QPOs may be equal to a few dynamical rotations. Using one of our models ($\mu0.5\theta5$) where 1–2 unstable tongues were observed, we demonstrate that QPO periods can be associated with 3–4 Keplerian rotations of the inner disc, which is in accord with their hypothesis. In models with a larger number of tongues, the picture of QPOs is more complex, though the strongest tongue may develop a QPO with the period of the inner disc.

## 5.2 Comparisons with other works

### 5.2.1 Donati et al. (2011, 2024)

Donati et al. (2011, 2024) measured the magnetic field in TW Hya in different observational sets, and obtain different strengths and
tilts of the dipole component. In more recent paper, (Donati et al. 2024, hereafter D24) they have shown that the magnetic field of TW Hya represents a slightly tilted dipole field with a strength of $B_d = 990-1190$ G. The left panels of fig. 5 of their paper show that the radial component of the field (its stronger part) is located at $\phi \approx 30^\circ - 40^\circ$ off the pole. The authors estimated the magnetospheric radius $r_m \approx 4.5_{-2.0}^{+1.0} R_*$, using the formulae of Bessolaz et al. (2008) ($B = 1.1$ kG and $\log(\dot{M}) = -8.65$). They concluded that matter flows from the inner disc to the star in a stable funnel stream and accretes close to the magnetic pole. Our simulations (see Table 3) show that at comparable values of $B_d = 1$ kG and similar accretion rate, $\dot{M} = 2 \times 10^{-9} \dot{M}_{\odot}$ yr$^{-1}$ the magnetospheric radius is $r_m \approx 4.9 R_*$ which is in a good agreement with that derived by D24. Our model shows that at these radii, accretion is unstable. However, in models with large magnetospheres (at $\mu = 1.5$ and 2.0), the unstable tongues are lifted above the inner magnetosphere (like in a stable regime) and form spots at high latitudes within $\phi \sim 30^\circ - 40^\circ$ off the magnetic pole. However, instead of one funnel stream, there are several streams and spots that reside at relatively high latitudes. More detailed observations are needed to distinguish between one or several accretion spots located at similar latitudes.

For estimates, D24 accepted parameters: $M_* = (0.8 \pm 0.01) M_{\odot}$, $R_* = (1.16 \pm 0.13) R_{\odot}$ (Baraffe et al. 2015) and corotation radius $r_{\text{cor}} = (7.9 \pm 0.3) R_*$. We performed special simulation runs using their corotation radius, $r_{\text{cor}} \approx 7.9 R_*$ and our model parameters $\mu = 1$ and 1.5. We obtain from simulations $r_m \approx 4.3 R_*$ at $\mu = 1$ and $r_m \approx 4.8 R_*$ at $\mu = 1.5$ (this radius is expected to be even larger if $\mu = 2$, see trend in Table 3). These radii are close to those obtained in our main models (where $r_{\text{cor}} = 9.6 R_*$). The ratios are $r_m/r_{\text{cor}} \approx 0.54$ and 0.61, respectively. They are slightly larger than in our main models and correspond to milder unstable regimes. In a milder regime, the unstable tongues are lifted above the inner magnetosphere at larger distances from the star, and spots form closer to the magnetic pole at $\phi \lesssim 30^\circ$ from the pole. These results are in line with observations of D24.

In earlier observations of TW Hya, Donati et al. (2011) derived predominantly octupolar magnetic field with a much smaller dipole component of a few hundred Gauss. At a weaker dipole component, the magnetospheric radius is smaller, and a star is expected to be in a strongly unstable regime. The light curve is expected to be more stochastic, and spots form further away from the magnetic pole. We do not possess the detailed light curve during their observational period. Future observations of the magnetic field structure with simultaneous photometric observations would be interesting for further understanding the dependence of the light curve on the structure of the dipole magnetic field.

### 5.2.2 Sicilia et al. (2023)

Sicilia et al. (2023) analysed several spectral lines that are expected to form near the footprint of the funnel flow (such as He I 5016 A and Fe II 5018) and found a periodic signal in radial velocities with a period of the star. They concluded that matter producing these spectral lines flows in a stable funnel stream, while stochastic photometric variability is produced by some spots that form in random places and have a random strength. We suggest that both phenomena can be explained by a model of accretion in a mildly unstable regime, in which spots form at a higher latitude and provide chaotic components in the light curve. At the same time, due to the tilt of the dipole and the inclination angle of the observer, the whole set of spots (located around the magnetic pole) rotates about

the rotational axis of the star with the angular velocity of the star and produces periodic components in spectra. We note that in fig. 5 of Sicilia et al. (2023; the folded phase-velocity plot for He I 5016 A and Fe II 5018 lines) there is a significant scatter of velocities, which may be connected to multiple funnel streams that fall near the magnetic pole but not in the same location.

### 5.2.3 GRAVITY Collaboration (2020)

GRAVITY Collaboration (2020) used interferometry to measure the size of the magnetosphere and radius of the TW Hya star (see details in their paper). They measured the radius of the Brγ emitting region as $R_{Br\gamma} = (3.49 \pm 0.20) R_*$ and suggested that this radius corresponds to the magnetospheric radius $r_m$. They measured the radius of the star as $R_* = (1.29 \pm 0.19) R_{\odot}$. At this radius (and other parameters: $M_* = 0.8 M_{\odot}$, $P_* = 3.566$ d), we obtain a corotation radius $r_{\text{cor}} = 7.06 R_*$ that is smaller than that accepted in our paper ($r_{\text{cor}} = 9.6 R_*$). We calculated an additional model taking $r_{\text{cor}}/R_* = 7.06$ and magnetospheric parameter $\mu = 1$. We observed an unstable regime of accretion. From simulations, we measured the magnetospheric radius $r_m \approx 3.9 R_*$ and the ratio $r_m/r_{\text{cor}} \approx 0.63$. This corresponds to a mildly unstable regime.

## 5.3 More complex fields

The magnetic field of TW Hya may have strong octupolar and higher order components of the field (Donati et al. 2011, 2024). Earlier 3D MHD simulations of models where the octupole component has been included show that initially, matter flows from the disc to the star along the dipole field, but closer to the star the octupole component redirects matter to a new position on the surface of the star determined by the octupole field (e.g. Long et al. 2008, 2011, 2012; Romanova et al. 2011; see also Gregory et al. 2010). In this case, the interchange instability at the disc-magnetosphere boundary will also produce unstable tongues and multiple flares in the variability curve. However, the location of hotspots will be determined by the octupolar field. In another research, the authors show evidence of accretion on to a low-latitude, almost equatorial spot using X-ray observations and argue that matter may accrete to one of the poles of the complex field (Argiroffi et al. 2017). It is possible if the dipole component is very weak (e.g. Gregory et al. 2010). Or if some unstable tongues penetrate deeper into a small dipole magnetosphere. A more detailed answer requires modelling accretion to a star with a dipole plus a more complex field, which can be done in the future.

## 5.4 Other possible origins of stochasticity

Stochastic-looking light curves can also originate in stable two-funnel magnetospheric accretion if matter flowing into the funnel has inhomogeneities of different sizes caused by turbulence. Robinson et al. (2021) used a one-dimensional accretion model from a turbulent disk, showing that it may explain stochastic accretion (see also Robinson et al. 2017). This model requires more detailed 3D simulations with well-resolved turbulence in the disc.

3D MHD simulations of MRI-driven turbulent matter on to a star with a very weak (dynamically unimportant) magnetic field (where the disc matter accretes on to the surface of the star in the boundary layer regime) show that the disc has turbulent cells of different sizes (see figs 3 and 4 from Romanova et al. 2012; hereafter R12). The accretion rate curve shows high-amplitude peaks associated with the accretion of the largest turbulent cells and more frequent small-amplitude peaks associated with smaller sized cells (see fig. 6
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Light curves in the model with $\mu = 1$ and tilts of the dipole $\theta = 10^{\circ}$ (left panel) and $\theta = 20^{\circ}$ (right panel) shown for different inclination angles of the observer: $i = 5^{\circ}$, $15^{\circ}$, and $30^{\circ}$. Fluxes are scaled for better visibility.</figcaption></figure>

from R12). In similar models but with larger, dynamically important magnetospheres, matter accretes predominantly through the large-scale azimuthally elongated turbulent cells if the tilt of the dipole is large $\theta = 30^{\circ}$ (see figs 7–9 from R12). In models with a small tilt of the dipole $\theta = 2^{\circ}$, turbulent matter penetrates the magnetosphere through interchange instability (see right panel of fig. 7 from R12).

The irregular photometric variability of CTTS can also be caused by variable circumstellar extinction. Such obscuration by dust is more typical for stars with a high inclination angle, where the line of sight intersects with the dusty disc winds or warp in the disc (e.g. Bouvier et al. 2007a). However, TW Hya is observed almost pole-on (e.g. Qi et al. 2004), and therefore, variability associated with dust extinction is expected to be insignificant.

## 5.5 Possible role of winds. Two-channel accretion

Many spectral lines of TW Hya show blueshifted absorption (e.g. Alencar & Batalha 2002; Lamzin et al. 2004; Herczeg et al. 2023; Wendeborn et al. 2024cc), which points to the presence of winds. Winds can originate at the surface of the star (e.g. Matt & Pudritz 2005) and may form polar winds (e.g. Cranmer 2009; Dupree et al. 2012) or at the disc–magnetosphere boundary (e.g. Shu et al. 1994; Goodson et al. 1997; Romanova et al. 2009; Zanni & Ferreira 2013; Lii et al. 2014). The variability in blueshifted absorption typically does not correlate with variability in spectral lines or with photometric variability (e.g. Alencar & Batalha 2002). Donati et al. (2024) notes that no periodicity was observed in the blue wing of the He 1083 nm line (see fig 8 in Donati et al. 2024), which means that it either comes from an erratic stellar wind or the inner disc. The origin of winds is yet to be understood.

Recent 3D MHD simulations of accretion on to a magnetized star confirmed accretion through instabilities (e.g. Takasao et al. 2022; Parfrey & Tchekhovskoy 2024; Zhu et al. 2024). Zhu et al. (2024) used the code with a high grid resolution (Cartesian with mesh refinement) and have shown that unstable filaments have a finer substructure which have not been seen in this work. This may explain the short-period QPOs observed by the TESS telescope. Takasao et al. (2022) report on observations of unstable accretion at a wider range of $r_m/r_{\text{cor}}$ compared with condition $r_m/r_{\text{cor}} \lesssim 0.71$ found by Blinova et al. (2016). They observed inflation and reconnection of the field lines connecting a star and the disc (e.g. Lovelace et al. 1995). Takasao et al. (2022) observed the formation of conical-shaped winds in the case of faster rotating stars, which may lead to large-scale

outflows. However, ejections are less ordered than in axisymmetric simulations (e.g. Zanni & Ferreira 2013; Lii et al. 2014). In cases of slowly rotating stars, Takasao et al. (2022) repot on the formation of slow turbulent winds which fail to live the system and accrete back to the star. These ‘failed’ winds may be another source of variability in CTTSs.

Simulations also point to two-channel accretion, where some matter accretes in several unstable tongues in the equatorial plane while some matter accretes from the top layers of the disc in funnel streams (Zhu et al. 2024). Earlier, Bachetti et al. (2010) studied two-channel accretion in 3D simulations and applied it to accreting millisecond pulsars (which represent a scaled version of CTTSs). They noticed that in the unstable regime moving spots at the stellar surface may produce two frequencies: the higher frequency is associated with unstable tongues, driven by the inner disc, and the lower frequency is caused by the moving spots, resulting from funnel streams originating at larger distances from the star. They used this model to explain two QPO frequencies observed in accreting millisecond pulsars. In our current models, we also see both accretion through instabilities and funnel accretion, which starts at the surface layers of the disc. The latter becomes more significant in models with milder unstable regime and larger sized magnetospheres.

## 5.6 Projection to other CTTSs

TW Hya is seen almost pole-on. In models, we can look at the star from different observer angles. We observed that at larger inclination angles, the light curve becomes more periodic or quasi-periodic. At the same time, the unstable accretion continues, providing QPO flares on different time-scales. For example, the light curve in the model $\mu = 1$ with the tilted angle of the dipole $\theta = 20^{\circ}$ looks stochastic at $i = 5^{\circ}$, but at $i = 15^{\circ}$, and $i = 30^{\circ}$, it becomes more quasi-periodic (see right panel of Fig. 17). On the other hand, at a small tilt of the dipole $\theta = 10^{\circ}$, the light curve looks stochastic at all above inclination angles (see left panel of the same plot). It becomes more periodic at larger inclination angles. Overall, we expect that the period of the star becomes more visible when the inclination angle increases because the whole set of spots is located around magnetic poles and rotates around the rotational axis together with the magnetic axis (see also Stauffer et al. 2014). At the same time, the stochastic component stays.

If a star is seen at high inclination, one can expect to observe stochastic occultations of the photosphere by dust lifted above the disc plane near the base of the accretion tongues whenever one passes
in front of our line of sight, leading to light curves dominated by aperiodic extinction events (e.g. McGinnis et al. 2015; Petrov et al. 2019).<fnref n="10" />

## 6 CONCLUSIONS

The conclusions of our work are the following:

(i) Wavelet analysis of light curves obtained with the TESS telescope in Sectors 9, 36, and 63 show multiple QPOs ranging from 8–9 d to less than 2 h. QPOs with periods close to the period of the star (3.56 d) are often present where the period varies in the range of 3 to 5 d. QPOs with longer periods of 7.5–8.8 d are always present. Multiple short-period QPOs are present and last for a shorter time.

(ii) A subset of ODYSSEUS light curves covering 70 d in 2022 (part of Epoch 2) was selected due to more frequent observations and analysed in detail. The light curve looks stochastic. Wavelet analysis shows multiple QPOs similar to those observed with TESS. QPOs with the stellar period are seen in both wavelet and Lomb–Scargle periodograms. No stellar period is seen in Epoch 1.

(iii) We developed 3D MHD models of a star with TW Hya parameters and several magnetosphere sizes from $3.2R_*$ to $5.6R_*$. In all models, matter accretes in the unstable regime, where unstable tongues form moving hotspots, which produce stochastic-looking light curves and QPOs with different quasi-periods. In some models, the light curves and wavelets are strikingly similar to those obtained from observations by TESS and ground-based telescopes.

(iv) In models with smaller tilts of the dipole magnetosphere $\theta = 5°$ and $10°$, variability is more stochastic, and the amplitude of QPO associated with the stellar period is either smaller or comparable to other QPOs. At larger tilts $\theta = 15° – 20°$, the QPO associated with a stellar period typically dominates in Fourier and wavelet spectra.

(v) The magnetic field of TW Hya varies from year to year (e.g. Donati et al. 2011, 2024). We suggest that the difference in light curves obtained in different years may be connected with variations in the strength and tilt of the dipole field.

(vi) In a model with a small magnetosphere ($r_m \approx 3.2R_*$), matter accretes in one or two funnel streams, which rotate with period of the inner disc. This quasi-period may be mistakenly accepted to be period of the star.

(vii) We show that persistent long-period QPOs can be caused by modulation of matter flow by density waves that form beyond the corotation radius.

(viii) Unstable tongues are stopped by the inner parts of the magnetosphere, form funnel-like streams, and hit the star at an angle of $\phi \sim 20° – 30°$ (from magnetic pole) in models with larger magnetospheres, and $\phi \sim 35° – 45°$, in models with smaller magnetospheres.

(ix) Light curves may have both: properties of the ordered magnetospheric accretion (due to matter flow around the closed magnetosphere) and stochasticity due to accretion of multiple tongues. In cases of a larger inclination angle of the observer, the light curves may become more ordered but with a stochastic component.

## ACKNOWLEDGEMENTS

The authors thank Dr. Michal Siwak for a careful reading of our paper and the insightful, helpful report. Resources supporting this work

were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and the NASA Center for Computational Sciences (NCCS) at Goddard Space Flight Center. MMR and RVEL were supported in part by the NSF grant AST-2009820. They thank Dr. Alexander Koldoba for earlier developed codes and Bez Thomas for technical support. CCE and JW were supported by NSFAST-2108446, HST AR-16129, and NASAADAP 80NSSC20K0451.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author (MMR).

## REFERENCES

* Alencar S. H. P., Batalha C., 2002, ApJ, 571, 378

* Alencar S. H. P. et al., 2010, A&A, 519, A88

* Andrews S. M. et al., 2016, ApJ, 820, L40

* Argiroffi C., Drake J. J., Bonito R., Orlando S., Peres G., Miceli M., 2017, A&A, 607, A14

* Armeni A. et al., 2024, A&A, 690, A225

* Arons J., Lea S. M., 1976, ApJ, 207, 914

* Arulanantham N. et al., 2023, ApJ, 944, 185

* Bachetti M., Romanova M. M., Kulkarni A., Burderi L., di Salvo T., 2010, MNRAS, 403, 1193

* Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42

* Batalha C., Batalha N. M., Alencar S. H. P., Lopes D. E., Duarte E. S., 2002, ApJ, 580, 343

* Bessolaz N., Zanni C., Ferreira J., Keppens R., Bouvier J., 2008, A&A, 478, 155

* Blinova A. A., Romanova M. M., Lovelace R. V. E., 2016, MNRAS, 459, 2354

* Bouvier J. et al., 1999, A&A, 349, 619

* Bouvier J. et al., 2007, A&A, 463, 1017

* Bouvier J., Alencar S. H. P., Harries T. J., Johns-Krull C. M., Romanova M. M., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. Univ. Arizona Press, Tucson, p. 479

* Calvet N., D'Alessio P., Hartmann L., Wilner D., Walsh A., Sitko M., 2002, ApJ, 568, 1008

* Cody A. M., Hillenbrand L. A., 2018, AJ, 156, 71

* Cody A. M. et al., 2014, AJ, 147, 82

* Cody A. M., Hillenbrand L. A., Rebull L. M., 2022, AJ, 163, 212

* Cranmer S. R., 2009, ApJ, 706, 824

* Donati J.-F., 2003, in Trujillo-Bueno J., Almeida J. S., eds, ASP Conf. Proc., Vol. 307, Solar Polarization. Astron. Soc. Pac., San Francisco, p. 41

* Donati J.-F., Landstreet J. D., 2009, Ann. Rev. Astron. Astrophys., 47, 333

* Donati J.-F. et al., 2011, MNRAS, 417, 472

* Donati J.-F. et al., 2020, MNRAS, 498, 5684

* Donati J.-F. et al., 2024, MNRAS, 531, 3256

* Dupree A. K. et al., 2012, ApJ, 750, 73

* Elsner R. F., Lamb F. K., 1977, ApJ, 215, 897

* Espaillat C. C., Robinson C., Grant S., Reynolds M., 2019, ApJ, 876, 121

* Espaillat C. C., Robinson C. E., Romanova M. M., Thanathibodee T., Wendeborn J., Calvet N., Reynolds M., Muzerolle J., 2021, Nature, 597, 41

* Espaillat C. C. et al., 2022, AJ, 163, 114

* Fischer W. J. et al., 2023, in Inutsuka S.-i., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 355

* France K., Schindhelm R., Bergin E. A., Roue E., Abgrall H., 2014, ApJ, 784, 127

* Goodson A. P., Bohm K.-H., Winglee R. M., 1997, ApJ, 489, 199

* Gravity Collaboration, 2020, Nature, 584, 547

* Gregory S. G., Jardine M., Gray C. G., Donati J.-F., 2010, Rep. Prog. Phys., 73, 126901

<footnote n="10">¹⁰Earlier, this idea was proposed by Bouvier et al. 1999 for stellar light obscuration by inner bending wave in the disc.</footnote>
Günther M. N. et al., 2020, AJ, 159, 16

Hartmann L., 2000, Accretion Processes in Star Formation, Vol. 32. Cambridge Univ. Press, Cambridge

Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135

Herbst W., Koret D. L., 1988, AJ, 96, 1949

Herbst W., Herbst D. K., Grossman E. J., Weinstein D., 1994, AJ, 108, 1906

Herbst W., Bailer-Jones C. A. L., Mundt R., Meisenheimer K., Wackermann R., 2002, A&A, 396, 513

Herczeg G. J., Linsky J. L., Valenti J. A., Johns-Krull C. M., Wood B. E., 2002, ApJ, 572, 310

Herczeg G. J. et al., 2023, ApJ, 956, 27

Huélamo N. et al., 2008, A&A, 489, L9

Hughes A. M., Wilner D. J., Calvet N., D'Alessio P., Claussen M. J., Hogerheijde M. R., 2007, ApJ, 664, 536

Ingleby L. et al., 2013, ApJ, 767, 112

Johns-Krull C. M., 2007, ApJ, 664, 975

Johnstone C. P., Jardine M., Gregory S. G., Donati J. F., Hussain G., 2014, MNRAS, 437, 3202

Kastner J. H., Huenemoerder D. P., Schulz N. S., Weintraub D. A., 1999, ApJ, 525, 837

Kastner J. H., Huenemoerder D. P., Schulz N. S., Canizares C. R., Weintraub D. A., 2002, ApJ, 567, 434

Koldoba A. V., Romanova M. M., Ustyugova G. V., Lovelace R. V. E., 2002, ApJ, 576, L53

Königl A., 1991, ApJ, 370, L39

Kulkarni A., Romanova M. M., 2013, MNRAS, 433, 3048

Kurosawa R., Romanova M. M., 2013, MNRAS, 431, 2673

Lamb F. K., Pethick C. J., Pines D., 1973, ApJ, 184, 271

Lamzin S. A., Kravtsova A. S., Romanova M. M., Batalha C., 2004, Astron. Lett., 30, 413

Lawson W. A., Crause L. A., 2005, MNRAS, 357, 1399

Lii P. S., Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E., 2014, MNRAS, 441, 86

Long M., Romanova M. M., Lovelace R. V. E., 2008, MNRAS, 386, 1274

Long M., Romanova M. M., Kulkarni A. M., Donati J.-F., 2011, MNRAS, 413, 1061

Long M., Romanova M. M., Lamb F. K., 2012, New Astron., 17, 232

Lovelace R. V. E., Romanova M. M., Bisnovatyi-Kogan G. S., 1995, MNRAS, 275, 244

Manara C. F., Testi L., Natta A., Rosotti G., Benisty M., Ercolano B., Ricci L., 2014, A&A, 568, A18

Matt S., Pudritz R. E., 2005, ApJ, 632, L135

McGinnis P. T. et al., 2015, A&A, 577, A11

Mekkaden M. V., 1998, A&A, 340, 135

Muzerolle J., Calvet N., Briceño C., Hartmann L., Hillenbrand L., 2000, ApJ, 535, L47

Nayak P. K., Narang M., Manoj P., Gorti U., Subramaniam A., George N., Mondal C., 2024, ApJ, 972, 19

Parfrey K., Tchekhovskoy A., 2024, ApJ, 975, 57

Petrov P. P. et al., 2019, MNRAS, 483, 132

Qi C. et al., 2004, ApJ, 616, L11

Rhee J. H. et al., 2007, ApJ, 660, 1556

Ricker G. R. et al., 2015, JATIS, 1, 014003

Robinson C. E., Espaillat C. C., 2019, ApJ, 874, 129

Robinson C. E., Owen J. E., Espaillat C. C., Adams F. C., 2017, ApJ, 838, 100

Robinson C. E., Espaillat C. C., Owen J. E., 2021, ApJ, 908, 16

Robinson C. E., Espaillat C. C., Rodriguez E., 2022, ApJ, 935, 54

Roman-Duval J. et al., 2020, RNAAS, 4, 205

Romanova M. M., Kulkarni A. K., 2009, MNRAS, 398, 1105

Romanova M. M., Lovelace R. V. E. 2006, ApJ, 645, L73

Romanova M. M., Owocki S. P., 2015, Space Sci. Rev., 191, 339

Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E, 2002, ApJ, 578, 420

Romanova M. M., Ustyugova G. V., Koldoba A. V., Wick J. V., Lovelace R. V. E, 2003, ApJ, 595, 1009

Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E, 2004, ApJ, 610, 920

Romanova M. M., Kulkarni A. K., Lovelace R. V. E., 2008, ApJ, 673, L171

Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E, 2009, MNRAS, 399, 1802

Romanova M. M., Long M., Lamb F. K., Kulkarni A. K., Donati J.-F., 2011, MNRAS, 411, 915

Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E, 2012, MNRAS, 421, 63

Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E, 2013, MNRAS, 430, 699

Romanova M. M., Koldoba A. V., Ustyugova G. V., Blinova A. A., Lai D., Lovelace R. V. E, 2021, MNRAS, 506, 372

Rucinski S. M., 1988, Information Bulletin on Variable Stars, No. 3146, #1

Rucinski S. M. et al., 2008, MNRAS, 391, 1913

Setiawan J., Henning Th., Launhardt R., Müller A., Weise P., Kürster M., 2008, Nature, 451, 38

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shu F., Najita J., Ostriker E., Shang H., 1994, ApJ, 429, 781

Sicilia-Aguilar A. et al., 2023, MNRAS, 526, 4885

Singh K. et al., 2024, ApJ, 968, 88

Siwak M., Rucinski S. M., Matthews J. M., Kuschnig R., Guenther D. B., Moffat A. F. J., Sasselov D., Weiss W. W., 2010, MNRAS, 408, 314

Siwak M. et al., 2011, MNRAS, 410, 2725

Siwak M. et al., 2014, MNRAS, 444, 327

Siwak M. et al., 2018, MNRAS, 478, 758

Sokal K. R., Deen C. P., Mace G. N., Lee J.-J., Oh H., Kim H., Kidder B. T., Jaffe D. T., 2018, ApJ, 853, 120

Stauffer J. et al., 2014, AJ, 147, 83

Stauffer J. et al., 2016, AJ, 151, 60

Takami M. et al., 2016, ApJ, 820, 139

Takasao S., Tomida K., Iwasaki K., Suzuki T. K., 2022, ApJ, 941, 73

Teague R., Bae J., Huang J., Bergin E. A., 2019, ApJ, 884, L56

Terquem C., Papaloizou J. C. B., 2000, A&A, 360, 1031

Venuti L. et al., 2019, A&A, 632, A46

Wendeborn J. et al., 2024a, ApJ, 970, 118

Wendeborn J. et al., 2024b, ApJ, 971, 96

Wendeborn J. et al., 2024c, ApJ, 972, 100

Zanni C., Ferreira J., 2013, A&A, 550, 99

Zhu Z., Stone J. M., Calvet N., 2024, MNRAS, 528, 2883

## APPENDIX A: DETAILS OF THE NUMERICAL MODEL

Reference parameters of the mass $M_0$, scale $R_0$, velocity $v_0$, and periods of rotation $P_0$ were described in Section 3.1. Other reference parameters depend on the magnetic field of the star $B_*$ (at the equator) and dimensionless parameter $\mu$, which determines the final size of the magnetosphere. We determine the reference magnetic field $B_0$ and magnetic moment $\mu_0 = B_0 R_0^3$ such that $\mu_0 = \mu_*/\mu$, where the magnetic moment of the star $\mu_* = B_* R_*^3$. Then $B_0 = \mu_0/R_0^3$. We determine the reference density from pressure balance at $R_0$: $\rho_0 = B_0 v_0^3/R_0^3$. The reference value for the accretion rate $\dot{M}_0 = \rho_0 v_0 R_0^2$, for energy flux (used to convert dimensionless energy fluxes shown in plots for light curves to dimensional ones) $\dot{E}_0 = \dot{M}_0 v_0^2$, and for energy flux per unit area (used to show energy flux distribution in hotspots on the surface of the star) $F_0 = \dot{E}_0/R_0^2$. Table A1 shows reference values for different $\mu$ and $B_*$.

3D MHD equations were solved numerically using a Godunov-type numerical scheme, where the seven-wave Roe-type approximate Riemann solver is used. The ‘cubed-sphere’ coordinate system rotates with the star (Koldoba et al. 2002). The energy equation is written in the form of entropy balance, and the equation of state is that of an ideal gas. Viscosity is modelled using the $\alpha$-model
<table><caption>**Table A1.** Reference values for models with different parameters $\mu$ and values of the equatorial magnetic field of the star $B_*$.</caption><thead><tr><th></th><th>Magn. field at the pole</th><th>B<sub>*</sub> (G)</th><th>1000</th><th>800</th><th>600</th></tr><tr><th></th><th>Equatorial field</th><th>B<sub>*</sub> (G)</th><th>500</th><th>400</th><th>300</th></tr></thead><tbody><tr><td rowspan="5">μ = 0.5</td><td>Reference density</td><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>3.24 × 10<sup>-12</sup></td><td>2.08 × 10<sup>-12</sup></td><td>1.17 × 10<sup>-12</sup></td></tr><tr><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>4.15 × 10<sup>-8</sup></td><td>2.66 × 10<sup>-8</sup></td><td>1.50 × 10<sup>-8</sup></td></tr><tr><td>Reference energy flux</td><td>E&#x0307;<sub>0</sub> (erg s<sup>-1</sup>)</td><td>1.50 × 10<sup>33</sup></td><td>9.58 × 10<sup>32</sup></td><td>5.39 × 10<sup>32</sup></td></tr><tr><td>Ref. flux per unit area</td><td>F<sub>0</sub> (erg s<sup>-1</sup>cm<sup>-2</sup>)</td><td>4.37 × 10<sup>10</sup></td><td>2.80 × 10<sup>10</sup></td><td>1.57 × 10<sup>10</sup></td></tr><tr><td>Reference density</td><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>8.11 × 10<sup>-13</sup></td><td>5.19 × 10<sup>-13</sup></td><td>2.92 × 10<sup>-13</sup></td></tr><tr><td rowspan="5">μ = 1.0</td><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>1.03 × 10<sup>-8</sup></td><td>6.65 × 10<sup>-9</sup></td><td>3.74 × 10<sup>-9</sup></td></tr><tr><td>Reference energy flux</td><td>E&#x0307;<sub>0</sub> (erg s<sup>-1</sup>)</td><td>3.74 × 10<sup>32</sup></td><td>2.39 × 10<sup>32</sup></td><td>1.35 × 10<sup>32</sup></td></tr><tr><td>Ref. flux per unit area</td><td>F<sub>0</sub> (erg s<sup>-1</sup>cm<sup>-2</sup>)</td><td>1.09 × 10<sup>10</sup></td><td>7.00 × 10<sup>9</sup></td><td>3.94 × 10<sup>9</sup></td></tr><tr><td>Reference density</td><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>3.60 × 10<sup>-13</sup></td><td>2.31 × 10<sup>-13</sup></td><td>1.30 × 10<sup>-13</sup></td></tr><tr><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>4.62 × 10<sup>-9</sup></td><td>2.96 × 10<sup>-9</sup></td><td>1.66 × 10<sup>-9</sup></td></tr><tr><td rowspan="5">μ = 1.5</td><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>1.66 × 10<sup>32</sup></td><td>1.06 × 10<sup>32</sup></td><td>5.99 × 10<sup>31</sup></td></tr><tr><td>Reference energy flux</td><td>E&#x0307;<sub>0</sub> (erg s<sup>-1</sup>)</td><td>4.86 × 10<sup>9</sup></td><td>3.11 × 10<sup>9</sup></td><td>1.75 × 10<sup>9</sup></td></tr><tr><td>Ref. flux per unit area</td><td>F<sub>0</sub> (erg s<sup>-1</sup>cm<sup>-2</sup>)</td><td>2.03 × 10<sup>-13</sup></td><td>1.30 × 10<sup>-13</sup></td><td>7.30 × 10<sup>-14</sup></td></tr><tr><td>Reference density</td><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>2.60 × 10<sup>-9</sup></td><td>1.66 × 10<sup>-9</sup></td><td>9.36 × 10<sup>-10</sup></td></tr><tr><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>9.35 × 10<sup>31</sup></td><td>5.99 × 10<sup>32</sup></td><td>3.99 × 10<sup>31</sup></td></tr><tr><td rowspan="5">μ = 2.0</td><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>2.70 × 10<sup>9</sup></td><td>1.75 × 10<sup>9</sup></td><td>1.75 × 10<sup>9</sup></td></tr><tr><td>Reference energy flux</td><td>E&#x0307;<sub>0</sub> (erg s<sup>-1</sup>)</td><td></td><td></td><td></td></tr><tr><td>Ref. flux per unit area</td><td>F<sub>0</sub> (erg s<sup>-1</sup>cm<sup>-2</sup>)</td><td></td><td></td><td></td></tr><tr><td>Reference density</td><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td></td><td></td><td></td></tr><tr><td>Reference accretion rate</td><td>M&#x0307;<sub>0</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Figure A1.** Same as in Fig. 7 but for the model $\mu 1.5\theta 5$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure A2. Same as in Fig. 7 but for the model $\mu1.5\theta20$.</figcaption></figure>

(Shakura & Sunyaev 1973), and is incorporated only into the disc so that it controls the accretion rate through the disc. We use a small $\alpha$-parameter $\alpha = 0.02$ in all simulation runs.

At both the inner and outer boundaries, most variables $A_j$ are taken to have free boundary conditions at both the inner and outer boundaries $\partial A_j/\partial r = 0$. At the stellar surface, accreting gas can cross the surface of the star without creating a disturbance in the flow. These conditions neglect the complex physics of interaction between

the accreting gas and the star. The cross-section of the inflowing matter at the distance of stellar radius is interpreted as a spot, and the energy flux distribution in the flow provides us with the energy distribution in hotspots. The magnetic field is frozen on to the surface of the star. That is, the normal component of the field, $B_n$, is fixed, while the other components of the magnetic field vary. At the outer boundary, matter flows freely out of the region. We forbid the back flow of matter from the outer boundary into the simulation region.
<figure><img src="image_21.png" /><figcaption>Figure A3. Same as in Fig. 7 but for the model $\mu=2\theta$.5.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure A4. Same as in Fig. 7 but for the model $\mu=2\theta=20$.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.