© The Authors 2026

1. Introduction

Classical T Tauri stars (CTTSs) are young, low-mass stars that accrete material from their surrounding circumstellar disks. These young star-disk systems are observable across a wide range of wavelengths, from X-rays to the millimeter, and typically exhibit strong surface magnetic fields of the order of kilogauss (kG) (Johns-Krull 2007; Flores et al. 2022). The stellar magnetic field plays an important role in these systems, controlling accretion from the disk and outflow through stellar and disk winds (Alencar et al. 2010; Hartmann et al. 2016).

According to magnetohydrodynamic (MHD) simulations (Romanova et al. 2002; Kurosawa & Romanova 2013), the stellar magnetic field mediates the star-disk interaction, truncating the disk where the magnetic pressure equals the dynamic pressure of the gas. The stellar magnetic field interacts with the inner circumstellar disk, producing accretion columns that channel gas onto the stellar surface. The accreting gas drags the dust accumulated at the base of the accretion funnel, lifting it above the midplane of the disk, and creating a warp at the disk edge (Bodman et al. 2017). When the star-disk system is seen at intermediate to high inclinations, the warp hides the stellar photosphere as the system rotates, causing brightness minima in the light curve of T Tauri stars (McGinnis et al. 2015). Stars that exhibit these characteristics are called dippers (Cody et al. 2014).

The advantage of observing dipper-like stars is that our line of sight passes through all the relevant structures of the system, including the disk wind, the inner disk region, the accretion funnel, and brightness features on the stellar surface. These systems allow us to investigate whether these structures are correlated as predicted by magnetospheric accretion and ejection models.

Bouvier et al. (1999, 2003) studied the dipper-like star AA Tau and originally proposed that its brightness variability might be caused by an inner disk warp at the base of the accretion column. Alencar et al. (2010) analyzed light curves of young stars from the NGC 2264 cluster obtained with the CoRoT satellite and observed that the photometric behavior of AA Tau was not an exception. They found that 28% of the CTTSs in NGC 2264 exhibited dipper-like variability, suggesting this phenomenon is common. In a subsequent observation campaign of NGC 2264 with CoRoT, Cody et al. (2014) significantly expanded the study of young stellar variability in CTTSs. They developed a detailed statistical analysis using metrics of periodicity (Q) and symmetry (M) to classify the observed light curves. These parameters defined a Q–M space suitable for the selection of dipper-like stars. McGinnis et al. (2015) further advanced the understanding of dipper stars by exploring the primary causes of photometric variability in CTTS within NGC 2264. They revealed that dipper-like light curves are predominantly observed in systems with high inner disk inclinations (i > 59°) with respect to our line of sight and that these light curves can transition from periodic to aperiodic and vice versa over a few years. This variability demonstrates that stars can transition between stable and unstable accretion regimes on these timescales (Blinova et al. 2016). More recently, Nagel et al. (2024) extended the inclination limit above which the dipper phenomenon can be observed to 45° by considering the dust survival and sublimation along the accretion funnel. These studies have significantly broadened our understanding of the dipper phenomenon in CTTSs, particularly in solar-mass systems. However, the extent to which the magnetospheric accretion model applies to very-low-mass stars still has some open questions. Studies showed that very-low-mass stars and young brown dwarfs can show signals of accretion (Hartmann et al. 2016; Manara et al. 2023). Jayawardhana et al. (2003) reported that brown dwarfs (i.e., young objects near or below the substellar limit with spectral types between M5 and M8) exhibit detectable accretion rates. Betti et al. (2023) demonstrated a correlation between accretion rate and stellar mass across a broad mass sample, ranging from 2M_⊙ to planetary-mass companions, indicating that accretion still occurs in very-low-mass objects. Despite evidence of accretion in very-low-mass objects, few studies have carried out in-depth analyses of a very-low-mass star using multi-technique, multi-epoch observations, studying the correlations between photometric and spectroscopic variability, the internal structure of the disk, and the role of large-scale magnetic fields in the accretion process. This gap represents one of the main motivations of our work, leading to a more detailed investigation into the applicability of the magnetospheric accretion model for mostly convective stars with masses below 0.5 M_⊙.

To investigate the dipper phenomenon in very-low-mass stars, we focus on JH 223 (Jones & Herbig 1979), also known as EPIC 248006676 or 2MASS J04404950+2551191. JH 223 is a dipper-like star (V = 15.50) that is a member of the Taurus star formation region (Kenyon & Hartmann 1995). Hartmann et al. (1991) classified the source as a weak-line T Tauri star, based on its Hα equivalent width EW_Hα. However, Hartmann et al. (2005) later observed a significant near-infrared excess in its spectral energy distribution, attributed to emission from a circumstellar disk, which suggested that JH 223 was an accreting star. Kraus & Hillenbrand (2007) reported the existence of an M6.5 substellar mass companion located at a separation of 2 arcsec (≃300 au). Herczeg & Hillenbrand (2014) derived a spectral type M2.8 for the brightest member and no detectable veiling in its low-resolution optical spectrum at 751 nm. From ALMA 1.3 mm observations, where both components of the system have been detected, Akeson et al. (2019) derived a primary disk mass of 2 × 10⁻⁴M_⊙. Rebull et al. (2020) reported a rotational period of 3.316 ± 0.003days¹ from a K2 light curve and classified the object as a quasiperiodic dipper. Roggero et al. (2021) further analyzed the K2 light curve, deriving a period of 3.31 ± 0.09days and a dip mean amplitude of 0.32 mag. Our initial motivation for selecting JH 223 was based on its K2 dipper light curve, which exhibits a remarkable stability in period and maximum brightness level over many rotations of the system. Its short rotation period has also enabled monitoring several rotational cycles over the course of a ground-based observational campaign. Additionally, the star was observed again years later by the TESS satellite, which revealed significant photometric variability over the years. Taken together, these datasets provide useful constraints on its variability and accretion processes.

The paper is structured as follows. In Sect. 2, we present the observations conducted on JH 223 over multiple epochs using several telescopes. These datasets formed the basis for the analyses presented in Sect. 3, where we detail our results. In Sect. 4, we discuss our results in the context of the star-disk interaction, as predicted by magnetospheric accretion models and simulations. Finally, in Sect. 5 we present our conclusions.

2. Observations

In this paper, we analyze multi-epoch observations of JH 223 spanning roughly six years, from 2017 to 2023. The data comprise ground-based photometry, photopolarimetry, high-resolution optical spectroscopy, and near-infrared spectropolarimetry obtained during this campaign, as well as complementary space-based photometry from K2 and TESS. A summary of the observations is available in Tables 1–3, while a detailed description of the instruments, observing sites, and data acquisition is given in Appendix A.

3. Results

3.1. Photometry

The K2 light curve of JH 223 is shown in Fig. 1 and exhibits a clear dipper behavior, with a constant maximum flux level interrupted by quasiperiodic luminosity dips (Cody et al. 2014). A periodogram analysis performed on this light curve returns a period of 3.31 ± 0.03 days, where the formal error is measured as the σ of the Gaussian fit of the periodogram peak in the frequency space. We modeled the light curve with a periodic kernel as a Gaussian process that yielded 3.32 ± 0.01 days. Although the period remains consistent within the error bars, we chose to adopt the more conservative uncertainty derived from the periodogram fit, using a period of 3.31 ± 0.03 days in the analyses presented in this study. The light curve folded in phase at this period is also shown in Fig.1.

Fig. 1. JH 223 K2 light curve shown as a function of Julian date (top) and rotational phase (bottom). The rotational phase was computed using a period of 3.31days, with the reference time JD0 = 2 457 817.63 chosen such that phase 0.5 corresponds to the photometric minimum. The color coding indicates to the Julian date of each observation.

Figure 2 shows the V-band light curve of JH 223 obtained during our campaign as a function of the Julian date and folded in phase with the photometric period of 3.31 days, as derived from the K2 observations. The amplitude of the V-band light curve is similar to that of the K2 light curve obtained more than 2 years earlier, with a maximum of m_V = 15.49, and the periodicity has clearly been recovered. The phased light curve exhibits a flat maximum brightness level and relatively deep dips (∼30%), whose depth and shape vary from one cycle to the next, as also seen in the K2 light curve. This type of light curve is characteristic of dippers and quite different from spot modulation, which usually produces smooth, sinusoidal-like shallower light curves that are stable over a timescale of a few weeks. Therefore, we ascribe the photometric variability of JH 223 to the recurrent occultation of the central star by circumstellar dust located close to corotation. A similar photometric behavior and periodicity can be seen in the R_c and I_c light curves of JH 223. The photometric periodicity and phase appear to have been retained over the full range of observations (i.e., over more than a year) from HJD 58530 (CrAO) to HJD 58905 (OACT), a time span corresponding to 113 rotational periods. The color-magnitude diagrams (CMDs) shown in Fig. 3 indicate that the system becomes redder when it gets fainter. Furthermore, the color slope is shallower than the reddening slope expected from ISM-type grains (Cardelli et al. 1989), which indicates that larger grains are responsible for the stellar occultations. For the 2019–2020 observational campaign, which includes contemporaneous ground-based photometry, photopolarimetry, spectroscopy, and spectropolarimetry, we adopted the following ephemeris,

$$\begin{array}{c}\hfill HJD(d)=2458770.2388+3.31E,\end{array}$$(1)

Fig. 2. JH 223 V-band light curve versus HJD (top) and rotational phase (bottom). The phased light curve uses a period of 3.31days and JD0 = 2 458 770.2388 (phase 0.5 at minimum). Arrows mark the epochs of spectroscopic and spectropolarimetric observations. CrAO and OACT data show larger scatter due to higher uncertainties.

Fig. 3. CMDs of (V − R)c vs. V (green) and (V − I)c vs. V (red), with (V − R)c shifted by +1.6mag for clarity. LCOGT (circles) and LT (squares) data are shown. Solid and dashed lines indicate least-squares fits and ISM-extinction slopes, respectively.

where E is the rotational cycle. For photometric datasets separated by more than one year from this campaign, we adopted different JD₀ values to keep the photometric minima around phase 0.5. This approach is required because the uncertainty in the rotational period accumulates over the years, leading to a loss in terms of the phase coherence.

The TESS satellite observed JH 223 during two epochs, in 2021 and 2023. Figures 4 and 5 show the light curves from these two observation periods. The 2021 light curve (sectors 43 and 44) shows smoother and more regular variability than the K2 and ground-based light curves, and it seems to be dominated by spots rather than dust obscuration, as compared to the other epochs of observation.

Fig. 4. Same as Fig. 1, but for the 2021 TESS light curve from sectors 43 and 44. The rotational phase was computed with the 3.31-day period and a reference time of JD0 = 2 459 471.7.

Fig. 5. Same as Fig. 1, but for the 2023 TESS light curves from sectors 70 (a) and 71 (b). The same period and reference time as in Fig. 4 were used.

The 2023 light curve presents two distinct behaviors in sectors 70 and 71, as can be seen in Fig. 5. In sector 70, the variations are complex, and a periodogram analysis suggests a period of 1.66days (i.e., half the period derived at other epochs). The latter is gradually recovered in sector 71, whose light curve presents similar patterns to the K2 light curve, with quasiperiodic brightness dips. The periodogram analysis of the TESS light curves returned period values consistent with the period from the K2 light curve: 3.3 ± 0.1days for the 2021 observations and 3.4 ± 0.2days for sector 71 of the 2023 observations. We still chose to use the period obtained from the K2 light curve in our analysis, since it is more precise, covering a larger number of rotational cycles, as compared to the TESS light curves. The phase-folded TESS light curves with the K2 adopted period are displayed in the bottom panels of Figs. 4 and 5.

3.2. Optical spectroscopy

3.2.1. Stellar parameters

We determined the stellar effective temperature (T_eff), logarithm of surface gravity (log g), projected rotational velocity (v sin i), and radial velocity (υ_rad). We fit the synthetic spectra to the GRACES data, using the spectral synthesis code ZEEMAN (Landstreet 1988; Folsom et al. 2016). This code computes synthetic spectra based on the MARCS atmosphere model (Plez 2008; Gustafsson et al. 2008) and spectral lines information from the Vienna Atomic Line Database (VALD, Ryabchikova et al. 2015). We assumed solar metallicity. The statistical significance of the fit is evaluated from a χ² minimization method based on the Levenberg-Marquardt algorithm (LMA).

We chose the spectrum collected on November 13, 2019 for the following analysis, as it presented the highest signal-to-noise ratio (S/N) among our observations. We ran the fitting procedure independently on nine windows within the 640nm to 870nm wavelength range, excluding the regions with molecular, telluric, and emission lines. Using multiple spectral windows helps to determine the stellar parameters more precisely, as different photospheric lines along the spectrum have different sensitivities to T_eff and log g. We adjusted one parameter at a time, namely, when fitting T_eff, we kept the other parameters fixed. This step was done iteratively for all parameters until we found a converged set of parameters that reproduced the window of interest. Figure 6 shows one spectral window used for stellar parameter determination overplotted with the corresponding synthetic spectra. We computed the average values obtained from the independent windows, excluding spectral windows that presented values larger than 1σ of the mean of each parameter, which provided the following optimal values for the parameters, T_eff = 3528 ± 213K, log g = 3.9 ± 0.1, v sin i = 14.6 ± 0.9km/s and υ_rad = 16.2 ± 0.4km/s. The radial velocity value listed in this section corresponds to the measurement obtained at a specific rotational phase, associated with the spectrum used for the determination of the stellar parameters. The adopted radial velocity of the star was actually derived from measurements taken at various rotational phases and is presented in Sect. 3.3.4. The uncertainty associated with each parameter was determined by evaluating the standard deviation of the values obtained from individual spectral windows. The final parameter values are presented in Table 4. The stellar parameters (T_eff, log g, and v sin i) determined with optical data agree within 1 σ with values obtained by Jönsson et al. (2020) in the near-infrared using APOGEE-2 data.

Fig. 6. Example of a spectral window used for the stellar parameter fitting. The corresponding synthetic spectrum (red line) is overplotted on the observed GRACES spectrum (blue line).

Using the best-fit parameters obtained for each spectral window from the GRACES spectra, we computed theoretical models for the Keck spectra, disregarding two windows that were unavailable as they fell within the inter-order gaps of the spectra. Keeping all parameters fixed except for the radial velocity, which was determined individually for each observation, we found that the resulting fits are consistent with the observed spectra. The radial velocity values for each observation were calculated by averaging the results obtained from each spectral region. The final mean radial velocity values obtained for each observation were 15.38 ± 0.24km/s, 16.17 ± 0.27km/s, and 14.74 ± 0.33km/s, respectively. The uncertainties were assessed based on the variation in χ² around the minimum value identified using the ZEEMAN code for different radial velocity values within each spectral window. The final error bar was determined as the average of the individual uncertainties.

We also computed the luminosity from the relationship between bolometric magnitude and stellar luminosity. We used the maximum magnitude on the V band obtained from the LCOGT observations, m_V = 15.49 (see Sect. 3.1), as the stellar magnitude. We further considered the extinction A_V = 1.2 ± 0.2 from Herczeg & Hillenbrand (2014), the distance d = 143 ± 3pc from Gaia EDR3 (Gaia Collaboration 2021), and we derived the bolometric correction BC_V = −1.91, interpolating for intermediate spectral types from the tables of Pecaut & Mamajek (2013). We then obtained the bolometric magnitude M_bol = 6.6 ± 0.2, resulting in a stellar luminosity L_★ = 0.18 ± 0.03L_⊙. The luminosity we derived agrees with previous values reported in the literature of 0.17 ± 0.01L_⊙ (Herczeg & Hillenbrand 2014) and 0.21 ± 0.03L_⊙ (Roggero et al. 2021).

Taking the values of T_eff and L_★ determined in this paper (Table 4), we obtained the stellar radius R_★ = 1.1 ± 0.2R_⊙ from the Stefan-Boltzmann law. This radius, combined with v sin i and the stellar rotational period, was used to derive the inclination of the star with respect to our line of sight, i = 57° ±14°. This inclination, within the error bars, is consistent both with the inner disk warp model discussed by McGinnis et al. (2015), which requires an inclination greater than 59° for the dipper phenomenon to be observed, and with the results obtained by Nagel et al. (2024), which predict that the dipper phenomenon can be observed in light curves at inclinations larger than 45°, taking into account dust sublimation along the accretion funnel.

Using pre-main sequence (PMS) evolutionary tracks from Baraffe et al. (2015), we derive the stellar mass of 0.4 ± 0.2M_⊙ and an age of ∼3 Myr from the location of JH 223 in the Hertzsprung-Russell diagram (HRD; see Fig. 14). The evolutionary model also indicates that JH 223 is fully convective. Table 4 summarizes the main stellar parameters of JH 223.

3.2.2. Veiling

The spectra of accreting T Tauri stars exhibit veiled photospheric lines compared to non-accreting systems. From the ultraviolet to the optical, veiling originates from the continuum emission of the accretion shock at the stellar surface (Hartmann et al. 2016).

To measure the optical veiling of JH 223, we selected four regions in the GRACES spectra between 534.5nm and 758.0nm. These regions were chosen for their high S/N and significant number of photospheric lines. We used the ESPaDOnS spectra of the weak-line T Tauri star TWA 7 as a template. TWA 7 has an age of 17 Myr, a mass of 0.62 ± 0.03M_⊙ and v sin i = 4.5 ± 0.2km/s (Nicholson et al. 2021). The spectral type of this star varies in the literature: according to Webb et al. (1999), TWA 7 has a spectral type M1, while Manara et al. (2013) quote an M2 spectral type and Herczeg & Hillenbrand (2014) list an M3.2 spectral type. We broadened the normalized spectrum of TWA 7 according to the projected rotational velocity of JH 223 (v sin i = 14.6 ± 0.9km/s) and corrected the shift due to the radial velocity difference between the two stars. Finally, we calculated the veiling (r), defined as the ratio of the continuum excess flux to the stellar flux (r = F_Excess/F_Phot), through a χ² minimization procedure to match the spectra of JH 223 with TWA 7 (Sousa et al. 2023).

Table 5 shows the mean veiling values of all the spectral regions between 534.5nm to 758.0nm, computed for each night of observation. The error was determined as the standard deviation of the veiling measurements obtained from the selected spectral regions. We derived the average veiling of 0.12 ± 0.05, where the error was calculated using the standard deviation of all measurements. We conclude that the star presents low optical veiling, which is constant in our observations within the errors, indicating a low mass accretion rate (see Sect. 3.2.4). Our low optical veiling values are consistent with those derived using the ZEEMAN code in Section 3.2.1, both indicating low veiling in the 600–750nm range. These results are also in agreement with those reported by Herczeg & Hillenbrand (2014), who found no veiling at optical wavelengths.

3.2.3. Circumstellar lines

We analyzed the Hα, Hβ, and He I 587.6nm lines in the GRACES and Keck spectra of JH 223. The Hα and Hβ lines are formed along the accretion funnel, while He I 587.6nm originates near the accretion shock, making them good tracers of accretion in the optical domain. We removed the photospheric contribution from the JH 223 spectra to obtain only the circumstellar component by subtracting from the JH 223 spectra a veiled and rotationally broadened mean spectrum of TWA 7 (see Sect. 3.3.1). Figure 7 shows the circumstellar profiles of the Hα, Hβ and He I lines. In Appendix B, we illustrate the emission lines before and after the removal of the photospheric contribution for each night of observation.

Fig. 7. Hα (top), Hβ (middle), and He I 587.6 nm (bottom) emission lines used to derive mass accretion rates. Labels indicate the rotational phase. GRACES and Keck data are shown as solid and dashed lines.

The Hα and Hβ lines exhibit a redshifted absorption component. This absorption reaches its maximum depth at phase 0.477, indicating that the accretion column crosses our line of sight at this phase. Furthermore, the maximum emission of the He I line at the same phase suggests that the accretion shock region also intersects our line of sight along with the accretion funnel. These results indicate that the accretion column and the accretion shock are spatially correlated with the inner disk warp, which obscures the star around phase 0.5, as observed in the light curves of JH 223.

3.2.4. Mass accretion rate

We estimated the mass accretion rate (Ṁ_acc) from the flux of the Hα and Hβ lines, which are considered to be good accretion tracers in optical spectra (Gullbring et al. 1998; Alcalá et al. 2017). Our spectra were not flux calibrated, so we calculated the line flux (F_line) as the product of the line equivalent width (EW_line) and the stellar continuum flux at the rotational phase of the spectral observation,

$$\begin{array}{c}\hfill F_{\mathrm{line}}=E{W}_{\mathrm{line}}{F}_0\times {10}^{-0.4(m-A)},\end{array}$$(2)

where F₀, m, and A are the reference flux of a zero magnitude star, the magnitude, and the extinction of JH 223 in the corresponding photometric band. We calculated EW_line by integrating the spectrum around the central wavelength of Hα and Hβ. The results obtained are shown in Table B.1.

For the Hα line, we considered the R band values: F₀ = 2.177 × 10⁻⁹ergs⁻¹cm⁻²Å (Bessell et al. 1998), A_R = 0.9 ± 0.2 (derived from the extinction relation A_R/A_V = 0.75, Cardelli et al. 1989), choosing the R band magnitude from LCOGT closest in rotational phase to the spectral observation. For this line we obtained EW_Hα ranging from 2.66Åto 4.56Å, resulting in line fluxes going from 3.06 × 10⁻¹⁴ergs⁻¹cm⁻² to 6.26 × 10⁻¹⁴ergs⁻¹cm⁻².

For the Hβ line, we used the V band values: F₀ = 3.631 × 10⁻⁹ergs⁻¹cm⁻²Å (Bessell et al. 1998), A_V = 1.2 ± 0.2 (Herczeg & Hillenbrand 2014), and chose the V band magnitude from LCOGT closest in rotational phase to the spectral observation. Knowing that the EW_Hβ ranges from 1.42Åto 2.64Å, we obtained line fluxes going from 0.91 × 10⁻¹⁴ergs⁻¹cm⁻² to 1.81 × 10⁻¹⁴ergs⁻¹cm⁻².

We derived the line luminosity from L_line = 4πd²F_line, using the stellar distance of d = 143 ± 3pc, and calculated the accretion luminosity (L_acc) from the empirical relations between L_line and L_acc provided by Alcalá et al. (2017). The accretion luminosity can be converted to mass accretion rate (Gullbring et al. 1998) via

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}={(1-\frac{R_{★}}{r_m})}^{-1}\frac{L_{\mathrm{acc}}{R}_{★}}{G{M}_{★}},\end{array}$$(3)

where R_★ and r_m are the stellar radius and the disk truncation radius, respectively. We used r_m = 5R_★, which is a typical value used in previous works (Gullbring et al. 1998; Bouvier et al. 2007b; Alcalá et al. 2017; Pouilly et al. 2020). Table B.1 shows the results of the mass accretion rate calculations for Hα and Hβ. We obtained a mean mass accretion rate of (5 ± 1)×10⁻¹¹M_⊙yr⁻¹ from Hα and (8 ± 2)×10⁻¹¹M_⊙yr⁻¹ from Hβ. The average mass accretion rate, computed using these values, was (7 ± 1)×10⁻¹¹ M_⊙yr⁻¹. Hartmann et al. (2016) and Betti et al. (2023) showed that there is a linear trend between mass accretion rate and stellar mass for young stars. When comparing their results and our mass accretion rate values, JH 223 presents a mass accretion rate lower than most young stars of similar mass.

We also calculated the mass accretion rates using the equivalent widths derived from the color index in the narrowband filters, CI_Hα, as described in Appendix A.2. The average mass accretion rate we obtained was Ṁ_acc = (5.5 ± 0.8) × 10⁻¹¹M_⊙yr⁻¹. These results align with the values calculated from the GRACES and Keck spectra, showing similar variability amplitude and average value. Furthermore, we observed a rotational modulation in the accretion luminosity, with a minimum occurring around phase 0.5 (see Fig. B.3). The modulation observed in L_acc indicates that the line emission region is obscured by the inner disk warp around phase 0.5. When the accretion luminosity is converted to Ṁ_acc, it would imply that the accretion rate is lower exactly when the main accretion stream, which creates the inner disk warp (see Sect. 4), is expected to be facing the observer. Therefore, we emphasize that a decrease in L_acc does not necessarily correspond to a real decrease in Ṁ_acc, but it instead reflects the obscuration of the accretion flow at this phase.

3.3. Near-infrared spectropolarimetry

3.3.1. Veiling

The nIR veiling may have a distinct origin compared to the optical veiling. It is caused by excess emission from the circumstellar disk, resulting from dust heating in its inner regions. This heating is associated with reprocessing stellar radiation and releasing energy during the accretion process (Chiang & Goldreich 1997; Sousa et al. 2023).

To measure the nIR veiling of JH 223, we applied the same procedures as described in Sect. 3.2.2, using the SPIRou spectra of JH 223. We selected the spectral region between 1158nm and 1167nm, in the Y band, to measure the veiling as it is close to the He I 1083nm line, which is one of the lines of interest in our analysis. As a template, we used the SPIRou spectra of the weak-line T Tauri star TWA 7.

Table 6 shows the veiling computed for each night of observation. We derived the average veiling of 0.10 ± 0.07, where the error was calculated using the standard deviation of all measurements. We did not detect any rotational modulation in the veiling variability. We conclude that our low veiling values are in agreement with the optical veiling results obtained in Sect. 3.2.2. We also measured the veiling in the spectral region from 2260nm to 2269nm, in the K band. The average veiling value of 0.28 ± 0.07 was larger than in the Y band, a behavior consistent with the results reported by Sousa et al. (2023), and the veiling values showed, once more, no rotational modulation.

3.3.2. Emission lines

We analyzed the He I 1083nm and Paβ lines present in the SPIRou spectra of JH 223. The He I line may originate in the accretion funnel, the accretion shock, and winds (Erkal et al. 2022), while Paβ is thought to be formed mostly along the accretion funnel. Therefore, these lines are good tracers in the infrared of accretion and mass loss (Sousa et al. 2021). We did not detect any other circumstellar line exhibiting signs of accretion.

We removed the photospheric contributions from the JH 223 spectra in the He I and Paβ line regions to obtain only the circumstellar contribution. Thus, we subtracted a veiled and rotationally broadened mean spectrum of TWA 7 from the JH 223 spectra.

The circumstellar Paβ line did not show any relevant signal, indicating that the absorption lines in the observed spectra were all photospheric (see Fig. B.2). The circumstellar He I line presents blueshifted and redshifted absorption components, as shown in Fig. 8. The blueshifted absorption component, around −40kms⁻¹, is commonly attributed to outflows projected on our line of sight, such as stellar and disk winds, which absorb photons emitted by the star and the accretion funnel (Edwards et al. 2006; Erkal et al. 2022). This component persists in all observations and shows small variability in depth and width. The redshifted absorption component traces gas in the accretion funnel, moving away from the observer, which absorbs photons from the shock region (Fischer et al. 2008). In the JH 223 spectra, this component exhibits large variations in velocity and depth. It becomes stronger at phases between 0.4–0.6, which correspond to the photometric minima and to the passage of the accretion funnel in our line of sight, than around phases 0.0–0.3, where it is hardly seen (see Fig. B.2). This indicates a spatial association between the funnel flow and the inner disk warp.

Fig. 8. Residual line profiles of HeI 1083nm. The colors represent different rotational phases.

3.3.3. Least-squares deconvolution

We analyzed the SPIRou spectra by applying the least-squares deconvolution (LSD) method (Donati et al. 1997) to atomic photospheric lines to derive weighted average Stokes I and V profiles for each night of observation. This method increases the S/N of the Stokes parameters (Donati et al. 1997). In particular, it is well suited for analyzing the Stokes V profiles, which usually present very low S/N values in the individual photospheric lines of T Tauri stars (Donati et al. 1997; Kochukhov et al. 2010).

Overall, LSD assumes that the intensity (polarization) spectrum results from a convolution between the weighted mean Stokes I (Stokes V) LSD profile and a mask of photospheric atomic lines. We used the VALD database and MARCS atmosphere models to generate the line mask, considering an effective temperature of T_eff = 3500K, surface gravity of log g = 4.0, and microturbulence velocity of v_mic = 2.0km/s. These values are compatible with the stellar parameters of JH 223 (see Sect. 3.2.1). We further selected only atomic lines with known Landé factors and a depth larger than 3% of the continuum level, resulting in a final mask with 2307 photospheric lines.

The LSD profiles were computed using the LSDpy code with a window of 200km/s, adopting for the normalization the average Landé factor of 1.2, a central wavelength of 1700nm, and a line depth of 0.15. In the next sections, we use the photospheric Stokes I and V LSD profiles, shown in Fig. C.2, to investigate the radial velocity and longitudinal magnetic field variability, and perform a Zeeman-Doppler imaging analysis of the spectropolarimetric data of JH 223.

3.3.4. Radial velocity

Radial velocity values were calculated as the first-order moment of the Stokes I LSD profiles, as in Donati et al. (2017), in the interval of ±30km/s with respect to the profile center. The radial velocities calculated for each night of observation are listed in Table 6 and shown in Fig. 9, as a function of the rotational phase using Eq. (1). We obtained radial velocities ranging from 14.2km/s to 16.2km/s, with an average value of 15.3 ± 0.6km/s. The radial velocity variation corresponds to distortions of the Stokes I LSD profiles due to the presence of spots at the stellar surface, while the mean radial velocity value is attributed to the true radial velocity of JH 223. We calculated the error of individual radial velocities propagating the uncertainties in the first-order moment equation of Stokes I. We obtained the error of the average value from the standard deviation of the individual measurements. The radial velocities obtained with the GRACES and Keck spectra (see Sect. 3.2.1) agree within 1σ with the individual radial velocity measurements from SPIRou (see Fig. 9).

Fig. 9. JH 223 radial velocities (top) and residuals (bottom). Radial velocities were derived from the first-order moment of the Stokes I LSD profiles of the SPIRou observations (black dots) and from spectral fitting of GRACES (red dot) and Keck spectra (green dots). The blue curve represents the sinusoidal fit with the rotational period fixed at 3.31days. Residuals (squares), computed as the difference between the RV values and the sinusoidal fit, have a standard deviation of 0.27km/s.

To verify whether the modulation of the radial velocity caused by spots on the stellar surface is associated with the period derived from the K2 light curve, we phase-folded the RV measurements using the ephemeris in Eq. (1). We then fit a sinusoidal curve to the phase-folded data, fixing the period to that of the light curve and optimizing the amplitude and phase offset parameters of the sine function. The fit showed good agreement between the RV modulation and the K2 light curve period, suggesting the stellar rotation period is consistent with the rotational period of the warp in the inner disk. This result indicates that the disk deformation responsible for the dips in the K2 light curve is located close to the disk’s co-rotation radius. Furthermore, the negligible phase offset of the sine function (0.02 ± 0.05) suggests that the main spot on the stellar surface crosses our line of sight around phase 0.5 (see Sect. 3.3.6 and Fig. 11), which is in agreement with the warp location. It is important to bear in mind that we cannot discard the possibility that the warp itself induces part of the RV variations, as it modulates the visible area of the stellar surface, much as a cool spot would do.

3.3.5. Longitudinal magnetic field

We estimated the longitudinal magnetic field component (B_l) from the first-order moment of the Stokes V LSD profiles (Kochukhov 2021) via

$$\begin{array}{c}\hfill B_{\ell }=-\frac{7.145\times {10}^5}{{\lambda }_0{g}_{\mathrm{eff}}}\frac{\int [v-{\upsilon }_{\mathrm{rad}}]V(v)\mathrm{d}v}{\int [1-I(v)]\mathrm{d}v},\end{array}$$(4)

where B_ℓ is expressed in G, λ₀ is the average wavelength in nm, and v and υ_rad denote the velocity at each point along the Stokes V profile and the radial velocity of the star, respectively, both in km/s. The average wavelength, λ₀, and the effective Landé factor, g_eff, are the same ones used in Sect. 3.3.3 to compute the LSD profile. I(v) and V(v) correspond to the Stokes I and V LSD profiles, respectively. We computed Equation (4) in the ±30km/s interval with respect to the center of the Stokes I LSD profile. We obtained a longitudinal magnetic field ranging from 0G to −190G with an average field value of −111 ± 56G. The values obtained for each observation are presented in Table 6 and shown as a function of the rotational phase in Fig. 10. We fitted a sine function to the B_l values, keeping the period of 3.31days fixed, and found a modulation in agreement with the rotational period of JH 223.

Fig. 10. Same as Fig. 9, but for the longitudinal magnetic field. Residuals have a standard deviation of 33 G.

3.3.6. Zeeman-Doppler imaging

We used the tomographic technique of Zeeman-Doppler imaging (ZDI, Semel 1989) to obtain information on brightness and magnetic field features at the surface of JH 223. This technique consists of tracing rotationally modulated signatures present in the time-series Stokes I and V LSD profiles to reconstruct the brightness distribution and large-scale magnetic field topology at the stellar surface (Donati et al. 1997, 2006). For this purpose, we employ the ZDIpy code written by Folsom et al. (2018), which has been largely used in the literature to model active stars (Cang et al. 2020; Pouilly et al. 2021; Bellotti et al. 2023).

The ZDI model treats the star as a spherical surface divided into a grid with elements of approximately equal area. We configured the surface of our model with 60 longitudinal rings, which correspond to approximately 2000 surface elements. The poles of the spherical coordinate system are aligned with the star’s rotation axis. The algorithm computes synthetic local Stokes I profiles for each surface element by employing a Voigt profile, weighted by its projected area, limb darkening, and local brightness. We derived the Voigt profile parameters (specifically, the Gaussian and Lorentzian widths) through a fitting procedure to the Stokes I LSD profile of the weak-line T Tauri star TWA7 (see Sect. 3.2.2), a slow rotator of the same spectral type as JH 223. The resultant fit provided Gaussian and Lorentzian widths of 9.7 km/s and 0.06 km/s, respectively. We retrieved the linear coefficient of limb darkening from the tables of Claret & Bloemen (2011), using the value of 0.2865 for TWA 7 and 0.2792 for JH 223. The local Stokes V profiles are derived using a weak field approximation, expressed as a function of the derivative of the Stokes I profile (Folsom et al. 2016), and the magnetic field is represented as a combination of spherical harmonics for the radial, azimuthal, and meridional components. The full disk-integrated line profiles, computed as the summation of local profiles at the visible surface Doppler-shifted by the line-of-sight projected velocity of the surface element, are used to reconstruct the surface brightness or the magnetic field at the stellar surface. The brightness map starts with a homogeneous brightness distribution, and the code iteratively adds dark and bright spots to the stellar surface to identify the optimal brightness distribution that reproduces the observed Stokes I profiles. To optimize the modeling of the large-scale magnetic field, the code adjusts the spherical harmonic coefficients of the radial, meridional, and azimuthal field components. This iterative process applies the principles of maximum entropy (Skilling & Bryan 1984) and χ² minimization as regularization criteria to the image reconstruction.

In our preliminary analysis, the brightness and magnetic field maps converged rapidly to ${\chi }_r^2\sim 1$, resulting in models with little information. Wade et al. (2000a,b) argued that the error bars of the LSD profiles are typically overestimated, bringing the need to aim for a target ${\chi }_r^2$ lower than 1 to compute the amount of information suitable for brightness and magnetic field maps. This highlights the need for a robust method to select the target ${\chi }_r^2$ that allows the ZDI model to extract information without overfitting the data. We proposed to apply a scaling factor to the error bars of the Stokes I and V profiles, in a similar way to that done by Wade et al. (2000b), following

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{new}}={\sigma }_{\mathrm{LSD}}\sqrt{{\chi }_{r,sf}^2},\end{array}$$(5)

where σ_LSD is the LSD error and ${\chi }_{r,sf}^2$ is the value that extracts information from the data without overfitting. As detailed in Appendix C.1, we obtained ${\chi }_{r,sf}^2=0.53$ for our dataset. We applied the correction to the error bars of all LSD profiles, which allowed us to target a unit ${\chi }_r^2$ in the following image reconstructions.

After applying the procedure described in the previous sections, we modeled the brightness and large-scale magnetic field topology on the surface of JH 223. Initially, we reconstructed the brightness distribution fitting the observed Stokes I LSD profiles (Fig. C.2, left). The reconstructed brightness map is displayed in Fig.11 with a spot and plage coverage of 18% of the stellar surface. The map shows a high contrast cool spot near the pole, located at about 60° latitude with respect to the equator, toward phases 0.40 − 0.45. Additionally, there are lower contrast features near the equator.

Fig. 11. Logarithmic brightness map of the surface of JH 223 in November, 2019. The star is shown in a flattened polar view down to a latitude of −30°, with the equator indicated by a solid line and latitudes of 60° and 30° by dashed lines. Outer radial ticks mark the phases of spectropolarimetric observations. Cool spots are shown in brown shades, while bright plages are shown in blue shades.

We then modeled the time series Stokes V profiles (Fig. C.2, right) to obtain the large-scale magnetic field topology, imposing the brightness distribution previously calculated. We truncated the spherical harmonics expansion describing the magnetic field at the order of ℓ_max = 5. Including terms with a higher order did not significantly change the reconstructed magnetic field, in line with what is expected for stars with low v sin i (Morin et al. 2010). The resulting maps for radial, azimuthal, and meridional components of the magnetic field are shown in Fig. 12. The large-scale magnetic field topology is predominantly poloidal, with 71% of the total magnetic energy, while the toroidal component represents 29% of the total magnetic energy. The poloidal field has 28% of axisymmetry, and the dipolar component is the most significant component, with 37% of the poloidal energy, tilted by 28° toward phase 0.38. The quadrupolar and octupolar components represent 6% and 14% of the poloidal energy, respectively. The average unsigned magnetic field strength across the stellar surface is 358G, with the maximum field intensity reaching 915G. The magnetic properties of JH 223 are summarized in Table 7.

Fig. 12. ZDI maps of the radial (top), meridional (middle), and azimuthal (bottom) components of the large-scale magnetic field at the surface of JH 223. Similar to Fig. 11, the star is represented in a flattened polar projection. Magnetic fluxes, indicated by the color bar, are expressed in gauss.

4. Discussion

JH 223 is a classical T Tauri star that exhibited typical dipper-like behavior when it was observed by K2 in 2017. The K2 light curve shows a constant maximum flux interrupted by minima that vary in width and depth from cycle to cycle. The Q and M parameters, calculated following the methodology of Cody et al. (2014) and Cody & Hillenbrand (2018), further support the classification of this light curve as a typical dipper, indicating a clear quasiperiodic variability (Q = 0.37) with positive asymmetry (M = 0.76). This quasiperiodic variability is similar to the photometric variability of known dipper stars, such as AA Tau (Bouvier et al. 2007a) and LkCa 15 (Alencar et al. 2018), and indicates the presence of circumstellar dust in the inner disk region that occults the star as the system rotates. These dust occultations are thought to come from an inner disk warp produced by the interaction of an inclined stellar dipole magnetic field component and the inner circumstellar disk, as seen in MHD simulations (Romanova et al. 2013). Our observing photometric, polarimetric, and spectroscopic campaign was carried out at the end of 2019, with the ground-based light curves obtained displaying the same characteristics as those seen in the prior K2 data. The CMDs show that the system becomes redder when fainter, and the optical linear polarimetric observations present an increase in the degree of polarization in the photometric minima, which gives support for a dust obscuration origin of the photometric variations.

JH 223 was later observed with the TESS satellite in 2021 and 2023. The 2021 TESS light curve (sectors 43 and 44; see Fig. 4) shows a much more regular variability pattern, reminiscent of spot-like light curves (see Cody & Hillenbrand 2018), but with a scalloped morphology, as though some dust still obscured the star at that epoch. Although the flux asymmetry for these sectors still corresponds to a dipper-like light curve with M = 0.5, the periodicity parameter Q = 0.04 suggests that a periodic physical process, such as spots, dominates the light curve. The influence of spots on the light curve variability is significant in 2021, imposing a strict periodicity that is typically not observed in dipper-like stars. This finding indicates that the presence of cool spots on the stellar surface, modeled with ZDI in 2019, is recurrent in JH 223. Additionally, there is a decrease of about 20% in the amplitude of the flux variability in the 2021 TESS light curve when compared to the other epochs, suggesting a lower amount of circumstellar dust at the warp location.

The 2023 TESS light curve variability seemed again to be due mostly to obscuration by dust (sectors 70 and 71, see Fig. 5). However, sector 70 presented an asymmetric (M = 0.56) and aperiodic (Q = 0.65) light curve, as though the circumstellar dust were randomly distributed in the inner disk region, as predicted by unstable+ordered and solely unstable MHD models (Blinova et al. 2016). Remarkably, the phased light curve of Sector 70 suggests two main obscuration events occurring at intermediate phases of ∼0.3 and ∼0.8; namely, in quadrature, compared to the main dip seen at phase 0.5 at other epochs. After about 25 days of observations, the photometric light curve in sector 71 gradually recovered the quasiperiodicity (Q = 0.30) and asymmetry (M = 0.73) that are typical of dipper-like stars, and looked again very similar to the K2 and our ground-based data. The morphological light curve variation in several timescales over the years seems to indicate that the dust distribution in the inner disk changes considerably in timescales of weeks to years. Blinova et al. (2016) MHD simulations showed the development of a pair of opposite accretion tongues in the ordered unstable regime, which, at times, seems to be perpendicular to the longitude of the magnetic pole (see their Fig.9). We might thus have witnessed the system transitioning from an ordered unstable regime during TESS Sector 70 (with two opposite tongues nearly perpendicular to the magnetic pole, producing the dips at phases 0.3 and 0.8) and going back to a dominant funnel flow in Sector 71, along with the associated disk warp responsible for the main dip at phase 0.5. Figure 13 shows the variability of the Q and M parameters for the light curves of JH 223 over six years of observations, overplotted with the parameters obtained for others CTTSs with confirmed disks in the Taurus star-forming region (Cody et al. 2022). Similar variations in light curve morphology of CTTSs on a timescale of 3 years were already pointed out by McGinnis et al. (2015) when analyzing CoRoT satellite light curves of 84 CTTSs in NGC 2264. They showed that about half of the dipper-like light curves became aperiodic and vice-versa on this timescale, hinting at the highly dynamical nature of the star disk interaction that shapes the dust distribution in the inner disk of accreting systems. This type of photometric behavior may actually be typical of low-mass accretion rate systems, as it was also reported for PDS 70 (Gaidos et al. 2024), a K7 CTTS with a transition disk and two orbiting exoplanets, which was observed with TESS on three different epochs over four years and showed light curves that varied from spot to dust dominated. PDS 70 also presented Hα and He I 1083 nm line profiles similar to those observed in JH 223 at certain epochs (Thanathibodee et al. 2020; Gaidos et al. 2024), showing that accretion and outflow are ongoing with about the same characteristics in both systems. As JH 223 is younger than PDS 70 and lacks detected planets, we may be witnessing the early stages of planet formation, potentially triggered by disk instabilities that could explain the observed dust density variability.

Fig. 13. Asymmetry (M) and quasiperiodicity (Q) parameters derived from the K2 and TESS light curves of JH 223. The gray points correspond to the Q and M parameters of others CTTSs analyzed by Cody et al. (2022) in the Taurus star-forming region.

Despite the morphological changes in the light curve, the periodicity of the photometric variability of the JH 223 TESS observations was maintained most of the time in the almost seven years of photometric data analyzed in this work. It is important to keep in mind that the period measured from a dipper light curve actually corresponds to the period of the disk radius where the dust overdensity is located, and may not be related to the stellar rotation period. However, in the case of JH 223 (as with many other dipper CTTSs), the inner disk warp at r_m = 6 ± 1R_★, calculated with Equation (6) of Bessolaz et al. (2008), occurs at the corotation radius between the star and the disk, r_cor = 6 ± 1R_★. The stellar rotation period (and, consequently, the location of the dust at the corotation radius) can also be verified through the radial velocity and longitudinal magnetic field variations, which come from cold spots at the stellar photosphere and the surface magnetic field, which also vary at the same period as the K2 light curve. The fact that the 2021 TESS light curve seems to be dominated by spots rather than dust variability, still showing the same periodicity as before, also points to the location of the circumstellar dust warp, when present, near the corotation radius. Additionally, the dust sublimation radius, calculated from Eq. 1 of Monnier & Millan-Gabet (2002) with T_sub = 1500K and Q_R = 1, lies at a distance of r_sub = 2.8 ± 0.3R_★, reinforcing that, at the disk warp location (6 ± 1R_★), the disk is sufficiently cold for dust to exist.

JH 223 presents a low mass accretion rate, with consistent measurements in the optical and in the nIR. The Hα and Hβ emission lines, as well as the He I 1083 nm line, show redshifted absorption around phase 0.5, indicating that the main accretion funnel is in our line of sight at about this phase. Furthermore, the maximum emission of the He I 587.6 nm line at the same phase suggests that the accretion shock region also intersects our line of sight along with the accretion funnel. These results coincide with the observations of the quasiperiodic luminosity dips in the contemporaneous LCOGT light curve and show the spatial connection between the accretion shock, the accretion column, and the dust warp in the inner disk.

The reconstructed large-scale magnetic field is complex and presents an inclined dipole component, which creates a major accretion column in each hemisphere, according to MHD simulations of star-disk interaction (Kurosawa & Romanova 2013). The dipole, with a maximum inclination around phase 0.4, is thought to give rise to the warp of gas and dust at the base of the accretion funnel near the disk, which periodically occults the stellar photosphere, with a maximum obscuration near phases 0.4–0.5, as observed. The ∼0.1 phase lag between the longitude of the magnetic pole and the maximum obscuration from the warp is predicted by analytical models of misaligned magnetic star-disk interaction to result from the delayed response of the inner disk to the magnetic perturbation (Terquem & Papaloizou 2000). These models also predict a trailing warp, which may account for the asymmetric shape of the main dip at the time of our observations, with a sharper ingress and a shallower egress.

According to our analysis, the characteristic dimming patterns observed in the light curves of JH 223 can be explained by the obscuration of the stellar surface by an inner disk warp, which originates from the interaction of an inclined dipole and the inner disk region. The disk warp is spatially connected to the accretion funnel, as evidenced by the appearance of redshifted absorption components in the Hα and He I 1083 nm lines, generated along the accretion funnel, which coincide in phase with the warp obscuration of the stellar surface. Other scenarios have been proposed in the literature to explain the dipper light curves, such as vertical hydrodynamic disk instabilities or dusty disk winds (Ansdell et al. 2020). Hydrodynamic instabilities can form vortices in the inner part of the disk, lifting dust above the disk midplane. These vortices are generated in regions of low ionization and weak coupling between the magnetic field and the gas (dead zones) (Lyra & Mac Low 2012; Flock et al. 2017), and for stars with an inclination between 70° and 80°, this mechanism can be a potential source of overdensity formation that periodically obscures the stellar photosphere (Ansdell et al. 2020). In the case of JH 223, the dust overdensity that obscures the star is located at the corotation radius of the disk, which coincides with the disk truncation radius, indicating that the deformation periodically observed is physically associated with the accretion column and, consequently, with the magnetic field, and therefore is not located in a dead zone region. The dusty disk wind scenario can also be excluded since the disk wind in JH 223, detected through a narrow blue-shifted absorption component in the He I line 1083 nm, does not vary with the stellar rotational period, remaining almost constant in our observations, which indicates that the winds mostly originate from a region far from the corotation radius.

Figure 14 shows JH 223 in an HRD that displays the large-scale magnetic characteristics of young low-mass stars, reconstructed with ZDI. Over the last decade, the large-scale magnetic properties of solar-mass stars were studied by several authors across the PMS. As can be seen in Fig. 14, their magnetic fields evolve from strong, axisymmetric, and mostly dipolar when the stars are fully convective, to fainter and complex fields as their internal structure changes, a radiative core develops, and they go toward the MS. However, JH 223 falls in a region of the HRD containing very few subsolar mass objects, which we do not have much information on with respect to the surface magnetic field structure. Although JH 223 is fully convective, as shown by V2247 Oph and V347 Aur, the only other PMS stars with known magnetic field topologies and masses less than 0.5 M_⊙, the three systems already show complex and rather faint surface magnetic fields at an early age, in contrast to solar-type stars. On the MS, very-low-mass, fully convective, M dwarfs present a dynamo bistability, with some stars showing complex and weak magnetic fields, while others present strong dipolar fields (Morin et al. 2011). It would be interesting to populate the same M dwarf region in the PMS to see if the bistability phenomenon holds for stars with the same internal structure but different ages and rotational properties.

Fig. 14. HRDs showing the large-scale magnetic properties of JH 223 and other classical and weak-line T Tauri stars. CTTSs are labeled by name. Symbol size scales with the average magnetic field strength, color with the poloidal magnetic energy fraction, and shape with the fraction of poloidal energy in axisymmetric modes. Evolutionary tracks (black lines) and isochrones (green dotted lines) from Baraffe et al. (2015) are shown, together with the fully convective boundary and the 40% radiative-core mass limit (dashed blue lines). Magnetic field data for the T Tauri stars are taken from Donati et al. (2014, 2015), Yu et al. (2017), Hill et al. (2017), Nicholson et al. (2018), Hill et al. (2019), Yu et al. (2019), Nicholson et al. (2021), Xiang et al. (2023), Donati et al. (2007, 2008, 2010a,b, 2011a,b, 2012, 2013, 2019, 2020a), Zaire et al. (2024), Donati et al. (2026).

5. Conclusion

We analyzed photometric, spectroscopic, and spectropolarimetric data of the young low-mass accreting star JH 223, aiming to investigate the stellar magnetic field topology and the star-disk interaction. The large-scale surface magnetic field is predominantly poloidal, with a dipolar component strong enough to truncate the accretion disk near the corotation radius. The observed dipper light curves of JH 223 are associated with the presence of an inner disk warp at the base of the accretion columns, obscuring the stellar surface periodically as the star rotates. The same periodicity is seen in the modulations of radial velocity and the longitudinal magnetic field. We also associated the accretion column, traced by the maximum redshifted absorption in the Hα and He I 1083 nm circumstellar lines, with the inner disk warp, as they both occur at similar rotational phases. The accretion process in JH 223 is dynamic, changing from an unstable to a stable regime on a timescale of a few weeks, as predicted by MHD models of the star-disk interaction. This study provides evidence that very-low-mass PMS stars with low mass accretion rates remain subject to magnetospheric accretion. Finally, our findings highlight the fact that classifications of individual stars based on photometric variability (e.g., dipper, quasiperiodic, or chaotic) only capture a transient state. Light curves represent a snapshot of the star’s condition at a given epoch, which can evolve rapidly, reinforcing the importance of multi-epoch monitoring to understand their underlying physical processes.

Acknowledgments

TPF, BZ and SHPA acknowledge financial support from CNPq, CAPES, and FAPEMIG. JB acknowledges funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 742095; SPIDI: Star-Planets-Inner Disk Interactions, https://www.spidi-eu.org). AB acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. AF and JAS acknowledge financial support from Large Grant INAF-2024 “Spectral Key features of Young stellar objects: Wind-Accretion LinKs Explored in the infraRed (SKYWALKER)”. We thank Dmitriy Blinov and Gina Panopoulou for providing the RoboPol data and for the discussions that enriched this work. This work is based on observations obtained at the Canada–France–Hawaii Telescope (CFHT), operated by the National Research Council of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawai‘i, including data acquired through the Gemini Remote Access to CFHT’s ESPaDOnS Spectrograph. Additional observations were obtained with the Liverpool Telescope, operated on La Palma by Liverpool John Moores University at the Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, with financial support from the UK Science and Technology Facilities Council, as well as from the Las Cumbres Observatory global telescope network, the 0.91-m telescope of the Osservatorio Astrofisico di Catania, the Crimean Astrophysical Observatory, and the Skinakas Observatory, supported by the European Social Fund (ESF) and National Resources under the Operational Programme Education and Lifelong Learning. CCP was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2019R1A6A1A10073437) This work made use of the LSDpy code (https://github.com/folsomcp/LSDpy) and the ZDIpy code (https://github.com/folsomcp/ZDIpy).

References

Appendix A: Observations

A.1. LCOGT photometry

Photometric observations were obtained at the Las Cumbres Observatory Global Network (LCOGT, Brown et al. 2013) with the 1m Sinistro telescopes from October 16 to November 30, 2019 (LCO 2019B-019, PI L. Rebull; CLN2019B-004, PI A. Bayo). We used SDSS/PanSTARRS g′, r′, and i′ filters with exposure times of 400s, 120s, and 80s, respectively. A total of 184 images centered on JH 223 were obtained over 47 days, yielding 46, 59, and 60 photometric measurements in the g′, r′, and i′ filters, respectively, after discarding some low-quality images. We used the photometric catalogs provided by LCOGT’s BANZAI reduction pipeline to retrieve the instrumental magnitudes of the target as well as those of nearby comparison and control stars. We chose JH 221, a nearby object of comparable magnitude as JH 223, as a comparison star. We checked JH 221 against XEST 07-OM-002 (2MASS J04404215+2552403), another nearby star of similar brightness. We found the two stars to be nonvariable over the observed time period and deduced photometric root-mean-squared (rms) uncertainties of 0.009, 0.014, and 0.010 mag in the g′, r′, and i′ filters, respectively. The g′, r′, and i′ differential light curves obtained for JH 223 using JH 221 as a comparison star were subsequently flux calibrated to within an rms accuracy of 0.02mag using 33 contemporaneous Cousins V, R_c, and I_c measurements of the same field obtained at Crimean Astronomical Observatory (CrAO).

A.2. OACT photometry

Photometric data were taken at the M. G. Fracastoro station (Serra La Nave, Mt. Etna, 1750m a.s.l.) of the Osservatorio Astrofisico di Catania (OACT, Italy) from October 21, 2019, to February 25, 2020, with the 0.91 m telescope. The CCD camera adopts a Kodak KAF1001E² 1k×1k chip that, with a focal reducer, covers a field of view (FOV) of about 11.5×11.5 arcminutes. We observed JH 223 with a set of broad band Bessel filters (V, R, I) as well as two narrow band filters, Hα₉ and Hα₁₈, centered at 656.8 nm and at 676.4 nm, whose full widths at half-maximum are 9 and 18 nm, respectively. We collected a total of 43, 70, and 73 useful images in V, R, and I bands, respectively, with exposure times of 30–40 s, 50–80 s, and 90–180 s, in the respective bands. Regarding the narrow band photometry, we collected 39 images in the Hα₉ filter, with exposure times of 180–300 s, and 78 images in the Hα₁₈ filter, with exposure times of 90–150 s. The latter were typically acquired just before and after the Hα₉ frames to take sky transparency variations into account.

The broad band V, R, and I photometry was calibrated using the following procedure. The stars in the standard area GD 71 (Landolt 2009) were observed in the nights with the best photometric conditions to calculate the zero points and transformation coefficients to the Johnson-Cousins system. These transformations were applied to the two non-variable stars in the FOV of JH 223, namely XEST 07-OM-002 and XEST07-OM-001 (2MASS J04403830+2554274), which were then used as comparisons to derive the magnitudes of JH 223 in each image³.

For details about the OACT data reduction and the narrow band Hα photometry, the reader is referred to Frasca et al. (2018). We remark that the CI_Hα = Hα₁₈ − Hα₉ is basically a measure of the intensity of the Hα emission in units of the continuum that can be converted into Hα equivalent width (EW_Hα), using the calibrations proposed by Frasca et al. (2018). Thanks to the contemporaneous GRACES and Keck high-resolution spectra, we have improved the aforementioned CI_Hα − EW_Hα calibrations by means of synthetic color indices measured on the latter ones, integrating them in the passbands of the OACT Hα filters.

A.3. CrAO photometry

Short (60, 15, and 5s in V, R_j, and I_j, respectively) and long (180s, 180s, and 120s in V, R_j, and I_j, respectively) exposure images were obtained from September 2 to November 19, 2019, at the Crimean Astronomical Observatory through Johnson’s V, R_j and I_j filters, using the AZT-11 1.25m telescope equipped with the CCD camera ProLine PL23042. CCD images were bias-subtracted and flatfield corrected following a standard procedure. On short exposures, we used XEST 07-OM-002 as a comparison star, for which we derived V = 16.83, (V_j − R_j) = 2.15, (V_j − I_j) = 2.79. The resulting photometric errors on long exposures typically amount to 0.09, 0.02, and 0.02mag in the V, R_j, and I_j filters, respectively. We finally converted the magnitude and colors of JH 223 obtained in Johnson’s system to Cousins’s, using the following relationships: (V_c − R_c) = (V_j − R_j)×0.71652 − 0.03206; (V_c − I_c) = (V_j − I_j)×0.7822 + 0.0065 (Landolt 1983). We collected a total of 50 useful V, R_j, and I_j images between September 02 and November 19, 2019.

A.4. Liverpool photometry

SDSS u′, g′, r′, and i′ photometry was performed with the 2m Liverpool Telescope at the Roque de los Muchachos Observatory (La Palma, Spain) using the IO:O camera equipped with the E2V detector. A total of 48 images were obtained from November 6 to 30, 2019, centered on JH 223 with a 10 arcmin square FOV, and with exposure times of 580s, 90s, 20s, and 10s in the u′, g′, r′, and i′ filters, respectively. The images were reduced by the IO:O pipeline, which includes bias subtraction, trimming of the overscan regions, and flat fielding. PSF photometry was performed with IRAF/DAOPHOT on the reduced images for JH 223, as well as for the two comparison stars XEST 07-OM-001 and XEST 07-OM-002. We estimated the photometric errors from the differential light curve computed between the 2 comparison stars, which yielded an rms error of 0.010mag in g′, r′, and i′ filters. The objects are too faint in the u band to yield reliable photometry, and measurements in this filter were thus discarded. We computed the differential light curve between JH 223 and XEST 07-OM-001, and calibrated it using contemporaneous LCOGT photometry.

A.5. K2 photometry

JH 223 (EPIC 248006676) was observed by the K2 satellite during Campaign 13, which took place over 80 days from March 8 to May 23, 2017, namely, about 2.5 years before our campaign. The observations were performed in a broadband filter (420–900nm) with measurements taken at a cadence of 30 minutes. The K2 light curve was reduced by A. M. Cody using a 3-pixel circular aperture with background subtraction and subsequent correction of pointing jitter, as described in Cody & Hillenbrand (2018) and Cody et al. (2022)⁴.

A.6. TESS photometry

The TESS satellite observed JH 223 during two epochs, from September 16 to November 05 in 2021 (sectors 43 and 44) and from September 20 to November 11 in 2023 (sectors 70 and 71), a few years after our campaign. The observations were obtained using the red-optical bandpass (600–1000nm) with a cadence of 2 minutes. In this work, we selected the reduction made by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016).

A.7. Skinakas observatory photopolarimetry

Observations were performed over two nights, on October 14 and 30, 2019, with the RoboPol photopolarimeter (Ramaprakash et al. 2019) mounted on the 1.3m telescope at Skinakas Observatory, Crete, gathering two measurements in the Cousins-R band and an additional one in the SDSS g filter. The target was placed in the center of the field of view, where a focal-plane mask produces shadows to reduce the background light. Polarized and unpolarized standard stars were observed on the nights of the target observations and were used in the calibration. The data were reduced using the RoboPol pipeline (King et al. 2014), which performs differential aperture photometry on the set of four target images produced by the instrument to measure the relative Stokes parameters Q and U. The results are listed in Table A.1.

A.8. Gemini/GRACES optical spectroscopy

The optical spectroscopic observations were obtained with the Gemini Remote Access to CFHT ESPaDOnS Spectrograph (GRACES), which consists of the ESPaDOnS high-resolution spectrograph (Donati 2003) at the Canada-France-Hawaii Telescope (CFHT) connected through two 270m optical fibers to the 8.1m Gemini North telescope (Chene et al. 2014). The observations were obtained over 4 nights between November 7, 2019, and January 20, 2020, with a resolution R ∼ 40000 in the wavelength range of 370nm to 950nm (proposal ID: GN-2019B-Q-101; PI: S. Alencar). The data reduction was done with the OPERA pipeline (Martioli et al. 2012), which includes bias subtraction, flat-field correction, and wavelength calibration. All spectra were continuum normalized.

A.9. Keck I optical spectroscopy

We also obtained optical spectroscopic observations with the Keck I Telescope using the High-Resolution Echelle Spectrometer (HIRES, Vogt et al. 1994), with a spectral resolution of R ∼ 48000, covering the wavelength range from 470nm to 930nm. A total of four spectra were acquired between November 29 and December 1, 2019. However, only three were included in the analysis, as one exhibited a low S/N. The data reduction was done using the MAKEE⁵ reduction pipeline written by T. Barlow.

A.10. CFHT/SPIRou near-infrared spectropolarimetry

JH 223 was observed for eight non-consecutive nights in November 2019 using the near-IR spectropolarimeter and velocimeter SPIRou at the CFHT (proposal ID: 19BF02; PI.: J. Bouvier). SPIRou is a stabilized high-resolution near-infrared (nIR) spectropolarimeter that covers wavelengths between 950nm and 2350nm with a resolution R ∼ 70000 (Donati et al. 2020b). SPIRou observations consist of sequences of four subexposures taken at distinct orientations of the quarter-wave Fresnel rhombs, where each subexposure (integrated over 903s for JH 223) includes two spectra with orthogonal polarization states (Donati et al. 1997; Bagnulo et al. 2009).

The data reduction was performed using the APERO pipeline software (Cook et al. 2022). In addition to the standard extraction, flat-field correction, and wavelength calibration, APERO performs the polarimetric data-reduction, combines subexposures to extract normalized unpolarized (Stokes I), circularly polarized (Stokes V), and null (Stokes N) spectra that are telluric-corrected following the method described in Artigau et al. (2022). The SPIRou data are publicly available from the Canadian Astronomy Data Center⁶.

Appendix B: Circumstellar lines from GRACES, Keck, and SPIRou observations

Figure B.1 shows the Hα, Hβ, and He I 587.6 nm emission lines from GRACES and Keck observations before (black) and after (red) the removal of the photospheric contribution. For the near-infrared circumstellar emission lines of He I 1083nm and Paβ, the lines profiles obtained from SPIRou observations are shown in Figure B.2.

Fig. B.1. Circumstellar emission lines of Hα (top row), Hβ (middle row), and He I 587.6nm (bottom row) from the GRACES (Cols. 3, 4, 6, and 7) and Keck (Cols. 1, 2, and 5) observations. The blue curves show the spectrum of TWA 7, used as a photospheric template. The black and red curves correspond to the spectrum of JH 223 before and after removal of the photospheric contribution, respectively. The photosphere-subtracted spectra are vertically shifted for clarity.

Fig. B.2. Similar to Fig. B.1, but for the He I 1083nm (top) and Paβ (bottom) emission lines from the SPIRou observations.

Table B.1 summarizes the equivalent widths, line fluxes, and line luminosities of the Hα and Hβ lines measured from the GRACES and Keck spectra. As described in Sect. 3.2.4, these quantities were used to derive the accretion luminosity and mass accretion rate of JH 223, which are also reported in Table B.1. In Fig.B.3, we illustrate the rotational modulation found in the accretion luminosity derived from the Hα line.

Fig. B.3. Accretion luminosity of Hα as a function of rotational phase. The yellow circles represent the accretion luminosities obtained from the equivalent widths derived from the narrowband color-index filters CIHα, while the blue circles correspond to the values derived from the GRACES and Keck spectra.

Appendix C: Details about the ZDI reconstruction

C.1. The χ² aiming criteria

We computed a series of ZDI reconstructions aiming at different values of ${\chi }_r^2$. We looked at the variation of the image entropy as a function of ${\chi }_r^2$ and applied the Kneedle algorithm (Satopaa et al. 2011), implemented in the Python package Kneed⁷, as an approach to search for ${\chi }_{r,sf}^2$. The Kneed algorithm identifies inflection points in curves characterized by negative or positive concavity in discrete datasets, using the mathematical definition of curvature for continuous functions. According to Satopaa et al. (2011), this method has shown greater effectiveness compared to other numerical approaches proposed to detect knees in discrete data, such as Angle-based and Menger Curvature (Basseville & Nikiforov 1993; Bollinger 2002; Zhao et al. 2008). Figure C.1 shows an example of the Kneedle algorithm applied to the reconstructed brightness images of JH 223. Through the Kneed code, we obtained the inflection point ${\chi }_{r,sf}^2=0.53$. We applied the correction to the error bars of all LSD profiles, which allowed us to target a unit ${\chi }_r^2$ in the following image reconstructions.

Fig. C.1. Entropy versus ${\chi }_r^2$ for several models of the reconstructed brightness images of JH 223. Each point corresponds to a converged ZDI solution. The vertical line indicates the inflection point of ${\chi }_{r,sf}^2=0.53$ that was determined using the Kneedle algorithm.

C.2. Best entropy model

Figure C.2 shows the maximum-entropy solution obtained by applying ZDI technique to the time-series of Stokes I and V LSD profiles computed from the SPIRou dataset.

Fig. C.2. ZDI fit of the Stokes I (left) and V (right) LSD profiles of JH 223. Observed (black) and maximum-entropy synthetic (red) profiles are shown. Dashed blue lines indicate the continua. Profiles are vertically shifted for clarity and labeled by rotational phase.

All Tables

All Figures

	Fig. 1. JH 223 K2 light curve shown as a function of Julian date (top) and rotational phase (bottom). The rotational phase was computed using a period of 3.31days, with the reference time JD0 = 2 457 817.63 chosen such that phase 0.5 corresponds to the photometric minimum. The color coding indicates to the Julian date of each observation.
In the text

	Fig. 2. JH 223 V-band light curve versus HJD (top) and rotational phase (bottom). The phased light curve uses a period of 3.31days and JD0 = 2 458 770.2388 (phase 0.5 at minimum). Arrows mark the epochs of spectroscopic and spectropolarimetric observations. CrAO and OACT data show larger scatter due to higher uncertainties.
In the text

	Fig. 3. CMDs of (V − R)c vs. V (green) and (V − I)c vs. V (red), with (V − R)c shifted by +1.6mag for clarity. LCOGT (circles) and LT (squares) data are shown. Solid and dashed lines indicate least-squares fits and ISM-extinction slopes, respectively.
In the text

	Fig. 4. Same as Fig. 1, but for the 2021 TESS light curve from sectors 43 and 44. The rotational phase was computed with the 3.31-day period and a reference time of JD0 = 2 459 471.7.
In the text

	Fig. 5. Same as Fig. 1, but for the 2023 TESS light curves from sectors 70 (a) and 71 (b). The same period and reference time as in Fig. 4 were used.
In the text

	Fig. 6. Example of a spectral window used for the stellar parameter fitting. The corresponding synthetic spectrum (red line) is overplotted on the observed GRACES spectrum (blue line).
In the text

	Fig. 7. Hα (top), Hβ (middle), and He I 587.6 nm (bottom) emission lines used to derive mass accretion rates. Labels indicate the rotational phase. GRACES and Keck data are shown as solid and dashed lines.
In the text

	Fig. 8. Residual line profiles of HeI 1083nm. The colors represent different rotational phases.
In the text

Fig. 9. JH 223 radial velocities (top) and residuals (bottom). Radial velocities were derived from the first-order moment of the Stokes I LSD profiles of the SPIRou observations (black dots) and from spectral fitting of GRACES (red dot) and Keck spectra (green dots). The blue curve represents the sinusoidal fit with the rotational period fixed at 3.31days. Residuals (squares), computed as the difference between the RV values and the sinusoidal fit, have a standard deviation of 0.27km/s.

In the text

	Fig. 10. Same as Fig. 9, but for the longitudinal magnetic field. Residuals have a standard deviation of 33 G.
In the text

	Fig. 11. Logarithmic brightness map of the surface of JH 223 in November, 2019. The star is shown in a flattened polar view down to a latitude of −30°, with the equator indicated by a solid line and latitudes of 60° and 30° by dashed lines. Outer radial ticks mark the phases of spectropolarimetric observations. Cool spots are shown in brown shades, while bright plages are shown in blue shades.
In the text

	Fig. 12. ZDI maps of the radial (top), meridional (middle), and azimuthal (bottom) components of the large-scale magnetic field at the surface of JH 223. Similar to Fig. 11, the star is represented in a flattened polar projection. Magnetic fluxes, indicated by the color bar, are expressed in gauss.
In the text

	Fig. 13. Asymmetry (M) and quasiperiodicity (Q) parameters derived from the K2 and TESS light curves of JH 223. The gray points correspond to the Q and M parameters of others CTTSs analyzed by Cody et al. (2022) in the Taurus star-forming region.
In the text

Fig. 14. HRDs showing the large-scale magnetic properties of JH 223 and other classical and weak-line T Tauri stars. CTTSs are labeled by name. Symbol size scales with the average magnetic field strength, color with the poloidal magnetic energy fraction, and shape with the fraction of poloidal energy in axisymmetric modes. Evolutionary tracks (black lines) and isochrones (green dotted lines) from Baraffe et al. (2015) are shown, together with the fully convective boundary and the 40% radiative-core mass limit (dashed blue lines). Magnetic field data for the T Tauri stars are taken from Donati et al. (2014, 2015), Yu et al. (2017), Hill et al. (2017), Nicholson et al. (2018), Hill et al. (2019), Yu et al. (2019), Nicholson et al. (2021), Xiang et al. (2023), Donati et al. (2007, 2008, 2010a,b, 2011a,b, 2012, 2013, 2019, 2020a), Zaire et al. (2024), Donati et al. (2026).

In the text

Fig. B.1. Circumstellar emission lines of Hα (top row), Hβ (middle row), and He I 587.6nm (bottom row) from the GRACES (Cols. 3, 4, 6, and 7) and Keck (Cols. 1, 2, and 5) observations. The blue curves show the spectrum of TWA 7, used as a photospheric template. The black and red curves correspond to the spectrum of JH 223 before and after removal of the photospheric contribution, respectively. The photosphere-subtracted spectra are vertically shifted for clarity.

In the text

	Fig. B.2. Similar to Fig. B.1, but for the He I 1083nm (top) and Paβ (bottom) emission lines from the SPIRou observations.
In the text

	Fig. B.3. Accretion luminosity of Hα as a function of rotational phase. The yellow circles represent the accretion luminosities obtained from the equivalent widths derived from the narrowband color-index filters CIHα, while the blue circles correspond to the values derived from the GRACES and Keck spectra.
In the text

	Fig. C.1. Entropy versus ${\chi }_r^2$ for several models of the reconstructed brightness images of JH 223. Each point corresponds to a converged ZDI solution. The vertical line indicates the inflection point of ${\chi }_{r,sf}^2=0.53$ that was determined using the Kneedle algorithm.
In the text

	Fig. C.2. ZDI fit of the Stokes I (left) and V (right) LSD profiles of JH 223. Observed (black) and maximum-entropy synthetic (red) profiles are shown. Dashed blue lines indicate the continua. Profiles are vertically shifted for clarity and labeled by rotational phase.
In the text