text:
# Flattening and truncation of stellar discs in edge-on spiral galaxies

## 1 Introduction

The study of the light and mass distribution of galactic discs is of fundamental importance for our understanding of the formation and evolution of spiral galaxies. It is well established that globally, the radial light distribution of galactic discs is exponential in nature (Freeman 1970), over a wide range of central surface brightness (de Jong 1996b; de Blok, van der Hulst, & Bothun 1995). This radial exponential distribution can be understood in the framework of the collapse theory of disc galaxy formation if the stellar mass-to-light ratio is approximately constant with radius (Fall & Efstathiou 1980; van der Kruit 1987; Dalcanton, Spergel, & Summers 1997). Perpendicular to the galactic plane, the light distribution of the disc is also approximately exponential (van der Kruit & Searle 1981a, b, 1982, hereafter KS1–3; Wainscoat, Freeman, & Hyland 1989; Shaw & Gilmore 1990; de Grijs, Peletier, & van der Kruit 1997). While the radial structure of the stellar disc is probably determined during disc galaxy formation, the vertical structure largely reflects the secular disc evolution (Wielen 1977; Bahcall 1984; Lacey 1984). Together, the radial and vertical disc structure determine the intrinsic disc flattening. Accordingly, disc flattening is governed by both disc galaxy formation and evolution and may offer additional insight into the physics of these.

The disc flattening can only be studied accurately in edge-on galaxies, where the radial and vertical distributions are essentially independent (KS1). Previous studies of large samples of edge-on disc galaxies pointed out that, on average, the disc light distribution slowly flattens from earlier to later morphological type (Guthrie 1992; Kudrya et al. 1994). These studies adopted a statistical approach, using isophotal axial ratios. A similar trend is observed in a sample of 47 edge-ons (de Grijs 1998), using the ratio of the exponential scalelength to the exponential scaleheight of the disc light (see also Schwarzkopf & Dettmar 2000). In this sample, a further comparison of the disc flattening with other global parameters, such as rotational velocity or H i content, may reveal additional trends concerning disc galaxy formation and evolution. Besides this, the flattening of the disc light can be used to estimate the disc contribution to the rotation curve (Bottema 1993). This allows one to study the dark matter content of spiral galaxies and provides an independent test of the so-called maximum disc hypothesis (van Albada et al. 1985; Freeman 1992).

A related study of the stellar kinematics of galactic discs by Bottema (1993) shows that more massive spirals tend to have larger stellar velocity dispersions. Together with the increase of the stellar velocity dispersion with age (Wielen 1977), this suggests that more massive discs are more dynamically evolved. If true, then more massive discs are also expected to have larger scaleheights. This can easily be verified observationally in edge-on spiral galaxies.

Another global parameter describing the radial light distribution of discs is its outer edge. This truncation of the disc light occurs at a radius $`R_{\mathrm{max}}`$, which varies from 3 to 5 scalelengths among disc galaxies (van der Kruit 2001). It is most easily found in edge-on spirals because of their higher surface brightness (KS1–3; Barteldrees & Dettmar 1994; Pohlen et al. 2000b; Schwarzkopf & Dettmar 2000; de Grijs, Kregel, & Wesson 2001). In less inclined spirals, for which azimuthally averaged radial light profiles are often studied, the non-axisymmetric component (e.g. due to spiral structure) tends to smooth out a truncation present in the old disc. This effect was first noted by van der Kruit (1988); for 16 face-on spirals, of which 15 did not show any sign of a truncation in the radial light profile, he found that the three outermost isophotes were much more closely spaced than the inner ones, providing clear evidence for truncations. The origin of the truncation is still unclear. If the truncation is also present in the disc mass distribution, it may be caused by the angular momentum distribution of the proto galaxy (van der Kruit 1987). This would imply that H i observed beyond the stellar disc (e.g. Broeils & Rhee 1997) has been accreted (cf. Bottema 1996). Another possibility is that the truncation corresponds to the radius at which the gas density drops below a threshold density necessary for star formation (Fall & Efstathiou 1980; KS3; Seiden, Schulman, & Elmegreen 1984; Kennicutt 1989).

In a previous study, de Grijs (1998) analysed BVIJK photometry of a statistically complete sample of nearby edge-on disc galaxies. He used profiles taken at suitably chosen positions to determine both the radial and vertical scale parameters (hereafter 1D method). The determination of the radial scalelengths was based on radial profiles, taken at some distance above the plane to minimize extinction. Still, the uncertainty in this analysis is considerable. First, spiral galaxies may have a surface brightness distribution which shows significant local deviations from assumed smooth model distributions (Shaw & Gilmore 1990; de Jong 1996a). Since the 1D method uses a small portion of the galaxy, it is susceptible to these local deviations. Secondly, the 1D method does not consider the presence of a radial truncation. Finally, an exponential was used as a fitting function, thereby neglecting the effect of line of sight projection. These omissions have probably caused a serious error in the disc scalelengths.

Here, we re-analyse the I-band photometry of a sub-sample of de Grijs’s sample, consisting of 34 spiral galaxies. In the first part of this paper, we use a two-dimensional (2D) least-squares algorithm to study the exponential scale parameters, and thereby the disc flattening. This 2D method minimizes the effects of a radial truncation and takes the line of sight projection into account. In addition, it allows for a bulge-disc decomposition and adopts a mask to minimize the effect of dust extinction. To assess the effect of dust extinction outside this mask, the B-band data are subjected to the same analysis. The 2D method was briefly introduced by de Grijs et al. (2001), where the results of a detailed study of the radial truncations of four spiral galaxies were presented. In the second part of this paper we report on a similar study of the radial truncations for our sample in the I band. Finally, the observed truncation radii are compared to those predicted by the scenarios proposed for the origin of the truncation. This paper is organized as follows. In Sect. 2 we present the sample. The 2D least-squares method is discussed and tested in Sect. 3. In Sect. 4, the results of its application to the I-band photometry are presented, together with the results of our study of the disc truncations. We further analyse and discuss the results in Sect. 5. Our main findings are summarized in Sect. 6.

## 2 The Sample

The parent sample (de Grijs 1998) contains 47 spirals and lenticulars and was selected from the ESO-LV catalogue (Lauberts & Valentijn 1989) according to the following criteria:

| (i) | inclination, i $``$ 87; |
| --- | --- |
| (ii) | blue diameter, D$`{}_{25}{}^{}{}_{}{}^{B}`$ $`>`$ 2.2; |
| (iii) | Hubble type ranging from S0–Sd; |
| (iv) | non-interacting |

From this sample we selected all regular spiral galaxies and excluded those that are clearly lopsided or warped. Important properties of the 34 selected spirals are summarized in Table 2, columns (1)–(6) ∗\*∗\*$``$All distance dependent parameters in this paper are calculated using $`H_0`$=75 km s-1 Mpc-1. We note that although the sub-sample covers a large range in maximum rotational velocities ($`v_{\mathrm{max}}`$ = 50–400 km s-1), it is dominated by spirals of intermediate- to late-type.

We checked the completeness of the sub-sample with the $`V/V_{\mathrm{max}}`$ test (Davies 1990). Based on the blue angular diameters given by the ESO-LV catalogue, the average $`V/V_{\mathrm{max}}`$ is $`0.58\pm 0.06`$. This is slightly larger than the average of 0.5 expected for a randomly distributed sample. In the remainder of this paper, one should keep in mind that our sample is of modest size and may not be statistically complete.

While our J and K-band data more closely resemble the underlying disc mass distribution than the optical BVI data (Rix & Rieke 1993), the quality of the former is poorer (de Grijs 1998) and will not be used. Of the optical passbands, the light in the I band is the least affected by dust extinction and the most dominated by the old stellar population (e.g. Charlot & Bruzual 1991). We will therefore primarily use those data to study the global structure of stellar discs. The data reduction and sky subtraction have been discussed by de Grijs (1997, 1998).

## 3 The two dimensional method

### 3.1 The model

The spiral galaxy luminosity density is assumed to be composed of a disc and a bulge component. The disc is assumed to be axisymmetric, transparent and, in first instance, without a radial truncation:

$$L(R,z)=L_0e^{R/h_\mathrm{R}}e^{z/h_\mathrm{z}},$$
(1)
where (R, z) are the cylindrical coordinates, $`L_0`$ is the central luminosity density and $`h_\mathrm{R}`$ and $`h_\mathrm{z}`$ are the radial scalelength and the vertical scaleheight (in terms of the vertical scale parameter $`z_0`$; $`h_\mathrm{z}`$ = 0.5 $`z_0`$). While other functionalities have been proposed for the vertical distribution, most notably the locally isothermal (KS1) and sech$`(z)`$ distributions (van der Kruit 1988), the exact choice is unimportant for our purpose (Sect. 3.2). The vertical exponential behaviour is taken to be independent of radius following de Grijs & Peletier (1997). They showed that, in the discs of most of the spirals in our sample (at least all late-type galaxies, which is about two thirds), the scaleheight does not increase with radius. It is further assumed that the galaxies are perfectly edge-on. As KS1 demonstrate, the error of the scale parameters due to this assumption is negligible as long as i $``$ 86 (see also de Grijs et al. 1997), which is the case for all the galaxies in our sample. With these assumptions the disc surface brightness distribution is then given by (KS1):

$$\mathrm{\Sigma }_{\mathrm{disc}}(R^{},z)=\mathrm{\Sigma }_0(R^{}/h_\mathrm{R})K_1(R^{}/h_\mathrm{R})e^{z/h_\mathrm{z}},$$
(2)
where R is the projected radius along the major axis, $`\mathrm{\Sigma }_0=2h_\mathrm{R}L_0`$ is the projected edge-on central surface brightness and $`K_1`$ is the modified Bessel function of the first order.

The surface brightness distribution of the bulge is modelled either by an exponential:

$$\mathrm{\Sigma }_{\mathrm{bulge}}(R^{},z)=5.360\mathrm{\Sigma }_\mathrm{e}e^{1.679r/r_\mathrm{e}},$$
(3)
or by an r1/4 law:

$$\mathrm{\Sigma }_{\mathrm{bulge}}(R^{},z)=\mathrm{\Sigma }_\mathrm{e}e^{7.669[(r/r_\mathrm{e})^{1/4}1]},$$
(4)
here $`r`$ is the projected radius; $`r=\sqrt{R^2+(z/q)^2}`$, $`r_\mathrm{e}`$ is the effective radius, $`\mathrm{\Sigma }_\mathrm{e}`$ the effective surface brightness, and q the bulge axial ratio (the ratio of the bulge minor to major axis diameter). The best-fitting bulge, in the reduced $`\chi ^2`$ sense, is adopted. To account for the effect of seeing, the model bulges are convolved with a two-dimensional Gaussian PSF. For the disc component the seeing effect can be neglected (Sect. 3.3). The central positions and position angles of the galaxies, determined by fitting ellipses to the I-band images (Sect. 4.1), are fixed parameters in the fit. Thus, the adopted model has at most six free parameters; $`\mathrm{\Sigma }_0`$, $`h_\mathrm{R}`$, $`h_\mathrm{z}`$, $`\mathrm{\Sigma }_\mathrm{e}`$, $`r_\mathrm{e}`$, and $`q`$.

### 3.2 The approach

Since the studies of van der Kruit & Searle (KS1–3), several techniques have been developed to model the global structure of stellar discs in edge-on spiral galaxies (Kylafis & Bahcall 1987; Shaw & Gilmore 1989; Xilouris et al. 1997). The application of these techniques to the region near the galactic planes was found to be problematic, mainly because of three reasons; (1) the presence of large amounts of dust, (2) the presence of additional components, e.g. in the form of spiral structure or a young thin disc, and (3) the limited accuracy, of about 1–2 degrees, with which the inclination can be determined.

Photometric studies of large nearby spirals which include a radiative transfer description, have pointed out that dust is more confined to the plane than the disc light; $`h_{z,\mathrm{dust}}0.5h_{z,\mathrm{stars}}`$ (Kylafis & Bahcall 1987; Xilouris et al. 1997, 1999). Considering the relatively small velocity dispersions of young stars (e.g. Wielen 1977), the same is true for young stellar populations. Hence, if we assume that the range in $`z/h_\mathrm{z}`$ affected by dust extinction and young populations is similar among galactic discs, it becomes possible to mask out that region in a systematic way. Indeed, the vertical $`IK`$ colour profiles of the galaxies in our sample, which trace the amount of dust extinction, show that in all cases the region $`|z|>1.5h_\mathrm{z}`$ is only marginally affected by dust extinction (de Grijs et al. 1997). Therefore, to solely study the old stellar population, we choose to mask out the region corresponding to $`|z|1.5h_\mathrm{z}`$. Since the differences among the various functionalities proposed for the vertical light distribution are small outside this mask (van der Kruit 1988, his fig. 1), the choice for the vertical distribution becomes unimportant. Still, we will further assess the presence of dust extinction outside the mask by (1) testing the 2D method on ‘dusty’ artificial images (Sect. 3.3) and (2) repeating the 2D fits in the B band to make a comparison with the I-band results (Sect. 4.1).

The 2D decomposition uses the non-linear least-squares technique developed by Marquardt (1963). The fit of Eqs. (24) is done in the linear regime, to the masked version of the sky-subtracted image (the data are not binned as in many previous methods, e.g. KS1–3, Shaw & Gilmore 1989). To ensure that the fitting algorithm minimizes the relative difference between our model and the data, each pixel is weighted inversely proportional to the value of the model surface brightness distribution of the disc at that pixel (cf. de Jong 1996a). As a consequence of this weighting, regions of the lowest signal-to-noise carry a large weight in the fit. To prevent this side effect, the data outside the region defined by the $`1\sigma `$ contour of the model image are excluded from the fit (the $`1\sigma `$ level is the standard deviation of the background). By using this criterion on the model instead of the data, there is no discrimination in favour of positive noise peaks. Since the adopted weighting scheme and the fitted region depend on the initial estimate of the model, an iterative application of the algorithm is required (see Sect. 3.3).

The light distribution of galactic discs may have a gradual or sharp radial steepening in the outer parts, indicating a possible truncation (e.g. KS1–3). If present, such a region will be excluded from the 2D fit by using an outer fitting boundary. This boundary will be assigned by eye, at a projected radius within the radius at which the radial light profile starts to deviate from that expected for an exponential disc. The error introduced by the subjectiveness of this procedure is estimated by performing several 2D fits in which the fitting boundary is varied by 10–20 percent.

### 3.3 Tests on artificial images

In order to assess the reliability of the method and the effectiveness of the mask used at $`|z|1.5h_\mathrm{z}`$, a series of tests was performed on artificial images. First, we created fifty artificial images using Eqs. (2) and (3). For the input parameters, random values were drawn from a uniform distribution covering the following ranges:

| *Disc* | *Bulge* |
| --- | --- |
| 16 $`<`$ $`\mu _0`$ $`<`$ 21 mag arcsec-2 | 17 $`<`$ $`\mu _\mathrm{e}`$ $`<`$ 22 mag arcsec-2 |
| 10 $`<`$ $`h_\mathrm{R}`$ $`<`$ 40 arcsec | 0.2 $`<`$ $`r_\mathrm{e}`$ $`<`$ 10 arcsec |
| 2 $`<`$ $`h_\mathrm{z}`$ $`<`$ 10 arcsec | 0.5 $`<`$ $`q`$ $`<`$ 1.0 |

where $`\mu _0`$ and $`\mu _e`$ are used to denote $`\mathrm{\Sigma }_0`$ and $`\mathrm{\Sigma }_e`$ in logarithmic units. To ensure the absence of artificial galaxies resembling bulge-dominated systems, the bulge-to-disc luminosity ratio was required to be less than two. Then, in order for the images to resemble the I-band observations, a sky background corresponding to $`\mu =20`$ mag arcsec-2 was added. The use of a different value for the background in each image is unnecessary, because sufficient contrast with the background is already provided by the variation in $`\mu _0`$ and $`\mu _e`$. The effect of seeing was mimicked by convolving the images with a circular Gaussian of 1.5 arcsec FWHM, which is typical for the observations. Finally, Gaussian read-out noise and Poisson noise were added and the background was subtracted.

For each artificial galaxy, initial estimates of the free parameters were generated by randomly offsetting the input parameters up to $`\pm 25`$ percent. The region corresponding to $`|z|1.5h_\mathrm{z}`$, was masked out. Then, the 2D fits were performed, while discarding data outside the $`1\sigma `$ contour of the model and redefining the mask and the pixel weights in each iteration. Generally, convergence was achieved within three iterations. In many of these 2D fits the bulge parameters did not converge. This is mainly due the use of the mask, causing the bulge contribution in the fitted region to be too small (less than 10 percent). However, for the same reason, the 2D fits could be safely repeated by fitting a disc-only model (cf. Fig. 2).

Figure 2 displays the test results for the disc parameters. It shows that the scale parameters are always recovered to within 10 percent, the standard deviation being about 2 percent for each parameter. As expected, the accuracy of the recovery decreases with decreasing central surface brightness. For example, the three outliers of Fig. 2b are all discs of low central surface brightness, $`\mu _0>20`$ mag arcsec-2. The amplitude of the deviations shown in Fig. 2 is mainly caused by the application of the mask; if no mask is applied the errors are less than 1 percent. Note that the application of the mask does not introduce any systematic error.

We repeated the tests after incorporating a strong dust lane. For the distribution of the dust, we chose a double exponential with $`h_{\mathrm{R},\mathrm{dust}}=h_{\mathrm{R},\mathrm{stars}}`$ and $`h_{\mathrm{z},\mathrm{dust}}=0.5h_{\mathrm{z},\mathrm{stars}}`$ neglecting scattering (see Kylafis & Bahcall 1987, their eq. (14)). Of course, in reality, dust shows clumpy and filamentary structures on small scales. However, the global nature of the fitting routine and the edge-on view largely average these out (Kuchinski et al. 1998), so that the exponential distribution is a good global approximation. Our neglect of scattering means that the effect of dust is slightly overestimated, since scattering would reduce the extinction. To resemble a strong dustlane, we chose a face-on central optical depth $`\tau _I(0)=1.0`$ (cf. Xilouris et al. 1999).

The effect of the dust on the scale parameters turned out to be moderate; the $`h_\mathrm{R}`$ and $`h_\mathrm{z}`$ were found to be on average 4 percent and 9 percent larger. The use of a smaller mask than that adopted here (or no mask at all) would mean that the effect becomes substantial, causing errors on the scale parameters larger than 10 percent.

The effect of a moderate stellar warp was investigated by simply introducing a vertical offset $`z_{\mathrm{warp}}(R,\varphi )`$ to the axisymmetric disc model (Eqn. 1) for $`R>R_{\mathrm{warp}}`$:

$`L(R,z,\varphi )`$ $`=`$ $`L(R)e^{|zz_{\mathrm{warp}}(R,\varphi )|/h_\mathrm{z}},`$
$`z_{\mathrm{warp}}(R,\varphi )`$ $`=`$ $`w[(RR_{\mathrm{warp}})/h_\mathrm{R}]\mathrm{cos}(\varphi )h_\mathrm{z},`$ (5)
where $`w`$ is the warping rate in units of scaleheight per scalelength, $`R_{\mathrm{warp}}`$ designates the onset of the warp, and the line of nodes of the warp is assumed to be straight and along the line of sight (such that the effect of the warp is the largest). For any reasonable warping, the effect on the derived scale parameters turned out to be small. For example, for a stellar warp with $`R_{\mathrm{warp}}`$ = 3 $`h_\mathrm{R}`$ and $`w=0.5`$, $`h_\mathrm{R}`$ and $`h_\mathrm{z}`$ are respectively under-and over-estimated by only 2 percent. This is mainly due to the fact that the inner parts of the edge-on disc are mostly unaffected by the warping.

We also studied the errors introduced by the fixed parameters. This was done by offsetting these parameters by their typical uncertainty and repeating the 2D fits. The following fixed parameters were investigated: sky subtraction – 1 percent error, seeing – 0.1′′error, central position – 1′′error, position angle – 1 degree error, and inclination – 3 degrees error. For example, for a disc with $`h_\mathrm{R}/h_\mathrm{z}=7`$ and an inclination of 87 degrees, our 2D method overestimates the scaleheight by 4 percent. The scalelength was found to be most sensitive to an error in the sky background, resulting in an error $``$ 4 percent. We estimate the total error due all systematic effects at $``$ 6 percent.

Further tests confirmed that the choice for the model used for the vertical light distribution is rather arbitrary when the region corresponding to $`|z|1.5h_\mathrm{z}`$ is discarded (as stated in Sect. 3.1). For example, when instead a vertical $`\mathrm{sech}(z/h_\mathrm{z})`$ distribution is used to create the artificial images, the scale parameters found by our 2D method differ from the input values by about 2 percent.

The algorithm does not take into account the effect of seeing on the disc component. This omission was tested by creating a number of artificial images, each with a different PSF, and applying the 2D fitting method. For the unlikely case that the seeing FWHM is comparable to the scaleheight of the disc, the error made by applying the method is only 3 percent. Since the scaleheight is much larger than the seeing FWHM for all galaxies in our sample, one can safely ignore this effect.

Finally, we investigated the effect of a truncation by applying the method to truncated artificial galaxies. These were constructed by including, in Eq. (1), $`L(R,z)=0.0`$ for $`R>R_{\mathrm{max}}`$ and numerically integrating along the line of sight. The difference in edge-on surface brightness between such a truncated disc and the infinite exponential disc increases with projected radius, causing the 2D method to underestimate the exponential scalelength. For galaxies with $`R_{\mathrm{max}}/h_\mathrm{R}4`$ this difference is small. For example, in a galaxy for which $`R_{\mathrm{max}}/h_\mathrm{R}`$ = 4, our method underestimates the scalelength by $``$ 6 percent. This error increases slightly for galaxies with smaller $`R_{\mathrm{max}}/h_\mathrm{R}`$, reaching $``$ 9 percent for $`R_{\mathrm{max}}/h_\mathrm{R}`$ = 3.

In summary, we conclude that the uncertainty in the scale parameters obtained with the 2D method is caused by several effects. In order of decreasing importance, these are: (1) the presence of a truncation, the error depending on the value of $`R_{\mathrm{max}}/h_\mathrm{R}`$, (2) residual dust extinction (error $`9`$ percent), (3) inaccurate sky subtraction (error $`4`$ percent), and (4) the application of the mask to $`|z|1.5h_\mathrm{z}`$ (error $`2`$ percent, see Fig. 2). For example, for discs with an $`R_{\mathrm{max}}/h_\mathrm{R}`$ of 4 (Sect. 5.2), the total measurement error in the scalelength is about 10 percent.

## 4 Results

### 4.1 The two dimensional fits

First, to prepare the images for the 2D fits, Hii regions, foreground stars and background objects were masked out using a combination of the sextractor algorithm (Bertin & Arnouts 1996) and manual editing in the galaxy’s immediate vicinity. The centre and major axis position angle of the galaxies were determined by fitting ellipses to the I-band images using the galphot package in iraf (Franx, Illingworth, & Heckman 1989) (see Table 12, appendix). Then, several profiles were extracted in the I band, both parallel and perpendicular to the major axis. The latter profiles were inspected for the effects of dust extinction. In most of these profiles, the turnover due to dust extinction lies below $`z=1.0h_\mathrm{z}`$ (see appendix, and de Grijs et al. 1997). However, in 16 galaxies the extinction clearly extends further out on one side of the galactic plane, e.g. in ESO 435-G14. This indicates that these galaxies are not exactly edge-on, and in these cases the affected side was excluded from the 2D fit.

Fits were made to each I-band image, taking the values determined in the 1D analysis (de Grijs 1998) as initial estimates. In 10 cases, the mask needed to be decreased from $`|z|1.5h_\mathrm{z}`$ to $`|z|1.0h_\mathrm{z}`$ in order to achieve convergence (Table 12, appendix). A bulge-disc decomposition was attempted in the 23 systems which clearly show an additional central component. As the best-fitting bulge model, either exponential (Eq. 3) or r1/4 (Eq. 4), the bulge-disc fit with the smallest reduced $`\chi ^2`$ was retained. However, in 14 galaxies the bulge parameters did not converge. This was not unexpected, since in many galaxies only a small portion of the bulge light is emitted in the region outside the mask. Still, two of the 14 systems do show an extended central component; ESO 321-G10 and ESO 416-G25. In these galaxies this component has a boxy or peanut-like shape, which probably caused the fit to diverge. To determine the disc parameters in these 14 cases, we assigned an inner radial fitting boundary, as in the ‘marking the disc’ method (Freeman 1970), and repeated the 2D fits with a disc component only. Tables 2 and 4 show the resulting parameters for the entire sample. The errors in the parameters are not the formal errors of the least-squares fit; these were in general less than a few percent and do not take into account the subjectiveness of the assignment of the radial fitting boundaries. Instead, the errors were estimated by comparing results from several similar fits in which the inner and outer radial fitting boundaries were varied by 10–20 percent. For each galaxy, the 2D fits are further illustrated in the appendix.

The 2D fits were repeated in the B band, using exactly the same fitting regions to further investigate the effect of dust extinction. For this, we calculated the scalelength ratios, $`h_\mathrm{R}^\mathrm{B}/h_\mathrm{R}^\mathrm{I}`$. A scalelength ratio is an indirect measure of the radial colour gradient of the disc. Since dust is known to be concentrated to the galactic centres and dust extinction increases towards bluer wavelengths, the presence of dust will result in a larger $`h_\mathrm{R}^\mathrm{B}/h_\mathrm{R}^\mathrm{I}`$. Consequently, when the observed scalelength ratio is compared to the ratio expected solely from radial stellar population changes, the effect of dust extinction can be inferred (Peletier et al. 1994; de Grijs 1998). The value of the scalelength ratio due to population changes is, however, somewhat uncertain. Peletier et al. (1994) argue that the scalelength ratio between the B and the I band due to population changes is in the range 1.05–1.1. From modelling of 6 nearby edge-on galaxies, which includes a radiative transfer description, Xilouris et al. (1997, 1999) find an average ratio of 1.25 $`\pm `$ 0.13 ($`1\sigma `$). Recently, from a study of 21 Sb galaxies of various inclination, Cunow (1999) derives a ratio of 1.13. For our sample, the average $`h_\mathrm{R}^\mathrm{B}/h_\mathrm{R}^\mathrm{I}`$ is $`1.17\pm 0.18(1\sigma )`$. This is clearly within the range of ratios attributed to radial population changes, and smaller than the ratios observed in less-inclined galaxies (Peletier et al. 1994). We conclude that the dust extinction in the fitted regions is small. This is in accordance with the results of Sect. 3.3.

In Fig. 4, the disc scale parameters are compared to those obtained from the 1D method (de Grijs 1998). For the scalelength, the differences are significant. The scalelengths determined with the 2D method are on average 13 percent larger. These differences arise mainly due to the combination of three effects. First, the fitted regions are different in both methods. According to Knapen & van der Kruit (1991), a difference of about 10 percent is not unusual in this respect. For example, the 1D method adopts a fixed radial outer fitting boundary at 4 scalelengths. However, many galaxies show a steepening of the radial surface brightness in this region, often caused by a truncation (Sect. 4.2). This effect causes the scalelength to be underestimated. Secondly, the bulge contribution is substantial in early-type spirals, even in the outer parts of the disc, also causing an underestimation of the scalelength. Thirdly, the 1D method adopts an exponential as a fitting model, instead of its projection (Eq. (2)). This effect leads to an overestimate of the scalelength by about 20 percent (Fig. 6).

For the vertical scaleheight, the differences between the two methods are mostly within 25 percent (Fig. 4b), except for 7 outliers for which the 2D method finds a scaleheight which is smaller by about 50 percent. These differences are caused by the differences in the fitting methods; the 2D method uses a two-dimensional fitting region and a bulge-disc decomposition, whereas in the 1D method the scaleheight is determined from profiles taken parallel to the major axis at various projected radii.

Our sample has 5 galaxies in common with the sample of Pohlen et al. (2000b). The model used by these authors is similar to that used here, but the method is different at two important points; (1) – no mask is applied to the region near the galactic plane and (2) – whereas with our method the truncation is largely avoided, their method is based upon the assumption that all discs show an infinitely sharp truncation. Table 6 lists both, the I-band scale parameters and truncation radii (Sect. 4.2), and the same parameters derived by Pohlen et al. (2000b) in the r band. The scale parameters agree well, except for the scalelengths of two galaxies, ESO 269-G15 and ESO 564-G27. For these, the determination of Pohlen et al. (2000b) is larger by about 30 percent, much larger than the difference expected between the I and the r-band from radial population changes alone. We note that for ESO 269-G15, our 2D fit was performed to the eastern side because it is not perfectly edge-on (see appendix). In ESO 564-G27, additional extended emission is present at large distances from the plane in the inner parts, indicating that the use of different fitting regions may be the cause of the discrepancy.

### 4.2 Fitting the truncation

The truncation, or cut-off, of a stellar disc can be inferred from a sharp drop in the radial light profile (e.g. KS1–3). The truncation radius is then obtained by extrapolating the radial light profile to zero flux. The error made in such an extrapolation is largely due to Poisson statistics, imperfect sky subtraction and foreground stars (KS1–3). While truncations are most easily observed in edge-on systems because of their higher surface brightness, there one views the result of a complicated line of sight integration across the truncated disc, which includes the effects of dust extinction and non-axisymmetric components such as spiral arms. Additionally, the truncation may not be axisymmetric, such that in an edge-on spiral the truncation becomes apparent at a different projected radius on either side (de Grijs et al. 2001). These considerations, as well as the variety of shapes of the radial light profiles seen in truncated edge-on spirals (e.g. Näslund & Jörsäter 1997; de Grijs et al. 2001), illustrate the difficulties faced in the construction of a model for the shape of the truncation from studies of edge-on spirals. For these reasons we did not include a truncation in the model adopted for the 2D fits (Sect. 3).

By applying an edge-fitting method similar to the one introduced by de Grijs et al. (2001), it is possible to objectively determine the truncation radius. The adopted method is as follows. A radial profile is constructed by first averaging vertically over the range $`1.0h_\mathrm{z}z2.0h_\mathrm{z}`$ on the least obscured side of the galactic plane. The choice for this range minimizes the effects of dust and young populations, while retaining a good signal-to-noise ratio. To further increase the signal-to-noise in the outer parts the profile is also averaged radially, by using a semi-logarithmic scheme (de Grijs et al. 2001). We made sure that this radial averaging does not smooth out the truncation by testing the edge-fitting method on artificial images (these were constructed as in Sect. 3.3 and included the infinitely sharp truncation). Then, the part of the profile containing flux from the galaxy is defined by clipping the vertically averaged profile at the 1$`\sigma `$ level ($`\sigma `$ being the standard deviation of the background in the vertically averaged profile), and using it as a mask (i.e. a method of conditional transfer). In this way, the outermost projected radius included in the profile is determined directly by the data quality. Following this, a series of trial fits of a straight line is performed in linear surface brightness at the edge of the profile (note that de Grijs et al. (2001) used an exponential function). In the first trial fit the four outermost points are included. Then, in each of the subsequent trial fits the inner fitting boundary is decreased by three points. Finally, only the fit with the smallest reduced $`\chi ^2`$ is retained and extrapolated to zero flux. The corresponding projected radius is then the adopted truncation radius, provided it satisfies the criterion outlined in the following.

An edge-fitting method was already applied, although in a somewhat different form, to four sample galaxies in B, V, and I band (de Grijs et al. 2001). There, as in all previous studies, the presence of a truncation is deduced from a qualitative comparison between the observed radial profile and the exponential disc model. In our sample, which is expected to contain many less obvious cases, a quantitative criterion is needed to determine its presence. To search for such a criterion, we performed tests of the edge-fitting method on artificial images. Not surprisingly, it was found that $`R_{\mathrm{max}}`$ is more easily recovered when the difference between the observed surface brightness and that expected for an untruncated exponential disc (Eq. (2)) is large. In particular, when the difference between the observed and expected surface brightness at the outermost radius of the profile is larger than about 2$`\sigma `$, then $`R_{\mathrm{max}}`$ could always be recovered ($`\sigma `$ is the standard deviation of the background in the original image). When the difference at the outermost radius is smaller than 2$`\sigma `$, which is caused by either a low central surface brightness and/or large $`R_{\mathrm{max}}/h_\mathrm{R}`$, then the recovery often fails. In what follows, a profile is said to be truncated if it satisfies this 2$`\sigma `$ criterion. The edge-fitting method is further illustrated in Fig. 8.

Of course, when measuring surface brightnesses in the outer parts of galactic discs, the accuracy of the background subtraction becomes very important. The main factors which limit this accuracy are large scale flatfielding errors and the presence of background sources and foreground stars; our small-scale flatfielding errors are within 0.5 percent (de Grijs 1997) and the CCDs were not affected by fringing. The background emission in the I-band images was originally estimated by fitting a plane to regions far away from the galaxy (de Grijs 1997). To check the influence of residual light due to objects within that region, we repeated these fits after heavily growing of our object masks. The results show that the earlier fits slightly overestimated the background, on average by 0.3 times the standard deviation in the background. Although this oversubtraction is small, much smaller than the amplitude necessary to produce artificial truncations (see below), it would affect the deduced values for $`R_{\mathrm{max}}`$ and we have therefore decided to apply the edge-fitting method to the new background-subtracted images. To estimate the amplitude of remaining variations due to large-scale flatfielding errors, we inspected the distributions of pixel values in 10x10 sized boxes at a large number of positions in the sky-subtracted background images. The standard deviation of the medians of these distributions was adopted as the error in the background determination.

The results of the application of the edge-fitting method to our sample are shown in Table 8 (and the appendix). To calculate the error on $`R_{\mathrm{max}}`$ two additional images were constructed for each galaxy by adding/subtracting the error in the background determination. The fits were repeated and the error on $`R_{\mathrm{max}}`$ was calculated by quadratically adding half the difference between the truncation radii found from the two modified images to the formal error of the original fit. For the galaxies for which no truncation could be determined, Table 10 shows a lower limit for $`R_{\mathrm{max}}`$.

To summarize, a truncation is found in 28 of the 68 profiles (counting each side of each galaxy separately). At the outermost point in these profiles, the differences between the exponential disc model and the data are larger than $`3\sigma `$, clearly satisfying our criterion. These differences are also much larger than the error in the background determination, i.e. the truncations are not produced by over-subtracting the sky background (Table 8 and the appendix). Of the other 40 profiles, about half had to be discarded because of the presence of foreground stars, background objects or a limiting field of view. In the remaining profiles, the difference between the data and the untruncated exponential disc at the outermost point is less than 2$`\sigma `$. The 28 truncations are found in 20 galaxies, 8 of which show a truncation on both sides. In all of these 8 cases the truncation appears symmetric; any asymmetries are less than 15 percent. Of the 12 galaxies in which a single truncation is seen, there are only 2 for which the profile at the other side of the galaxy could be studied as well; ESO 435-G14 and ESO 460-G31. Looking directly at the images of these, we note that their discs are not perfectly edge-on and show a signature of strong spiral arms which may be causing the single-sided truncation (cf. de Grijs et al. 2001). Galaxies which are not exactly edge-on and show signs of spiral structure are marked with an asterisk in Table 8. As an alternative possibility, the single truncations may be caused by a non-axisymmetric or lopsided stellar disc in which the truncation occurs at a smaller projected radius (= higher surface brightness) on one side.

To assess the effect of dust extinction, the fits were repeated using the region below one scaleheight. The resulting $`R_{\mathrm{max}}`$ are consistent with the values of Table 8, indicating that dust extinction near the edge of the disc is not very important. Observations of dust tracers such as CO emission and mm/submm continuum in the outer parts of spiral galaxies confirm this picture (Combes & Becquaert 1997; Neininger et al. 1996).

Other effects which may cause an apparent truncation are warping or flaring of the stellar disc. Although our sample contains no spirals which are strongly warped (Sect. 2), artificial images of warped discs (Sect. 3.3) show that a moderate warping can still introduce a break in the radial profile similar to the signature of a truncation. More importantly, away from the major axis the radial profiles show a strong asymmetry; a break at a small galactocentric radius on one side and an upturn followed by a break at a large galactocentric radius on the other. This causes a large difference between the apparent truncation radii of both sides. While similar upturns in a few of our sample galaxies may be a sign of warping (or the effect of spiral arms), the galaxies having a double-sided truncation show no evidence for such a strong asymmetry. To address disc flaring we also investigated artificial images of flared discs viewed edge-on. These which were constructed as in Sect. 3.3, but included a scaleheight increasing linearly with radius. It turned out that a very strong and localized thickening in the outer parts is required to produce a feature resembling a radial truncation (a rate larger than 1 $`h_\mathrm{z}`$ per $`h_\mathrm{R}`$). For our sample galaxies, and late-type spirals in general, such thickening is not observed (de Grijs & Peletier 1997; Morrison et al. 1994; Fry et al. 1999). We conclude that warping and flaring are very unlikely to have caused the truncations.

Again, we can compare our results to those obtained by Pohlen et al. (2000b), who estimated the truncation radius by visually comparing radial profiles to a truncated model (Table 6). For these 5 galaxies the results for the truncation radius are consistent. Thus, the truncation is seen in different data-sets using different methods, indicating that the truncations are real. This is further emphasized by the tight relation observed between $`h_\mathrm{R}`$ and $`R_{\mathrm{max}}`$ (Sect. 5.2).

Finally, we would like to emphasize that the fraction of galaxies in which a truncation is found, about 60 percent, is a lower limit. Clearly, the detection of a truncation depends on the edge-on central surface brightness, on the $`R_{\mathrm{max}}/h_\mathrm{R}`$ of the disc, on the limiting magnitude reached by the observation and, to a lesser extent, on the shape of the truncation. Since these elements are different in each image, the selection effect is different in each case.

## 5 Discussion

### 5.1 The disc flattening

We further investigate the disc scale parameters obtained in Sect. 4.1. In Figs. 10a & d the distributions of the I-band scalelengths and scaleheights are plotted. The corrected distributions were obtained by applying a volume correction to the diameter-limited distributions (van der Kruit 1987; Davies 1990; de Jong 1996b) and normalizing to the number of galaxies in the sample. The diameter selection limit implies that at small physical sizes the sampled volumes are too small to contain any edge-on spirals (from our sample the density of spirals with inclination $``$ 87  is only about 2 x 10-3 Mpc-3). This selection effect causes the low number of small spirals seen in the scalelength distribution. Small scalelength spirals are known to dominate the population of spirals (van der Kruit 1987; de Jong 1996b). Recognizing the fact that these small spirals are underrepresented, the shape of the volume-corrected distribution is consistent with the results of van der Kruit (1987) and de Jong (1996b). Interestingly, the distribution of scaleheight shows a similar shape, including a lack of small scaleheight galaxies (Fig. 10d). This lack is also due to the selection against intrinsically small galaxies because scaleheight and scalelength are observationally related (see below). Still, the smaller scaleheight discs are the most numerous; the volume-corrected distribution of scaleheights implies that about 90 percent of the spirals has $`h_\mathrm{z}<0.6`$ kpc.

For 33 spirals with known maximum rotational velocities (Table 2), both scale parameters are plotted as a function of maximum rotational velocity in Fig. 10b & e. This shows clearly that the disc scalelength increases with maximum rotational velocity. A Spearman rank correlation test results in a coefficient of 0.75, corresponding to a confidence level greater than 99 percent. Since the maximum rotational velocity is an indicator of total mass, this merely shows that more massive spirals tend to have larger stellar discs. More interestingly, the disc scaleheight shows a similar behaviour (Spearman correlation coefficient 0.74). Hence, more massive spiral galaxies have both larger and thicker discs. The most extreme example of this is ESO 509-G19. This huge spiral has an extremely large maximum rotational velocity of 392 km s-1 and a disc scalelength and scaleheight of 17 kpc and 1.7 kpc respectively. A similar correlation is seen between the scale parameters and the total H i mass (Fig. 10c & f), with the notable exception of ESO 437-G62 (= NGC 3390). This outlier has a remarkably small amount of H i and a thick disc, and is classified as an S0–a in the ESO-LV catalogue. For spirals however, more massive H i layers tend to be associated with larger and thicker stellar discs. We observe no clear trend with morphological type.

Observations of the stellar kinematics in 12 spirals by Bottema (1993) have revealed a related trend; more massive spirals tend to have larger stellar velocity dispersions. A simple straight line fit, adopted from his fig. 6, gives:

$$\sigma _\mathrm{z}|_{\mathrm{R}=0}=(0.29\pm 0.10)v_{\mathrm{max}},$$
(6)
where $`\sigma _\mathrm{z}|_{\mathrm{R}=0}`$ is the central vertical velocity dispersion of the disc stars. This vertical velocity dispersion is directly related to the scaleheight of the disc mass distribution through the vertical Jeans equation. When the vertical velocity dispersion is integrated in the vertical direction (see van der Kruit 1988), one has:

$$\sigma _\mathrm{z}|_{\mathrm{R}=0}=\sqrt{CG\mathrm{\Sigma }_0h_\mathrm{z}},$$
(7)
where $`\mathrm{\Sigma }_0`$ is the central disc surface density and $`C`$ is a constant which depends on the exact choice of the vertical density distribution; for an exponential distribution $`C`$ = $`3\pi /2`$, for a sech-dependence $`C`$ = 1.7051 $`\pi `$ and for the locally isothermal distribution $`C`$ = 1. If we assume that the scaleheight of the disc mass distribution is proportional to that of the old stellar population, Eq. (7) predicts that an increase in the vertical velocity dispersion is associated with an increase in the scaleheight of the light of the old population. Furthermore, by eliminating $`\sigma _z|_{R=0}`$ between Eq. (6) and (7), the scaleheight of the disc light is then also expected to be related to the maximum rotational velocity; $`h_zv_{\mathrm{max}}^2/\sqrt{\mathrm{\Sigma }_0}`$. As an example, we illustrate this in Fig. 10e (dashed line) for vertically exponential discs with face-on central surface brightness $`\mu _B=21.65`$ mag arcsec-2 (Freeman 1970) and constant mass-to-light ratio $`M/L_\mathrm{B}`$ = 2 $`(\mathrm{M}/\mathrm{L})_{}`$. While the disc mass-to-light ratios in and among spiral galaxies are still poorly known, at least the shape of the predicted relation for a constant mass-to-light ratio is roughly similar to that of the observations.

In Fig. 12 we investigate the ratio of the I-band scale parameters, which is a direct measure of the intrinsic flattening of the disc light. The volume corrected distribution of the disc flattening (Fig. 12a) is remarkably narrow; about 70 percent of spirals seem to have a flattening between 6 and 8. The average flattening of the diameter-limited distribution is $`h_\mathrm{R}/h_\mathrm{z}=8.5\pm 2.9`$ (1$`\sigma `$), while that of the volume corrected distribution is $`h_\mathrm{R}/h_\mathrm{z}=7.3\pm 2.2`$ (1$`\sigma `$). Figure 12b reveals a clear correlation between $`h_\mathrm{R}`$ and $`h_\mathrm{z}`$ (Spearman correlation coefficient 0.84). Considering the errors, the spread in the $`h_\mathrm{R}`$$`h_\mathrm{z}`$ plane is real, suggesting that the physical relation between $`h_\mathrm{R}`$ and $`h_\mathrm{z}`$ cannot be described by a simple linear functionality expected for discs of constant flattening. This is even more clear in Fig. 12c & d. The range and maximum of the disc flattening seem to increase with both maximum rotational velocity and total H i mass. A similar behaviour is seen with absolute magnitude $`M_\mathrm{I}^0`$ (not shown). The disc flattening shows a similar, but weaker, trend with Hubble type (de Grijs 1998; Schwarzkopf & Dettmar 2000), in accordance with the results on the disc axial ratio obtained by Guthrie (1992) and Kudrya et al. (1994). On the other hand, the minimum disc flattening seems to be independent of all of these parameters. At present, there is no clear physical explanation for these trends.

The disc flattening can be used to estimate the contribution of the disc to the rotation curve. For a self-gravitating exponential thin disc, there is a relation between maximum rotational velocity ($`v_{\mathrm{max}}^{\mathrm{disc}}`$) and disc scalelength (Freeman 1970):

$$v_{\mathrm{max}}^{\mathrm{disc}}=0.88\sqrt{\pi G\mathrm{\Sigma }_0h_\mathrm{R}}.$$
(8)
When the central surface density is substituted from Eq. (7), a simple relation between the maximum rotation, vertical velocity dispersion and flattening of the disc is obtained (Bottema 1993):

$$v_{\mathrm{max}}^{\mathrm{disc}}=(0.69\pm 0.03)\sigma _\mathrm{z}|_{\mathrm{R}=0}\sqrt{\frac{h_\mathrm{R}}{h_\mathrm{z}}},$$
(9)
where the constant accounts for density laws ranging from an exponential to a sech-dependence, and $`h_\mathrm{R}/h_\mathrm{z}`$ is the flattening of the disc mass. This equation states that a flatter disc with the same vertical stellar velocity dispersion or ‘temperature’ will be more massive. By eliminating $`\sigma _\mathrm{z}`$ between this equation and the empirical relation Eq. (6), one can estimate the contribution of the disc to the rotation curve:

$$v_{\mathrm{max}}^{\mathrm{disc}}/v_{\mathrm{max}}=(0.21\pm 0.08)\sqrt{\frac{h_\mathrm{R}}{h_\mathrm{z}}}.$$
(10)
By taking $`h_\mathrm{R}/h_\mathrm{z}`$ = 10 for the average flattening of the disc mass and using a locally isothermal vertical distribution, Bottema (1993) argued that the disc supplies on average 63 percent of the observed maximum rotation, i.e. $`v_{\mathrm{max}}^{\mathrm{disc}}/v_{\mathrm{max}}=0.63`$. Of course, this argument can be used for our sample, where we measured the flattening of the disc light in the I band. Assuming that the flattening of the disc mass equals that of the disc light from the old stellar population, and using the volume corrected average of the flattening of $`7.3\pm 2.2`$, we find $`v_{\mathrm{max}}^{\mathrm{disc}}/v_{\mathrm{max}}=0.57\pm 0.22`$. This confirms the result of Bottema (1993) and suggests that in a typical spiral galaxy, the actual disc mass-to-light ratio is about one third of that expected in a maximum disc situation (van Albada et al. 1985; Freeman 1992). However, the result is not entirely unambiguous. Given the error in the determination, the value of the disc contribution to the rotation could be about 80 percent. This would mean the disc can still be close to maximum, since, to avoid hollow halo cores and allow for a bulge contribution, a maximum disc has $`v_{\mathrm{max}}^{\mathrm{disc}}/v_{\mathrm{max}}0.9`$. The relatively large error in our determination stems from the scatter in the empirical relation Eq. (6), a relation which is still ill-constrained due to the small number of spirals for which stellar velocity dispersions are known.

### 5.2 The disc truncation

The truncation radius is an important global parameter, for which we need to find an explanation in theories of disc galaxy formation and/or evolution. In Sect. 4.2, we showed that the light distribution of at least 20 out of the 34 galaxies in our sample is radially truncated. These galaxies display a remarkably tight relation between the truncation radius and disc scalelength (Fig. 14a) (Spearman correlation coefficient 0.95). By taking the average $`R_{\mathrm{max}}`$ of both sides in case of a double-sided truncation, we obtain an average $`R_{\mathrm{max}}/h_\mathrm{R}=3.6\pm 0.6`$ (1$`\sigma `$) in the I band. Previous studies have also revealed the relation between the truncation radius and disc scalelength, but yielded different values for the ratio $`R_{\mathrm{max}}/h_\mathrm{R}`$. For example, for 7 nearby spirals KS3 find a mean value $`R_{\mathrm{max}}/h_\mathrm{R}=4.2\pm 0.5`$, Pohlen, Dettmar, & Lütticke (2000a) find $`2.9\pm 0.7`$ from 31 galaxies, whereas Schwarzkopf & Dettmar (2000) obtain $`3.7\pm 0.9`$ from 61 galaxies (although the latter two samples have several galaxies in common). These differences may be partly due to the use of different passbands; the use of bluer passbands will result in larger scalelengths, and hence leads to smaller $`R_{\mathrm{max}}/h_\mathrm{R}`$.

We further investigate the distribution of $`R_{\mathrm{max}}/h_\mathrm{R}`$ in Fig. 14b by plotting $`R_{\mathrm{max}}/h_\mathrm{R}`$ as a function of scalelength. For those galaxies for which no truncation was found, a lower limit is shown (arrows). The plot seems to reveal an increase in $`R_{\mathrm{max}}/h_\mathrm{R}`$ towards small scalelengths; the average ratio for the spirals with $`h_\mathrm{R}<4`$ kpc is 4.4, more than one standard deviation higher than that of the entire sample. This increase of $`R_{\mathrm{max}}/h_\mathrm{R}`$ may be related to its decrease found at very large scalelengths by Pohlen et al. (2000a). However, the reality of this feature is not clear considering the modest sample size and the selection effect against galaxies of small physical sizes. In addition, the ability to find a truncation depends, among other things, on the edge-on central surface brightness (Sect. 4.2). We are, for example, still missing low surface brightness galaxies, which may have entirely different $`R_{\mathrm{max}}/h_\mathrm{R}`$.

While our view of the distribution of $`R_{\mathrm{max}}/h_\mathrm{R}`$ is certainly obscured by these complex selection effects, additional information can be obtained from face-on samples. We have estimated a minimum $`R_{\mathrm{max}}`$ for the spirals in the face-on sample of de Jong & van der Kruit (1994). By taking the lowest contour in their R-band contour plots, we find an average $`\left(R_{\mathrm{max}}/h_\mathrm{R}\right)_{\mathrm{min}}=4.0\pm 1.1`$ (1$`\sigma `$). When combined with the results from our truncation analysis, this suggests that the ratio of truncation radius to disc scalelength in small scalelength spirals is at least four. This is important, because small spirals are the most numerous in the local Universe (van der Kruit 1987; de Jong 1996b, Sect. 5.1). Note that this also implies that the truncation is more difficult to detect in small scalelength spirals.

We can compare the results on the truncation radii to the predictions of the collapse model of disc galaxy formation (White & Rees 1978; Fall & Efstathiou 1980; Dalcanton et al. 1997). The collapse theory is compelling because its predictions are in general agreement with the basic structure and rotation curves of disc galaxies, over the full range of central surface brightness (Dalcanton et al. 1997; Mo, Mao, & White 1998). If it is further assumed that proto galaxies are sharp-edged then the collapse theory also predicts the outermost radius of the baryonic disc, which corresponds to the material with the highest specific angular momentum. To quantify that prediction and make a comparison with the observed truncation radii of stellar discs, we calculated model surface density profiles over a range in total mass, $`M_{\mathrm{tot}}`$ = 1010–1013 M$`_{}`$, and spin parameter, $`\lambda `$ = 0.01–0.28, using the method of Dalcanton et al. (1997). In the calculations, we used a fixed baryonic mass fraction $`\mathrm{F}=0.10`$, and a Hernquist halo at the onset of baryon-decoupling (as in Dalcanton et al. 1997). From the calculated surface density profiles, the scalelengths were derived using a method similar to the ‘marking-the-disc’ method (Freeman 1970) and the outermost radii were obtained by taking the radius at which the density drops to zero. The result is plotted in Fig. 14b as lines of constant mass. Both the outermost radius and the scalelength of the baryonic proto-disc increase with the mass and angular momentum of the proto galaxy, such that their ratio remains roughly constant at 3–4. Taking a different or non-constant baryon mass fraction does not significantly change this result. The predicted ratio is somewhat smaller than the ratio of 4.5, predicted by van der Kruit (1987) based on a comparison of the angular momentum distribution of an exponential disc to that of a uniformly rotating, uniform sphere with $`\lambda =0.07`$.

The predictions of the collapse model of disc galaxy formation roughly coincide with the observed truncation radii of stellar discs, although they do not show an increase of $`R_{\mathrm{max}}/h_\mathrm{R}`$ toward small scalelengths. It is clear from Fig. 14b that some galaxies fall above or below the predicted range, which may be attributed perhaps to environmental effects such as accretion and mergers. At any rate, in general the ratio $`R_{\mathrm{max}}/h_\mathrm{R}`$ predicted by the collapse model is similar to that seen in the stellar luminosity density. This may indicate that the present radial structure of the stellar disc still largely traces that of the baryonic proto-disc. If true, this result implies that H i present beyond the stellar disc (Broeils & Rhee 1997) is of a more recent origin (van der Kruit 1987). The result does not seem to be in accordance with the radial colour gradients observed in nearby face-on spirals. These indicate a stellar population change, with both the average stellar age and metallicity decreasing towards larger radii (de Jong 1996c). However, whether this shows that spiral galaxies have grown since their formation is largely an open question, for it is unclear if the radial colour gradient applies to the entire stellar disc or is caused by a young population in the plane of the galaxy.

It should be noted that our observations compare equally well to the value predicted by the theory of stochastic self-propagating star formation, for which, assuming a flat rotation curve, $`R_{\mathrm{max}}/h_\mathrm{R}`$ = 4 (Seiden et al. 1984). The other possible scenario, in which the truncation of the stellar disc is caused by a threshold on star formation due to its dynamical stability (Fall & Efstathiou 1980; Kennicutt 1989), does not make a direct prediction for $`R_{\mathrm{max}}`$. However, in a disc galaxy formation model similar to that of Dalcanton et al. (1997), van den Bosch (2001) has incorporated a star formation recipe which includes a threshold surface density based on Toomre’s (1964) criterion for local stability. For models of small bulge-to-disc ratio and with maximum rotational velocities exceeding 100 km s-1, comparable to the galaxies considered here, the presence of this threshold produces truncations in the stellar discs with $`R_{\mathrm{max}}/h_\mathrm{R}`$ = 2–4 (his fig. 7). This is in good agreement with the observations and also corresponds to the outermost radius of the proto-disc according to the collapse model investigated above. In particular, this latter coincidence seems to be caused by the fact that for massive late-type spirals the model of van den Bosch (2001) predicts that star formation occurs essentially throughout the entire gas disc, such that the present day truncation radius of the stellar disc is only marginally smaller than the outermost radius of the gaseous disc. We tentatively conclude that although the origin of the truncation of the stellar discs is still unclear, and may well be caused by a threshold on star formation, the truncation radius mainly reflects the material of the highest specific angular momentum which has collapsed during formation, at least for the intermediate to late-type galaxies considered here. To test whether the stellar truncation radius corresponds to the predicted threshold radius on star formation would require a detailed study of H i and H2 surface density profiles, as first suggested by Kennicutt (1989).

## 6 Conclusions

We have analysed I-band photometry of 34 edge-on spiral galaxies, taken from the sample of de Grijs (1998), to study the flattening and truncation of stellar discs. The exponential scale parameters of the disc light were derived using a 2D least-squares fitting method. This procedure includes a bulge-disc decomposition and adopts a mask in the region $`|z|1.5h_\mathrm{z}`$ to systematically exclude the region affected by dust extinction and young populations. It is better suited to extract the disc scalelength than the earlier 1D analysis (de Grijs 1998). We used an edge-fitting method, similar to the one introduced by de Grijs et al. (2001), to study the disc truncation. This method uses a simple straight line as a fitting function and a quantitative criterion to establish the presence of a truncation.

We find a clear increase in the scaleheight of the stellar disc as a function of maximum rotational velocity and total H i mass. This is in general accordance with observations of the stellar kinematics in spiral galaxies (Bottema 1993); larger discs are clearly more dynamically evolved and thicker. Both, the maximum and the range of the flattening of stellar discs seem to increase with maximum rotation and total H i mass. We used the average volume corrected flattening of the disc light, $`h_\mathrm{R}/h_\mathrm{z}=7.3`$, to estimate the disc contribution to the rotation curve, with the assumption that the flattening of the disc mass equals that of the disc light. This resulted in a disc contribution of 57$`\pm `$22 percent, which still – barely – allows for a maximum disc situation.

At least 20 galaxies in our sample (60 percent) are truncated, displaying a tight relation between scalelength and truncation radius. The average truncation radius corresponds to $`3.6\pm 0.6`$ I-band scalelengths. In addition, there seems to be a small increase in $`R_{\mathrm{max}}/h_\mathrm{R}`$ towards galaxies of small scalelength. Small scalelength spirals, which are the most numerous spirals in the local Universe, have an $`R_{\mathrm{max}}/h_\mathrm{R}`$ of at least four. A comparison of this ratio to that predicted by theories of disc galaxy formation seems to indicate that for the intermediate and late-type galaxies considered here the truncation radius mainly reflects the material of the highest specific angular momentum which has collapsed during formation. These scenarios are inconclusive with respect to the question of the origin of the truncation, and imply that a structure analysis of the truncation alone is insufficient to solve this issue.

It is clear that many questions are left unanswered. In an attempt to partly remedy this situation, we are currently studying the stellar and H i kinematics of the galaxies in our sample. This will, for instance, further improve the statistics on the empirical relation between maximum rotation and the stellar velocity dispersion, which could significantly advance our knowledge on the dark matter content of spiral galaxies. Still, studies of larger samples are needed to further constrain many of the relations found here, e.g. the trend of the disc flattening with disc mass. These future studies should preferably use deep, wide field, near infrared photometry. By targeting face-on as well as edge-on spirals, it becomes possible to more accurately study the disc truncation as well as its shape.

## Acknowledgments

We would like to thank Reynier Peletier and the anonymous referee for useful suggestions. MK would like to express his sincere thanks to Roelof Bottema for his sharp comments on an earlier version of this paper. RdeG acknowledges partial support from NASA grants NAG 5-3428 and NAG 5-6403, and hospitality at the Kapteyn Astronomical Institute on 2 visits. This work is based on observations obtained at the European Southern Observatory, La Silla, Chile. We have made use of the LEDA database (http://leda.univ-lyon1.fr).

## Appendix

We further illustrate the results of the 2D fits (Sect. 4.1) and the truncation analysis (Sect. 4.2). For each galaxy, information is presented in three panels.

The panel at the top left shows the surface brightness distribution in the I band, after editing and rotation (Sect. 4.1). The surface brightness at the lowest contour corresponds to the 3$`\sigma `$ noise level of the background. The 3$`\sigma `$ level is marked by the solid arrow on the magnitude scale of the lower left panel. The contour interval is 0.5 mag arcsec-2, except for ESO 033-G22, ESO 041-G09 and ESO 138-G14, where we chose an interval of 1.0 mag arcsec-2 for clarity. The orientation of each image is given in Table 12, columns (5) & (6). The inner and outer radial fitting boundaries and the mask (Sect. 4.1) are marked by the vertical and horizontal dotted lines respectively. The dashed contour marks the boundary outside of which the data were not included in the 2D fit (Sect. 3.2).

The bottom left panel shows, on the same scale as the top left image, two profiles taken parallel to the major axis. In each case the best-fitting model is shown for comparison (dashed lines). The bottom profile (solid line) is the average of two profiles. These were taken at a distance of 1.5 $`h_\mathrm{z}`$ on each side of the galactic plane and averaged over a range of 1/3 $`h_\mathrm{z}`$ to yield an acceptable signal-to-noise. When one side of the galaxy was fitted, only that side is used to construct this profile. The inner and outer radial fitting boundaries are indicated by the horizontal dotted lines. The top profile in this panel is the one used in the truncation analysis. It is shifted upwards by 1.5 magnitudes for clarity. If a truncation is present then the truncation fit (solid line) is also shown. The dashed–dotted lines indicate an envelope corresponding to $`\pm 3`$ times the background error. The open arrow at the magnitude scale, also shifted by 1.5 mag, indicates the standard deviation of the background in the vertically averaged profile used to clip the light profile (see Sect. 4.2 for further details).

The panel at the bottom right shows two profiles taken perpendicular to the major axis. The corresponding projected radii at which these were taken are mentioned in Table 12, columns (7) & (8). In each case, the model is shown for comparison (dashed lines). Again, each of the profiles is an average of two profiles, one taken on each side of the galaxy after averaging over a range of 1/3 $`h_\mathrm{R}`$. The boundary of the mask is indicated by the dotted lines.


format: text

text:
# Two-dimensional Site-Bond Percolation as an Example of Self-Averaging System.

## I Introduction

The influence of disorder on a phase transition is one of the important problems in the theory of phase transitions. In experimental measurements, the thermodynamic properties of one (or several) large-sized sample (with respect to the number of molecules) are usually studied. Therefore, it is important to know, whether a single large-sized sample with quenched realization of impurities can represent the properties of a model? The self-averaging properties of the system provide an answer to this question. If some quantity is self-averaging, the measurement for a single large sample gives a reasonable value for all samples of such size. If the quantity is not self-averaging, an increase in the system size does not make the measurement for a single sample representative.

Let us consider some statistical model with impurities on the $`d`$-dimensional lattice with a linear size $`L`$, in which we average some quantity $`X(\omega )`$ for a certain sample, where $`\omega `$ is the impurity realization from some impurity ensemble $`\mathrm{\Omega }`$. Here we assume that $`X(\omega )`$ is the exact thermal average for impurity realization $`\omega `$. Let us denote by $`[X]`$ the average over all impurity realizations of ensemble $`\mathrm{\Omega }`$: $`[X]=\underset{\omega \mathrm{\Omega }}{}P(\omega )X(\omega )`$, where $`P(\omega )`$ is the probability of impurity realization $`\omega `$. Then, variance $`V_X=\underset{\omega \mathrm{\Omega }}{}P(\omega )X^2(\omega )[X]^2`$. Let us define the relative variance of $`X`$: $`R_X=\frac{V(x)}{[X]^2}`$. This quantity characterizes the size dependent properties of the model. If $`R_X(L)L^x,x=d`$ the quantity $`X`$ is said to be strongly self-averaging. If $`x<0`$, $`X`$ said to be weak self-averaging, and if $`\underset{L\mathrm{}}{lim}R_X(L)\mathrm{const}`$, than $`X`$ is non-selfaveraging. We will now discuss self-averaging of a disordered system; therefore we examine the averaging over all impurity realizations, instead of self-averaging of a pure system, when we study the averaging over the whole statistical ensemble (for example, over all possible spin configuration for a spin models).

The Harris criterion states  that the weak randomness does not change the critical behavior of the $`d`$-dimensional second-order phase-transition model if the specific heat index $`\alpha <0`$, which corresponds to the pure system correlation length index $`\nu >2/d`$. As was first mentioned by Brout , far for critiality where the system size is much greater than the correlation length $`L\xi `$, all additive thermodynamical quantities are strongly self-averaging. The self-averaging properties of a disordered system near the critiality obeys Harris-Aharony (HA) criterion see  and . This criterion states that if randomness is irrelevant the system is governed by a pure fixed point and the relative variance $`R_XL^{\frac{\alpha _p}{\nu _p}}`$, where $`\alpha _p`$ and $`\nu _p`$ are the critical exponents of a pure system. If however, the system is governed by a random fixed point, then $`\underset{L\mathrm{}}{lim}R_X=\mathrm{const}`$. This criterion explains the numerical results of many works . Wiseman and Domany in  show that the Ashkin-Teller model with $`\alpha <0`$ is weakly self-averaging.

Two-dimensional percolation, which can be treated as a $`q`$-state Potts model with $`q1`$, is one of the most popular two-dimensional statistical models with the second-order phase transition. The two-dimensional percolation critical index $`\alpha =\frac{1}{2}<0`$, and, in accordance with the HA criterion we can expect that this model exhibits weak self-averaging properties. We study numerically the relative variances $`R_M`$ and $`R_\chi `$ of the probability $`M`$ of a site to belong the maximum cluster (the analog of magnetization) and the mean finite cluster size $`\chi `$ (the analog of the magnetic susceptibility). It was shown, that two-dimensional site-bound percolation, where bonds play the role of impurities and sites play the role of statistical ensemble over which the averaging is performed, exhibits self-averaging properties. The article is arranged as follows. In Section 2 we describe the site-bond percolation model and two types of impurities distribution; in Section 3 we discuss the self-averaging criterion and its phenomenological derivation; in Section 4 we present the numerical results and Section 5 is the conclusion.

## II Two-dimensional site-bond percolation model and canonical and grand-canonical impurity distributions

Site-bond percolation is a combination of purely site percolation and purely bond percolation. In purely site percolation, one investigates clusters of occupied sites. In purely bond percolation, one investigates clusters of sites, connected by occupied bonds. In site-bond percolation one, investigates clusters of occupied sites connected by occupied bonds. Each site is occupied with probability $`p_s`$, and each bond is occupied with probability $`p_b`$. In critical point $`(p_s^{},p_b^{})`$ the correlation length becomes infinite. Yanuka and Engelman  proposed approximate formula for a critical curve in $`(p_s,p_b)`$ plane

$$\frac{\mathrm{log}(p_s)}{\mathrm{log}(p_s^{})}+\frac{\mathrm{log}(p_b)}{\mathrm{log}(p_b^{})}=1$$

The two-dimensional site-bond percolation belongs the same universality class as the two-dimensional purely site percolation. Let us define by $`n_C`$ the mean number of clusters of size $`C`$ per lattice site. Then, the probability of a site to belong to the maximum cluster is

$$M=n_{C_{max}}C_{max}$$
(1)
and the mean finite cluster size

$$\chi =\frac{\underset{C,CC_{max}}{}n_CC^2}{\underset{C,CC_{max}}{}n_CC}$$
(2)
Near the critical point $`(p_s^{},p_b^{})`$ this quantities scales as follows

$$M(p_s,p_b=p_b^{})(p_sp_s^{})^\beta ,p_s>p_s^{},\beta =\frac{1}{18}$$
(3)
and

$$\chi (p_s,p_b=p_b^{})(p_sp_s^{})^\gamma ,\gamma =\frac{43}{18}$$
(4)
The correlation length $`\xi (p_s,p_b=p_b^{})(p_sp_s^{})^\nu `$ scaling index is $`\nu =\frac{4}{3}`$.

As mentioned above, each bond is occupied with probability $`p_b`$ and is empty with probability $`1p_b`$; therefore the total number of occupied bonds on the lattice fluctuates. This method of generating impurity configurations is known as Grand Canonical (GC) by analogy with Grand Canonical statistical ensemble with a fluctuating number of particles. However we can fix the number of occupied bonds $`N_b`$ and then distribute them randomly on the lattice. This method of generating impurity configurations is known as Canonical (C) by analogy with the Canonical statistical ensemble with fixing number of particles.

## III Self-Averaging criterion

Let us explain the HA criterion via simple phenomenological considerations . In this section, we use the $`T`$ temperature as a parameter of the model but keeping in mind, that in the case of percolation, the role of this parameter is played by the site concentration $`p_s`$.

We characterize every sample (impurity realization) $`\omega `$ with size $`L`$ by a pseudo-critical temperature $`T^{}(\omega ,L)`$. $`T^{}(\omega ,L)`$ fluctuates about its mean value and is averaged over all impurity configurations $`T^{}=[T^{}(\omega ,L)]`$. We introduce the reduced temperature for each sample

$$\dot{t}_\omega =\frac{TT^{}(\omega ,L)}{T^{}}$$
(5)
In the vicinity of the critical point, the quantity $`X`$ scales as

$$X_\omega (T,L)=L^\rho \stackrel{~}{Q}_\omega (\dot{t}_\omega L^{y_t})$$
(6)
Here, $`\rho `$ is the exponent characterizing the behavior of $`[X]`$ at $`T^{}`$ and the thermal scaling index $`y_t=\frac{1}{\nu }`$ is assumed to be universal for all samples. The form of function $`\stackrel{~}{Q}_\omega `$ is assumed to be sample dependent. We assume that (this relation is numerically checked in the next section)

$$\left(\delta T^{}\right)^2\left(\delta \dot{t}_\omega \right)^2L^d$$
(7)
So the variance of the argument of the function $`\stackrel{~}{Q}_\omega `$ scales as

$$\left(\delta \dot{t}_\omega L^{y_t}\right)^2L^{d+\frac{2}{\nu }}$$
(8)
We keep in mind that the scaling relation  $`\alpha =2\nu d`$. So, the variance $`V_x`$ of quantity $`X`$ at the critical point scales as

$$V_XL^{2\rho }\left(\delta \dot{t}_\omega L^{y_t}\right)^2L^{2\rho d+\frac{2}{\nu }}=L^{2\rho +\frac{2d\nu }{\nu }}=L^{2\rho +\frac{\alpha }{\nu }}$$
(9)
and the relative variance ($`[X]L^\rho `$) is

$$R_X=\frac{V_X}{[X]^2}=L^{2\rho +\frac{\alpha }{\nu }2\rho }=L^{\frac{\alpha }{\nu }}$$
(10)
We thus obtain the Harris-Aharony criterion: if the scaling index $`\alpha =const<0`$, then the model is self-averaging and the relative variance for quantity $`X`$ is $`R_XL^{\frac{\alpha }{\nu }}`$ from simple scaling relations.

## IV Numerical results

In our computation, the bond play the role of impurities; therefore, we fix the bond concentration $`p_b^{}=0.875`$ and consider the site concentration $`p_s`$ as parameter of our model.

1. we generate (by C and GC methods) the bond configuration $`\omega `$.
2. for this bond configuration we generate a set of site configurations (for site we use only GC method); For each site configuration we, calculate the probability of a site belonging to the maximum cluster and the mean finite cluster size;
3. We average these quantities over site realizations and obtain mean values $`M(p_s,p_b)`$ and $`\chi (p_s,p_b)`$. We split the set of site configurations into ten series to evaluate the numerical inaccuracy $`\mathrm{\Delta }M`$ and $`\mathrm{\Delta }\chi `$
4. we perform steps 1.–3. for another bond configuration.

The $`\chi _\omega (p_s)`$ dependence for three different bond realizations $`(\omega =1,2,3)`$ is shown in Fig.References. Let us calculate the pseudocritical site concentration $`p_s^{}(\omega )`$ for each bond realization $`\omega `$. We assume, that the mean finite cluster size $`\chi _\omega (p_s)`$ has maximum at the pseudocritical point $`p_s^{}(\omega )`$, and therefore, we approximate the data for $`\chi _\omega `$ near the maximum of $`\chi _\omega `$ by the parabola $`\chi _\omega (p_s)ab\left(p_sp_s^{}(\omega )\right)^2`$ and treat $`p_s^{}(\omega )`$ as the pseudocritical concentration for bond realization $`\omega `$. The locations of pseudocritical points for bond realizations $`\omega =1,2,3`$ is shown in Fig.References by vertical lines.

Let us find the critical point $`p_s^{}`$. We calculate the mean critical site concentration $`[p_s^{}(L)]`$ averaged over $`100`$ bond realizations created by GC and C methods as a function of lattice size $`L`$. These numerical data and results of approximation

$$\begin{array}{cc}\hfill \mathrm{C}\mathrm{method}:& [p_s^{}(L)]_C0.6519(5)0.34(10)L^{0.94(9)}\hfill \\ \hfill \mathrm{GC}\mathrm{method}:& [p_s^{}(L)]_{GC}0.6511(4)0.39(15)L^{1.02(11)}\hfill \end{array}$$
(11)
is shown in Fig.References. We take the concentration $`p_s^{}=0.6515`$ (averaged by C and GC methods) as critical. The variance of pseudocritical concentration behaves as $`(\delta p_s^{})^2L^2`$ (Fig.References), as we assumed in the previous section.

Let us investigate the self-averaging properties of $`M`$ and $`\chi `$ at the critical point $`(p_b^{}=0.875,p_s^{}=0.6515)`$.

Now we describe computational procedure that we use to calculate the relative variance. Here, we follow the method described . First, we note that, for each impurity realization $`\omega `$, instead of the exact value of quantity $`X_\omega `$ we get some value $`\overline{X}_\omega `$ averaged over site configurations with a numerical error

$$\left(\delta \overline{X}_\omega \right)^2=\frac{\sigma _{p_b,\omega }}{N_s/\tau _\omega }$$
(12)
Here, $`N_s`$ is the number of site configurations (the length of Monte-Carlo run) and $`\tau _\omega `$ is autocorrelation time. To calculate $`\left(\delta \overline{X}_\omega \right)^2`$ we split the MC sequence of site configurations of lenth $`N_s`$ into 10 subsequences and treat each subsequence as independent. We define by $`[\mathrm{}]`$ the averaging over impurity – bond configurations. Thus, the error of $`[\overline{X}_\omega ]`$ averaged over $`N_b`$ bond configurations is

$$\left(\delta [\overline{X}_\omega ]\right)^2=\frac{1}{N_b1}\underset{\omega \mathrm{\Omega }}{}\left([\overline{X}_\omega ^2][\overline{X}_\omega ]^2\right)$$
(13)
The total total error $`\left(\delta [\overline{X}_\omega ]\right)^2`$ has two contributing terms: from the sample to sample fluctuation of exact $`X_\omega `$ about $`[X_\omega ]`$ and from the fluctuations of $`\overline{X}_\omega `$, averaged over the finite number of spin configurations, about $`X_\omega `$ for each bond realization $`\omega `$.

$$\left(\delta [\overline{X}_\omega ]\right)^2=\frac{V_X}{N_b}+\left[\frac{\sigma _{p_b,\omega }}{N_bN_s/\tau _\omega }\right]$$
(14)
Here, $`N_b`$ is the number of bond configurations. We can calculate $`\left[\frac{\sigma _{p_b,\omega }}{N_bN_s/\tau _\omega }\right]=\frac{1}{N_b}\left[\left(\delta \overline{X}_\omega \right)^2\right]`$ by averaging the error $`\left(\delta \overline{X}_\omega \right)^2`$ over the bond configurations $`\omega \mathrm{\Omega }`$. So, we can express the relative variance via $`\left(\delta [\overline{X}_\omega ]\right)^2`$ and $`\left[(\delta \overline{X}_\omega )^2\right]`$

$$R_X=\frac{V_X}{[\overline{X}_\omega ]^2}=\frac{1}{[\overline{X}_\omega ]^2}\left(N_b\left(\delta [\overline{X}_\omega ]\right)^2\left[(\delta \overline{X}_\omega )^2\right]\right)$$
(15)
The numerical data for relative variance $`R_M`$ and $`R_\chi `$ computed by the C and GC methods in accordance (15), are shown in Fig.References and References respectively. We can see that points lies on the straight lines in log-log scale. Thus, we approximate the data by power function $`aL^b`$ on the interval $`L[48,128]`$. The results of approximation are plotted in Fig.References and References and placed below. As we might expect, for the GC method the relative variance is greater, than for the C-method because of the fluctuation of bonds number.

$$\mathrm{C}:R_M0.0041(3)L^{0.504(19)}$$

$$\mathrm{GC}:R_M0.0147(15)L^{0.50(2)}$$

$$\mathrm{C}:R_\chi 0.0243(2)L^{0.53(2)}$$

$$\mathrm{GC}:R_\chi 0.0259(5)L^{0.51(5)}$$

We see, that in excellent agreement with HA criterion, the relative variance of measured quantities obeys the power low dependence $`RL^{\frac{1}{2}}=L^{\frac{\alpha }{\nu }}`$. As a result, we can state that this model is self-averaging.

## V Conclusion

We have shown numerically that the two-dimensional site-bond percolation exhibit self-averaging properties at critical point. In our numerical experiments bonds play the role of the quenched disorder, and the sites play role of the statistical ensemble. We have found find that relative variance scales $`R_ML^{\frac{1}{2}}`$, $`R_\chi L^{\frac{1}{2}}`$, where $`\frac{1}{2}=\frac{\alpha }{\nu }`$. We assume that we can consider the sites as a quenched disorder and bonds, as a statistical ensemble, and the results will be the same.

We can also expect for the site percolation that, if we ”freeze” the state of some sites and average over other sites, we get the same weakly self-averaging behavior with respect this ”frozen” sites, which play the role of quenched disorder. Of cause, the same is valid for bond percolation.

I am grateful to the Joint SuperComputer Center RAN (www.jscc.ru) for providing computational resources.

Figure captions

Figure 1. The mean finite-cluster size $`\chi (p)`$ for three realization of impurity configurations. The critical concentrations $`p_s^{}`$ for each realization $`\omega `$, are shown by vertical lines.

Figure 2. The mean pseudocritical site concentration $`[p_s^{}]`$ as a function of lattice size $`L`$ for bond realizations, created by C and GC methods, the results of approximation by a power function of $`L`$, and the extrapolated critical concentration $`p_s^{}=0.6515`$

Figure 3. The variance $`\left(\delta p_s^{}\right)^2`$ of critical concentration as a function of lattice size $`L`$ for C (crosses) and GC (triangles) distributions, and the results of approximation by function $`aL^b`$.

Figure 4. The relative variance $`R_M`$ of probability $`M`$ of a site belonging to the maximum cluster as a function of lattice size $`L`$ for GC (crosses) and C (triangles) distributions and the results of approximation by function $`aL^b`$.

Figure 5. The relative variance $`R_\chi `$ of the mean finite-cluster size $`\chi `$ as a function of lattice size $`L`$ C (crosses) and GC (triangles) distributions and the results of approximation by function $`aL^b`$.