© ESO 2020

1 Introduction

Z CMa is a binary intermediate-mass star known for its strong photometric and spectroscopic variability (Covino et al. 1984; Shevchenko et al. 1999). The two components are separated by 0.1″, or about 100 au (Barth et al. 1994), for a distance of 1 kpc (Shevchenko et al. 1999; Kaltcheva & Hilditch 2000). The SE companion is classified as a FU Orionis (FUor) object (Hartmann et al. 1989; Hessman et al. 1991), while the NW component is an embedded, variable intermediate-mass star, initially labeled as an infrared companion (Koresko et al. 1991). The NW component was detected later in the optical (Barth et al. 1994; Thiebaut et al. 1995) and classified as a 16 M_⊙, B0 III star based on its de-reddened photometric colors (van den Ancker et al. 2004), although its radius and luminosity are disputed (e.g., Monnier et al. 2005). A B0 spectral type would place the NW component of Z CMa among the earliest-type stars known to have a disk (e.g., see the reviews by Zinnecker & Yorke 2007; Beltrán & de Wit 2016) and to undergo variable accretion in a way not dissimilar to lower-mass objects, although disks similar to those of solar-type stars are also increasingly found around massive stars (e.g., Bik & Thi 2004; Alonso-Albi et al. 2009; Kraus et al. 2010; Fedriani et al. 2020).

Z CMa has a complex light curve with anomalous variability (Shevchenko et al. 1999). Highly complex, short-timescale, aperiodic variations were confirmed by MOST (Microvariability and Oscillations of Stars; Siwak et al. 2013). Although historically the main variability was attributed to FUor, further observations, including polarization and scattered light, have suggested that the NW component is responsible for the most dramatic changes (Lamzin et al. 1998; Szeifert et al. 2010) and also dominates the polarization observed in the system (Fischer et al. 1998). The photometric evolution of Z CMa NW is uncertain since the first optical observations. Early works focussed on the FUor, but there are several reports of anomalous behavior that may be related to the NW component. Z CMa was always presented as an atypical FUor, in part for its long rising time (Hartmann et al. 1989; Hessman et al. 1991; Hartmann & Kenyon 1996), but also because of its bumpy light curve. Several authors found the IR object to be brighter than expected in the optical since the 1990s, although Thiebaut et al. (1995) argued that this could be due to scattered light of the central object on the wall of a low-inclination bipolar cavity. The same conclusion of the NW component being optically brighter was reached by Barth et al. (1994) from observationsobtained three years after Thiebaut et al. (1995). Subsequent interferometric observations revealed a ~4 au dusty ring around Z CMa NW (Monnier et al. 2005) and a larger mid-infrared structure (~68 mas × 41 mas, ~40–70 au) dominated by Z CMa NW but extended toward the FUor companion (Monnier et al. 2009). The FUor companion also has a ring that is similar in size (Millan-Gabet et al. 2006).

The classification of the FUor object emphasized its double-peaked absorption lines, a telltale of strongly accreting sources (Hartmann & Kenyon 1996). The spectra revealed a radial velocity ~30 km s⁻¹, and double-peaked absorption lines consistent with a self-luminous disk 6 au in size and with maximum velocity ~120 km s⁻¹ (Hartmann et al. 1989). In addition, Z CMa showed lines with strong P Cygni profiles (Hartmann et al. 1989; van den Ancker et al. 2004)that were enhanced during the episodes of anomalous, bumpy variability (Hessman et al. 1991). Covino et al. (1984) pointed out significant spectral variations since observations began in the 1920s, and several authors have suggested that the origin of the strong line emission lies in the NW component (Garcia et al. 1999; Benisty et al. 2010; Szeifert et al. 2010; Bonnefoy et al. 2017).

There is little information on the accretion mechanisms and accretion evolution in intermediate-mass stars. Unlike T Tauri stars, Herbig Be (HBe) stars are not expected to have magnetospheric accretion as their magnetic fields are too weak (Alecian et al. 2013) as well as very small magnetospheres (Cauley & Johns-Krull 2015). This would lead to a different accretion mechanism compared to solar-type stars (Cauley & Johns-Krull 2016), such as a boundary layer (Popham et al. 1993; Vink et al. 2002; Eisner et al. 2004; Mendigutía et al. 2011; Wichittanakom et al. 2020), a hot inner disk scenario (e.g., Fairlamb et al. 2015; Mendigutía et al. 2015), or non-axisymmetric radial flows (Mendigutía et al. 2017). Therefore, exploring the small-scale accretion structures is especially important in objects such as Z CMa NW.

A new outburst started in January 2008 (Grankin & Artemenko 2009), ending in October 2009. High spatial resolution observations confirmed that the intermediate-mass Z CMa NW was responsible for the increased brightness, and also resolved the near-infrared (near-IR) spectra of both components, revealing the expected characteristic lines of FUor objects for the SE companion, and a large number of emission lines reminiscent of EXor variables for Z CMa NW (Bonnefoy et al. 2017). An enhanced bipolar wind with velocities up to 700 km s⁻¹ was also detected during the outburst (Benisty et al. 2010; Szeifert et al. 2010), which is a further sign that variability episodes are related to increased accretion.

In this paper, we present the spectroscopic analysis of the 2008-2009 outburst and return to quiescence of Z CMa NW, exploring the physical structure and properties in the wind and the innermost disk with the aim of understanding accretion in high-mass stars. Observations, data reduction, and emission line classification are presented in Sect. 2. The analysis of the emission lines is presented in Sect. 3, while the implications for the outbursts of massive versus low-mass stars are discussed in Sect. 4. Our results are summarized in Sect. 5.

2 Observations and data reduction

2.1 Spectroscopic data

Two series of spectra were obtained for Z CMa at the Observatoire de Haute-Provence (OHP) using the SOPHIE spectrograph (Perruchot et al. 2008). The first set of two spectra was obtained on 2008 December 12 and 19 during the outburst. An additional set of three spectra was obtained with the same setup from 2009 October 16 to 19 as the system had returned to quiescence. OHP/SOPHIE spectra are automatically reduced online at the telescope to provide a continuum normalized 1D spectrum ranging from 382 to 693 nm at a spectral resolution of 39 000. With an integration time of 1 h, and depending on weather conditions, the spectra have a signal-to-noise ratio (S/N) ranging from 40 to 150 at 600 nm. Another series of spectra were obtained with the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS; Donati 2003) in spectroscopic mode at the Canada-France-Hawaii Telescope (CFHT) from 2008 January 8 to 14 as the system was reaching the peak of the outburst. A total of six spectra were obtained over seven nights, covering from 370 to 1000 nm at a spectral resolution of 65 000. An exposure time of 4800 s yielded a S/N of 500–600 at 600 nm. Spectra were reduced with the Libre-Esprit software (Donati et al. 1997) to provide continuum normalized 1D spectra. Table 1 summarizes the available spectroscopic data.

Our data are not spatially resolved, so the spectra contain emission from both Z CMa NW and from the FUor companion. We assume that the FUor contribution can be neglected during outburst because the luminosity is at least 8 times higher during outburst as in quiescence, the FUor outburst seems to have ended (or substantially weakened) after 1995 (see Appendix A), and thewind component has been estimated to be at least 20 times stronger for Z CMa NW in outburst as for the FUor (Antoniucci et al. 2016). The quiescence spectra do not reveal any features typical of a FUor. Even though there is a larger relative contribution from the FUor companion, the spectra resemble those of a typical, strongly accreting intermediate-mass star. The FUor features such as double-peaked absorption components that were present before (e.g., Hartmann et al. 1989; Hartmann & Kenyon 1996) are no longer observed, indicating the evolution of the FUor outburst in the last ~20 yr. We thus suspect that the quiescence data may also be dominated by Z CMa NW, although a significant contribution from the companion cannot be excluded.

2.2 Photometric observations

Z CMa was observed by the All Sky Automated Survey (ASAS; Pojmanksy 2003) between December 2000 and December 2009 in the V filter. The photometry is calculated using five different apertures. We chose the smallest (2 pix, where the native pixel is 15″¹) to minimize the amount of nebular emission included, although the differences between apertures are negligible, especially considering the variability range. Only the quality A data are used, containing 637 photometric data points spread over nine years (Fig. 1). The ASAS data include the full outburst, rising from quiescence levels by more than two magnitudes and coming back to quiescence, and offer a good baseline to explore the status of the object at the time of the spectroscopic observations. The ASAS data also reveal that, although the 2009 outburst is significantly longer and very bright, the source has suffered several bursts, rising to nearly the same magnitude for shorter times, since their records started.

The historical data from the American Association of Variable Stars Observers (AAVSO; see Appendix A) show that the smooth FUor outburst profile is punctured by brief bursts (similar to the bumps cited by Hessman et al. 1991) since the early 1980s. The overall shape of the light curve suggests that the NW component has become increasingly active and likely brighter and less extincted in the last ~40 yr. The FUor outburst had faded substantially by 1995. More recent AAVSO observations show that the bumpy behavior recorded by ASAS continues until the present day (see Fig. A.1), so that the outbursts described in this work are recurrent.

Analysis of the ASAS data does not reveal any significant periodic signature on short (days–months) or long (months–years) timescales. A low-significance 249 d quasi-period can be derived from the whole V data (ASAS plus the AAVSO, see Appendix A), suggesting that outbursts happen roughly every 249 days. Such a period is similar to that of a body in Keplerian rotation at a distance 2.0 × (M_*/16 M_⊙^)1∕3 au. Regular perturbations of the disk by companions may trigger outbursts (Lodato & Clarke 2004), although the results have low significance and need follow up for confirmation. The ASAS data also show a brief dip or occultation event that happened just before the first spectrum was acquired. Although the dip was not completely over at the time of the first observation, we do not find any significant spectroscopic signature linked to this event.

Fig. 1 Full light curve (left) from ASAS and zoom around the times of the spectroscopic observations (right). The times at which thespectra were taken are indicated as blue for outburst spectra and orange for quiescence spectra. A complete light curve is given in Fig. A.1.

2.3 Emission line classification

The emission lines were visually selected and classified using the NIST database (Ralchenko et al. 2010; Kramida et al. 2018) and the list of lines observed in Z CMa and other outbursting objects (van den Ancker et al. 2004; Sicilia-Aguilar et al. 2012, 2015, 2017). The narrow neutral metallic emission lines in quiescence reveal a radial velocity ~27 ± 3 km s⁻¹, which is consistent with previous estimates (~30 km s⁻¹; Hartmann et al. 1989). This radial velocity was used to revise the classification of the remaining emission lines, in particular, those that were weak or not observed in other young stars. In total, we estimated more than 1000 emission lines in the optical outburst spectra, but because of low S/N and blends, we only resolved 497 of these, among which 26 lines cannot be identified in the literature. The complete line list is presented in Appendix B (Table B.1). The table contains all identified lines plus blends, together with excitation potentials, transition probabilities, profile types, strength, and whether they have been observed in outburst and/or quiescence.

The lines include typical H, He, and Ca II emission plus a large number of Fe I, Fe II, Ti I, Ti II, Si II, and other metallic lines commonly found in young stars (Joy 1945; Hamann & Persson 1992). Nearly all lines correspond to permitted transitions and the only forbidden lines that can be clearly identified² are the [O I] lines at λλ6300, 6363 Å and the [S II] lines at λλ6717, 6731 Å. Manyof these are common to both quiescence and outburst, but their profiles are different. Several examples of lines with various profiles are shown in Fig. 2. Although most of the lines are broad in both quiescence and outburst, the outburst lines are significantly broader; some such as Hα have line wings up to ±1000 km s⁻¹. A P Cygni profile is most common among the strong, high-energy lines, such as H and Fe II, in both outburst and quiescence; 89 lines have this kind of profile. Very high energy lines, such as He I and Si II, have extreme P Cygni profiles in outburst, which are entirely dominated by blueshifted absorption and with a very small emission component (see Fig. 2), being thus the only lines that show the typical profiles of FUor objects. No such lines are observed in quiescence, so they are likely related to the outbursting Z CMa NW and not to the FUor. The blueshifted absorptions become much weaker (or disappear completely, especially inthe most energetic lines) in the quiescence spectra. This suggests a powerful hot wind developing during outburst that is likely related to the increased accretion episode. Most of the outburst blueshifted absorption features have several components and show a very high degree of variability, as we discuss in Sect. 3.2.

We detected many low-energy, neutral metallic lines in emission during both quiescence and outburst, especially Fe I. In quiescence, most of the Fe I, Na I, and Ni I lines are narrow (FWHM < 30 km s⁻¹) or absent, while in outburst, 91 of these lines show a box-like or double-peaked profile as expected for Keplerian disk emission (Horne & Marsh 1986) and similar to the disk emission profiles observed in HAeBe stars in CO (e.g., Hein Bertelsen et al. 2016). The double-peaked profiles are very similar to those observed by Hartmann et al. (1989) in FUor objects, but they are in emission instead of in absorption. Although viscous dissipation often dominates disk heating in strongly accreting systems (producing a disk that is hotter in the midplane than in the surface; Hartmann et al. 1989), this is a signature of a temperature inversion in the Z CMa NW disk, which remains hotter in its upper layers as seen in irradiated disks (Calvet et al. 1992). Interferometry and polarimetry have also revealed structures that are consistent with disks in some HBe stars (Eisner et al. 2004; Ababakr et al. 2017), which could be similar to what we trace in the neutral metallic lines. There are differences in the disk-like profiles, considering line strength (or transition probability) and some asymmetries in the blue- versus redshifted parts. The former is expected since lines with various strengths saturate at different heights over the disk and stronger lines can be produced over larger, less well-defined regions (and may thus include other components, including low-velocity ones) compared to weaker lines, resulting in box-like rather than double-peak profiles (Ferguson 1997). The latter could be related to asymmetries in the disk and/or wind- or accretion-related absorption components, which is discussed in detail in Sects. 3.4 and 3.5.

Redshifted absorptions are typically a clear indication of infall in the spectra of young stars. We did not observe clear redshifted absorptions in the Z CMa spectra, although there are several lines with redshifted asymmetries that could be consistent with weak absorption components at various velocities. Redshifted absorption asymmetries were observed at ~200 km s⁻¹ in the higher Balmer lines, IR H lines, the O I line at 8446 Å and the Ca II IR triplet (Fig. 2). As we discuss in Sect. 3.6, Fe I lines also have weak redshifted asymmetries in outburst and quiescence that could be related to absorption components, although the velocities are significantly lower.

Regarding forbidden line emission, the [O I] line is the most evident because of its strength and lack of blends with other features. The [O I] lines have two components with typical shock profiles (Fig. 3): one component slightly blueshifted with respect to the source velocity and a second, high-velocity component at approximately −420 km s⁻¹. The high-velocity component is consistent with the high-velocity peak of the jet identified by Whelan et al. (2010), who also detected the SE source jet at velocities −300 to −100 km s⁻¹ (not evidentin our spectra). [S II] emission at 6731 and 6717 Å is also observed. Its high-velocity component at −400 km s⁻¹ is detected in quiescence, but it is negligible in outburst, which could be a contrast effect combined with blending with other broad lines in the outburst spectrum. Additional [N II] emission and some [Fe II] lines may be also present, but because they are heavily blended with other lines (see Table B.1), it is hard to confirm their presence and to reveal any details about their velocity components.

The Li I resonance line, which is characteristic of young T Tauri stars (White & Basri 2003), is also detected (Fig. 4), albeit mostly in emission, as observed for V1118 Ori in its 2005 outburst (Herbig 2008). In outburst, the line has a profile that is very similar to that of EX Lupi in outburst with broad line wings in emission and a shallow absorption at the stellar rest velocity (Sicilia-Aguilar et al. 2012); this could suggest a similar origin. It is nevertheless very hard to produce Li I emission under normal conditions and abundances (e.g., Shore & De Gennaro Aquino 2020), so we explored other possibilities. The profile is clearly different from disk- or box-like profiles, which have distinct peaks or, for box-like lines, a central peak instead of a central absorption. As in the case of EX Lupi, Li I is the only line with this type of profile. There is potential contamination by nearby Fe I and Fe II lines and by the high-velocity component of the [S II] line at 6717 Å, which is likely the cause of the emission feature observed in quiescence, although the wavelength of the weak absorption feature is more consistent with Li I. It is unlikely that [S II] is responsible for the feature observed at 6708 Å during outburst because the high-velocity component of the 6731 Å [S II] line is negligible. The Fe II line is clearly moreenergetic than the typical Fe II emission observed, and thus unlikely to contribute. The Fe I line is only slightlymore energetic than observed Fe I emission lines, so a contribution is hard to rule out, although it would be blueshifted. In addition, a V I line could be also present at 6708.07 Å and has similar energetics to the V I line observed at 6643.786 Å, which has a disk-like profile in outburst and an emission profile in quiescence. Nevertheless, the lack of consistency between the profiles of the tentative V I identifications makes the association uncertain. Li I absorption is usually not seen in stars earlier than F0 (Zappala 1972), which could hint at an origin in very hot and dense circumstellar material or inner hot disk atmosphere for the small absorption feature.

Finally, a handful of absorption lines are observed, mostly in quiescence, but some of these lines are also observed in outburst (e.g., 5781, 5798, 6614 Å), which are likely diffuse interstellar bands (DIBs; Herbig 1995; Hobbs et al. 2008). There is no clear identification of photospheric absorption lines in outburst nor quiescence.

Fig. 2 Examples of lines observed for Z CMa. The lines are normalized to the local continuum. In this and the following figures, the quiescence data are plotted in various shades of orange, while the outburst spectra are shown in blue. The dashed line shows the radial velocity of the system (27 km s−1) and the dotted line indicates the continuum level. First row: some H Balmer lines and an IR H line. The quiescence spectra of Hα and Hβ are scaled by × 0.5 for better display. Second row: some of the Fe I lines with disk-like profiles and a K I resonance line. Third and fourth rows: some lines with various types of P Cygni profiles. The quiescence spectra do not cover the range beyond 6940 Å.

Fig. 3 Upper: forbidden [O I] emission. Lower: forbidden [S II] emission. The dotted lines indicate the continuum level, and the dashed lines represent the radial velocity of 27 km s−1. The data from the outburst are plotted in various shades of blue, while the quiescence data are shown in orange. The extra lines observed in the 6363 Å spectrum correspond to unrelated Fe I and Fe II lines, and the strong absorption is due to Fe II 6369 Å. The narrow absorption lines are telluric.

Fig. 4 Zoom around the Li I 6708 Å line. The dotted black line indicates the continuum level, the dashed black lines indicate the Li I transitions (6707.76 and 6707.91 Å), the vertical red dotted lines indicate the Fe I 6707.43 Å and Fe II 6708.88 Å lines, and the green dotted line indicates the V I line at 6708.07 Å; all of these are shifted to match a radial velocity of 27 km s−1. The emission feature in quiescence is likely due to the high-velocity component of [S II].

3 Analysis

In this section, we constrain the properties of the various components of the star-disk system using the emission and absorption lines. First, we revise the stellar parameters and accretion in Sect. 3.1, followed by a discussion of the wind components in Sect. 3.2. Then, we use several methods to extract information about the physical conditions and velocities from the observed emission lines. The details of the methods are given in Appendix C. There are many uncertainties in the structure of a complex object such as Z CMa NW, so we use various techniques that work under different assumptions. These include Saha’s equation, which assumes local thermodynamic equilibrium (LTE), and ratios from lines from the same upper level, which do not depend on ionization equilibrium (Beristain et al. 1998). We finally use velocity brightness decomposition for disk-like profiles (Acke & van den Ancker 2006) to introduce further constraints derived from the line profile on the velocity, temperature, and density structure of the system, as well as to reveal where the observations depart from the simplified models. Blended lines and those that may belong to several species and/or transitions are excluded from the analysis. In addition, non-LTE effects and line pumping may affect line ratios; for instance, due to strong UV emission lines directly feeding the upper level of a line. Lacking UV data, but knowing that lines in the UV, as in the optical, may be very broad, we used the exhaustive line list from Herczeg et al. (2005) for the young star RU Lupi to exclude any line that could be pumped by UV lines with ±1000 km s⁻¹ wings.

Fig. 5 Extinction in the Gaia G band and Gaia BP-RP around a 0.2 degree radius field toward Z CMa. The rise of extinction is fully consistent with the location of the Z CMa cloud at 1 kpc.

3.1 Constraining the stellar properties and accretion

Because of the disagreements regarding the spectral type and quiescence luminosity of the star (van den Ancker et al. 2004; Monnier et al. 2005), our first step was to revise these assumptions. First, we used Gaia DR2 data to revise the distance. Z CMa itself is not a reliable Gaia source because it is too bright, is a binary, and is surrounded by nebulosity. Nevertheless, a rise in extinction at a certain distance can be used to identify the location of its cloud (Green et al. 2015; Sicilia-Aguilar et al. 2017). Examining the Gaia DR2 stellar extinction for stars 0.2degrees around Z CMa reveals a clear rise around 1 kpc (Fig. 5), confirming the previously assumed distance to the Z CMa cloud, and thus the stellar luminosity. For a quiescence luminosity of 2400 L_⊙ (Thiebaut et al. 1995) and assuming that the system has an age in the 0.3–3 Myr range, the CMD 3.3 stellar evolutionary tracks³ (Bressan et al. 2012; Marigo et al. 2017; Pastorelli et al. 2019) suggest a mass between 6 and 8 M_⊙. The temperaturewould be consistent with a B-type star (16000–22 000 K), albeit later than B0. The star could be colder (≿10 000 K) but similarly massive if it is younger than 0.3 Myr. Assuming that half of the observed luminosity in quiescence is due to accretion, we still get a mass between 6 M_⊙ and 8 M_⊙ and a temperature between 13 000 and 20 000 K (for ages 0.3–2 Myr) or 8000 and 13 000 K (for ages younger than 0.3 Myr). So while the object is definitely an intermediate-mass star, the mass is likely lower and the spectral type is uncertain.

We thus adopted the quiescence luminosity of L_* = 2400 L_⊙ and give all results as a function of stellar mass. An effective temperature of 29 200 K (for a B0 star) results in a stellar radius of about 1.93 R_⊙ in quiescence, but the radius would be larger if T_eff is lower, especially if the star is very young. During outburst, the luminosity increased by a factor of 8. Assuming a blackbody-like scenario, the effective temperature for this increased luminosity could reach 49 100 K if the emission originated from the star. Nevertheless, we expect a range of temperatures for the star and innermost disk and accretion structures. Despite all these uncertainties, if we assume that the extra energy due to accretion during outburst is dissipated in the innermost regions (star and innermost disk), the temperature of the regions outside this area depends on the total luminosity rather than on the actual temperatureand radius of the star or the size of the emitting region. The temperature required for silicate dust sublimation (~1500 K) is reached at ~3.4 au during quiescence (consistent with the dusty ring at 4 au resolved by Monnier et al. 2005) and at ~9.6 au in outburst. These distances are lower limits to the inner dusty rim in the disk; its radius could be larger if there is in situ disk heating. Large amount of silicates may have been vaporized in outburst. The temperatures at 1 au are ~2800 K in quiescence and ~4700 K in outburst. Orbital velocities are 120$\times {(M_{*}/16M_{\odot })}^{1/2}$ km s⁻¹ at 1 au, and 38$\times {(M_{*}/16M_{\odot })}^{1/2}$ km s⁻¹ at the outburst dust sublimation radius (9.6 au).

The accretion rate (Ṁ) can be estimated from the increase in luminosity, as the luminosity in outburst is significantly higher than the quiescence luminosity of the object. For T Tauri stars with magnetospheric accretion, the accretion luminosity is $$L_{\text{acc}}=\frac{G{M}_{*}\dot{M}}{R_{*}}(1-\frac{R_{*}}{R_{\text{in}}}),$$(1)

where M_* and R_* are the stellar mass and radius, G is the gravitational constant, and R_in is the radius from which the material is accreted (usually, the magnetospheric or corotation radius in T Tauri stars; e.g., Gullbring et al. 1998; Sicilia-Aguilar et al. 2010; Fairlamb et al. 2015). We assume that the infall factor (1 − R_*∕R_in) is about 0.5, which can be attained from material infalling from a minimum distance of 2 stellar radii and is also equivalent to the result for a boundary layer (Bertout et al. 1988). We obtain Ṁ ~ 10⁻⁴ × (16 M_⊙∕M_*) M_⊙ yr⁻¹, which is similar to the accretion rates of FUors (Hartmann & Kenyon 1996). In general, B stars are fast rotators (Abt et al. 2002; Müller et al. 2011), which results in small, if any, magnetospheres (Tambovtseva et al. 2016). A magnetosphere extending to several times the stellar radius would result in a lower accretion rate, while detailed boundary layer fits to B stars tend to produce higher accretion estimates (Wichittanakom et al. 2020).

The feasibility of free-falling material feeding accretion can be also estimated looking at the free-fall velocity, v_ff, $$v_{ff}={(\frac{2G{M}_{*}}{R_{*}})}^{1/2}{(1-\frac{R_{*}}{R_{\text{in}}})}^{1/2}.$$(2)

For Z CMa and an infall factor of 0.5, we obtain a free-fall velocity ~1200$\times {(M_{*}/16M_{\odot })}^{1/2}$ km s⁻¹. This is similar to the velocities observed in the Hα line wings (up to ~ 1000 km s⁻¹), but Hα can be strongly affected by Stark broadening. The strongest metallic lines have maximum velocities of the range of 500–600 km s⁻¹ and any redshifted absorption features appear at lower velocities (see Sect. 2.3). Therefore, whether velocities can be explained by free-fall needs further discussion after other possibilities have been considered (Sect. 4).

Some initial constraints on the density can be derived from the estimated accretion rate. Assuming that accreted material moves in at a typical velocity v over one or more accretion channels that cover a fraction f of the stellar surface, we can write $$\dot{M}=4\pi R_{*}^{2}fv\rho =4\pi R_{*}^{2}fv\mu m_{H}n,$$(3)

where ρ is the density of the accreting material, which can be written in terms of the mean atomic weight (μ = 1.36 is a typical value), the mass of the hydrogen atom (m_H), and the number density (n; e.g., Sicilia-Aguilar et al. 2012). Taking a typical velocity v = 500 km s⁻¹, we obtain a density n ~1.6e14∕f cm⁻³. For any small value of f, these densities are about two orders of magnitude larger than observed in outbursting T Tauri stars. Nevertheless, if the material is brought onto the star via an extended structure involving part of the innermost disk, we could have f > 1 and lower densities, as found, for instance, for ASASSN- 13db (Sicilia-Aguilar et al. 2017, but note that ASASSN-13db has magnetospheric accretion).

Temperature estimates also suggest that the energy is deposited over a region larger than the stellar surface. If the accretion luminosity was released in a stellar spot at the stellar surface, the temperature of the spot could be unphysical, reaching (47e3/$\sqrt{f}$) K, where f is the fraction of the stellar surface covered by the spot. Assuming a boundary layer scenario, the boundary layer temperature T_BL can be written as $$T_{\text{BL}}={(\frac{L_{\text{acc}}}{4\pi R_{*}\delta \sigma })}^{1/4},$$(4)

where σ is Stefan’s constant, δ is the scale of the boundary layer, and the rest of symbols retain their previous meanings (Popham et al. 1993; Blondel & Djie 2006; Mendigutía 2020). The value ofδ is calculatedas $$\delta ={\delta }_0{R}_{*}{\left[1-{\left(\frac{{\Omega }_{*}}{{\Omega }_K}\right)}^2\right]}^{-1/3},$$(5)

where δ₀ is a scale parameter, with typical values between 0.01 and 0.5 (Popham et al. 1993; Blondel & Djie 2006); Ω_* is the angular velocity of the star; and Ω_K is the Keplerian angular velocity. For a rotation rate of 25% of the breakup velocity (typical for B stars; Abt et al. 2002), δ₀ = 0.5 would still result in T_BL ~ 56 000 K, therefore the requirement of a structure that extends over several stellar radii remains. This is also consistent with the very extended boundary layers, often hard to distinguish from the disk, expected in objects with very high accretion rates (Popham et al. 1993).

3.2 Wind variability in outburst and quiescence

Variable and complex, strong winds are typical of quiescent early-type stars (Prinja & Howarth 1986). The many emission lines with P Cygni profiles are a sign of a powerful wind that gets stronger and more complex during outburst, which is typical of accretion outbursts and similar to what is observed in previous events (Covino et al. 1984; van den Ancker et al. 2004). At least three absorption components (at −900, −500 and −100 km s⁻¹) were reported by van den Ancker et al. (2004) during the 1999 outburst; the blueshifted absorption extended up to −1000 km s⁻¹ in the Balmer H series. According to Antoniucci et al. (2016), the wind from Z CMa NW is at least 20 times stronger than the FUor wind for velocities −90 to −400 km s⁻¹. This is confirmed by our quiescence spectra, which show much weaker absorption components despite having a higher proportional contribution from the FUor, and justifies our approach of treating the observed outburst lines as coming from Z CMa NW⁴. The analysis of the quiescence wind is thus limited to pointing out the differences with outburst, since the quiescence data have an unknown contribution of Z CMa SE.

In our 2008/2009 data, the absorptions in the strongest Balmer lines are highly saturated and the wind is stronger than in the past, consistent with the 2 mag rise experienced in the 2009 outburst compared, for instance, to the 1 mag rise in 1999. The wind affects preferentially the high-energy lines, which is a sign of high wind temperature. The absorption features have complex and highly variable, high-velocity profiles, up to −800 km s⁻¹ in certain lines. The velocities of the fast and intermediate-velocity winds (−100 to −500 km s⁻¹) are similar to escape velocities in the innermost regions of the system (<0.1–0.3 au), which is consistent with stellar winds; these winds are associated with innermost accretion structures or initially poorly collimated launching regions that become the large-scale collimated jet. Since hotter temperatures are more likely to originate closer to thestar, the high-velocity wind is the best candidate for a stellar or innermost-disk wind, while the low-velocity wind could originate in a disk wind.

Inner disk winds can explain the deep and broad absorptions observed in FUors (Milliner et al. 2019), which are similar to what is observed for Z CMa in the higher energy lines. He I and Si II are peaked toward the fastest wind components at > − 400 km s⁻¹ (see Fig. 2), indicating a very high temperature. He I and energetic lines may come from a more direct view (hotter, denser, saturated) than the Fe II and Ti II lines, which can also explain the differences in wind velocity and maximum redshifted velocities in the emission part of the PCygni profiles. The O I 8446 Å line shows fast and slow components similar to Fe II lines, such as the Fe II 5018/4923 Å multiplet, maybe due to pumping by UV lines (e.g., H₂ C-X0-2Q(10) line; Herczeg et al. 2005). The absorption components observed in the optical lines are also in agreement with the large-scale jet observed around Z CMa NW, which have velocities from −450 to −600 km s⁻¹ (Whelan et al. 2010; Antoniucci et al. 2016). Ti II lines, which typically trace material at lower temperatures and densities, are dominated by the intermediate and slow wind.

The redshifted emission peaks of most strong ionized lines also suggest a slower wind around −20 km s⁻¹. The velocity of the slow wind is comparable to the escape velocity at ~10–30 au, which is beyond the disk region detected in emission lines, but the orbital period at this distance is too long to expect significant day-to-day modulations, unlike what is observed.

Variations in the velocity of the absorption components are observed in all lines with non-saturated PC profiles, including those that have only a low-velocity absorption component, so that all wind components are variable. The day-to-day variations of the wind components could be caused by real variability, rotation, and/or small-scale occultation events in the innermost regions. Rapid variability in a high-velocity wind component has been also observed for ASASSN-13db in outburst, caused by rotation of an accretion-powered, non-axisymmetric wind (Sicilia-Aguilar et al. 2017). A difference in inclination angle can cause a variation in the absorption versus emission parts of the wind (e.g., Milliner et al. 2019), therefore we may observe material at multiple inclinations depending on the temperature structure.

The variability of the wind shows a clear time evolution. This is particularly evident for the six spectra taken on consecutive nights between 2009 January 7 and 13 (Fig. 6). To examine these spectra, the various absorption components were locally fitted by Gaussians, considering their centers as the representative velocity of each individual wind component. This fit is a simple approximation, not taking into account other effects (e.g., optical depth, radiative transfer). Figure 7 shows that the Fe II and Ti II absorption profiles vary in parallel, even if their central velocities differ, which may arise from different line opacities (Sicilia-Aguilar et al. 2015). If we consider the Fe II 5316 Å and Ti II 4468 Å lines as examples, the Fe II very fast, broad wind component (with velocities between −300/−380 km s⁻¹) is the most variable and does not display a well-defined trend (besides not being observed in the lower energy Ti II lines). The intermediate-velocity wind component at −95/−120 km s⁻¹ and the slower wind component at −47/−55 km s⁻¹ are correlated for both lines (correlation coefficient 0.90, with false-alarm-probability 0.04 using a Spearman rank test) and vary smoothly with time, revealing the development and evolution of each individual component during the six observations.

Two scenarios are possible. The first would be rotational modulation of the wind (either in a stellar or disk wind), as it has seen in the multiple system GW Ori (Fang et al. 2014) and in the outbursting star FU Orionis (Powell et al. 2012) and ASASSN-13db (Sicilia-Aguilar et al. 2017). The wind modulation in GW Ori is caused by a companion, but any process resulting in a non-axisymmetric wind associated with either the star or the disk can produce the desired effect. The wind modulation in FU Orionis is periodic, which cannot be confirmed in this work owing to the few spectra available. The second scenario would result from the propagation of a clump of matter along the wind. The radial velocities could vary owing to changes in their angle with respectto the line of sight, as also occurs in the rotation scenario, or to changes in the physical velocity (e.g., in an accelerating wind). Both scenarios produce a smooth velocity evolution (as it is observed), but an accelerated wind would move toward larger velocities, while rotation would result in changes in the angle with respect to the line of sight. The variability ranges proportional to their velocity; the fast component changes by ~25 km s⁻¹ and the slow wind component changes by ~8 km s⁻¹ during thesesix day (see Fig. 7). This is consistent with rotation and projection effects as observed in rapidly variable, winding winds in certain B-type stars (Prinja & Howarth 1986, 1988), where a rotating/spiraling, non-axisymmetric clumpy wind, moving not only radially, but also azimuthally, causes differences in opacity along the line of sight and thus rapid velocity variations. The trend extends to the very fast wind component observed in Fe II, although the error in its central velocity is large because the absorption features are very broad. Testing whether the system shows any periodic or quasi-periodic behavior in the wind velocity and strength can help to distinguish these scenarios.

The conclusions of this section are that Z CMa NW develops a hot, fast wind in outburst, which is likely of stellar or innermost disk origin, as well as colder components that could be related to the disk. The wind is variable on timescales of days, and the intermediate and strong, slow wind components are clearly correlated, betraying a small-scale structure and a clumpy nature in a rotating, non-axisymmetric wind.

Fig. 6 Examples of the time variability observed for the Ti II 4468 Å line (left, blue), the Fe II 4923 Å line (middle, green), and the Fe II 5316 Å line (right, red). In each case, the 8 subfigures display the individual spectra as they were obtained in time. The dotted lines indicate the continuum level and the radial velocity of the source.

Fig. 7 Day-to-day velocity variation of the wind components in 2009 January, classified according to the central velocities as very fast (top, observed clearly only for Fe II), fast (middle), and slow (lower). Results for the Fe II 5316 (red), Fe II 4923 (green), and Ti II 4468 Å (blue) lines, which have good S/N and no signs of contaminations, are shown. The velocities are given as measured, without correcting for the radial velocity of the system. The correlations between the components observed for both intermediate velocity and slow wind lines, despite some shifts in velocity, are shown. The fast wind component appears to split into several components for some dates for the Fe II 4923 Å.

3.3 Saha’s equation constraints on temperature and density in outburst and quiescence

We can use the observed emission lines to put constraints on the density and temperature of the emitting material, assuming LTE. We use Saha’s and Boltzmann’s equations (Appendix C.1) to derive expected intensities for a range of temperatures and densities and compare these with the observations. For velocity-resolved lines, it is possible to distinguish various gaseous structures considering their velocities. For Z CMa NW, where different lines are related to different physical structures (e.g. winds, hot inner disk, accretion structures), it is crucial to classify the lines according to their profiles. Most Fe II lines have P Cygni profiles, while Fe I lines have disk- or box-like profiles, a sign that they are not produced by the same structure and cannot be treated together. Many lines are saturated and probably originate in optically thick regions, so a proper fit would require full radiative transfer. As a result, the only cases for which we can apply this method are weak Fe I and Fe II lines with similar box- or disk-like profiles (see Table 2). There are only two clean Fe II lines with box-like profiles, λλ 5627.497, 6084.103Å. These lines have transition probabilities on the order of 1e3 s⁻¹ and are very weak⁵. Among the clean Fe I lines, only four have low enough transition probabilities for the assumption of optically thin lines to be justified.

The lines were extracted from the continuum and scaled to the same arbitrary-unit scale. Considering a reasonable range of temperatures and electron densities, we can estimate the relative line intensity expected depending on the physical conditions. Leaving the total density of Fe I as an unknown, the relative strength of Fe II lines with respect to Fe I lines is given by Saha’s equation (Eq. (C.1)). Then, for each transition we consider the atomic parameters so that the total line strength is given by the transition probability (A_ki) divided by the wavelength and multiplied by the Boltzmann distribution (Eq. (C.2)). We estimated theoretical line ratios for temperatures from 1000 to 12 500 K and electron densities 10³ –10¹⁵ cm⁻³. Saha’s equation has an intrinsic degeneracy between the electron density and the temperature, which means that higher temperatures with higher electron densities may reproduce the same line ratios, but we add further constraints in the following sections. The spectra are not flux-calibrated, and the theoretical intensity is calculated as a ratio of the intensities depending on the ratio of ionized to neutral elements, so the lack of scale invalidates the χ² statistics (e.g., Andrae et al. 2010). Therefore, the best-fitting density is obtained by minimizing the standard deviation in the ratio of the observed to theoretical line intensities, taking the median of the ratio of the observed to theoretical lines as a scaling factor. The contribution of each line is weighted taking into account the uncertainties in its atomic parameters using the NIST data qualityflag. To assess the uncertainties in the results, the exercise was repeated for each one of the eight outburst observationsand for the average outburst spectrum. There are no significant differences between dates, thus the variations are due to noise.

The exercise can be repeated for the lines observed in quiescence, remembering that there is no evidence that quiescence lines originate in the disk. The data from 2009-10-18 suggest higher temperatures and densities than for the other quiescence spectra, which could be related to an increased accretion rate in either of the two Z CMa components. Although the Hα emission line is strong on that particular day, the noise level is also higher and there are no other significant differences between all three quiescence spectra. Therefore, random noise cannot be fully excluded and the results from the averaged quiescence spectra maybe more representative of the typical quiescence state.

The explored range of electron densities and temperatures and the best-fitting values are shown in Fig. 8 and Table 3. One warning regarding the interpretation of these results is that the deviations from the best-fit models for different groups of lines are not random. While neutral lines tend to favor a relatively low temperature (on the order of 2000–3000 K), to reproduce the ionized lines at the same level we tend to require higher temperatures and densities (~7600 K) and the best fit tends to the lower temperature values because the atomic parameters of the neutral lines have smaller uncertainties. There are no perfect fits and, in particular, the FeI line at 5506 Å and the FeII line at 6084 Å, which has a particularly large uncertainty in its atomic parameters, are usually badly fitted (within a factor of ~10) by any combination of temperatures and electron densities that reproduces the remaining lines. As observed in EX Lupi (Sicilia-Aguilar et al. 2015), this may indicate multiple components with various densities and temperatures, although for Z CMa there are too few weak Fe II lines to attempt a more complex model fit. Despite the uncertainties and degeneracy, the best fits suggest a relatively low temperature for the formation of the observed disk-shaped lines. The constraints on the density are weaker, although the best-fitting values tend toward the lower end as well (see Fig. 8). Although we cannot use the Ti I/Ti II lines for this exercise because all the Ti II lines are heavily affected by the winds, we note that the temperature and density ranges obtained from Fe I/Fe II are consistent with the strongest Ti I lines having rounded or narrow profiles. These lines are very sensitive to temperature, disappearing for temperatures more than 3000 K at the observed best-fit densities. The fact that no disk-like profiles are observed for Ti I lines suggests that they may come from regions in the disk at larger radii (and thus colder) than the Fe I lines. Therefore, this section adds an additional proof of the inner disk emission being extended, which agrees with the range of temperatures required to explain both the Fe I and Fe II emission (from 2000–3000 to 7000 K); it also agrees with the variety of observed disk-like profiles. The narrow quiescence lines may be linked to the accretion structures in quiescence.

Fig. 8 Result of Saha’s equation analysis for the lines with disk-like profiles and low transition probabilities for the observationstaken during the outburst phase (upper panels) and for the same lines as observed during quiescence (lower row). The panels show the data corresponding to individual days in chronological order; last panel: result of fitting the average spectrum. The color scheme reveals the relative variations in the line intensity ratio for each combination of temperature and electron density (such that 0.1 would correspond to 10% average deviations). The white contour indicates the regions for which the relative variations of the line ratio are up to 3 times the minimum value observed, which is considered as our best-fit region; a white star denotes the best fit.

3.4 Line ratio constraints on the disk and wind properties

Here, we follow Beristain et al. (1998) to derive the optical depth of the line-emitting disk material using lines that share the same upper level. This avoids having to deal with the way the upper level is populated, so that the resulting relative intensities depend only on the atomic parameters, temperature, and line opacity. For velocity-resolved lines, we can apply the Sobolev approximation to write down the line opacity, which results in the line strength depending only on the atomic parameters, temperature, and column density-velocity gradient (nd v∕dl, where n is the density, v is the velocity, and l is the spatial scale). The details of the procedure are given in Appendix C.1. We first apply this method to the disk-like emission lines, and then examine the results for the lines with P Cygni profiles.

3.4.1 Line ratio constraints on the inner disk

There are 13 groups of lines from the same upper level among the disk-like lines (Table 4). Nevertheless, only 3 groups remain for which the lines are not blended, are close enough to minimize the effects of the continuum variations, and have high-quality NIST atomic data⁶. A fourth pair, Fe I λλ 8674.7, 8838.4 Å, also satisfies the above conditions, but the presumably stronger line is observed to be weaker. This is probably caused by the 25% uncertainty in the transition probabilities plus the fact that both lines have a very similar strength. Considering the uncertainties, the best line pair is Fe I 5572/5624 Å⁷ (10% uncertainty or class B, compared to 25% uncertainty for the other two pairs). The smaller difference between the transition probabilitiesfor the Fe I 6393/6492 Å pair makes it the least constrained. All three line pairs are also detected in quiescence, but the Fe I 5572/5624 Å pair in quiescence is too weak to be used.

To derive the ratios, we normalized every line individually. First, the global shape of the spectrum has to be taken into account, which for lines that have similar wavelengths can be done assuming a linear fit between the two regions. We imposed a minimum of 5σ emission above the noise level to obtain the ratio, where σ corresponds to the standard deviation of the flux measured in a line-free nearby region, to avoid unphysical values due to noise. We analyzed the outburst spectra individually, together with the averaged spectrum, to estimate the noise and variability. For comparison, we also analyzed the quiescence data in the same manner, which is discussed in Sect. 3.6.

The results are shown in Fig. 9 and continued in Appendix C. All figures show the velocity-by-velocity line ratio, followed by a comparison with models covering a range of temperature and velocity-density gradient going from the conditions where both lines are optically thin, to the strong line saturating, and finally, to the limit where the two lines are saturated (Beristain et al. 1998). Leaving aside the Fe I 5572/5624 Å pair that is close to saturation in the line wings as well as the center, the central part (closer to zero velocity) appears to be optically thinner than the high-velocity wings in theoutburst lines. For a disk scenario, this can be explained if the density decreases toward larger radii (lower velocities), which is reasonable. For the quiescence data, the emission is close to being optically thin in the weak line for both the Fe I 6256/6191 Å and the Fe I 6462/6393 Å pairs. This is significantly different from what has been observed in EX Lupi (Sicilia-Aguilar et al. 2012, 2015) and DR Tau (Beristain et al. 1998), where narrow lines are optically thicker than the broad components and can be traced back to dense, hot spots on the stellar surface⁸, and suggests a disk origin.

We can combine all of the lines to derive a better constraint on temperature and density (see Fig. 10). The temperatureis not too well constrained from this exercise, but most of the lines agree with column density-velocity gradients (nd l∕dv) on the order of 10¹⁰ –10¹² cm⁻³s, higher than in quiescence. The weakness of the lines, the uncertainty of some of the atomic parameters, and the potential breakdown of the Sobolev approximation for lower velocity gradients may cause part of the temperature mismatch between line pairs. If we combine this result with Saha’s equation from the previous section, both converge toward larger densities and lower (~4000–5000 K) temperatures, which is also consistent with the origin in a disk and the values derived from the lines with lower uncertainties.

A limitation of this procedure is that lines with various strengths (and line profiles) trace material originating in slightly different locations, which adds a further uncertainty. The line pairs examined to trace the disk emission during outburst have typical ratios in between the optically thin and optically thick limits. This means that the weak line is usually optically thin, while the strong line is optically thick. The optically thin line probes material deeper into the disk, which explains why lines with low transition probabilities look more disk-like (instead of box-like) than lines with high transition probabilities. In addition, narrow-line components (originating in more compact but relatively hot and dense regions, such as hot spots at the stellar surface or in the accretion shock region) would be stronger for the lines with high transition probability because low transition probability lines would be more susceptible to collisional deexcitation in these environments. Such departures from disk-like profile related to optical depth are also observed in CO lines (e.g., see Panić et al. 2008; Hein Bertelsen et al. 2016; Carmona et al. 2017), contributing to the box-like appearance of optically thick lines. As in such cases, Fe I lines with various optical depths can give information concerning the vertical temperature structure. In the case of the strong lines, the box-like features would not result from optically thick material in a flat disk (e.g., as in the examples from Horne & Marsh 1986), but rather from material coming from various layers and distances in the optically thick, flared disk. In addition, contributions from a disk wind can affect the double-peaked profiles in optically thick lines, and can, together with optical depth, result in box-like profiles (Murray & Chiang 1996, 1997; Ferguson 1997). Wind can quickly overwhelm disk (and magnetospheric) emission due to its larger volume (Tambovtseva et al. 2016), which is not observed here.

Most of the observed disk lines are systematically brighter on their blue side. This blue-red asymmetry could be due to inclination and flaring – the red and the blue parts may appear to have different inclinations; (Hein Bertelsen et al. 2016) – a density enhancement on the blue side, an off-center disk, or contamination by other contributors to the line emission such as accretion structures or the disk wind, which does not need to be uniform. We further examine the location of the emitting material based on the velocities in Sect. 3.5, but the results are consistent. A blue asymmetry, likely related to wind, is also observed, albeit in absorption, in the strongest disk-like lines in FU Ori (Donati et al. 2005).

Column density-velocity gradients can be used to estimate the densities and constrain the values derived from Saha’s equation. Considering that the maximum dv∕dl is given by approximately v∕r where v and r are the Keplerian velocity and radius, the best-fitting column density-velocity gradient of 1e12 cm⁻³s suggests a particle density at least 3e4 cm⁻³. This value is clearly lower than our estimate for the accretion structures (Sect. 3.1), which is consistent taking into account that the disk is likely less dense than the accretion channels, especially at larger radii. Together with the temperature derived from this procedure (~4000–5000 K), this section confirms the agreement of the neutral line emission with an origin in an extended disk structure around the star as well as a density (optical depth), radial differences between the high-velocity disk (line wings, formed closer in), and the lower-density, lower-velocity regions (related to the outer disk).

Fig. 9 Analysis of the line ratio for Fe I transitions with common upper levels (continued in Appendix C). Upper left panels: line emission in outburst and upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. The gray area shows the region for which the line emission is >5σ (in the velocity panels) and the 1σ error in the temperature and density planes. Quiescence lines are often very weak and thus uncertain.

Fig. 10 Result of the line ratio analysis for neutral lines with disk profiles in outburst (left) and narrow profiles in quiescence (right).

3.4.2 Line ratio constraints on the wind absorption and redshifted emission components

The same exercise of line ratios was done for the pairs of lines with wind signatures that share an upper common level. We find a total of 15 pairs or triplets of lines with common upper levels and P Cygni profiles (see Table 5), but only 5 pairs remain if we exclude pairs. For the line ratio analysis, we checked both the line ratio in the emission part (measured as positive over the continuum) and in the absorption part (measured as negative under the continuum). The velocity space was divided in the region in which both lines have significant (>5σ) emission and the regions in which both lines have significant absorption. The ratio was only calculated where both lines were either in emission or in absorption. The wind was separated in two components, the fast and the slow wind, but since the emission part of the lines is very broad, it is very likely that the slow wind is contaminated by emission.This explains the unphysical ratios observed at some velocities. The line pairs show no significant absorption during quiescence, a signature of the wind becoming optically thin, as expected from an accretion-powered wind.

Figure 11 shows the results. The wind components cover a large range of temperatures and densities. We can put togetherthe column density-velocity gradient and temperature ranges favored by each line ratio to get a better constraint for both the emission and fast wind (Fig. 12). Since the wind absorptions are very weak for all the line pairs, derived wind properties are strongly affected by the noise. Weak lines with low transition probabilities have very weak or even absent absorption, a signature of being optically thin. This constrains the maximum density allowed for a given wind temperature and the range of densities that this method can probe. For the emission component of the line, both the Fe II line pairs and Ti II line pairs are consistent with a temperature in the range 6000–7000 K, although the Ti II lines favor a column density-velocity gradient value three orders of magnitude lower. Given that the line profiles for Fe II and Ti II are not identical, differences in the physical location are expected. Only the two Ti II pairs give significant results for the fast wind component, suggesting a much lower column density-velocity gradient and a very high temperature (>2e4 K). Since the uncertainties in the atomic parameters for Ti II are ~50%, the final results have uncertainties on the order of 70%, which would be compatible with a lower bound temperature of ~ 10⁴ K. Such large temperatures are more consistent with a stellar, rather than a disk wind, origin. The relative values for the Ti II lines suggest that the column density-velocity gradient is between 1 and 2 orders of magnitude lower for the wind than for the redshifted (potential infall-related) part.

The emission part of the line has a narrow, central component in the weakest ionized lines, probably because the strong lines with large transition probabilities are highly saturated. The lack of narrow components in low transition probability neutral lines is likely related to temperature: the originating accretion structure, wind, or hot spots in the proximities of the star are likely too hot to allow for strong Fe I emission. If the temperatures are low enough to ensure that the metals are all once-ionized and hydrogen remains neutral, the electron density would be approximately 0.1% of the hydrogen number density. Typical redshifted velocities in the metallic lines are around 100 km s⁻¹, which is likely too low to be explained by free fall around an intermediate-mass star, but compatible with, for example, rotation at a few stellar radii. The 10¹³ cm⁻³ s⁻¹ column density velocity gradient constrained by the line ratios would then result in electron densities on the order of 10⁵, which require temperatures of at least 3500 K to avoid Fe I emission (according to Saha’s equation). The temperatures from line ratios are about a factor of 2 higher (see below), explaining the lack of Fe I emission with accretion or wind signatures and the scarcity of Fe II lines with disk profiles.

The temperatures of the redshifted wings derived from the line pairs are around 7000 K, which are low for an origin in the accretion columns in the proximity of a luminous, intermediate-mass star. This suggests that the lines originate in a more distant region, inconsistent with magnetospheric accretion; infall velocity onto an intermediate-mass star would rapidly exceed the 100–200 km s⁻¹ ranges observed in the line wings. A boundary layer scenario also fails to explain the velocities in the inner disk, which are expected to be significantly lower than observed in disk-like lines (Popham et al. 1993), although increased velocities due to turbulence may reach a few 100 km s⁻¹ (Bertout et al. 1988). As in the disk analysis, because optical spectra trace relatively low-energy lines, we have a selection bias toward lower temperatures. This could explain their origin not too close to the star, so we may be indeed missing the hottest, highest velocity emission in most of these line pairs. Much higher velocities, more consistent with infall on an intermediate-mass star, are seen in the very extended wings of the H Balmer lines. A hotter and higher-velocity region could be explored further with spectroscopic observations covering metallic UV lines.

Summarizing, line ratios reveal a highly structured wind and a complex origin for the redshifted emission components. Both cover a broad range of densities and temperatures that can be probed by lines with various excitation energies. The wind temperatures are consistent with a stellar (or innermost disk) origin, while the redshifted emission could originate at a few stellar radii. The wind becomes optically thinner in quiescence, which is consistent with being powered by accretion.

Fig. 11 Analysis of the line ratio for the FeII and TiII lines with PCygni profiles from transitions with common upper levels (continued in Appendix C). Upper left panels: line emission in outburst and upper right panel: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. In the velocity panels, the red, cyan, and gray areas show the regions with significant redshifted emission and slow and fast wind absorption, respectively. The shaded colored areas in the temperature and density planes show the regions that are compatible with the line ratios observed in the emission and the two wind component parts with their 1σ error. See text for details.

Fig. 12 Result of the line ratio analysis ion pairs of lines with a common upper level for the emission and wind absorption parts (denoted as “wind” in the legend). Some of the Ti II lines have sharp cuts in the n dv∕d l vs. T plane due tounphysical values that arise in regions of low S/N (see Appendix C). The results have to be regarded with care when approaching such regions.

3.5 Disk structure derived from the velocity analysis

We used the brightness decomposition technique from Acke & van den Ancker (2006) to study the velocity structure of the disk-like lines. The details of the procedure are given in Appendix C.2. For lines with disk-like profiles, the highest velocities are caused by material in the innermost disk. Starting with the highest velocity, we fit a ring that can reproduce the observed flux and subtracted it from the line profile, repeating the procedure until the line is entirely fitted or we reach a velocity in which emission from narrow, non-disk components dominates. The emission scales as a blackbody with temperature given by a power law of the radius. The final emission is obtained by summing all rings. All the observed lines are best fit using the temperature power law for a flared, irradiated disk model (D’Alessio et al. 1999). Given that the accretion during the outburst is very high⁹, to account for an irradiated disk at few astronomical units we need the bulk of the accretion energy to be released closer to the star than the region that produces the disk-like lines, at a fraction of an astronomic unit or less, or radii smaller than those inferred from the lines with disk profiles, justifying our discussion in Sect. 3.1. Some viscous heating is likely present, but because of its r⁻³ dependence with the disk radius (D’Alessio et al. 1998; Woitke 2015), any viscous heating at the distances at which we observe the double-peaked lines would be strongly suppressed, so that irradiation (by the central star, accretion spots, or even the viscous heating in the innermost part of the disk) dominates.

Our models assume i = 30 degrees, but the inclination is unknown. Since the observed velocity is v_obs = v_keplersini, the radius derived assuming Keplerian rotation is proportional to sin²i. If the stellarmass was lower or the disk was closer to edge-on, the emission zone would move toward smaller radii, constraining even more the locations at which the main source of outburst energy is released. The results are easily scalable, noting that the temperatures, based on irradiation/luminosity, would not change. Nevertheless, the vertical scale height and flaring could be affected if the gravity pull is different. Putting together the effects of inclination and stellar mass, the scaled radius (r_s) for an inclination i and a stellar mass of M_* would be $$r_{\text{s}}=r\frac{{\text{sin}}^2(i)}{{\text{sin}}^2({30}^{°})}\frac{M_{*}}{16M_{\odot }},$$(6)

where r is the value derived for i = 30 degrees and M_* = 16 M_⊙ shown in the figures.

The main limitations of our analysis are the presence of non-Keplerian velocity components (e.g., arising from a boundary layer at low velocities, or from an accretion column or a wind at high velocity) and the fact that different lines (and different velocities) trace material at different depths within the disk. As in previous sections, the results from weak lines are the most reliable. The final result is obtained by fitting the median data, while individual fits to each spectrum are used to determine the uncertainties. For the noisiest lines (Fe I 6358 Å and Ca I 6439 Å), some of the individual datasets were disregarded. A summary of the fitted lines is shown in Table 6.

The results of the velocity decomposition are shown in Fig. 13. The central radial velocity of the lines is ~20–35 km s⁻¹ (average 30 ± 5 km s⁻¹), which is consistent with but not identical to the 27 ± 3 km s⁻¹ radial velocity. The line decomposition algorithm defines the center of the line by symmetry, so that deviations in the line symmetry would alter the line radial velocity. The decomposition reveals that the bulk of disk emission originates at radii r ~ 0.5–3 × (M_*/16 M_⊙) au, and for most lines (especially, those with well-defined profiles and high S/N), the emission peaks at r ~ 1–2 × (M_*/16 M_⊙) au; this is consistent with the radius at which the temperature originated by irradiation is expected to be high enough to result in significant Fe I emission.

The blue and red sides are not symmetric: most of the lines are dominated by the blue component. This suggests that the gaseous disk may be either non-axisymmetric (as it had been observed in the metallic emission lines of EX Lupi in outburst Sicilia-Aguilar et al. 2012) and/or not fully centered around the star. Stronger lines, which are likely to originate in the disk surface layers, tend to be more asymmetric and more variable, suggesting that the outer layers of the disk have a more complex structure than the midplane. The disks of HAeBe stars often have asymmetries in the outer layers at all radii (Benisty et al. 2015; Laws et al. 2020; Kluska et al. 2020), so this is not unexpected. Hot spots could produce extra emission at low velocity, although we do not observe any rotational modulation in the narrow components. Many of the lines have high-velocity wings or tails (in the range of −75 to −90 km s⁻¹) that are clearly distinct from the main disk profile. Although in some cases the wings are seen in both the blue and red part of the spectrum (e.g., Fe I 6400 Å), in most of them the blue high-velocity wing is stronger (or even the only one seen; e.g., Fe I 5624 Å) and certain linesare particularly asymmetric on their blue sides (e.g., Fe I 6750 Å, Ni I 6767 Å). The stronger Fe I lines with higher A_ki have stronger blue wings¹⁰. The high-velocity region could suggest a discontinuous, asymmetric disk structure rotating faster and closer to the star (weak because of its smaller area). Alternatively, the high-velocity tails and blue/red asymmetries could result from absorptions in the range of −60 to −40 km s⁻¹ related to the slow wind component (in the blue) or to infall absorption (in the red, at positive velocities). Nevertheless, the lack of significant variability in these high-velocity tails contradicts what is observed in the wind components, and is more in agreement with an origin in a stable structure (with a longer period) such as a radially and azimuthally asymmetric disk.

Besides the difference in line profile with the transition probability (weak lines look more disk-like than strong ones), there are no significant correlations between excitation potentials, A_ki, and the brightness decomposition. Correlations between the velocity, the full width at half maximum (FWHM), or the peak-to-peak distance, and the atomic parameters are uncertain and noise-dominated. Part of the noise may be caused by the larger variability observed in the line wings of the strong lines that affect their FWHM. Leaving aside the changes that could be attributed to depth and inclination described before, this homogeneity could be a sign that the emitting region is very well constrained and ring-like (Hein Bertelsen et al. 2016).

To summarize, this section confirms the previous results of neutral lines in an irradiated disk-like structure, extending 0.5–3 × (M_*/16 M_⊙) au, and presenting some asymmetries in the surface and the disk midplane as well as radial and azimuthal variations. The double-peaked lines arise from very large radii compared to the stellar radius (and to any reasonable value of the stellar magnetosphere or boundary layer), so these lines cannot provide much information on the accretion mechanism. The temperature ranges inferred from the disk-like lines in the previous sections are not too different from those observed in FUors (i.e., similar to photospheres of F- or G-type stars Hartmann & Kenyon 1996), although this could result from Fe I emission being strongly suppressed at much higher temperatures. Observations of line profiles in the UV during an outburst would be needed to obtain information on the innermost zones around the star and its temperature and velocity structure. As observed in FUors (Hartmann & Kenyon 1996), we also see that the lines at longer wavelengths are less disk-like. This could respond to the same effect observed in FUors where lower temperatures, corresponding to further away distances in the disk, have different line profiles. In this respect, it would be useful to check the profiles of near-IR lines, but unfortunately the available near-IR spectroscopy does not resolve the line profiles of neutral atomic lines (Bonnefoy et al. 2017).

3.6 Return to quiescence

Z CMa NW in quiescence may not be much brighter than its companion, so the three quiescence spectra may contain a combinationof Z CMa SE (likely no longer in outburst; Appendix A) and the quiescent Z CMa NW. In the quiescence spectra, the number of emission lines is still high despite being lower than in outburst, but the line profiles are dominated by redshifted emission components instead of blueshifted absorption. The lines are weaker, but due to the changes in the continuum and in the wind absorption, equivalent widths tend to be higher in quiescence. There is no significant variability in the quiescence data (see Fig. 14, the apparent variability in faint atomic lines is at the noise level). This could be due to the short time baseline (2009-10-16 to 2009-10-19), but since the outburst spectra show variations on consecutive days, the rapid variations observed in outburst clearly subsided in quiescence.

The maximum velocity observed in the Hα red wing decreases from ~1000 km s⁻¹ in outburst to ~750 km s⁻¹ in quiescence, and the wind components faster than −400 km s⁻¹ are absent from the metallic lines in quiescence. High-energy lines such as Si II are absent in quiescence, and in the He I lines the emission component is stronger and the absorption features are shallower in quiescence (Fig. 2). The wind during the quiescence stage is not as clumpy as in outburst, but dominated by a single, more symmetric, wind component instead of several non-axisymmetric components. The intermediate-velocity wind (v_wind > −100∕−200 km s⁻¹) is still observed in the metallic (including Fe II and Ti II), H I, and He I lines in quiescence, but the strong slow wind (v_wind ~−50∕−100 km s⁻¹) disappears in quiescence (Fig. 2). This is not unreasonable if the slow wind originates in the disk and the amount of matter transported through the disk decreases in quiescence, affecting the wind mass and optical depth, as is often found in accretion-driven outbursts. The low-excitation material also becomes less optically thick during quiescence, which causes most of the weak lines to become very faint or disappear altogether (see Fig. 9). The pair Fe I 5624/5572 Å is too weak in quiescence to provide physical results, leaving the temperature largely unconstrained (Fig. 12). Optically thinner circumstellar material in quiescence suggests that the disk accretion rate has substantially decreased after the accretion burst, as has also been observed in CO lines in other outbursting stars (Banzatti & Pontoppidan 2015).

The decrease in optical thickness and column density-velocity gradient in quiescence reveals that, even if the quiescence lines were dominated by the FUor component, the structures accreting onto the FUor are significantly less dense than those of the outbursting Z CMa NW. Differences in inclination may change the observed line widths and depths, but unless accretion covers a significantly larger fraction of the FUor, the accretion rate of the FUor inferred from the densities would be lower, by a couple of orders of magnitude, than the estimated value of 10⁻⁴ M_⊙ yr⁻¹ for Z CMaNW in outburst. Since 10⁻⁴–10⁻⁵ M_⊙ yr⁻¹ is a typical accretion rate for a FUor (Hartmann & Kenyon 1996), this is an additional sign that, as of 2009, Z CMa SE had returned to quiescence, although it could be still a strong accretor. We also checked for further signatures of redshifted absorption, characteristic of nonspherical infall or magnetospheric accretion, in the less-blended strong Ti II lines observed in quiescence (Fig. 15). The asymmetric redshifted side of the line is suggestive of redshifted absorption, but no consistent trend is observed and the main signature of the Ti II quiescence lines is a strong, narrow wind component at −130 km s⁻¹.

The [O I] components at ~0 km s⁻¹ and ~ − 400 km s⁻¹ seem stronger in quiescence owing to the lower continuum levels, but their velocities do not change between quiescence and outburst (see Fig. 3). The high-velocity component becomes visible in [S II] in quiescence. The profiles, which are only clearly resolved in quiescence because of contrast, are typical of Herbig-Haro objects. The high-velocity, forbidden line components could be related to the main wind (~ −400 km s⁻¹) and could correspond to a shock in the wind emission. The second, at ~ 0 km s⁻¹, could be related to a shock along the disk or extended envelope, but the origin is uncertain since Z CMa SE also drives shocks (Bonnefoy et al. 2017).

Many low-excitation lines with box- or disk-like profiles in outburst are detected in quiescence, but they are narrow and do not have double-peaked profiles. Since in quiescence the irradiated disk temperature is expected to decrease rapidly, we would expect the disk emission to move to faster velocities as the disk cools down, but the emitting area (and thus the line strength) would also decrease rapidly. Strong neutral lines are often still strong and have complex profiles, suggesting collisional de-excitation in a still-high-density environment. Some of the strongest Fe I lines have asymmetric profiles (e.g., Fe I 6191 and 6393 Å lines in Fig. 9), dominated by blueshifted emission. This redshifted asymmetry could be indicative of infall or nonspherical accretion, although the velocities are very low (~ −20 km s⁻¹). The typical FWHM of 20–30 km s⁻¹ for the metallic quiescence lines would correspond to Keplerian rotation at ~16 × (M_*/16 M_⊙) au, at which the disk will be too cold in quiescence. Their wings extend out to ± 50 km s⁻¹ and could be related to the innermost disk (3–5 × (M_*/16 M_⊙) au), but line width is very small compared to infall or rotation velocities around an intermediate-mass star. The lower luminosity in quiescence suggests line formation closer to the star. If the lines were associated with small and localized hot spots on the stellar surface, we would expect to see some rotational modulation, unless the star is seen close to pole-on. The lack of rotational modulation in the narrow quiescence emission lines suggests that the material is either at a polar, always-visible spot, or extended in a rather uniform way around the star. This could be expected for unstable accretion (Kurosawa & Romanova 2013) or if the star accreted through an axisymmetric structure or boundary layer. Stellar rotation could account for the line broadening if the star is a slow or highly inclined rotator. Further observations over a more extended period of time would be needed to clarify this issue and, in particular, to determine the longer-term variability during quiescence that may betray its spatial location. Summarizing, quiescence brings a dramatic decrease of mass accretion and wind density, triggering the disappearance of the low-velocity disk wind and non-axisymmetric wind components.

Fig. 13 Brightness profile decomposition for the disk-like lines (continued in Appendix C). For each line, upper panel: normalized, median line emission (blue), the individual data points (colored dots), and the model fit. The black solid line is the best fit to the median data, while the dashed lines are the fits to the individual datasets that are used to compute the errors. The zero velocity (measured with respect to the symmetry axis of the line) is shown as a dotted line. The individual ring fits are shown by thin red lines, and the residuals are shown in gray. Since we impose a limit of 35 au to avoid fitting potential narrow emission within the line, the center of the line is not well fitted. Lower panel: for each line, the red and blue thick lines show the best fit to the median spectrum for the red and the blue sides, while the shaded areas show the range spanned by the fits to the individual spectra. Since the noise in the single spectra is higher than in the median spectrum, there are no significant results at extreme velocities for some of the individual datasets. The vertical dotted line indicates the place where direct irradiation results in a temperature of 2100 K, which is independent on the model assumptions and stellar mass. All models are calculated for an inclination of 30 degrees and a stellar mass of 16 M⊙. A different inclination i would change the derived radii by sin2i/sin2(30), and a different stellar mass M* would change them by a factor of M*/16 M⊙, but the relative brightnesses would remain the same.

Fig. 14 Comparison of the variability observed in several lines in outburst and quiescence. The average flux is plotted in blue, while the standard deviation is shown in black with gray shadowing. The zero velocity (dotted vertical line) is given with respect to the stellar radial velocity. Note: the Fe I 6400 Å line in quiescence may contain some contribution of a nearby forbidden [Fe I] line.

Fig. 15 Line profiles of the strong Ti II as observed in quiescence. The data for the various dates are all plotted together, given that there is no significant variability. The vertical line indicates the 27 km s−1 radial velocity, and the dotted horizontal line indicates the continuum level.

4 Outbursts in high-mass versus low-mass stars

One of the most surprising results of our analysis is that the spectral behavior of Z CMa NW in outburst and quiescence is qualitatively not too different from what has been observed in other EXors (Herbig 2007, 2008; Sicilia-Aguilar et al. 2012, 2017; Holoien et al. 2014; Audard et al. 2014). From historical photometric evidence (see Appendix A), the outbursts of Z CMa NW are not different from the outbursts of lower-mass EXors, repeating in time every few years and showing various patterns in time duration and intensity (Herbig et al. 2001; Herbig 2007; Audard et al. 2014). In both EX Lupi and Z CMa NW, the accretion-related wind becomes stronger during outburst, and the line wings of essentially all lines move toward greater (positive and negative) velocities. Emission originating from the inner disk is also detected in both objects, although in EX Lupi the emitting material within the inner disk is highly non-axisymmetric (e.g., a single accretion column or disk spiral) and very close-in, which results in dramatic day-to-day velocity variations in the broad components of the lines (Sicilia-Aguilar et al. 2012). Non-axisymmetric winds have also been observed in the M5 star ASASSN-13db during outburst (Sicilia-Aguilar et al. 2017). Evidence of changes in the inner disk content after bursts has been documented in EX Lupi (Banzatti &Pontoppidan 2015), which is consistent with the decrease in optical thickness in the line-emitting region observed for Z CMa NW after outburst. The larger luminosity of this outburst agrees with the idea of a continuum of bursting stars, not only in the magnitude of the burst (Contreras Peña et al. 2017), but also in stellar mass.

The emission lines observed in Z CMa NW, EX Lupi, and ASASSN-13db have statistically indistinguishable distribution of transition probabilities. The only significant difference arises between the excitation potentials of the lines observed in Z CMa during outburst and ASASSN-13db, thus rejecting at the 1% level the hypothesis that both samples are drawn from the same distribution,according to a double-sided Kolmogorov-Smirnov test. This difference in excitation potentials likely results from the large difference in stellar temperature (B-type vs. M5). The similarities point to a rather uniform range of physical conditions for the materialin the inner disk and accretion structures. However, because of the differences in the masses and effective temperatures of the three stars, the regions in which this emission originates ranges from a few stellar radii for ASASSN-13db (Sicilia-Aguilar et al. 2017), to 0.1–0.2 au for EX Lupi (Sicilia-Aguilar et al. 2012), and to 0.5–3 × (M_*/16 M_⊙) au for Z CMa NW.

The accretion scenario for Z CMa NW can be also explored through its emission and absorption lines. Although contributions from different origins (e.g., nonspherical winds) may affect the extremely redshifted Hα wings, thereare some small redshifted absorption asymmetries in the Ti II and Fe I lines (in the range of 50–150 km s⁻¹) and in the high-velocity wings of accretion-related lines. These asymmetries are mainly observed during outburst in the higher Balmer lines and in H IR lines (see Fig. 2, small absorption components appear at approximately +200 and +400 km s⁻¹ in Hβ and ~250 km s⁻¹ in Hγ). These could indicate a high-velocity infall scenario (e.g., Koenigl 1991; Edwards et al. 1994; Muzerolle et al. 1998) rather than a classical boundary layer, where the velocities in the boundary would be much reduced (e.g., Popham et al. 1993). Nonetheless, these velocities suggest a stellar mass much lower than 16 M_⊙, even if the magnetosphere was very small. Using as the maximum infall velocity that observed in the Hα line wings (1000 km s⁻¹) this would result in a 2.7 R_* or 5.2 R_⊙, which is on the small side for magnetospheric accretion in intermediate-mass stars (e.g., Mendigutía et al. 2017), albeit not unexpected in a fast-rotating star with a weak magnetic field. A variation in stellar mass by a factor of two would bring the infall radius to the classical infall radius for a low-mass stars (Gullbring et al. 1998).

Nevertheless, the corresponding field to ensure a magnetospheric cavity extending as far as the corotation radius for R_cor = 2.5 R_* would be on the upper limit for HAeBe stars. The minimum magnetic field B_min required to support magnetospheric accretion can be estimated as $$B_{\text{min}}=1.1{(\frac{M_{*}}{M_{\odot }})}^{2/3}{(\frac{\dot{M}}{{10}^{-7}{M}_{\odot }{\text{yr}}^{-1}})}^{23/40}{(\frac{R_{*}}{R_{\odot }})}^{-3}{(\frac{P}{1d})}^{29/24}\text{kG}$$(7)

(see Collier Cameron & Campbell 1993; Johns-Krull et al. 1999; Mendigutía 2020) as a function of the rotational period of the star, P, which is a priori not known. A corotation radius of 2 × R_* requires a stellar rotation period of 0.22 d (for a 16 M_⊙ star), down to 0.44 d (for a 4 M_⊙ star), but the minimum magnetic field needed to support the accretion rate observed in outburst via magnetospheric accretion would be at least 8 kG and up to 16kG for a 4 M_⊙ star, since a lower stellar mass also means a higher accretion rate for the same accretion luminosity. This is much larger than the existing estimates for the Z CMa pair (1.2 kG; Szeifert et al. 2010; Hubrig et al. 2015) and the typical fields of HAeBe stars. On the other hand, assuming that in quiescence the accretion rate drops to the typical HAeBe rate of 10⁻⁷M_⊙ yr⁻¹, a magnetic field of 0.2 kG would be able to support magnetospheric accretion for a 16 M_⊙ star, and a lower-mass star would require a weaker field.

A more detailed estimate of the equatorial magnetic field required could be derived following (Bessolaz et al. 2008) $$\begin{array}{ll}{B}_{*}=\hfill & 140\times 2^{-7/4}{m}_{\text{s}}^{-1/2}{\left(\frac{R_{\text{t}}}{R_{*}}\right)}^{7/4}{\left(\frac{\dot{M}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1./2.}{\left(\frac{M_{*}}{0.8M_{\odot }}\right)}^{1/4}\hfill \\ \hfill & \times {\left(\frac{R_{*}}{2R_{\odot }}\right)}^{-5/4}G,\hfill \end{array}$$(8)

where R_t is the truncation radius, m_s is the Mach number (~0.45; Bessolaz et al. 2008), and the rest of symbols retain their previous meaning. This approach results in an even higher magnetic field during outburst (46–65 kG) for a star with a mass between 4–16 M_⊙, which is completely unfeasible. The result for the quiescence accretion rate of 1e-7 M_⊙ yr⁻¹ would range from 1–1.5 kG for stars with masses between 4–16 M_⊙, which although feasible for the more magnetic HAeBe stars, remains questionable until such a high magnetic field has been confirmed. Therefore, although magnetospheric accretion is highly unlikely in outburst, which would be a significant difference compared to lower-mass EXors, it may be possible in quiescence.

5 Summary and conclusions

We present optical spectroscopy during the outburst and post-outburst phase of Z CMa NW, using the wealth of emission lines observed to explore the causes of the outburst and the properties of the accreting star. We confirm that the outburst of Z CMa NW was related to an increased accretion episode reaching Ṁ ~ 10⁻⁴ M_⊙ yr⁻¹. This accretion rate is significantly higher than the current accretion rate of the FUor companion Z CMa SE, for which photometric and spectroscopic evidence suggest a return to quiescence by year 2000. The historical light curve (including AAVSO data since the 1930s) and past spectroscopic activity (e.g., Hessman et al. 1991) suggest that the outbursts of Z CMa NW are repetitive, as in low-mass EXors. We use several techniques – including velocity decomposition to track the spatial location, Saha’s equation, and line ratios of lines from the same upper level – to estimate the density and temperatureof the various regions, which include a non-axisymmetric, multicomponent wind and extended disk emission.

The optical spectra strongly resemble EXor outburst spectra for lower-mass stars, suggesting a continuum of accreting behaviors for disked stars with all masses. The spectra of Z CMa NW in outburst reveal lines with multiple origins. A double-peaked line emission originates in a hot disk (detected in neutral atomic lines, especially Fe I). Complex wind absorption with several velocity components is most evident in H, He I, and in the higher-energy metallic ionized linessuch as Fe II, Si II, and Ti II. There is also some evidence of redshifted asymmetries that could be related to infall or accretion structures, which are detected preferentially in H and single-ionized elements such as Fe II, and Ca II in outburst, and in the asymmetry of neutral lines in quiescence. This is not dissimilar to what is seenin lower-mass stars, although the velocities also suggest that the stellar mass is likely lower than 16 M_⊙. Age and luminosity considerations suggest a value around 6–8 M_⊙.

The magnetosphere – if present – is expected to be very small, and it is highly unlikely that the magnetic field is strong enough to support magnetospheric accretion during outburst, although it may be possible in quiescence if the star belongs to the most magnetic HAeBe stars. This is significantly different from lower-mass EXors, and challenges the role of amagnetosphere in the triggering of some outbursts (e.g., as suggested for EX Lupi; Sicilia-Aguilar et al. 2015). We do not observe a boundary layer in the classical sense either, but rather extended accretion structures. We also observe a more distant Keplerian, irradiated disk. The energetics of the optical emission lines may impose a limitation to the temperature we can reach, and UV lines may be needed to trace such a hotter, closer structure.

The disk emission observed in Fe I and Ca I lines is best fitted with an irradiated, flared disk originating in the inner ~0.5–3 × (M_*/16 M_⊙) au regions, which have temperatures in the 2000–7000 K range and densities that decrease toward lower velocities (larger radii). Weak Fe II lines are also consistent with disk emission, although they trace higher temperature regions. For the disk to be dominated by irradiation, it is necessary that the bulk of energy released in the accretion burst is deposited at a small radius (e.g., in the innermost disk, accretion structures, or on stellar hot spots). There is no evidence of the temperature inversion associated with viscous heating in strongly accreting systems, but this could happen if the part of the irradiated disk that we observe in the Fe I lines is too distant from the hot, dense, accretion-powered inner disk. The disk-like features of higher energy and ionized lines from the hotter, innermost regions in the disk may be masked by the strong emission and absorption profiles produced by the hot stellar wind. We observe a decrease in column density-velocity gradient per velocity bin, and for the redshifted components of the lines, decreases by at least one order of magnitude between outburst and quiescence, a sign of a lower mass transport through the disk.

The blueshifted absorption features are consistent with an accretion-powered wind. The fastest and slowest wind components disappear in quiescence, when the wind density decreases. The temperature in the fast wind is higher than in the slow wind, betraying an origin very close to the star, while the slow component could be a disk wind. The changes in velocity and depth of the absorption are consistent with a clumpy, non-axisymmetric wind with spatially variable opacity, modulated by rotation. This is similar to the winding winds observed in B-type stars (Prinja & Howarth 1988), but the behavior disappears in quiescence and thus is linked to the outburst rather than to the spectral type. Our small number of observations does not allow us to determine whether the velocity changes display any periodicity, but the presence of a non-axisymmetric wind is an additional sign that accretion is non-axisymmetric as well.

This study, together with our previous work on low-mass stars (Sicilia-Aguilar et al. 2012, 2015, 2017), demonstrates that optical spectroscopy can be used to trace the structure, radial distance, and physical conditions of accretion in stars with a large range of mass. Combined with historical data, it shows that a similar outbursting behavior can be observed through a very large range of stellar masses and spectral types (from M5 to B-type), pointing to common outburst mechanisms and similarities in accretion and inner disk transport within this whole range, despite potential differences in the structure and size of the magnetosphere. The limited energy ranges covered by the lines that can be observed in optical spectra also point out the need for high-resolution UV and IR spectroscopy data to study the full energetic range of this kind of objects.

Acknowledgments

We thank Lynne Hillenbrand, Andrew Collier Cameron, Lee Hartmann, Peter Woitke, and Ignacio Mendigutía for inspiring discussion, as well as the anonymous referee for valuable comments. A.S.A. is partly supported by the STFC grant number ST/S000399/1 (‘‘The Planet-Disk Connection: Accretion, Disk Structure, and Planet Formation”) and by the 2017 Visiting Teacher Campaign at the University of Grenoble. We thank the former OHP Director, Michel Boer, for granting discretionary time to perform these observations. We thank Anne Eggenberger for performing the second set of spectroscopic observations at OHP during quiescence. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742095; SPIDI: Star-Planets-Inner Disk-Interactions, http://spidi-eu.org, and grant agreement No 740651 NewWorlds). We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research, and in particular, observers from the Royal Astronomical Society of New Zealand (RASNZ).

Appendix A Variability: 90 yr of photometric data

The AAVSO contains data on Z CMa since the 1930s. These data provide a great opportunity to examine the behavior of the Z CMa pairover nearly 90 yr. The AAVSO confirms that the light curve is anomalous compared to other FUors (Hartmann et al. 1989; Hessman et al. 1991), since the documented magnitude rise spans only 2 mags and seem to have progressed very slowly for the first 30 yr. The light curves are shown in Fig. A.1, suggesting that the FUor outburst probably ended by 1995. Since then, the outbursts of Z CMa NW have been recurrent and increasing in brightness, although both the AAVSO data and early mentions of small-scale variability over the FUor curve (Hessman et al. 1991) suggest that the outbursts may have been occurring at least since the early 1980s, masked in part by the FUor outburst of the companion. Thus the P Cygni profiles observed to dominate the spectrum during bumps (Hessman et al. 1991) may have been caused by the HBe star, and not by its FUor companion. In any case, the bursts appear to have increased in intensity since the 2000s, and deserve some future follow-up.

A generalized Lomb-Scargle periodogram (GLSP; Scargle 1982; Horne & Baliunas 1986; Zechmeister, & Kürster 2009) reveals a low-significance quasi-periodof 249 ± 20 d (Fig. A.2) when the ASAS and AAVSO V data (covering 14 yr) are taken into account. The modulation appears to change inphase, which suggests a quasi-periodic behavior in a component that changes at about the same rate as the periodicity is observed. Using a red noise model (as done in Sicilia-Aguilar et al. 2020), the false-alarm-probability of the modulation is ~1%. Although this potential period needs to be confirmed in the coming years, it could be related to repeated disturbances in the innermost disk causing the bursts.

As seen in Fig. 1, the two first spectra were taken as the star emerged from an occultation-like event, in which the rising magnitude experienced a drop of more than 0.5 magnitudes and then recovered; the whole event lasted about 30 days. Unfortunately, the spectra were taken when the star had nearly recovered, so any occultation effects may be mild. These two first datasets also seem to show distinct wind features, such as a less complex profile more similar to the quiescence profile. This could be caused by either a larger contribution of the FUor spectrum if Z CMa NW is partly occulted or suggest that part of the wind-originating location could have been occulted by the material that caused the observed dip. The high temperatures inferred for the fast wind are consistent with an origin closer to the star, which would be also easier to be hidden by an occultation event than, for instance, a more extended disk wind. The effect is mostly observed in Fe II lines (5018, 4923, 4233 Å), which are also some of the most variable lines among those with P Cygni profiles. Nevertheless, since the first spectrum has stronger wind absorption than the second spectrum, the wind is by itself very variable andthese observations were taken about a month before the following observations, what we observe in the wind may be just part of the usual wind variability unrelated to the occultation event. Figure 1 reveals that similar dips during previous outbursts of Z CMa NW have occurred. They are more evident during the last ten years, after the FUor outburst had dimmed and during the time when the activity of Z CMa NW seems to have increased, and could be related to variable extinction suggested by Blondel & Djie (2006). Future follow-up may help to understand the origin of the dips and to locate the occulting material and the regions it affects.

Fig. A.1 Light curves from AAVSO for Z CMa since 1930s (upper) and from 1999 to 2018 (lower). Only verified observations are included. The visual and V filters are plotted together, since they are similar for the current purpose of examining the timescales of the outbursts. As a comparison, the ASAS V filter is also plotted (red), which is in excellent agreement with the amateur V data.

The line identification was done using the NIST database, except for lines denoted with WS, for which the line identification and atomic parameters have been taken from Wahlgren & Shore (S. Shore, priv. comm.) and their observationsof the symbiotic star V407 Cyg outburst. We compare the observed lines with the lists of van den Ancker et al. (2004) from observations of Z CMa between the late 1990s and the early 2000s (including weaker outbursts) and the lines observed in two outbursting stars, EX Lupi (Sicilia-Aguilar et al. 2012, 2015) and ASASSN-13db (Sicilia-Aguilar et al. 2017), for which a similar analysis is published. The observed profiles were also used to classify a line, or at least to identify the dominant component in case of blends or multiple options, when more than one transition was possible. As mentioned in the text, strong Fe II lines (and also Cr II and Si II lines) have P Cygni profiles. Weak Fe II lines with low transition probabilities are similar to Fe I lines. A few forbidden lines are also observed; [O I] 6300/6363 Å is the strongest, followed by [S II]. Some redshifted outflow emission could be present, butblends make this possibility uncertain at best. Some [Fe II] lines are detected, but they tend to be very weak and, in most cases, they appear blended with or very close to other very strong lines so that their profiles are hard to compare and the detections remain tentative. Some of the [Fe II] lines may also be observed in quiescence.

Fig. A.2 Data from ASAS and AAVSO V wrapped according to a 249.1 d period. The data are colored depending on the epoch to enable visualization.

Appendix B Complete list of emission lines

Appendix C Methodology

This appendix contains the details of the different methods used to analyze the emission and absorption lines that are applied in Sect. 3.

C.1 Deriving physical properties from the spectral lines

To constrain the physical conditions of lines originating from the same location, we compared the strength of the ionized and neutral lines to put constraints on the density and temperature using Saha’s equation, $$\frac{n_{j+1}{n}_{\text{e}}}{n_j}={(\frac{2\pi m_{\text{e}}kT}{h^2})}^{3/2}\frac{2U_{j+1}(T)}{U_j(T)}{e}^{-{\chi }_{\text{I}}/kT},$$(C.1)

where n_j+1 and n_j are the densities of the ionized and neutral species, respectively; χ_I is the ionization potential; n_e is the electron density; m_e is the electron mass; k is the Boltzmann constant; h is the Planck constant; T is the temperature; and U_x (T) is the partition function for the x-ionized ion at a temperature T (Mihalas 1978). The main limitations of this approach are non-LTE effects and optical depth and self-absorption (as observed in EX Lupi, Sicilia-Aguilar et al. 2012). We expect at least some of the lines to be optically thick, which is a problem when using line ratios to derive the physical properties. Nevertheless, the non-detection (or detection) of low-energy lines is a very powerful way to set a limit to the temperature for reasonable densities, or vice versa. For this exercise, we used the atomic parameters and partition functions for different temperatures from the NIST spectral line database¹¹ (Ralchenko et al. 2010; Kramida et al. 2018) and made sure that the line profiles of the compared transitions are very similar in velocity and line profile, to minimize the risk of including emission produced by different structures.

Fig. C.1 Analysis of the line ratio for Fe I transitions with upper common levels (Fig. 9, continued). Upper left panels: line emission in outburst; upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. The gray area shows the region for which the line emission is >5σ (in the velocity panels), and the 1σ error in the temperature and density planes.

A further constrain is obtained from the ratio of the transitions of the same ion. As a first approach, if we ignore radiative transfer effects and assume that the lines are optically thin, the levels are populated according to Boltzmann statistics, so that $$\frac{n_i}{n_0}=\frac{g_i}{g_0}{e}^{-Ei/kT},$$(C.2)

where n_i and n₀ are the densities in the levels i and the ground state, respectively; g_i and g₀ are their statistical weights; and E_i is the energy of the level with respect to the ground state.

A second calculation can be performed considering the ratio of two lines that share a common upper level. In such a case, the way the upper level is populated does not matter, and the line ratio depends only on the optical thickness of each line (Beristain et al. 1998). If we estimate the ratio of the weakest to the strongest line, the optically thin limit (for which both lines are optically thin) depends on the line properties alone. While the strongest line is saturated, the ratio increases until it reaches a limit once the weakest line is also saturated. In particular, for isotropic high-velocity gradients, we can use the Sobolev approximation to treat the escape probability (Beristain et al. 1998) so that the line ratio can be written as $$r_{w/s}=\frac{A_{ki}^w{\lambda }_s}{A_{kj}^s{\lambda }_w}\frac{1-e^{-{\tau }_w}}{{\tau }_w}\frac{{\tau }_s}{1-e^{-{\tau }_s}},$$(C.3)

where A_ul is the Einstein coefficient for the spontaneous transition from levels u to l, λ is the wavelength, and τ is the opacity, and the indices w and s mark the weak and the strong line, respectively. The opacity in the Sobolev approximation for the transition between levels k (up) and i (low) can be written as $$\tau =\frac{hc}{4\pi }\frac{n_i{B}_{ik}-n_j{B}_{ki}}{\text{d}v/\text{d}l}\approx \frac{hc}{4\pi }{B}_{ik}{g}_i\frac{e^{-E_i/kT}}{U(T)}\frac{n\text{d}l}{\text{d}v},$$(C.4)

where B_ik and B_ki are the Einstein coefficients for stimulated absorption and emission and the approximation can be obtained by neglecting the population in the upper level and assuming thermal equilibrium (Mihalas 1978; Beristain et al. 1998). The velocity gradient is written dv∕dl and E_i is the energy of the lower level. Rearranging Eqs. (C.3) and (C.4), we can write the line ratio in terms of the temperature and modified column density or column density-velocity gradient $\mathcal{N}=n\text{d}l/\text{d}v$.

Using nearby lines (with wavelengths within a ≤500 Å range) is preferred to avoid effects of underlying variable continuum. We note that the main uncertainty related to the atomic parameters comes from the B_ij parameters, thus the main uncertainties lie in the derived line opacities. This exercise was done for several velocity ranges along the line profile (see Sect. 3.4). Some lines are shown in the main text, with the remaining lines shown in Figs. C.1 and C.2.

Fig. C.2 Analysis of the line ratio for the FeII and TiII lines with PCygni profiles from transitions with upper common levels (Fig. 11, continued). Upper left panels: line emission in outburst; upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. In the velocity panels, the red area, cyan,and the gray area show the regions with significant redshifted emission, and slow and fast wind absorption, respectively. The shaded colored areas in the temperature and density planes show the regions that are compatible with the line ratios observed in the emission and the two wind component parts with their 1σ error.

C.2 Disk analysis by brightness decomposition

To explore the location of the material producing the double-peaked emission lines, we used a procedure similar to Acke & van den Ancker (2006) to convert line intensity in each velocity bin into a relative brightness per unit area at each velocity. Assuming Keplerian rotation, we can convert the velocities into radial distances. One main limitation comes from gaseous disks being partly supported by their own gas pressure, so the rotation of hot atomic gas is sub-Keplerian. The thermal deviation from Keplerian rotation depends on the temperature, and this quantity ranges from ~11 km s⁻¹ for a temperature of 10⁴ K, to ~3.5 km s⁻¹ for a temperature of 10³ K. This is a small fraction of the total velocity span observed in the lines, justifying the approximation. We fit the ESPaDOnS data only, since the OHP data are too noisy for the fit.

For the brightness decomposition, we assume a disk inclination of 30 degrees and the outburst stellar values (R_* = 1.93 R_⊙, M_* = 16 M_⊙, L_* = 8 × 2400 L_⊙). A different inclination or the stellar mass can be easily rescaled in the brightness decomposition plots, as long as the luminosity is kept constant¹². With these values we can estimate an effective temperature T_eff, assuming that the stellar radius does not change during outburst with respect to the quiescent value, that reproduces the observed luminosity. We can also estimate the temperature at various distances from the star, which depends on the heating mechanism and disk structure. We assume that the temperature of the innermost disk rim (the part producing the highest observed velocities) is set by direct irradiation, and for the rest of the disk it follows a power law with radius. Although the temperature estimate would be smaller if the emitting zone is larger (e.g., if the innermost disk is the emitting region), at larger distances the only quantity that matters is the total luminosity (or $R_{\text{emit}}^2{T}_{\text{emit}}^4$) as long as the emitting region is located at a smaller radius than the detected disk emission¹³.

The temperature power law is poorly known, since we can have various sources of heating (e.g., irradiation, viscous dissipation) and gaseous disks are flared and have a vertical temperature structure. We tested several temperatureversus radius power laws. These include the Chiang & Goldreich (1997) structure for an irradiated flat disk, for which the temperature in the disk is taken to be a power law of radius with exponent −3/4, which is also the same power law than if we assume a self-luminous disk powered by viscous dissipation (e.g., Hartmann 1998); and the flared irradiated disk model, in which the temperature varies as a power law with an exponent −3/7 versus radius (D’Alessio et al. 1999). The result of using a steeper versus a shallower power law is that the strongest relative brightness moves to larger radii when the temperature decreases faster (i.e., if the temperature is lower, a larger area is needed to reproduce the observed flux), preserving the general shape of the profile and the blue/red (a)symmetry. Figure C.3 shows some examples of the effect in the derived brightness profiles depending on the power law assumed. In general, the residuals of the flareddisk models are smaller than those of the flat disk fits, thus we conclude that an irradiated flared disk model can fit better the velocity profiles observed, and this is what we use through Sect. 3.5.

Fig. C.3 Examples of how the temperature power law affects the radius derived from the brightness profile. For each line, the left models show the results from a flat disk (power-law exponent −3/4; Chiang & Goldreich 1997), while the right models show the results assuming a flared disk (power-law exponent −3/7; D’Alessio et al. 1999). The difference is larger for lines with higher transition probabilities (Aki = 9.3e6 s−1 for FeI 6400 Å,Aki = 5.6e4 s−1 for FeI 6462 Å) that are likely to be formed higher over the disk midplane. Colors and symbols as in Fig. 13.

Finding Fe I emission in LTE is unlikely at temperatures below a thousand degrees. For a flat, irradiated disk and the stellar luminosity in outburst, T = 1000 K is reachedat about 20 au. Even though other sources of heating (e.g., thermal dissipation due to accretion) are likely to operate in this strongly accreting system, we impose a cutoff at 35 au. Beyond this cutoff we deem that the emission becomes unreliable, especially for the strongest lines that often have an emission core at low velocities that cannot be separated from the disk emission¹⁴. We also estimated the radial velocity of the line using the 80% level of the double peak of the line to estimate the line center, which typically produces radial velocity values between 20–35 km s⁻¹, in agreement with the narrow quiescence lines. A highly asymmetric line and/or eccentric or asymmetric disk could lead to inaccurate values of the radial velocity.

Following Acke & van den Ancker (2006), we decomposed the normalized line profile starting from the highest velocity at which the flux is ≥3σ. We assume that the emission at each velocity originates in a narrow ring located at the position with the same Keplerian orbital velocity. We generated simple disk model emission, assuming that each element of the velocity ring emits a line where the width is due to thermal broadening. The emission of the ring is proportional to the blackbody emission at the relevant temperature and to the area of the emitting ring. We constructed the emission of the ring by summing over all angles, breaking the ring into a 70 point grid in radius and azimuthal angle. The ratio of the observed to theoretical flux at each velocity point gives us the brightness associated with that ring. The resulting ring is subtracted from the total observed line, and the procedure is repeated for the next highest velocity until the line has been completely fitted or until we reach the 35 au cutoff.

The comparison between the blue and the red sides can provide information on the axisymmetry of the disk. In addition, this simple model allows to check whether there is emission at locations more distant than those expected to have a temperature of up to 2100 K (similar to the lower limit we obtained from Saha’s equation) for a directly irradiated disk. This would suggest that the disk is not (or not only) directly irradiated, but also heated by other mechanisms such as accretion, but in general the emission is very low at larger radii. The analysis is only valid if the line is not contaminated by nearby features and the continuum is relatively smooth. This results in only a subset of the lines with disk-like profiles being appropriate for the brightness decomposition, which is discussed in Sect. 3.5. A few examples of the brightness decomposition are shown in the main text; the rest of the lines can be found in Fig. C.4.

Fig. C.4 Brightness profile decomposition for the disk-like lines not shown in the main text (Fig. 13, continued). For each line, upper panel: normalized, median line emission (blue), individual data points (colored dots), and model fit. The black solid line represents the best fit to the median data, while the dashed lines indicate the fits to the individual datasets that are used to compute the errors. The zero velocity (measured with respect to the symmetry axis of the line) is shown as a dotted line. The individual ring fits are shown by thin red lines, and the residuals are shown in gray. Since we impose a limit of 35 au to avoid fitting potential narrow emission within the line, the center of the line is not well fitted. Lower panel for each line, the red and blue thick lines show the best fit to the median spectrum for the red and blue sides, while the shaded areas show the range spanned by the fits to the individual spectra. Because the noise in the single spectra is higher than in the median spectrum, there are no significant results at extreme velocities for some of the individual datasets. The vertical dotted line indicates the place where direct irradiation results in a temperature of 2100 K, which is independent of the model assumptions and stellar mass. All models are calculated for an inclination of 30 degrees and a stellar mass of 16 M⊙. A different inclination i would change the derived radii by sin2i/sin2(30), and a different stellar mass M* would change the derived radii by a factor of M*/16 M⊙, but the relative brightnesses would remain the same.

References

All Tables

All Figures

	Fig. 1 Full light curve (left) from ASAS and zoom around the times of the spectroscopic observations (right). The times at which thespectra were taken are indicated as blue for outburst spectra and orange for quiescence spectra. A complete light curve is given in Fig. A.1.
In the text

Fig. 2 Examples of lines observed for Z CMa. The lines are normalized to the local continuum. In this and the following figures, the quiescence data are plotted in various shades of orange, while the outburst spectra are shown in blue. The dashed line shows the radial velocity of the system (27 km s−1) and the dotted line indicates the continuum level. First row: some H Balmer lines and an IR H line. The quiescence spectra of Hα and Hβ are scaled by × 0.5 for better display. Second row: some of the Fe I lines with disk-like profiles and a K I resonance line. Third and fourth rows: some lines with various types of P Cygni profiles. The quiescence spectra do not cover the range beyond 6940 Å.

In the text

Fig. 3 Upper: forbidden [O I] emission. Lower: forbidden [S II] emission. The dotted lines indicate the continuum level, and the dashed lines represent the radial velocity of 27 km s−1. The data from the outburst are plotted in various shades of blue, while the quiescence data are shown in orange. The extra lines observed in the 6363 Å spectrum correspond to unrelated Fe I and Fe II lines, and the strong absorption is due to Fe II 6369 Å. The narrow absorption lines are telluric.

In the text

Fig. 4 Zoom around the Li I 6708 Å line. The dotted black line indicates the continuum level, the dashed black lines indicate the Li I transitions (6707.76 and 6707.91 Å), the vertical red dotted lines indicate the Fe I 6707.43 Å and Fe II 6708.88 Å lines, and the green dotted line indicates the V I line at 6708.07 Å; all of these are shifted to match a radial velocity of 27 km s−1. The emission feature in quiescence is likely due to the high-velocity component of [S II].

In the text

	Fig. 5 Extinction in the Gaia G band and Gaia BP-RP around a 0.2 degree radius field toward Z CMa. The rise of extinction is fully consistent with the location of the Z CMa cloud at 1 kpc.
In the text

	Fig. 6 Examples of the time variability observed for the Ti II 4468 Å line (left, blue), the Fe II 4923 Å line (middle, green), and the Fe II 5316 Å line (right, red). In each case, the 8 subfigures display the individual spectra as they were obtained in time. The dotted lines indicate the continuum level and the radial velocity of the source.
In the text

Fig. 7 Day-to-day velocity variation of the wind components in 2009 January, classified according to the central velocities as very fast (top, observed clearly only for Fe II), fast (middle), and slow (lower). Results for the Fe II 5316 (red), Fe II 4923 (green), and Ti II 4468 Å (blue) lines, which have good S/N and no signs of contaminations, are shown. The velocities are given as measured, without correcting for the radial velocity of the system. The correlations between the components observed for both intermediate velocity and slow wind lines, despite some shifts in velocity, are shown. The fast wind component appears to split into several components for some dates for the Fe II 4923 Å.

In the text

Fig. 8 Result of Saha’s equation analysis for the lines with disk-like profiles and low transition probabilities for the observationstaken during the outburst phase (upper panels) and for the same lines as observed during quiescence (lower row). The panels show the data corresponding to individual days in chronological order; last panel: result of fitting the average spectrum. The color scheme reveals the relative variations in the line intensity ratio for each combination of temperature and electron density (such that 0.1 would correspond to 10% average deviations). The white contour indicates the regions for which the relative variations of the line ratio are up to 3 times the minimum value observed, which is considered as our best-fit region; a white star denotes the best fit.

In the text

Fig. 9 Analysis of the line ratio for Fe I transitions with common upper levels (continued in Appendix C). Upper left panels: line emission in outburst and upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. The gray area shows the region for which the line emission is >5σ (in the velocity panels) and the 1σ error in the temperature and density planes. Quiescence lines are often very weak and thus uncertain.

In the text

	Fig. 10 Result of the line ratio analysis for neutral lines with disk profiles in outburst (left) and narrow profiles in quiescence (right).
In the text

Fig. 11 Analysis of the line ratio for the FeII and TiII lines with PCygni profiles from transitions with common upper levels (continued in Appendix C). Upper left panels: line emission in outburst and upper right panel: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. In the velocity panels, the red, cyan, and gray areas show the regions with significant redshifted emission and slow and fast wind absorption, respectively. The shaded colored areas in the temperature and density planes show the regions that are compatible with the line ratios observed in the emission and the two wind component parts with their 1σ error. See text for details.

In the text

	Fig. 12 Result of the line ratio analysis ion pairs of lines with a common upper level for the emission and wind absorption parts (denoted as “wind” in the legend). Some of the Ti II lines have sharp cuts in the n dv∕d l vs. T plane due tounphysical values that arise in regions of low S/N (see Appendix C). The results have to be regarded with care when approaching such regions.
In the text

Fig. 13 Brightness profile decomposition for the disk-like lines (continued in Appendix C). For each line, upper panel: normalized, median line emission (blue), the individual data points (colored dots), and the model fit. The black solid line is the best fit to the median data, while the dashed lines are the fits to the individual datasets that are used to compute the errors. The zero velocity (measured with respect to the symmetry axis of the line) is shown as a dotted line. The individual ring fits are shown by thin red lines, and the residuals are shown in gray. Since we impose a limit of 35 au to avoid fitting potential narrow emission within the line, the center of the line is not well fitted. Lower panel: for each line, the red and blue thick lines show the best fit to the median spectrum for the red and the blue sides, while the shaded areas show the range spanned by the fits to the individual spectra. Since the noise in the single spectra is higher than in the median spectrum, there are no significant results at extreme velocities for some of the individual datasets. The vertical dotted line indicates the place where direct irradiation results in a temperature of 2100 K, which is independent on the model assumptions and stellar mass. All models are calculated for an inclination of 30 degrees and a stellar mass of 16 M⊙. A different inclination i would change the derived radii by sin2i/sin2(30), and a different stellar mass M* would change them by a factor of M*/16 M⊙, but the relative brightnesses would remain the same.

In the text

	Fig. 14 Comparison of the variability observed in several lines in outburst and quiescence. The average flux is plotted in blue, while the standard deviation is shown in black with gray shadowing. The zero velocity (dotted vertical line) is given with respect to the stellar radial velocity. Note: the Fe I 6400 Å line in quiescence may contain some contribution of a nearby forbidden [Fe I] line.
In the text

	Fig. 15 Line profiles of the strong Ti II as observed in quiescence. The data for the various dates are all plotted together, given that there is no significant variability. The vertical line indicates the 27 km s−1 radial velocity, and the dotted horizontal line indicates the continuum level.
In the text

	Fig. A.1 Light curves from AAVSO for Z CMa since 1930s (upper) and from 1999 to 2018 (lower). Only verified observations are included. The visual and V filters are plotted together, since they are similar for the current purpose of examining the timescales of the outbursts. As a comparison, the ASAS V filter is also plotted (red), which is in excellent agreement with the amateur V data.
In the text

	Fig. A.2 Data from ASAS and AAVSO V wrapped according to a 249.1 d period. The data are colored depending on the epoch to enable visualization.
In the text

Fig. C.1 Analysis of the line ratio for Fe I transitions with upper common levels (Fig. 9, continued). Upper left panels: line emission in outburst; upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. The gray area shows the region for which the line emission is >5σ (in the velocity panels), and the 1σ error in the temperature and density planes.

In the text

Fig. C.2 Analysis of the line ratio for the FeII and TiII lines with PCygni profiles from transitions with upper common levels (Fig. 11, continued). Upper left panels: line emission in outburst; upper right panels: quiescence data. The ratios observed are compared to the theoretical calculations for a range of temperatures and column density-velocity gradients in the lower panels. In the velocity panels, the red area, cyan,and the gray area show the regions with significant redshifted emission, and slow and fast wind absorption, respectively. The shaded colored areas in the temperature and density planes show the regions that are compatible with the line ratios observed in the emission and the two wind component parts with their 1σ error.

In the text

Fig. C.3 Examples of how the temperature power law affects the radius derived from the brightness profile. For each line, the left models show the results from a flat disk (power-law exponent −3/4; Chiang & Goldreich 1997), while the right models show the results assuming a flared disk (power-law exponent −3/7; D’Alessio et al. 1999). The difference is larger for lines with higher transition probabilities (Aki = 9.3e6 s−1 for FeI 6400 Å,Aki = 5.6e4 s−1 for FeI 6462 Å) that are likely to be formed higher over the disk midplane. Colors and symbols as in Fig. 13.

In the text

Fig. C.4 Brightness profile decomposition for the disk-like lines not shown in the main text (Fig. 13, continued). For each line, upper panel: normalized, median line emission (blue), individual data points (colored dots), and model fit. The black solid line represents the best fit to the median data, while the dashed lines indicate the fits to the individual datasets that are used to compute the errors. The zero velocity (measured with respect to the symmetry axis of the line) is shown as a dotted line. The individual ring fits are shown by thin red lines, and the residuals are shown in gray. Since we impose a limit of 35 au to avoid fitting potential narrow emission within the line, the center of the line is not well fitted. Lower panel for each line, the red and blue thick lines show the best fit to the median spectrum for the red and blue sides, while the shaded areas show the range spanned by the fits to the individual spectra. Because the noise in the single spectra is higher than in the median spectrum, there are no significant results at extreme velocities for some of the individual datasets. The vertical dotted line indicates the place where direct irradiation results in a temperature of 2100 K, which is independent of the model assumptions and stellar mass. All models are calculated for an inclination of 30 degrees and a stellar mass of 16 M⊙. A different inclination i would change the derived radii by sin2i/sin2(30), and a different stellar mass M* would change the derived radii by a factor of M*/16 M⊙, but the relative brightnesses would remain the same.

In the text