© The Authors 2025

1 Introduction

Protoplanetary discs, known to be the nurseries of planets, play a crucial role in determining the formation conditions and early evolution of planetary systems. The creation of close-in planets, similar to our own, is particularly influenced by the innermost regions of the disc (< 1 AU) and so it is here that valuable insights can be gained into the formation process (Lissauer et al. 2011). It is also in this region that material is accreted onto the central star (e.g. Hartmann et al. 2016) and ejected through disc winds and outflows (Frank et al. 2014), processes that play an important role in the evolution of the star and its circumstellar disc. Observational evidence has routinely shown that the outer regions of the disc are highly structured (e.g. Long et al. 2018; Andrews 2020; Benisty et al. 2023); however, the challenging nature of studying the inner disc has still left the question of whether structures are present in the inner 0.5 au. These scales can only be marginally resolved with available observations (Andrews et al. 2016); however, what is known through photometric and spectroscopic surveys is that the region is highly complex and variable.

Time variability of the continuum emission, as well as of the accretion and wind-tracing lines, can offer a unique insight into the inner disc properties (Sicilia-Aguilar et al. 2020; Campbell-White et al. 2021; Armeni et al. 2023; Sicilia-Aguilar et al. 2023). For example, the variable light curve of AA Tau was interpreted to result from the occultation of the star by an inner disc warp induced by a misaligned magnetic field (Bouvier et al. 2003). Large photometric campaigns further revealed a class of young stars - ~30% of the young stellar object (YSO) population - the so-called dippers, with multiple brightness dips in their light curves (Alencar et al. 2010; Cody et al. 2014; Cody & Hillenbrand 2018; Roggero et al. 2021; Campbell-White et al. 2021). Some objects have very regular, small variations, likely due to spots on the stellar surface.

Other objects show very irregular and deep dips (accounting for up to 20-50% of the photospheric flux), likely due to variable obscuration by circumstellar dust (McGinnis et al. 2015). Such events, whose origin is a matter of debate, are more likely related to transient changes of the disc itself, or of the amount of dust lifted by winds (Bans & Königl 2012; Rodenkirch, P. J. & Dullemond, C. P. 2022), or entrained in the accretion columns. Since in several cases these irregular dips occur in discs that have low to moderate inclination at large scales (≤ 10-100 au) as probed by the Atacama Large millime-tre/submillimetre array (ALMA), they may trace the presence of highly misaligned inner discs (Ansdell et al. 2020; Scaringi et al. 2016). Further evidence can be found in the shadows left by the misaligned inner disc onto the outer disc as observed in scattered light (Benisty et al. 2017, 2018). In the case of the dipper RXJ1604-2130, the shadows vary in location and amplitude (Pinilla et al. 2018), with a short timescale that suggests that the obscuring dust is located very close to the star (Sicilia-Aguilar et al. 2020).

Considering the range of possible variability mechanisms, dip amplitudes, and frequency of stellar spot evidence, irregular dippers are therefore ideal targets to investigate the complex phenomena occurring in the inner disc. This can be done by comparing properties at different epochs, as the star suffers different degrees of obscuration. Time series of flux-calibrated spectra that simultaneously cover a broad range of wavelengths can provide unique information on accretion and wind properties, as well as on the location and properties of the obscuring grains (Schneider et al. 2018).

This paper presents the results of a survey of 16 irregular dippers in the star-forming regions Upper Scorpius, observed with X-Shooter mounted on the ESO’s Very Large Telescope (VLT), focussing on the analysis of the continuum variability. The paper is structured as follows. In Section 2, we present the VLT/X-Shooter observations and the sample properties. In Sect. 3 we describe the method of determining the stellar photosphere used to compare the spectra of the objects in a dip. Section 4 shows the first results of the survey along with the measured dust properties of the dips. Section 5 investigates the near-infrared (NIR) emission of the objects, while 6 looks at the line veiling of the dips to search for any connection to variable accretion. In Section 7, we discuss the findings and their significance for such objects, which is followed by our conclusions in Section 8.

2 Sample and observations

2.1 Sample properties

The sample comprises of 16 objects from the Upper Scorpius star forming region. They are classified as irregular dippers by Ansdell et al. (2016, 2018); Ansdell et al. (2020) and Hedges et al. (2018). We selected objects with frequent and deep dips (median dip depth of 20%), characteristics of irregular dips. The overall sample represents ~35% of the known dipper population in the region.

The selected objects are all low mass, late spectral type (SpT) stars, mostly K and M, with just one G-type star. Their stellar parameters are shown in Table 1. This table reports also the disc inclination as measured by ALMA (Carpenter et al. 2025; Ansdell et al. 2020; Barenfeld et al. 2017). In general, the selected objects have relatively high inclination angles (i) varying from 50° to 60°. Only two objects have measurements of i ≤ 20° (namely J155836 and J160421). These objects have been reported to have misaligned discs (Pérez et al. 2018; Benisty et al. 2018; Nealon et al. 2020; Mayama et al. 2018). See Appendix E for more details.

2.2 VLT/X-Shooter observations and data reduction

The data were taken with the X-Shooter spectrograph (Vernet et al. 2011). This work combines spectra from a collection of observing programmes, spanning from 2016 to 2021 (see Table F.1). The first programmes (Id. 097.C-0378, 0101.C-0866 and 0103.C-0887; PI Manara, observations in 2016, 2018 and 2019, respectively) account for several epochs of four objects (J155836, J160900, J160421, and J160207). The most recent programme (Id. 0105.C-0513(A); PI Manara) obtained in 2021 contains two to four spectra on additional 12 objects completing the dataset. The time interval between observations vary from more than a year to just one day. Therefore, for all stars, we obtained 2-4 spectra, with no specific strategy in terms of time separation. This method of stochastic observations maximises the probability of catching an object in a dipping state while also allowing the combination of archival and new data. The same observing strategy was used in all programmes. More details on the observations and data reduction for the data from 2016 and 2018 can be found in Manara et al. (2020), where a detailed study of the accretion properties was completed. Here the analysis was done using only the first epoch for each target.

X-Shooter can obtain simultaneous observations in its three arms, (UVB, VIS, and NIR) covering the wavelength range 3002500 nm, allowing for a broad wavelength coverage of the dipper variability. The spectra were obtained with narrow (1″.0, −0″.4, −0″.4) slit widths, leading to spectral resolutions of ≳5000 for the UVB, 18 000 for the VIS, and 10 000 for the NIR arms. Additional spectra with wide (5^,') slit widths were taken immediately after for flux calibration purposes. This was to ensure no flux losses due to varying seeing conditions. The slits were aligned at the parallactic angle to minimise atmospheric dispersion. The combination of the data leads to a database containing a total of 16 sources with two to four epochs per source (a total of 46 spectra).

Standard data reduction was completed using the X-Shooter pipeline v3.6.1 (Modigliani et al. 2010) with EsoReflex v2.11.5. The image reduction and analysis facility (IRAF) was used to extract the spectra in each arm (see Manara et al. 2020 for more details). After standard flux-calibration procedures, the spectra taken with the narrow slit widths were scaled to the flux calibrated spectra obtained with the wide slit set-up. The wavelength calibration from the pipeline was refined by using the Li I (6707.856 Å) absorption line characteristic of the photosphere of late type low mass stars (Campbell-White et al. 2023). Any remaining noise in the continuum was simply cut through high threshold (17σ) sigma clipping. No telluric correction was performed, regions effected by atmospheric interference were cut from the spectra. The final sample of reduced flux calibrated spectral energy distributions are shown in Appendix A.

3 The stellar photosphere

The reduced and flux-calibrated spectra (Appendix A) show the characteristic spectral energy distributions (SED) of Class II objects with little to moderate infrared excess in this wavelength coverage. The majority of objects have multi-epoch SEDs showing significant differences in their continuum levels. This is the dipping behaviour seen from a spectroscopic point of view. In some cases, the spectra also show different amounts of NIR excess between epochs.

In order to characterise the dipping events and probe their nature, the ‘intrinsic’ spectrum of the object (outside dimming events) needs to be identified and characterised. This is done by defining the inherent photospheric flux of the source, along with its spectral type. In the following we describe how the source spectral type and intrinsic continuum flux have been derived.

3.1 Spectral type

The spectral type of the objects in the sample was derived by comparing the observed spectra to a collection of Class III spectra following Manara et al. (2013b, 2017); Claes et al. (2024). As described in this work, Class III objects are preferred over synthetic spectra as the latter do not account for the high chromospheric activity characteristic of low gravity objects such as pre-main sequence stars. The Class III sample was selected from Claes et al. (2024)¹. As in our case, these spectra were observed with X-Shooter, and they were reduced and analysed in the same fashion. To find the best match between templates and observed spectra we considered the position and depth of the molecular bands in the visual at 705-706 nm, 708 nm-709 nm, and 712713 nm in the case of M-type stars. For earlier spectral types, the determination of the SpT is based on a Mg-line index (Herczeg & Hillenbrand 2014; Claes et al. 2024). The overall continuum shape and presence of a NIR excess was also taken into account when determining the spectral type. To determine the uncertainty in the SpT, the molecular bands and previously mentioned spectral index of the same object were derived at different epochs. On average, the depth of the molecular bands used in the spectral type determination of M-type stars vary by <5% in six out of the eight objects, and by <15% in the remaining two. For earlier spectral types, variations of <10% are found. In all cases, this translates into an uncertainty of the SpT of less than two subclasses.

The derived spectral types are reported in Table 2, together with the adopted Class III template. The differences with previous estimates are always within three subclasses. It should be noted that in contrast with our study, previous spectral type determinations of our objects may have been performed from single epoch observations, and therefore, local variations in the extinction - causing the dipping events, not related to the interstellar medium (ISM) - were not taken into account. This can account for the reported differences.

Using the method described above, suitable templates were found for all our sources with the exception of two objects: J162647 and J161544. Both objects show a very red SED, with J161544 being potentially a misclassified Class I object. Both objects are excluded from further analysis.

3.2 Normalisation of the class III template and interstellar extinction

Once the SpT is determined, the template spectra need to be calibrated to the intrinsic photospheric flux. Multi-epoch observations of dippers show that their fluxes decrease multiple times from a maximum flux level that is seldom overreached. Unfortunately, the sparse and sporadic time coverage of our data do not allow us to determine the maximum (or quiescent) flux level of our sources. However, for 13 out of 16 objects, Gaia offers time-sequenced photometric measurements (Gaia Collaboration 2016, 2023). We can make use of this multi-epoch Gaia photometry to determine the quiescent flux level of our sources. Furthermore, as Gaia optical observations are free from excess emission due to the presence of, for example, a disc, one can assume the quiescent flux level as the intrinsic photospheric flux of our sources. This assumes that the observed variation in the flux level of our spectra at optical wavelengths is due to a variable extinction due to the sum of line-of-sight extinction, that should be constant in time, and of a local, unknown contribution, variable in magnitude and wavelength dependence, which causes the dips.

To obtain the photospheric flux level of our sources and the ISM reddening, synthetic G, BP and RP magnitudes were estimated from our Class III templates. The templates were reddened by the ISM extinction using Cardelli et al. (1989) with R_V = 3.1. The ISM extinction was varied until there was a good match between the BP and RP magnitudes normalised to the Gaia G band magnitude at quiescence, and the reddened synthetic magnitudes of our Class III template. It should be noted that normalising to the broad band, Gaia G filter, also corrects for the distance to the objects. Figure B.2 plots the difference between the observed Gaia magnitudes and the synthetic magnitudes of the adopted photospheric flux. The derived ISM extinction values are included in Table 2. The derived values agree with those reported in the literature within typically ∆A_V = 0.5 mag. The only exception is J160421 where a difference of 0.8 mag is found.

For two of the three objects with no multi-epoch Gaia photometry (J162647 and J161544), the photospheric flux calibration was done in a similar fashion but using the available photometric measurements from the literature. These were compared with the available single epoch G-band photometry to determine the brightness state of the source. It should be noted that for these two sources, the resulting photosphere has a larger uncertainty than the rest of the sample, as the value of A_V cannot be accurately constrained. Finally, no representative Gaia photometry was reported for J161411, and therefore the photospheric level of the source could not be determined. This source was excluded from further analysis.

The reconstructed, reddened intrinsic continuum photospheric emission of each star is shown in Fig. A.1 alongside the observed spectra. Following the discussion above, the total sample included in this work contains 13 objects.

4 Dust properties as traced by the dips

In order to probe the continuum variations from epoch to epoch, the ratio of the ISM corrected observed spectrum to the photospheric one, R_λ, was estimated (see Fig. C.1). The ratios were computed from the observed spectra smoothed with equally spaced bins in wavelength (95 Å), and sigma-clipped to remove any effects of emission lines. To better visualise the results, Fig. 1 shows the value of this ratio at λ = 550 nm (R₅₅₀) computed for each multi-epoch observation. As shown in the figure, obscuring events are variable from source to source and epoch to epoch for the same source, with values generally ranging from R₅₅₀~1 (i.e. no dip) to R₅₅₀~0.5. Only 7 events (belonging to 4 different objects) show R₅₅₀ values ≤ 0.5. In some cases, dimming events are clustered around a specific value (e.g. J160418). Fig. 1 also highlights that some spectral types (e.g. K3, M1) have a wide spread of R550values indicating no relation between spectral type and the dip ratios. On a separate note, there is no trend seen when comparing the time separation between observations and values of R550 . This confirms that our observing strategy samples randomly the dipping history of the objects, and, also, that these are indeed irregular dippers.

The most likely scenario for the presence of the observed dimming events is the presence of local dust, increasing the extinction towards the source, and/or the presence of cold spots on the stellar surface. Both of these scenarios will produce a decrease in the observed flux. Cold spots are generally shortlived and appear regularly in time, at least for a period, in phase with the stellar rotation. In order to relate the observed dimming events to the presence of cold spots, we analysed K2 curves (Howell et al. 2014; Huber et al. 2016) of several of our objects in order to look for regular variability and assess the degree of flux variations. Based on our analysis, we find that regular dips, likely associated with cold spots, are unlikely to reduce the photospheric flux by more than ~15%. Therefore, dimming events above this limit are most likely due to variable local dust obscuration. In this section, we will focus on probing the properties of the dust causing the dimming events. Therefore, we will limit our study to events with R₅₅₀ < 0.85, avoiding significant contamination from cold spot contributions. Taking this limit into account, our analysis will focus on 20 dimming events from a total of 11 sources (marked with a * in Table 1).

Fig. 1 Ratios of the observed to photospheric flux at 550 nm (R550) for each star (see colour legend) are shown as a function of the spectral type. When not visible, the errors are smaller than the size of the dots. Dotted and dashed grey lines represent R550thresholds of 1.0, 0.85 respectively. No relation between spectral type and R550is observed.

4.1 Wavelength dependence of the optical depth

The wavelength dependence of R_λ provides very important clues on the properties of the local layer of dust causing the dimming events. From R_λ, it is possible to estimate the optical depth of the obscuring layer of dust as $${\tau }_{\lambda }=-\ln R_{\lambda }.$$(1)

A rough estimate of the optical depth dependence on the wavelength can be estimated by computing β, which represents the ratio of the optical depths at two reference wavelengths λ₁ and λ₂: $$\beta =\frac{\log ({\tau }_2/{\tau }_1)}{\log ({\lambda }_2/{\lambda }_1)}.$$(2)

We choose to focus our measurement in the visual range of the spectrum (λ₁ ≈ 400 nm and λ₂ ≈ 900 nm). In doing so we avoid contamination from the lower signal-to-noise ratio (S/N) in the UVB, and the hot dust emission that dominates in the NIR.

The values of τ₁, τ₂ were computed from the weighted mean of the three nearest values at λ₁, λ₂ from the binned τ_λ values, as can be seen in Fig. D.1. The uncertainties on the values of τ₁,τ₂ are derived from the noise of the spectra in each bin of R_λ. These are carried through, using the standard method of error propagation, in computing τ_λ, the weighted means, and ultimately the value of β. Note that τ_λ and β refer to the local dust only as ISM extinction correction has already been applied as explained in Section 3.2 and reported in Table 2. Figure D.1 shows τ_λ for the 20 dimming events associated with the presence of dust. Note that τ_λ and β refer to the local dust only as ISM extinction correction has already been applied as explained in Section 3.2 and reported in Table 2. Figure 2 shows the values of β versus the optical depth at 550 nm, τ₅₅₀. As a reference, the β value corresponding to the ISM extinction taken from Cardelli et al. (1989) assuming R_V = 3.1 is also shown in the plot.

In practice, Fig. 2 shows a crude measure of the dependence of τ on λ (τ ∝ λ^−β), and thus provides information about how the dust extinction (i.e. effective scattering and absorption opacities; k_λ) varies between λ₁, λ₂ for each dip. As shown in the figure, β values are always lower, and in many cases significantly lower, than that expected from ISM reddening, except for two dimming events (for objects J160207 and J160205). It should be noted that five dips show negative β values implying that the dust extinction increases with increasing wavelength. This is the opposite behaviour expected from ISM extinction, and it shows a more extreme behaviour than that expected for large grain size distributions (i.e. grey extinction, β ~ 0; Draine 2003).

Fig. 2 β (Eq. (2)) values as a function of τ550 for objects with R550<0.85. Multiple dips of the same object are colour coded according their shortened 2MASS ID (see Table 1). Dashed blue line represents the value expected for ISM reddening under typical conditions, assuming the Cardelli extinction law (Cardelli et al. 1989) with RV=3.1.

4.2 Dust grain models

To further interpret the dependence of τ on wavelength and gain some insight into the nature of the dust causing the dimming events, dust opacity models were computed using the Optool code (Dominik et al. 2021). The code performs opacity calculations for various dust grain properties and conditions. It should be noted that the analysis presented here is not intended as an exhaustive study of the dust properties of the disc, but as a tool to give rough estimates of the dust properties producing the dimming events and the relative dust property variations between epochs.

We consider grains composed of composite silicate materials (astrosilicates; Draine (2003) with a size distribution n(a) ∝ a^−q between a minimum and a maximum size a_min and a_max. The grains are assumed to be compact and spherical. We compute absorption (k_abs) and scattering (k_sca) opacities, and the phase function (< g >) as a function of wavelength. Note that we primarily consider single scattering when computing k_sca and < g >. However, considering pure absorption, k_abs only, can be interpreted as the extreme case of multiple scattering. For a first comparison to the observations, we characterise each curve by the corresponding value of the quantity β, defined as in Sect. 4.1. We vary a_min and a_max between 0.03 μm and 100 μm and fix q = 3. Figure 3 (first and second panel from top) shows the values of β for absorption and scattering opacities separately. In reality, the actual attenuation of the stellar radiation due to the intervening dust is a complex function not only of the dust properties but also of the geometry of the system (Lakeland & Naylor 2022). Two additional examples are provided in the lower panels of Fig. 3, showing the values of β for extinction (k_ext = k_abs + k_sca), and for the case where the contribution of scattering is weighted by the phase function (κ_eff = κ_abs + κ_sca(1 - < g >); Birnstiel et al. 2018).

As shown in Fig. 3, distributions with both a_min and a_max < 0.3 μm result in high values of β (≳ 1.4). If absorption dominates, a wavelength-independent attenuation (β ~ 0), as detected in several spectra, can only be reproduced if the smallest grains (i.e. a_min ≳ 1-2 μm) are eliminated. The need to significantly increase a_min remains if (single) scattering is included, even if less extreme. In the case of extinction, it is possible, for example, to have β ~ 0 with a_min ≳ 0.2, as long as a_max ≳ 1-2 μm. The very interesting cases of negative β, as measured in some spectra, can only be explained if scattering plays a significant role in the attenuation, and only if the grain size distribution is very narrow, limited to the range from few tenths to ~1 μm (dark blue areas in Fig. 3). Note that, because of the range of wavelengths used in the definition of β (i.e. 400-900 μm), the results are in all cases not sensitive to values of a_min ≪ 0.3 μm and a_max ≫ 10 μm, respectively.

Although the trends discussed above are likely to be qualitatively correct, similar values of β can be obtained by different sets of parameters even for a specific grain model. In addition, it is difficult to account for the effect that the uncertainties on β have in the interpretation of Fig. 3. Moreover, it should be noted that β provides a very crude description of the actual wavelength dependence of the observed optical depth, as it is based on two reference wavelengths only (see Equation (2)). With these considerations in mind, we examined each of the 20 τ(λ) curves seen in Fig. D.1. For each τ(λ), we derived the best fitting values of the grain size distribution parameters over the whole range between 400 nm and 900 nm, assuming, for simplicity, that the observed optical depth is due to extinction. As we are interested in the possibility that small grains (e.g., <1 μm) do exist in the obscuring dust, we explore the full range of a_min ≥ 0.03 μm, decreasing q when required. This has the effect of making the grain size distribution more top heavy. The results are shown in Fig. D.1 and summarised in Table 3. Note that we do a visual fit only as we consider that the unknown geometry and the uncertainties in the actual grain composition will make any statistical estimate of the parameters’ errors meaningless.

Inspection of Figs. 3, D.1, and Table 3 suggest that the observed τ(λ) curves can be classified into three different groups, based on their wavelength dependence and matching grain size distributions. Figure 4 shows the three characteristic examples of this behaviour: decreasing τ(λ) with increasing wavelength (top panel); very flat τ(λ) (middle panel); curves not well described by a power-law in the Logτ-Logλ plane (bottom panel). The curves shown in the top panel (i.e. decreasing τ(λ) with increasing wavelength) can be reproduced by extinction curves with q = 3, a_min ~ 0.03 μm, and a_max ranging from tenths to micron size particles. In some cases (e.g. J160205), a_max is also small, and the grains have a size distribution similar to the ISM one. The curves with very flat τ_λ (β ~ 0.0) (Fig. 4; middle panel) are consistent with very small values of a_min ~ 0.03 μm, provided that a_max is very large, of the order of a few tenths of micron or more. Finally, the example of curves not well described by a power-law in the Logτ-Logλ plane (Fig. 4; bottom panel) often results in β < 0, and they require a narrow size distribution of grain sizes, in the tenth of micron range. Note that these kinds of curves require significant contributions of scattering to the total opacity. In this case, the albedo (ω = k_sca/k_ext) at λ = 550 nm increases from 0.60 to 0.86 from the example shown in the middle to the bottom panel.

In general, very few dips show evidence of ISM-type grains. Typically, the dips sample in this study require flat τ(λ) curves, typical of grain size distributions with a_max much larger than in the ISM, reaching in several cases, sizes of >10 μm. The minimum size is often consistent with ISM values, but in all cases the grain size distribution is shifted to much larger sizes, providing strong evidence that the grains in the occulting region come from a population that has undertaken significant growth.

Fig. 3 β colour maps from dust opacity models. Each figure represents a grid of dust size distributions varying minimum and maximum grain sizes (n(a) ∝ a−q with q=3.0). Coloured by the values of β derived from the opacities for (a) pure absorption (βabs), (b)scattering (βsca), (c) extinction (βext), and (d) effective extinction (βeff; see text). Contours are marked by black lines at β = 0.0 and βism = 1.38). Pattern effects are artefacts of the interpolation used.

Fig. 4 Optical depth as function of wavelength for three representative dips. The observed values are shown by grey dots with error bars. The red dashed line shows the extinction opacity (absorption + scattering) for the grain size distribution that best fit the data in the wavelength interval 400900 nm (limits shown by the vertical grey dotted lines), normalised to the observed value at 700 nm. The values of amin, amax, q are given on top of each panel; the corresponding grain size distributions (n(a)) are displayed in the inset plots.

5 Near-infrared emission

So far, our analysis has been focussed on the visual range of wavelengths, where the observed flux variations are not affected by the innermost disc hot dust emission. However, an analysis of the NIR excess can potentially shed further light on the dimming events and in particular on the location of the obscuring dust.

In order to quantify the NIR emission and how it varies among dips, the excess emission was calculated for all objects for which a photospheric level was determined. We define eK to represent the amount of excess flux (F_obs, F_phot) at a noisefree region in the K-band of their SEDs (2120-2150 nm). This is calculated as $${ϵ}_K=\frac{F_{obs}}{F_{phot}}-1.$$(3)

Discs containing little to no NIR emission in the K-band will therefore have ε_K ~ 0, while strong discs will have values ≥0.2.

Figure 5 (top panel) shows the computed ε_K values for all the dimming events as a function of R₅₅₀. At ε_K = 0, a conservative uncertainty of ±20% was taken into account (grey area in Fig. 5). All other errors were derived from the RMS of the continuum spectra. As shown in the top panel, most of our measurements are located at 0.0 ≤ ε_K ≤ 0.2, and R₅₅₀ ≥ 0.75 (i.e. τ₅₅₀ ≥ 0.28), indicating very little excess emission at NIR wavelengths, and none to little difference between the photospheric flux and the single epoch measurement (i.e. none or small dimming events). This is in agreement with the SEDs of these objects (see Fig. A.1). Only seven epochs show ε_K >  0.2, suggesting stronger disc emission. This is again in agreement with the SEDs of these objects (namely J155836, J162218, and J160421) as shown in Fig. A.1. Interestingly, eight epochs from four different objects (J160647, J160418, J160727, and J160207) show ε_K <  0. This would suggest the unlikely scenario in which the photospheric flux contribution is higher than the observed total flux.

Based on these results, and in order to get hints on the location of the obscuring dust, the excess emission was recalculated by correcting the observed flux by the obscuring dust (local extinction). This was only done for those objects showing dimming events caused by obscuring dust, i.e. with R₅₅₀ < 0.85 as explained in Sect. 4. Correcting ε_K for the local extinction is equivalent to assuming that the obscuring dust is occulting both the photosphere as well as the hot dust disc emission (assumed to be constant for simplicity). The local extinction was estimated by extrapolating the optical depths of the dust grain distribution (see Sect. 4.2) to the NIR region. The observed flux was then multiplied by e^τ_λ. The results are shown in the bottom panel of Fig. 5.

As expected, most of the dimming events have now higher ε_K values after correcting for the local NIR local extinction with only one epoch remaining with ε_K < 0. It should be noted, however, that for this object (J160727) no multi-epoch Gaia data were available and thus the estimated photospheric level is highly uncertain (see Sect. 3). After correction, five epochs show ε_K ~ 0, while seven dips now show ε_K > 0.3 indicating a stronger infrared excess emission. Fig. 5, the dimming events with R₅₅₀ <0.6 (i.e. higher optical depths) are those showing the greatest changes in ε_K, moving from ε_K < 0 to ε_K ~ 0 (with the exception of one epoch as discussed above). In these cases, it is possible that both the photosphere and the NIR disc emission are obscured by the same dust causing the dimming event. Although any conclusions on the location of the obscuring dust are difficult to make given the generally low ε_K values. Interestingly, two objects (namely J155836 and J160207) show epochs with large ε_K changes of 0.57 and 0.30. These variations imply a degree of intrinsic K-band variability. Although due to the small statistics it is difficult to draw further conclusions and ascertain whether these changes are related to those seen in the optical. Similarly, it is also difficult to shed further insights onto the location of the obscuring dust for the remaining dimming events as ε_K variations are likely due to a combination of several factors such as the wavelength dependence of the optical depth, and the strength of the NIR excess.

Fig. 5 Excess emission at 2.1 μm, εK shown as a function of R550. The top panel plots the ratio between the observed flux and the intrinsic photospheric one (${ϵ}_K=\frac{F_{obs,K}}{F_{ph,K}}-1$). The bottom panel shows the ratio between the observed flux corrected for the extinction derived at optical wavelengths and the photospheric flux (${ϵ}_K=\frac{F_{obs,K}{e}^{\tau }}{F_{ph,K}}-1$). The grey area covers a ±20% region around εK = 0. Vertical error bars on εK are smaller than the size of the points.

6 Accretion variability

YSOs are also known to undergo short term changes in their levels of accretion, adding a continuum emission to the photospheric one at all wavelengths. The result is a general increase in the observed continuum emission, stronger in the UV, but also clearly detected in the visual, where it causes the photospheric absorption lines to look shallower (veiling). The variable level of local extinction in dippers complicates the typical approach of measuring accretion properties in T Tauri stars from the UV continuum excess (Manara et al. 2013a, 2020; Claes et al. 2024). Measurement of the properties such as the accretion luminosity are therefore beyond the scope of this work. However, to investigate if the observed variability can be explained by changing emission levels, the ratio of the observed Equivalent Widths (EWs) of the LiI(6707.8Å) and Ca I(4226.7Å) lines with respect to each object’s brightest epoch was calculated (see Fig. 6). These lines are well-known proxies of veiling variations. In particular, the Li I 6707.8 Å line is minimally affected by accretion line emission (Campbell-White et al. 2023), while the Ca I is a diagnostic of veiling in the blue part of the spectrum (Herczeg & Hillenbrand 2008). If the photospheric line does not vary, the ratio of the EW between the two spectra should be proportional to the change of the continuum flux at the wavelength of the line. We find in all cases that in general the variations of the EW of both the LiI and CaI lines are <10%. Therefore, our sample show little evidence to suggest any of the observed variability can be explained by sudden accretion-linked emission changes leading to line veiling.

7 Discussion

In this work, we presented the most complete sample of irregular dippers in the Upper Scorpius star-forming region to date. In our analysis, we employed a novel approach which consisted of deriving the spectral type of our sources using a sample of Class III templates, and retrieving their quiescent photospheric level from multi-band Gaia photometric observations. This method allowed us to derive the intrinsic photospheric spectrum of 13 out of 16 objects. Unlike previous studies where the study of the dimming events is based on continuous high-cadence photometric light curves (Cody et al. 2014; Hedges et al. 2018; Roggero et al. 2021; Capistrant et al. 2022), here the relative variability is measured with respect to the retrieved photospheric level. This approach therefore relies upon the accuracy in reproducing the true stellar photosphere. We have shown that we are able to reproduce the Gaia observations with a high degree of precision using class III templates (see Fig. B.2). As already mentioned (Sect. 3) there is an uncertainty of less than two subclasses on the spectral type determination. However, we have found that this difference does not significantly change the derived values of R_λ,τ_λ,β whose uncertainties are dominated by the noise of the spectra. A side product of this procedure is that the derivation of the interstellar extinction was possible. Accurately identifying the true photosphere, spectral type, and A_V is critical for the studies of such variable objects, as the dimming events can lead to inaccurate measurements of key physical properties such as stellar and accretion luminosities, dust masses, etc. We have shown, for example, that there can be discrepancies between the spectral type (up to three subclasses) and A_V (up to 0.8 mag) measured in this study and those found in literature (Table 2).

X-Shooter observations allowed us to probe variations in the continuum spectrum from ~400 nm up to ~2.4 μm and simultaneously look for accretion variability. Our study showed that our objects show very little variation in their line veiling and therefore the observed photometric changes are not due to sudden accretion variability shocking the stellar surface. Furthermore, our analysis concluded that cold spots are also unlikely to reduce the photospheric flux by more that ~15%. Therefore, dimming events above this limit are believed be due to local dust obscuration, consistent with a number previous studies of irregular dippers (Bouvier et al. 1999; Alencar et al. 2010; McGinnis et al. 2015; Ansdell et al. 2016; Bodman et al. 2017; Schneider et al. 2018; Sitko et al. 2023).

In this regard, 11 out of 13 objects show an observed flux lower than that of the photosphere at optical wavelengths, indicating the sources were undergoing a dip. The magnitude of the dips varied between sources and observed epochs, ranging from a few percent to ~75% of the photospheric flux at 550 nm. While the spread in the degree of variability is common to other studies (e.g. Morales-Calderón et al. 2011; Ansdell et al. 2016; Hedges et al. 2018; Alencar et al. 2010; Roggero et al. 2021; Venuti et al. 2021), the deepest dips seen in our sample exceed previously recorded ones, which typically showed flux variations of up to 50%. In particular, our sample show six dips across three targets (namely J160727, J160418 and J160647) that show dips in excess of 50% of the quiescent photospheric flux level at 550 nm. Interestingly, our sample does not show any correlation between the disc inclination and spectral type with the strength of the dimming events. One might expect, for example, that deeper dips occur in systems that are more highly inclined as the line of sight intersects more closely the surface of the disc (Bouvier et al. 1999; Lakeland & Naylor 2022), or that they correlate with higher optical depths. The lack of evidence for such correlations highlights the unpredictability of these objects.

The analysis of the dependence of the dip optical depths on wavelength shows a variety of behaviours, ranging from very flat (low β values) to rather steep (high β) in the optical wavelength range (400-900 nm). Under the interpretation that the dimming events are due to obscuring dust, our results show a very dynamical inner disc region. They can have changing dust properties for different dimming events belonging to the same source. Under the assumption of a typical MRN power-law size distribution (n(a) ∝ a^−q) composed of solid astrosilicate grains, our results show that in most cases a top heavy dust grain size distribution is required. In all cases grain sizes larger than those expected in the ISM, are required in order to reproduce the corresponding optical depth wavelength dependence in the 400-900 nm range. There is a large variation in a_max, with no correlation with the depth of the dip. In some cases, our models can only provide lower limits for the maximum grain size with a_max ≳ 10-20 μm, and in others we require grain growth for a_min (up to tenths μm). Our analysis also shows that in some dips (e.g. J155836) scattering likely plays a significant role, as in several cases the optical depth dependence on wavelength cannot be reproduced by absorption only (see Fig. 3). While our models do test the extreme case of multiple scatterings (pure absorption), we focus on understanding the contribution to the extinction from single scattering events only. We find that in many cases the albedo of the local extinction is different from that of typical ISM conditions, with an increased contribution of scattering. Under the assumption of classical extinction (k_abs + k_sca), it is the scattering that dominates above the absorption allowing for the unusual τ_λ curves to be possible (see Fig. 3, τ_λ plot). This arises from the transition of Rayleigh to Mie scattering regimes. This increased scattering effect is a result of the changing grain sizes and the fact that the dips are typically optically thin, instead of any geometrical effects. The scattering contribution is important to consider when observing such targets, and has been previously linked with the dipping mechanisms (Bouvier et al. 1999, 2003; Dodin et al. 2019). Interestingly, in these cases the dips are linked to very narrow grain size distributions with both a_min and a_max of the order of tenths of microns. Given the typical timescale of the dips, the obscuring dust is understood to be coming from the inner disc regions (<<1 au; Morales-Calderón et al. 2011). The fact that dust grains of sizes up to tens of μm exist in these regions suggests some level of grain growth has occurred (Nagel & Bouvier 2020; Sitko et al. 2023).

Regarding the NIR emission at 2 μm, by assuming the same dust model as those derived at optical wavelengths and propagating into the NIR we were able to give some constraints on the location of the obscuring dust. Although the short time scale variability of dippers is generally linked to the closest-in location of the obscuring dust (e.g. Morales-Calderón et al. 2011), here we were able show that in some cases it may be possible that both the photosphere and region emitting the 2 μm continuum emission are obscured by the same dust. In most objects, the extinction-corrected 2 μm continuum emission does not vary between epochs. However, in two cases (J155836, J160207), there remains a difference in excess emission between two dips. This suggests an intrinsically varying NIR emission for these objects. However, in general low NIR excess emissions combined with the range of possible uncertainties make it difficult to draw firm conclusions.

Although similar results regarding the location and spread in dust grain size distributions, including for some objects in our sample, have also been reported in literature (see e.g. Bouvier et al. 1999, 2003; Schneider et al. 2018; Sitko et al. 2023), in this study, we show the presence of relatively quick changes in the grain size distribution in the inner disc. With time-scales ranging from a few weeks to several months, the dust grain size distributions can change from epoch to epoch. The distributions in some cases are top-heavy and others they have narrow sub-micron dust grain distributions dominated by scattering. The changes are in general too quick to be explained by simple radial drag of different dust grain populations migrating from outer to inner regions of the disc. This points towards some mechanism capable of accessing various grain sizes from differing depths in a vertically stratified disc structure. Recent ALMA and SPHERE observations of protoplanetary discs show indeed the presence of different grain size populations at similar spatial scales (e.g. Andrews et al. 2018; Garufi et al. 2024). For instance, SPHERE scattered light images detect stellar radiation that is scattered off small dust grains at the disc surface, whereas at the same radii, ALMA dust continuum emission reveals the presence of much larger dust grains at deeper vertical depths (e.g. IM Lup Andrews et al. 2018; Avenhaus et al. 2018).

One mechanism capable of lifting such dust grains to significant heights above the disc plane to obscure the star, and sometimes also the region emitting at 2 μm, independent of the disc inclination, is the so-called dusty disc wind. In recent years, several studies have investigated the entrainment of dust grains in MHD and photoevaporative winds (Owen et al. 2011; Miyake et al. 2016; Giacalone et al. 2019; Booth & Clarke 2021). These models predict different regions from which these winds will be launched as well as different sizes of dust grains that can be lifted at those radii. In general, models predict a maximum grain size that can be lifted from a certain distance from the central source. Some of these models also include a vertical dust grain size stratification with grain sizes at 1 au ranging from submicron to 100 μm (see e.g. Miyake et al. 2016). For a vertical cut in the disc, at a certain distance from the source, the size of the grains that can be lifted depends on the mass flux of the disc wind (Miyake et al. 2016; Giacalone et al. 2019). Therefore, variations in the mass flux of the disc wind might explain the changes observed in the grain size distributions needed to reproduce the observed wavelength dependence of the optical depth in the optical. Dusty disc winds have been previously used to explain the dimming events in RW Aur (Dodin et al. 2019) as well as scattered light observations of the disc surface of IM Lup (Avenhaus et al. 2018), and in general linked with YSO variability (Ellerbroek et al. 2014; Lakeland & Naylor 2022; Grinin et al. 2023; Sitko et al. 2023).

Alternatively, dust entrainment along magnetospheric accretion columns could also explain the presence of the obscuring dust at high scale heights from the disc surface. Stable and unstable accretion columns have been shown to be able to lift dust grains above the disc surface following the magnetic field lines (Romanova et al. 2003; Nagel et al. 2024). The evaporation of the dust grains approaching the central star may explain the different levels of optical depths encountered. While this configuration tends to favour occultation of systems at high inclinations, varying dust survival along the accretion columns might be able to explain lower inclinations down to 45° (Nagel et al. 2024). In this context, unstable accretion columns would be particularity favourable for irregular dippers given their short timescales and lack of periodicity, and was believed to be the case for four systems (J160418, J160205, J160727 and J160900) in this sample (Bodman et al. 2017). Stellar occultation by disc warps in the inner regions, arising from such accretion events, is also a possibility (Romanova et al. 2013; Sicilia-Aguilar et al. 2020; McGinnis et al. 2015). Instabilities arising from such mechanisms have also been shown to affect the dust grain growth (Aly et al. 2024). Given that a number of the objects in this sample have been linked to this accretion driven variability, we cannot rule it out as a possible dip driving mechanism.

Fig. 6 Equivalent width (EW) ratios between epochs of the Li I 670.78 nm line (top panel) and Ca I 422.67 nm (bottom panel) photospheric absorption lines. Ratios are calculated with respect to the brightest epoch. Results show that in all cases EW changes are within 10%. Dotted black lines are at ratios equal to 1, signifying no change in the EW.

8 Summary and conclusions

In this paper, we presented the results of a programme aimed at studying the properties of a number of irregular dippers in the Upper Scorpius star-forming region. The sample includes 16 objects, each of them was observed with X-Shooter two to four times, with time intervals ranging from days to years, for a total of 46 spectra. Our main results are summarised below:

• By combining multi-epoch Gaia photometry with the spectral type derived from Class III templates, we find evidence that for 11 out of 13 objects the observed emission was lower than that of the photosphere, that is the star was undergoing a dip;

• The magnitude of the dips varies between epochs even for the same star ranging from few percent to up to ~75% of the photospheric flux at 550 nm;

• Assuming that the dips are caused by dust obscuration, we find that in most cases the grains causing the dips favour a top-heavy size distribution (i.e. most of the dust mass is in grains significantly bigger than those in the ISM). We find in some cases lower limits on the maximum grain size to be tens of microns. However, for a few cases, scattering also plays a significant role. In this case, a narrow grain distribution of sub-micron dust grains is needed;

• The properties of the obscuring dust changes from dip to dip and from source to source, with no evident relationship between the properties of the dip and stellar properties such as spectral type or disc inclination;

• By analysing the continuum emission at 2 μm it was shown that most of our objects do not show significant NIR excess;

• By assuming the same dust model as at optical wavelengths, the extinction at 2 μm was estimated to constrain the location of the obscuring dust. Although a detailed modelling of the dust properties is beyond the scope of this paper, under this simple assumption it was possible to infer that in a few sources the same dust may be obscuring both the source and the 2 μm emitting region of the disc. A few sources also show hints of 2 μm continuum variability.

In general, the picture that comes out of our results is that of very variable phenomena, not just in the depth of the dips but also in terms of the properties of the dust that cause them. This is true not only between different stars, but also for the same star, that can have dips due to grains of different properties, with no obvious trend that links that to the magnitude of the dip. These changes in the dust properties of different dips also point towards a highly dynamic dust growth and/or transport process in the inner disc. To further answer the questions on dipper variability, simultaneous spectroscopic and photometric observations covering full dimming events would be highly beneficial. Additionally, scattered light observations could inform of the presence of the dust substructures we believe to causing the dips.

Acknowledgements

We thank the anonymous referee for their careful consideration their comments that helped improve the clarity of our work. A.E. acknowledges funding from Taighde Éireann - Research Ireland under Grant number GOIPG/2023/4396 and the UCD Physics Scholarship in Research and Teaching (SIRAT). CFM is funded by the European Union (ERC, WANDA, 101039452). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. A.E. acknowledges support from the ESO Early-Career Scientific Visitor Programme when part of this work was carried out. MB has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration programme.

Data availability

This project was based on observations collected at the European Southern Observatory under ESO programmes 0105.C-0513(A), 097.C-0378(A), 0101.C- 0866(A), and 0103.C-0887(B). Reduced spectra can be found on the ESO Science Archive Facility with DOI: https://doi.eso.org/10.18727/archive/71.

References

Appendix A Spectral energy distributions

Figure A.1 shows the observed spectra for all the targets where the photospheric spectrum could be derived (see sect. 3).

Fig. A.1 Observed spectra for each epoch are plotted alongside the ISM-reddened photosphere where applicable.

Fig. A.2 Observed spectra for objects for which a photosphere could not be derived.

Appendix B Gaia photometry

We show in this section (fig. B.1 the observed photometric broad-band (G) magnitudes available from the Gaia Variability programme (Gaia Collaboration 2016, 2023). In each panel, the dots show the observed value. The grey area highlight the range of magnitudes between the maximum and 85% of it. The dashed line plots the median of the G values included in the region (in magenta), that we assume to be representative of the photosphere of the object.

Fig. B.1 G-band magnitudes observed by the Gaia Variability programme (Gaia Collaboration 2016, 2023) for all the objects included in the Gaia dataset. In each panel, the observed magnitudes are shown by dots. The shaded region defines the range between the minimum and 85% of it (i.e. the range of the brightest values). The dotted line shows the median of the points included in this range (in magenta), and is assumed to be the G-band magnitude of the stellar photosphere.

The procedure followed to derive the ISM extinction and the photospheric spectrum for each object is described in sect. 3. Fig. B.2 shows for each object the difference between the reddened synthetic magnitude of the adopted photosphere and the observed values. The difference is always smaller than 5%.

Fig. B.2 Difference in magnitude between the synthetic magnitudes of the reddened adopted photosphere and the Gaia observed values for the BP and RP filters. By definition, G-band difference between the reddened adopted photosphere and the observed magnitude is 0. Dashed lines mark the ±5% interval around ∆m = 0.

Appendix C Ratios

Fig. C.1 shows the values of the ratio of the observed to the ISM-reddened photospheric fluxes as function of the wavelength.

Fig. C.1 Ratio Rλ of the observed over the ISM-reddened photospheric flux versus λ. The ratios are computed from the spectra binned in equal intervals of 95Å and σ-clipped.

Appendix D Grain properties

Fig. D.1 shows the wavelength-dependence of the optical depth for the 20 dips studied in Sect. 4.2.

Fig. D.1 τλ as a function of λ for all the 20 dips analysed in Sect.4. The name of the star and the epoch number are given on the secondary axis. The observed values are shown by black dots with error bars. The red coloured lines show the results for the grain size distribution that best fit the data in the wavelength interval 400-900 nm for the extinction opacities, each normalised to the observed value at 550 nm. The values of amin, amax, q are given on top of each panel; the corresponding grain size distribution is displayed in the insert. The vertical grey dotted lines represent the boundaries of λ = 400, 900 nm. Three representative cases are also shown in the main paper.

Appendix E Individual dippers

Of the objects in this study, two are particularly noteworthy for being well studied in the context of YSO variability, J155836 and J160421.

Appendix E.1 J155836

J155836, often referred to as HD 143006 in the literature, is a K0 T Tauri star (Manara et al. 2020) (M_* ~ 1.5 M_⊙) with shallow, irregular photometric dips (Ansdell et al. 2018). Its outer disc is highly structured, with a series of rings and gaps, and a bright arc exterior to the rings apparent at high angular resolution in sub-millimetre images (Andrews et al. 2018; Pérez et al. 2018). In the infrared scattered light, the disc appears to show a large-scale asymmetry, with the eastern side of the disc brighter than the western side (Benisty et al. 2018), indicating a shadow cast by a potentially misaligned inner disc region, inaccessible to direct imaging observations. Interestingly, the innermost ~7 au of the disc appear depleted of large dust grains (emitting in the sub-mm) (Jennings et al. 2022), and the ring beyond the cavity could be misaligned by ~40° with respect to the outer disc (Pérez et al. 2018) possibly due to a strongly inclined binary companion within the cavity (Ballabio et al. 2021). Nevertheless, the cavity is not fully depleted of material as the star is accreting (Manara et al. 2020) and the emission of small micron-sized dust at the dust sublimation edge is spatially resolved with the Very Large Telescope Interferometer (Lazareff et al. 2017).

Appendix E.2 J160421

In contrast, J160421, a K3 star, shows negligible accretion (Manara et al. 2020). It hosts a nearly face-on transition disc with a large depleted cavity while being one of the most massive discs in Upper Sco (Barenfeld et al. 2017). Its outer disc resolves into a single ring (Mayama et al. 2012, 2018) that has clear indication of being shadowed by a misaligned inner disc. Unlike the case of J155836, the shadows appear azimuthally narrow (Mayama et al. 2012; Pinilla et al. 2018) and highly variable in both their location and morphology within very short timescales of a day. Some of the properties of the shadows can be explained with a broken circumbinary disc, and a precessing, strongly misaligned, inner disc (Nealon et al. 2020). Additional support for the presence of a misaligned inner disc responsible for the shadows are found in the warped kinematics of the gas within the cavity and in the shadow counterparts in the sub-mm (Mayama et al. 2018). Davies (2019) analysed projected rotational velocities and determined that the stellar rotation axis inclination angle strongly differs from the one of the outer disc, indicating that the inner disc and the star are well aligned. This system was qualified as a dipper from its K2 light curve, and shows aperiodic dimming events of ~1 mag in depth, superposed onto a ~5.1 days periodic signal (Ansdell et al. 2016, 2018). Time-resolved multi-wavelength photometric studies recover the ~5 days quasi-periodicity but show variability in the eclipse depths and frequency, indicating an irregularly shaped inner disc, whose mass is completely accreted (and then replenished) in days to weeks timescales (Sicilia-Aguilar et al. 2020).

Appendix F Log of the observations

All Tables

All Figures

	Fig. 1 Ratios of the observed to photospheric flux at 550 nm (R550) for each star (see colour legend) are shown as a function of the spectral type. When not visible, the errors are smaller than the size of the dots. Dotted and dashed grey lines represent R550thresholds of 1.0, 0.85 respectively. No relation between spectral type and R550is observed.
In the text

	Fig. 2 β (Eq. (2)) values as a function of τ550 for objects with R550<0.85. Multiple dips of the same object are colour coded according their shortened 2MASS ID (see Table 1). Dashed blue line represents the value expected for ISM reddening under typical conditions, assuming the Cardelli extinction law (Cardelli et al. 1989) with RV=3.1.
In the text

Fig. 3 β colour maps from dust opacity models. Each figure represents a grid of dust size distributions varying minimum and maximum grain sizes (n(a) ∝ a−q with q=3.0). Coloured by the values of β derived from the opacities for (a) pure absorption (βabs), (b)scattering (βsca), (c) extinction (βext), and (d) effective extinction (βeff; see text). Contours are marked by black lines at β = 0.0 and βism = 1.38). Pattern effects are artefacts of the interpolation used.

In the text

Fig. 4 Optical depth as function of wavelength for three representative dips. The observed values are shown by grey dots with error bars. The red dashed line shows the extinction opacity (absorption + scattering) for the grain size distribution that best fit the data in the wavelength interval 400900 nm (limits shown by the vertical grey dotted lines), normalised to the observed value at 700 nm. The values of amin, amax, q are given on top of each panel; the corresponding grain size distributions (n(a)) are displayed in the inset plots.

In the text

Fig. 5 Excess emission at 2.1 μm, εK shown as a function of R550. The top panel plots the ratio between the observed flux and the intrinsic photospheric one (${ϵ}_K=\frac{F_{obs,K}}{F_{ph,K}}-1$). The bottom panel shows the ratio between the observed flux corrected for the extinction derived at optical wavelengths and the photospheric flux (${ϵ}_K=\frac{F_{obs,K}{e}^{\tau }}{F_{ph,K}}-1$). The grey area covers a ±20% region around εK = 0. Vertical error bars on εK are smaller than the size of the points.

In the text

	Fig. 6 Equivalent width (EW) ratios between epochs of the Li I 670.78 nm line (top panel) and Ca I 422.67 nm (bottom panel) photospheric absorption lines. Ratios are calculated with respect to the brightest epoch. Results show that in all cases EW changes are within 10%. Dotted black lines are at ratios equal to 1, signifying no change in the EW.
In the text

	Fig. A.1 Observed spectra for each epoch are plotted alongside the ISM-reddened photosphere where applicable.
In the text

	Fig. A.2 Observed spectra for objects for which a photosphere could not be derived.
In the text

Fig. B.1 G-band magnitudes observed by the Gaia Variability programme (Gaia Collaboration 2016, 2023) for all the objects included in the Gaia dataset. In each panel, the observed magnitudes are shown by dots. The shaded region defines the range between the minimum and 85% of it (i.e. the range of the brightest values). The dotted line shows the median of the points included in this range (in magenta), and is assumed to be the G-band magnitude of the stellar photosphere.

In the text

	Fig. B.2 Difference in magnitude between the synthetic magnitudes of the reddened adopted photosphere and the Gaia observed values for the BP and RP filters. By definition, G-band difference between the reddened adopted photosphere and the observed magnitude is 0. Dashed lines mark the ±5% interval around ∆m = 0.
In the text

	Fig. C.1 Ratio Rλ of the observed over the ISM-reddened photospheric flux versus λ. The ratios are computed from the spectra binned in equal intervals of 95Å and σ-clipped.
In the text

Fig. D.1 τλ as a function of λ for all the 20 dips analysed in Sect.4. The name of the star and the epoch number are given on the secondary axis. The observed values are shown by black dots with error bars. The red coloured lines show the results for the grain size distribution that best fit the data in the wavelength interval 400-900 nm for the extinction opacities, each normalised to the observed value at 550 nm. The values of amin, amax, q are given on top of each panel; the corresponding grain size distribution is displayed in the insert. The vertical grey dotted lines represent the boundaries of λ = 400, 900 nm. Three representative cases are also shown in the main paper.

In the text