© ESO 2020

1 Introduction

Dippers, also called AA Tau-type stars, are young stars with lightcurves characterized by aperiodic dimmings, or eclipsing events consistent with variable extinction by circumstellar material (Bouvier et al. 1999). Dippers are particularly interesting sources, because the occultations of the star by the disk material offer an opportunity to explore the structure and composition of the innermost disk (Schneider et al. 2018) that are inaccessible for most other sources.

The solar-mass star RX J1604.3-2130A¹ is a member of the 5–11 Myr old Upper Scorpius Association (Preibisch & Zinnecker 1999; Pecaut et al. 2012; Pecaut & Mamajek 2016) that has been identified as a dipper (Ansdell et al. 2016). It possesses one of the brightest disks detected in the region, which also has a large inner cavity (Carpenter et al. 2006, 2009, 2014; Dahm & Carpenter 2009; Mathews et al. 2012; Zhang et al. 2014; Dong et al. 2017), and is therefore considered a transition disk. The resolved outer disk is nearly face-on, with an estimated inclination of 6° (Zhang et al. 2014; Dong et al. 2017). The disk contains a substantial amount of gas (Mathews et al. 2013), although the cavity shows significant CO depletion (Dong et al. 2017; Mayama et al. 2018).

The disk gap could be the result of planetary formation or the presence of a yet-undetected stellar companion. A low-mass (M2) companion, which is itself a binary with a stellar companion at 0.082 arcsec (Köhler et al. 2000), is found at 16″ (RX J1604.3-2130B). Considering the Gaia DR2 parallax (6.662 ± 0.057 arcsec; Gaia Collaboration 2016, 2018), available through VizieR (Gaia Collaboration 2018), its distance is 150 ± 1 pc, which implies that the potential companions are located at about 2400 au and 12 au (from B), respectively. Nearby companions of RX J1604.3-2130A at >22 au have been excluded down to 2–3 M_J (Kraus et al. 2008; Canovas et al. 2017), but there is no record of objects in the innermost region of the gap.

The disk was observed with VLT/SPHERE, revealing a close-to-face-on ring in scattered light of about 65 au in radius (Pinilla et al. 2015, 2018). Dark dents on the scattered light image of the disk rim suggested shadows cast by a highly misaligned inner disk (as it has been observed in other systems; Marino et al. 2015; Benisty et al. 2017, 2018). The origin of the shadows in an inclined inner disk was recently confirmed by multi-epoch scattered light observationsbetween 2016 and 2017, which show that the position angle of the dips varies only slightly (with PA ~ 83.7 ± 13.7° and PA ~ 265.9 ± 13.0°, respectively Pinilla et al. 2018), although the morphology of the shadows is quite variable on timescales of days. The presence of a misaligneddisk has also been suggested from ALMA gas observations (Mayama et al. 2018) and is also in good agreement with the eclipsing activity observed, but a classical, smooth inner disk cannot explain the optical and scattered light variability observed.

In this paper, we use ground-based optical and near-IR photometry from REM/La Silla, together with K2 data and archival photometry and spectroscopy of RX J1604.3-2130A, to explore the causes behind the observed eclipses, and to constrain the structure of the innermost disk and its variability timescales. Observations are presented in Sect. 2. The periodicity of the lightcurves and variability causes in terms of eclipses, accretion, rotation, and variations in the inner disk structure are analyzed in Sect. 3. The discussion and conclusions are presented in Sects. 4 and 5.

2 Observations and data reduction

2.1 REM g′ r′i′z′ JHK observations

RX J1604.3-2130A was observed with the Rapid Eye Mount (REM) 60 cm telescope² in La Silla, as part of a DDT proposal. With its two instruments, ROS2 and REMIR, REM can obtain nearly-simultaneous images in the Sloan g′, r′, i′, z′, and the IR JHK filters over a 10 × 10 arcmin² field. The observations took place at irregular intervals during approximately five months, from 2018 May 09 to 2018 October 01 (MJD 58 247.21–58 392.08). The observations were repeated at intervals of between half an hour and few days, with a period of high-cadence on 58 255.09 (2018 May 16–17), during which the images were taken about every two minutes. For the optical bands, a total of 338 images were obtained, among which 168 of them belong to the high-cadence dataset. The total exposure time for the optical images was 5s. Data reduction was performed with the automated REM pipeline, and the images were aligned using Astrometry.net. The IR observations were obtained during the same epochs, but due to technical issues, only 42 datasets were completed for J, 43 for K, and 142 for H. Typically, a total of five dithered 3s exposures were taken in each case, although some of them have less dithers if the imaging failed (especially on bad nights). Sky images were also acquired, although some of the final images have less dithers due to the lack of quality of some of them. The images were reduced, combined, sky-subtracted, aligned using the automated REM pipeline, and aligned using Astrometry.net.

Aperture photometry was performed using iraf task noao.digiphot.apphot. For the optical data, a relative calibration was performed for each filter by comparing all observations to the data taken on MJD 58 346.021 (one of the best nights, based on seeing/FWHM). An iterative process based on the median and standard deviation of the magnitude difference (Sicilia-Aguilar et al. 2008, 2017) was used on all stars in the field to identify the non-variable comparison stars. The main issue is that there are few comparison stars in the field, and all of them are fainter than the target star. This is due to RX J1604.3-2130A brightness (Gaia G 11.87 mag) requiring low exposure times to avoid saturation. We thus imposed quality limits on the calibration, rejecting all those for which four or fewer comparison stars were identified, and those where the calibration had magnitude-dependent offsets or very large errors. Typically, we found between five and 15 comparison stars for g′ and r′, and between 10 and 27 comparison stars for i′ and z′, spanning a range of two to four magnitudes around (but mostly fainter than) the object. The nights for which the calibration fails are typicallythose with poor seeing and poor weather conditions that make it hard to detect enough field stars. This results in 205 dates for which all four optical filters are complete. The final data (relative magnitudes) are listed in Table 1. The uncertainties provided in the tables and figures include the photometric uncertainty and the uncertainty in the relative calibration. In general, the uncertainty in the relative calibration dominates the value in cases with few comparison stars.

The absolute calibration of the g′ r′ i′ z′ data was done using griz data from PAN-STARRS³ (Chambers et al. 2016; Flewelling et al. 2016; Magnier et al. 2016). The procedure for the absolute calibration was similar to the relative calibration, comparing the data from MJD 58 346.021 with those of PAN-STARRS. The g′, r′ calibration is quite robust (2 and 4% uncertainties, respectively), but for i′, there is a large amount of scatter, and for z′, only four reference stars could be identified. In addition, we find that there are color terms for all filters except g′. The color terms are particularly large for i′ and z′. Since there are no stars as bright as RX J1604.3-2130A in the field, and very few non-variable stars in any magnitude range, the absolute calibration is very uncertain, and the errors could be >50% in i′ and z′. The magnitudes in the reference night MJD 58 346.021 are thus g′ = 12.47 ± 0.02 mag, r′ = 11.01 ± 0.04: mag, i′ = 11.7: mag, and z′ = 11.0: mag. Because the magnitudes in r′,i′, and z′ are estimated from much fainter stars, we label them as uncertain (:) and treat them as merely indicative, and we focus the discussion on relative REM magnitudes and relative color variations. In total, there are 261 datapoints in g′, r′ and i′, and 270 in r′ observations.

For the JHK data, we followed a similar procedure. The data were calibrated against those of MJD 58 363.04 (one of the best nights for which all three IR filters were obtained), which was also calibrated against 2MASS data (Skrutskie et al. 2006). The REM and 2MASS JH filters are essentially identical, while the relation between the standard Johnson K filter and the 2MASS K_s is K = K_s + 0.044 mag (Bessell 2005). The systematic errors in the absolute calibration are 4% for J and H and 2% for K. There are no evident color terms between the 2MASS and REM filters, but because the data quality of the JHK images is generally worse than in the optical, there are fewer stars for comparison. Following the same quality criteria, we only have complete JHK data on 10 epochs, although the H lightcurve is much more complete (73 points), while J and K only have eight. The final results are listed in Table 2.

The final REM lightcurve is displayed in Fig. 1. The optical data show the typical dimming events described in Ansdell et al. (2016). The JHK data follows the behavior observed in the optical, although they are more scarce. The high-cadence data reveal the star emerging from one of the eclipses. We find that the eclipses observed by REM are deeper than previously reported (0.57 mag based on K2 data; Ansdell et al. 2016). Even though our filters are significantly narrower than the K2 filter, the maximum depth varies between 0.4 and 1.8 mag in g′, 0.3 and 1.5 mag in r′, 0.3 and 1.5 mag in i′, and 0.2 and 1.2 mag in z′. The variations in H (the only IR filter for which we have enough eclipse data) are up to 0.2–0.7 mag. There are also smaller variability events, but we detect at least 10 deep eclipses during our observations (all of which are recovered inmultiple bands), in addition to other shallower ones similar to those reported by Ansdell et al. (2016). The decrease in depth with increasing wavelength suggests extinction events, which we explore in Sect. 3.

2.2 Other optical lightcurves: K2 and CSS

RX J1604.3-2130A was observed by K2 as part of the Ecliptic Plane Input Catalog (EPIC; Huber et al. 2016) as source EPIC 204 638 512. The K2 data was obtained from the Mikulski Archive for Space Telescopes (MAST⁴) at the Space Telescope Science Institute⁵. Although some K2 data are affected by interlopers, for RX J1604.3-2130A the MAST K2SFF public lightcurves have been validated so that there is no need of further corrections (Ansdell et al. 2016). The data were acquired between 2014 August 23 and 2014 November 10 (MJD 56 892.78–56 971.55), thus about four years before the REM data. The K2 data are essentially uniformly sampled, with a sampling period of 30 min, and are presented as relative fluxes, thus having values around one out of eclipse. Figure 2 displays the complete K2 lightcurve, which shows the irregularly-shaped dimming events up to 0.57 mag reported by Ansdell et al. (2016). We also note that, besides the dips, there are several cases where a sudden brightness increase is observed. These may be stellar flares and are discussed in Appendix A.

RX J1604.3-2130A has been observed by the Catalina Sky Survey (CSS; Drake et al. 2009), Data Release 2⁶. The CSS archivecontains 293 photometry points distributed over nearly eight years, from 2005-08-01 to 2013-07-22 (MJD 53 583.454–56 495.535: see Fig. 2). The object ID in the survey is SSS_J160421.7–213 028. The sampling is very sparse compared to the relevant timescales, but it shows the same behavior detected in the REM and K2 data, with sudden dimmings that, in some cases, go down by nearly 1.4 magnitudes in V and some periods of relative stability. Some of the data are very close to the saturation limit (11 mag), to the extent that the highest magnitudes may be uncertain, but the eclipse data are well below saturation. Although no further information can be obtained from these data regarding periodicity, they essentially confirm the behavior observed and the fact that the eclipse depths are highly variable and persistent. Note that the Catalina V filter has a non-negligible color term⁷. The color term is stronger for very red objects, so it likely affects the eclipse depth. Considering the color variations observed with REM and the typical colors for a K3-type star (Kenyon & Hartmann 1995), the maximum eclipse depth in V (Cousins system) may be shallower by up to 0.3–0.4 mag with respect to the value observed in the CSS lightcurve.

Fig. 1 REM lightcurve for RX J1604.3-2130A. Upper panel: full lightcurve observed, while lower panel: zoom in the high-cadence data. All magnitudes shown are relative magnitudes evaluated with respect to the σ-clipped average magnitude of all datapoints in each band, marked as zero point with a dotted line.

Fig. 2 Top: K2 lightcurve for RX J1604.3-2130. Typical uncertainties are smaller than dots. Bottom: Catalina Survey DR2 lightcurve. Since Catalina data cover the epochs of the Spitzer, WISE, and HIRES spectroscopy observations, we have marked them in the figure as vertical lines (see text). In both cases, the average flux (magnitude) is shown as a dotted line.

2.3 Archival optical spectroscopy

With the aim of understanding the causes of variability, we also need to constrain rotation and accretion, which are two of the major causes leading to magnitude fluctuations observed in young stars. We thus study archival high-resolution spectroscopy in the analysis. RX J1604.3-2130A was observed six times with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) between June 2006 and April 2010, available through the Keck Observatory Archive (KOA⁸). Exposure times, coverage, and resolution varied and are listed in Table 3. The data were reduced using the automated MAKEE pipeline⁹. The automated reduction includes bias and flat field correction and calibration using a ThAr lamp. The long-slit data were used to extract and subtract the sky spectrum. The spectra were extracted in vacuum wavelength and subsequently transformed using PyAstronomy¹⁰ routine vactoair2. No flux calibration was performed.

Besides looking for variability signatures, the spectrum with the best S/N (MJD 55 287.614) was also used to confirm the projected rotational velocity (v sin i) of the object. We selected the region between 5500 and 5800 Å, which is relatively devoid of both accretion- and activity-related emission lines and telluric lines (Curcio et al. 1964), and we measured the rotational and radial velocity by cross-correlating the object spectrum with three different rotational standards with similar spectral types¹¹ that had been observed with Keck under similar conditions. These included HD 114386 (K3, also used by Dahm et al. 2012), HD 10780 (K0), and HD 151541 (K1).

PyAstronomy task rotBroad was used to create artificially broadened templates, and crosscorRV was used to obtain the cross-correlation. The location of the cross-correlation peak was used to determine the radial velocity (v_rad), and the width of the cross-correlation function was compared with that of the broadened templates to obtain the rotational velocity v sin i. We obtained v_rad = −6.8 ± 0.1 km s⁻¹ and v sin i = 16.2 ± 0.6 km s⁻¹. Both are in good agreement with Dahm et al. (2012), and confirm that the star is a relatively fast rotator compared to young stars with similar spectral types (e.g. Sicilia-Aguilar et al. 2005; Weise et al. 2010), especially if we take into account that the system is accreting (Dahm et al. 2012). Although the rest of datasets are significantly worse in quality, the results derived are consistent (see a summary in Appendix B).

2.4 Further archival data and stellar parameters

RX J1604.3-2130A has been repeatedly observed in the mid-IR by WISE and Spitzer (Carpenter et al. 2006; Luhman & Mamajek 2012; Esplin et al. 2018). The star has some signs of intriguing mid-IR variability, changing from photospheric colors as observed with Spitzer on MJD 53 820, to clear mid-IR excess as observed with WISE after MJD 55 249 (Luhman & Mamajek 2012). We thus include the Spitzer and WISE data in the discussion, taking the IRAC photometry values reported by Carpenter et al. (2006) and the available AllWISE Multiepoch Photometry¹² (Wright et al. 2010) lightcurve data provided by IRSA¹³.

To examine the disk structure and behavior, we use the multiwavelength data available within VizieR to construct the spectral energy distribution (SED¹⁴; see Table C.1 for wavelengths, fluxes and references). The available data include SDSS, Gaia, POSS-II, hipparcos, and SkyMapper optical data; 2MASS, POSS-II, UKIDSS, and Vista near-IR data; and Spitzer, Akari, IRAS and IRAS mid- and far-IR data, together with the 880 μm datapoint from ALMA. The magnitude data are converted to fluxes using the same calibrations as VizieR (Bessell & Brett 1988; Fukugita et al. 1996; Cohen et al. 2003).

Finally, in the whole discussion, we use the most recent estimates of the stellar parameters derived from X-shooter spectroscopy (see Manara et al., in prep.), which give a K3 spectral type, T_eff = 4730 K, L_* = 0.90 L_⊙ ¹⁵ (so the stellar radius is R_* = 1.4 R_⊙), and M_* = 1.24 M_⊙ (using the spectral type calibration and tracks from Luhman et al. 2003; Baraffe et al. 2015, respectively) and estimates an accretion rate of 3e-11 M_⊙ yr⁻¹ at a time when the star was in a relatively bright state. Note that the accretion rate is on the limit of what can be detected with X-shooter (the object is classified as a potential accretor by Dahm et al. 2012, based on its weak accretion features), which adds uncertainty to the measured value, although the line profiles seen with HIRES show clear accretion. Using these stellar parameters, we can also estimate that the dust sublimation radius (for T = 1500–1000 K) is located at about 0.06–0.15 au (~10–22 R_*).

3 Analysis

3.1 Periodicity analysis

A period of about 5 d (albeit very uncertain) has been suggested as the rotational period of RX J1604.3-2130A from K2 data (Ansdell et al. 2016; Rebull et al. 2018). Here we revisit the periodicity in the K2 lightcurve to examine whether it is most likely due to rotation, or related to the obscuration events. The lightcurve is extremely irregular, suggesting variations in both the phase, the period, and the amplitude of the modulations and the presence of correlated, non-Gaussian noise. Therefore, we take two approaches to search for periodical signatures: generalized Lomb-Scargle periodograms (GLSP; Scargle 1982; Horne & Baliunas 1986; Zechmeister & Kürster 2009) and wavelet analysis (Torrence & Compo 1998; Liu et al. 2007).

Simple GLSP fail when applied to quasi-periods, so, for the first approach, we use stacked GLSP (SGLSP; Mortier & Collier Cameron 2017), where the data are filtered around each single date to study periodicity only over a limited number of days. Repeating the exercise over time, changes in the period and phase can be tracked. Since the data distribution is highly non-Gaussian and strongly correlated, a red noise model is needed to assess the significance of the signatures. The red noise model was derived from the correlation between consecutive datapoints, parameterized by α, the slope of the correlation between one datapoint and the next. For K2 data, we find α = 0.98. We then simulated data following a similar distribution where each value depends on the previous one (as given by α) plus a stochastic component with values drawn from random numbers following the standard deviation of the observed data. For each K2 (S)GLSP, we computed 1000 red noise simulations with the same number of datapoints, the same sampling rate, and the same red noise model but without any periodic signature and used their periodograms to derive the confidence intervals.

A GLSP including all the available data suggests some signals with periods 2.5 d (which could correspond to half the 5 d period), 5 d, and 9 d, but none of them is highly significant (see Fig. 3). The 5 d periodic signature is reported in the literature as a rotational period, so we first examined the data belonging to the out-of-eclipse phase, and the data belonging to the eclipses,separately. Defining the “out-of-eclipse phase” is not easy in a lightcurve that does not show a clear difference between the “in-eclipse” and “off-eclipse” parts (see Fig. 4). We thus cleaned the data using a σ-clipping algorithm to calculate the mean and standard deviation, and then removed all points that are beyond 3σ from this value. The analysis of the off-eclipse data revealed that the periodic signatures are stronger when the full dataset is considered, and thus the (quasi-)periodicity is strongly linked to the eclipses, and not only to rotation.

We then stacked the data for different numbers of days around each point and calculated SGLSP in intervals of ± 10 d, ± 20 d, ± 30 d, and ± 40 d. The number of days was selected to be large enough to detect the potential periods observed in the full collection of data, up to the limit of ± 40 d that includes essentially the whole dataset, and thus converges to the full-data GLSP. Figure 5 displays the results. The data reveals a 5 d period in the stacked ±10 d and ±20 d diagrams, which progressively dilutes when more data are added. The period is not always present, and the peak changes between 4.8 and 5.5 d, which makes it more plausibly a quasi-period related to a rapidly variable phenomenon than a typical rotational period, evenif both may be connected. The 2.5 d period is also present, although it has a lower significance, except during the time of increased eclipse activity (approximately, from mission day 2085 to 2105), when eclipses occur at a higher rate. The 9 d period is very broad and not well-defined.

For the second periodicity estimate, we use wavelet analysis based on a Morlet function with ω₀ = 6, which typically offers the best results for complex datasets (Torrence & Compo 1998). The Morlet wavelet Ψ(η) (Grinsted et al. 2004) is very similar to a sinusoidal function tappered by a Gaussian, written as $$\Psi (\eta )=\frac{e^{i{\omega }_0\eta }}{{\pi }^{1/4}}{e}^{-{\eta }^2/2}.$$(1)

Here, the time dependency is wrapped in the dimensionless parameter η, which takes values that are multiple of powers of two assigned through the dimensionless time series¹⁶. In essence, the SGLSP and wavelet analysis are quite similar, with the main difference being that SGLSP uses a square passband to filter the data around a certain date, plus a collection of sinusoidal functions, while for wavelet analysis, the wavelet functions (Eq. (1)) play the roles of both the filter and the (complex) periodic function.

The wavelet analysis was performed using the Python PyCWT spectral analysis module¹⁷ (based on Torrence & Compo 1998; Grinsted et al. 2004; Liu et al. 2007). The procedure requires the data to be uniformly sampled, which is the case for K2 observations. The few small gaps and inhomogeneities in the data were filled by interpolation. While this may have an effect on the shortest timescales sampled, it does not affect the final significant periods, which are in the range of days. The significance was estimated considering a model for red noise estimated through a Lag-1 autocorrelation of the original data (Torrence & Compo 1998; Liu et al. 2007). The results of the wavelet analysis are plotted in Fig. 6. The wavelet analysis recovers the results of the SGLSP and shows the same trends of quasi-periodicity, with significant signatures in the range of 2.4–2.6 d, 4.9–5.5 d, and, to a lesser extent, 9 d. As for the SGLSP analysis, not all the periods are recovered on all epochs and there is a drift in phase and periodicity that suggests changes in the structure that causes the eclipses.

To help visualize the periodicity, we phase-wrapped the data for the various potential periods. This exercise reveals a clear modulation for a period of 5.02 ± 0.12 d (corresponding to the main peak of the GLSP, see Fig. 7) and, to a lesser extent, for a period of half this value. Any other period either fails the visual check, or appears to be spurious (e.g., being an integer multiple of the sampling rate). The overall shape of the curve suggests a 5 d quasi-periodicity for the eclipses, and at the same time reveals a clear but rapidly-changing modulation in the flux observed out of eclipse. The tip of the phased lightcurve has an “M” shape that could be consistent with rotational modulation in a star with two non-identical cold spots (see Fig. 7, top right). Nevertheless, we also observe that the modulation suffers significant changes from one period to the next, which is unfeasible for typical long-lived, cold stellar spots. The eclipses are clearly associated to this 5 d period, although they also change in depth from period to period and show a drift in phase during the observed epochs (Fig. 7 bottom). The epochs of increased eclipse activity correspond to the times when the 2.5 d period is stronger. Having ruled out relatively stable structures (e.g., spots) for the variability, the observed lightcurve requires something that changes on the timescale of the rotational period, such as significant variations in the obscuring material in the innermost disk. The various possibilities are discussed in the following sections.

The REM data do not offer the same kind of time coverage and photometric stability, and thus the periodicity that can be inferred from them has low-significance. Moreover, the correlated errors in the REM magnitude are very complex, as the “redness” of the noise strongly depends on the observed cadence and the magnitude (given the typical cadence, it is unlikely to find many consecutive points on eclipse so that low magnitudes tend to be followed by quite uncorrelated ones, while “out-of-eclipse” points are often followed by a measurement with very similar value). This makes it very hard to assess the significance of the GLSP, and makes a wavelet analysis impossible. Nevertheless, the same rough behavior is observed, and wrapping the lightcurve reveals a dominant periodicity around 4.9 d, which appears correlated with the extinction events (see Fig. 8). There is no high-cadence periodic signature.

Fig. 3 Top: GLSP for complete K2 dataset. Bottom: GLSP forσ-clipped dataset (out-of-eclipse data). The significance levels are estimated according to a red noise model with correlation parameter α = 0.98, and the same uncertainty distribution and sampling as observed in the data (see text), so that the observed peaks are not significant.

Fig. 4 Histogram of K2 data by magnitude. Note that there is no clear separation between “on-eclipse” and “off-eclipse” parts, except in very deep eclipses (relative flux <0.75, green line). The red dotted line marks the separation from the σ-clipping filter. Selecting down to a different level (e.g., 0.983, orange line, where the distribution appears to slightly flatten out) does not introduce any significant change.

Fig. 5 Stacked GLSP for full K2 data (top three panels, stacking at intervals of ± 10 d, ± 20 d, ± 30 d) and σ-clipped data (out-of-eclipse data; bottom two panels, stacking at intervals of ±10 d and ± 20 d). The most significant periods detected in the individual GLSP are marked with vertical lines (see text for discussion). In each panel, x-axis shows periods, and y-axis shows date (in K2 mission days) around which we consider the time interval to estimate the SGLSP. Dates for which no data are available are left blank. Color scale is set so that purple is equivalent to 95% significance and dark blue is equivalent to 90% significance for a red noise model with the same number of points, similarly distributed (see text).

Fig. 6 Wavelet power spectrum for K2 data on RX J1604.3-2130. Left panel: wavelet spectrum in time. Color scale is set between minimum and maximum power, and global significance limits are marked as thick black (95% confidence) and thin black (90% confidence) contours, calculated for a red noise data model with α = 0.98. Regions where edge effects could be significant are masked out. Horizontal lines mark positions of the GLSP peaks found with K2 and REM (at 8.9 d, 5.5 d, 5.0 d, 4.9 d, 2.59 d, 2.39 d). Right panel: global wavelet power together with the significance, the noise spectrum, and the scaled GLSP for comparison.

Fig. 7 K2 data wrapped according to a period of 5.015 d. Data are color-coded according to date to better display points thatwere taken close in time. Upper left panel: entire K2 data wrapped, showing the phase twice to aid the eye. Upper right panel: zoom in on the “out-of-eclipse” part of the K2 data (colored points, selected via σ-clipping), together with the average and standard deviation in 10 intervals in phase space (black points; the errorbars correspond to the standard deviation within each bin) and an interpolation curve to trace the shape of the variations. Bottom panel: full K2 data plotted against mission day, together with interpolated modulation curve of top-right panel, to show how the phase of the eclipses drifts at certain times during the observation campaign, although 5 d periodicity remains visible throughout the full dataset.

3.2 REM data, extinction, and the inner-disk structure

The multi-band REM data allows us to study the color variation during the eclipses for the first time. The JHK data could provide a good insight into the properties of the obscuring matter, but the only dates for which we have complete data do not reveal substantial variability, and the only filter for which we cover a significant number of points and that shows variability is H. The optical data taken around the same dates as JK reveal that the object was essentially out of eclipse. The observed JHK colors are consistent with the colors of a pre-main-sequence star without significant near-IR excess (Kenyon & Hartmann 1995).

We thus concentrate on examining the optical and H-band data. Figure 9 shows the color variation as a function of the magnitude variation for several combinations of bands. Considering the standard interstellar extinction laws (R_V = 3.1–5.5, see, for instance, Schlegel et al. 1998; Stoughton et al. 2002; Cardelli et al. 1989), the slopes of the color variations up to Δg′ ~ 1 mag are quite consistent with the standard extinction vectors. There are some places where the color variation becomes suddenly smaller, for instance around Δg′ ~ 1 mag and (especially in the g′ vs. r′−z′ diagram), around Δg′ ~ 0.2 mag. These could be due to scattering shifting the colors towards bluer regions as the eclipse progresses (as has been observed in UXors, e.g. Grinin 1988). There are further changes in the slope as the eclipse progresses, with the deeper eclipses being better fit with an extinction law with higher R_V (which is, in general, attributed to high extinction clouds with more processed material; Cardelli et al. 1989), and the very deep eclipses having an offset with respect to the standard trend and showing a flatter reddening curve.

The color variation towards flatter reddening curves is only observed for eclipses with Δg′ ~ 0.9–1.1 mag. In those cases, the reddening continues into the deepest eclipses (for which Δg′ ≥ 0.9) with a slope that cannot be explained with the standard R_V = 3.1–5.5. The slope variation is likely related to strongly processed dust grains, although the change is different in different filters, maybe indicating an effect of grain size (Eiroa et al. 2002). Higher density or differences in dust properties and/or scattering within the clumpy material of the disk could also contribute.

There are also some shifts parallel to the extinction vectors visible at about all magnitudes. Small variations in the stellar luminosity (e.g. due to accretion) and/or scattering on longer timescales may also contribute to shifting some of the data taken on different dates parallel to the standard extinction vectors as seen in Fig. 9, even though small stellar luminosity and accretion changes are not expected to produce a significant change in the observed colors. It should be noted that a large fraction of the shallower eclipse points belong to the high-cadence observations, and thus were taken very close in time, so that their variation is smooth.

Leaving aside the UXor behavior towards bluer color, the changes in the color slope consistent with reddening could be caused by extinction by more opaque material in the deepest parts of the eclipse, and variations in the sizes of the grains, depending on the height in the innermost disk (e.g., McGinnis et al. 2015), due, for instance, to differential settling (e.g., D’Alessio et al. 2006; Laibe et al. 2014). The dust content around a star is constrained by the dust sublimation temperature (~1500 K, with some variations in the 1000–2000 K range depending on disk structure, density, and composition; Isella & Natta 2005; Kama et al. 2009), and typical observations can be well-fit with inner rims at T = 800–1200 K (McClure et al. 2013). Significant grain processing already happens at these (and much lower) temperatures (Tielens et al. 2005), so that the grains are likely different from plain ISM silicates. Grain growth is also generalized in disks, and larger grains tend to produce gray extinction with a lower color dependency at optical and near-IR bands (Miyake & Nakagawa 1993; Eiroa et al. 2002). Large grains are often involved in the best-fitting models for inner disk walls (McClure et al. 2013). In addition, there are other effects, such as a dusty wind (e.g., as observed in RW Aur; Bozhinova et al. 2016a), which could also cause obscurations. The Hα profile of RX J1604.3-2130A shows some blueshifted absorption (Manara et al., in prep., see also, Sect. 3.3), but since the accretion rate is low, accretion-related winds are expected to be weak and carry significantly less mass than the accretion flows.

The high-cadence data offers us a chance to explore what happens during a single eclipse on a timescale where the variations (accretion, luminosity) are likely negligible (Fig. 10). The data show a smooth transition from eclipse to maximum, although the high-cadence eclipse does not go as deep as those where a significant color offset is observed. The high-cadence data are fully consistent with a mild UXor behavior around Δg′ ~ 0.2 mag, as well as increasing dust extinction with standard extinction laws for thin ISM dust (R_V = 3.3). For the rest of the eclipses, including those happening at half of the period, the coverage is very scarce, but wrapping the data does not reveal a significant difference between the colors of different eclipses, nor between the colors of the full-period vs. half-period eclipses (see Fig. 8, bottom), other than the differences observed between shallow and deep eclipses.

We can then use the observed extinction to estimate the amount of material needed to produce the eclipses. The typical depth of the extinction event in g′ is 1.2 mag. Assuming standard interstellar extinction (Schlegel et al. 1998), this is equivalent to A_V = 1 mag, or N_H = 1.8e21 cm⁻² (Predehl & Schmitt 1995). For an average particle weight of 1.36 m_H (for solar metallicity neutral gas, e.g., see Mihalas 1978), we obtain 0.004 g cm⁻². It should benoted that because the metal vs. hydrogen content in the inner disk may be different (either higher or lower, e.g. Panić et al. 2009; Riviere-Marichalar et al. 2013), there is some uncertainty in this parameter. The stellar parameters suggest a dust destruction radius at 0.06–0.15 au. If this mass were distributed in a ring of radius 0.06 au (0.15 au) and at a height comparable to the stellar radius, the total mass associated with the structure would be about 2e-3 (5e-3) M_Ceres, including gas and dust. Nevertheless, we can derive a better estimate of the thickness of the structure from the observed shadows in the outer disk, which span about 20° in average (see sketch in Fig. 11). For an obscuring structure at 0.06 au (0.15 au) distance, this would mean a size of about 0.02 au (0.05 au) or 3 R_* (8 R_*) as a function of the stellar radius. The total as plus dust mass would be in the range of 1–6% of the mass of Ceres. This means that the innermost disk ring does not need to be very massive in order to reproduce the observed behavior. An extinction law for dark nebulae would result in a slightly larger mass, although for reasonable values, the main uncertainty in our estimate remains the uncertainty in the size of the innermost disk and the gas to dust ratio.

Pinilla et al. (2018) revealed flux variations in the scattered light J-band observations up to 0.4–0.6 with respect to the non-obscured flux. Standard extinction laws suggest A_g ~ 2.6–4.3 mag, about a factor of 1.4–2.3 higher than in the deepest optical eclipses, which would result in a similarly larger mass. There are several possibilities to explain this difference. Differences in alignment between the star-inner disk and the inner disk-outer disk could result in different column densities viewed on each line of sight (see Fig. 11). In addition, the fact that the shadow observations and the optical photometry are not simultaneous, and the dust content is known to be variable, plus the possibility that the dust does not follow a standard extinction law, especially in the deepest eclipses, may also play a role. Another possibility would be if the inner disk does not cover the full stellar disk along our line of sight, or if the disk does not generally cover the surface of the star, but a localized warp or blob on it does. In fact, the observed differences in extinction (observed from the lightcurve A_g ~ 0.4–1.8 mag, expected from shadows A_g ~ 2.6–4.3 mag) can be explained with a variable stellar disk coverage by the eclipsing material. Considering the density estimated from the extinction a variation in coverage ranging from 0% out of eclipse, to between 30% for shallow eclipses, down to 100% at the deepest ones, could explain the observed lightcurve. We highlight that the optical variability depends on various poorly-constrained parameters, such as the relation between infrared and optical extinction (i.e., dust properties), whether the disk is occulting the star towards the equator or towards the poles (due to limb darkening), and whether the star has additional causes of variability (e.g., hot and cold spots).

The observed small changes and displacements of the shadows are consistent with a “wavy” or warped disk (Pinilla et al. 2018), which is also consistent with other observations (Grinin et al. 2008; McGinnis et al. 2015), models of self-shadowed UXor disks (Dullemond et al. 2003), and the non-axisymmetric or clumpy structures observed in gas and dust in the inner disks of some young stars (e.g., Sicilia-Aguilar et al. 2012; Siwak et al. 2014; Scholz et al. 2019). A warp at the point where the disk and the star are connected at the basis of the accretion column (see, e.g., the models of Alencar et al. 2012, 2018; McGinnis et al. 2015; Bodman et al. 2017) could also explain the observations, as shown in the cartoon in Fig. 11. The variations of the position of the shadows (~20° in average; Pinilla et al. 2018), together with the estimated radius of the disk, suggest that the innermost disk is as wavy as it is thick. A wavy/warped disk that deviates from a flat structure in vertical scale by about a third of the disk radius (to account for the angle variation of the shadows) would also explain why the star is not always eclipsed, while the shadows on the outer disk are nearly ubiquitous, despite variability.

A 5 d orbital period corresponds roughly to the distance where the disk would have a temperature of around 1540 K (~0.06 au or 9.4 R_*). The precise location and temperature depend strongly on the dust properties, the structure of the inner disk rim, and the stellar parameters (especially the stellar mass derived from evolutionary tracks). But, in any case, the corotation radius is compatible with the dust sublimation radius. With this in mind, the quasi-periodicity of the eclipses discussed before suggests that the structure causing them is rapidly changing on timescales comparable to the Keplerian period of the material in the innermost disk (days), so even consecutive eclipses have different depths and lightcurve profiles. Disk precession may alter the inclination of the inner disk in time, but their timescales are typically much longer. For the general relativistic precession, the timescale is proportional to ${(c/v_{\text{K}})}^2$ times the orbital period, where c is the speed of light and v_K is the Keplerian orbital velocity. Kozai or secular resonances may also cause precession, but they are usually weaker than the relativistic effect (Kozai 1962; Ford et al. 2000) and require very massive and close-in companions that would have been spectroscopically detected.

The quasiperiodicity on short timescales and color variations suggest that the disk is highly asymmetric, with warps or clumps that are denser than the rest. It is also likely being externally modified/fed due to viscous matter transport, which can explain the sudden periods of intense dimming and concatenated eclipses, such as the one observed between days 2085 and 2105 in the K2 data (see Fig. 2). Taking into account the accretion rate ~3e-11 M_⊙ yr⁻¹ and assuming that the rate of transport of matter in the inner disk is similar, the mass of the innermost disk required by the eclipses is comparable to what accretion transport could provide in between two weeks’ and two months’ time. This means that the inner disk is filling up (and draining) on relatively short timescales, and therefore, significant changes could be observed on 5 d periods.

To summarize this section, the color variability confirms that the eclipses are consistent with extinction by dust with properties ranging from ISM dust to more processed grains, located in an irregular, warped disk at the corotation radius. The total mass depends on the dust properties and on the size of the disk, but ~1% M_Ceres of total mass (including gas and dust) is enough to explain the observations. Considering the accretion rate, it is not unexpected that the dust distribution in the innermost disk changes on short timescales, since a significant fraction of the obscuring matter in the innermost disk will be fed to the star (or drained from the innermost disk) on each rotational period, explaining the rapid variability.

Fig. 8 Top panels: REM data, wrapped for 4.9 d period to reveal the presence of periodic signatures. Bottom panels: REM colors, wrapped for 4.9 d period. All magnitudes are given relative to the median magnitudes in each filter.

Fig. 9 Relative color-magnitude variations in optical and H-band REM data. Dotted lines indicate the extinction vector for several extinction laws (RV = 3.3, 5.0 and 8.0, derived following Schlegel et al. 1998). High-cadence data are plotted in blue, while remaining data are shown in black. All magnitudes are relative to the median value in each filter. Black arrow on the left of the first panel shows the correlation direction for cases where g′ is also used to calculate the color, showing 3× the average uncertainty.

Fig. 10 Color variations in dip observed with high-cadence data. We see that the eclipse is not deep enough to show significant variations in the extinction law, although some of the UXor behavior is seen, especially in i′ − z′.

Fig. 11 Not-to-scale sketch showing inferred structure of the system. The star and its magnetosphere are displayed in the center, the red arrow indicates the rotation axis of the star, magnetosphere, and corotating inner disk. Inner disk is deformed in regions where the magnetosphere is joining in, which causes warps that cross our line of sight when the disk rotates. Theouter and inner disk are highly inclined in relation to each other, so there are always shadows along the line where the plane of the inner disk and the plane of the outer disk cross. It should be noted that the inner disk is expected to be more wobbly and irregular than shown here.

3.3 Rotation, accretion, and variability in the inner disk

Rotation rates suggested a very high inclination for the disk system around RX J1604.3-2130A (Davies 2019), despite its outer disk being nearly face-on (6°; Zhang et al. 2014; Dong et al. 2017). If the 5 d modulation observed was caused by rotation, we can use the observed rotational velocity (16.2 ± 0.6 km s⁻¹) to infer the inclination of the system. The maximum period that this velocity range allows is 4.2–4.5 d. A 5 d period would require a stellar radius ~10–15% larger than estimated, or larger uncertainties in the v sin i (e.g. due to differences in limb darkening). In any case, the star cannot be aligned with the outer disk, because it would rotate at breakup velocity, but it should be closer to equator-on, and thus rather aligned with the inner disk. If the 2.4–2.5 d period were a rotation period, the inclination angle would be about 38° (still far from aligned with the outer disk). Nevertheless, such a low angle would be inconsistent with the lower limit of 61° of Davies (2019), besides the fact that the shape of the lightcurve agrees with a 5 d rotational period in a star with two asymmetric eclipsing structures, rather than with a 2.5 d period.

Constructing a model for the out-of-eclipse lightcurve is complicated, because the sampling of the REM photometry is not high enough to show the detailed color variation during a single rotational period, and there are significant changes every few days (see Fig. 8). A cold spot model is not sufficient, since spot-related variability is usually of the order of 0.1 mag (Grankin et al. 2007) and can by no means explain the observed eclipses nor shadows, but it could explain the M-shaped part of the K2 data. The flux modulation at phase ~0.5 (see Fig. 7) can be obtained with spots 200–500 K cooler than the stellar photosphere and spot coverage between 0.03 and 0.36 (as in, e.g., Bozhinova et al. 2016b), all of them reasonable parameters, but the rapid variability from period to period is hard to explain. For instance, a smaller scale eclipse could also produce the same effect, and since the amount of matter responsible for this would be very small compared to what is observed to flow around the star from its accretion rate, rapid variations are more plausible. In addition, there is more evidence of the disk-star connection being dynamic and rapidly variable than in the case of stellar spots (e.g. Fonseca et al. 2014), and the color variations (Figs. 8, 9, and 10) suggest that occultations are the dominant process.

All the datasets (CSS, K2, and REM) show brief periods of increased eclipsing activity on timescales shorter than 5 d (usually, about 2.5 d). These could be times at which an increased accretion rate throughout the disk triggers accretion onto the other side of the star and feeds a secondary warp, causing additional extinction events. Modeling this situation would require, at a minimum, simultaneous high-cadence multi-color photometry and spectroscopy, and it is thus beyond the scope of this paper.

The archival HIRES data are consistent with the low accretion rate estimates from X-shooter spectra (Manara et al., in prep.), not showing the emission lines characteristic of strong accretors (Hamann & Persson 1992, for instance, no He I emission and Ca II IR narrow emission lines in the center of strong absorption components). There is also no evidence of spectroscopic binarity (see Appendix B). Hα and Hβ have strongly variable intensity and line profiles, with timescales of months to years. The features are better viewed after photospheric subtraction, using a photospheric template derived from the standards HD 114386 (which has been observed in the Hβ region) and HD 151541 (for which the availabledata cover the Hα region; see Fig. 12). After photospheric subtraction, Hβ emission becomes evident, together with evidence of both redshifted and (in Hα) blueshifted absorption components. Photospheric subtraction also reveal further complex absorption profiles in other lines, such as Na I D. Since there is no detailed day-to-day high-resolution spectroscopic followup, it is hard to assess whether the variations in the line profile are due to variations in the accretion rate or changes in the orientation of the accretion flows with respect to the observer. Nevertheless, the changes in equivalent width and 10% velocities suggest that the accretion rate is variable, with up to two orders of magnitude variability between the maximum and the minimum widths (using the 10% Hα width; Natta et al. 2004), with the 3e-11 M_⊙ yr⁻¹ value from Manara et al., in prep. being an intermediate rate. Accretion is very weak (consistent with no accretion) in the spectra from MJD 53 902.

For both the Hα and Hβ lines, the redshifted absorption features are dominant, and on two of the dates (MJD 54 642 and MJD 55 287) show a YY Ori or inverse PCygni profile (IPC) with an absorption component that goes below the continuum. This suggests that the accretion column must have been very close to viewed along the line of sight on these dates. Accreting along the line of sight could also mean that the gaseous material could be obscuring the star (as has been suggested in other objects, e.g., TW Hya; Siwak et al. 2014). Obscuration by dust could happen in a warp induced in the place where the accretion column is attached to the disk (McGinnis et al. 2015; Alencar et al. 2018), although unfortunately none of the spectra has simultaneous photometry. The spectrum taken on MJD 54 642, which shows a mild IPC profile in Hα, is the one closest to a photometric datapoint from CSS corresponding to a deep eclipse. However, the HIRES data were taken 24 h after the photometry, and since typical eclipses last less than 24 h, it is not possible to assume that the star would have been in eclipse. Moreover, phase shifts are expected between the very-close-in gas emitting Hα and eclipses associated to dust at the corotation radius. On MJD 55 287, the IPC is clearly visible in Hα, and very strong in Hβ, a typical signature of a highly inclined system (Alencar et al. 2012, 2018; Donati et al. 2019). This, together with the observed eclipses, is a sign that RX J1604.3-2130A may be very similar to other highly inclined systems such as Lk Ca15, AA Tau, and V354 Mon (Bouvier et al. 2007; Alencar et al. 2018; Donati et al. 2019; Fonseca et al. 2014).

Since the analysis of the extinction suggests an inner disk that is routinely drained on a short timescale due to accretion, we explore whether the rate at which the magnitude varies during the eclipse is consistent with material transported due to accretion. A change in magnitude vs. time is equivalent to a change in column density over time, using the conversion between extinction and column density, as in Sect. 3.2. Assuming that this extinction event covers the stellar surface uniformly, the obscuring mass involved per time can be obtained by multiplying by the area of the stellar disk. This value can then be transformed into a approximate “projected” accretion rate that can be compared to the accretion rate measured by spectroscopy. Because the REM data has only very scarce coverage of each dip, we need to do the exercise with the K2 data, although one of the main limitations is the lack of color information.

For K2, we calculate the change in magnitude between each two points, i and i + 1, using the flux ratio between these two points to estimate the variation in magnitude as − 2.5 log₁₀(f_i+1∕f_i). Dividing by the time interval, we can obtain the change in magnitude vs. time (Fig. 13) and transform it into an approximate accretion rate as explained above. Using the conversion between K2 extinction and A_V (A_V = 0.4 A_K2¹⁸; Rodrigues et al. 2014), we obtain a typical change of 1 mag day⁻¹. With the relation between A_V and column density (Predehl & Schmitt 1995), we derive an approximate accretion rate of 6e-12 M_⊙ yr⁻¹. The large majority (99.1%) of the observed points fall between the (−5, +5) mag day⁻¹ interval, and these would correspond to accretion rates up to 3e-11 M_⊙ yr⁻¹, fully consistent with the accretion rate estimates from X-shooter observations (Manara et al., in prep.). Note that these “changes” do not need to be accretion rate variations, since we are only taking into account matter along the line of sight, and that the estimates are also subject to the same uncertainties as the inner disk mass estimates, namely dust properties and gas to dust fraction. It is also important to keep in mind that the obscuring matter contains dust, and must then be located at least at the dust destruction radius near the star-disk connection, while what is observed in Hα corresponds to hot gas likely closer to the stellar photosphere. Because of this, changes in the dust are not expected to be observed as immediate changes in Hα, and a phase delay is very likely to occur, as it is also observed between accretion-related emission lines with different energetics (Sicilia-Aguilar et al. 2015), or between veiling and line emission, or optical and X-ray accretion signatures (Dupree et al. 2012).

Fig. 12 Photosphere-subtracted Hα (top) and Hβ (bottom) spectra observed with HIRES/Keck (for various epochs) and X-shooter (corresponding to the time when the accretion rate and the stellar properties were estimated; Manara et al., in prep.). Zero levels for flux and radial velocity are marked by dotted lines. A rotationally-broadened template photosphere has been subtracted to show absorption and emission features.

Fig. 13 Time variation for K2 observations. Upper left panel: change in magnitude per day as a function of the observing date. Gray lines mark the cases in which the K2 relative flux drops below 0.94 (the limiting value for eclipses considered, see text). Upper right panel: histogram of the same plot. Lower panels: zoom-in of the magnitude change vs. time for several of the eclipses.

Fig. 14 WISE lightcurve during the two periods of observation. The W4 magnitude, significantly brighter than the rest, is shifted by2 mags for easier viewing.

3.4 The long-term variability of the inner disk

The disk around RX J1604.3-2130A has been classified as transitional (Carpenter et al. 2006; Mathews et al. 2012; Esplin et al. 2018). Luhman & Mamajek (2012) pointed out that the mid-IR observations with Spitzer/IRAC were consistent with a bare photosphere, while WISE data revealed a clear IR excess. Considering the fluxes provided by Carpenter et al. (2006) and the IRAC zeropoints¹⁹, we find that the source changed by about 1.47 mag at 4.5 μm between the Spitzer and the WISE observations²⁰, which are roughly separated by four years in time (from MJD 53 820 to MJD 55 249). The AllWISE lightcurve, on the other hand, reveals only mild mid-IR variability (as expected if the cause is extinction) during the half year interval, during which RX J1604.3-2130A was observed by WISE (see Fig. 14) with amplitudes of 0.17 mag (W1), 0.13 mag (W2), 0.07 mag (W3), and 0.24 mag (W4). It is interesting that while W1, W2, and W3 follow roughly the same pattern and vary in parallel, W4 behavesdifferently, although the uncertainties are also larger. At longer wavelengths, the extinction becomes negligible, so large changes in the flux are most likely related to changes in the disk structure. Wavelengths around 22 μm trace material at considerably larger distances, and could be dominated by the emission of the outer disk, less variable on short timescales. Note that the dramatic (>1 mag) IR variability affects only the IRAC bands, since W4 and MIPS 24 μm roughly agree despite being observed at times where there is substantial difference for the 3–10 μm fluxes. Therefore, although the strong variability in the near-IR is reminiscent of the “seesaw” behavior observed in some disks (e.g., Espaillat et al. 2011; Flaherty et al. 2012), we note that the situation here is quite different, because the mid-IR fluxes are not variable. While seesaw behavior in (usually, pre-transitional) disks can be explained by changes in the vertical scale without much modification to the contents of the disk, here, the observed near-IR photospheric colors need a very strong depletion of warm dust.

The CSS data reveal decreased eclipse activity soon after the Spitzer observations were obtained, which points towards the intriguing possibility that the inner disk may have been depleted of dust around that date. The CSS data closer in time to the Spitzer observations (56 datapoints in total within 10–20 days) do not show any eclipse, although there are no observations for nearly half a year before the Spitzer observations, and the first data point was taken more than 40 days after Spitzer. The sparse sampling may have missed the eclipses, since the CSS observations were obtained at irregular intervals of a few days during 2 months. However, with the sampling frequency and given that altogether, 35% of the CSS observations appear to have been taken during eclipses, it is significant that not a single eclipse would have occurred near the time of the Spitzer observations, suggesting that the eclipse activity may have indeed been lower at the time. In addition, the HIRES spectra taken during the CSS stable phase after the Spitzer observations show a significantly weaker Hα feature, consistent with decreased accretion, and maybe an additional signature of a lack of material in the innermost disk.

The strong Spitzer vs. WISE mid-IR variability is very similar to what has been observed for GW Ori, a triple stellar system with a circumtriple disk and an unstable or “leaky” dust filter that leads to remarkable changes in the SED due to variations in the innermost disk dust content on timescales of decades (Fang et al. 2014). Similarly to GW Ori, we can make an order of magnitude estimate of the amount of dust needed in the inner disk to account for the mid-IR variability. This can be compared with the accretion rate and with the feeding rate needed to support the variability of the innermost disk structure deduced from the extinction events, assuming matter flows through the disk in a stable way. Nevertheless, studying the SED of RX J1604.3-2130A is more complicated, because the data are non-simultaneous, and, unlike GW Ori, the star is highly variable in the optical. To constrain the disk properties, we need to disentangle what is caused by extinction from what may be caused from variations in the innermost disk structure. We construct the SED using VizieR multiwavelength data (see Appendix C) plus our REM photometry. We consider the brightest magnitudes observed in each optical filter as stellar photosphere (out-of-eclipse data). The out-of-eclipse data are corrected for A_V = 1 mag (Preibisch & Zinnecker 1999) using standard color relations (Cardelli et al. 1989; Stoughton et al. 2002; Bessell 2005). The data that are fainter than the out-of-eclipse values are not corrected for extinction and are not used for the luminosity fit. The SED reveals that the JHK archival data were likely taken during one of the obscuration episodes (see Fig. 15), so our REM observations are brighter than previous ones and more consistent with a K3 photosphere. The CSS observations suggest that the star was in a bright state near the Spitzer observations, and the WISE observations do not seem seriously affected by the strong optical variability at the time. Thus it is reasonable to assume that the main changes at wavelengths >4 μm are caused by variations in the disk structure, and thus a model can be constructed assuming the disk illuminated by the non-extincted stellar photosphere.

A very simple inner disk model (not taking the outer disk into account) can be constructed to estimate a lower-limit to the disk mass. Assuming that the innermost disk is optically thin in the IR and dominated by a single temperature, the flux observed at a given frequency ν, F_ν, can be written as $$F_{\nu }=\Omega B_{\nu }(T_{\text{ring}})k_{\nu }\Sigma ,$$(2)

where Ω is the solid angle subtended by the structure, B_ν(T_ring) is the black body function evaluated for the dominant ring temperature, T_ring, k_ν is the dust mass absorption coefficient (cm² g⁻¹), and Σ is the surface density of the ring. Since Ω depends on the area of the disk (with an inclination factor included), this means that F_ν is proportional to the mass of the ring. Taking k_ν to be a power law of the frequency with exponent β = 1.9 and k_{70 μm} = 118 cm² g⁻¹ (for small interstellar grains without ice mantles; Roccatagliata et al. 2013), we find that the mid-IR emission can be relatively well fit with a dominant ring temperature of 850 K, and a dust mass of 3e-3 M_Ceres, or 4e-4 M_Ceres for a higher temperature of 1270 K. Considering the corotation temperature, 1540 K, the amount of dust required would decrease to 2e-4 M_Ceres. If the gas-to-dust ratio takes the usual value of 100, this would mean a mass between 2–30% M_Ceres, including gas and dust. The result is strongly dependent on the ring temperature, and because a single temperature cannot fit all the observed datapoints equally well, we expect the inner disk to cover a range of temperatures. Further variations in the dust properties (highly processed silicates or large grains have a different emissivity), and whether the disk is optically thin (the mass is a lower limit) also play a role here. Nevertheless, this exercise reveals that, as for GW Ori (Fang et al. 2014), a small change in the dust content of the inner disk can be responsible for a very remarkable change in the mid-IR fluxes. For the observed accretion rate, this amount of mass is comparable to what can be accreted in a few months to a few years, so piling up the extra material in ~4 yr time could be done with a mismatch between disk accretion and stellar accretion happening at some place between the outer disk and the stellar magnetosphere.

Since the assumptions of having an optically thin inner disk and constant temperature are a poor approximation for a relatively dense and likely extended protoplanetary disk, we also modeled the disk using the radiative transfer code, RADMC-2D (Dullemond & Dominik 2004; Dullemond 2011). Because the SED is so uncertain in terms of extinction and variability, we just aim to estimate how much mass would be needed to reproduce the variability in the innermost disk, and the misinclination of the disks is not taken into account. As a photospheric model, we use a MARCS photosphere (Gustafsson et al. 2008) scaled to the observed luminosity. We explored various dust distributions within the nearly-cleared inner hole of the disk, and set up the scattered light ring as the maximum outer disk inner radius. The total mass needed depends on the grain properties. We assume standard amorphous silicates and carbon (in a proportion 3:1) with typical sizes between 0.1 μm to 1cm (following a collisional distribution with exponent −3.5), with opacities derived from the Jena database (Jäger et al. 2003). Including a dust surface density of the order of 1e-8 g cm⁻² at 0.06 au, decreasing with the disk radius as r⁻² or steeper (so that most of the mass is concentrated in the innermost few au), is enough to go from photospheric fluxes to the observed ones, and the total mass involved would be of the order of 3–5×10⁻⁵ M_Ceres if we assume agas-to-dust ratio of 100. It should be noted that the RADMC model underestimates the near-IR fluxes, so an inner wall (not included in the current model), and thus a larger mass, are likely. As a final point, although the far-IR part of the SED is highly uncertain due to the large beam and unknown variability, these simple models also reveal that the 22–24 μm flux is dominated by the outer disk inner wall, which explains the lack of variability between W4 and MIPS 24 μm.

The model confirms that the innermost disk structural changes required to alter the IR SED appearance as observed could be caused by small accretion variations between the outer and the inner disk, without needing any more dramatic processes. Although the SED can be fitted with an inner disk that is clearly devoid of matter compared to the outer disk, there has to be enough material close to the star to support the observed accretion, although the present data does not allow us to set further constraints.

From this exercise, despite all uncertainties involved in the analysis of the non-simultaneous SED of a highly variable object, the conclusion we reach is that the variability observed in the innermost disk is consistent with matter transport at rates similar to the observed accretion rate onto the star. The timescales of the variations (~4 yr) imply that the material cannot be located much further away than a couple of au from the star. The mass in the inner disk would be in any case many orders of magnitude lower than the disk mass, which is estimated to be 0.018 M_⊙ in the RADMC models to account for the 880 μm flux. The time constraint imposed by the observed variability would be consistent with the orbital time of an undetected companion at a couple of au. Dusty material moving at few au should definitely emit in the mid-IR, so that future time-resolved mid-IR data over several years may help to explore this scenario.

Fig. 15 SED with archival data for RX J1604.3-2130A and simple disk models. The open dots show all the available data obtained from VizieR and not corrected for extinction. REM data are shown for the filters with negligible color terms. Blue dots show the selected brightest magnitudes for each band, corrected by the nominal AV = 1.0 mag (Preibisch & Zinnecker 1999). Three very simple models are shown for comparison for the purpose of estimating the mass that could be associated with inner disk emission (see text for details): two constructed using RADMC radiative transfer code, and two assume a grey body model for material located in a ring with temperature 850 K and dust mass 0.003 MCeres, a temperature 1270 K, dust mass 4e-4 MCeres, and a temperature of 1540 K with dust mass 2e-4 MCeres, respectively.

4 Discussion

We now use all the previous information to trace a complete picture of the RX J1604.3-2130A system and to investigate the physical mechanism(s) behind the variability observed. The lightcurve presents dramatic changes on very short timescales, including phase shifts, lightcurve profile variations from period to period, and short-lived period drifts, even if the system eventually reverts to the 5 d period. Such changes are too rapid to be explained by long-lived, cold stellar spots. Moreover, the dusty composition of the occulting material requires the structures to be located beyond the dust destruction radius, for instance, in a clumpy or warped disk. The eclipsing material would be located at the inner disk at the corotation radius, where star-disk interactions are expected to be highly dynamical (Fonseca et al. 2014; McGinnis et al. 2015), although current stellar parameters can be only reconciled with a rotational period if the radius is increased by ~10–15%. The observations are consistent with one major and one minor warps on a dusty disk that is highly inclined with respect to the outer disk. The warps may be associated to quasi-stable accretion columns (as in Alencar et al. 2012, 2018; Bodman et al. 2017). This highly inclined, wobbly/warped disk can also explain the shadows observed in the outer disk (Pinilla et al. 2018). The dust observations are in agreement with ALMA gas observations, suggesting that the innermost gaseous disk is also misaligned with respect to the outer disk (Mayama et al. 2018).

Magneto hydro-dynamic (MHD) models usually distinguish between two accretion scenarios: stable (with accretion columns that are fixed to the star for at least several rotational periods) and unstable (where accretion proceeds through several rapidly-changing fingers due to e.g., Rayleigh–Taylor instability between the innermost disk and the star, see Kurosawa & Romanova 2013). Unstable accretion can explain quick changes from one rotational period to the next. Nevertheless, although we observe a drift in phase along the K2 observations, the timing of the eclipses is not chaotic (Fig. 7), but rather consistent with two structures on either side of the star. These structures do change from one rotation period to the next, but they are more stable than Rayleigh–Taylor fingers distributed over the stellar surface. The well-defined redshifted absorption features in some of the Hα and Hβ spectra are also characteristic of a system with stable accretion columns viewed close to the line of sight. Thus, one possibility would be RX J1604.3-2130A having intermediate characteristics between stable and unstable accretion, such as relatively stable accretion structures attached to a particular longitude and latitude on the stellar surface (e.g., due to a localized, dipolar magnetic field) but locally unstable (e.g., due to Rayleigh–Taylor instability as proposed by Kurosawa & Romanova 2013, or any other localized instability). If the accretion columns are relatively stable and locked to the disk at corotation, the most intense eclipses would be related to the rotational period. Unstable accretion along an otherwise-well-defined column could produce changes from one rotational period to the next and changes in the vertical structure of the inner disk.

The period of increased eclipsing activity (from approximately mission day 2085 to 2105) is consistent with the picture of relatively unstable accretion through two well-defined structures. The periodicity is dominated by the 2.5 d signature during this epoch (see Figs. 5 and 6). This suggests that a change in the inner disk mass (deeper and more frequent eclipses) results in accretion instability on the side of the star that is usually quiescent. Triggering accretion by Rayleigh–Taylor instability does not need a change in the stellar magnetic field, but could result from a change on the viscosity and/or disk accretion rate (Kurosawa & Romanova 2013). Additional material flowing inwards from the outer disk could increase the mass and density of the inner disk from time to time. Some of the eclipses observed during this time are particularly deep, which would be also in agreement with an increased amount of matter in the inner disk.

Irregular feeding of the innermost disk could therefore be a way to trigger both the 2.5 d eclipses and to keep the star in an unstable state, where Rayleigh–Taylor instability changes the properties of the accretion columns on timescales of days. It is likely that the accretion rate would be higher if the star is accreting through both sides. Nevertheless, an increase in accretion by a factor of few in a star that does not have a particularly high accretion rate (so the accretion luminosity is small compared to the stellar luminosity) and that has complex extinction is hard to measure, besides, the two spots cannot be observed simultaneously in a star that is nearly equator-on. Detailed time-resolved spectroscopy could help to disentangle rotational modulations from the effect of increased accretion (Sicilia-Aguilar et al. 2015).

Most of the line profile variations observed with HIRES can be reproduced by rotational modulations, although both the Hα variability and the periods of increased eclipsing activity suggest that the accretion rate is variable on longer timescales. Moreover, the times at which the star was observed to have photospheric mid-IR colors are also coincidental with the times at which Hα is weakest, suggesting lower accretion in agreement with the picture that the time-variability observed in the inner disk depends on the mass transport on a larger radial scale. Future observations of similar events are required for confirmation.

Massive bodies in the inner disk could be responsible for small variations in the matter flow that could increase the pressure over the stellar magnetosphere and promote accretion along the two magnetically-active regions on both sides of the star. A highly asymmetric magnetic field distributed over the surface of the star may be also taking part of the regulation so that accretion only proceeds through intense magnetic field regions where the star can efficiently lock into the innermost disk, not necessarily matching the rate of accretion through the disk and onto the star. If the corotation radius is located at the distance where the Keplerian period is 5 d, the magnetosphere will be quite large (~9–10 R_*) compared to what is usually assumed for young stars. This may be an additional signature of a strong magnetic field, maybe similar to LkCa 15 (Donati et al. 2019).

Despite the rough agreement between the accretion rate needed to transport matter through the inner disk and the mass required to produce the shadows and eclipses, an important problem remains: dealing with the angular momentum transport when matter from a face-on outer disk is transferred to a nearly edge-on inner disk. Companions between the outer and inner disk may help in the process, but accretion of angular momentum from the outer disk would also tend to change the orientation of the very low-mass inner disk. A strong stellar magnetic field ensuring disk locking may force the material in a particular orientation, but probably only very close-in. Therefore, the structure of the inner disk and of the material in the cavity may be tridimensional, highly complex, and highly dynamic, including bridges and streamers, as has been suggested for another dipper, AA Tau (Loomis et al. 2017). Near- and mid-IR observations following the object for several years could help to constrain the structure and matter flow within the cavity. Due to the location of the outer disk, a companion aligned with the outer disk could help to explain the formation history of the system, but it would not help with the change in angular momentum for material that is transported from the outer to inner disk. Nevertheless, massive companions (>2–3 M_J) at >22 au have been ruled out (Kraus et al. 2008; Canovas et al. 2017), which may be a problem to explain matter and angular momentum transfer from the outer to the inner disk. RX J1604.3-2130B appears to be located too far away from the disk (Köhler et al. 2000) to have a significant effect, unless it has a highly eccentric orbit. The Gaia DR2 parallax and proper motions for RX J1604.3-2130A (ϖ = 6.662 ± 0.057 mas, −12.33 ± 0.10 mas yr⁻¹ and −23.83 ± 0.05 mas yr⁻¹) and RX J1604.3-2130B (ϖ = 6.79 ± 0.10, −12.64 ± 0.18 mas yr⁻¹ and −24.73 ± 0.09 mas yr⁻¹) are consistent with a common origin, but do not constrain the orbital properties. A complex formation history, with two different protostellar collapse and accretion episodes (one originating the inner disk plus the star, and the second producing the outer disk), could explain the formation. Studying such formation scenarios, as well as the transfer of matter from the outer to the inner disk, would require a better knowledge of the detailed disk structure and companions, in addition to hydrodynamical simulations.

5 Summary and conclusions

This study reveals the power of time-resolve data together with years of archival multi-wavelength data in the task of disentangling the properties of very complex systems such as RX J1603.4-2130A. Our main results are summarized below.

The observed eclipses have a quasi-periodicity of 5 d, consistent with the rotational period of the star, and can be explained as extinction by dusty structures, probably warps at the points where the accretion columns leave the disk. The eclipsing dust is located at the dust destruction radius in corotation with the star. The corotation radius is of the order of 9–10 R_*, suggesting a very extended magnetosphere. The location of the warps suggests that the star is accreting on a relatively stable regime, although occasional instabilities, maybe induced by variations in the amount of matter in the inner disk, may trigger enhanced accretion through a secondary spot, also resulting in more frequent eclipses.

The amount of material required to produce the optical eclipses (to the order of few percent of M_Ceres, including gas and dust) is comparable to what is brought in by accretion in a few days to a few weeks (depending on the inner disk structure and gas to dust ratio, and on the variable accretion rate), which explains the rapid variability. Leaving aside the uncertainties in dust properties, structure shape, and inclination, the amount of material needed to produce the IR shadows is also roughly consistent with the mass needed to reproduce the optical eclipses and rapid variability.

As observed for GW Ori (Fang et al. 2014), the dusty (and likely gaseous) contents of the inner disk are variable in timescales of years. These long-term variations could be produced by small mismatches of the accretion rate between the outer and the inner disk, and can explain the lack of eclipsing activity that coincides with the lack of near-IR excess observed in the past. Future simultaneous optical and mid-IR data would be required to confirm the transport of material throughout the inner disk. In particular, simultaneous time-resolved photometry and spectroscopy may help us to pin down potential phase shifts between Hα (close to the star) and the eclipses (close to the dust destruction radius) that can reveal the innermost disk and magnetosphere structure. Longer term followup (years) in the optical and IR could be used to understand transport from the outer to the inner disk.

Together, the time-resolved photometry and spectroscopy also confirm that the star is not aligned with the outer disk, but highly inclined and rather aligned with the low-mass, highly variable, inner disk. Systems with highly inclined disks with respect to the stellar rotation also pose a problem to the standard picture of protostellar collapse and disk formation mechanisms, so a deeper understanding of them may help us to understand the angular momentum transfer and the effect of initial conditions at protostar or cluster levels on the formation, and subsequent evolution, of protoplanetary disks.

Acknowledgements

We thank the referee for the comments and suggestions that helped to clarify this paper, and S. Matsumura,P. Palmer, and M.H. Kristiansen for useful comments and discussion. A.S.A. is partly supported by the STFC grant number ST/S000399/1 (The Planet-Disk Connection: Accretion, Disk Structure, and Planet Formation). C.F.M. acknowledges an ESO fellowship and has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 823 823 (RISE DUSTBUSTERS). C.F.M. was partly supported by the Deutsche Forschungs-Gemeinschaft (DFG, German Research Foundation) – Ref no. FOR 2634/1 TE 1024/1-1. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming disks). P.P. acknowledges support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. J.d.B. acknowledges funding from the European Research Council under ERC Starting Grant agreement 678 194 (FALCONER). J.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742 095; SPIDI: Star-Planets-Inner Disk-Interactions. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. This research hasmade use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research hasmade use of the VizieR catalogue access tool, CDS, Strasbourg, France (Ochsenbein et al. 2000). This work makes use of the CSS survey, which is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the US National Science Foundation under grants AST-0909182 and AST-1 313 422. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration This research made use of Astropy²¹ a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), and PyAstronomy²².

Appendix A Flares observed in the K2 data

Despite the mild accretion, there are signatures of activity. There are several brief increases in brightness in the K2 data that could be stellar flares (see Table A.1). All the cases were checked against contamination from Solar System objects, which can mimic the flaring behavior (Rappaport et al. 2019), using SkyBot (Berthier et al. 2006, 2016) as done by Szabó et al. (2015), and no contamination was detected (Martti Holst Kristiansen, priv. comm.). In five cases, the increased flux is detected in at least three datapoints, revealing a profile consistent with stellar flares (Doyle et al. 2018; Gershberg & Shakhovskaia 1983), with a sharp rise and slower decline and duration between 1 and 3 hr (see Fig. A.1). Four more potential events are detected in fewer than three observations (consistent with short flares, Doyle et al. 2018, but also harder to confirm) and/or are weaker. One among them happens during an eclipse, so their classification is more uncertain. One of them (on MJD 59 021.611483) has extreme intensity, although it does not show the usual flare profile (rapid rise, slower decay), but rises over two observations and falls of from the peak within half an hour. The two largest flares happen at the beginning and at the end of the period between 2091 and 2010 d, when the star also presents concatenated eclipses, but more data would be needed to tell whether this is significant. One of the flaring events may be composed of at least two consecutive flares, as has been observed in some M-type stars (Doyle et al. 2018). No flares are detected in the REM data, likely due to the sparse sampling.

Fig. A.1 Observed flares.

Appendix B Summary of the kinematics from optical spectroscopy

Table B.1 lists the rotational and radial velocities derived for all the individual Keck spectra. The rotational and radial velocities are consistent throughout all the spectra, although the quality of the data from MJD 55 287.614 are better than the others. There is no evidence of radial velocity variability within the observed errors. Constraining the possible radial velocity companions with sparse-sampled data in an object with variable accretion and quasi-stable accretion columns along the line of sight would require a detailed analysis to distinguish accretion and gas absorption from the disk from those of potential companions (e.g. see Mora et al. 2002; Sicilia-Aguilar et al. 2015), which is beyond the scope of this paper. Note that the spectra from 54 689.304 and 55 256.556 are noisier than the rest in the regions used for the velocity estimates, and thus the results are highly uncertain.

Appendix C Archival SED photometry

Table C.1 contains the archival photometry from VizieR that we selected as being most likely representative for the out-of-eclipse SED. We chose the brightest magnitudes observed in each filter, adding also those that are unlikely to change due to extinction and on short timescales (mid-IR and longer wavelengths). Due to the non-homogeneity of the sample and the lack of detailed knowledge of the full spectrum variability, some of the selected points may have been taken during eclipse phases.

References

All Tables

All Figures

	Fig. 1 REM lightcurve for RX J1604.3-2130A. Upper panel: full lightcurve observed, while lower panel: zoom in the high-cadence data. All magnitudes shown are relative magnitudes evaluated with respect to the σ-clipped average magnitude of all datapoints in each band, marked as zero point with a dotted line.
In the text

	Fig. 2 Top: K2 lightcurve for RX J1604.3-2130. Typical uncertainties are smaller than dots. Bottom: Catalina Survey DR2 lightcurve. Since Catalina data cover the epochs of the Spitzer, WISE, and HIRES spectroscopy observations, we have marked them in the figure as vertical lines (see text). In both cases, the average flux (magnitude) is shown as a dotted line.
In the text

	Fig. 3 Top: GLSP for complete K2 dataset. Bottom: GLSP forσ-clipped dataset (out-of-eclipse data). The significance levels are estimated according to a red noise model with correlation parameter α = 0.98, and the same uncertainty distribution and sampling as observed in the data (see text), so that the observed peaks are not significant.
In the text

Fig. 4 Histogram of K2 data by magnitude. Note that there is no clear separation between “on-eclipse” and “off-eclipse” parts, except in very deep eclipses (relative flux <0.75, green line). The red dotted line marks the separation from the σ-clipping filter. Selecting down to a different level (e.g., 0.983, orange line, where the distribution appears to slightly flatten out) does not introduce any significant change.

In the text

Fig. 5 Stacked GLSP for full K2 data (top three panels, stacking at intervals of ± 10 d, ± 20 d, ± 30 d) and σ-clipped data (out-of-eclipse data; bottom two panels, stacking at intervals of ±10 d and ± 20 d). The most significant periods detected in the individual GLSP are marked with vertical lines (see text for discussion). In each panel, x-axis shows periods, and y-axis shows date (in K2 mission days) around which we consider the time interval to estimate the SGLSP. Dates for which no data are available are left blank. Color scale is set so that purple is equivalent to 95% significance and dark blue is equivalent to 90% significance for a red noise model with the same number of points, similarly distributed (see text).

In the text

Fig. 6 Wavelet power spectrum for K2 data on RX J1604.3-2130. Left panel: wavelet spectrum in time. Color scale is set between minimum and maximum power, and global significance limits are marked as thick black (95% confidence) and thin black (90% confidence) contours, calculated for a red noise data model with α = 0.98. Regions where edge effects could be significant are masked out. Horizontal lines mark positions of the GLSP peaks found with K2 and REM (at 8.9 d, 5.5 d, 5.0 d, 4.9 d, 2.59 d, 2.39 d). Right panel: global wavelet power together with the significance, the noise spectrum, and the scaled GLSP for comparison.

In the text

Fig. 7 K2 data wrapped according to a period of 5.015 d. Data are color-coded according to date to better display points thatwere taken close in time. Upper left panel: entire K2 data wrapped, showing the phase twice to aid the eye. Upper right panel: zoom in on the “out-of-eclipse” part of the K2 data (colored points, selected via σ-clipping), together with the average and standard deviation in 10 intervals in phase space (black points; the errorbars correspond to the standard deviation within each bin) and an interpolation curve to trace the shape of the variations. Bottom panel: full K2 data plotted against mission day, together with interpolated modulation curve of top-right panel, to show how the phase of the eclipses drifts at certain times during the observation campaign, although 5 d periodicity remains visible throughout the full dataset.

In the text

	Fig. 8 Top panels: REM data, wrapped for 4.9 d period to reveal the presence of periodic signatures. Bottom panels: REM colors, wrapped for 4.9 d period. All magnitudes are given relative to the median magnitudes in each filter.
In the text

Fig. 9 Relative color-magnitude variations in optical and H-band REM data. Dotted lines indicate the extinction vector for several extinction laws (RV = 3.3, 5.0 and 8.0, derived following Schlegel et al. 1998). High-cadence data are plotted in blue, while remaining data are shown in black. All magnitudes are relative to the median value in each filter. Black arrow on the left of the first panel shows the correlation direction for cases where g′ is also used to calculate the color, showing 3× the average uncertainty.

In the text

	Fig. 10 Color variations in dip observed with high-cadence data. We see that the eclipse is not deep enough to show significant variations in the extinction law, although some of the UXor behavior is seen, especially in i′ − z′.
In the text

Fig. 11 Not-to-scale sketch showing inferred structure of the system. The star and its magnetosphere are displayed in the center, the red arrow indicates the rotation axis of the star, magnetosphere, and corotating inner disk. Inner disk is deformed in regions where the magnetosphere is joining in, which causes warps that cross our line of sight when the disk rotates. Theouter and inner disk are highly inclined in relation to each other, so there are always shadows along the line where the plane of the inner disk and the plane of the outer disk cross. It should be noted that the inner disk is expected to be more wobbly and irregular than shown here.

In the text

Fig. 12 Photosphere-subtracted Hα (top) and Hβ (bottom) spectra observed with HIRES/Keck (for various epochs) and X-shooter (corresponding to the time when the accretion rate and the stellar properties were estimated; Manara et al., in prep.). Zero levels for flux and radial velocity are marked by dotted lines. A rotationally-broadened template photosphere has been subtracted to show absorption and emission features.

In the text

	Fig. 13 Time variation for K2 observations. Upper left panel: change in magnitude per day as a function of the observing date. Gray lines mark the cases in which the K2 relative flux drops below 0.94 (the limiting value for eclipses considered, see text). Upper right panel: histogram of the same plot. Lower panels: zoom-in of the magnitude change vs. time for several of the eclipses.
In the text

	Fig. 14 WISE lightcurve during the two periods of observation. The W4 magnitude, significantly brighter than the rest, is shifted by2 mags for easier viewing.
In the text

Fig. 15 SED with archival data for RX J1604.3-2130A and simple disk models. The open dots show all the available data obtained from VizieR and not corrected for extinction. REM data are shown for the filters with negligible color terms. Blue dots show the selected brightest magnitudes for each band, corrected by the nominal AV = 1.0 mag (Preibisch & Zinnecker 1999). Three very simple models are shown for comparison for the purpose of estimating the mass that could be associated with inner disk emission (see text for details): two constructed using RADMC radiative transfer code, and two assume a grey body model for material located in a ring with temperature 850 K and dust mass 0.003 MCeres, a temperature 1270 K, dust mass 4e-4 MCeres, and a temperature of 1540 K with dust mass 2e-4 MCeres, respectively.

In the text

	Fig. A.1 Observed flares.
In the text