2.5D globaldisk oscillation models of the Be shell star ζ Tauri I. Spectroscopic and polarimetric analysis
Key Words.:
Methods: numerical – Radiative Transfer – Stars: emissionline, Be – Stars: individual: Tau – Stars: kinematics and dynamicsAbstract
Context:A large number of Be stars exhibit intensity variations of their violet and red emission peaks in their H i lines observed in emission. This is the socalled phenomenon, usually explained by the precession of a onearmed spiral density perturbation in the circumstellar disk. That globaldisk oscillation scenario was confirmed, both observationally and theoretically, in the previous series of two papers analyzing the Be shell star Tauri. The vertically averaged (2D) globaldisk oscillation model used at the time was able to reproduce the variations observed in H, as well as the spatially resolved interferometric data from AMBER/VLTI. Unfortunately, that model failed to reproduce the phase of Br 15 and the amplitude of the polarization variation, suggesting that the inner disk structure predicted by the model was incorrect.
Aims:The first aim of the present paper is to quantify the temporal variations of the shellline characteristics of Tauri. The second aim is to better understand the physics underlying the phenomenon by modeling the shellline variations together with the and polarimetric variations. The third aim is to test a new 2.5D disk oscillation model, which solves the set of equations that describe the 3D perturbed disk structure but keeps only the equatorial (i.e., 2D) component of the solution. This approximation was adopted to allow comparisons with the previous 2D model, and as a first step toward a future 3D model.
Methods:We carried out an extensive analysis of Tauri’s spectroscopic variations by measuring various quantities characterizing its Balmer line profiles: red and violet emission peak intensities (for H, H, and Br 15), depth and asymmetry of the shell absorption (for H, H, and H), and the respective position (i.e., radial velocity) of each component. We attempted to model the observed variations by implementing in the radiative transfer code HDUST the perturbed disk structure computed with a recently developed 2.5D globaldisk oscillation model.
Results:The observational analysis indicates that the peak separation and the position of the shell absorption both exhibit variations following the variations and, thus, may provide good diagnostic tools of the globaldisk oscillation phenomenon. The shell absorption seems to become slightly shallower close to the maximum, but the scarcity of the data does not allow the exact pattern to be identified. The asymmetry of the shell absorption does not seem to correlate with the cycle; no significant variations of this parameter are observed, except during certain periods where H and H exhibit perturbed emission profiles. The origin of these socalled triplepeak phases remains unknown. On the theoretical side, the new 2.5D formalism appears to improve the agreement with the observed variations of H and Br 15, under the proviso that a large value of the viscosity parameter, , be adopted. It remains challenging for the models to reproduce consistently the amplitude and the average level of the polarization data. The 2D formalism provides a better match to the peak separation, although the variation amplitude predicted by both the 2D and 2.5D models is smaller than the observed value. Shellline variations are difficult for the models to reproduce, whatever formalism is adopted.
Conclusions:
1 Introduction
Following Struve’s (1931) picture, shell stars are classical Be stars observed with their disks edgeon. The circumstellar disks lie between the observer and the central star (whereas the disks do not project on the photosphere in the case of classical Be stars) and this results in distinct phenomenology. For instance, their spectra exhibit sharp absorption profiles in Balmer lines and in lines from heavier elements that are ionized once. These socalled shell lines are expected to be formed in the innermost section of the rotating disk. Their narrow profiles are the consequence of the extremely small differential velocity projected in the line of sight toward the star (where the rotation is perpendicular to the line of sight).
In the simplest case of an axisymmetric, unperturbed disk, both emission and shellline profiles would be timeinvariant. This is certainly the case for some wellknown objects like CMi for the classical Be stars or Aqr for the shell stars (Hanuschik et al., 1996; Rivinius et al., 2006). However, some shell stars exhibit large variations in the shape, intensity, and radial velocity of their central core, manifesting the existence of timedependent asymmetries and nonazimuthal velocity fields in the disk. As such, shell lines are excellent diagnostic tools to probe Be disk dynamics, as was proven by the analyses performed by Štefl et al. (2012) on 48 Lib and by Mon et al. (2013) on EW Lac.
The shell star Tauri (HR 1910, HD 37202) has been widely studied in the past. The first reason is probably because it is a bright object, whose ideal position in the sky makes it observable from both the northern and southern hemispheres, thus providing a rich dataset. The second and certainly most important reason is because Tauri exhibits various peculiarities that make it the perfect laboratory for understanding the physics of Be stars and for testing theoretical models. Its most important feature is the socalled variation, i.e., the often cyclic fluctuation of the ratio between the violet and red emission peak intensities, in H i lines observed in emission. These variations are thought to be the manifestation of a globaldisk oscillation (Kato, 1983), provoking important modifications to the density and velocity structures in the disk.
In Štefl et al. (2009, hereafter ZT I), a large pool of spectroscopic and spectrointerferometric data covering three complete cycles was gathered and analyzed. In Carciofi et al. (2009, ZT II), a vertically averaged (hereafter 2D) globaldisk oscillation model of Okazaki (1997) and Papaloizou et al. (1992) was combined with the radiative transfer code HDUST to model this dataset. They successfully reproduced the observed cycle of H, as well as the spatially resolved AMBER/VLTI interferometric data, providing the first observational and theoretical confirmation of the existence of global density waves. However, the models shown in ZT II were not able to reproduce the cycle of Br 15 and its polarimetric variations. These two observables originate in the inner disk; it was thus hypothesized at the time that the model geometry in these regions was incorrect. A possible solution emerged with the 3D disk oscillation model of Ogilvie (2008), who has shown that accounting for the vertical structure of the perturbation strongly modifies the properties of the solution for the perturbation pattern in the disk.
There are multiple objectives of the present paper. The first is to quantify the shellline variations of Tauri that have been, until now, only considered in first order, global quantities. We thus propose a finer grained quantitative description of various aspects of the line profile, not limited to the emission components alone. By doing so, we can assess the diagnostic potential of various line features for probing the disk structure and dynamics. The second objective is to model these variations and test whether or not the parameters determined in ZT II are consistent with the observational diagnostics introduced here. The final objective is to test the model of Ogilvie (2008) in a modified 2.5D version and compare its predictions with the former 2D formalism.
In Sect. 2, we briefly present the observational data and give some details about the various line characteristics we measured. Section 3 is dedicated to the quantitative description of the observed temporal shellline variations and their relation with the cycles. In Sect. 4, globaldisk oscillation models will be compared to observations. Before concluding in Sect. 6, the structure of the models and the effect of some model parameters on the observables will be discussed in Sect. 5.
2 Observational data and analysis method
We focus our analysis on the hydrogen lines across the optical and nearIR ranges because these are the lines our numerical tools can model (see next section). The observational material is the same as that used in ZT I and ZT II. It consists of around 300 spectra, combining observations of our own as well as archival data from ESO and from the BeSS database (Neiner et al., 2007). The dataset includes three complete cycles of approximately four years each observed between 1997 and 2010. An inspection of the BeSS data observed after 2010 shows that the strength of the disk emission has drastically decreased and at the same time the cycle has ceased so that since then.
Figure 1 (top panel) shows the curve for H and H between 1997 and 2010, highlighting the quasicyclic nature of the phenomenon. The regions in gray correspond to the triplepeak phases (see below). Although the three cycles (labeled I, II, and III) have slightly different lengths, we followed the definition of ZT I for the phase , taking the origin at the very first maximum (MJD 50414) and an average cycle length of 1429 days.
Figure 2 presents a sample of Balmer lines observed at different phases; the vertical dotted lines mark the rest wavelength (null radial velocity). The most striking features are the following:

: This phase corresponds to the maximum. The violet emission peak is extremely strong in H and H while it is much weaker, but still observable, in H and H. The red emission component is absent from the last two lines. The central shell absorptions are all redshifted.

: These profiles are observed at a phase where , moving toward the minimum. Red emission peaks are a bit stronger than the violet ones in H and H. A faint emission is observed on the violet side of H and H. The shell profiles are all centered on the rest velocity ( km s) and are deeper than in any other phase.

: This phase corresponds to the minimum; the observed features are qualitatively the opposite to the ones described for : the red emission peaks are stronger than the violet ones and the shell absorptions are blueshifted with respect to the rest velocity. Additionally, the shell components appear broader and much shallower. We note that all the lines exhibit a somewhat perturbed profile on their red side. For instance, H presents a central depression in its red emission component and, at nearly the same radial velocity, H, H, and H present a peak. These specific features will be discussed in Sect. 5.3.

: This phase illustrates very well the difficulty of measuring the ratio and the peak radial velocities during the transition from to (triplepeak phase), as H exhibits three emission peaks of comparable intensity. Meanwhile, a small dip appears in the H and H blue emission peaks. The central shell absorptions are all redshifted and are shallower than during previous phases.
This quick overview provides a qualitative description of the phenomenon. We performed a more quantitative approach by measuring various line profile parameters for each spectrum, illustrated on an H profile in Fig. 1 (bottom panel). The values plotted below for H and Br 15 are those from ZT I, they were completed by new measurements on H, H, and H. All the quantities were derived by fitting a Gaussian curve to the top ( bottom) of the violet and red emission peaks ( central absorption). In the case of a triplepeaked emission profile, the Gaussian was applied to the red component as a whole without considering the additional absorption, making the measurements more uncertain in these specific cases.

V/R ratio: The ratio can be defined in two different ways, the first considering the peak height above the continuum level (normalized to unity in the present study), the second considering the peak height above the zero flux level. That second more straightforward definition was adopted here. Hence, the values correspond to the ratio between () and ().

Depth of the central absorption: The depth of the central core () was defined as the difference between the local continuum level (equal to unity) and the minimum value at the bottom of the line.

Position/Radial velocity: The respective heliocentric radial velocities of the central absorption ( at ), violet ( at ), and red ( at ) emission peaks were measured.

Asymmetry: Based on the recent work of Mon et al. (2013), we introduced a parameter to quantify the degree of asymmetry of the shell absorption. This parameter is defined as the ratio between the width of the blue wing at half maximum () measured from the line center and the FWHM of the central absorption (). A value of 0.5 corresponds to a symmetric profile. Higher values indicate that the blue edge is broader (i.e., farther from the line center) than the red edge; smaller values indicate that the red edge is the broader one.
3 Observed spectroscopic variations
3.1 Emission line features
The measured variations were presented in Fig. 1 (upper panel) and repeated in Fig. 3a as reference for the other observables discussed below. Only H (black) and H (red) are plotted as only they exhibit emission peaks pronounced enough to allow meaningful values to be measured for entire cycles.
A secondary maximum, whose position coincides with the triplepeak phase, can be identified for H in Cycle II and more marginally in Cycle III. Unfortunately, the lack of measurements during these periods does not allow a secondary maximum for H to be identified. The pattern observed in H is of smaller amplitude and occurs slightly earlier than that of the H line. Time lags between the variations of different emission lines (Balmer, Brackett, even He i and heavier elements) were highlighted in ZT I and Wisniewski et al. (2007). Similar features were also found in EW Lac (Kogure & Suzuki, 1984; Mon et al., 2013). Since the emission lines probe different radial zones of the unperturbed circumstellar disk, the time lag between their respective cycles is evidence for the existence of a spiral density wave and, as will be shown in Sect. 4, an excellent diagnostic tool for the globaldisk oscillation model.
Despite the large scatter in the measurements, a sawtooth pattern can be recognized for the peak separation (Fig. 3b, black solid lines). This is especially visible in Cycles I and II, where the variation seems to follow the variation, with a possible phase shift. In Cycle III, because of the even more extreme scatter during the triplepeak phase, the distribution looks like a double sawtooth, with the second peak exhibiting a maximum separation larger than the first.
It was mentioned previously that the largescale variations vanished after 2010. Some hints of a weakening density perturbation are found in Figs. 3a and 3b. The first piece of evidence is the decreasing amplitude, from cycle to cycle, of the variations. The second is the increasing peak separation at the maximum; the larger velocity difference between the two peaks is the likely consequence of a disk progressively moving toward a lower density state by which the onearmed oscillation moves out of resonance.
3.2 Shell line features
The violet and red emission peaks are not the only ones to experience velocity shifts; the central shell component exhibits oscillations of its radial velocity apparently synchronized with the variations (Fig. 3c). When the emission peaks are the farthest apart (i.e., at the maximum), the central core reaches its maximum redshift (v km s and km s at the first and second maximum, respectively). Afterward, as the emission peaks move closer to each other and the ratio decreases, the central core progressively shifts blueward. It reaches its maximum blueshift at the minimum, before shifting back toward red. An extremely large blueshift (down to km s) occurs during the triplepeak phase of Cycle III. One can wonder whether these velocity shifts are not just artifacts resulting from the fillingin effect of the lines on their blue side (at maximum) or on their red side (at minimum). This might be conceivable for lines experiencing very strong variations, like H, but less obvious for H and H in which the fillingin effect is much weaker. In addition, the similar velocity shifts that are observed in lines from heavier elements that do not experience any variations, like SiII and HeI (see Fig. 3 of Štefl et al. 2009), favors the scenario in which the central absorption does experience a shift from its rest velocity, thus confirming that the phenomenon is not only a matter of density structure, but also that of associated velocity field.
All the shell lines see their depth (Fig. 3d) dropping very quickly during the maximum before going back rapidly to their usual level. They all exhibit exceptionally shallow profiles during the triplepeak phase of Cycle III, suggesting that the origin of the triple peak may actually be an emission phenomenon filling in the central absorption, instead of an extra absorption as hypothesized by ZT I.
We can barely distinguish an overall consistent cyclic behavior for the shell line asymmetry (Fig. 3e). The mean asymmetry parameter is 0.545 for H (red), 0.562 for H (blue), and 0.565 for H (green), which indicates that all the lines are on average slightly broader on their blue side. This blueward asymmetry seems to increase progressively for higher members of the Balmer series; however, the difference from one line to another is extremely small. Intriguingly, the asymmetry of H and H are strongly and differently affected during the triplepeak phase of Cycle III, the former becoming broader on its blue side and the latter on its red side; H stays surprisingly unperturbed during that period.
We conclude that the depth and position of the shell absorption offer promising diagnostic potential for the phenomenon, whereas the shellline asymmetry does not, as it seemingly does not correlate with the cycle. This last statement should nonetheless be moderated by the fact that shell profiles are, on average, broader on their blue edge, which should be accounted for while testing our models as it is another hint of the spiral pattern.
4 Modeling the data
4.1 Tools and method
The computer code HDUST (Carciofi & Bjorkman, 2006) was used to model the observed features described above. It is designed to solve the coupled problems of radiative transfer, radiative equilibrium, and statistical equilibrium in a nonlocal thermodynamical equilibrium (NLTE) regime. A Monte Carlo approach is adopted: once emitted by the photosphere, the path of each photon across the stellar envelope is followed. The radiative equilibrium is ensured because the absorbed photons are reemitted with a direction and a wavelength determined by the local conditions. Both continuum and line processes are included in order to compute the opacity and the emissivity of the gas.
The highly variable nature of Tauri requires a methodic approach and the combination of a disk oscillation model with HDUST. The modeling campaign from ZT II was performed in two steps. The first step consisted in fitting the average observed properties (optical and infrared SEDs, continuum polarization H peak separation) in order to fix the stellar and average disk parameters. The theoretical background describing the structure of the unperturbed disk (see Appendix A) assumes a steadystate viscous decretion disk in hydrostatic balance in the vertical direction and truncated by the binary at the tidal radius (18.6 for Tauri). The second step adopted the 2D globaldisk oscillation model of Okazaki (1997) in order to reproduce the timedependent properties of Tauri: on the unperturbed state was superposed a linear perturbation in the form of a normal mode of frequency .
In this paper we adopt the model proposed by Ogilvie (2008) for describing the dynamical structure of the disk, and combine it with HDUST for computing line profiles and other observables. The model of Ogilvie (2008) includes the hitherto ignored component (in cylindrical coordinates , , ) in the calculations. This 3D formalism requires the perturbed quantities to be expanded in the direction by means of Hermite polynomials (see Appendix B). A positive aspect of these 3D models is that, in contrast to their 2D counterparts, they are able to produce prograde confined perturbation modes, as expected on both theoretical (Papaloizou et al., 1992; Oktariani & Okazaki, 2009) and observational grounds (Telting et al., 1994; Vakili et al., 1998), without using any weak line force (Chen & Marlborough, 1994).
Since this new model predicts a complex structure for the circumstellar disk, we adopted in this work a 2.5D approach, resulting from the integration of the 3D perturbed disk structure in the direction. This simplified approach is a first step toward a planned full 3D treatment. It allows a meaningful comparison between the new theory and the previous 2D formalism. In practice, the method consists first in solving the 3D equations (B.19 to B.21) describing the perturbed disk, the solutions of which are expressed by equations B.24 to B.27. Afterward, by adopting a null value for the secondorder terms ( and h), the vertical components are neglected and only the zerothorder perturbed quantities, which are independent of the vertical coordinate, remain. As will be seen below, this 2.5D model already brings various changes to the model structure and, consequently, in the resulting observables. It qualitatively improves our modeling capacities, thus building a useful bridge between the previous 2D and the upcoming 3D models.
4.2 Presentation of the models
Our reference model (Model 1) was computed with the 2D formalism and the stellar, orbital, and geometric parameters constrained in ZT II (their Table 1). Ogilvie (2008) noted that 3D models are naturally more confined than their 2D counterparts. The degree of confinement of the perturbation mode is defined by the viscosity parameter , first introduced by Shakura & Sunyaev (1973). To allow a direct comparison between the two formalisms, Model 2 was run with the same parameters as Model 1 and a 2.5D structure for the perturbed disk. In order to test the effect of the viscosity parameter, various 2.5D models with different values of were computed. The last model discussed below (Model 3) was computed with the same parameters as the other two and (instead of 0.4). It was chosen because it provides a density wave confinement, materialized by the amplitude of its variations, comparable to the one predicted by the reference model (Sect. 4.3 and Fig. 5). The respective disk structure of each model will be discussed further in Sect. 5.1; their main differences are summarized in Table 1.
To make the comparison with our models easier, the observational data were converted in phase based on the ephemerides measured in ZT I. The phase is usually defined such that at the maximum. To simulate the temporal variation the observables experience while the density wave precesses around the star, for a given model we compute the emerging spectrum for several azimuthal coordinates . The phase is related to the model phase by
where is the phase shift applied to each model in order to match the observational phase (see Table 1).
An illustration of the Balmer line profiles predicted by all three models is given in Fig. 4, at phases close to those presented in Fig. 2. The features are qualitatively the same as those described in Sect. 3. The most obvious difference is the absence of triplepeaked profiles in the models, suggesting that whatever the physical process underlying this phenomenon is, it is neither included in nor predicted by our current models. This idea is also supported by the very similar central absorption depths predicted by the models on phases and , contrasting with the large differences observed in Fig. 2 between phases and . The same quantities presented in Fig. 1 (bottom panel) were measured on the model spectra, and compared to the observational ones in Figs. 5  7.
parameters  

Stellar and disk parameters  Formalism  Weak line force  Viscosity parameter  
Model 1  ZT II  2D  280  0.006  5.7410(/)  0.4 
Model 2  ZT II  2.5D  35  0.008    0.4 
Model 3  ZT II  2.5D  55  0.009    0.8 
4.3 Emission line features
The properties of the emission lines are essential diagnostic tools for our models because they probe the physical conditions throughout the disk. The variations and peak separation predicted by our models are plotted in Fig. 5 and compared to the observed ones. We chose H to illustrate our results because the larger dataset allows a better comparison. Nevertheless, all the comments below apply identically to H. Model 1 (from ZT II) provides a reasonable match to the variations. Switching to the 2.5D formalism and keeping all the parameters fixed (Model 2) results in a more confined density wave, as demonstrated by a curve whose shape and amplitude do not match the observations. This was fixed by increasing the viscosity parameter (Model 3). The resulting shape is quite different, however, as Model 3 reaches the minimum sooner than Model 1 and presents a sort of plateau between and . During the phases, the ability of Model 1 and Model 3 to reproduce the observations does not seem qualitatively different. During the phases, Model 3 predicts a sharper profile, more in agreement with the observations. As expected from the line profiles plotted in Fig. 4, the models are not able to explain the triplepeak phases (marked as grayshaded areas in this and the following plots for reference).
Reproducing the peak separation variations is more problematic. All the models exhibit variations whose amplitude is smaller than what is observed. The 2D model (Model 1) is the one that matches best the sawtooth pattern observed in Cycle I and Cycle II. The 2.5D models predict more perturbed variations, especially Model 3, which exhibits a secondary minimum that looks like the double sawtooth pattern observed in Cycle III but with a much smaller amplitude.
4.4 Shell line features
While the emission lines sample the disk as a whole, the region probed by the shell lines is restrained to the line of sight. The diagnostics they provide should not be neglected, as their properties (absorption depth, central position, asymmetry) reflect the structure and dynamics of the disk in the observer’s direction. A comparison between our models and the observed shellline properties is proposed in Fig. 6, for H.
Looking first at the position of the central absorption (top panel), there is an obvious discrepancy between the model predictions and the observations. First, the models are displaced in phase while for the same value of they match the phase of H perfectly. Second, the average position (similar for each model 47 km s) is at a higher velocity than the observed ( 11 km s). Finally, the amplitude of the variations is clearly far too small in the models.
If we exclude the triplepeak phases, the models exhibit variations in the shell absorption depth whose amplitude is compatible with the observations. The scarcity of the observational data does not allow the variation pattern predicted by the models to be validated.
The results for the line asymmetry (bottom panel) also highlight strong challenges for the models; they predict largeamplitude (between 0.2 and 0.8) cyclic variations, while the observations show very little variation (amplitude 0.1) and no sign of cyclic behavior. Model 3 predicts shell lines that are on average slightly broader on their blue side (average asymmetry parameter ), as might be expected from the results in Sect. 3.2, whereas Model 1 predicts more symmetric shell line profiles (average asymmetry parameter ) and Model 2 shell lines slightly broader on their red side (average asymmetry parameter ). The differences between the three models are marginal, however, considering the typical error on the asymmetry parameter ( for the spectral resolution of our models).
Overall, these results highlight a limitation of our models in their current form. The 2D model was shown by ZT II to offer a good description of the spectrointerferometric data observed by VLTI/AMBER, which are sensitive to the onthesky distribution of the disk intensity. This and the good overall fit to the variations suggest that the models reasonably reproduce the density distribution of gas across the disk. However, the line features of Tauri, considered in such detail for the first time, are influenced not only by the density distribution, but also by the disk kinematics. In particular, the shell features are sensitive to the lineofsight projected velocity of the disk. Therefore, the general failure of both the 2D model by Okazaki (1997) and Papaloizou et al. (1992) and the present 2.5D model based on the works of Ogilvie (2008) to reproduce them suggests that both models are predicting the wrong kinematic structure for the disk.
4.5 IR lines and polarization
Two important issues were highlighted in ZT II concerning the 2D globaldisk oscillation model. The first was its inability to reproduce the phase of the cycle of Br 15. The model reproduced the observed amplitude, but had to be artificially shifted by 0.2 in phase in order to match the observations. The second was about the polarization: the model predicted variations whose amplitude was much larger than the observed one. Since both the IR lines and the polarization originate from the innermost regions of Be disks, it was concluded in ZT II that the 2D model was unable to predict the correct disk structure in these regions. This was supported by the work of Ogilvie (2008), who showed how much stronger the confinement of the oscillation modes is when a 3D perturbation is adopted.
Figure 7 (top panel) compares the three aforementioned models with the observed Br 15 cycle (data from ZT I). Models 1 and 3 both reproduce the amplitude of the variations, but Model 3 clearly does a better job at reproducing the phase. Model 2 obviously presents the worse fit quality, as was the case for the H cycle. It was evoked previously that the time lag, or equivalently the phase shift, between the maxima of different emission lines is an important diagnostic tool for the globaldisk oscillation model. The ability of the 2.5D formalism to consistently fit the variations of H and Br 15 represents a great step in our understanding of the phenomenon, as it puts some constraints on the properties of the spiral arm perturbation, further discussed in Sect. 5.
We now consider the continuum polarization (Fig. 7, bottom panel). The observations show an average level of 1.45 and a standard deviation of 0.07. Model 2 is the model providing the smallest polarization level (1.12 in average). Model 1 and Model 3 exhibit a higher level (1.28 and 1.26, respectively) but remain below the observations. This result tends to favor less confined models, but we see in Sect. 5.2 that other factors can produce an average level comparable to the observed one. Regarding the amplitude of the variations, Model 2 matches the observations better, while Model 1 and Model 3 predict too large large (0.14) and too small variations (0.05), respectively. To summarize, none of these models provides a satisfactory fit to the polarization data, as they can either match the average level or the amplitude, but not the two simultaneously. On a more qualitative note, the shape of the polarization variations predicted by the 2.5D models present two maxima per cycle, as expected on observational grounds, while the 2D formalism presents a centered, singlemaximum pattern.
Model  H (R)  (g cm)  (log)  (km s)  v (km s)  Pol(V)  

Model 1  0.208  5.9010  85  0.4  0.495  218.5  24.1  1.284  0.279 
Model 2  0.208  5.9010  85  0.4  0.209  230.8  23.8  1.122  0.199 
Model 3  0.208  5.9010  85  0.8  0.384  215.3  41.2  1.256  0.098 
0.208  7.0810  85  0.8  0.348  209.2  34.7  1.269  0.097  
0.208  8.8510  85  0.8  0.314  202.7  22.9  1.291  0.119  
0.208  1.1810  85  0.8  0.317  196.3  18.1  1.314  0.146  
i  0.208  5.9010  83  0.8  0.319  209.6  27.3  1.357  0.114 
i  0.208  5.9010  87  0.8  0.447  220.9  61.6  1.129  0.113 
H  0.222  5.9010  85  0.8  0.397  216.8  44.3  1.304  0.107 
H  0.236  5.9010  85  0.8  0.416  218.9  46.0  1.344  0.112 
H  0.249  5.9010  85  0.8  0.426  220.4  48.1  1.390  0.115 
5 Discussion
5.1 Structure of the solution
In the previous sections, we showed that switching from the 2D to the 2.5D formalism and increasing the viscosity parameter, , imply large modifications on the observables computed by our models. In this section we discuss how these modifications reflect important changes in the model disk structure, especially regarding the shape of the spiral density wave.
An illustration of these structural differences is given in Fig. 8, where the disk density perturbation pattern is represented for the three models detailed previously. The shading indicates the amplitude and sign of the perturbation, from density deficit (lighter colors) to density enhancements (darker colors). A simple visual comparison between Model 1 and Model 2 (left and center panels, respectively) shows qualitatively how the strength of the confinement of the density wave changes while moving from the 2D to the 2.5D formalism: in the former, the spiral arms extend up to while in the latter they are restrained within .
It is only by adopting a high value of that the 2.5D formalism produces density perturbations able to spread at larger radii (Model 3, right panel). That parameter was introduced by Shakura & Sunyaev (1973) to characterize the efficiency of the angular momentum transport in a turbulent, differentially rotating medium. Higher values of imply a more efficient diffusion of the density wave, explaining why the perturbation spread farther in Model 3. The results of Shakura & Sunyaev (1973) suggested that cannot be greater than 1 in the case of turbulent viscosity, so adopting a value as large as 0.8 in our models can almost be seen as a limit case (see, e.g., the recent determination of for the disk of 28 CMa made by Carciofi et al. 2012).
The spiral arm perturbation can be characterized not only by its extent/confinement, but also by how tight it is wound around the star. On that score, it can be seen in Fig. 8 that the 2.5D models exhibit more coiled spiral patterns than the 2D model. Keeping in mind that ZT II defined as the azimuth of the minimum of the density perturbation pattern at the base of the disk, the different coil from one model to another is what explains the large phase shift seen between Model 1 () and Models 2 and 3 ( and , respectively). In addition, the tighter coil of the spiral arm of Model 3, compared to Model 1, explains why the first does a better job at fitting the phase of Br 15 than the second.
5.2 Diagnostic potential
Overall, the previous results indicate that the 2.5D formalism provides a better representation of the globaldisk oscillation phenomenon, at least regarding the intensity diagnostics, such as the curves of H and Br 15. It remains less conclusive regarding the kinematics diagnostics, such as the shell features.
Various model parameters determined in ZT II were kept fixed for the analysis presented above. Among these parameters, some may strongly affect the onthesky projected disk density, namely the base density of the disk, its vertical scale height, and its inclination angle with respect to the observer. The impact of these three quantities on the observables are now studied to test whether the values determined in ZT II remain the best ones in the 2.5D formalism framework.
Model 3 (2.5D with ) will be taken as reference here; starting from that model, we made the (unperturbed) base density , the scale height and the inclination angle vary. The model parameters are compiled in Table 2 and the results pictured in Fig. 9 for three observables: the H variations (top row), peak separation (middle row), and the continuum polarization (bottom row).
We consider first the unperturbed base density (left column, blue), i.e., the density (expressed in g cm) in the equatorial plane of the disk and at its innermost radius. The models presented in Sect. 4 were computed with = 7.910 g cm; we increased that value by 20, 50, and 100%. Doing so tends to progressively reduce the amplitude of both variation and peak separation, hence degrading the fit to these two observables. On the other hand, higher densities tend to match the polarimetric data slightly better, although without reaching the observed level. As highlighted by Haubois et al. (2014, their Fig. 4), the relation between polarization and density is not straightforward as it also depends strongly on the effective temperature of the central star. For tenuous disks, the polarization level increases linearly with the density, whatever the effective temperature is. For denser disks, the growing contribution of multiplescattering and prescattering absorption leads to a nonlinear variation of the polarization with base density. Zeta Tauri seems to lie in a regime where raising its base density by a large amount does not affect its polarization level significantly.
In the isothermal viscous decretion disk model (Carciofi & Bjorkman, 2006), the radial density falloff is a power law with a slope of and the vertical density structure is a Gaussian with a scale height given by
where,
is the scale height at the base of the disk, the sound speed, the critical velocity, and for isothermal disks. In ZT II, was fixed to 0.208 by adopting a gas temperature of about 70% of the stellar effective temperature on the equatorial plane. For our tests, we increased up to 0.249 , corresponding to a gas temperature equal to the effective temperature. The effect on the observables can be summarized as follows: (i) the curve appears little affected; (ii) the maximum peak separation slightly increases, but still remains well below the observed values; (iii) the average polarization increases to a level comparable to the observations while the amplitude of the variations remains almost unchanged, and thus is a bit too small.
The inclination angle was determined in ZT II through the fits to SED and continuum polarization. It was shown at the time that a difference as small as 2 greatly affects the polarization. This is confirmed in Fig. 9 (right column, red): the average polarization level decreases ( increases) significantly for slightly higher ( lower) inclination, with respect to the reference model. A lower inclination angle (, dotted line) seems more consistent with the observed polarization level, but it slightly degrades the fit to both the variations and peak separation. Instead, a higher inclination angle (, dashed line) provides larger values of the maximum peak separation, almost comparable to the observed maximum, but then the polarization level becomes far too low.
In summary, it appears possible to improve the match with the polarization data by tuning the three aforementioned quantities, but this implies degrading the fit to the observed variations. It is also possible to improve the match with the observed peak separation, at the cost of making the match with the polarization data even worse. From those considerations, the base density and inclination angle determined in ZT II remain the best parameters to fit the quantities. This was not unexpected since they were determined independently and without any former consideration for the disk oscillation formalism.
The scale height of the disk has been explored here for the first time. Its strong impact on the polarization level together with its extremely light impact on both H variations and peak separation might indicate that its influence is limited to the innermost regions of the disk. Based on the polarization diagnostics, it appears that values of higher than those adopted in ZT II are favored.
5.3 The H profile: triplepeak or doubledip?
In the two previous sections, we pointed out the various unusual features observed during the triplepeak phases. The origin (i.e., the physical process giving rise to those profiles) and nature (three distinct emission components selfabsorption overlapping with an emission component) of the triplepeaked profiles are puzzling issues that, together with the shellline profile perturbations highlighted here, deserve some discussion. Various systems, like Be/Xray binaries (Moritani et al., 2011), exhibit spectroscopic features that look similar but whose origin might differ from those observed in Tauri. The discussion below likely applies to objects like Tauri, namely Beshell binaries undergoing variations, and does not exclude another origin for triplepeaked features observed in other Be stars.
The spectroscopic diagnostics provide somewhat contradictory clues regarding the origin of the triplepeaked features. For instance, ZT I noted that other nonhydrogen emission lines in the visual range do not experience a triplepeak phase. On the one hand, the absence of triplepeaked profiles in optically thin lines like O i 8446, which probes the Hforming region due to Ly resonance pumping (Andrillat et al., 1988), suggest that its origin is probably not an actual local density enhancement. Considering the shell nature of Tauri, a selfabsorption located at the base of the disk and overlapping with the emission component might be invoked. On the other hand, the absence of triplepeaked profiles in optically thick lines, like Fe ii 6318, suggests that the feature originates from the outer disk, since those lines are supposed to be formed closer to the star than H and H. This would support a binary origin, as the companion would more likely affect the outermost regions of the disk.
Observable  Probed Region  Related Parameter(s)  Favored Model 

(H)  Entire Disk  , Base Density, Inclination Angle  2.5D 
Peak Separation (H)  Entire Disk  , Base Density, Inclination Angle  2D 
Shell Lines  Line of Sight     
Polarization  Inner Disk  , Base Density, Inclination Angle, Disk Scale height  2.5D 
(Br 15)  Inner Disk  2.5D 
The results presented in Sects. 3.1 and 3.2 support the idea that, whatever mechanism is creating the triplepeaked profiles,
it should be effective both in the outer disk, from where the emission lines originate, and in the inner disk, from where the shell
lines originate. On that score, the triplepeak phase of Cycle III is particularly interesting: it is the shortest one,
Identifying the nature of the triplepeaked profiles might help us to better understand their origin. Traditionally, disturbed H profiles like the ones displayed by Tauri around = 0.75 have been called triplepeaked profiles, for obvious reasons. However, the density modulation induced by a onearmed disk oscillation is very smooth (see Fig. 8). It provides little inspiration for the delineation of three distinct lineemitting regions. The general success of the model even for Tauri suggests that the presence of the companion does not change this significantly.
However, restricting the search for an interpretation to absorption could lead to an increasing number of additional free parameters for our models. For instance, there could be selfabsorption in a highdensity region or the companion star could excite a warping mode of the disk, which would require not only an extension of our code to 3D, but also an expansion of its dynamical basis. However, there is an ansatz, which does not add any new parameters: the doubledip phenomenon, which it will be discussed here, is most prominent in H.
Because it is almost a resonance line, H absorption forms along the line of sight to the central star over a very large range in distance, larger than for any other Balmer line. If the effect of the onearmed disk oscillation is pictured as a spiral pattern of density enhancement (Fig. 8), it is conceivable that the line of sight intersects the increaseddensity zone of the spiral twice, at about the same azimuth (modulo 360) but different radii. Of the resulting two discrete absorption components, one would be roughly shared with the highorder Balmer lines forming in the inner disk while the latter would exhibit at most traces of the second absorption in H. Because the spiral is onearmed, double dips only occur once per V/R cycle. This description seems to be in qualitative agreement with the observations of Tau. In this scenario, double dips would always occur at roughly the same phase in all stars, namely when the outer Habsorbing end of the spiral crosses the line of sight, which corresponds to the middle of the ascending branch of cycles ().
The apparent restriction of the doubledip phenomenon to shell stars with a companion (Rivinius et al., 2006; Štefl et al., 2007) could imply that the companion not only truncates the disk but also changes its dynamics. For instance, the disk eigenmodes could be affected or the viscosity could be increased, leading to a more strongly wound spiral pattern (see Fig. 8). Cycletocycle differences in the appearance of doubledip line profiles, as observed in Tau, may supply additional useful interpretative constraints. A forthcoming paper will more specifically investigate the physical basis, if any, of this picture and explore more quantitatively the diagnostic value of doubledip shell absorptions.
6 Conclusions
We undertook the analysis of the Be shell star Tauri, introducing new criteria to quantify the variations of its spectroscopic characteristics. The quantitative analysis of its observational data was complemented with a detailed modeling, adopting a new formalism based on the results of Ogilvie (2008). The main observational results are the following:

As expected, the variations are accompanied by variations of the peak separation; the latter exhibiting a sawtooth pattern identifiable despite the large scatter among the observational points.

Shell lines exhibit marginal variations in depth and more obvious variations in velocity. Their apparent synchronism with the variations make them good diagnostic tools for our models. Surprisingly, the asymmetry of the shell lines does not vary much. The average value for the asymmetry parameter indicates that shell lines are generally broader on their blue side.

Between 1997 and 2010, Tauri experienced three complete cycles. Each cycle exhibits a period during which the emission profiles of H and H are perturbed, making the measurement of the ratio and peak separation uncertain. The origin and nature of this socalled triplepeak phase has been discussed and a possible explanation, the doubledip, has been suggested. The existence of such features only in Beshell stars with a companion suggest that the physical mechanism creating these profiles may be triggered (or enhanced) by the companion itself.
From our modeling efforts, we reached the following conclusions (also summarized in Table 3):

The new 2.5D globaldisk oscillation formalism is generally preferred over the previous 2D model. In particular, the 2.5D formalism solves the Br 15 phase problem found by Carciofi et al. (2009) and improves the agreement with the observed continuum polarization if a larger disk scale height is adopted. This is an encouraging first step toward our forthcoming 3D formalism.

Our models (both 2D and 2.5D) fail to reproduce most of the shellline properties. They predict variations that are too small for the position of the shell absorption component and large asymmetry variations that are not observed. They were thus not used as diagnostic tools for constraining model parameters, but they highlighted a current limitation of our models that we shall explore in the near future.

As expected by the theory, for a given value of the viscosity parameter , 2.5D models produce much more confined spiral density perturbations than their 2D counterparts. Since reproducing the variations requires less confined perturbation modes, larger values of must be adopted while computing 2.5D models. This finding is in agreement with the recent determination of a large () for 28 CMa by Carciofi et al. (2012).

The effect on the observables of some model input parameters has been assessed. The base density and inclination angle determined in ZT II remain valid within the 2.5D formalism. The disk vertical scale height has been explored here for the first time; the tests realized on various observables, the polarization in particular, tend to favor larger values for that parameter.
The present paper is the first of a series. The second is already in preparation and aims at applying the 2.5D model to interferometric data observed with AMBER/VLTI and MIRC/CHARA. As mentioned earlier, the 2.5D formalism discussed here is a useful step toward the 3D treatment, not the final model. The third paper of the series will be dedicated to the 3D modeling of Tauri, gathering all the data presented in the first two papers. The results presented here suggest that the density distribution across the disk is already correctly reproduced by the 2.5D model, but not the kinematic structure. Therefore, the biggest challenge the 3D model will have to face is reproducing the amplitude of the lineofsight velocity variations, and providing a consistent picture of both the density and kinematic structure of perturbed Be disks.
Acknowledgements.
The authors are very grateful to the anonymous referee for the careful review of the manuscript and the useful comments. This work was supported by the FAPESP grant N 2011/179305 to C.E. and made large use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), the purchase of which was made possible by the Brazilian agency FAPESP (grant N 2009/540064) and the INCTA. C.E. and A.C.C. are grateful to the technical computing staff of the astronomy department of the IAG/USP for their continuous help. A.C.C. acknowledges support from CNPq grant N 307076/20121. A.T.O. is grateful to support from the JSPS grant N 24540235.Appendix A Unperturbed disk state
Before we consider the equations describing the 3D formulation of Ogilvie (2008) for the onearmed oscillation model, we shall note the fundamental equations describing the unperturbed disk. The model adopted here is the same as that detailed in ZT II (their Sect. 2). We assume a steady, axisymmetric viscous decretion disk, for which the surface density is given by (Bjorkman & Carciofi 2005)
(1) 
where is the horizontal distance from the center of the star, is the average mean molecular weight, is the mass of the hydrogen atom, is the viscosity parameter of Shakura & Sunyaev (1973), is the average gas kinetic temperature, is the stellar decretion rate, and is an arbitrary integration constant.
The massdecretion rate is related to the radial velocity via the equation of continuity (Carciofi & Bjorkman 2008). It is given by
(2) 
where is the verticallyaveraged, radial velocity.
In binaries with low eccentricities, the Be disk is expected to be truncated by the tidal torques from the companion (Okazaki & Negueruela 2001). For a circular binary with a small mass ratio (), such as Tau, the truncation takes place at the tidal radius, which is approximately given by (Whitehurst & King 1991), with being the Roche radius approximately given by (Eggleton 1983)
(3) 
where is the semimajor axis.
We also assume that the unperturbed disk is in hydrostatic balance in the vertical direction. Then, the unperturbed density is given by
(4) 
where is the disk density at the midplane (), and the disk scale height is given by
(5) 
where is the disk scale height at the inner radius .
We can find the radial dependence of the midplane density by integrating in the direction:
(6) 
Appendix B 3D Structure of oscillations in truncated decretion disks
On the unperturbed state described in Appendix A, we superpose a linear perturbation in the form of a normal mode of frequency , which varies as . For simplicity, we assume the perturbation to be isothermal. The linearized perturbed equations are then obtained as follows:
(9) 
(10) 
(11) 
(12)  
Here is the isothermal sound speed and () the perturbed velocity; and are respectively the unperturbed and perturbed enthalpy, where is the density perturbation.
Following the formulation made by Ogilvie (2008, see also ), in order to take the 3D effect into account, we must expand perturbed quantities in the direction in terms of Hermite polynomials:
(13) 
(14) 
(15) 
(16) 
(Ogilvie 2008; see also Okazaki et al. 1987), where is the Hermite polynomial defined by
(17) 
with and being a dimensionless vertical coordinate and a nonnegative integer, respectively. We note that for .
Replacing the perturbed quantities by their Hermite polynomials expansion, Eqs. (B.1) to (B.4) become
(18) 
(19) 
(20) 
(21) 
We note that Eqs. (18)(21) are reduced to Eqs. (48)(51) of Ogilvie (2008)
provided that the effects of the radial velocity and viscosity are neglected and is
replaced with
Still following Ogilvie (2008), we close the system of resulting equations by assuming for . Using the approximations , , , and , the basic equations for linear perturbations in viscous decretion disks are then
(22) 
(23) 
(24) 
(25) 
(26) 
Eliminating and from Eqs. (22)(26), we obtain a set of three equations for , , and :
(27) 
(28) 
(29) 
Equations (B.19), (B.20), and (B.21) correspond respectively to Eqs. (16), (17), and (15) from ZT II. They are solved with the same boundary conditions, i.e., a rigid wall at the inner disk radius and at the outer disk radius, being the Lagrangian perturbation of pressure.
We note that the terms ( in isothermal disks) in Eq. (27) and in Eq. (28) result from the 3D effect and provides an important contribution to the confinement of the modes. The other terms are the same as in the basic equations used in ZT II.
Once Eqs. (27)(29) are solved, the secondorder terms, and , are obtained from Eqs. (25) and (26),
(30) 
(31) 
where . Finally, under the closure approximation adopted here, the perturbed density and velocity are given by
(32)  
(33)  
(34)  
(35)  
where the unperturbed density is given by Eq. (4) and the density perturbation is related to the enthalpy perturbation by and .
Equations (B.24) to (B.27) provide the 3D density and kinematic structure of the perturbed disk. In the framework of the 2.5D formalism, the secondorder terms, namely and , are taken as equal to zero. Consequently, the component of the velocity field () and the dependent terms of the perturbed density () are neglected. The final solution is thus independent of the vertical coordinate and includes zerothorder terms only.
Footnotes
 It is worth noticing that Cycle III is also the shortest cycle, with 1230 days (see Table 2 of ZT I).
 This results from the difference between our normal mode form, , and that of Ogilvie (2008),
References
 Andrillat, Y., Jaschek, M., & Jaschek, C. 1988, A&AS, 72, 129
 Bjorkman, J. E. & Carciofi, A. C. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 337, The Nature and Evolution of Disks Around Hot Stars, ed. R. Ignace & K. G. Gayley, 75–+
 Carciofi, A. C. & Bjorkman, J. E. 2006, ApJ, 639, 1081
 Carciofi, A. C. & Bjorkman, J. E. 2008, ApJ, 684, 1374
 Carciofi, A. C., Bjorkman, J. E., Otero, S. A., et al. 2012, ApJ, 744, L15
 Carciofi, A. C., Okazaki, A. T., Le Bouquin, J.B., et al. 2009, A&A, 504, 915
 Chen, H. & Marlborough, J. M. 1994, ApJ, 427, 1005
 Eggleton, P. P. 1983, ApJ, 268, 368
 Hanuschik, R. W., Hummel, W., Sutorius, E., Dietle, O., & Thimm, G. 1996, A&AS, 116, 309
 Haubois, X., Mota, B. C., Carciofi, A. C., et al. 2014, ApJ, 785, 12
 Kato, S. 1983, PASJ, 35, 249
 Kogure, T. & Suzuki, M. 1984, PASJ, 36, 191
 Mon, M., Suzuki, M., Moritani, Y., & Kogure, T. 2013, PASJ, 65, 77
 Moritani, Y., Nogami, D., Okazaki, A. T., et al. 2011, PASJ, 63, L25
 Neiner, C., de Batz, B., Mekkas, A., Cochard, F., & Martayan, C. 2007, in SF2A2007: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. J. Bouvier, A. Chalabaev, & C. Charbonnel, 538
 Ogilvie, G. I. 2008, MNRAS, 388, 1372
 Okazaki, A. T. 1997, A&A, 318, 548
 Okazaki, A. T., Kato, S., & Fukue, J. 1987, PASJ, 39, 457
 Okazaki, A. T. & Negueruela, I. 2001, A&A, 377, 161
 Oktariani, F. & Okazaki, A. T. 2009, PASJ, 61, 57
 Papaloizou, J. C., Savonije, G. J., & Henrichs, H. F. 1992, A&A, 265, L45
 Rivinius, T., Štefl, S., & Baade, D. 2006, A&A, 459, 137
 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
 Struve, O. 1931, ApJ, 73, 94
 Telting, J. H., Heemskerk, M. H. M., Henrichs, H. F., & Savonije, G. J. 1994, A&A, 288, 558
 Štefl, S., Le Bouquin, J.B., Carciofi, A. C., et al. 2012, A&A, 540, A76
 Štefl, S., Okazaki, A. T., Rivinius, T., & Baade, D. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 361, Active OBStars: Laboratories for Stellare and Circumstellar Physics, ed. A. T. Okazaki, S. P. Owocki, & S. Stefl, 274
 Štefl, S., Rivinius, T., Carciofi, A. C., et al. 2009, A&A, 504, 929
 Vakili, F., Mourard, D., Stee, P., et al. 1998, A&A, 335, 261
 Whitehurst, R. & King, A. 1991, MNRAS, 249, 25
 Wisniewski, J. P., Kowalski, A. F., Bjorkman, K. S., Bjorkman, J. E., & Carciofi, A. C. 2007, ApJ, 656, L21