text:
# The Balmer Wavelength Range of BP Tauri

## 1 Introduction

The International Ultraviolet Explorer (IUE) was operational for almost 20 years and took over 107000 images of a wide array of objects. In particular, it observed about 130 T-Tauri stars (TTSs) in the low resolution ($`\lambda 6\mathrm{\AA }`$), long wavelength range (from $`\lambda =1850\mathrm{\AA }`$ to $`\lambda =3350\mathrm{\AA }`$). This dataset constitutes the most complete set of observations of TTSs in the Balmer range.

Now that the IUE Final Archive has been completed and all the data is in the public domain111http://archive.stsci.edu/iue/ we have undertaken the project of exploring all this dataset. In this paper we will examine the observations of BP Tauri (also known as HBC 32 and HD 281934). A complete discussion of the full dataset will appear elsewhere.

We have chosen this star in particular because of the large number of observations by IUE: there are 61 low resolution spectra. They will allow us to study correlations between different physical aspects of the T-Tauri phenomenon over a long period of time.

BP Tauri is a single (Ghez, Neugebauer & Matthews 1993), Classical T-Tauri star (CTTS) located in L1495 (Lynds 1962) in the Taurus-Auriga complex of dark clouds at a distance of 140 pc (Elias 1978). The star has a constant radial velocity close to that of its molecular surroundings and a projected equatorial velocity of $`v\mathrm{sin}i<10`$ km/s (Hartmann et al. 1986; Basri & Batalha 1990). The spectral type of BP Tauri corresponds to a K5-K7 dwarf with a strong absorption in the lithium resonance line (Duncan 1991; Hartigan et al. 1989). The mass is $`M=0.50.8M_{}`$ (Gulbring 1994; Gullbring et al. 1998).

Photometric observations (Vrba et al. 1986; Bouvier, Bertout & Bouchet 1988; Simon, Vrba & Herbst 1990 - SVH90 - and the observations summarized by Rydgren et al. 1984) have shown that BP Tauri has a large infrared emission and is variable in all optical bands, with a range in B of up to 3 magnitudes (Herbig & Bell 1988). Rapid fluctuations on time-scales of minutes have been observed by Schneeberger, Worden & Africano (1979).

These fluctuations make it difficult to extract a possible photometric period. No period at all (Bouvier et al. 1988, Simon et al. 1990) and periods of 6.1 days (Simon et al. 1990), 7.6 days (Vrba et al. 1986), 7.7 days (Osterloh, Thommes & Kania 1996), 8.3 days (Richter et al. 1992), and $`10`$ days (Gómez de Castro & Franqueira 1997 - GdCF97) have been proposed. The fact that BP Tauri grows redder as it becomes fainter has been used to explain this variability in the period in terms of a hot spot on the surface of the star (see for example SVH90). In a small fraction of CTTSs, spots cooler than the photosphere have also been observed, but no evidence of cool spots exists for BP Tauri (Gullbring et al. 1996). Variations in the measured period of BP Tauri might correspond to the hot spot (or spots) moving to different latitudes in a differentially rotating star. An inclination of $`i30^o50^o`$ for BP Tauri has been derived on the basis of the measured $`v\mathrm{sin}i`$ and the longest measured periods (SVH90; Gullbring 1994). The variability in the photometric period can also be explained by considering a beat frequency between the Keplerian frequency of clumps of material circling the star and the rotation rate of the stellar surface (Bouvier et al. 1999; Smith, Bonnell & Lewis 1995; Smith, Lewis & Bonnell 1995).

Measurements of spot sizes in this and other TTSs give values that vary between a fraction of a percent to almost 40% the visible stellar surface (Vrba et al. 1986; Vrba et al. 1993; Herbst, Herbst & Grossman 1994; Gullbring 1994; Bouvier et al. 1995; Fernández & Eiroa 1996).

As has been observed in many CTTSs, the spectrum of BP Tauri appears veiled by variable continuum excess emission. This is the origin of the photometric variability. The optical veiling has been observed to be between 30 and 75% the value of the photospheric flux (Basri & Batalha 1990; Valenti, Basri & Johns 1993; Hartigan, Edwards & Ghandour 1995; Gullbring et al. 1998).

Some spectral lines (notably, the Balmer series and the Mg II line) show evidence of blueshifted and redshifted absorption (see for example, Johns & Basri 1995ab; GdCF97). It has also been noted that line emission and the continuum level are sometimes correlated (GdCF97; SVH90). We will explore one of these correlations in Section 3.

These observations can be framed in terms of the current paradigm for CTTS. As we understand it today, a CTTS consists of a late-type pre-main-sequence star with a magnetic field, surrounded by an accreting, dusty disk. At the distance from the star at which the torque due to the magnetic field equals the viscous torque in the disk, the disk is truncated. Some material from the disk will then be captured by the magnetic field (see for example Königl 1991) and some will escape as a wind (Shu et al. 1994). Collisional coupling between ions and neutrals is such that even the neutral gas component follows the magnetic field lines (Martin, 1996). At the stellar surface, the material will have supersonic velocity. The simplest models assume that this velocity should be the free-fall velocity ($`v_{ff}300\mathrm{km}\mathrm{s}^1`$, Calvet & Gullbring 1998) but more complicated models slow down the material to a fraction ot this velocity (Ostriker & Shu 1995). Nevertheless, a shock will occur some distance above the star. The shocked gas will heat up the surface of the star and a hot spot will be seen. Emission from this hot spot is thought to be responsible for the visible and ultraviolet excess observed in CTTS spectra (Valenti et al. 1993; Calvet & Gullbring 1998; Lamzin 1998). Furthermore, emission from the whole accretion flow is believed to be responsible for other features of the TTS spectra, in particular, the strong emission in the Balmer lines (Hartmann, Hewett & Calvet 1994; Johns-Krull & Basri 1997) and some atomic lines (Batalha et al. 1996; Beristain, Edwards & Kwan 1998). In this model redshifted absorption features would correspond to material falling into the accretion spot. Blueshifted absorption features would correspond to absorption by an optically thin wind. It is unclear at this point what fraction of the excess emission (if any) is due to the wind.

In the context of the magnetospheric accretion model, the ultraviolet is a particularly interesting range. Emission from the underlying star is low compared to that of the veiling continuum. This means that, even though is not possible to measure the veiling directly (as it is possible in the visible range), conclusions derived from the ultraviolet analysis are not very sensitive to errors in the spectral type.

In this paper we will explore the IUE dataset for BP Tauri. Relatively little work has been done on the individual IUE observations and their relation to the accretion paradigm. Here we assume the excess ultraviolet continuum emission to be caused by accretion. We start with a discussion of the raw data. A more complete discussion (including many more stars) will appear soon, in a paper dealing with the full low-resolution IUE dataset. In Section 3, we will use our dataset to study the correlation between the UV continuum and the flux in the Mg II blend. To understand the IUE data to the fullest we need a physical model. Therefore in Section 4 we will use the Calvet & Gullbring (1998) (CG98) description of the accretion shock to interpret our spectra in terms of accretion rates. Section 5 and 6 contain the analysis and the conclusions.

## 2 Observations

We searched the IUE Merged Log for observations of TTS taken with the long wavelength cameras (LWP and LWR, from $`\lambda =1850\mathrm{\AA }`$ to $`\lambda =3350\mathrm{\AA }`$) and obtained spectra for a total of 131 different TTS. For BP Tauri we found 61 low resolution ($`\lambda 6\mathrm{\AA }`$) spectra. A summary of the observations is shown in Table 1. The oldest spectra comes from April 1979 and the most recent from January 1992. There are clusters of observations taken in October 1986 (SVH90) and January 1992 (GdCF97).

All the spectra used in this paper have been processed by the New Spectra Image Processing System (NEWSIPS). For an explanation of this system and the differences with IUESIPS, see Nichols & Linksky (1996). For each observation the spectrum has been sky-subtracted by using the background in the camera. NEWSIPS spectra come with an estimate of the error at each wavelength. This error is derived empirically for each camera by measuring the scatter in the flux around the mean in the background regions of several hundred images (Garhart et al. 1997). As such, it is an estimate of the instrumental, statistical and background noise for each spectrum. NEWSIPS also provides a “quality” flag at each wavelength of each spectrum. This flag indicates certain special events, such as a saturated pixel, cosmic ray or the presence of a calibration mark.

Figure 1a. shows the mean and the normalized variance (Johns & Basri 1995b) of all spectra. The feature at $`\lambda 3050\mathrm{\AA }`$ is a reseaux mark. The central line is the blending of the Mg II h & k lines. Figure 1a shows that the edges of the spectra vary the most, whereas the center is very constant. The noise of the spectra also increases towards the edges (due to the reduced sensitivity of the Vidicon detectors used by IUE) but not enough to account for the value of the variance.

Figure 1b shows two extreme manifestations of the spectra (not corrected for reddening). The lower one (lwp18976, solid line) shows strong emission of C II\] and Fe II multiplets. For a more exhaustive description of the spectral features present, see Imhoff & Appenzeller (1987) and Brown, de M. Ferraz & Jordan (1984). In the upper one (lwp09256, dash-dot line), those emissions are very difficult to distinguish. These spectral changes are typical of CTTS and are probably due to inhomogeneities in the accretion flow and/or changes in the orientation of the star as it rotates. The temperature of the upper spectrum is almost a thousand Kelvin higher than that of the lower spectrum (Table 1).

We have de-reddened all the spectra using $`A_V=0.51`$ magnitudes (Gullbring et al. 1998). This reddening determination is based on the shape of the optical spectrum and not only on broad band photometry. It is therefore more reliable than previous values of $`A_V`$.

In addition to its science objects, IUE also observed a suite of standard stars. We have compared spectral characteristics of BP Tauri with two of these standards (Table 1). The color temperature of the BP Tauri spectra was evaluated using the continuum fitting procedures described in Sections 4.2 and 5.3. HBC 399 (also known as V 827 Tau) was the only WTTS observed by IUE with enough signal-to-noise to measure its color temperature. It is classified as K7, M0 and has an H$`\alpha `$ equivalent width of $`1.8\AA `$ (Herbig & Bell 1988). We take the reddening to be $`A_V=0.28`$ magnitudes (Kenyon & Hartmann 1995). IUE also observed main-sequence stars including 61 Cyg B, with spectral type K7V. According to the SIMBAD database its parallax is $`0.28`$ arc-sec, which means that the reddening is negligible. As Table 1 shows, the continuum level and the flux in the Mg II line in BP Tauri are one and three orders of magnitude larger than those in HBC 399 and 61 Cyg B respectively. This observation confirms the fact that the intrinsic photospheric continuum of the star is invisible in the Balmer range. As shown in Section 4 this ultraviolet excess can be understood as being due to the accretion process.

## 3 The continuum and the Mg II blend

As mentioned in the introduction, relationships between line emission and continuum level have been published by a number of authors. The study of these relationships is interesting because they illuminate what fraction of the observed line emission is due to the accretion process, as opposed to the wind or the stellar photosphere. For BP Tauri in the ultraviolet, two main articles reaching opposite conclusions have been published. These articles use subsets of the entire IUE dataset for BP Tauri.

From six spectra taken in 1986, SVH90 reached the conclusion that variations in the flux of Mg II blend can be correlated with continuum variations. They also looked at O I ($`1300\mathrm{\AA }`$), Si II ($`1810\mathrm{\AA }`$) and C II ($`1335\mathrm{\AA }`$) and reached the same conclusion. This would seem to indicate that the bulk of the emission in those lines comes from a process directly related to the formation of the accretion spot. From a set of twenty-two spectra taken in 1992, GdCF97 reached the conclusion that the Mg II blend, Si II and C IV ($`1550\mathrm{\AA }`$) light curves are dominated by short timescale fluctuations, whereas O I and He II ($`1640\mathrm{\AA }`$) are well correlated with variations in the UV continuum. Mg II, Si II and C IV are collisionaly excited lines, but O I and He II can be excited by recombination process.

The IUE data utilized in this article includes all the data used by SVH90 and GdCF97 and much more and therefore it can resolve the issue of whether or not the Mg II line flux correlates with the UV continuum. Following GdCF97, we have measured the UV continuum variations in a window of $`50\mathrm{\AA }`$ width centered at $`2900\mathrm{\AA }`$. This region of the spectra stands between two groups of Fe II lines (Brown et al, 1984) and therefore it samples closely the real continuum of the star.

Disregarding those spectra for which the continuum or the Mg II blend are affected by “extraordinary” events (defined as having a “quality” flag different from zero), we have 43 measurements. The dispersion in the continuum flux measured in this way, is $`m_{cont.}1`$ magnitudes. The dispersion in the flux of the line is $`m_{\mathrm{MgII}}0.8`$ magnitudes. Figure 2 shows that, indeed, no correlation exists between the flux in the blend and the UV continuum flux.

The difficulty of interpreting low-resolution observations of the Mg II line is compounded by the fact that high resolution observations, as described by GdCF97 show that the line has a strong wind component. If part of the line is created in the accretion flow and part in an extended wind, a measurement of the low resolution flux is sampling to different regions of the matter flow (see Section 5.3).

We have also compared the color temperature of the continuum (with the continuum as found in Section 4.2) and the amount of continuum flux (Figure 3). The errors in the temperature are $`10\%`$. As Figure 3 shows, no correlation exists between color temperature and flux in the continuum. This implies that there is no correlation between the flux in the line and the color temperature. We defer further discussion of this Figure until Section 5.3.

## 4 Modeling the ultraviolet continuum

In order to interpret the spectra, a physical model is needed; we have decided to use the work of CG98. A full description of the model can be found there. We describe only those aspects of it directly relevant to our data. In this section we will find the accretion rate, the energy flux and the spot size indicated by each of our spectra.

### 4.1 The Model

CG98 construct a model in which the observed spectral excess in CTTS is the sum of optically thick emission coming from the heated photosphere, and optically thin emission coming from the pre- and postshock regions. The shock is assumed to be parallel to the surface of the star. The heated postshock material radiates mainly in X-rays. Half of that radiation is reprocessed into optical and UV wavelengths in the photosphere; the other half is absorbed in the preshock region and reprocessed in optical and UV wavelengths. A fraction of this emission will be re-emitted towards the star and will contribute to heating the photosphere, whereas the remainder will travel away from the star.

The amount of observed radiation at each wavelength depends on two factors. The first factor is the strength of the potential depth of the star, M/R, where M is the mass and R is the radius of the star. This ratio determines the temperature of the postshock gas and therefore, everything else being equal, the hardness of the incident radiation. We have assumed $`\mathrm{M}/\mathrm{R}=0.25\mathrm{M}_{}/\mathrm{R}_{}`$ for BP Tauri (Gullbring et al. 1998). The second factor is the kinetic energy flux carried by the accretion flow. From mass conservation we have that the density of the accretion flow is $`\rho =\dot{M}/(A_{intspot}v_{ff})`$, where $`\dot{M}`$ is the accretion rate, $`v_{ff}`$ is the free-fall velocity of the material. $`A_{intspot}`$ is the intrinsic area of the spot, equal to the observed area divided by $`\mathrm{cos}\theta `$, the inclination of the spot with respect to the line of sight. To consider inclination effects explicitly, we write $`\rho =\dot{M}\mathrm{cos}\theta /(A_{spot}v_{ff})`$, where now $`A_{spot}`$ is the observed spot area. We then write the kinetic energy flux (or the particle flux) as:

$$=\frac{1}{2}\rho v_{ff}^{}{}_{}{}^{3}=\mathrm{9.8\; 10}^{10}(\frac{\dot{M}\mathrm{cos}\theta }{10^8M_{}yr^1})(\frac{M}{0.5M_{}})(\frac{R}{2R_{}})^3(\frac{f}{0.01})^1\mathrm{ergs}\mathrm{cm}^2\mathrm{s}^1$$
(1)
in which $`v_{ff}300\mathrm{km}\mathrm{sec}^1`$, the free-fall velocity from a truncation radius equal to $`\mathrm{R}_\mathrm{T}=5\mathrm{R}_{}`$ (Meyer, Calvet & Hillenbrand 1997), $`f`$ is the apparent fractional size of the accretion region or spot, $`A_{spot}=f4\pi R_{}^2`$ and $`\mathrm{cos}\theta `$ is the inclination of the spot with respect to the line of sight. The factor of $`\mathrm{cos}\theta `$ is obtained by assuming that the star has only one spot. If the region producing the emission is not one spot but a set of spots or a deformed ring (Mahdavi & Kenyon 1998) $`\mathrm{cos}\theta `$ should be interpreted as the cosine of the mean inclination of the region responsible for the emission. In the following analysis we do not find the intrinsic accretion rate but $`\dot{M}\mathrm{cos}\theta `$, the “projected” accretion rate. Similarly, the value that is found for $`f`$, the spot size, corresponds to the observed spot size and not the intrisic one. Both are related by $`f=f_{int}\mathrm{cos}\theta `$.

CG98 assumes that the truncation radius $`\mathrm{R}_\mathrm{T}`$ remains constant. Modeling by a number of authors (e.g. Königl 1991, Ghosh 1995, Ostriker & Shu 1995) suggests that $`\mathrm{R}_\mathrm{T}\alpha \dot{\mathrm{M}}^{2/7}`$. This result can be obtained in a number of ways, but the simplest is by realizing that it is the only dimensionally correct dependence between $`\mathrm{R}_\mathrm{T}`$ and $`\dot{\mathrm{M}}`$. Furthermore, in most models the truncation radius affects the size of the spot in such a way that a large $`\mathrm{R}_\mathrm{T}`$ implies a small spot (see Section 5.6). If we were to assume a varying truncation radius, the factor $`\zeta =1R_{}/R_T`$ would appear multiplying the accretion rate. As the truncation radius depends only weakly on the accretion rate, we will assume $`\zeta `$ to be constant .

CG98 divides the emission accretion zone into three regions: Shock/postshock, preshock, and heated atmosphere. In the shock and postshock region a number of simplifying assumptions are made to solve the mass, momentum, and energy conservation equations: the flow is assumed to be plane-parallel and strictly along field lines (i.e., one-dimensional) perpendicular to the stellar surface. Energy transport by conduction is neglected, as are instabilities on the flow and effects of the magnetic field on the dynamics of the gas. For the boundary conditions the shock is assumed to be strong. The temperature structure towards the stellar surface is calculated until a minimum is reached. The emission from the accretion column is calculated with CLOUDY.

The preshock region is assumed to be illuminated by the radiation field provided by the postshock calculation. The small size of this region compared to the stellar radius is used to justify an assumption of constant density through it.

The stellar atmosphere is assumed to be semi-infinite, plane-parallel, and in radiative equilibrium. The emission from it is the sum of the unperturbed stellar photosphere, plus reprocessed incident radiation. For $`\mathrm{T}_{}=4000`$ K, $`\mathrm{M}=0.5\mathrm{M}_{}`$, $`\mathrm{R}=2\mathrm{R}_{}`$ and $`\mathrm{log}=11.5`$, the effective temperature of the continuum at $`\tau _{Ross}=2/3`$ is then $`\mathrm{T}_{\mathrm{eff}}6000\mathrm{K}`$. The output of the model is a grid of the excess spectra each of which is characterized by a different value of $`\mathrm{log}`$.

Through CG98’s model, the excess Balmer continuum is shown to be due to optically thick emission from the heated atmosphere and optically thin emission from the pre- and postshock regions. The former dominates at high energy fluxes ($`\mathrm{log}>11`$) with the total spectrum becoming more like blackbody emission. The optically thin emission becomes more important as the energy flux decreases. The Paschen continuum is always dominated by optically thick photospheric emission. We refer the reader to the original paper for further details.

### 4.2 Fitting the continuum

To use CG98 models, the first step is to find the UV continuum. We have only used IUE spectra in which there are regions with S/N$`>`$2.5. As this region differs for each image, the useful wavelength coverage of each spectrum varies slightly, but in general, the useful range is reduced from the nominal IUE coverage to $`\lambda 2400\mathrm{\AA }`$ to $`\lambda 3100\mathrm{\AA }`$. The signal-to-noise constraint reduces the number of useful spectra to 45. We have developed an automated procedure that eliminates from each spectrum the regions with the resonant Mg II blend and the Fe II uv 1 blend (2771.3 to 2821.3 Å and 2576 to 2636 Å, respectively). For the remaining data, the bottom of the flux in a $`\lambda 10\mathrm{\AA }`$ wide box is found. This defines the continuum. The error in the continuum determination is found from the variations in the flux within this box. The continuum points found with this procedure are binned until we have 20 data points.

As Figure 1b shows, there are certain regions of the spectrum that sample the apparent continuum. Our procedure makes sure that the continuum spline fit goes through certain points like $`2300\mathrm{\AA }`$, $`2650\mathrm{\AA }`$, $`2900\mathrm{\AA }`$ and $`3000\mathrm{\AA }`$.

### 4.3 The Temperature of the Underlying Star

To obtain the spectrum of the excess emission we must first characterize the spectrum of the underlying K7 star. As shown in Table 1 the color temperature of the standard stars of the same spectral type as BP Tauri are lower than the canonical 4000 K of a K7 star. This is because the UV photosphere is optically thicker than the visible one and therefore it is situated close to the temperature minimum of the star. If we assume that BP Tauri is intrinsically like HBC 399, the color temperature of its photosphere should be $`3200`$K. It has been shown that the ultraviolet color temperature of early G stars is $`80\%`$ of their effective temperature (Haisch & Basri 1985). It is also shown by CG98 that the ratio of Balmer to Paschen continua for a $`\mathrm{T}_{\mathrm{eff}}=4000\mathrm{K}`$ star is better reproduced by assuming that the underlying star has a temperature of 3580 K in the Paschen continuum. They attribute this to the possibility that accretion is occurring into a cool spot. What the temperature of this spot is in the Balmer range cannot be estimated without further modeling.

Therefore, for the purposes of calculating the excess emission, we assume that the underlying star has a temperature of $`3580`$K. This estimate is consistent with the color temperature of the standards and allows us to use CG98 models. We will examine the consequences of this assumption in Section 5.4.

### 4.4 Finding the Energy Flux and the Spot Size

Once the underlying star has been subtracted, the resulting observed excess emission is compared with model spectra given by CG98. To perform this comparison, we obtain the mean value and the slope (assuming a straight line) of the logarithm of the observed spectrum. These two parameters have zero covariance and so any correlation between them is not due to the fitting procedure.

In CG98 models, the slope of the predicted excess emission for a fixed value of $`\mathrm{M}/\mathrm{R}`$ depends only on $`\mathrm{log}`$. By fitting the continuum we are able to associate a value of $`\mathrm{log}`$ with each observed spectrum. When our values for $`\mathrm{log}`$ extend beyond their published grid, we have linearly extrapolated their models. Once the value of $`\mathrm{log}`$ is found, we find the relative size of the spot $`f`$ by comparing the mean value of the theoretical excess spectrum to that of the observed excess spectrum.

In this way we obtain the values for $`\mathrm{log}`$ and $`f`$ quoted in table 2. Errors in $`f`$ are $`\sigma _f/f10`$% and are a combination of the errors in the value of the flux at each wavelength, and errors in the determination of the continuum. Errors in $`\mathrm{log}`$ are $`\sigma _{\mathrm{log}}/\mathrm{log}1`$%. With Equation 1 we obtain $`\dot{M}\mathrm{cos}\theta `$, the accretion rate uncorrected for projection effects. With this model it is impossible to obtain the “unprojected” accretion rate $`\dot{M}`$, unless $`\mathrm{cos}\theta `$ is known: $``$ is independent of $`\theta `$ but the observed spot area $`f`$ is proportional to $`\mathrm{cos}\theta `$ and so their product is proportional to $`\dot{M}\mathrm{cos}\theta `$. The errors in $`\dot{M}\mathrm{cos}\theta `$ are $`\sigma _{\dot{M}\mathrm{cos}\theta }/\dot{M}\mathrm{cos}\theta 30\%`$. All errors include an error in $`\mathrm{A}_\mathrm{V}`$ of 0.05 magnitudes.

## 5 Analysis

The average size of $`f`$ is $`\overline{f}=\mathrm{4.4\; 10}^3`$ with a dispersion of $`\sigma =\mathrm{5\; 10}^3`$. The average of $`\mathrm{log}`$ is $`\overline{\mathrm{log}}=11.7`$, and the dispersion is $`\sigma =0.2`$ (see Figure 4a and 4b). It should be noted that the value of $`f`$ is a lower limit to the real spot size, given that what one observes is the projected spot area. Also, because of the period changes mentioned in the introduction, we may not be observing the same spot (or set of spots) all the time, but their size and distribution over time. Figure 4c shows the distribution of $`\dot{M}\mathrm{cos}\theta `$. The average is $`\overline{\dot{M}\mathrm{cos}\theta }=\mathrm{1.6\; 10}^8M_{}\mathrm{yr}^1`$, and the dispersion is $`\sigma \mathrm{0.8\; 10}^8M_{}/yr`$. The dispersion of $`f`$ is larger than the observational errors, but the dispersions in $`\mathrm{log}`$ and $`\dot{M}\mathrm{cos}\theta `$ are consistent with the errors.

As mentioned in the introduction, a wide range of spot sizes have been identified by other authors. We find spot sizes ranging from 0.06 to 3% the size of the visible star disk. The average accretion rate we observe agrees with the most recent determinations for BP Tauri, based on single optical spectra: $`\dot{M}2\mathrm{3\; 10}^8M_{}\mathrm{yr}^1`$ (Gullbring et al. 1998, CG98).

The logarithm of the energy flux, $`\mathrm{log}`$, does not depend on the inclination angle $`\theta `$. Therefore, the variations observed in Figure 4b are due to variations in $`\dot{M}`$ and/or the intrinsic size of the accretion spot. The dispersion in $`\mathrm{log}`$ implies a fractional variation of $`50\%`$ in $`\dot{M}/f_{int}`$, where $`f_{int}`$ is the intrinsic (i.e. deprojected, $`f_{int}=f/\mathrm{cos}\theta `$) spot size. This dispersion is consistent with the errors in $``$. Variations in $`\dot{M}`$ cannot be decoupled from those in $`f_{int}`$. What can be said is that changes in the accretion rate $`\dot{M}`$ or $`f_{int}`$ are less than 50% during the time spanned by the observations. This upper limit puts constraints on the inhomogeneities of the disk (that may contribute to the variations in the accretion rate) and in the geometry of the stellar field.

### 5.1 Variable Extinction

Recently, Bouvier et al. (1999) reported on photometric observations of the CTTS AA Tauri that can be interpreted as periodic occultations of the stellar photosphere by a warp in the accretion disk. One of the effects of this warp is to act as a source of variable reddening. This mechanism might also play a role in the IUE observations of BP Tauri. To find out if that is the case, we performed the following analysis. The procedure of finding which model corresponds with an observed spectrum has three free parameters: reddening, $`f`$ and $`\mathrm{log}`$. In this paper we assume a fixed reddening and find the values of $`f`$ and $`\mathrm{log}`$ that match the observed spectrum to the models. We can also assume a fixed spot size and find the values of $`\mathrm{log}`$ and reddening that match the observed spectra. We did so, assuming $`f=\mathrm{4\; 10}^3`$. For the cases in which this procedure succeeded in finding a $`\mathrm{log}`$ and $`A_V`$ that fit the data, we obtained a large range of values for $`A_V`$, from 0.4 to 10 magnitudes, but most of them bigger than 1. The largest value that we found in the literature for a reddening determination of BP Tauri was 0.85 (Valenti et al. 1993). Therefore, most of the reddenings determined in this way do not fit the optical data. Photometric observations in various bands simultaneously have also failed to detect color changes in BP Tauri analogous to those reported by Bouvier et al. (1999) for AA Tauri. We therefore conclude that the changes in $`\mathrm{log}`$ and $`f`$ are not due to changes in reddening.

### 5.2 The accretion luminosity versus the continuum flux

Figure 5 shows a tight correlation between $`\dot{M}\mathrm{cos}\theta `$ and the continuum flux, which implies that the observed (projected) accretion rate determines the brightness of the UV continuum. We can convert the observed UV flux to luminosity and the model accretion rate to accretion luminosity. We obtain

$$\mathrm{log}(L_{acc}/L_{})=(1.1\pm 0.1)\mathrm{log}(L_{2900\AA }/L_{})(0.5\pm 0.05)$$
(2)
A similar relationship has been found by GC98 and Gullbring et al. 1997, who observe a very clear correlation between the accretion luminosity and the brightness in the U-band, based on individual observations of various CTTSs in the optical. The empirical relation they find is

$$\mathrm{log}(L_{acc}/L_{})=1.09\mathrm{log}(L_U/L_{})+0.98$$
(3)
in which $`L_{acc}`$ is the accretion luminosity (which is proportional to the accretion rate) and $`L_U`$ is the unreddened luminosity of the star in the U-band once the intrinsic (photospheric) luminosity has been subtracted.

An independent way to verify the accretion rates obtained here is to use Equation 3 as to predict accretion luminosity and then compare it with the accretion luminosities from Table 4. Given the rapid photometric variations of CTTSs, one needs to use concurrent $`L_{acc}`$ and $`L_U`$ to test this relationship. We searched the photometric database maintained by William Herbst222http://www.astro.wesleyan.edu/ bill/research/ttauri.html#ov2 described in Herbst et al. (1994) and found two measurements of the U-magnitude close enough to the IUE measurements: RJD=5615.51 (IUE at 5615.85) and RJD=5620.99 (IUE at 5620.97), which imply $`\mathrm{log}(L_U/L_{})=1.88`$ and $`\mathrm{log}(L_U/L_{})=2.05`$ respectively. These determinations assume $`L_{BPTauri}=0.925L_{}`$ (Gullbring et al. 1998). Using Equation 3, the predicted values are $`\mathrm{log}(L_{acc}/L_{})=1.072`$ and $`\mathrm{log}(L_{acc}/L_{})=1.25`$ respectively. We observe $`\mathrm{log}(L_{acc}/L_{})=1.101`$ and $`\mathrm{log}(L_{acc}/L_{})=1.303`$, by interpolating Figure 5. This is consistent with Equation 3 within the errors of our data. Admittedly, two points out of 45 is a very limited sample, but the fact that they obey Equation 3 lends support to our determinations of the accretion rate.

### 5.3 The Line Flux and the Continuum Temperature

Figure 3 shows that no correlation exists between the continuum temperature and the flux in the continuum. Within the context of CG98 model it is possible to understand this result. The Balmer continuum is the sum of optically thick emission from the heated photosphere and optically thin emission from the pre-shock region. Only when the value of $``$ is very high does the observed color temperature approach the temperature of the heated photosphere. The flux in the continuum at $`2900\mathrm{\AA }`$ depends, again in the context of CG98’s model, not only on the spot size, but also on the value of $``$. Two spectra with the same accretion spot size but different $``$ will have different values of continuum flux. Therefore, the color temperature is a measure of two different things, the spot size and the energy flux.

We do not observe any correlation of the Mg II line flux with $`\mathrm{log}`$, the accretion rate $`\dot{M}\mathrm{cos}\theta `$, or the spot size $`f`$.

It has been argued by Hartmann (1990) that strong emission lines (like H$`\alpha `$ and Mg II) in the spectra of CTTS are formed in an extended envelope. His argument is based on the similar widths of the strong lines. Table 1 shows that the Mg II line is very bright in BP Tauri, compared with that of non-CTTS of the same spectral type. There have also been unsuccessful attempts to model the Mg II line blend as purely chromospheric (e.g. Calvet, Basri & Kuhi 1984). All these arguments point to the conclusion that the Mg II line blend is a consequence of the accretion process. The lack of correlation with the accretion rate maybe the result of time delay: if the Mg II line is created in the funnel flow one would expect a delay of $`5`$ days between a spectral manifestation in Mg II and the landing on material on the star. An analogous delay is expected for the wind component of the line.

On the other hand GdCF97 suggest that magnetospheric flaring may give origin to increased Mg II emission. Such flaring is present in the model by Shu et al. (1994). If this is the case, it is unclear whether or not any time delay would correlate the line flux with the accretion rate. The flaring of the magnetosphere maybe be causally unrelated to the process of accretion. If this is the case we may be seeing a Mg II line corresponding to more than one source: a wind component, a funnel component and a magnetospheric component.

As the model used here to find the accretion rate is sensitive only to regions below the accretion shock, the lack of correlation of the Mg II line with the accretion rate is telling us that the dominant component of the line is not created in the region below the accretion shock.

### 5.4 The Accretion rate and the spot size

Figure 6a shows that there is a correlation between the spot size, $`f`$, and the accretion rate, $`\dot{M}\mathrm{cos}\theta `$. The slope of the points in Figure 6a is such that $`\dot{M}\mathrm{cos}\theta \alpha f^{0.5\pm 0.05}`$.

The correlation between $`\dot{M}\mathrm{cos}\theta `$ and $`f`$ implies (via Equation 1) a correlation between $`\mathrm{log}`$ and $`f`$, as shown in Figure 6b. Here $`\alpha f^{0.5}`$. Assuming a constant accretion rate of $`\dot{M}\mathrm{cos}\theta =\mathrm{1.6\; 10}^8M_{}\mathrm{yr}^1`$, we obtain the solid line shown in Figure 6b, which shows that smaller spots have under-average $`\dot{M}\mathrm{cos}\theta `$ and larger spots have over-average $`\dot{M}\mathrm{cos}\theta `$. This merely reflects the fact that the observed (projected) accretion rate is not constant.

As Equation 2 shows, the observed (in the detector) ultraviolet flux, in $`\mathrm{ergs}/\mathrm{sec}\mathrm{cm}^2`$ is such that $`F_{UV}\alpha \dot{M}\mathrm{cos}\theta `$. This implies, from Figure 6a, that $`F_{UV}\alpha f^{0.5}`$. The flux on the surface of the star is proportional to the detected flux divided by the emitting area. Therefore, at the stellar surface, $`F_{UV}^{sur}\alpha f^{0.5}`$ and the particle flux and the ultraviolet flux are proportional to each other.

The proportionality between particle and surface UV flux is implicit in CG98 models because a large particle flux produces a large UV flux. What needs an explanation is the physical origin of only one of the correlations ($`F_{UV}^{sur}\alpha f^{0.5}`$, $`\alpha f^{0.5}`$ or $`\dot{M}\mathrm{cos}\theta \alpha f^{0.5}`$). Once one is explained, the others follow.

CG98 observe values of energy flux and spot sizes suggesting that, for increasing mass accretion rate, material will fall onto the star over a larger surface area, while the energy flux $``$ remains more or less constant. They reach this conclusion by single optical observations of 15 CTTSs (Figure 6c). We also find that a bigger mass accretion rate implies a bigger spot size, but the energy flux that we calculate does not remain constant. While there is a tendency in CG98’s data towards larger accretion spots with larger area (as implied by the constant $``$), this tendency is not strong. Most of their stars lie on the locus defined by the BP Tauri observations and therefore, their data could be construed as confirming the analysis done here. However, their stars are not a homogeneous population, as they encompass a range of radii and masses.

The correlations from Figure 6 are somewhat sensitive to the estimate of the temperature of the underlying star. A higher temperature star will produce smaller inferred accretion spots with higher $`\mathrm{log}`$. For the 3580 K, the temperature that we have chosen here, the ratio between the observed UV continuum flux and the photospheric flux is $`30`$. For 4000 K this ratio $`10`$. For the 4000K case, the accretion spots are 0.3 times the size of the spots obtained with the 3580K case. The values of $`\mathrm{log}`$ and $`\dot{M}\mathrm{cos}\theta `$ are 1.04 times and 0.6 times the values of the 3580K case. The shape of the correlations from Figure 6, including slopes and overall values, does not change, within the errors in the calculated quantities. In other words, even though a given spot will change size when the assumed temperature of the underlying star changes, the overall relationship between spot size, energy flux, and accretion rate does not change.

### 5.5 Simulating the Correlation between Accretion Rate and Spot Size

A possible explanation for the correlation from Figure 6a lies in the geometry of the accretion region. As the spot size, $`f`$, and the accretion rate, $`\dot{M}\mathrm{cos}\theta `$, are both proportional to the geometric factor, $`\mathrm{cos}\theta `$, one would expect a correlation between them. When the projected size of the accretion spot is large (large $`f`$ due to a small inclination angle) we should see large projected accretion rates. The larger range of the spot sizes ($`1.6`$ dex) compared to the accretion rates ($`1.2`$ dex) could be explained by different intrinsic distributions.

To decide whether or not the correlation of Figure 6a is due to geometrical effects we have perfomed a simulation. We assume that the situation is one in which we are looking at an opaque disk from an inclination angle $`i`$. The spot has the colatitude angle $`\theta _{col}`$. As all the spectra in our data set show ultraviolet excess, one always see the spot. This gives the constraint that $`i\pi /2\theta `$. Let $`\widehat{𝐢}`$ be the unit vector pointing from the center of the star to the observer. Then $`\widehat{𝐢}=(\mathrm{sin}(i),0,\mathrm{cos}(i))`$. Let $`\widehat{𝐧}`$ be the unit vector perpendicular to the spot area. Then $`\widehat{𝐧}=(\mathrm{sin}(\theta _{col})\mathrm{cos}(\omega t),\mathrm{sin}(\theta _{col})\mathrm{sin}(\omega t),\mathrm{cos}(\theta _{col})`$, where $`\omega =2\pi /8`$ rad per day, the angular frequency of the star, assumed to have a period of 8 days. Therefore $`t`$ maps the phase in which we are observing the star. The cosine of the angle between the spot and the observer (which we have called here $`\theta `$) is given by $`\mathrm{cos}\theta =\widehat{𝐢}\widehat{𝐧}=\mathrm{sin}(i)\mathrm{sin}(\theta _{col})\mathrm{cos}(\omega t)+\mathrm{cos}(i)\mathrm{cos}(\theta _{col})`$. We then generate uniform random points between 0 and 8 for $`t`$. This procedure should produce the observed distribution of $`\mathrm{cos}\theta `$.

We assume that $`\dot{M}`$ and $`f_{int}=f/\mathrm{cos}\theta `$ have normal distributions. The simulation therefore has 6 parameters: the means and sigmas of the normal distributions, the colatitude angle of the spot and the inclination of the system. For a given set of parameters, we randomly choose 2000 points for $`\dot{M}`$, $`f_{int}`$ and $`\mathrm{cos}\theta `$ and multiply $`\dot{M}`$ by $`\mathrm{cos}\theta `$ and $`f_{int}`$ by the same $`\mathrm{cos}\theta `$. The values of $`f_{int}\mathrm{cos}\theta `$ and $`\dot{M}\mathrm{cos}\theta `$ are constained to be larger than the minima of their distributions. In other words, we are assuming that the accretion rate and the intrinsic size are uncorrelated and that any correlation comes from the $`\mathrm{cos}\theta `$ that multiplies both quantities.

We do not find any combination of parameters that will satisfy all the constraints of the problem: Figures 4a, 4b and 6a. All the simulations that produce a correlation between $`\dot{M}\mathrm{cos}\theta `$ and $`f`$ give a slope close to 1 for the log-log plot of the two quantities. The width of the gaussians determine the scatter of this correlation. This means that the basic hypothesis of the model (that $`\dot{M}`$ and $`f_{int}=f/\mathrm{cos}\theta `$ are uncorrelated gaussians) is not verified by the simulation.

Dendy, Helander & Tagger (1988) have considered the possibility that the intrinsic distribution of accretion rates should be a power law. They arrive to this conclusion based on models of self-organized criticality. The range of accretion rates found here is not large enough to test this hypothesis.

While it is clear that a full understanding of Figure 6a should include the effects of changing inclinations of the spot as the star rotates, the fundamental explanation is not geometry. We have to look for other mechanisms.

It is possible that some of the inferred changes in energy flux are due to “limb brightening”. Assuming that spot size does not change for short cycles, a small spot size would imply a high spot inclination, and therefore we would observe higher regions of the spot (farther away from the center of the star), which the models predict have higher temperatures (Lamzin 1998). Therefore, these higher regions will look like regions of high energy flux. One would therefore expect a correlation between $`\mathrm{log}`$ and $`f`$, as observed here. If “limb brightening” were indeed an important effect, the slope of the spectra would not measure the value of $``$ but the temperature of that layer, which would change for different orientations. The CG98 models assume that the spot is observed from directly from above, and therefore it is difficult to account for this effect. Resolving this issue will have to wait for more detailed models.

### 5.6 Theoretical Expectations of the Relationship between Spot Sizes and Accretion Rates

The total size of the accretion region is determined by geometry of the stellar magnetic field. An accretion spot of variable size would therefore imply that the field itself is changing. So far, the most sophisticated idea of the accretion process is the so-called ’X-wind’ model (Shu et al. 1994). This model has been developed for a magnetic dipole aligned with the disk, even though most authors agree that such geometry is an oversimplification of the more complicated real geometry of the field. In the X-wind model, the fractional coverage of the accretion spot with respect to the area of the star (the quantity we have called here $`f`$) is given by $`f=\sqrt{1R_{}/R_X}\sqrt{11.5R_{}/R_X}`$ (Ostriker & Shu 1995), where $`R_X`$ is the radius of the ’X-point’, that coincides to a few percent with the corotation radius of the star and with the truncation radius $`R_T`$. Furthermore, $`R_X\alpha \dot{M}^{2/7}`$. For $`R_{}/R_X0.2`$ as assumed in this paper, $`f6\%`$. Let’s assume the accretion rate increases by a factor of 2, from the value that produces by $`R_{}/R_X0.2`$. $`R_X`$ decrease bys a factor of 0.8 and then $`R_{}/R_X0.24`$. Then $`f7\%`$. According to Figure 6a, $`f`$ should increase by a factor of 4 when the accretion rate increases by a factor of 2. The increase in spot size with accretion rate is in the right direction but is not fast enough to explain the correlation of Figure 6a.

The model by Mahdavi & Kenyon (1998) uses a tilted dipole and consider the accretion as occurring onto a thin ring on the surface of the star. For $`R_{}/R_T0.20.3`$ and a dipole inclined by 60 degrees with respect to the plane of the accretion disk, this model predicts $`f24\%`$, and $`f\alpha \dot{M}^{2/7}`$ which is again a very slow increase in the right direction.

It is also clear that the problem remains if we include explicitly the varying truncation radius in the correlation from Figure 6a: $`f\alpha (\dot{M}\mathrm{cos}\theta )^2(1R_{}/R_T)^2`$. The increase in the accretion rate dominates over the slow decrease in the factor $`\zeta =(1R_{}/R_T)`$.

The root of the problem is probably in the assumption of a magnetic dipole as the magnetic field of the star. Even though this is a good assumption at large distances from the star it is clear that it is an oversimplification if one wishes to explain the geometry of the accretion spots, which depends on the field close to the star.

In Ostriker & Shu (1995) and Mahdavi & Kenyon (1998), the predicted accretion spots are on the upper tail of our distribution. This has also been observed by CG98, which conclude that the accretion process has to be inhomogeneous. If not all the available area is covered at the same time, the observed accretion spot will be a fraction of the available area.

Independently of these considerations one could argue that, according to Figure 6a, matter loading on the field lines is not very efficient. We see that the accretion rate increases only with the linear dimensions (e.g. the radius or the circumference) of the accretion spot: $`\dot{M}f^{0.5}`$. So a bigger area produces a more diluted accretion, and therefore, a weaker surface flux $`F_{UV}^{sur}`$. Another way of saying this is that the spot size is very sensitive to changes in accretion rate.

The models that we use here are not complete: they assume a shock parallel to the surface of the star, where the real shock may be oblique. They also assume free-fall velocity, which, as Ostriker & Shu (1995) show, is not the most appropiate hypothesis. The inclusion of these two effects would increase the observed spot size, but it is unlikely that they will change the slope of Figure 6a. To rule out this possibility, more complete simulations are needed.

In summary, the IUE data on BP Tauri puts a strong constraint in any model that seeks to explain the size of accretion spots. Accurate modeling should consider the geometry of the observations (i.e. inclination angle, spot latitude), a more realistic geometry for the magnetic field, oblique shocks, the real velocity of the material falling onto the star, and possibly the effects of limb brightening. There may also be an unmodelled efficiency factor that determines the way in which flux tubes are actually filled with material, in such a way that the accretion rate scales as the radius of the tube.

## 6 Conclusions

We have analyzed all the observations of BP Tauri taken by the International Ultraviolet Explorer in the low resolution ($`\lambda 6\mathrm{\AA }`$), long wavelength range (from $`\lambda =1850\mathrm{\AA }`$ to $`\lambda =3350\mathrm{\AA }`$). We found 61 useful spectra.

Compared with spectra of standard stars, all the BP Tauri spectra show strong ultraviolet excess and color temperatures around $`8000`$ K. Using 43 spectra, we observe variability in the ultraviolet continuum of $`m_{cont.}1`$ magnitude and variability in the Mg II line flux of $`m_{\mathrm{MgII}}0.8`$ magnitudes. We observe no correlation between continuum flux and Mg II line flux. This result resolves a controversy in the literature regarding the relationship of the continuum flux and Mg II line flux in BP Tauri. We conclude that the bulk of the Mg II in BP Tauri is not emitted below the accretion shock. Therefore, it may be emitted by the accretion column, the stellar wind and/or magnetospheric flaring. We do not observe any correlation between color temperature and the average value of the continuum.

Within the context of the magnetospheric accretion paradigm, we have used models by CG98 to obtain the energy flux into the star due to the accretion process, and the spot size for each of the spectra of BP Tauri. We assume a constant inner truncation radius for the disk of $`5R_{}`$. With 45 spectra, we find an average energy flux of $`\mathrm{5\; 10}^{11}\mathrm{ergs}\mathrm{cm}^2\mathrm{sec}^1`$ and an average projected spot size of $`\mathrm{4.4\; 10}^3`$ times the area of one hemisphere of the star. The latter result coincides with previous observational determinations of spot sizes. From the energy flux and the spot size we obtain the accretion rate times $`\mathrm{cos}\theta `$, where $`\theta `$ is the average inclination angle of the region that produces the emission. The average of this projected accretion rate is $`\mathrm{1.6\; 10}^8M_{}\mathrm{yr}^1\pm 30\%`$, consistent with values found by other researchers. This values are weakly sensitive to the value of the underlying star used obtain the continuum excess. No correlation is observed between the flux in Mg II line and the energy flux, the accretion rate or spot size, affirming the conclusion from the previous paragraph.

The nature of the data and the models is such that we cannot say anything about intrinsic variations in the accretion rate or the real (i.e., deprojected) spot size, but we conclude that the dispersion in their ratio is $`\dot{M}/f50\%`$, consistent with the errors.

We observe a correlation between the ultraviolet flux and the accretion rate, similar to correlations found by other researchers between the U-band flux and the accretion rate: $`\mathrm{log}(L_{acc}/L_{})=(1.1\pm 0.1)\mathrm{log}(L_{2900\mathrm{\AA }}/L_{})(0.5\pm 0.05)`$. This correlation is implicit in the models by CG98. Also, we observe a correlation between the accretion rate $`\dot{M}\mathrm{cos}\theta `$ and the projected spot size $`f`$: $`\dot{M}\mathrm{cos}\theta \alpha f^{0.5\pm 0.05}`$. Together, the two correlations imply that the UV flux in the surface of the star is proportional to the particle flux.

These results are significant, in spite of the large errors involved in the determinations of each quantity. These correlations stand in contrast to results obtained by Calvet & Gullbring (1998). Using individual optical observations of 15 CTTSs, they concluded that the accretion rate increases linearly with spot size. A resolution of this discrepancy will have to wait until we have applied the magnetospheric accretion models to other stars.

Using a simulation, we conclude that the geometry of the system determined by the observer, the rotation axis and the inclination of the star can not explain the correlation between accretion rates and spot sizes. One must look at other more fundamental explanations.

Models of magnetospheric accretion by Ostriker & Shu (1995) and Mahdavi & Kenyon (1998) do predict that the spot size should increase with the accretion rate , but in both models, the predicted increase is much slower than what we found here ($`f\alpha \dot{M}^2`$).

The correlation between accretion rate and spot size could be due to an unmodelled efficiency factor that determines how the spot is filled. Conversely, it indicates that the spot size is very sensitive to accretion rate. Future modeling should take into account the possibility of a non-dipole field for the star, oblique shocks, limb brightening and the effects of the geometry determined by the observer. The correlation between accretion rate and spot size puts a strong constraint in these models.

We would like to thank Nuria Calvet and Erik Gullbring for providing us with their models. Dr. Calvet in particular spent quite a bit of time helping us to understand the nature and strengths of their models. We are grateful to F. Shu and J. Bouvier for suggestions that led to the improvement of the text. Thanks are also due to Debi Howell-Ardila who edited the manuscript for language. We acknowledge the support of the NASA ADP program through grant # NAG5-3471.

Figure 1: (a) Average of all the available spectra (solid line) and the normalized variance (dotted line). The feature near $`3050\mathrm{\AA }`$ is a reseaux mark. (b) Two sample spectra (lwp18976, solid line, lwp09256, dash-dot line) uncorrected by reddening. Notice the strong variation in slope and flux level. The position of some spectral features is indicated (Brown et al. 1984).

Figure 2: Mg II line flux vs. UV continuum flux. The continuum flux has been calculated from a box $`50\mathrm{\AA }`$ wide centered at $`2900\mathrm{\AA }`$. No correlation between the two quantities is observed, which means that the MgII line is not sensitive to the characteristics of the accretion spot.

Figure 3: Average of the logarithm of the flux versus the color temperature. Both have been calculated from the continuum fitting procedure (Section 4.2). No correlation is observed between the two quantities.

Figure 4: (a) Histogram for the calculated values of $`f`$. $`\overline{f}=\mathrm{4.4\; 10}^3`$ with a dispersion of $`\sigma =\mathrm{5\; 10}^3`$. The errors are $`10\%`$. The bar indicates the error at $`f=0.0095`$. (b) Histogram for the calculated values of $`\mathrm{log}`$, the logarithm of the energy flux into the star; $`\overline{\mathrm{log}}=11.7`$ and the dispersion is $`\sigma =0.2`$. The error is $`1`$%, as indicated by the bar centered in 11.4. With $`\mathrm{log}`$ and $`f`$ we calculate $`\dot{M}\mathrm{cos}\theta `$, the accretion rate uncorrected for projection effects. As shown in (c), $`\overline{\dot{M}\mathrm{cos}\theta }=\mathrm{1.6\; 10}^8M_{}\mathrm{yr}^1`$ and the dispersion is $`\sigma \mathrm{0.8\; 10}^8M_{}/yr`$. The errors are $`30`$%, as indicated by the bar for $`\dot{M}\mathrm{cos}\theta =\mathrm{1.3\; 10}^8M_{}\mathrm{yr}^1`$.

Figure 5: The projected accretion rate vs. the UV continuum flux. The UV continuum flux is taken from Table 1. This relation is similar to previously found relationships between the accretion luminosity and the U-band flux.

Figure 6: (a) The logarithm of the accretion rate vs. $`\mathrm{log}(f)`$. The boxes indicate the extent of the errors. The correlation is such that $`\dot{M}\mathrm{cos}(\theta )\alpha f^{0.5\pm 0.05}`$. (b) The correlation between $`\mathrm{log}`$ and $`\mathrm{log}(f)`$ implied by (a) via Equation 1: $`\alpha f^{0.5}`$. The solid line is what one would obtain from Equation 1, assuming a constant $`\dot{M}\mathrm{cos}\theta =\mathrm{1.6\; 10}^8M_{}\mathrm{yr}^1`$. (c) Same as (a) without the error boxes. The names indicate the positions of stars analyzed by CG98 using single optical observations (i.e. AA means AA Tauri, DK means DK Tauri and so on for all stars except GM Aur and UY Aur). In their data, some tendency can be seen to larger accretion spots with larger accretion rates.


format: text

text:
# Near-field EM wave scattering from random self-affine fractal metal surfaces: spectral dependence of local field enhancements and their statistics in connection with SERS

## I Introduction

Recently, single molecule probing by means of Surface Enhanced Raman Scattering (SERS) has been reported both on Ag single nanoparticles  and on Ag colloidal aggregates; the latter exploits the extremely large near infrared (NIR) Raman scattering cross sections of dye molecules. Bearing in mind how inefficient normal spontaneous Raman scattering is, enhancement factors of $`10^{14}`$ or larger are required to achieve single molecule detection.

Typically, SERS is known to yield Raman signals enhanced by a factor $`𝒢10^6`$ with respect to those of conventional Raman scattering. Two mechanisms are responsible for such enhancement factors: the surface-roughness-induced intensification of the electromagnetic (EM) field both at the pump frequency and at the Raman-shifted frequency (EM mechanism), and the charge-transfer mechanism. The former mechanism is widely accepted to be the most relevant one from the quantitative standpoint, providing gains of $`𝒢_{EM}>10^4`$ in most experimental configurations. Extensive theoretical work has been devoted to the explanation of the EM mechanism (cf. e.g. the reviews Refs. ), and the consensus is that what underlies such EM field enhancement (FE) factors is the roughness-induced excitation of surface-plasmon polaritons  (SPP), either propagating along a continuous surface (extended SPP), confined within metal particles (particle plasmon resonances), or even confined due to Anderson localization (localized SPP).

In light of the SERS enhancement factors estimated for single molecule detection it is evident that, in addition to the well known average SERS enhancement factors, extremely large EM fields must appear in the vicinity of SERS substrates, even if the charge-transfer mechanism is also known to be especially intense, as in Ref. . This is supported by the observation, through photon scanning tunneling microscopy (PSTM), of very intense and narrow EM modes (called hot spots) on rough metal surfaces  or rough metal-dielectric films. Interestingly, these rough metal surfaces used as SERS substrates posses, in some cases, scaling properties within a sufficiently wide range of scales (physical fractality). The substrates can present self-similarity, as the widely employed colloidal aggregates,  or self-affinity, as in the case of deposited colloids or evaporated or etched rough surfaces and thin films.

Therefore, inasmuch as the quantitative evaluation of the surface EM field is central to the SERS effect, knowledge of the EM scattering process for surface models as realistic as possible is obviously needed. In recent years, the theoretical efforts have been directed towards either describing through approximate methods realistic surface models, or using the full EM theory to study simplistic surface models,  though introducing increasingly complex properties.

In this paper we study the near EM field scattered in the vicinity of rough, one-dimensional self-affine fractal surfaces of Ag, Au, and Cu, with the aim of determining the appearance of strong local optical excitations (hot spots) and characterizing them with regard to their spatial and spectral width, their polarization, and their excitation spectra; in addition to that, the global optical response of such fractal surfaces will be studied through the statistical properties of the surface EM fields. Both local and global responses are discussed in light of the influence on SERS. For this purpose, we make use of numerical simulation calculations based on the Green’s theorem integral equation formulation, rigorous from the classical EM standpoint. Unlike recent work also for self-affine fractals (though sensitively less rougher) based strictly on magnetic field calculations, the magnitude that naturally arises in this formulation when applied to 1D surfaces and $`p`$ polarization (the one relevant for its light-SPP coupling selectivity), we fully characterize here the surface and near electric field components (crucial in SERS and other nonlinear optical effects) through simple expressions in terms of the magnetic field and its normal derivative. The necessary details of the theoretical formulation are given in Sec. II. The local optical excitations are studied in Sec. III and their statistical properties in Sec. IV, leaving for Sec. V the conclusions of this work.

## II SCATTERING FORMULATION

### A Surface integral equations

The scattering geometry is depicted in Fig. 1. A rough metal surface $`z=\zeta (x)`$ is the substrate onto which molecules are adsorbed in SERS typical experimental configurations. The semi-infinite metal occupying the lower half-space \[$`z\zeta (x)`$\] is characterized by an isotropic, homogeneous, frequency-dependent dielectric function $`ϵ^<(\omega )`$. From the medium of incidence, characterized by a frequency-dependent dielectric function $`ϵ^>(\omega )`$, a monochromatic, linearly polarized incident beam of frequency $`\omega `$ impinges on the interface at an angle $`\theta _0`$, measured counterclockwise with respect to the positive $`z`$ axis. The polarization is defined as shown in Fig. 1: the magnetic (respectively, electric) field is perpendicular to the $`xz`$ plane for $`p`$ (respectively, $`s`$) polarization, also known as the transverse magnetic (respectively, transverse electric) one.

We restrict the analysis to 1D surfaces (invariant along the $`y`$ direction), without any loss of generality as far as the physics underlying the SERS EM mechanism is concerned.  It has been shown  that the latter condition simplifies considerably the formulation based on the integral equations resulting from the application of Green’s second integral theorem (with the help of the Sommerfeld radiation condition). In such circumstances, the starting 3D vectorial problem can be cast into a 2D scalar one, where the unknown is the $`y`$-component of either the magnetic field \[$`H_y^{(p)}(𝐫,\omega )`$, with $`𝐫(x,z)`$\] for $`p`$ polarization, or the electric field \[$`E_y^{(s)}(𝐫,\omega )`$\] for $`s`$ polarization; no depolarization takes place when pure $`s`$ or $`p`$ polarized fields are incident. This simplification is very convenient from the analytical and numerical point of view, and yields straightforwardly the far field scattered intensity. We point out that, if the electric (respectively, magnetic) field is needed for $`p`$ (respectively, $`s`$) polarization, it can be calculated on the basis of Maxwell equations, as we shall see below.

Let us focus on the electric field calculation for the case of $`p`$ polarization. Evidently this is the most relevant one for our SERS problem (only $`p`$ polarized light can excite in the present configuration the SPP responsible for the EM field enhancements), and also to other interesting problems such as second harmonic generation on metal surfaces  or PSTM studies.  As mentioned above, the integral equation formulation is simplified if written in terms of the magnetic field amplitude. Our monochromatic incident field of frequency $`\omega `$ is a Gaussian beam of half-width $`W`$ in the form:

$`H_y^{(p,i)}(x,z|\omega )=\mathrm{exp}\{ık_ϵ(x\mathrm{sin}\theta _0z\mathrm{cos}\theta _0)`$ (3)
$`\times [1+w(x,z)]\left\}\mathrm{exp}\right[{\displaystyle \frac{(x\mathrm{cos}\theta _0+z\mathrm{sin}\theta _0)^2}{W^2}}],`$
$`w(x,z)={\displaystyle \frac{1}{k_ϵ^2W^2}}\left[{\displaystyle \frac{2}{W^2}}(x\mathrm{cos}\theta _0+z\mathrm{sin}\theta _0)^21\right]`$ , (4)
where $`k_ϵ=n_c^>\omega /c`$ and $`n_c^>=\sqrt{ϵ^>}`$. From now on, since a time harmonic dependence $`e^{ı\omega t}`$ is assumed, the functional dependence on frequency will be omitted unless necessary for the sake of clarity. The surface integral equations that fully describe the EM linear scattering problem for $`p`$ polarization, in the geometry of Fig. 1, are

$`H_y^{(p,i)}(𝐫)+{\displaystyle \frac{1}{4\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}[H_y^{(p,>)}(𝐫^{}){\displaystyle \frac{G^>(𝐫,𝐫^{})}{n^{}}}`$ (6)
$`G^>(𝐫,𝐫^{}){\displaystyle \frac{H_y^{(p,>)}(𝐫^{})}{n^{}}}]`$ $`=H_y^{(p,>)}(𝐫),z>\zeta (x)`$ (8)
$`=0,z<\zeta (x),`$
$`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}[H_y^{(p,<)}(𝐫^{}){\displaystyle \frac{G^<(𝐫,𝐫^{})}{n^{}}}`$ (9)
$`G^<(𝐫,𝐫^{}){\displaystyle \frac{H_y^{(p,<)}(𝐫^{})}{n^{}}}]`$ $`=0,z>\zeta (x)`$ (11)
$`=H_y^{(p,<)}(𝐫),z<\zeta (x),`$
where $`H_y^{(p,>)}(𝐫)`$ and $`H_y^{(p,<)}(𝐫)`$ are the magnetic field in the upper ($`z>\zeta `$) and lower ($`z<\zeta `$) semi-infinite half-spaces, and the normal derivative is defined as $`/n(\widehat{𝐧})`$, with $`\widehat{𝐧}\gamma ^1(\zeta ^{}(x),0,1)`$ and $`\gamma =(1+(\zeta ^{}(x))^2)^{1/2}`$. The 2D Green’s function $`G`$ is given by the zeroth-order Hankel function of the first kind $`H_0^{(1)}`$.

The four integral equations (II A) fully describe the scattering problem for $`p`$ polarization in terms of the $`y`$ component of the magnetic field. Analogous integral equations can be obtained for $`s`$ polarization dealing with the $`y`$ component of the electric field. In order to solve for the surface field and its normal derivative, defined as the functions $`H(x)`$ and $`\gamma ^1L(x)`$, two of the integral equations (note that they are not independent), typically Eqs. (8) and (11), are used as extended boundary conditions, leading to two coupled integral equations once one invokes the continuity conditions across the interface:

$`H(x)=`$ $`H_y^{(>)}(𝐫)_{z=\zeta ^{(+)}(x)}=H_y^{<)}(𝐫)_{z=\zeta ^{()}(x)},`$ (13)
$`\gamma ^1L(x)=`$ $`\left[{\displaystyle \frac{H_y^{(>)}(𝐫)}{n}}\right]_{z=\zeta ^{(+)}(x)}`$ (15)
$`={\displaystyle \frac{ϵ^>}{ϵ^<}}\left[{\displaystyle \frac{H_y^{(<)}(𝐫)}{n}}\right]_{z=\zeta ^{()}(x)},`$
with $`\zeta ^{(\pm )}(x)=lim_{\epsilon 0}(\zeta (x)\pm \epsilon )`$. The resulting system of integral equations can be numerically solved upon converting it into a system of linear equations through a quadrature scheme, the unknowns being $`H(x)`$ and $`L(x)`$. Then Eqs. (8) and (11) permit to calculate the scattered magnetic field in the upper incident medium and inside the metal, respectively.

But what if the magnitude of interest is the electric field? This is indeed the situation in SERS where the surface electric field locally excites the molecule vibrations that produce the Raman-shifted radiation that is detected. In Refs. , the EM field enhancement factor has been defined as the normalized magnetic field intensity:

$`\sigma _H(\omega )`$ $`=`$ $`{\displaystyle \frac{H_y^{(p)}}{H_y^{(p,i)}^2}}.`$ (16)
Even if the enhancement factor thus defined closely resembles the correct total electric field enhancement factor, we are evidently losing information about the different electric field components, in turn relevant to the SERS polarization selectivity.

### B $`p`$ polarization: Electric field

In order to obtain the electric field components from the $`y`$ component of the magnetic field, use can be made of the Maxwell equation

$$\times 𝐇=ı\frac{\omega }{c}ϵ𝐄.$$
(17)
In the incident medium, Eq. (8) provides the only non-zero component of the magnetic field. Use of Maxwell’s equation (17) leads to the following electric field components:

$`E_x^{(p,>)}(𝐫)=E_x^{(p,i)}(𝐫)ı{\displaystyle \frac{c}{4\pi \omega ϵ^>}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}[H_y^{(p,>)}(𝐫^{}){\displaystyle \frac{}{}}`$ (20)
$`\times {\displaystyle \frac{^2G^>(𝐫,𝐫^{})}{zn^{}}}{\displaystyle \frac{G^>(𝐫,𝐫^{})}{z}}{\displaystyle \frac{H_y^{(p,>)}(𝐫^{})}{n^{}}}]`$
$`E_y^{(p,>)}(𝐫)=0`$ (21)
$`E_z^{(p,>)}(𝐫)=E_z^{(p,i)}(𝐫)+ı{\displaystyle \frac{c}{4\pi \omega ϵ^>}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}[H_y^{(p,>)}(𝐫^{}){\displaystyle \frac{}{}}`$ (23)
$`\times {\displaystyle \frac{^2G^>(𝐫,𝐫^{})}{xn^{}}}{\displaystyle \frac{G^>(𝐫,𝐫^{})}{x}}{\displaystyle \frac{H_y^{(p,>)}(𝐫^{})}{n^{}}}]`$
These equations can be rewritten in terms of the source functions $`H(x)`$ and $`L(x)`$ as follows:

$`E_x^{(p,>)}(𝐫)=E_x^{(p,i)}(𝐫){\displaystyle \frac{\omega }{4c}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}\{H(x^{})`$ (28)
$`\times [{\displaystyle \frac{z\zeta (x^{})}{𝐫𝐫^{}^2}}(𝐧(𝐫𝐫^{}))H_2^{(1)}(k_ϵ𝐫𝐫^{})`$
$`{\displaystyle \frac{1}{\gamma ^{}k_ϵ𝐫𝐫^{}}}H_1^{(1)}(k_ϵ𝐫𝐫^{})]`$
$`L(x^{}){\displaystyle \frac{z\zeta (x^{})}{\gamma ^{}k_ϵ𝐫𝐫^{}}}H_1^{(1)}(k_ϵ𝐫𝐫^{})\}`$
$`E_y^{(p,>)}(𝐫)=0`$ (29)
$`E_z^{(p,>)}(𝐫)=E_z^{(p,i)}(𝐫){\displaystyle \frac{\omega }{4c}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\gamma ^{}dx^{}\{H(x^{})`$ (33)
$`\times [{\displaystyle \frac{xx^{}}{𝐫𝐫^{}^2}}(𝐧(𝐫𝐫^{}))H_2^{(1)}(k_ϵ𝐫𝐫^{})`$
$`{\displaystyle \frac{\zeta ^{}(x^{})}{\gamma ^{}k_ϵ𝐫𝐫^{}}}H_1^{(1)}(k_ϵ𝐫𝐫^{})]`$
$`+L(x^{}){\displaystyle \frac{xx^{}}{\gamma ^{}k_ϵ𝐫𝐫^{}}}H_1^{(1)}(k_ϵ𝐫𝐫^{})\}`$
where the explicit form of the Green’s function has been taken into account, leading to the appearance of 1st and 2nd order Hankel functions of the first kind $`H_1^{(1)},H_2^{(1)}`$. For the Gaussian incident field given by Eq. (II A), the electric field components are

$`E_x^{(p,i)}(𝐫)=`$ $`{\displaystyle \frac{ı}{n_c^>}}H_y^{(p,i)}(𝐫)[{\displaystyle \frac{}{}}ı\mathrm{cos}\theta _0(1+w(x,z))`$ (37)
$`(ı{\displaystyle \frac{4}{k_ϵ^2W^4}}(x\mathrm{sin}\theta _0z\mathrm{cos}\theta _0)`$
$`{\displaystyle \frac{2}{k_ϵW^2}})\mathrm{sin}\theta _0(x\mathrm{cos}\theta _0+z\mathrm{sin}\theta _0)]`$
$`E_y^{(p,i)}(𝐫)=`$ $`0`$ (38)
$`E_z^{(p,i)}(𝐫)=`$ $`{\displaystyle \frac{ı}{n_c^>}}H_y^{(p,i)}(𝐫)[ı\mathrm{sin}\theta _0(1+w(x,z))`$ (41)
$`+(ı{\displaystyle \frac{4}{k_ϵ^2W^4}}(x\mathrm{sin}\theta _0z\mathrm{cos}\theta _0)`$
$`{\displaystyle \frac{2}{k_ϵW^2}})\mathrm{cos}\theta _0(x\mathrm{cos}\theta _0+z\mathrm{sin}\theta _0)]`$
Equations (II B) and (II B) provide the electric field components in the incident medium of the resulting $`p`$-polarized EM field, incident plus scattered from the rough surface. The scattered electric field involves an additional surface integral in terms of the source functions, previously obtained (numerically) from the above mentioned coupled integral equations. Analogous expressions, not shown here, for the corresponding electric field components inside the metal can be obtained from Eq. (11). On the other hand, recall that a similar procedure can be straightforwardly developed to yield the magnetic field components in the case of $`s`$-polarized EM waves as surface integrals in terms of the surface electric field ($`y`$ component) and its normal derivative.

### C $`p`$ polarization: Normal and tangential surface electric field

It should be pointed out that when trying to evaluate the electric field on the surface, or even very close to it, from Eqs. (II B), non-integrable singularities appear associated with the Green’s functions derivatives for vanishing arguments. Use of expressions (II B) for the evaluation of the electric field close to the surface will produce unphysical results. To deal properly with this situation, more care should have been taken in doing the derivatives of the integral (8) describing the magnetic field, whose integrand already exhibits singularities, though integrable. A simple way to work around this problem consists of evaluating the electric field at the surface itself, and we have found very simple relations connecting the normal and tangential components of the electric field (see Fig. 1) with the surface magnetic field and its normal derivative \[cf. Eqs. (II A)\]. These are

$`E_n^{(p,>)}(x)`$ $`=`$ $`{\displaystyle \frac{ıc}{\omega ϵ^>}}\gamma ^1{\displaystyle \frac{dH(x)}{dx}}`$ (43)
$`E_t^{(p,>)}(x)`$ $`=`$ $`{\displaystyle \frac{ıc}{\omega ϵ^>}}\gamma ^1L(x).`$ (44)
These expressions are extremely useful, for they facilitate considerably our study of the SERS EM mechanism.

Thus, taking advantage of the one-dimensional scattering geometry (which, although is not general, is not inappropriate to study the SERS EM mechanism ), we only have to deal with the y-component of the magnetic field to obtain a simplified, basically exact solution to the scattering problem from the classical EM viewpoint, at the excitation frequency. The numerical solution of the resulting integral equations yields the surface magnetic field and its normal derivative as the main results. The drawback of working with the magnetic field, when the quantity of interest is the electric field, is avoided by expressions (II C), which allows us to obtain the surface electric field with the only additional algebra of calculating a spatial derivative.

We now properly define the electric field enhancement factors for either a single component or the total field as the normalized intensities:

$`\sigma _\alpha (\omega )`$ $`=`$ $`{\displaystyle \frac{E_\alpha ^{(p,>)}^2}{E^{(p,i)}^2}},`$ (45)
$`\sigma (\omega )`$ $`=`$ $`{\displaystyle \frac{E^{(p,>)}^2}{E^{(p,i)}^2}},`$ (46)
with $`\alpha =n,t,x,z`$ and $`E=𝐄^2=E_n^2+E_t^2=E_x^2+E_z^2`$.

### D Numerical implementation

The numerical procedure has been implicitly outlined above; further details have been given in Ref. . Self-affine random fractal surfaces numerically generated by means of Voss’ fractional Brownian motion algorithm  are studied. This kind of fractals exhibit self-affine scaling properties in a broad spatial range, and have properties that resemble those of some SERS substrates, such as cold-deposited metal films or etched metal surfaces. In the numerical calculations, surface realizations of length $`L=10.29\mu `$m, consisting of $`N_p=n_iN`$ sampling points obtained by introducing $`n_i=4,6,8`$ or 10 cubic-splined interpolating points into a sequence of $`N=201`$ points extracted from each generated fractal profile with $`N_f=1024`$ points; note the considerably larger sample density with respect to that of Ref. . The statistical properties of the physical quantities of interest will be calculated on the basis of Monte Carlo simulations for an ensemble of fractal realizations.

## III Local Field Enhancement: Hot Spots

We now turn to the investigation of the occurrence of very large near EM field enhancements. Particularly, we will concentrate on self-affine fractals with Hurst exponent $`H=0.1`$ (namely, local fractal dimension $`D_f=2H=1.9`$), which have been shown in Ref.  to give rise to large surface magnetic fields. The lower scale cutoff $`\xi _L50`$ nm has been chosen to resemble that of SERS substrates. The upper scale cutoff, typically $`\xi _L50\mu `$m, is considerably larger than the illuminated area $`L`$, and this is in turn sufficiently (in order to avoid finite length effects) larger than the incoming wavelength (0.4 $`\mu `$m$`<\lambda <1.3\mu `$m). Thus physical scaling is meaningful for the relevant interval of this scattering problem. The effect of further reducing the lower scale cutoff will be investigated elsewhere; in this regard, it should be recalled that the minimum scale relevant to the far field pattern has been studied for Koch fractals.

### A Near Field Intensity Maps

In Fig. 2, the intensity (on a logarithmic scale) of the electric and magnetic near field in the vicinity of a self-affine Ag surface with $`D=1.9`$ and rms height $`\delta =514.5`$ nm, in a particular region (of about $`1\times 1\mu `$m2), is shown for normal incidence with light of wavelength $`\lambda =514.5\mu `$m; in addition, the intensities of the two different components ($`x`$ and $`z`$) of the electric field are separately shown. Before analyzing the results, some comments are in order with regard to the numerical calculations. Whereas the electric and magnetic fields in vacuum away from the interface are given by the integral equations (II B) and (8), respectively, and similarly for the EM field inside silver, their corresponding values on the interface are directly obtained from the source functions through Eqs. (II A) and (II C). As mentioned in the preceding section, Eqs. (8) and (II B) exhibit singularities upon approaching the surface, so that they are not accurate at points very close to the surface, typically within distances smaller than the surface sampling interval. Thus the EM field intensity in Fig. 2 at distances from the surface smaller than $`4L/N_p`$ for the magnetic field and $`6L/N_p`$ for the electric field are obtained from the weighted values on the two closest sampling points on the surface, taking explicitly into account for points inside silver the continuity conditions for the magnetic \[Eq. (II A)\] and electric field components. Despite that, some slight (not inherent) mismatch might still appear when entering into the surface field area, mostly in the intensities of the electric field components.

It is, of course, expected that such continuity conditions across the interface could be roughly observed in the calculations even if we were not considering the above explicit matching. On the one hand, the continuity of the tangential component of the magnetic field is neatly seen in Fig. 2(d); on the other hand, the continuity of the tangential electric field is appreciable in Fig. 2(b) and (c) through the continuity of the $`x`$ (respectively, $`z`$) component of the electric field at locally flat (respectively, vertical) parts of the rough surface, whereas the discontinuity of the normal component of the electric field (continuity of the normal component of the displacement vector) is inferred from the discontinuity of the $`x`$ (respectively, $`z`$) component of the electric field at locally vertical (respectively, flat) areas. Incidentally, note also that the EM field inside silver decays very rapidly as expected from the Ag skin depth $`d=(c/\omega )(ϵ^<)^{1/2}27`$ nm (cf. Ref.  for the Ag dielectric constant).

It is evident from Fig. 2 that the maximum local EM fields are located right on top of the Ag surface, whereupon some particularly bright spots appear. Thus we next plot in Fig. 3 the intensities of the surface EM fields (including electric tangential and normal components) for the surface area shown in Fig. 2, including the surface profile. Very narrow peaks surrounded by dark areas are observed, whose widths are well below half wavelength of the SPP ($`\lambda _{SPP}/2=243`$ nm). Some of these peaks can be considered as optical excitations (hot spots) where very large local FE’s occur, such that non linear optical processes would be strongly enhanced therein. In particular, the largest in Fig. 3, shown in the inset, is of the order of $`\sigma 610^2`$. Note that our calculations identify the local electric field component that is responsible for such FE: the normal component. In Fig. 4, the near electric field in the vicinity of this hot spot (zooming in Fig. 2(a)) is given; interestingly, it is associated with a surface peak.

It has been experimentally shown by near-field microscopy that such optical excitations rapidly disappear upon changing the frequency of the incident radiation. Our rigorous calculations corroborate those experimental observations, as seen in Fig. 5, where the surface electric field intensity is plotted for several incident wavelengths close to $`\lambda =514.5`$ nm. For the wavelengths $`\lambda =`$495.9 and 539.1 nm, the high intensity spot is still clearly visible, though less bright. It fades away, however, for larger frequency shifts, becoming barely visible for $`\lambda >563.6`$ nm or $`\lambda <476.9`$ nm (about 10$`\%`$ frequency shift).

The hot spot shown in Fig. 4 is strongly polarized along the normal to the surface. This is extremely relevant to SERS spectroscopy, since it might impose selection rules to the vibrational modes of the adsorbed molecule. Is it possible to find hot spots with different polarizations? Only in the case of silver at wavelengths close to the surface plasma wavelength, we have found certain spots strongly polarized along the tangential direction too, though weaker than those normally polarized. These tangentially polarized hot spots exhibit FE factors not larger than $`\sigma _t10^2`$, and are typically located in regions presenting larger values of $`\sigma _n`$: this will be discussed elsewhere.

The occurrence of local optical excitations has been studied for the same self-affine surface profile illuminated with different excitation wavelengths, and also for metals such as Au and Cu. Although not shown here, similar normally polarized hot spots are found on a broad spectral range on all the self-affine surfaces of Ag, Au, and Cu. We now discuss some of the characteristics of these hot spots.

It has been argued  that these hot spots are due to Anderson localization of SPP. Theoretical works based on a dipolar model, also demonstrate the possibility of creating strongly confined and intense excitations on self-affine fractal surfaces; in fact, in the case of random metal-dielectric films, Anderson localization of surface plasmon modes is predicted. Our numerical calculations, not subject to dipolar (and quasi-static) restrictions, indeed reveal the existence of these kind of optical excitations and, although compatible with the possibility of them being due to Anderson localization of SPP, do not permit to draw further conclusions in this respect. Our scattering geometry involving the interaction between a propagating beam of light and a metal surface does not lend itself well for the characterization of the SPP Anderson localization phenomenon. To that end, the study of the propagation and transmission of SPP through rough surfaces would be more adequate. Only the fact that light can couple into these, possibly localized, SPP modes through the roughness can be inferred from our calculations and from the typical PSTM configurations.

It is also worth pointing out that the extremely rapid decay and widening of the optical excitations in the near field (see Figs. 2 and 4) implies that PSTM images taken at certain distance from the surface, leaving aside the rounding effects of the tip, will manifest themselves as much wider and weaker optical excitations, and this seems to be the case.  Direct probing of hot spots, on the other hand, could be carried out by the nonlinear effects of physi- or chemi-sorbed molecules.

### B Spectral dependence

In order to analyze the polarization and spectral dependencies of the optical excitations, we present in Fig. 6 the maximum local FE values found at Ag, Au, and Cu self-affine surfaces with $`D=1.9`$ and $`\delta =514.5`$ nm, obtained from numerical calculations of the surface EM field for an ensemble of $`N_r=60`$ realizations generated as mentioned above (only the data from the central half of each realization are used). The results for weaker self-affine surfaces ($`\delta =102.9`$ nm) used in Ref.  to compute magnetic FE’s are also shown. In addition, the case of having water as incident medium (solvent) has been analyzed, though in Fig. 6 only the results for H2O/Ag are shown. Several remarks are in order with regard to Fig. 6.

Very large FE’s appear for a wide spectral range covering the visible and entering into the NIR. In the red and NIR parts of the spectrum, all three metals being studied behave similarly, giving rise to hot spots exhibiting strongly enhanced electric field intensities normal to the surface, the tangential component tending to vanish. This behavior can be understood in accordance with the spectral evolution of the Ag, Au, and Cu dielectric constants, all showing increasingly large negative real parts, tending to the perfectly conducting limit $`ϵ\mathrm{}`$ that predicts vanishing tangential electric fields. For wavelengths $`\lambda 1240`$ nm, however, FE’s are expected to slowly decrease as the surface is “seen” by the incoming radiation of increasing wavelength as increasingly flatter. Other calculations, not shown here, indicate that this is the case for $`\lambda >2\mu `$m. In fact, this decay is observed at lower wavelengths (within the spectral interval covered by Fig. 6) for the fractal surface with $`\delta =102.9`$ nm.

The optical responses of Au and Cu manifest significant differences with respect to that of Ag in the blue part of the spectrum. The onset of interband transitions, which takes place in Ag at $`\lambda 300`$ nm unlike in Au and Cu (slightly below 600 nm), makes the difference inasmuch as such transitions constitute a strong absorption mechanism. Consequently, FE’s should be significantly reduced for wavelengths below the onset threshold, as is evident in Fig. 6 for Au and Cu below $`\lambda 600`$ nm, but not seen for Ag since the lower wavelength considered in Fig. 6 is above the Ag onset threshold. Moreover, silver surfaces at small incoming wavelengths approaching the surface plasmon wavelength (but above the onset of interband transitions) present strong local optical excitations tangentially polarized, as mentioned above. These tangential-electric hot spots can lead to local FE nearly comparable to those corresponding to the normal-electric ones at such wavelengths, although more than an order of magnitude weaker than those obtained at larger wavelengths \[see Fig. 6(b) and (c)\]. This has important implications in SERS, since Ag substrates of the kind studied here, when illuminated at wavelengths $`\lambda <600`$ nm, could enhance the Raman signal coming from a vibrational mode of the molecule sensitive to the tangential electric field, as well as those sensitive to the normal electric field (typically established as predominant according to SERS selection rules ).

Finally, note that using water as solvent does not introduce significant changes in the qualitative and quantitative behavior of the maximum local FE. In addition, we would like to point out that the magnetic FE follows qualitatively (and almost quantitatively) the normal electric FE.

## IV Statistical properties of the surface FE

In this section we study the statistical properties of the FE’s occurring on the surface of self-affine fractal profiles with fractal dimension and rms roughness as in the preceding section. These properties are obtained from Monte Carlo numerical simulations results performed as described in Sec. II D. Preliminary results based on magnetic field calculations for weaker fractal surfaces have been presented in Ref. .

### A Probability Density Function

In Sec. III we have found, by direct observation of the calculated near-field excited in the neighborhood of the interface, that very large fields can be excited at the surface. These enhanced excitations are coupled to the surface by the surface roughness and, thus, depend strongly on its properties. The question then arises as to how probable these values are and, in turn, how the EM field intensity is distributed over the fractal surface. In Fig. 7, we show the PDF of the surface EM FE (including separately tangential and normal components) for self-affine fractal surfaces with $`D=1.9`$ and $`\delta =`$514.5 nm for two incoming wavelengths $`\lambda =`$ 514.5 and 1064 nm, averaging over various angles of incidence. For the sake of comparison, the results for fractal surfaces with smaller rms height $`\delta =`$102.9 nm, and also with both smaller fractal dimension $`D=1.2`$ and rms height $`\delta =`$102.9 nm are included (the latter as used in Ref.  for magnetic FE calculations).

For the smoother surface, the resulting PDF is a narrow distribution centered at the surface EM FE value for a flat metal surface, as expected for such a weakly rough surface and in agreement with Refs. , wherein similar results were shown for the intensity of the magnetic field (typically, $`\sigma _H|1+R|^2`$, $`R`$ being the corresponding Fresnel reflection coefficient). Recall that the electric field components on flat metal surfaces follow $`\sigma _t\mathrm{cos}^2\theta _0|1+R|^2`$ and $`\sigma _n\mathrm{sin}^2\theta _0|1R|^2`$, so that the PDF distributions for $`D=1.2`$ in Fig. 7 are only significant for FE’s approximately in between the minimum and maximum expected values for the different $`\theta _0`$.

For rougher surfaces, however, the surface EM no longer resembles the flat surface result, presenting alternating dark and bright regions, and giving rise with increasing roughness parameters to very bright hot spots surrounded by large dark regions. The corresponding PDF becomes wider, turning into a slowly decaying function that is maximum at zero and exhibits a long tail for large FE values (see Ref.  for the magnetic FE PDF for the self-affine fractal with $`D=1.9`$ and $`\delta =`$102.9 nm). Upon comparing Figs. 7 (a) and (b), it is evident that moderately large $`\sigma _t`$’s become feasible at $`\lambda =`$514.5 nm as well as very large $`\sigma _n`$’s, whereas only the $`\sigma _n`$’s are expected to be intense at $`\lambda =`$1064 nm.

### B Average and Fluctuations

As a result of the change in the surface EM field PDF for increasing surface roughness parameters, the moments of the distribution are also modified. Particularly relevant are the average and the statistics of the fluctuations, as they can also provide some information about the global response of larger surfaces (of the order of centimeters) under broad beam illumination.

In Fig. 8 we present the spectral dependence of the mean FE for the same self-affine fractal surfaces whose local FE where shown in Fig. 6. In fact, the qualitative behavior of the mean FE does not differ substantially from that exhibited by the maximum local FE. Basically, there is a broad excitation spectra for the rougher fractals (slightly narrower for the smoother fractal), covering the visible and the NIR (at least up to $`\lambda =2\mu `$m). The behavior is similar for Ag, Au, and Cu and, in all cases, the electric field is predominantly normal to the surface. The blue part of the excitation spectra reveals, on the other hand, a rapid decrease for Au and Cu associated with the onset of interband transitions, whereas large normal FE’s can still be found in that spectral region for Ag fractals as well as an increase of the tangential electric field upon approaching the surface plasma wavelength (but still above the threshold wavelength of interband transitions). The latter blue tangential electric FE increase is slightly larger when water rather than air constitutes the propagating medium. On the other hand, it should be emphasized that the rougher fractal surfaces used in Fig. 8 give rise to an estimated SERS FE factor $`𝒢_{EM}10^5`$, in good agreement with the phenomenological factor experimentally induced.

The absorption spectra are shown in Fig. 9 for the sake of comparison. It is evident that the absorption spectrum does not resemble the qualitative behavior of the excitation spectrum in Fig. 8. Therefore, for this kind of fractal surfaces yielding wide excitation and absorption spectra, the maximum absorption region as experimentally obtained from the absorption spectra cannot be straightforwardly related to the optimum excitation wavelength. Furthermore, it should be emphasized that strong absorption can even be associated with very low surface EM fields (and thus the substrates being SERS inactive), as is the case of Au and Cu self-affine fractals in the blue spectral region. In other substrate configurations, nonetheless, the contrary might be the case, and absorption bands can be used to identify excitation bands resulting in strong surface FE (substrates becoming SERS active), as in rough surfaces presenting surface shape plasmon resonances, such as colloidal aggregates, or in gratings diffracting into propagating SPP. In summary, one has to carefully interpret absorption spectra when using such information to determine the appropriate SERS (or any other surface optical nonlinear effect) excitation frequency.

Finally, we show in Fig. 10 the spectral dependence of the FE fluctuations. In accordance with the previously discussed FE PDF widening with increasing surface roughness, it is obvious that the rougher the surface, the larger the fluctuations. And for sufficiently rough surfaces, the fluctuations can be even larger than the average, as seen upon comparing Fig. 10 with Fig. 8. Actually, the FE fluctuations rather than the average provide a good estimate of how probable and how bright hot spots are. Indeed, the spectral dependence of the FE fluctuations in Fig. 10 closely follows that of the maximum local FE shown above in Fig. 6.

## V Conclusions

By means of a rigorous Green’s theorem integral equation formalism, we have studied the occurrence of strong local optical excitations (hot spots) on self-affine fractal surfaces of Ag, Au, and Cu. The statistics of the surface field fluctuations that produce these strong excitations have also been studied. The formalism exploits the scalar character of the resulting integral equations for one-dimensional surfaces illuminated with linearly $`s`$\- or $`p`$-polarized light, by treating the problem in terms of the electric or magnetic field, respectively. In the case of $`p`$ polarization, which is the relevant one in our problem due to the SPP excitation selectivity, we have calculated the electric field from the resulting only nonzero component of the magnetic field and its normal derivative on the surface. The problem is studied numerically by means of Monte Carlo simulations of the interaction of light with self-affine metal fractals whose profiles were obtained from the trace of a fractional Brownian motion. The appearance of hot spots and their statistics have been determined for a broad spectral range of the incoming light (400 nm $`<\lambda <`$ 1300 nm).

We have found hot spots on self-affine fractals with fractal dimension $`D=1.9`$ and rms height $`\delta =`$514.5 nm. These hot spots constitute very strong and narrow (considerably narrower than half of the SPP wavelength) surface EM field excitations, with very selective excitation spectra (both temporally and spatially). Typically, they give rise to local FE strongly polarized along the normal to the rough surface. The largest ones we have found yield local SERS FE factors of $`𝒢_{EM}10^7`$, and appear for a wide range of incoming wavelengths covering the visible and NIR up to $`\lambda 2\mu `$m. Interestingly, our results reveal that weaker, tangentially polarized hot spots can be found in Ag fractals for small excitation wavelengths (blue or smaller).

All these features are not incompatible with the suggestion that Anderson localization of SPP is the underlying physical mechanism responsible for such optical excitations (e.g. the exponential decay of the SPP transmission versus rough surface length could certainly provide some indication of such mechanism ), although no direct evidence of this can be obtained from our scattering geometry and calculations.

The PDF of the surface EM field for those self-affine fractals exhibiting hot spots is a slowly decaying function with a significant tail for large surface FE’s. It differs substantially from that for smooth surfaces for which the PDF is a narrow distribution centered at the value of the EM field expected on a flat metal surface.

The mean FE acquires considerable values in a broad spectral region. For the three metals considered, in most of the visible and NIR ($`\lambda <2\mu `$m) excitation regions, the component responsible for this enhancement is the electric field normal to the surface. We have found that, for these self-affine fractal surfaces the average SERS FE factors are $`𝒢_{EM}10^5`$.

We have also analyzed the spectral dependence of the surface FE fluctuations. Such fluctuations are indeed very large for the self-affine fractals that give rise to hot spots, and present a qualitative behavior similar to that of the maximum local FE in the vicinity of the hot spots. This is an interesting property that could be used in, e.g., PSTM studies to identify samples with the potential capability of yielding large optical excitations: even if no hot spots are found in the region being scanned, the calculated fluctuations of the resulting intensity map could provide a statistical account on the probability of finding hot spots (simpler than calculating the total PDF for which much more data from a larger scanning area would be required).

With regard to the quantitative aspects of SERS, the maximum local enhancement factors are still below those that could be deduced from experimental works on single molecule detection  and from approximate theoretical calculations.  This, however, is not entirely surprising, given the differences in the type of SERS substrates being considered. In fact, the typical SERS spectroscopy enhancement factors are fairly similar to those found in this work.

We would like to emphasize that our formulation is exact within the classical EM framework, at least as far as the linear (direct) field enhancement factor is concerned. Thus, our calculations should provide a truthful picture of the linear optical response of self-affine fractal metal substrates. Further work is of course needed to test the perhaps naive assumption that the enhancement factor at the Raman-shifted frequency is identical to that obtained at the excitation frequency. Also, more work is required to study the effects of lower scaling cutoffs in the generation of the fractal surfaces, as these spatial frequencies might contribute to build up the enhancement factors.

Finally, we mention that work involving rigorous calculations of the kind presented here is also in progress for the study of the field enhancements produced on self-similar substrates, such as those found in colloidal aggregates, and for the study surfaces covered by a monolayer of Raman active molecules (Langmuir-Blodgett films).

###### Acknowledgements.

This work was supported by the Spanish Dirección General de Enseñanza Superior e Investigación Científica y Técnica, through grant No. PB97-1221. We also thank the Mexican-Spanish CONACYT-CSIC program for partial travel support.