α Virginis (Spica, HD 116658) is a double-lined spectroscopic binary consisting of two early B-type stars in a short-period (∼4 days) eccentric orbit. The B1 III?IV primary was discovered to be a β Cephei-type star in early observations, with a pulsation period of 0.1738 days in both the light curve (Shobbrook et al. 1969 ) and radial velocity (RV) variations (Smak 1970 ). Although the light variations became undetectable a short time after they were first detected (Lomb 1978 ; Sterken et al. 1986 ), the line-profile variability persists (Smith 1985a , 1985b ; Riddle 2000 ). Its most striking characteristic is the presence of discrete absorption and emission features traveling from the blue toward the red wing of the absorption line.

The first detection of this type of periodic line-profile variability was made by Walker et al. ( 1979 ) in the rotationally broadened O9.5V star ζ Ophiuchi. In their seminal study, Vogt & Penrod ( 1983 ) showed that this variability could be explained by non-radial pulsations (NRPs). Subsequent observations show a prevalence of this type of line-profile variability (cf. Baade 1984 ; Gies & Kullavanijaya 1988 ; Reid et al. 1993 ; Fullerton et al. 1996 ; Rivinius et al. 2001 ; Uytterhoeven et al. 2001 ). The numerous models based on the NRP mechanism are reviewed by Townsend ( 1997a , 1997b ), some of which incorporate the effects of rapid rotation. A general feature of these models is to assume a prescription for the surface velocity field of the pulsating star, from which the perturbed photospheric absorption line profiles are computed.

A very different approach was used by Moreno et al. ( 2005 ), which calculates the surface velocity field from first principles. This calculated (but not assumed) field is then projected along the line of sight to the observer to produce photospheric absorption line profiles. The one basic approximation in this model is to assume that only the external layer of the star oscillates in response to the forces in the system, while the interior region is assumed to rotate rigidly. This one-layer approximation is essentially equivalent to the assumption that the surface layer behavior is primarily controlled by the binary companion perturbations.

The advantages of this method are threefold: (1) it makes no a priori assumption regarding the mathematical formulation of the tidal-flow structure since the velocity field v is derived from first principles; (2) the method is not limited to slow stellar rotation rates nor to small orbital eccentricities; and (3) it is computationally inexpensive. It is not clear to what extent the response of the inner layers may affect the line-profile calculation because the amplitudes of induced oscillations decrease sharply in deeper layers (Dolginov & Smel'Chakova 1992 ). In this paper, we confront the results obtained from the Moreno et al. ( 2005 ) method with very high quality observational data of Spica. We show that the azimuthal velocity components of the tidal perturbations (the "tidal flows") play a dominant role in determining the line-profile variability as produced using this theoretical framework.

The structure of this paper is as follows: In Sections 2 and 3 , we describe the observations and the characteristics of Spica's spectrum and variability; in Section 4 , we present the results of the model calculations of Spica's surface velocity field and line-profile variability; the conclusions are presented in Section 5 ; and details of the data reduction process and a discussion of the local line-profile effects and size of the perturbations on the stellar surface are given in the Appendices .

The observations were performed in queue mode at the Canada France Hawaii 3.6 m telescope (CFHT) with the ESPaDOnS spectropolarimeter. ESPaDOnS is a fiber-fed cross-dispersed echelle covering 3700?10480 A in a single exposure with 40 spectral orders on the EEV1 detector with a nominal spectral resolution of R = 68,000 in polarimetric mode. The instrument has a dedicated reduction package, Libre-Esprit, that automatically processes the data. This script does typical flat fielding, bias subtraction, wavelength calibration using both calibration lamp and Fabry-Perot frames, and optimal spectral extraction (see Donati et al. 1997 ). Our integration times were 6 s per exposure in polarimetric mode. This mode uses a Wollaston prism to produce two separate orthogonally polarized spectral orders on the detector. A polarimetric sequence consists of two groups of four individual exposures at changing wave-plate angles. Each group of four exposures is processed by Libre-Esprit into average intensity and polarization spectra. The doubling of the spectral orders when combined with the intrinsically broadened spatial point-spread function caused by the image slicer gives a peak signal-to-noise ratio (S/N) above 1500 in individual spectral pixels without saturation. Thus, a polarimetric sequence consisting of eight exposures (16 individual spectra) spatially broadened by the image slicer gives very high precision at full spectral resolution.

On each night during 2008 March 20?March 26, we obtained one full polarimetric sequence. On March 27 and 28, two complete sequences were obtained. Table 1 lists the date and time of the exposures, the number of images combined to produce the spectrum, the signal to noise for various lines in this combined spectrum, the Julian date, and finally the corresponding orbital phase. The orbital phase was computed using P and T ₀ listed in Table 2 from Herbison-Evans et al. ( 1971 ). A phase of zero corresponds to periastron passage.

ESPaDOnS is optimized for a broad wavelength coverage, having a high efficiency from roughly 4000 to 8000 A. Libre-Esprit automatically outputs an estimate of the signal to noise for each spectral order based on the average flux in each order. Figure 1 shows the wavelength and temporal dependence of the S/N for our data set. The peak signal to noise occurs in order 17 around 5150 A. The ratio is above half the peak value over the 4000?8500 A range. Table 1 shows that the precision obtained for all the lines presented in this paper are similar. The ratio only changes by 35% from 4270 to 6660 A.

Zoom In Zoom Out Reset image size

All cross-dispersed echelle spectrographs have multiple spectral orders with wavelength-dependent order overlap. Typically shorter wavelength orders overlap significantly while longer wavelength orders leave order gaps. The spectral orders must be individually continuum normalized before being spectrally merged to account for the wavelength-dependent sensitivity. The Libre-Esprit package produces both continuum-normalized and unnormalized spectra. In this package, each individual spectral order is rectified using a high-order polynomial response function to the continuum regions. However, this polynomial fit was not accurate enough for this work. The continuum normalization was done on each unnormalized individual exposure after orders were spectrally merged. This merging and continuum misfitting can be problematic for accurate spectroscopic analysis. The analysis presented in this paper avoids lines falling near order boundaries and in regions where the order overlap is incomplete. This continuum fitting provides sufficient accuracy for this work (see Appendix A and ESPaDOnS instrument description on the Web).

We have created IDL scripts to deal with the non-uniform overlapping spectral regions by binning the intensity values into essentially constant wavelength intervals of Δλ = 0.1 A. This wavelength interval is typically four times the individual spectral-pixel size. These scripts result in a monotonic and nearly regular wavelength coverage as well as increased signal to noise per wavelength bin resulting from the 2?12-pixel averaging. The properties of the ESPaDOnS spectral coverage and a detailed description of the effects of this script are contained in Appendix A . As an example, Figure 2 shows a plot of the March 20 merged and rectified spectrum where the very high signal-to-noise profiles can be seen.

Zoom In Zoom Out Reset image size

The spectrum of Spica is very similar to that described by Riddle ( 2000 ). It contains two sets of photospheric absorptions lines, each typical of an early B-type spectrum. The components arising in the primary star, m ₁ are significantly stronger than the lines arising in the secondary component, m ₂ . Figure 2 displays a fragment of the spectrum at two orbital phases close to elongations, where the two sets of lines are clearly resolved. The m ₂ components are indicated in this figure, along with the identifications of the strongest absorption lines.

The most noteworthy features of Spica's spectra are the following. (1) The width of m ₁ 's lines is consistent with the value of v sin? i = 161 km s ^?1 reported by Smith ( 1985a ), while that of m ₂ is, indeed, significantly smaller. (2) We are able to clearly identify He? ii 4685.7. Although its blue wing is blended with Si? iii 4683.80, the relatively small shift (?20 km s ^?1 ) of the centroid of the feature with respect to m ₁ 's reference frame indicates that He? ii is the dominant contributor to this blend. If both Si? iii and He? ii contributed equally, the expected shift would be ?64 km s ^?1 . We also confirm the absence of He? ii 4541.6, thus concurring with Riddle ( 2000 ) and J. Aufdenberg et al. (2009, in preparation) that m ₁ is a B0.5 star. (3) All weak lines present in m ₁ 's spectrum display the discrete narrow features, while lines having larger opacities, such as the H-Balmer series and the He? i lines do not. We do not measure any significant difference in the equivalent widths of m ₂ lines at the two opposite elongation phases illustrated in Figure 2 .

Zoom In Zoom Out Reset image size

3.1.?Radial Velocities

The first step in the line-profile analysis consists of transforming the wavelength scale to a velocity scale in the frame of reference of m ₁ , which requires knowledge of its orbital motion. There are currently at least two determinations of the stellar and orbital parameters of Spica, the first from Herbison-Evans et al. ( 1971 ) and the second from J. Aufdenberg et al. (2009, in preparation). The parameters we used are listed in Table 2 in two columns, the first derived from the results from Herbison-Evans et al. (labeled "H-E 1971") and the second derived from those of Aufdenberg et al. (labeled "Aufden 2009").

Table 2. ?Spica Parameters

Parameter	Set "HE"	Set "Aufden"
m 1 ( M ? )	10.9 ± 0.9(1)	10.25 ± 0.68(4)
m 2 ( M ? )	6.8 ± 0.7(1)	6.97 ± 0.46(4)
R 1 ( R ? )	8.1 ± 0.5(1)	7.40 ± 0.57(4) a
R 2 ( R ? )	...	3.64 ± 0.28(4) a
P (days)	4.01459(1)	4.0145898
T 0 (JD)	2440678.09(1)	2440678.09
e	0.146(1)	0.067 ± 0.014(3)
i (deg)	66 ± 2(1)	54 ± 6(4)
ω (deg) at T 0	138 ± 15	140 ± 10(4)
Aps. period (yr)	124 ± 11(1)	135 ± 15(4)
v 1 sin? i (km s ?1 )	161 ± 2(2)	161 ± 2(2)
v rot1 (km s ?1 )	176 ± 5	199 ± 5
v 2 sin? i (km s ?1 )	70 ± 5(2)	70 ± 5(2)
v rot2 (km s ?1 )	77 ± 6	87 ± 6
β 0 ( m 1 )	1.3	1.88 ± 0.19
β 0 ( m 2 )	...	1.67 ± 0.5

Notes. This table lists the parameters of the Spica binary system from Herbison-Evans et al. 1971 (H-E 1971) and J. Aufdenberg et al. 2009, in preparation (Aufden 2009). ^a These are the radii at the poles. References. (1) Herbison-Evans et al. 1971 ; (2) Smith 1985a ; (3) Riddle 2000 ; (4) J. Aufdenberg et al. 2009, in preparation.

Download table as:? ASCII Typeset image

For the analysis presented in this paper, we chose the lines of Si? iii 4552.62 and He? ii 4685.7 (+Si? iii ) for the following reasons: (1) they are relatively isolated from other neighboring lines; (2) they are relatively free from contamination by other atomic transitions; and (3) the contribution from m ₂ is clearly visible in Si? iii , but absent from He? ii . We also note that the equivalent widths of these Si? iii lines are insensitive to effective temperature.

The strong line-profile variability in Spica requires that the measurement of spectral lines be performed with great caution. We exemplify the problem in Figure 3 with the Si? iii 4552.62 A line observed on March 26. The location of the line centroid is indicated for measurements performed at different intensity levels; as the measured intensity level approaches the continuum, the centroid of the line becomes increasingly more negative. The RV measured near the minimum of the line differs by 32 km s ^?1 from the RV measured at the continuum level. This particular example is very illustrative because the contribution from m ₂ , clearly visible in the plot, cannot contribute in any manner to the shape of m ₁ 's profile. Thus, it is very difficult at this stage to decide what portion of the line profile most reliably represents orbital motion and this situation introduces a high degree of uncertainty in the determination of the stellar and orbital parameters.

Given this uncertainty, we measured Si? iii and He? ii 4685.7 A at as many intensity levels as possible, without reaching levels where m ₂ 's absorption is evident. The average value of these measurements is listed in Table 3 , together with an uncertainty that corresponds to 1/2 of the total range in RV values measured for the line. The resulting RVs are plotted in Figure 4 . Also plotted in this figure are the RV curves corresponding to the J. Aufdenberg et al. (2009, in preparation) parameters (solid) for three different values of ω, the argument of periastron, which defines the location of line of apsides with respect to the plane of the sky. J. Aufdenberg et al. (2009, in preparation) give ω = 140° ± 10° at T ₀ = JD 2440678.09, as well as the apsidal period P _aps = 135 ± 15 years. Thus, for the current observations, ω ₂₀₀₈ lies between 248° and 257°. The measured RVs indicate that ω ₂₀₀₈ = 255°, as illustrated in Figure 4 . Since this lies within the range of acceptable values, we adopt ω ₂₀₀₈ = 255° for this study.

Zoom In Zoom Out Reset image size

Table 3. ?Radial Velocities

Day	JD	Phase	Si? iii		He? ii		He? i	m 2	Na? i
Mar			RV	Err	RV	Err	RV	EW	ISM
20	546.022	0.377	94	3	64	4	?115	0.086	?7
21	547.100	0.652	?77	14	?91	1	...	...	?8
22	548.031	0.878	?87	14	?118	10	179	0.11	?8
23	549.054	0.132	73	6	44	3	?135	0.12	?9
24	550.098	0.393	84	3	43	2	...	...	?10
25	551.115	0.652	?79	9	?96	1	...	...	?8
26	551.899	0.841	?106	16	?141	8	193	0.10	?10
27	552.980	0.110	68	3	39	4	...	...	?9
28	553.970	0.357	103	2	71	2	?139	0.11	?8

Notes. This table shows the radial velocities (RVs) measured for each observation with respect to laboratory wavelengths. The columns show the 2008 March observation date, Julian date (JD 2454000), orbital phase of observation, RV and uncertainty (err) of the primary star for Si? iii (4552.2 A) and He? ii (4685.7 A) lines, the secondary companion He? i (6678.15 A) RV, and equivalent width in angstroms and finally the Na? i D ISM line RV.

Download table as:? ASCII Typeset image

3.2.?Line-profile Variability

Line-profile variability in Spica was first analyzed in detail by Smith ( 1985a , 1985b ) in a closely spaced set of data obtained over the complete orbital cycle. He described the variability in m ₁ in terms of the presence of "bumps" that deform the shape of the photospheric absorption line profile from that which it would have if its surface were rotating as a rigid body. The individual "bumps" were found to travel across the line profile, from blue to red, at rates of ∼12?17 km s ^?1 per hour. The analysis of the pattern of moving bumps seen in the Si? iii lines at 4552, 4567, and 4574 A led him to conclude that they could be accounted for by four NRP modes. Two of the modes had periods of ∼1/12 P and 1/2 P , the latter being associated to "the spectroscopic equivalent of the ellipsoidal light variability" which is due to the tidal distortion.

Our limited data set is unable to fully address the short-timescale variability since we only have a single spectroscopic observation set of eight exposures taken over less than 5 minutes on most nights. The two nights with 16 exposures (27th and 28th) only cover 10 minutes. However, we do have two pairs of nights (March 27 and 23, 25 and 21) where there are small orbital-phase differences between observations corresponding to a few hours of stellar rotation. These pairs allow the characterization of changes in the "bump" structure over short timescales, though with one complete orbit actually in between the observations. The night-to-night variability can be easily addressed with such high-quality data.

Figure 5 illustrates the Si? iii 4552.62 A line in the nine different orbital phases covered by our observations, plotted on a velocity scale after applying the correction for orbital motion as listed in Table 3 . The variability in this line is typical of that seen in all other weak lines. Note that the absorption from m ₂ can be clearly followed in these spectra. In the first five profiles (? = 0.1?0.4, from bottom to top), the absorption line originating in the secondary component, m ₂ , is observed to move toward maximum negative velocity and then return toward the rest-frame velocity. In the last two profiles (? = 0.8472, 0.8838), m ₂ 's contribution may be seen near its maximum positive velocities. In the two profiles at ? = 0.652, the contributions from m ₁ and m ₂ overlap.

Zoom In Zoom Out Reset image size

We have marked in Figure 5 the location of the discrete absorptions that are resolved in our data and whose RV variations could be measured. In general, at least five such troughs may be identified, although at times some of these may blend together leading to the appearance of only three troughs (for example, at orbital phases 0.3623, 0.3776, and 0.3987). If we assume that the troughs retain their identity between two consecutive orbital phases, ⁵ then we find that between ? = 0.1167 and 0.1387, the rate by which they travel from blue to red is Δ V _bump /Δ t = 22 km s ^?1 per hour, and between ? = 0.8472 and 0.8838, Δ V _bump /Δ t = 8 km s ^?1 per hour. For ? = 0.35?0.39, one of the bumps appears to remain stationary, the second moves slightly "redward," and the third bump seems to split into at least two absorption-like features.

Thus, the traveling bumps described by Smith ( 1985b ) are clearly visible in Figure 5 . However, only the troughs for ? ∼ 0.1 in Figure 5 travel at an average displacement speed similar to the ∼17 km s ^?1 per hour speed that may be inferred from Smith ( 1985a ; Figures 1 ? 4 ). At the other orbital phases covered by our data, the velocity is significantly slower. These speeds, of course, are the result of the actual motions on the stellar surface combined with the projection of the stellar rotation velocity along the line of sight to the observer. Further insight into the phenomenon requires the use of a theoretical model.

The stability on ∼2 hr timescales may be fully appreciated in Figure 6 . In three of the four comparisons of Figure 6 , there is excellent agreement and the identity of individual bumps is almost entirely preserved. This is despite the fact that the observations are separated by a full orbital cycle and that diverse variability sources could possibly be present. The most striking difference appears in the third comparison, corresponding to the line profiles from March 21 and 25 obtained at nearly the same orbital phase, 0.6515 and 0.6523. These observations have the closest phase overlap of all four shown in the figure and yet they present the most significant line-profile differences. This indicates that the variability is not fully orbital-phase locked.

Zoom In Zoom Out Reset image size

Figure 6 also provides a view of the manner in which the absorption line wings rise from the core to the continuum level. There is a marked difference between the shape of the line wings at different orbital phases. For example, at orbital phases ∼0.3?0.4, the red wing gives the profile a more "boxy" shape than at other phases, changes that may be quantified using the kurtosis , which describes the degree to which a distribution is centrally peaked or flat. The bottom panel of Figure 4 shows the phase-dependent behavior of the kurtosis, confirming the difference that is apparent from Figure 6 in the line shapes. The kurtosis was computed over a wavelength range that avoids the contribution from m ₂ , and the error bars correspond to an uncertainty ±0.2 A on the definition of each line wing.

In summary, two categories of line-profile variability are present in our data: (1) there are traveling bumps that retain their identity over timescales of at least a day, though some may change significantly and (2) there is a varying shape of the line wings which go from being very extended to "box like" on a nightly timescale. We will show below that both of these characteristics are a consequence of surface flows that can be produced by the tidal interaction between the two components.

A system is in synchronous rotation when the orbital angular velocity Ω equals the angular velocity of rotation ω ₀ . In eccentric orbits, the degree of synchronicity varies with orbital phase. We use periastron passage as the reference point for defining the synchronicity parameter, β _per . ⁶ When at any orbital phase β ≠ 1, non-radial oscillations are excited, driven by the tidal interactions. We refer to the azimuthal components of the forced oscillations as "tidal flows."

In the following paragraphs, we summarize the general characteristics of a simple model developed to gain some insight into the effects that the tidal flows introduce in the photospheric line profiles. A more detailed description is provided in Moreno & Koenigsberger ( 1999 ) and Moreno et al. ( 2005 ).

Consider a binary system in which the light of the primary star with mass m ₁ and radius R ₁ dominates the shape and intensity of a photospheric spectral line. The main body of the star, below the surface layer of thickness Δ R is assumed to have a constant spin angular velocity ω ₀ with respect to an inertial frame. The binary companion has mass m ₂ and an instantaneous orbital angular velocity Ω in its (in general) elliptic orbit around m ₁ . The axis of stellar rotation is assumed to be perpendicular to the orbital plane. We define β ₀ = ω ₀ /Ω ₀ , where Ω ₀ is the orbital angular velocity at periastron. β ₀ specifies whether the stellar rotation and the orbital periods are synchronized (β ₀ = 1) or whether the rotation is super- or subsynchronous (β ₀ > 1 or β ₀ < 1, respectively). If β ₀ ≠ 1.000, then tidal flows on the stellar surface are induced and non-radial oscillations are excited. Our model is based on the fact that the surface layer is the one most strongly affected by the external gravitational potential of the companion star (see, for example, Dolginov & Smel'Chakova 1992 ; Eggleton et al. 1998 ) and we therefore compute the response of this layer to all the forces in the system, while the rest of the star below it is treated as a rigidly rotating body. We find a posteriori that the dominant contribution to the photospheric line-profile variability is produced by the azimuthal component of the velocity field and not the radial component. In our model, the azimuthal flow structure is governed by the tidal forcing and appears to be very weakly dependent on the radial component of the velocity field.

A non-uniform Lagrangian grid is constructed over the stellar surface by specifying the number of surface elements along the equatorial belt and the number of latitudes for the computation. The number of surface elements for each latitude is defined so that in the unperturbed star all elements have similar sizes in the longitudinal coordinate. For each surface element in the grid, an equation of motion is specified in which the following forces are included: the gravitational field of m ₁ and m ₂ ; gas pressure; Coriolis force; centrifugal force and the viscous shear forces associated with the interface of the thin surface layer with the inner stellar body; and the interface between two adjacent surface elements. The surface layer may be treated either in the polytropic or in the isothermal approximation. The cases described below correspond to calculations with polytropic indices of 1.5. Cases run with a different polytropic index or in the isothermal approximation yield similar global characteristics. The simultaneous solution of the equations of motion for all surface elements yields values of the radial, v _r , and azimuthal, v _φ , velocity fields over the stellar surface. Azimuth angle, φ, is measured from the line joining the centers of m ₁ and m ₂ , in the non-inertial reference frame centered on m ₁ , and rotating with the instantaneous orbital angular velocity Ω. Polar angle, θ, is measured from the pole of m ₁ to its equator.

The model is not restricted to slow rotation nor to small eccentricities. The only a priori assumptions built into the calculation are (1) only the thin outer layer responds to the tidal perturbation, (2) the stellar rotation axis is perpendicular to the orbital plane, and (3) the response of the outer layer is either polytropic or isothermal. It must be kept in mind, however, that because of the single-layer approximation, we compute only the motion of the center of mass of each surface element, thus ignoring the details of the inner motions within these elements, and thus the buoyancy force is neglected.

4.1.?Input Parameters

4.2.?The Azimuthal Velocity Field

We use the term "tidal flows" to refer to the large-scale azimuthal velocity patterns that are induced by the gravitational field of the binary companion. An intuitively simple way to grasp the essence of the idea of tidal flows is given by the description of Tassoul ( 1987 ): "Assume that at some initial instant the primary is rotating rigidly with a constant supersynchronous angular velocity. If there were no companion, the system would remain forever axisymmetric. Because of the presence of the secondary, however, a fluid particle moving on the free surface of the primary will experience small accelerations and decelerations in the azimuthal direction." In the simplest of cases, the tidal-flow pattern may be characterized in terms of four large contiguous zones over which the acceleration alternates between positive and negative values. This defines four quadrants over the stellar surface, the first of which always lies near the subbinary longitude, and in which the acceleration is negative for supersynchronous rotation. Figure 7 is an example of such velocity perturbations obtained from the analytical solution given by Scharlemann ( 1981 ) for a binary system with m ₁ = 10.9 M _? , m ₂ = 6.8 M _? , R ₁ = 8.1? R _? , and P = 4.01?days. The figure shows a plot of the azimuthal velocity versus azimuth angle for different latitudes, starting near the equator. Because his formulation is valid only for circular orbits, near-synchronous rotation and inviscid fluids, we used β = 1.001 and ν = 0.04 R ² _? day ^?1 , the minimum allowed viscosity by the code for this case. The dashed curves in this figure illustrate the results from the calculations of our code for a Δ R / R ₁ = 0.025 surface layer and the same orbital and stellar parameters. The one-layer velocity perturbations are very similar in shape and in amplitude to those of Scharleman's model at latitudes lower than ∼44°. At intermediate latitudes (∼44°?70°) the shape is still the same, although the one-layer amplitudes are systematically larger.

Zoom In Zoom Out Reset image size

The velocity field becomes significantly more complex when eccentricity and faster stellar rotation rates are introduced, as illustrated in Figure 8 . There we show the azimuthal velocity field for the GGD6 one-layer model calculations for one of the parameter sets similar to those of Spica. Note that the pattern still consists of two maxima and minima, but their shape can no longer be described as simply sinusoidal and, in addition, a great deal of substructure is apparent. Furthermore, because the orbit is eccentric, the pattern changes from one orbital phase to another. The maxima appear around azimuth angle ∼90°?150° and ∼270°?330°, and correspond to material that is flowing in the direction of the stellar rotation. The minima correspond to surface material that, with respect to the rigidly rotating inner region, is flowing in the opposite direction; with respect to the external observer, this material is flowing slower than its corresponding rotation speed.

Zoom In Zoom Out Reset image size

These "tidal flows" are represented for our model run GGD in Figure 9 (top left) with light and dark colors, corresponding to positive and negative perturbation speeds, respectively. This figure also includes maps of the radial perturbations in both velocity, v _r (middle), and in the stellar radius, R _φ (right).

Zoom In Zoom Out Reset image size

It is important to note that the "tidal bulge" is usually characterized as a region over which the maximum radial extent of the star is attained. Our model suggests that, although on average the tidal bulge region is more extended, the detailed structure of the tidal bulge for an eccentric, supersynchronously rotating star such as Spica consists of small regions of alternating large and small radial extent. These are accompanied by similar oscillations in the velocity perturbations. Hence, we conclude that even for a relatively small eccentricity, the "dynamical tide" plays a major role and cannot be neglected. In Figure 10 , we illustrate the azimuthal velocity field evolving over an orbit on the star's hemisphere that would be visible to a distant observer. The velocity fields are shown for four different orbital phases: ? = 0 (periastron), 0.25, 0.5 (apastron), and 0.75, and assuming a longitude of periastron of ω = 255°.

Zoom In Zoom Out Reset image size

A second important point to note is that the pattern of flow velocities is relatively fixed with respect to the coordinate axis that rotates with the binary system; it is not fixed with respect to the coordinate axis that is rotating with m ₁ . Thus, this is conceptually very different from the idea of "spots" that have relatively fixed locations on the stellar surface and rotate in and out of our field of view with the rotation angular velocity. Rather, the actual pattern changes with the orbital period; but given that there are two maxima and two minima per cycle, this leads to variability on timescales 1/2 P .

4.3.?Line-profile Variability

Once the surface velocity perturbations have been computed, the projection of the velocity field along the line of sight to the observer provides the Doppler shift needed to produce the rotationally broadened and tidally perturbed line profile that may be compared with observations. The equations needed to perform the projection along the line of sight to the observer are described in Moreno et al. ( 2005 ). The Doppler shifts derived from these projections are then applied to the "local" line profiles, i.e., the shape of the spectral line that is assumed to correspond to the emergent line radiation from each surface element.

Ideally, the local line profile should be derived from a model atmospheres computation. In our current calculation, we adopt a Gaussian local line profile for the following reasons.

First, we have investigated the effect of various local line-profile shapes on the shape of the total perturbed integrated profile and found the variation to be a very small. We used CMFGEN model atmosphere profiles (Hillier & Miller 1998 ) for T _eff = 25,000 K and log( g ) = 3.75 (kindly provided by John Hillier) as well as Gaussian and Voigt profile fits to these CMFGEN profiles. The CMFGEN model provides emergent flux (local line profiles) for a set of locations across the stellar disk (μ). These model atmosphere profiles are then fitted with Gaussian and Voigt functions at every μ so an accurate integrated line profile can be made using different local line-profile shapes at every μ. For the Voigt function fits, both parameters are varied to minimize the deviation between Voigt and CMFGEN profiles. Typical values for the damping parameter were 0.27?0.48. This illustrates how the variation in local line-profile shape influences the total perturbed line profile and is detailed in Appendices B ? D . All local profiles give essentially identical rotationally broadened line widths. The width of the total line profile is rotational velocity dominated and the intrinsic line width and shape play little role. The local line-profile shape is somewhat more significant when computing the detailed shape of the discrete absorption-like features. A wider local profile will tend to smooth out these bumps, increasing their width and decreasing their amplitude, but the bump velocity and number of identifiable bumps remain unchanged for all the cases we examined. It is interesting to note that the CMFGEN model spectra produce slightly narrower local line profiles for the Si? iii line (smaller wings, shallower core) than do the Gaussian or Voigt profiles and would slightly increase the relative contrast of the bumps. We have investigated this effect using our computed tidal velocity fields as well as simulated velocity fields. These local line-profile effects were found to be very small and can be easily neglected for the velocity fields and perturbed line profiles calculated here. See the Appendix for more details.

Second, it is important to note that given the dynamical conditions on the stellar surface, non-uniform temperature, and microthermal speed distributions are most likely present. In addition, the shearing flows may lead to radial temperature gradients that differ from point to point over the stellar surface. Without a full three-dimensional radiative transfer treatment, the impact of these effects on the shape of the local line profiles is difficult to predict and constitutes a source of systematic error. Thus, adopting local line profiles from a one-dimensional model atmosphere is subject to a wide range of systematic errors and is as much of an approximation as is the use of a Gaussian local line profile.

Table 4 summarizes the range in input parameters that was explored in the nearly 200 model calculations run for Spica as well as some specific examples. A general feature of all the calculations is that the line profile undergoes a strong change in its shape on a day-to-day timescale, transitioning from a Gaussian-like profile to a more "boxy" shape. These changes may be associated with the "equilibrium tide" component of the gravitational perturbation. As pointed out by Smith ( 1985a ), it is the spectroscopic analog of the ellipsoidal variability in the light curves and produces a period of 1/2 the orbital period.

The change in the line profiles associated with this effect can be quantified using the kurtosis. The bottom panel in Figure 4 shows the variation of the kurtosis in the computed line profiles for the model Case GGD over the orbital cycle. As expected from the above discussion, there are two maxima and two minima in this curve, corresponding to the transitions into and out of the "boxy" line-profile shape. In order to compare the theoretical kurtosis with the observational data, we also measured the kurtosis in the Si? iii 4552.62 A line of our nine spectra. Here, however, the wavelength interval measured was chosen to avoid the contribution from the binary companion. The trends in these measurements are plotted in Figure 4 as the filled-in triangles. The error bars in the data correspond to the uncertainty introduced by increasing or decreasing the wavelength interval by 0.1 A. Thus, we conclude that our model adequately reproduces the variability trends in the observational line profiles showing that tidal flows are an important source of non-RVs.

It is interesting to note that the line-profile variability introduces a systematic error in the measurement of the line centroid, as illustrated in the middle panel of Figure 4 . There we plot the difference between the centroid measured on the perturbed line profile and the unperturbed profile at the same orbital phase. For the cases considered here, the largest error occurs in the orbital-phase interval 0.6?0.7.

A comparison of the observations with examples of the predicted line profiles is shown in Figure 11 , illustrating that the model reproduces the orbital-phase-dependent variations in both the general shape and in the number and approximate location of the discrete absorptions. It also reproduces the "blue-to-red" motion of the discrete features as illustrated in the montage of all the theoretical line profiles illustrated in Appendix B . Also interesting to note is the presence in our computed profiles of the "red spike" found by Smith ( 1985b ) and attributed to the horizontal motions. Finally, we note that since the profiles for all orbital phases are computed sequentially during the same run, the comparison of their time-dependent changes (from one orbital phase to the next) with that of the data constitutes yet another constraint. Figure 11 illustrates that the model calculation is consistent with the data also in this respect.

Zoom In Zoom Out Reset image size

The one important inconsistency between the model line-profile calculation shown in Figure 11 and the data is that we used ω = 209° for the input parameter corresponding to the argument of periastron, instead of the value ω = 255° implied by the RV curve. We wish to stress, however, that the theoretical profiles presented in Figure 11 are not the product of a fitting process; these example profiles result directly from the tidal interaction computation using one of the sets of stellar and orbital parameters that have been derived for the Spica system. Small modifications in some of the input parameters can have a significant impact on the details of the line profiles. A more in-depth discussion of the effects on the line profiles due to different input parameters will be presented elsewhere. However, it is important here to illustrate the effects produced by the two free parameters of the model, ν and dR / R ₁ . Figure 12 shows that the large-scale orbital-phase-dependent variability is generally much stronger than the variations produced by changing ν and dR / R ₁ . Hence, once the stellar and the binary parameters are constrained, the number of free parameters that enter into our calculation are greatly reduced. Finally, Figure 13 illustrates the perturbed line profiles computed with the full velocity field (as in Figures 11 and 12 ) compared with the profiles computed with only the horizontal velocity component, Δ v _φ (dashes). These profiles are nearly identical, showing that the radial component of the velocity field contributes very little toward the total line-profile variability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we have presented very high quality observations of line-profile variability. The combination of high spectral resolution with very high precision allows a much more detailed comparison of observation and theory. Though we do not have a continuous sequence of spectra, we presented some observations that show remarkable similarities after a full orbit while others show significant differences between orbits.

We find that the Moreno et al. ( 2005 ) approach to computing line-profile variability from tidal flows yields theoretical profiles that can be very similar to observational data. Specifically, the theoretical profiles display two characteristics of Spica's line-profile variability: (1) the "boxy" appearance of the line that is observed twice per orbital cycle and that is associated with the "equilibrium tide" and (2) the superposed narrow absorption-like features that travel from the blue wing to the red wing on timescales of ∼1 day and that may be associated with the "dynamical tides." This is also the "spike" reported in Smith ( 1985b ). Both of these characteristics are caused primarily by the horizontal components of the surface velocity field and thus, we conclude that tidal flows are the dominant contributor to Spica's line profile. Because the simple one-layer model cannot be expected to yield information on the internal structure of the star, the relatively good coincidence between the ab initio calculation and the observations indicates that in a close binary system such as Spica, the tidal forcing by the companion is the dominant factor determining the surface perturbations. The azimuthal velocity fields that are calculated from tidal interactions have the potential to explain much of the observed variability.

Given the fact that our computation is performed from first principles with no built-in assumptions regarding the surface velocity field and heavy physical constraints on our input parameters, we consider the coincidence between the theoretical and observational line profiles to be very encouraging. We do use a single-layer fluid approximation but all the dynamics are controlled from the fundamental equations of motion. The structure of the fluid flow is not constrained to fit any specific mathematical form. The velocity field, subject to the single-layer approximation, is computed directly from the equations of motion. This illustrates the importance of properly treating surface flows (the horizontal component of the surface motions) when analyzing observed line-profile variability.

There are several aspects of this model that can be readily explored. The strength of the narrow spectral features (bumps) are not ideally reproduced. The feedback from shear energy dissipation into surface temperature and local line-profile shape should be included to test the significance of dissipation in computing the predicted profiles. Ideally, model atmosphere line profiles and limb darkening should be used in the general scheme, although a proper treatment requires the use of a three-dimensional model atmosphere calculation. In the meantime, however, we have shown the shape of the local line profiles to be of little significance for this particular system. Also, incorporating a more realistic stellar structure into the calculation should allow exploration of the effects produced by the interplay between the tidal forcing and the normal modes of oscillation.

We express our gratitude to the referee, Myron Smith, for comments and suggestions that helped improve the paper significantly, to John Hillier for providing the CMFGEN model spectra, Ben Brown for guidance in producing maps, and Jason Aufdenberg for providing preliminary results on Spica's fundamental parameters. This program was partially supported by the NSF AST-0123390 grant, the University of Hawaii and the AirForce Research Labs (AFRL), as well as UNAM/PAPIIT IN-106708 and CONACYT 48929. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. These observations were reduced with the dedicated software package Libre-Esprit made available by J. -F. Donati. This program used the Simbad data base operated by CDS, Strasbourg, France as well as the Markwardt and Coyote IDL Libraries available on the web. The NASA Astrophysics Data System (ADS) was frequently utilized.

The nonlinear and irregular wavelength sampling of the ESPaDOnS spectropolarimeter led us to develop routines to linearize the output spectra. Libre-Esprit produces wavelength-calibrated spectra for every possible spectral pixel in the focal plane. The spectral orders are listed in Table 5 . As with all cross-dispersed echelle spectrographs, the dispersion, order overlap and sampling are all functions of wavelength. An additional complication is that the ccd pixels correspond to 2.6 km s ^?1 but the spectra are output by Libre-Esprit at 1.8 km s ^?1 resolution, 30% smaller than a pixel. These are called ccd and spectral pixels, respectively. Figure 14 (a) shows the order overlap as a function of wavelength for a typical ESPaDOnS exposure in our data set. The order overlap increases from in overlap from 20 A at 4000 A to over 60 A around 7000 A before falling strongly to become a more than 40 A gap in wavelength coverage by 10000 A The dispersion, shown in Figure 14 (b), runs from 0.02 A per pixel to 0.06 A per pixel in proportion to the wavelength increase. This leads to an essentially constant spectral resolution, shown in Figure 14 (c) as roughly R ∼ 170,000. As the full width at half-maximum of most calibration-lamp lines is over 2 pixels, the spectral resolution is around 70,000.

Zoom In Zoom Out Reset image size

We desired a monotonic wavelength grid that assigns one intensity to one wavelength and combines all pixels covering the same wavelength to achieve higher S/N. The processing scripts we wrote calculate the number of spectral pixels to bin to achieve 0.12 A per pixel coverage for all wavelengths. This is shown in Figure 14 (d). Since the order overlap is significant, in some cases 30% of an order, there is a bi-modal distribution in the number of pixels binned. For the part of an order with overlapping coverage, a linear wavelength coverage gives roughly twice the number of pixels binned. Beyond roughly 8000 A, the order overlap disappears and the dispersion reaches 0.05 A per pixel resulting in a nearly regular 2-pixel binning factor.

After applying the rebinning routines, the dispersion remains constant at 0.12 A per pixel as seen in Figure 15 (a). This results in a spectral resolution that is wavelength dependent, going from 30,000 to 80,000 per spectral-pixel at 1 pixel per resolution element. This is shown in Figure 15 (b). As stated above, the physical spectral resolution of the raw output is roughly 70,000 with over two spectral pixels per resolution element. This new binned output undersamples 2:1 at 4000 A increasing to 1:1 sampling at 7000 A and rising with longer wavelengths. The signal to noise will be increased roughly as the square-root of the number of pixels averaged as all our observations are in the high signal-to-noise regime. The signal-to-noise improvement factors are shown in Figure 15 (c) for both order overlap and no overlap. Libre-Esprit outputs an estimate of the signal to noise for each order as part of the reduction routine. An example of this is seen in the lower curve of Figure 15 (d). After the regularization of the wavelength array, the signal to noise is increased. Using the modest improvement with no order overlap, the resulting signal to noise is shown as the upper curve in Figure 15 (d). Though these routines result in a wavelength-dependent resolution (with constant dispersion), the regular, linear, monotonic wavelength sampling allows for easier processing of large numbers of observations along with generally higher signal to noise.

Zoom In Zoom Out Reset image size

We also performed some experiments on the accuracy of the signal-to-noise predictions given by Libre-Esprit. In the very high signal-to-noise regime, systematic effects caused by flat-fielding issues, nonlinearity or instrument instabilities can become significant. The Libre-Esprit reduction logs output two estimates of the signal to noise based on the flux. There is an estimate for each "spectral pixel" and "ccd pixel." The difference between these pixels is the resolution of 2.6 km s ^?1 per physical ccd pixel compared to 1.8 km s ^?1 per spectral pixel as resampled in Libre-Esprit.

The noise properties of the ESPaDOnS data in the high signal-to-noise regime can be best seen in our best observation set: March 27. These properties are measured directly on the Libre-Esprit output data without applying our wavelength-regularizing scripts. Figure 16 shows four separate spectra of Spica for order number 44 centered at 5150 A. There are two complete polarization measurements which each consist of two intensity spectra output by Libre-Esprit. Each individual spectrum there is calculated from four individual exposures each with 2 orders. Thus, each spectrum in Figure 16 (a) represents an average of eight individual spectra. The Libre-Esprit package gives S/N for each plotted spectrum of roughly 1350/1600 per spectral and ccd pixel. The wavelength solution is derived independently for each four-exposure polarization measurement. The shift between the wavelength solutions is typically 0?2 pixels. In Figure 16 (b), the change in wavelength solution is shown to be a small fraction of a pixel across this spectral order. The illumination pattern on the CCD changes only slightly between individual exposure sets.

Zoom In Zoom Out Reset image size

The noise in the spectra can be estimated by subtracting a boxcar-smoothed spectrum from a raw spectrum. As the full width at half-maximum of spectral features is roughly 2.2 pixels, any wavelength variation present at the one-spectral-pixel level is surely noise. Figure 16 (c) shows residual noise from each spectrum after subtracting a 3-pixel boxcar smooth. This highlights the noise and removes spectral structure. This noise spectrum can then be cross-correlated to determine the level of any systematic errors present. Figure 16 (d) shows the cross-correlation between each of the four noise spectra and the first noise spectra (the top spectrum in box c). The top curve represents an autocorrelation and a clear signature from the boxcar smooth is seen. The lower three cross-correlation curves also show this same signature, but spectrally offset by 1 or 2 pixels and at significantly lower amplitude. The wavelength solutions differ by less than 10% of a pixel across an order and are generally stable. However, a 1?2 pixel shift between exposures has been recorded by Libre-Esprit between the four wavelength solutions. This correlation between the noise spectra is likely caused by structure on the CCD.

In this appendix, we will describe some general properties of the model calculations. Figure 17 shows a typical model run. The code is allowed to run for roughly 30 orbital cycles to damp out any perturbations from initial conditions. The perturbed velocity field is then calculated at each time step, as is its projection along the line of sight to the observer, from which a total line profile is output. The right-hand panel of Figure 17 shows the average line profile over the orbital cycle as the dark solid line. Several individual spectra from throughout the orbit are also overplotted as the dashed and dotted lines. Once the average line profile is subtracted from each individual measurement, a grayscale plot of the deviations can be constructed to show how individual surface perturbations move across a line profile and preserve their identity over an orbital cycle. The left-hand panel is dominated by the main tidal bulge?the light and dark patches that are most prominent on the line wings (largest velocities). However, the smaller-scale surface perturbations can be more easily seen near line center as the individual diagonal streaks.

Zoom In Zoom Out Reset image size

In order to test the accuracy and utility of using Gaussian local line profiles, we performed an investigation as to the effect of using different shape local line profiles in computing the integrated rotationally broadened, tidally perturbed profiles. We obtained CMFGEN model atmosphere profiles with T _eff = 25, 000 K and log( g ) = 3.75 provided by John Hillier. Since the CMFGEN model provides emergent flux (local line profiles) for a set of projected locations on the stellar disk (μ) this database of profiles represents our starting point to simulate local line-profile effects. We now have a model local line-profile shape for any point on the stellar surface.

These model atmosphere profiles are then fit with both Gaussian and Voigt functions at every μ. These fits are straightforward using IDL's curve-fitting routines: Gaussfit and Curvefit. Once we have CMFGEN, Gaussian and Voigt local line profiles for every μ an accurate integrated line profile can be made to simulate an observation of a tidally perturbed star. The projected velocity is known for every surface element as is the value for μ and emergent flux at every surface element. An example of these local line profiles is shown in Figure 18 .

Zoom In Zoom Out Reset image size

To illustrate the effects of local line-profile shapes on the perturbations caused by macroscopic surface flows, we computed a perturbed line profile for a velocity field using all three of these local line-profile shapes. A simulated perturbing velocity field is created as a combination of spherical harmonic functions and random fluctuations with a maximum amplitude of ±15 km s ^?1 . This creates a highly perturbed line profile with known properties that can clearly demonstrate the effect of local line-profile shapes. Figure 19 shows a perturbed velocity field in the top panel as well as integrated line profiles in the bottom panel. The integrated line profiles are calculated using the perturbed stellar velocity field. The CMFGEN (solid), Gaussian (dot-dashed) and Voigt (dash) local Si? iii line profiles are used to generate an integrated profiles seen in the bottom panel. There is essentially no difference between integrated line profiles even though significant differences in local line-profile shapes are apparent.

Zoom In Zoom Out Reset image size

As a third test of local line-profile shape effects on integrated line profiles, we adapted our tidal-flow code to compute integrated line profiles using a more boxy local line-profile shape $-\sim e^{-x^4}$ instead of $\sim e^{-x^2}$. As expected, the difference in local line-profile shape was negligible. Figure 20 shows integrated line profiles computed for Case GGD3 parameters. The more boxy local line profiles are used for the solid curve and the Gaussian local line profiles are used for the dotted curve. There is a very small difference in the depth of the perturbations but the overall morphology is very similar.

Zoom In Zoom Out Reset image size

We consider the use of Gaussian local line profiles to be entirely sufficient to demonstrate this tidal-flow model for integrated line profiles. The integrated line-profile intensities differ at the 0.05% level for various local line-profile shapes. This is considered to be entirely below any systematic errors or errors caused by uncertainties in the physical parameters used in calculating the properties of the binary systems.

We feel it is worth pointing out that the amplitude and size of the "bumps" fairly tightly constrain the size and amplitude of the corresponding velocity fluctuations on the stellar surface. In order to produce the bumps we observe, one simply needs a velocity perturbation of order 10?30 km s ^?1 with a few "patches" of velocity perturbation on the stellar surface for each observed bump. Figure 21 illustrates this idea.

Zoom In Zoom Out Reset image size

We created radially projected velocity perturbations on the stellar surface with a random number generator scaled to a peak amplitude of ±15 km s ^?1 . The surface of the star was perturbed with this random field scaled to have an increasing number of "granules" on the stellar surface. This illustrates the effect of decreasing the size scale of velocity perturbations on the integrated line profiles. Our experiments were run with spatial scales of 0, 2, 5, 10, 20, 50, 100, 200, and 400 individual perturbations across an equatorial slice.

The top panel of Figure 21 shows first the rotational broadening of 160 km s ^?1 in the top left followed by these perturbed velocity fields. Once these perturbations are incorporated with the rotational broadening, the corresponding line profiles are calculated. The bottom panel of Figure 21 shows these perturbed line profiles. A fairly simple trend can be seen. A small number of perturbed surface patches (2?5) results in small bulk-shape deviations of the line profile. When the number of surface perturbation patches reaches 10?20, line profiles resembling those observed are seen. As the number of surface perturbations increases, more bumps appear at increasingly smaller amplitudes until (by 400 equatorial patches) there is essentially no visually apparent difference between perturbed and unperturbed line profiles. This example illustrates a simple constraint on the perturbed velocity field coming directly from the observed amplitude and number of bumps?relatively large-scale surface perturbations (a few to several patches in an equatorial slice) with an amplitude that is a fraction of the rotational velocity can cause the observed bumps. Since the structure of the observed perturbations (the line-profile shape) is more or less similar from one orbit to the next, the corresponding velocity perturbations are clearly not random.