format: text

text:
# 1 Introduction

## 1 Introduction

Optical emission lines from pre-main-sequence stars are observed as clues to understanding the acceleration of jets and mass accretion on to the stars. The spatial resolution of most techniques is insufficient to resolve their dynamics directly, while high-resolution spectroscopy has provided a variety of information for these mechanisms. H$`\alpha `$ emission shows a rich variety of profiles due to an interplay between emission and absorption features (e.g., Reipurth, Pedrosa, & Lago 1996), and various models have been discussed to explain their profiles using accelerating stellar wind (e.g., Decampli 1981), decelerating stellar wind (Mitskevich, Natta, & Grinin 1993), disc wind (Calvet, Hartmann, & Hewett 1992), and infalling envelope (e.g., Calvet & Hartmann 1992). The forbidden lines often have two blue-shifted velocity components, which have provided constraints on geometry and mechanism of jet and/or wind (Eislöffel et al. 2000, and references therein). The lack of red-shifted components of these line supports the existence of circumstellar discs (Eislöffel et al. 2000, and references therein). In addition to high-resolution spectroscopy, some techniques have been tried in the last decade to obtain spatial information. Speckle observations reveal elongated structure of the H$`\alpha `$ emission with an angular scale of 0.09 arcsec towards T Tau (Devaney et al. 1995). The $`Hubble`$ $`Space`$ $`Telescope`$ has provided fine structure of jets towards young stellar objects (Ray et al. 1996; Burrows et al. 1996) although the bandwidths of its filters are often insufficient to suppress the stellar continuum from pre-main-sequence stars. Ground-based long-slit spectroscopy is effective in measuring the typical angular scale of a few forbidden lines suppressing the contamination of the strong continuum emission (e.g., Hirth, Mundt, and Solf 1997).

Here we present the results of spectro-astrometric observations of RU Lupi. The spectro-astrometry technique studies the small scale structure by observing the relative position of the centroid of the point-spread function on scales of a few milli-arcsec (Bailey 1998a, b). The target, RU Lupi, is one of the most active T-Tauri stars in the southern sky. The H$`\alpha `$ emission has broad wings up to 900 km s-1 (Reipurth, Pedrosa, & Lago 1996), and the equivalent width of the emission line is up to 200 Å, which is one of the largest among T-Tauri star (Giovannelli et al. 1995). Forbidden lines of \[SII\] and \[OI\] are also present in the spectrum of RU Lupi, and at least the \[OI\] lines have an extended blue wing up to 300 km s-1 and a redder peak up to 30 km s-1 (Gahm, Lago, & Penston 1981; Hamann & Persson 1992). In this paper, the detail of the observations and the data reduction are described in §2, the results are shown in §3, and the origin of H$`\alpha `$ emission and the forbidden lines are discussed in §4 together with the existence of a circumstellar disc with a gap. Throughout this paper, the distance of 140 pc to the Lupus cloud is adopted based on Hughes, Hartigan, & Clampitt (1993).

## 2 Observations & Data Reduction

Observations were carried out on 1996 Aug 25, 1997 Jun 27, and 1999 Jul 2 at the 3.9-m Anglo-Australian Telescope using the RGO spectrograph. The configuration with a 1-arcsec slit width, a 1200 line mm-1 grating, and the 82-cm camera provided a spectral resolution ($`\lambda /\mathrm{\Delta }\lambda `$) of 7000. The first two observations used a Tektronix 1024$`\times `$1024 thinned CCD. The 1999 Jul 2 observation used a 2048$`\times `$4096 deep depletion CCD from MIT Lincoln Labs. The pixel scale was 0.23 arcsec with the Tektronix CCD and 0.15 arcsec with the MITLL CCD, which provides good sampling of the seeing profile (1–1.5 arcsec) at each wavelength. The spectra were obtained at four slit position angles (0, 90, 180, and 270). The flat fields were made by combining many exposures with the spectrograph illuminated by a tungsten lamp. Wavelength calibrations were made by observations of a CuAr lamp.

The data were reduced using the FIGARO package. After subtracting the bias level and flat-fielding, the position spectrum was determined by fitting the seeing profile at each wavelength with a Gaussian function. About 1-2$`\times 10^4`$ photons at each wavelength provides a typical accuracy of 2-5 milli-arcsec at the continuum level. Any instrumental effect in the position spectra were eliminated with high accuracy by subtracting those with opposite position angles (0 $`{}_{}{}^{}180^{}`$, or 90 $`{}_{}{}^{}270^{}`$), since these will remain constant on rotation whereas any real structure in the source will reverse its sign. Effects which are removed in this way include any curvature or optical distortion introduced by the spectrograph, any misalignment of the spectrum with the CCD columns, and any effects due to departure of the CCD pixels from a regular grid, imperfect flat fielding, or charge transfer deficiencies in the CCD. The derived position spectra have wavelength coverages of 6520-6750 Å for the former two runs and 6160-6780 A for the latest run, which include H$`\alpha `$, \[SII\] 6716/6731 A, \[NII\] 6583 A, and \[OI\] 6300 Å. Each position spectrum has an arbitrary zero point, which we adjust to correspond to the continuum position. In addition to the position spectra, the intensity spectra were obtained by subtracting the bias, flat-fielding, subtracting the adjacent sky and extracting bright columns on the CCD. In all the spectra, the radial motion of the star is calibrated using the Li 6707.815Å absorption line on the stellar atmosphere. The measured radial velocity in the local standard rest frame is $`+8\pm 2`$ km s-1.

Similar methods were previously used by Hirth, Mundt, and Solf (1994, 1997) to study forbidden emission lines in pre-main sequence stars. However, the techniques described here enable us to achieve substantially higher accuracy in the relative positions.

## 3 Results

### 3.1 Positional displacement and intensity profile of H$`\alpha `$, \[OI\] and \[SII\] emission

Fig. 1 shows the raw spectra of the intensity and the positional displacement. In the intensity spectra, a number of emission lines including H$`\alpha `$, HeI, \[SII\], \[OI\], \[NII\], FeI, FeII, and SiII are detected. The wavelengths and equivalent widths without atmospheric correction are summarized in Table 1.

Among these lines, the H$`\alpha `$, \[SII\] 6731/6716Å, \[NII\] 6583 Å, and \[OI\] 6300 Å lines show displacements from the continuum position in the position spectra. The positions shown in Fig. 1 are the centroid of the combined line and continuum emission of the star, that is, the centroid of their positions weighted by their intensities. To get the true position of the emission line region the observed position must be corrected for the effect of the underlying continuum by dividing the positional displacement in the raw spectra by $`I_{\lambda (\mathrm{line})}/(I_{\lambda (\mathrm{line})}+I_{\lambda (\mathrm{conn})})`$. Fig. 2 shows the positional displacement of H$`\alpha `$, \[OI\] 6300 Å, and \[SII\] 6731/6716 Å observed on 1999 Jul 2 after this calibration. Typical uncertainties for the displacements in Fig. 2 are 2, 15, 40, and 100 milli-arcseconds for H$`\alpha `$, \[OI\] 6300 Å, \[SII\] 6731 Å and 6716 Å, respectively. For H$`\alpha `$, higher accuracy was obtained than that at the continuum level because of its brightness while the accuracy of the other forbidden lines are much less than the continuum because of their lower intensity. The line profiles shown in Fig. 2 are similar to those previously presented by Giovannelli et al. (1995), Reipurth, Pedrosa, & Lago (1996), and Hamann (1994).

Fig. 2 shows that the positional displacement of H$`\alpha `$ emission extends towards the south-west and the north-east over angular scales of up to 20-30 milli-arcsec. The intensity profile of H$`\alpha `$ emission is symmetric except for a redshifted absorption dip at +50 km s-1 from the star, and the blue-shifted wing corresponds to the displacement at the south-west while the red-shifted wing corresponds to that at the north-east. At these wings, the displacement monotonically increases with the relative velocity from the star in the velocity ranges of $`50`$ to $`200`$ km s-1 and $`+200`$ to $`+300`$ km s-1. On the other hand, the displacement is distributed within 5 milli-arcsec from the star in a lower velocity range between $`50`$ and $`+200`$ km s-1, and does not show clear correlation with the velocity in Fig. 2.

The positional displacement of the three forbidden lines extends only towards the south-west, the same direction as those of the blue-shifted H$`\alpha `$ wing, and that of HH55 which lies 3 arcmin from RU Lupi (cf. Krautter, Reipurth, & Eichendorf 1984). The displacements are larger than that of H$`\alpha `$ emission by a factor of 10-30, and smallest in \[OI\] and largest in \[SII\] 6716Å: up to 200, 400, and 600 milli-arcsecs for \[OI\] 6300Å, \[SII\] 6731 Å, and 6716Å, respectively. These line profiles show two peaks with typical velocities between $`250`$ and $`100`$ km s-1 for the high velocity component (HVC) between $`100`$ and $`0`$ km s-1 for the low velocity component (LVC). The HVC is more extended spatially than in the LVC as in other T Tauri stars. The typical angular scales of the LVC are 30-40 milli-arcsec in \[OI\] 6300 Å and 150 milli-arcsecs in \[SII\] 6731 Å. Fig. 2 does not show a clear difference in the displacement between the two components in \[SII\] 6716Å.

The relation between the velocity and the radial distance of the displacement for H$`\alpha `$, \[OI\] 6300Å, and \[SII\] 6731Å, is plotted in Fig. 3. A higher signal to noise ratio of the positional displacements has been obtained by binning along the wavelength. This figure shows that the velocity monotonically increases with the distance for all the lines, and the gradient for H$`\alpha `$ emission is much larger than those for the forbidden lines. Fig. 3 also shows slight asymmetry of the velocity field between the blue-shifted and red-shifted wings of H$`\alpha `$ emission: the positional displacement is prominent in the figure beyond the velocity of 100 km s-1 for the blue-shifted wing and 200 km s-1 for the red-shifted wing. In the forbidden lines, the displayed velocity gradient is larger in the HVC ($`\mathrm{\Delta }`$V $`>`$ 100 km s-1) and smaller in the LVC ($`\mathrm{\Delta }`$V$`<`$100 km s-1). On the other hand, the LVC in \[OI\] reveals no velocity gradient while those in the \[SII\] 6731 Å line show slight positive gradients. The position of the LVC is offset from the star even at the zero velocity: 30 and 70 milli-arcsecs for \[OI\] and \[SII\] 6731 Å, respectively.

The observed tendencies for the line profiles and the positional displacement of the HVC agree with the previous observations of other T Tauri stars. Ground-based long slit spectroscopy shows the positional displacement of the HVCs of T Tauri stars in the Taurus-Auriga cloud (Hirth, Mundt, and Solf 1997), which lies at a similar distance (140 pc) to that of the Lupus cloud. The HVC and LVC in these objects have typical velocities $`50`$ to $`150`$ km s-1 and $`5`$ to $`20`$ km s-1, respectively, which are similar to those of RU Lupi. The measured displacements of the HVC are 0.2 and 0.6 arcsecs for \[OI\] and \[SII\] 6731Å, respectively, which are similar to those seen in our observations. The \[SII\] lines have much lower critical density than the \[OI\] line (typically 103 and 106 cm-3, respectively), and the observed positional difference between these lines indicates that the outer region has a lower electron density. Such a distribution of the electron density could be due to diverging stream lines of the outflow, as suggested by Hirth, Mundt, and Solf (1997).

### 3.2 Time variation of H$`\alpha `$ and \[SII\] emission

Fig. 4 shows the time variation of the intensity profile and the positional displacement for H$`\alpha `$, \[SII\] 6731Å, and \[SII\] 6716Å emission. Since photometric observations were not made on the same dates, each line profile is displayed in a unit of the continuum intensity on each date. The positional displacements are displayed after removal of the contamination of the continuum, as in Fig. 2.

In the case of H$`\alpha `$, time variation was detected for both the line profile and the position spectra. The intensity profiles observed on 1996 Aug 25 and 1997 Jun 27 are nearly the same at the wings. On the other hand, at the peak with the velocity of $`100`$ to $`+100`$ km s-1, the profile on the latter date has lower intensity relative to the continuum than that on the former date. On 1999 Jul 2, the relative intensity of the H$`\alpha `$ to the continuum is slightly lower than the others except that at the absorption dip. In the same observing run, the absorption dip is shallower and relative intensity is higher at the top than those observed in the other runs. Clear change of the position spectra of H$`\alpha `$ was observed during 1997-1999 with the velocity ranges of $`200`$ to $`100`$ km s-1 and $`+100`$ to $`+300`$ km s-1. The positional displacements of these wings are larger on the latest date than in the others by a factor of 2-3. The repeatability of the measured positional displacement between $`100`$ and $`+100`$ km s-1 confirms that the displacement is not due to slight oscillation of the tracking of the telescope, since any other pre-main sequence stars observed in 1999 do not show the same displacements as observed in this object.

The relative intensity of the \[SII\] lines to the continuum is smallest on the earliest date, and largest on the latest date at the whole velocity range. On the other hand, the peak velocity and the line widths of the HVC and LVC did not change throughout the observations. To determine the difference of the time variation between the HVC and LVC, we measure the equivalent widths of the HVC and LVC respectively, and these are tabulated in Table 2. For both the \[SII\] 6716Å and 6731Å, the LVC shows a larger increase in the equivalent width than the HVC: the difference of the relative intensities between 1996 and 1999 is 40-50 and 10-30 percent for the LVC and HVC, respectively. Table 2 also shows the equivalent width ratio of the two lines of 0.53-0.60 and 0.37-0.40 for the HVC and LVC, respectively, and these ratios indicate electron densities of 4-6$`\times 10^3`$ and more than $`10^4`$-$`10^5`$ cm-3 assuming a electron temperature of $`10^4`$ K (cf. Osterbrock 1989). Such a difference of the electron density between the HVC and LVC was also observed in other T Tauri stars by Hamann (1994). On the other hand, Fig 4. shows no systematic time variation in the position spectra of the two \[SII\] lines. The displayed variation would indicate motion of gas, although higher signal-to-noise ratio is required to comfirm the existence of the variation.

In addition to the time variation of the emission line fluxes, variation of the continuum flux can affect the intensity profiles shown in Fig. 4. Its contribution is not clear owing to the lack of photometric data. Previous observations show the time variation of the continuum which may be caused by obscuration by clumpy clouds (Gahm et al. 1974; Giovannelli et al. 1995). In addition, Table 1 shows time variation of the equivalent width for the Li absorption, which suggests variation of the veiling continuum. The intensity profiles of H$`\alpha `$ and the two \[SII\] lines show different time variation, and it cannot be simply explained by variation of the continuum. The ratios of the equivalent width between the LVC and the Li absorption are nearly constant over time (1.4-1.6 and 0.5-0.6 for \[SII\]6731Å/Li and \[SII\]6716Å/Li, respectively), and suggest that the time variation of the LVC is due to the variation of the veiling continuum.

### 3.3 Comparison with the other lines profiles

The profiles of the other lines are plotted in Fig. 5. Most of the permitted lines including FeI, FeII, SiII have triangular profiles with peak velocities of $`50`$ to $`+50`$ km s-1, and full width zero intensities of less than 250 km s-1. The width of these lines is much narrower than that of the H$`\alpha `$, and the velocity ranges of these lines do not overlap that of the H$`\alpha `$ wings: thus, these lines are not related to the H$`\alpha `$ wings observed in 1999. On the other hand, the profiles of \[OI\] 6364 Å and \[NII\] 6583 Å are much more similar to those of \[OI\] 6300 Å and the two \[SII\] lines than that of H$`\alpha `$. The profiles of \[OI\] 6364 Å, which shares the upper transition level with \[OI\] 6300 Å, is nearly the same as that of the latter line shown in Fig. 2. No clear signature of the LVC is found in the \[NII\] profile as in the case of other T Tauri stars (e.g., Hartigan, Edwards, and Ghandour 1995).

The HeI 6678 Å line shows different profiles from the others, and a prominent time variation was observed. In 1996 and 1997, the line had nearly the same triangular profile with a full width zero intensity of about 400 km s-1. On the other hand, in 1999, the line width was reduced to be about a half without changing the peak velocity and the relative intensity to the continuum at the peak. Two components with different line widths exist in the latter profile: a narrow component (NC) with a full width half maximum of about 50 km s-1, and a broad component (BC) with a full width zero intensity of more than 200 km s-1. Such coexistence of two components has been often observed in HeI 5876Å, CaII, and FeII lines of other T Tauri stars (e.g., Batalha et al. 1996, Muzerolle, Hartmann, & Calvet 1998, Beristain, Edwards, & Kwan 1998), and both the components are considered to be associated with mass accretion in the vicinity of the star (e.g., Najita et al. 2000, and references therein). As shown in Fig. 4, H$`\alpha `$ profiles have a red absorption dip indicative of infalling material to the star. The H$`\alpha `$ and HeI profiles observing in 1999 are characterized by a shallower redshifted absorption dip and a weaker broad component, respectively, and such relation between the H$`\alpha `$ dip and the HeI broad component suggests that the latter component arises from gas accreting to the star. No prominent emission associated with the H$`\alpha `$ wings was found in the broad component of HeI in 1999: thus, we conclude that the H$`\alpha `$ wings have a different origin from that of the broad component of the HeI line.

## 4 Discussion

### 4.1 Origin of H$`\alpha `$ emission

Positional displacement as measured by spectro-astrometry should be sensitive to any asymmetric structure about the continuum source. In the case of H$`\alpha `$ emission from pre-main-sequence stars, binary companions are responsible for the displacement as shown by Bailey (1998a). However, the feature seen in RU Lupi is quite unlike those normally seen in binary systems which are not symmetric about the line centre. Binary stars with the same spectrum and different radial velocities could produce the observed positional displacement, although such motion at the observed direction does not agree with the orbital motion perpendicular to the outflow. No binary companion has been detected by $`the`$ $`Hubble`$ $`Space`$ $`Telescope`$ (Bernacca et al. 1995) or infrared speckle observations (Ghez et al. 1997). Thus, it is likely that the detected feature reflects the distribution of circumstellar matter around the star.

A variety of mechanisms have been previously proposed to explain the observed H$`\alpha `$ profiles from pre-main-sequence stars. Magnetically-driven stellar winds which accelerate with distance have been suggested for Balmer emission including H$`\alpha `$ (e.g., Decampli 1981; Hartmann, Edwards, & Avrett 1982; Lago 1984; Natta, Giovanardi, & Palla 1988), while a decelerating ”stochastic” stellar wind was proposed by Mitskevich, Natta, & Grinin (1993). Calvet & Hartmann (1992) and Hartmann, Hewett & Calvet (1994) argue that the H$`\alpha `$ emission is produced in infalling envelopes. Edwards et al. (1994) presented for 15 T Tauri stars a set of high resolution profiles including Balmer, HeI, and NaI lines and showed that the blue/red asymmetry of upper Balmer lines agrees well with magnetospheric accretion models. On the other hand, Reipurth, Pedrosa, and Lago (1996) compared their H$`\alpha `$ profiles with various models, and suggested that a single model cannot account for the variety of the H$`\alpha `$ profiles. In addition to outflowing and inflowing motion, Keplerian rotation would also be responsible for the shape of the H$`\alpha `$ profile.

Among these candidates, we consider bipolar outflow as a plausible explanation for the observed displacement for the following reasons: (1) the positional displacement of the blue-shifted component is in the same direction as those of forbidden lines, which are considered to trace outflow, and (2) the other permitted lines, which probe mass accretion, do not have such broad wings as H$`\alpha `$, thus these wings should have different origin (cf. §3-3). No clear evidence for inflowing motion has been detected in the position spectra, consistent with previous models which propose that the H$`\alpha `$ emission arises from a region within less than ten stellar radii (e.g., Calvet & Hartmann 1992; Hartmann, Hewett & Calvet 1994). Rotation cannot explain the observed displacement for the following reasons: (1) the measured velocity increases with positional displacement while Keplerian rotation should behave in an opposite fashion, and (2) the star needs to be more than 100 M to bound the gas with the measured positional displacement and velocity (3 au and 200 km s-1, respectively) while previous research suggests the mass of the star is 0.5-2 M (Gahm et al. 1974; Lamzin et al. 1996).

Figs 2 and 3 show that the velocity of the outflowing gas increases with positional displacement, a tendency which agrees with that of magnetically-driven wind models which predict that the flow accelerates with distance (e.g., Decampli 1981; Hartmann, Edwards, & Avrett 1982; Lago 1984). Most of the models adopt a spherically symmetric line forming region not exceeding ten stellar radii, and predict that the acceleration of the flow ends within this scale. On the other hand, the observed position spectra of RU Lupi indicate bipolar geometry and a velocity gradient extending to a few au (Figs 2 and 3). The latter results would indicate slower acceleration of the outflow than those predicted by models, however, models with more appropriate geometry and direct prediction of the positional displacement are necessary for detailed comparison. Actually, H$`\alpha `$ emission from pre-main-sequence stars should be optically thick (cf. Reipurth, Pedrosa, & Lago 1996), thus the observed velocity field would be highly sensitive to the geometry which affects radiative transfer.

The time variation of the positional displacement and the intensity profile for the H$`\alpha `$ wings are similar to that of HeI in the following points: there is no clear difference between 1996 and 1997, and a remarkable difference between 1997 and 1999. Since the HeI emission is related to the accreting gas, as shown in §3-3, such similar time variation suggests that the accreting gas also contributes to the H$`\alpha `$ flux at least in 1996 and 1997. Even if the positional displacement of the bipolar outflow is constant over time, the net positional displacement can change as the flux of the accreting gas at the centre varies. Transient time variation of the displacement in high-velocity wings could be explained by direct motion of the outflowing ionized gas, although we reject this possibility since line profiles with the same velocities only show slight change of line-to-continuum ratio at the same velocities.

Here we attempt to explain the observed positional displacement with a steady bipolar outflow and time-variable flux from inflow. To remove contamination from inflow on the total H$`\alpha `$ flux, and reproduce the position spectra of the outflow we make the the following assumptions: (1) the H$`\alpha `$ emission from infalling gas arises from the star and cannot be spatially resolved by our observations, (2) inflow contributes to the observed H$`\alpha `$ profiles only in 1996 and 1997, and (3) the time variation of the Li equivalent widths shown in Table 1 is due to the veiling continuum. The first assumption provides the observed net positional displacement as follows:

$$\mathrm{\Delta }x_\lambda =\mathrm{\Delta }x_{\lambda \mathrm{out}}\frac{I_{\lambda (\mathrm{out})}}{I_{\lambda (\mathrm{out})}+I_{\lambda (\mathrm{in})}}$$
(1)
where $`\mathrm{\Delta }x_\lambda `$ is the net displacement shown in Fig. 4, $`\mathrm{\Delta }x_{\lambda \mathrm{out}}`$ is the positional displacement caused by the bipolar outflow, $`I_{\lambda (\mathrm{out})}`$ and $`I_{\lambda (\mathrm{in})}`$, are the absolute H$`\alpha `$ intensity from the outflow and inflow, respectively. This assumption is consistent with our observations, and the previous inflow models described above. From the second assumption, the intensities from the outflow and inflow are obtained as follows:

$$I_{\lambda (\mathrm{out})}=I_{\lambda (\mathrm{H}\alpha 99)},I_{\lambda (\mathrm{in})}=I_{\lambda (\mathrm{H}\alpha )}I_{\lambda (\mathrm{H}\alpha 99)},$$
(2)
where $`I_{\lambda (\mathrm{H}\alpha )}`$ is the absolute H$`\alpha `$ intensity in a certain observing run, and $`I_{\lambda (\mathrm{H}\alpha 99)}`$ is that observed in 1999. This assumption is plausible for the H$`\alpha `$ wings since those observed in 1999 are not associated with the other permitted lines indicative of mass accretion as described in 3-3. The third assumption provides the absolute H$`\alpha `$ intensities from the observed parameters as follows:

$$\frac{I_{\lambda (\mathrm{H}\alpha )}}{I_{\lambda (\mathrm{H}\alpha 99)}}=\frac{i_{\lambda (\mathrm{H}\alpha )}I_{\lambda (\mathrm{conn})}}{i_{\lambda (\mathrm{H}\alpha 99)}I_{\lambda (\mathrm{conn99})}}=\frac{i_{\lambda (\mathrm{H}\alpha )}}{i_{\lambda (\mathrm{H}\alpha )}}\frac{W_{\lambda (\mathrm{Li99})}}{W_{\lambda (\mathrm{Li})}},$$

where $`i_{\lambda (\mathrm{H}\alpha )}`$ and $`I_{\lambda (\mathrm{conn})}`$ are the relative H$`\alpha `$ intensity to the continuum and the absolute intensity of the continuum in a certain observing run, respectively, and $`i_{\lambda (\mathrm{H}\alpha 99)}`$ and $`I_{\lambda (\mathrm{conn99})}`$ are those observed in 1999. $`W_{\lambda (\mathrm{Li})}`$ and $`W_{\lambda (\mathrm{Li99})}`$ are the equivalent widths of the Li absorption in a certain observing run and in 1999, respectively. By combining equations (1)-(3), the positional displacement of the outflow is obtained from the observed parameters as follows:

$$\mathrm{\Delta }x_{\lambda (\mathrm{out})}=\mathrm{\Delta }x_\lambda \frac{i_{\lambda (\mathrm{H}\alpha )}}{i_{\lambda (\mathrm{H}\alpha 99)}}\frac{W_{\lambda (\mathrm{Li99})}}{W_{\lambda (\mathrm{Li})}},$$
(3)
Fig. 6 shows the reproduced position spectra using equation (4). The spectra for 1999 are identical with the observed one shown in Fig. 4, since we assume that the outflow dominates the observed intensity in this observing run. These spectra agree well with each other, thus support the idea that the time variation of the positional displacement can be explained by the combination of steady bipolar outflow and time-variable flux from the accreting gas.

The contribution of the inflow to the total H$`\alpha `$ flux could be estimated from equations (2) and (3): according to these equations, 32 and 28 percents of the total H$`\alpha `$ flux would arise from the accreting gas in 1996 and 1997, respectively. These values would be lower limits, since the inflow with low velocities would contribute to the H$`\alpha `$ emission in 1999. The permitted lines related to mass accretion were detected with a typical velocity range between $`150`$ and $`+150`$ km s-1 as shown in §3-3, and there is no reason to assume that H$`\alpha `$ emission does not arises from the same region.

### 4.2 Origin of the forbidden lines

Optical forbidden emission lines including \[OI\] 6300Å and \[SII\] 6716/6731Å provide information on the outflows from T Tauri stars. These line profiles often have two blue-shifted velocity components: a high-velocity component (HVC) and low-velocity component (LVC) which have typical velocities of $``$50 to $``$150 km s-1 and $`5`$ to $`20`$ km s-1, respectively (Hirth, Mundt, and Solf 1997). The HVC is spatially resolved by ground-based long-slit spectroscopy in the nearest star-forming regions (e.g., Solf & Böhm 1993; Hirth, Mundt, and Solf 1997). Hirth, Mundt, and Solf (1997) have shown that the typical angular scale of the HVC is 0.2–0.6 arcsec from their large sample of data mainly in the Taurus-Auriga cloud. The high speed and narrow linewidth of the HVC indicates that the velocity vectors of the contributing gas particles are nearly parallel to one another, i.e., suggesting a well-collimated flow or jet (Kwan & Tademaru 1988). On the other hand, the low-velocity component (LVC) has smaller angular scales in the same sample (0.1-0.2 arcsec or less), and its origin is less understood. Geometric effects of a single outflow are proposed to explain the existence of the two velocity components (e.g., Edwards et al. 1987; Hartmann & Raymond 1989), while two distinct mechanisms, jet and disc wind, are proposed for each velocity component by Kwan and Tademaru (1988, 1995). Hamann (1994) measured the electron density, temperature, and the ionization state between the HVC and LVC in many T Tauri stars, and the measured difference of the physical conditions suggest distinct origins for these two components as proposed by Kwan and Tademaru (1988, 1995). Similar results were also presented by Hartigan, Edwards, & Ghandour (1995). The velocity dispersion and flow velocity are comparable in the LVC, suggesting moderate collimation such as that due to a disc wind (e.g., Solf & Böhm 1993). Most pre-main-sequence stars do not have red-shifted components to their forbidden lines, and this is explained by obscuration of these components by circumstellar discs (e.g., Eislöffel et al. 2000, and references therein).

The \[OI\] and \[SII\] lines from RU Lupi have two blue-shifted components, and our results also support the existence of two distinct mechanisms. The \[SII\] intensity ratio indicates different electron densities between the HVC and LVC as shown in §3-2 (4-6$`\times 10^3`$ and more than $`10^4`$-$`10^5`$ cm-3 assuming a electron temperature of $`10^4`$ K), a tendency also observed in other T Tauri stars by Hamann (1994). In addition, the intensity of each component exhibits different time variations as described in §3-2. The coexistence of two distinct mechanisms can easily explain these differences.

Furthermore, the velocity field in the LVC shown in Fig. 3 agrees with the disc wind hypothesis. Fig. 3 shows that the position of the LVC is offset from the star even at zero velocity and nearly constant over all velocities in \[OI\] emission, as described in §3-2. The former cannot be explained by acceleration of outflow, while both tendencies can easily be explained by moderate collimation and/or rotation which would be expected in a disc wind. In the case of \[SII\] 6731Å, the positional displacement in the LVC increases with the velocity. Such discrepancy between the \[OI\] and \[SII\] emission can be explained by blending of the HVC and LVC near their boundary. The profiles of the forbidden lines in Fig. 2 show that the two components are partially blended, thus contamination from the HVC could increase the positional displacement in the LVC near their boundary. It is likely that this effect is more prominent for the \[SII\] line, since the lower critical density of this line allows for larger contamination from the HVC which has a lower electron density (cf. §3-2).

### 4.3 Circumstellar disc and disc gap

As shown in Fig. 2, the forbidden lines have only blue-shifted components, similar to many other T Tauri stars, a fact usually explained by obscuration of the red-shifted flow by a circumstellar disc (e.g., Eislöffel et al. 2000, and references therein). Infrared excess and a bright millimetre continuum suggest the existence of a circumstellar disc around RU Lupi (e.g., Giovannelli et al. 1995; Carballo et al. 1992; Nürnberger et al. 1997), supporting the explanation for the lack of red-shifted components.

On the other hand, H$`\alpha `$ emission has both blue- and red-shifted components. The measured positional displacement for the H$`\alpha `$ emission is smaller than any velocity component of \[OI\] and \[SII\] emission. The detection of the red-shifted H$`\alpha `$ outflow can be explained by a gap in the disc which allows this emission to be seen. The presence of disc holes and gaps is often inferred around other pre-main-sequence stars by infrared spectral energy distributions (SEDs) (e.g., Mathieu et al. 1991; Marsh & Mahoney 1992; Jensen & Mathieu 1997), and directly observed in near-infrared imaging (e.g., Roddier et al. 1996; Silber et al. 2000) and at radio wavelengths (e.g., Dutrey, Guilloteau, & Simon 1994). Fig. 3 suggests that the projected radius of the hole or gap is 20-30 milli-arcsec, which corresponds to 3-4 au, assuming symmetric distribution of the blue-shifted and unseen red-shifted components of the forbidden lines.

To investigate the presence of a disc gap, we plot the infrared SED in Fig. 7. The flux at $`JHKLM`$-bands and 12, 25, 60, 100 $`\mu `$m were obtained from Giovannelli et al. (1995) and Carballo et al. (1992). Giovannelli et al. (1995) show the near-infrared flux is highly variable with time, thus the data-set with the largest and smallest flux were selected from their data and plotted in the figure. In Fig. 7, a spectrum obtained via the $`Infrared`$ $`Space`$ $`Observatory`$ (ISO) data archive is also plotted to interpolate the flux in the mid-infrared. Lamzin et al. (1996) claim that the mid-to-far infrared flux obtained by Carballo et al. (1992) would be overestimated, since another source discovered by Cohen & Schwartz (1987) is included in the same aperture. However, the consistency of their data with the ISO spectrum shows that the contribution from this source is much less than the total flux: the aperture of the ISO Short-Wave Spectrometer (SWS) with a size of 20$`\times `$33 arcsec2 does not include this source which lies about 3 arcmin away from the target.

Fig. 7 shows two peaks in the SED which arises from dust. The near-infrared peak at 1.65 $`\mu `$m is longer than the peak for the star (1.0 $`\mu `$m) with an effective temperature of 3900 K (Nürnberger et al. 1997; Lamzin et al. 1996). Previous authors show that the selective extinction towards RU Lupi is 0.3-1.3 AV (Hughs et al. 1994; Giovannelli et al. 1995). Such values can shift the peak of the stellar continuum by less than 0.2 $`\mu `$m suggesting that the near-infrared excess from the stellar photosphere represents hot dust from the accretion disk. Such a contribution by the accretion disk to the near-infrared SED is also supported by the near-infrared colours (J-H and H-K of 0.84-0.92 and 0.72-0.77, respectively) and diagrams shown by Hughes et al. (1994), Greene & Meyer (1995), and Itoh et al. (1999). On the other hand, the mid-to-far infrared SED shows another excess which is due to cooler dust with a temperature of less than 100 K. The shallow dip at 4-15 $`\mu `$m suggests a lack of dusty material at temperatures of 200-900 K. These temperatures correspond to the radiative temperature at a radius of 0.1-2 au from the star, assuming the total luminosity from the star and accreting gas of 5 L based on Lamzin et al. (1996). These spatial scales agree well with that of the suggested gap in Figs 2 and 3.

Theoretical work has shown that disc gaps and holes could be induced by binary companions or young planets (e.g., Lubow & Artimowicz 2000). In the case of RU Lupi, no other pre-main-sequence companions have been detected by $`Hubble`$ $`Space`$ $`Telescope`$ (Bernacca et al. 1995), infrared speckle observations (Ghez et al. 1997), or spectro-astrometry which is sensitive to the presence of binaries (cf. Bailey 1998a). Therefore, the observed disc gap suggests the existence of an unseen companion. A young planet could be one of possible candidates for the following reasons: the presence of the fragmented dusty clouds around RU Lupi are suggested by variability at optical to infrared wavelengths (Gahm et al. 1974; Giovannelli et al.1995), and Gahm et al. (1974) claim that such clouds would indicate ongoing planetary formation around this object. In addition, the observed disc gap has a similar size to that of the orbit of the Jovian planets in our solar system, and Takeuchi, Miyama, & Lin (1996) suggest that a proto-Jupiter could induce a disc gap with a similar size to its orbit. Their numerical calculations show that a Saturnian and Jovian mass protoplanet can open a gap in 103-104 and 102-103 yrs, respectively, timescales which are much shorter than the age of RU Lupi (106 yrs - Lamzin et al. 1996).

## 5 Conclusions

Spectro-astrometric observations have been used to determine the spatial distribution of optical emission lines from RU Lupi. Positional displacements were detected in H$`\alpha `$, \[OI\] 6300 Å, \[SII\] 6731/6716 Å, and \[NII\] 6583 Å with accuracies down to a few milli-arcsec. The intensity profiles were also obtained for these lines together with HeI, FeI, FeII, SiII, and another \[OI\] line.

The positional displacement of H$`\alpha `$ emission extends towards the south-west and north-east with angular scales of 20-30 milli-arcsecs. The blue-shifted wing in the intensity profile corresponds to the displacement at the south-west while the red-shifted wing corresponds to that at the north-east. We conclude that the positional displacement of the H$`\alpha `$ emission is due to bipolar outflow for the following reasons: (1) the distribution of the positional displacement is almost symmetric between the blue-shifted and red-shifted components, (2) the displacement of the blue wing aligns well with that of the forbidden lines which are considered to trace outflowing gas. The velocity of the outflow increases with the distance, suggesting the flow is magnetically-driven. On the other hand, the time variation of the intensity profile and the positional displacement can be explained well by the variation of H$`\alpha `$ flux from the accreting gas at the star. We estimate the contribution of the accreting gas to the total flux to be more than 32 and 28 percent in 1996 and 1997, respectively.

The positional displacement of \[OI\] and \[SII\] emission extends towards the south-west, in the same direction as HH55 which lies 3’ from the star. The intensity profiles have two blue-shifted components, as seen in many other pre-main-sequence stars, and the positional displacements of hundreds of milli-arcsecs and down to 30 milli-arcsecs were detected for the high-velocity components (HVC) and the low-velocity components (LVC), respectively. The differences of the derived electron density and time variation between the two components suggest that these components have distinct origins, e.g., jet and disc wind. The positional displacement in the LVC is consistent with the LVC originating from a disc wind.

The forbidden lines have only blue-shifted components, as do those in many other T Tauri stars, and this is usually explained by the obscuration of the red-shifted flow by a circumstellar disc. On the other hand, the H$`\alpha `$ outflow has both blue-shifted and red-shifted components. Such a difference suggests the presence of a disc gap with an outer radius of 3-4 au. The infrared spectral energy distribution is consistent with the presence of a gap on this scale. This gap could be induced by a unseen companion such as a young planet.

We would like to thank the staff at the Anglo-Australian Observatory for their help and support during the observations. We also thank ISO staff for their assistance in using their archive system, and the referee (T.P. Ray) for his valuable comments. MT thanks PPARC for support through a PDRA.

References

Bailey J., 1998a, MNRAS, 301, 161
Bailey J., 1998b, SPIE Proceeing, 3355, 932
Batalha C.C., Stout-Batalha N.M., Basri G., Terra M.A.O., 1996, ApJS, 103, 211
Beristain G., Edwards S., Kwan J., 1998, ApJ, 499, 828
Bernacca P.L., Lattanzi M.G., Porro, I., Neuhäuser R., Bucciarelli B., 1995, A&A, 299, 933
Burrows C.J., Stapelfeldt K.R., Watson A.M., Krist J.E., Ballester G.E., Clarke J.T., Crisp D., Gallagher J.S.III, Griffiths R.E., Hester J.J., Hoessel J.G., Holtzman J.A., Mould J.R., Scowen P.A., Trauger J.T., Westphal J.A. 1996, ApJ, 473, 437Calvet N., Hartmann L., 1992, ApJ, 386, 239
Calvet N.,Hartmann L., Hewett R., 1992, ApJ, 386, 229
Carballo R., Wesselius P.R., Whittet D.C.B., 1992, A&A, 262, 106
Cohen M., Schwartz R.D., 1987, ApJ, 316, 311
Decampli W.M., 1981, ApJ, 244, 124
Devaney M.N., Thiébaut E., Foy R., Blazit A., Bonneau D., Bouvier J., de Batz B., Thom Ch., 1995, A&A, 300, 181
Dutrey A., Guilloteau S., Simon M., 1994, A&A, 286, 149
Edwards S., Cabrit S., Strom S.E., Heyer I., Strom K.M., Anderson E., 1987, ApJ, 321, 473
Edwards S., Hartigan P., Ghandour L., Andruis C., 1994, AJ, 108, 1056
Eislöffel J., Mundt R., Ray T.P., Rodríguez L.F, 2000, Protostars and Planets IV, 815
Gahm G.F., Nordth H.L., Olofsson S.G., Carlborg N.C.J., 1974, A&A, 33, 399
Gahm G.F., Lago M.T.V.T., and Penston M.V., 1981, MNRAS, 195, 59
Ghez, A.M. McCarthy D.W., Patience J.L., Beck T.L., 1997, ApJ, 481, 378
Giovannelli F., Vittone A.A., Rossi C., Errico L., Bisnovatyi-Kogan G.S., Kurt V.G., Lamzin S.A., Larionov M., Sheffer E.K., Sidorenkov V.N. 1995, 114, 341Greene T.P., Meyer M.R., 1995, ApJ, 450, 233
Hamann F., 1994, ApJS, 93, 485
Hamann F. Persson S.E., 1992, ApJS, 82, 247
Hartigan P., Edwards S., Ghandour L., 1995, ApJ, 452, 736
Hartmann L., Hewett R., Calvet N., 1994, ApJ, 426, 669
Hartmann L., Edwards S., Avrett E., 1982, ApJ, 261, 279
Hartmann L. Raymond J.C., 1989, ApJ, 337, 903
Hirth G.A., Mundt R., Solf J., 1994, A&A, 285, 929
Hirth G.A., Mundt R., Solf J., 1997, A&AS, 126, 437
Hughes J., Hartigan P., Clampitt L., 1993, AJ, 105, 571
Hughes J., Hartigan P., Krautter J., Kelemen J., 1994, AJ, 108, 1071
Itoh Y., Tamura M., Nakajima T., 1999, AJ, 117, 1471
Jensen E.L.N. Mathieu R.D., 1997, AJ, 114, 301
Krautter J., Reipurth B., Eichendorf W., 1984, A&A, 133, 169
Kwan J. Tademaru E., 1988, ApJ, 332, L41
Kwan J. Tademaru E., 1995, ApJ, 454, 382
Lamzin S.A., Bisnovatyi-Kogan G.S., Errico L., Giovannelli F., Katysheva N.A., Rossi C., Vittone A.A. 1996, A&A, 306, 877Lago M.T.V.T, 1984, MNRAS, 210, 323
Lago M.T.V.T, Penston M.V., 1982, MNRAS, 198, 429
Lubow S., Artimowicz P., 2000, Protostars and Planets IV, 731
Marsh K.A. Mahoney M.J., 1992, ApJ, 395, L115
Mathieu R.D., Adams F.C., Latham D.W., 1991, AJ, 101, 2184
Mitskevich A.S., Natta A., Grinin V.P., 1993, ApJ, 404,751
Muzerolle J., Hartmann L., Calvet N., 1998, AJ, 116, 455
Najita J., Edwards S., Barsi G., and Carr J., 2000, Protostars and Planets IV, 457 Natta A., Giovanardi C., Palla F., 1988, ApJ, 332, 921
Nürnberger D., Chini R., and Zinnecker H., 1997, A&A, 324, 1036
Osterbrock D.E., 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Mill Valley: Univ. Science Books)
Osterbrock D.E., Tran H.D., Veilleux S., 1992, ApJ, 389, 305
Ray T.P., Mundt R., Dyson J.E., Falle S.A.E.G., Raga A.C., 1996, ApJ, 468, L103
Reipurth B., Pedrosa A., Lago M.T.V.T., 1996, A&AS, 120, 229
Roddier C., Roddier F., Northcott M.J., Graves J.E., Jim K. 1996, ApJ, 463, 326
Silber J.M., Gledhill T.M., Duchêne G., Ménard F. 2000, ApJL, 536, 89
Solf J., Böhm K.H., 1993, ApJ, 410, L31
Takeuchi T., Miyama S.M., Lin D.N.C., 1996, ApJ, 460, 832


text:
# Dynamical scaling of the DNA unzipping transition

## Abstract

We report studies of the equilibrium and the dynamics of a general set of lattice models which capture the essence of the force-induced or mechanical DNA unzipping transition. Besides yielding the whole equilibrium phase diagram in the force vs temperature plane, which reveals the presence of an interesting re-entrant unzipping transition for low $`T`$, these models enable us to characterize the dynamics of the process starting from a non-equilibrium initial condition. The thermal melting of the DNA strands displays a model dependent time evolution. On the contrary, our results suggest that the dynamical mechanism for the unzipping by force is very robust and the scaling behaviour does not depend on the details of the description we adopt.

The replication of DNA is a correlated process involving many proteins and other molecules working at different points in space and time. An understanding of the nature and origin of this correlation is expected to shed light on this complex mechanism. It has recently been shown that the force induced unzipping of DNA is a genuine phase transition different from the thermal melting transition of DNA. It was then hypothesized  that the initiation of replication at the origins along the DNA, e.g, by dnaA for E.Coli or by the “origin recognition complex” (ORC) in eukaryotes is like this unzipping near the critical threshold (with dnaA or ORC acting as the force-inducing agent) and the resulting correlation during unzipping leads the co-operativity required for replication.

In contrast to real biological situations, techniques like laser tweezers, atomic force microscopes (AFM) etc have been used to study DNA by pulling at one end. This has led to strand separation by force. In particular, AFM experiments reported hysteresis in the unzipping process, indicating the presence of a first order transition. These mechanical unzipping experiments have opened up new ways of thinking about DNA, just as similar stretching experiments of DNA showed the possibility of several structures other than the most prevalent B-DNA. The activities of polymerases, topoisomerase etc on single stranded DNA have now been analyzed in terms of the force they exert or the force applied against them. What needs to be investigated, to mimic the biological situation, is the coupling between the opening of the strands and the subsequent events during replication. Such a study involves the dynamics of the unzipping process.

The purpose of this paper is to define a set of simpler models, in the spirit of Poland and Sheraga, for which the unzipping transition can be studied exactly. On the basis of this, the dynamics can be understood. The proposed lattice models (bubble models: b-models) incorporate the mutual-avoidance (hard-core repulsion) of the strands (and also self-avoidance). A further simplification is obtained by suppressing the bubbles along the chains, thereby defining a “fork model” or “Y-model”. The phase diagram of the equilibrium system displays in both cases a re-entrant region at low $`T`$: for a finite range of forces the molecule gets unzipped by decreasing the temperature. The dynamics of both the b- and the Y-models in the various phases and on the phase boundary are then studied, by starting from a zig-zag non-equilibrium bound state as the initial condition. We find that in all regimes above or below the critical line in the phase diagram, the time evolutions of the order parameters follow dynamical scaling laws.

The models and their equilibrium phase diagram: We model the two strands of DNA by two directed self and mutually avoiding walks. In two dimensions, on the square lattice (see Fig.1), this means that the two walks are forced to follow the positive direction of the diagonal axis $`(1,1)`$ (i.e. the coordinate along the direction $`(1,1)`$, to be called the $`z`$ direction, always increases). The force is then acting along the direction perpendicular to direction $`z`$, i.e., $`(1,1)`$: this transverse direction is called the $`x`$-direction. By measuring the $`x`$ separation in unit of the elementary square diagonal, we say that two complementary monomers are in contact when $`x=1`$: a (positive) binding energy $`ϵ`$ is gained for each contact. Due to the geometrical properties of the lattice, all these contacts contributing to the energy involve monomers labelled by the same $`z`$-coordinate, as one would require for base pairing in DNA (see Fig. 1).

By fixing $`ϵ=+1`$ throughout the calculations, all the thermodynamic properties of the system depend on the temperature $`T=\beta ^1`$ and on the force $`f`$. We have set $`k_B=1`$. The partition function of the two $`n`$-step chains, one starting at $`(0,1)`$ and the other at $`(1,0)`$ is

$$𝒵_n(\beta ,f)=\underset{x1}{}d_n(x)\mathrm{exp}\left(\beta fx\right)e^{n\beta },$$
(1)
where $`d_n(x)`$ represents the fixed distance partition function, i.e., the sum over all interacting pairs of directed chains whose last monomers are at distance $`x`$ and $``$ is the free energy density. The quantity $`d_n(x)`$ whose $`\beta `$ dependence has been omitted, obeys simple recursion relations. In $`d=2`$, on the square lattice, the recursion relation is

$`d_{n+1}(x)`$ $`=`$ $`\left[2d_n(x)+d_n(x+1)+d_n(x1)\right]`$ (3)
$`\times \left[1+\left(e^\beta 1\right)\delta _{x,1}\right],`$
valid for $`x1`$ with the conditions $`d_0(x)=\delta _{x,1}`$ and $`d_n(0)=0`$ for all $`n`$. Note that this kind of equation appears also in other models of DNA and in studies of random walks adsorption. The model can be asymptotically solved by locating the singularity closest to the origin of its related generating function

$$𝒢(\beta ,f)=\underset{n=0}{\overset{\mathrm{}}{}}z^n𝒵_n(\beta ,f).$$
(4)
A phase transition is indicated by any singular change in the location of the closest singularity.

With some simple calculations a critical line, $`f_c(T)`$, separating the zipped from the unzipped phase can be obtained

$$f_c(T)=T\mathrm{cosh}^1\left[\frac{1}{2}\frac{1}{\sqrt{1e^\beta }1+e^\beta }1\right],$$
(5)
and it has been plotted in Fig. 2.

At $`f=0`$ the critical melting temperature is $`T_m=\frac{1}{\mathrm{log}4/3}`$, the critical force at $`T=0`$ is $`f_0=1`$ whereas for $`f>f_{\mathrm{max}}=1.358806\mathrm{}`$ the system is always in the unzipped state. Similar results can be obtained for the physically relevant three dimensional case and for higher dimensions. The average fraction of contacts $`\mathrm{\Theta }`$ is zero (not zero), when $`f>f_c(T)`$ ($`f<f_c(T)`$) while the average end-to-end distance $`x=N\frac{}{f}`$, where $`N`$ is the number of base pairs in each chain, and the average number of “liberated” monomers (i.e. from the last contact to the end), $`m`$ (See Fig. 1) have the opposite behavior. The transition line is first order at any point with $`f_c(T)0`$ in $`d=2,3,4,5`$ and everywhere for $`d>5`$.

We notice that there is a re-entrance in the $`fT`$ phase diagram. For $`f_{\mathrm{max}}>f>f_0`$ the usual denaturation transition is present but if the temperature is further lowered the two strands separate again through a “cold unzipping”.

A further simplification to this model can be obtained by suppressing the bubbles along the chain, i.e. by considering only conformations having the first $`Nm`$ monomers bounded, whereas the remaining $`m`$ are separated in a Y-like conformation.

This Y-model will be extremely useful to study the unzipping dynamics and it presents a phase diagram similar to the one previously obtained but, for example, with $`T_m=\frac{1}{\mathrm{log}2}`$, $`f_0=1`$ and $`f_{\mathrm{max}}=1.282143\mathrm{}`$ in d=2.

Dynamics. We now consider the dynamics of the models previously introduced, the b-model of directed walks and the simplified Y-model. In both cases, we start from a non-equilibrium “zero-temperature” initial condition with the two chains zipped in a zig-zag configuration, and let the system evolve at a temperature $`T`$ and under a force $`f`$, with $`T`$ and $`f`$ chosen so that the equilibrium state is either on or above the critical line. The simpler 2-dimensional case will be presented first and the generalization to higher dimensions will be discussed later.

The five regimes considered are marked A-E in Fig. 2. Numerically, a Monte-Carlo dynamics is used to monitor the time evolution of the order parameters $`m`$ and $`x`$ previously introduced (see also Fig. 1). For the dynamics, we selected only local physical moves, so that the model should be the lattice counterpart of the Rouse model in the continuum. In all cases we find that, far from saturation, the order parameters evolve as power laws as functions of time $`t`$: in particular we can define two exponents $`\theta _1`$ and $`\theta _2`$ as

$$m(t)t^{\theta _1},\mathrm{and}x(t)x(0)t^{\theta _2},$$
(6)
Notice that $`m`$ and $`x`$ in Eq. (6) refer to the ensemble average values, but average signs have been omitted for simplicity of notation. We in fact find the following dynamical scaling laws to hold:

$$\frac{m(t)}{N^{d_1}}=G_m\left(\frac{t}{N^{z_1}}\right),\frac{x(t)x(0)}{N^{d_2}}=G_x\left(\frac{t}{N^{z_2}}\right),$$
(7)
where $`N`$ is the length of the chains and $`G_{m,x}`$ are two scaling functions. Eq. 7 also defines the exponents $`d_{1,2}`$ and $`z_{1,2}`$ for the two variables. Note that $`\theta _1=\frac{d_1}{z_1}`$, $`\theta _2=\frac{d_2}{z_2}`$ and that $`d_{1,2}`$ can be obtained through equilibrium considerations as one requires $`m(t)\stackrel{t\mathrm{}}{}N^{d_1}`$ and $`x(t)\stackrel{t\mathrm{}}{}N^{d_2}`$. The crossover to the equilibrium behavior is described by the “dynamic” exponents $`z_{1,2}`$. The possibility of the two quantities having different relaxations is kept open by two different exponents.

The values of the exponents in Eqs. (6) and (7) obtained from simulations, for both the Y-model and the b-model, are shown in Table I. They were obtained by collapsing Monte-Carlo data according to eq. (7) and by using a recently proposed search algorithm  (see Fig. 3)

It is possible to explain the exponents for the Y-model found numerically. There are two relevant mechanisms which can drive the separation of the strands: the unbinding of bases at the bifurcation point and the stretching of the strands at the extremes. The combination of these two processes controls the overall behavior in different regimes.

Regime A: $`f=0`$, $`T>T_m`$: Above the critical temperature the dominance of the entropy implies that at every time step one base pair breaks, yielding a linear behavior with $`\theta _1=1`$ and $`d_1=1`$. Also $`x(t)`$ tends to increase, up to its equilibrium value $`N^{\frac{1}{2}}`$: this is reflected in the equilibrium probability distribution $`P_{eq}(x)`$ which displays an upward derivative at $`x=1`$. This suggests (see ) that the dynamics of this quantity should be in the same universality class of the one-dimensional Kardar-Parisi-Zhang equation, and so $`\theta _2=\frac{1}{3}`$, $`d_2=\frac{1}{2}`$, and $`z_2=\frac{3}{2}`$.

Regime B: $`f=0`$, $`T=T_m`$: In this regime $`\theta _1=\frac{1}{2}=z_1^1`$ because at criticality the probabilities to increase and to decrease $`m`$ are expected to be equal, so that $`m`$ performs, roughly speaking, a random walk in time with reflecting boundaries at $`m=0,N`$. Also for the end-to-end distance, steps toward larger or smaller values of $`x`$ are equally probable, and this means that the equation obeyed by $`x`$ should stay in the same universality class as the one-dimensional Edward-Wilkinson equation, and therefore $`\theta _2=\frac{1}{4}`$ and $`z_2=2`$.

Regime C: $`f>0`$,$`T>T_m`$: The strands microscopically tend to stretch along the direction of the pulling force. However, once we have pulled the two chains up to an end-to-end distance $`x`$, to increase $`x`$ further by one unit we first need to move all the stretched part, which would take a time typically of order $`x`$. In other words, one has $`x(t+t_0)x(t)+\frac{t_0}{x}`$ and we expect the dynamical exponent $`\theta _2`$ to be $`\frac{1}{2}`$ (of course $`\theta _1=1`$ as before and $`d_1=d_2=1`$).

Regime D: $`T<T_m`$,$`f>f_c(T)`$: Here the only microscopic mechanism for opening the fork is through the applied force: the strands must stretch completely in the vicinity of the bifurcation point and only at this point will the fork liberate one more monomer, because otherwise the opening is energetically very unlikely. Thus we expect that $`xm`$ and, using arguments as done for regime $`C`$, $`\theta _1=\theta _2=\frac{1}{2}`$, $`d_1=d_2=1`$. In ref. it was found, in the mean field approach and in a model resembling our Y-model, that the time necessary to unzip the two strands completely is of the order of $`N^2`$. This is in agreement with our analysis, but works only in this regime.

Regime E: $`T<T_m`$,$`f=f_c(T)`$: At criticality, one expects that the cost for unzipping and zipping is the same (the equilibrium probability distribution of having $`m`$ monomers unzipped or an end-to-end distance equal to $`x`$ is flat), therefore $`x(t+t_0)=x(t)\pm \frac{t_0}{x}`$ with equal probability. Therefore, the end-to-end distance makes a random walk in the rescaled time $`\frac{t}{x}`$ so that $`x(\frac{t}{x})^{\frac{1}{2}}`$ and $`\theta _2(=\theta _1)=\frac{1}{3}`$ and $`d_1=d_2=1`$ since at coexistence there is a finite fraction of liberated monomers.

Another way of obtaining $`\theta _1=\frac{1}{3}`$ is to demand that a kink liberated at the fork needs to diffuse out of the end before the next one is released. In other words, the rate of change of $`m`$ is determined by the diffusion of a kink over a distance $`m`$. The latter time-scale being of order $`m^2`$, we expect $`dm/dtm^2`$ which gives $`mt^{\frac{1}{3}}`$.

This simple model deserves some observations. Remarkably, it displays two different time scales in its dynamics, mirrored in the difference of the exponents $`z_{1,2}`$ in regimes A and C. The first time scale ($`tN^{z_1}`$) quantifies the time necessary for the unbinding (or unzipping) of the bases while the second ($`tN^{z_2}`$) gives the time needed to open (and to stretch whenever $`f0`$) the two chains up to their equilibrium end-to-end distance. At $`T<T_m`$ the two processes are virtually the same, because the unbinding (or unzipping) is dragged by the stretching. However, above $`T_m`$, the processes decouple and the unbinding gets faster, being controlled by the temperature, and therefore $`z_1<z_2`$.

Furthermore, in the numerical calculations we found large sample-to-sample fluctuations, thereby requiring a huge number of runs to get good averages of $`m(t)`$ and $`x(t)`$. This is due to the long time correlation that exists in the system, which keeps samples with different initial histories far apart for any $`t`$.

Turning to the b-model with bubbles, amazingly the dynamical exponents in regimes B, E and D with $`TT_m`$ are the same as found in the Y-model. This establishes that at $`T<T_m`$ not only for statics, as we saw in the previous section, but also for the dynamics, bubbles are not relevant in the scaling properties (for $`T=T_m`$ the equality is due to another reason ). At $`TT_m`$, on the other hand, the opening of bubbles heavily affects the base unpairing process, unlike the Y-model case where bubbles are forbidden. The length of the unzipped part in the present case now can change by $`\pm l(t)`$, where $`l(t)`$ is the typical length of bubbles, and the motion of the fork point can by no means be approximated by a simple random walk ( and so $`\theta _1`$ changes as shown in Table I). The quantity $`x(t)`$ instead has a dynamics in the b-model similar to the fork case and indeed $`\theta _2`$ is the same for both the models in all regimes. We show in Fig.3 the collapse leading to $`\theta _1=\frac{3}{4}`$ in regime $`A`$.

An important question is that of the dependence of these results on dimensionality. As for the Y-model, the arguments we gave suggest that there be no $`d`$-dependence. For the b-model, instead, this should be true only in regimes D and E for $`T<T_m`$, where the Y-model gives the exact result; at $`TT_m`$, on the other hand, bubbles play a dominant role and so we expect a dependence on $`d`$. We confirmed this picture with some calculations on a simpler model which should be in the same universality class of the one under study: that of a single random walk, pinned at the origin by an attractive interaction and subject to a stretching external force. In this system $`m`$ is defined as the number of monomers from the last visit to the origin to the end of the walk. For $`T>T_m`$, in regimes $`A`$ and $`C`$, our calculations show that the exponent $`\theta _1`$ increases as dimension increases, apparently with no upper critical dimension. Just at criticality at zero force, instead, we find that the exponent $`\theta _1`$ is very close to $`\frac{1}{2}`$ in any $`d`$. The emerging picture of robust results for $`T<T_m`$ and model-dependent dynamics for $`TT_m`$ would be preserved even if, in the original models, the directedness constraint is relaxed (as is the case for the statics).

The arguments presented so far could be generalized to add other ingredients of dynamics also. An example is the nonlocal effects in dynamics, as is the case e.g. in the Zimm model. In ref., the author suggests that, in the regime we call $`D`$, in the mean field approximation, non-locality could speed up the opening so that $`\theta _{1,2}`$ would be $`\frac{1}{1+u}`$ instead of $`\frac{1}{2}`$ ($`u`$ is the exponent characterizing the length dependence of the mobility as defined, e.g., for the Zimm model, see  for the notation and  for a review). In regime $`E`$, we can combine our arguments with the same reasoning to get $`\theta _1=\theta _2=\frac{1}{2+u}`$. As we expect, the Rouse model results are obtained with $`u=1`$.

Lastly, the above analysis can also be extended to binary disorder in the contact potential, which is of course a realistic feature to be included in the model, i.e. the energy of the contact of the $`i`$-th base may be $`ϵ_i=ϵ+\mathrm{\Delta }_i`$, $`\mathrm{\Delta }_i`$ being a random variable with binary distribution and zero mean. The Y-model offers a good starting point for the study of the effects of heterogeneity, because the critical line of the quenched model can be proved to be the same as that of the pure model with energy $`ϵ`$. If we call $`(m)`$ the free energy density of a configuration with $`m`$ (out of $`N`$) unzipped monomers, we expect on general grounds that $`P(m,t)`$, the probability of having $`m`$ monomers unzipped at time $`t`$, will obey a master equation with transition rates $`W_{\pm ,m}`$ depending on the realization of disorder, i.e. we expect that $`W_{\pm ,m}\mathrm{min}\{e^{\beta N\frac{d(m)}{dm}},1\}`$, with $`N\beta \frac{d(m)}{dm}`$ $`=+\frac{\mathrm{\Delta }_{Nm}}{T}+\frac{ϵ}{T}\mathrm{log}(1+\mathrm{cosh}(\frac{f}{T}))`$ , which can be seen to consist of a zero-mean random “noise” ($`\frac{\mathrm{\Delta }_{Nm}}{T}`$) plus a “bias”. We thus believe that the dynamics of the Y-model can be mapped onto Temkin’s model of a random walker in a random environment  (provided that at $`T<T_m`$ we rescale time as described above in regimes $`D`$ and $`E`$ ), so that, following Ref. , there will be a region around the critical line, at $`T<T_m`$, in which $`m,xt^{\frac{1}{2}}`$ (and just on the critical line at $`f0`$ one is tempted to expect $`m,x\left(\mathrm{log}\left(\frac{t}{x}\right)\right)^2`$). The quantities $`m,x`$ are now quenched averages over realizations of disorder. The scenario we propose is sketched in Fig.4. The curve bounding the region where, at $`T<T_m`$, disorder should be important for the dynamics has been found by applying the criterion of Ref.  to the system with the above transition rates. This curve coincides with the critical line of the static model with annealed disordered. Preliminary runs at very low $`T`$ suggest that as $`T0`$, for a given realization of disorder and with any force between $`ϵ_1`$ and $`ϵ_2`$ (the two binding energies), the unzipping will take place as far as the first more attractive base pair is found. Besides, the probability that the monomer of index $`m_0`$ will ever be unzipped decays exponentially: $`P(m_0)2^{m_0}`$. This is consistent with the mapping suggested above. A detailed study of the effect of disorder will be done elsewhere.

Conclusion. We obtained the phase diagram of a lattice model for the unzipping transition of a double stranded DNA by a force. This model incorporates both the mutual and self-avoiding nature of the two strands and the equilibrium problem can be solved exactly. A still simpler model has been defined by suppressing the bubbles that are important for thermal denaturation. This simpler version is also exactly solvable and is shown to retain the basic features of the b-model. Thanks to the exact knowledge of the phase boundary, we have been able to investigate the dynamics of unzipping from a non-equilibrium bound state both on and away from the phase boundary. The dynamics shows scaling behaviors in different regimes of the phase diagram. These scalings in most cases could be understood from the plausible mechanisms of unzipping and denaturation as discussed in the previous sections. We end with two notable features. In the case of the b-model with bubbles, the unzipping at high temperatures remains a puzzling issue, especially the dynamic exponent $`z_1=\frac{4}{3}`$ with or without force. We believe that this has to do with the statistics of bubbles, though no satisfactory answer could be found. This remains an open problem. Lastly, on a positive note, the unzipping dynamics on the phase boundary in the presence of a force is found to be distinctly different from the thermal denaturation at zero force. Whether a real biological system takes advantage of these differences to distinguish the unzipped region of DNA from a fluctuation-induced bubble formation remains to be probed.

Acknowledgments: This work was supported by MURST(COFIN 99). DM also acknowledges INFM funding. We thank H. Dobbs and A. Trovato for discussions.


text:
# Introduction

## Introduction

Understanding the complex relaxation behavior in supercooled liquids has been a field of much research interest in recent times. In this regard the response of a system to an energy fluctuation namely the frequency dependence of specific heat of a supercooled liquid has been investigated by a number of authors. Generally, specific heat is a property that is usually linked to the thermodynamic property of a system . The pioneering experiments done by Birge and Nagel  have studied the dynamic response in glassy systems namely glycerol and propylene and obtained interesting dynamical response behavior, expressed in terms of a frequency dependent specific heat . In an experiment usually the frequency dependent product of thermal conductivity($`\kappa `$) and specific heat, $`\kappa c_p(\omega )`$ is measured. However, in the temperature range $`190^oK220^oK`$, over which we are interested here, $`\kappa `$ has weak frequency dependence, and hence the dynamics observed in the product $`\kappa c_p(\omega )`$ is entirely due to the frequency dependence of the specific heat. The theoretical modeling for the frequency dependence of specific heat in a supercooled liquid has been studied by various authors ,. New internal mode was proposed to be present in the supercooled liquid to explain the frequency dependent response. In a simple analysis, Zwangig however had argued  that the frequency dependence of the specific heat can be obtained without introducing any such internal mode. This work showed that what is measured as the frequency dependent specific heat is actually related to that of the longitudinal viscosity in the liquid. In this model the dynamics of fluctuations around the equilibrium was studied in terms of a simple set of slow variables of Hydrodynamics of fluids. These equations of motion used were the conservation laws of mass, momentum and energy in the system. The resulting formula for the specific heat is equivalent to linking of structural relaxation in a supercooled liquid to the frequency dependence of the specific heat. In the present work we take the data from the Specific heat measurement  and extract the frequency dependence of the viscosity as will be required from such a formulation proposed in Ref. . We then address the question if this value of the longitudinal viscosity will indeed be self consistent with independent measurements on the structural properties. The basic idea is to consider the frequency dependence of the specific heat solely in terms of the structural relaxation in the supercooled liquid.

For structural relaxation \- the microscopic model for the liquid dynamics namely the Mode Coupling Theory (MCT) has been studied by a number of authors in recent years . A simple application  of the MCT to fit the specific heat data indicated a very large exponent is required to match the experimental data to power law divergence. In the present work we also use the MCT as a model for structural relaxation and obtain the corresponding frequency dependence of the longitudinal viscosity. We then use it to predict the behavior of the specific heat with frequency as predicted from the theory proposed in Ref . In the microscopic model of the Mode coupling Theory the wave number dependence of the longitudinal viscosity is obtained using proper input for the structure factor of the liquid. For this purpose standard results from the integral equations for simple liquids are used for the structure factor. The frequency dependent specific heat is then computed for different wave numbers. Thus the effect of the response to heat fluctuations can be computed over different length and time scales in the present approach. While this extends the theory with scope of further comparison, the main goal of the present work is to test if the frequency dependence of the specific heat can be understood solely in terms of the structural relaxation and if the two sets of measurements agree in a self-consistent manner. The paper is organized as follows. In the next section, Sec. II, we consider the schematic model for the time dependence of the viscosity and, in Sec. III, we compare the theoretical results with the experimental observations. A wavenumber dependent calculation for the specific heat is presented in the Sec.VI. In Sec. V, we present the Mode Coupling results for the specific heat over different length-scales and temperatures. In the last section we discuss the results.

## II. Frequency Dependent Specific Heat

Since we are concerned here with the dynamic properties of a supercooled liquid, an obvious choice is to consider a hydrodynamic model for the system. To start with, we write down the linearized hydrodynamic equations,

$$\frac{}{t}\delta \rho +.\stackrel{}{g}=0$$
(1)
$$\frac{g_i}{t}+_iP\eta ^2\stackrel{}{g}_i(\frac{1}{3}\eta +\zeta )(.\stackrel{}{g}_i)=0$$
(2)
$$\rho _oc_v\frac{}{t}\delta T+\kappa ^2T+\frac{T_o}{\rho _o}\left(\frac{P}{T}\right)_\rho .\stackrel{}{g}_i=0,$$
(3)
governing the time evolution of fluctuations of conserved variables- mass density($`\rho `$), momentum density($`\stackrel{}{g}`$) and the temperature(energy)$`T`$. Here $`\rho _o`$ and $`T_o`$ represent equilibrium density and temperature respectively and $`\delta \rho `$ and $`\delta T`$ are the fluctuations from the equilibrium values. $`c_v`$ is the specific heat per unit mass at constant volume, $`\eta `$ and $`\zeta `$ are the shear and bulk viscosities respectively. The viscosity coefficients here are divided by the density. The fluctuation of the pressure $`P`$ around the equilibrium value can be expanded to the lowest order in density and temperature as,

$$\delta P=\left(\frac{P}{\rho }\right)_T\delta \rho +\left(\frac{P}{T}\right)_\rho \delta T,$$
(4)
where we have assumed that the change of the Pressure functional with the density function at the equilibrium can be replaced by the equilibrium partial derivative - replaced by the corresponding thermodynamic quantity. Using the above equations the energy conservation equation (3) reduces to the Fourier heat law for thermal fluctuations,

$$\iota \omega \delta T=\mu (\omega )^2\delta T,$$
(5)
with the frequency dependent thermal diffusivity $`\mu (\omega )`$ defined in terms of the specific heat $`c_p`$ as,

$$\mu (\omega )=\frac{\kappa }{\rho _oc_p(\omega )}.$$
(6)
The specific heat $`c_p(\omega )`$ is expressed in the form,

$$c_p(\omega )=c_v+(c_pc_v)\frac{K_T(0)}{K_T(\omega )}.$$
(7)
The quantity $`K_T(\omega )`$ is called the generalized bulk modulus and is given by

$$\frac{K_T(\omega )}{K_0}=1+\iota \omega \mathrm{\Gamma }(\omega ),$$
(8)
and is expressed in terms of the reduced form $`\mathrm{\Gamma }(\omega )=\eta _l(\omega )/c_o^2`$ of the frequency dependent longitudinal viscosity $`\eta _l(\omega )`$. In equation (8), $`K_0`$ is the $`\omega =0`$ limit of $`K_T(\omega )`$. Obviously for the liquid state with finite viscosity the zero frequency limit of $`K_T(\omega )`$ relates to the thermodynamic property of the supercooled liquid. The sound speed $`c_o`$ is given by

$$c_{o}^{}{}_{}{}^{2}=\left(\frac{\delta P}{\delta \rho }\right)_T=\frac{K_0}{\rho _o}.$$
(9)
A frequency dependent longitudinal modulus $`M(\omega )`$, the inverse of compliance, is defined along a similar line as,

$$\frac{M(\omega )}{K_0}=\gamma +\iota \omega \mathrm{\Gamma }(\omega ),$$
(10)
where $`\gamma `$ is the ratio of the long time limit of the specific heat $`c_p(\omega =0)`$ to $`c_v`$. Equation (7) is the key formula used in this paper for testing the idea of modeling the frequency dependence in the specific heat solely in terms of the structural relaxation. In obtaining equation (5) one also needs to assume that the following self-consistent relation holds,

$$\mathrm{\Delta }(\omega )(\gamma 1)\frac{\omega }{\overline{M}(\omega )[\omega +i\nu \overline{M}(\omega )]}<<1.$$
(11)
Here we have expressed $`M`$ in the dimensionless form as, $`\overline{M}(\omega )=M(\omega )/K_0`$. $`\nu =c_0^2/\mu _o`$, with $`\mu _o=\kappa /(\rho _oc_v)`$ is the bare thermal diffusivity. We test the validity of the assumption (11) in the frequency range where the analysis with respect to experimental data is made.

## III. Comparison of Experimental Data

In this section we test the self-consistency in expressing the frequency dependent data on specific heat and Structural relaxation. For the supercooled liquids the relaxation over longest time scales, i.e. the $`\alpha `$-relaxations, follows the stretched exponential behavior,

$$\eta (t)=\eta _oexp\left[\left(\frac{t}{\tau }\right)^\beta \right]$$
(12)
where $`\eta _o`$ is the amplitude and $`\beta `$ is the stretching parameter which defines the degree of deviation from the exponential decay. In fitting the specific heat data we use the dimensionless form for the specific heat ratio,

$$c_p(\omega )=c_v\left[1+\frac{(\gamma 1)}{1+\iota \omega \mathrm{\Gamma }(\omega )}\right],$$
(13)
which reduces the formula in a dimensionless form. We fit the specific heat data of Ref.  to the formula (13) using a simple stretched exponential (12) relaxation function. In Fig. 1(a) and 1(b) we show the respective fitting of the experimental data for the real and imaginary parts of $`c_p(\omega )`$ in the supercooled glycerol for three different temperatures, T = 201.4$`{}_{}{}^{o}K`$, 203.9$`{}_{}{}^{o}K`$, 211.4$`{}_{}{}^{o}K`$. The arrows in Fig. 1(b) indicate the peak positions in the imaginary parts of the viscosities at the corresponding temperatures. In calculating the specific heat, the three parameters $`\mathrm{\Gamma }(t=0))`$, the relaxation time $`\tau `$ and the stretching exponent $`\beta `$ are used as the free parameters. $`\gamma `$=1.86. Using the best fit values of the parameters with the specific heat data we compute the structural properties of the liquid given by Modulus $`M(\omega )`$ defined in eqn. (10) and the longitudinal viscosity. The resulting behavior for these quantities are compared with the experimental results as shown in Figures 2 and 3.

In Figure 2, we show the viscosity $`\eta `$ in the zero frequency limit in units of $`K_o\tau _o`$ where $`\tau _o`$ is the unit of time used. In the Inset we show the corresponding experimental data  for the viscosity. The experimental data shown here is over a much wider temperature range (317$`{}_{}{}^{o}K190^oK`$) \- both theoretical and experimental data agree with the Vogel Fulcher fit ( $`\eta exp[A/(TT_o)]`$ for $`T_o=128^o`$K and $`A=2480^o`$K shown in both the figure and the inset as solid lines. The zero frequency modulus $`K_o`$ is roughly temperature independent  over the range considered here. We notice that the viscosity increases by four orders of magnitude as the temperature is decreased over a small range($`200^oK220^oK`$) near the glass transition temperature($`T_g=`$190$`{}_{}{}^{o}K`$). We define a normalized longitudinal modulus

$$\stackrel{~}{M}(\omega )=\frac{M(\omega )M(0)}{M(\mathrm{})M(0)}.$$
(14)
In Fig. 3 (a) and (b) we show the real and imaginary parts of $`\stackrel{~}{M}`$ respectively denoted by $`M^{}(\omega )`$ and imaginary $`M^{\prime \prime }(\omega )`$ against the frequency for the three different temperatures, T=203.9$`{}_{}{}^{o}K`$(solid), 211.4$`{}_{}{}^{o}K`$ (dashed), 221.5$`{}_{}{}^{o}K`$ (dotted) lines. The frequency in each case is expressed in terms of the ratio with the corresponding peak position $`(\omega _p)`$ in the imaginary part. The corresponding results from measurements on the modulus $`\stackrel{~}{M}`$  are also shown with filled circle. In Fig. 4 we show the plot of the peak positions as found in the fitting with different temperatures. As the temperature is decreased, peak in the imaginary part of the specific heat shifts towards the lower temperatures, signifying the slower relaxations in the system. The solid line indicate V-F fit with $`T_o=128^o`$K. In Fig 5 we show the frequency dependent Specific heat and the viscosity function ( in the inset ) at the same temperature. The peaks appear nearly at the same position on the frequency scale for the two quantities showing that the dominant time scales are same in the two cases. Finally we test the validity of the assumption (11) that is crucial in reaching the Fourier heat law (5) with the frequency dependent specific heat - defined above. For this we calculate both the real$`(\mathrm{\Delta }^{}(\omega ))`$ and imaginary $`(\mathrm{\Delta }^{\prime \prime }(\omega ))`$ parts of $`(\mathrm{\Delta }(\omega ))`$ for the supercooled glycerol. In Fig. 6 we plot both the real and imaginary parts of $`\mathrm{\Delta }(\omega )`$ on the frequency range over which specific heat(frequency dependent) is observed. These figures clearly show that the quantity $`\mathrm{\Delta }(\omega )`$ is much smaller as compared to unity over this frequency range. This substantiates the assumption made in the previous section to reach the Fourier heat law in a generalized sense.

## IV. Wavevector Dependence of Specific Heat

In the previous section, we studied the specific heat and other quantities like longitudinal viscosity and modulus using a schematic model to show the self-consistency of the relation between structural relaxation and the frequency dependence of specific heat. Here we consider the wavevector and frequency dependent specific heat in a liquid. Starting from the generalized hydrodynamic equations for the conserved densities in q-space, we obtain an equation,

$$\iota \omega \rho _oc_v\delta T(q,\omega )=q^2\kappa \delta T(q,\omega )+\frac{\iota q^2\omega T_o}{(\rho _o\omega ^2q^2K_T(q,\omega )}\left(\frac{P}{T}\right)_\rho ^2\delta T(q,\omega )$$
(15)
which describes dynamics of the energy fluctuations over different length and time scales. $`K_T(q,\omega )`$ is the wave vector and frequency dependent bulk modulus given by,

$$\frac{K_T(q,\omega )}{K_T(q)}=1+\iota \omega \mathrm{\Gamma }(q,\omega )$$
(16)
where $`\mathrm{\Gamma }(q,\omega )`$ is the wavevector and frequency dependent longitudinal viscosity devided by the square of the speed of sound $`c_s^2(q)=K_T(q)/\rho _o`$. $`K_T(q)`$ is the zero frequency limit of $`K_T(q,\omega )`$. The energy equation (15) reduces to the wavevector dependent Fourier heat equation,

$$\iota \omega \delta T(q,\omega )=q^2\chi (q,\omega )\delta T(q,\omega )$$
(17)
where $`\chi (q,\omega )=\kappa /(\rho _oc_p(q,\omega ))`$ is thermal diffusivity and $`c_p(q,\omega )`$ is the $`q`$-dependent specific heat given by,

$$c_p(q,\omega )=c_v\left[1+(\gamma _q1)\frac{1}{1+\iota \omega \mathrm{\Gamma }(q,\omega )}\right]$$
(18)
and $`\gamma _q`$ is the ratio $`c_p(q)/c_v`$. Here in obtaining the Fourier heat equation (17), we have assumed that the quantity,

$$(\gamma _q1)\frac{\omega }{\overline{M}(q,\omega )[\omega +\iota \nu (q)\overline{M}(q,\omega )}<<1$$
(19)
where $`\nu (q)=c_s^2(q)/\mu _o`$ and $`\overline{M}(q,\omega )=\gamma (q)+\iota \omega \mathrm{\Gamma }(q,\omega )`$. is much smaller as compared to unity. In the long wavelength limit this quantity reduces to $`\mathrm{\Delta }(\omega )`$ given by Eq. (11).

## V. Results from the Mode coupling Theory

In this section, we consider the time dependent longitudinal viscosity obtained from a microscopic theory of Statistical Mechanics, instead of taking inputs from experimental results as was done in the section III. We predict structural aspects i.e. the wave vector dependence in the specific heat at different frequencies.

In the simplest form the self-consistent mode coupling theory predicts a sharp transition of the supercooled liquid to nonergodic phase. In later versions it was shown that due to coupling of density fluctuations with currents this sharp transition is eliminated - the full model with the cutoff mechanism included is termed as the extended model. In the work we consider the extended model where the cutoff function is adjusted to obtain agreement with the viscosity of the supercooled Liquid to the results obtained from Computer simulations. The details of the model and the scheme for computation of the density correlation function using the proper cutoff function is presented elsewhere. We consider a one component Lennard-Jones system for computing the structural relaxation properties using the MCT. The temperature $`T^{}`$ and density $`\rho ^{}`$ are expressed in the standard units of $`ϵ/K_B`$ and $`\sigma ^3`$ respectively. $`ϵ`$ is the unit of energy in a L-J system and $`\sigma `$ is the diameter of a particle.

In computing the dynamical behavior of the density correlation function we estimate the cutoff parameters of the theory so that the shear viscosity obtained from the self consistent MCT agrees with computer simulation results. For the simulation results on one component model we use the recent results of Rucco et. al.  using special techniques that avoid the typical problem of crystallization in one component systems. From the self consistent results for the density correlation functions we compute the mode coupling integrals for the longitudinal part of the memory function related to the decay of the density correlation functions.

The longitudinal viscosity in the zero wave number limit is shown in Fig. 7 for the temperature range around $`T_c`$. The longitudinal viscosity shown for the temperature range around $`T_c/T`$ less than 1, follows the power law behavior. The viscosity diverges with exponent equal to 1.9 around $`T_c`$ and for lower temperatures ( $`T_c/T`$ 1 ) the behavior follows a Vogel Fulcher form. This is usual with results of extended MCT . We then use this frequency dependent memory function or the longitudinal viscosity to compute the corresponding frequency dependent specific heat. In the Mode-Coupling approximation, the normalized longitudinal viscosity, $`\mathrm{\Gamma }(q,\omega )`$, is given by,

$$\mathrm{\Gamma }(q,\omega )=\frac{1}{2n}S(q)\frac{d\stackrel{}{k}}{(2\pi )^3}[\frac{\widehat{q}.\stackrel{}{k}}{q}c(k)+\frac{\widehat{q}.(\stackrel{}{q}\stackrel{}{k})}{q}c(|\stackrel{}{q}\stackrel{}{k}|)]^2\psi (k,t)\psi (|\stackrel{}{q}\stackrel{}{k},t)$$
(20)
where $`\psi (q,t)`$ and $`S(q)`$ represent density correlation function and the structure factor respectively. $`c(q)`$ is the direct correlation function of the system. $`\widehat{q}`$ denotes the unit vector along $`\stackrel{}{q}`$ and $`n`$ is the number density. Using the above expression for the longitudinal viscosity in Eq. (18), we calculate the $`q`$-dependent specific heat in the supercooled liquid.

In figure 8(a) and (b) we show respectively the real and imaginary part of $`c_p^{}(q,\omega )`$ vs. the frequency. This is shown here for three different wave vectors, $`q\sigma `$=0, 7.05(peak of the structure factor) and 30 (upper cutoff taken for the k-integral) at temperature $`T^{}`$=.559. The inset of the corresponding figures shows the secondary peak predicted for fast processes at very high frequency window - representing the so called $`\beta `$-processes in the supercooled liquid. The peak in the imaginary part shifts to lower frequency with lowering of temperature. In Fig. 9(a) we show the variation of the peak position with temperature of the liquid. The solid line in the figure shows Vogel-Fulcher fit with $`T_o`$=.014. In order to indicate the structural dependence we also show in Fig. 9(b), the dependence of the peak position on the wave number $`q`$ at a fixed temperature $`T^{}`$=.559 . The peak frequency signifying the dominant time scale for relaxation at different wavenumbers follows the nature of the structure factor. It shows a minimum at $`q`$-value which corresponds to the peak in the structure factor of the liquid. Successive minima in the figure correspond to the other less pronounced maxima in the structure factor.

## VI. Discussion

In section III, we considered two types of experimental measurements on the supercooled liquids, respectively related to the energy fluctuations and the structural relaxations. This was done to check the consistency of formula for specific heat obtained from simple analysis of the Hydrodynamic equations. The comparisons done in section III indicates that the frequency dependent specific heat can be understood in terms of the structural relaxation data in terms of the analysis proposed in Ref. .

Subsequently we apply the standard forms of the self-consistent Mode Coupling Theory to compute the frequency dependence of the viscosity and compute those for the specific heat. Here we use the formula in the terms of the Generalized Hydrodynamics, extending the model to large $`k`$, or small wavelengths. We find that the dispersion in the specific heat decreases as we go to the smaller length scales(higher $`q^{}`$) with corresponding increase in spectrum width. This demonstrates the fact that at short length scales, the relaxation is fast and here the memory effects reflecting cooperativity are not strong.

The peak position ($`\omega _p`$) in the imaginary part of the $`c_p^{}(q,\omega )`$ shifts to higher frequency with the increase of the corresponding wavevector $`q^{}`$. It however reaches to a minimum frequency at the structure factor peak. Finally since the MCT relates to the two step relaxation process in supercooled liquids, there is a corresponding implication on the specific heat curve predicting a peak at very high frequency in the specific heat. This is shown in the inset of figure 8(b). It is a consequence of secondary relaxation in the supercooled liquid. Due to the constraints on the MCT at very low temperatures, we could not study the thermal response of the system close to the glass transition temperature $`T_g`$. However, as is shown in the Fig. 9(a), the main peak in the specific heat moves towards the smaller frequencies with decreasing temperature, thus at the temperatures very close to the glass transition one can expect the two peaks to lie further apart from each other. We have ignored here effects of nonlinearities in the energy equation . This can produce frequency dependence on other transport coefficients like thermal conductivity as well. However, observation of such a peak will further strengthen the validity of the simple analysis presented here in energy transport in terms of structural relaxation behavior.

## Acknowledgement

We are greatful to Prof. N. O. Berge for providing data on specific heat measurements. UH acknowledges the financial support from the University Grant Commission, India.

## Figure Captions

Fig. 1(a): Theory fit to the experimental data(filled circles) for the real part $`c_p(\omega )`$ of the specific heat ($`c_p(\omega )`$) at three temperatures $`T`$= 201.4$`{}_{}{}^{o}K`$(continuous line), 203.9$`{}_{}{}^{o}K`$(long dashed) and 211.4$`{}_{}{}^{o}K`$(dotted). $`\omega ^{}=\omega \tau _o`$, where $`\tau _o`$ is the units of time used(see text).

Fig. 1(b): imaginary part $`c^_p(\omega )`$ of specific heat corresponding to the real part shown in Fig. 1(a). arrows along the frequency axis indicate the peak position in the imaginary part of the corresponding viscosity.

Fig. 2: temperature variation of viscosity in supercooled glycerol as obtained from the specific heat fitting. Viscosity is given in units of $`K_o\tau _o`$. experimental results for the viscosity are shown in the inset. The continuous line( both in the main figure and the inset) is the Vogal-Fulcher, $`\eta =\eta _oexp[A/(TT_c)]`$, fit with $`\eta _o`$=1.08$`\times 10^4P`$, $`A`$= 2480$`{}_{}{}^{o}K`$ and $`T_o`$=128$`{}_{}{}^{o}K`$.

Fig. 3(a): real part, $`M^{}(\omega )`$, of the normalized longitudinal modulus $`\stackrel{~}{M}(\omega )`$(see text) is plotted at three different temperatures T=203.9$`{}_{}{}^{o}K`$(solid line), 211.4$`{}_{}{}^{o}K`$(dashed line) and 221.5$`{}_{}{}^{o}K`$(dotted line). dots are the experimental results of ref. . frequency axis is scaled with respect to the peak values $`\omega _p`$ for the three different temperatures.

Fig. 3(b): Imaginary parts of the normalized longitudinal modulus $`\stackrel{~}{M}(\omega )`$ corresponding to the real parts shown in the figure 3(a).

Fig. 4: peak position $`(\omega _p)`$ in the imaginary part of the specific heat is plotted as a function of the temperature$`(T)`$. continuous line is the VF fit :$`\omega _p=\omega _oexp[A/(TT_c)]`$; with $`\omega _o=1.0\times 10^{15}`$ Hz, $`E`$=2559.35$`{}_{}{}^{o}K`$ and $`T_o=`$ 128.22$`{}_{}{}^{o}K`$.

Fig. 5: Imaginary part of the specific heat, $`c_p^{\prime \prime }(\omega )`$ at T=214$`{}_{}{}^{o}K`$. Arrow along the frequency axis at 2.45 indicates the peak position in the imaginary part of the corresponding viscosity shown in the inset.

Fig. 6 real and imaginary parts of $`\mathrm{\Delta }(\omega )`$(see text), for supercooled glycerol, is plotted at $`T=201.4K`$.

Fig. 7: mode-coupling viscosity is plotted as a function of temperature. it shows a vogul-Fulcher behaviour, $`\eta _oexp[A/(TT_o)]`$ for $`T<T_c`$ with $`T_o`$=.023 while for higher temperatures$`(T>T_c)`$ it follows a power low behaviour with exponent 1.9 . Arrow along the temperature axis at $`T_C/T`$=1.36 indicates the power-low divergence.

Fig. 8(a): MCT results for the real part of the normalized specific heat $`c_p^{}(q,\omega )`$ = $`(c_p(q,\omega )c_v)/(c_p(q)c_v)`$ for three wave vectors $`q`$=0 (dotted), 7.05(continuous) and 30(dashed) at $`T^{}`$=.559. Here frequency $`\omega ^{}`$ is in the units of the inverse of Lenard-Jones time, $`\tau =\sqrt{\frac{m\sigma ^2}{ϵ}}`$.

Fig. 8(b): imaginary part of $`\stackrel{~}{c}_p(q,\omega )`$ from the MCT corresponding to the real parts shown in Fig 8(a). Inset shows the secondary peaks predicted by the MCT for the same three $`q`$-values.

Fig. 9(a): variation of the peak position ($`\omega _p`$) in the specific heat $`\stackrel{~}{c}_p(q,\omega )`$ with temperature for $`q`$=0. along the y-axis, we have shown $`\omega _p^{}=\omega \tau \times 10^2`$.

Fig. 9(b): variation of the peak frequency $`(\omega _p)`$ with wave vector $`q^{}`$. $`\omega _p`$ reaches to minimum at $`q^{}`$=7.05 at which the structure factor shows a maximum. $`\omega _p^{}=\omega \tau \times 10^2`$.