format: text

text:
# Structure of Disk Dominated Galaxies I. Bulge/Disk Parameters, Simulations, and Secular Evolution

## 1 Introduction

Stellar density distributions provide important constraints for bulge and disk formation models. Historically, astronomers have embraced the $`r^{1/4}`$ brightness “law” (de Vaucouleurs, 1948) and exponential brightness profile111The exponential nature of galaxy disk profiles emerges naturally in analytical models of disk formation (e.g. Lin & Pringle 1987; Dalcanton, Spergel, & Summers 1997; Ferguson & Clarke 2001). (de Vaucouleurs 1959a; Freeman 1970) to model the light distribution of the galaxy bulge and disk, respectively222It is important to remind ourselves from the onset that bulge-to-disk decompositions, and inward extrapolations of the disk into the central bulge and/or bar, may have no physical (or dynamical) basis. They provide a convenient description of the light distribution of a galaxy’s components that are otherwise dynamically coupled. The effective integrals of motions are likely similar for all the co-spatial components, though kinematically distinct bulges (counter-rotating nuclei) are known to exist.. Departures from the standard de Vaucouleurs profile in the inner light distribution of early- and late-type spirals have however been demonstrated in a number of early studies (de Vaucouleurs, 1959; van Houten, 1961; Burstein, 1979), including the Milky Way (Kent, Dame, & Fazio, 1991). Andredakis & Sanders (1994), de Jong (1996a), Courteau, de Jong, & Broeils (1996), and Carollo (1999) later used small samples of high-quality surface brightness (SB) profiles to establish the exponential profile as a better match to late-type disk bulges; thus SB profiles of most late-type spirals are best modeled by a double-exponential fit to the bulge and disk.

A broader analysis suggests a range of bulge shapes from early- to late-type spirals (Andredakis et al., 1995; de Jong, 1996a; Courteau et al., 1996; Graham, 2001). Most of these analyses rely on the modeling of a generalized surface density function such as that proposed by Sérsic (1968);

$$I(r)=I_0\mathrm{exp}\left\{\left(\frac{r}{r_0}\right)^{1/\mathrm{n}}\right\}$$
(1)
or, in magnitudes,

$$\mu (r)=\mu _0+2.5\mathrm{log}(e)\left\{\left(\frac{r}{r_0}\right)^{1/\mathrm{n}}\right\}.$$
(2)
where $`\mu _0`$ ($`I_0`$) is the central surface brightness (intensity), $`r_0`$ is a scaling radius, and the exponent $`1/n`$ is a shape parameter that describes the amount of curvature in the profile. For $`n=1`$ or 4 one recovers a pure exponential or the de Vaucouleurs $`r^{1/4}`$ profile respectively.

Collectively, the works above suggest that the bulge shape parameter $`n`$ correlates with absolute luminosity and half-light radius, such that bigger, brighter systems have larger values of $`n`$. This result was extended to brightest cluster galaxies by Graham et al. (1996). Courteau et al. (1996) also demonstrated a tight correlation between the bulge and disk exponential scale lengths, for all spiral types, with $`h_b/h_d=0.1\pm 0.05`$ (where $`h=r_0`$ and $`n=1`$ in Eq. 1). The exponential nature of late-type galaxy bulges and the correlation between bulge and disk scale lengths was interpreted by Courteau et al. (1996) as evidence for regulated bulge formation by redistribution of disk material to the galaxy center by a bar-like perturbation. We will return to this important constraint for secular evolution models in § 5.3.

This study focuses on the development of a reliable set of observables and constraints for structure formation models. An important goal is to measure the range of the Sérsic $`n`$ parameter for virialized disk systems. The analyses described above are reproduced and expanded upon with the largest multi-band survey of its kind to date and a clearer understanding of model limitations than previously attained. We aim to characterize and quantify the intrinsic structural properties of the bulge and disk and the extent of their variation with wavelength. These characterizations are made through reliable modeling of bulge and disk parameters from SB profile decompositions. Multi-wavelength information also provides insight about structural variations within and among galaxies due to dust and stellar population effects. While some of these issues have been addressed before, there remains a number of significant measurement uncertainties and technical limitations which we now investigate thoroughly.

This paper is organized as follows: a brief description of the database is given in § 2 and in § 3 we discuss our B/D decomposition algorithms (1D and 2D) and the simulations to test the reliability of our technique. For the readers interested mostly in final profile decompositions and results, a summary of the simulation results and guidelines is given in § 3.4. Actual B/D decompositions of galaxy SB profiles are presented in § 4, followed by a discussion and interpretation of the results in terms of secular evolution models in § 5. A discussion on the nature of Freeman Type-II profiles is also presented in § 5. We conclude with future directions in § 6. Two appendices present (A) a discussion of the functional form for the Sérsic coefficient $`b_n`$, and (B) decomposition results for our Type-I profiles.

## 2 The Data

Our structural analysis of galaxy luminosity profiles is based on the catalog of multi-band images of late-type spiral galaxies by Courteau, Holtzman, & MacArthur (2002; hereafter Paper II). It consists of over 1000 deep B, V, R, and H images of 322 nearby bright late-type spiral galaxies. The data were collected between 1992 and 1996 at Lowell Observatory and Kitt Peak National Observatory (KPNO). A full description of the sample selection, observations, and reductions is presented in Paper II. A summary is given below.

The galaxy sample was selected from the Uppsala General Catalogue (UGC, Nilson 1973) with the following criteria:

* Predominently late Hubble types
* Zwicky magnitude $`m_B`$ 15.5
* Blue Galactic extinction $`A_B=4\times E(BV)`$ 0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$5 (Burstein & Heiles 1984)
* Inclination bins covering face-on ($`i`$ 6°), intermediate (50$`{}_{}{}^{}<i<`$ 60), and edge-on ($`i`$ 78°) projections
* Blue major axis $`3\mathrm{}`$ .

This catalog is not complete in any sense of the (much abused) term. The diameter limit was constrained primarily by the field of view of the infrared cameras in use at KPNO (IRIM and COB) and Lowell Observatory (OSIRIS) in 1992–1996 and the requirement for blank areas in the field of view for sky subtraction. Additionally, peculiar and interacting galaxies (e.g. no visible tidal tails) were excluded to ensure that the sample consisted only of isolated disk dominated galaxies. Barred galaxies, as classified in the UGC, were not excluded per se but only a handful were observed. For the present analysis, we use a sub-sample of 121 galaxies with face-on and intermediate inclinations only, for a total of 523 images333Radial brightness profiles cannot be measured for fully edge-on galaxies. Our images of edge-on systems will be used in a forthcoming analysis of stellar and dust scale heights and truncation radii in spiral disks.. The distribution of Hubble types in our reduced sample is: 2 Sab, 26 Sb, 19 Sbc, 38 Sc, 25 Scd, 11 Sd.

All distances are corrected to the reference frame of the Local Standard of Rest (Courteau & van den Bergh, 1999), and we use H$`_\text{0}`$ = 70 km sec$`^{\text{-1}}`$ Mpc$`^{\text{-1}}`$. The survey effective depth is $`cz5500\mathrm{km}\mathrm{s}^1`$ or $`80`$ Mpc.

### 2.1 Observations and Basic Reductions

All optical BVR images were obtained from 1992 to 1994 at Lowell Observatory with a TI 800$`\times `$800 chip (scale = 0$`\stackrel{}{\mathrm{.}}`$5/pix) on the Perkins 72″ telescope. The infrared H-band images were acquired from 1993 to 1995 at KPNO with the 2-meter and 4-meter telescopes equipped with either a HgCdTe (IRIM) or an InSb (COB) 256$`\times `$256 array (1$`\stackrel{}{\mathrm{.}}`$09/pix and 0$`\stackrel{}{\mathrm{.}}`$5/pix respectively), and from 1995 to 1996 with the OSIRIS imager (1$`\stackrel{}{\mathrm{.}}`$49/pix) mounted on the Lowell 72″ telescope. The exposure times were: 300s in R, 400s in V, 1500s in B, and on-target integration of 1200s in H. Landolt (1992) standards covering a wide range of airmasses and colors were observed each night at Lowell Observatory, giving a photometric accuracy of $``$2% for the optical passbands. UKIRT standards (Guarnieri, Dixon, & Longmore 1991) observed each night yielded H-band photometric calibrations good to $``$3%. Stars and defects were edited from the images prior to further analysis.

The typical seeing full-width at half-maximum (FWHM) at Lowell and KPNO was 2$`\stackrel{}{\mathrm{.}}`$0 with typical standard deviations of $``$20% (optical) and $``$35% (IR) per image. These measurements were computed as the mean of the FWHMs of all non-saturated stars measured automatically on each image frame; typically 5 to 40 measurements per frame were used.

We measure mean sky levels (for a 5–6 day old moon) of $`B=21.9\pm 0.8,V=21.2\pm 0.5,R=20.6\pm 0.5`$, and $`H=14.1\pm 1.2`$ mag arcsec-2 . Typical systematic errors in the sky measurement, computed from 4 or 5 sky boxes suitably located between the galaxy and the edge of the frame, are $`0.51.0\%`$ in the optical and $`0.0050.01\%`$ in the IR.

Azimuthally-averaged SB profiles were extracted for all the galaxies using ellipse fitting with a fixed center. To ensure a homogeneous computation of structural parameters and color gradients, we use the isophotal maps from the R-band to determine the SB profiles in BVH. Even though dust effects can still play a role at 7000 Å, the R-band was adopted for our isophotal templates as it has the most stable sky and deepest profiles. We allowed a variable position angle and ellipticity at each isophote, but a comparison with SB profiles extracted using concentric isophotal fits demonstrated that our results do not depend on the fitting technique. Further information about profile extraction and CCD surface photometry can be found in Courteau (1996a) and Paper II. We trace SB profiles to $``$26 mag arcsec-2 in optical bands and $``$22 mag arcsec-2 at H-band. These levels correspond to a surface brightness error of $``$0.12 mag arcsec-2 .

### 2.2 Surface Brightness Corrections

The observed surface brightness of a galaxy can change when viewed at different inclination angles, depending on the distribution of a galaxy’s interstellar medium and its opacity. Surface brightnesses are also affected by Galactic foreground extinction and redshift dimming. We account for the latter effects but defer any treatment of internal extinction, which vary greatly from author to author, to Paper II. Our conclusions do not depend on the exact values of the central and effective surface brightnesses of galaxies.

We correct for Galactic foreground extinction using the reddening values, $`A_\lambda `$, of Schlegel et al. (1998) and assuming an RV = 3.1 extinction curve (e.g. Cardelli et al. 1989),

$$\mu _{c,Gal}^\lambda =\mu _{obs}^\lambda A_\lambda .$$
(3)
We correct surface brightnesses for the $`(1+z)^3`$ cosmological redshift dimming (per unit frequency interval) as

$$\mu _{c,z}^\lambda =\mu _{obs}^\lambda 7.5\mathrm{log}(1+z).$$
(4)
The final correction to the observed surface brightnesses is thus,

$$\mu _c^\lambda =\mu _{obs}^\lambda A_\lambda 7.5\mathrm{log}(1+z).$$
(5)
Examples of the types and extent of the SB profiles for typical late-type spiral galaxies in our sample are shown in Fig. 1. For Type-I disks (Freeman, 1970), the inner profile always lies above the surface brightness of the inward extrapolation of the outer disk, whereas Type-II systems have a portion of their brightness profiles lying below the inward disk extrapolation. We define a Transition case for luminosity profiles that change from Type-II at optical wavelengths to Type-I in the infrared. Many galaxies classified as Type-II show a weakening of the inner profile dip at longer wavelengths and, in this sense there is no clear distinction between the Type-II and Transition galaxies. Likely interpretations for Type-II profiles are discussed in § refsubsec:typeII.

## 3 Simulations of Bulge-to-Disk Decompositions

In order to measure galaxy structural parameters, we have developed two independent algorithms to decompose the galaxy 1D and 2D light distributions into bulge and disk components. These programs allow for a generalized Sérsic bulge, an exponential disk, and a central bar for 2D images. There are several issues involved with accurate decompositions, particularly with the measurement of bulge parameters, including; the sensitivity of final results to starting guesses, effects of statistical and systematic errors in sky brightness and seeing estimates, choice of fit baseline, etc. We explore these in great detail below using both 1D and 2D analyses to determine the robustness of our codes and the reliability of our final solutions. Because projected surface brightness profiles contain fewer data points than full 2D images, we can create 1D simulations faster than 2D models. Thus our most extensive tests rely on 1D simulations, which are shown to be fully consistent with 2D simulations when considering axisymmetric features.

### 3.1 1D and 2D Algorithms

Our brightness profile (1D) bulge-to-disk (B/D) decomposition algorithm was initially developed by Broeils & Courteau (1997) and subsequently improved by LM. This program reduces 1D projected galaxy luminosity profiles into bulge and disk components simultaneously using a non-linear Levenberg-Marquardt least-squares (NLLS; see § 15.5 in Press et al. (1992)) fit to the logarithmic intensities (i.e. magnitude units). Random SB errors are accounted for in the (data$``$model) minimization, whereas systematic errors such as uncertainties in the sky background and determination of the image mean PSF are accounted for separately in a series of experiments designed to calibrate their effects. Seeing effects in our model galaxies are accounted for by convolving the theoretical bulge-disk surface brightness profiles and images with a radially symmetric Gaussian Point Spread Function (PSF), of the form

$$I_s(r)=\sigma ^2e^{r^2/2\sigma ^2}_0^{\mathrm{}}I_{total}(x)I_0(xr/\sigma ^2)e^{x^2/2\sigma ^2}x𝑑x$$
(6)
where $`I_{total}(x)`$ is the intrinsic surface brightness profile, $`\sigma `$ is the dispersion of the Gaussian PSF and $`I_0`$ is the zero-order modified Bessel function of the first kind (see also Trujillo et al. (2001) for a study of the Moffat PSF.)

The 2D decomposition program is based on the same NLLS technique as above but uses the full 2D image in intensity units instead of a logarithmic radial surface brightness profile. While computationally more intensive than its 1D analogue, the 2D decomposition should yield more physically meaningful results since the azimuthal information is lost in 1D profiles. Byun & Freeman (1995), de Jong (1996a), and Simard et al. (2002) have discussed the merits of the 2D approach, such as greater ability to recover true parameters (based on simulations), and the potential to model non-axisymmetric features such as bars, rings, and spiral arms. The need for the implementation and testing of a robust 2D B/D decomposition package is thus obvious, but we find that 1D decompositions compare favorably for reliability and predictive power provided high S/N 1D radial profiles are used. Note that neither 1D nor 2D decompositions are impervious to dust extinction effects. Extinction effects are lessened at H-band, but can still be significant in disk bulges and spiral arms. A proper recovery of the true stellar density profile would require a full 3D radiation transfer treatment, and such an analysis is beyond the scope of this work.

### 3.2 Methodology

A fundamental aspect of profile decompositions is the choice of fitting functions. The disk light is modeled with the usual exponential function,

$$I_d(r)=I_0\mathrm{exp}\left\{\frac{r}{h}\right\}$$
(7)
or, in magnitudes,

$$\mu _d(r)=\mu _0+2.5\mathrm{log}(e)\left\{\frac{r}{h}\right\}$$
(8)
where $`\mu _02.5\mathrm{log}I_0`$ and $`h`$ are the disk central surface brightness (CSB) and scale length respectively, and $`r`$ is the galactocentric radius measured along the major axis. In the 2D decompositions, the computation of the radius at each pixel requires two additional parameters: the position angle (PA) of the disk major axis on the sky and the disk ellipticity, $`\epsilon =1b/a`$, where $`a`$ and $`b`$ are the major and minor axes of the disk respectively). To test for the shape of the bulge luminosity profiles we adopt the generalized formulation of Sérsic (Eqs. 1 & 2).

It has become customary to express the disk parameters in terms of scale length and CSB ($`h`$ and $`\mu _0`$), while the bulge parameters are expressed in terms of effective parameters ($`r_e`$ and $`\mu _e`$). We adopt this formalism, thus parameters with subscript $`e`$ refer to the bulge. Eq. 1 can be re-written as:

$$I_b(r)=I_e\mathrm{exp}\left\{b_n\left[\left(\frac{r}{r_e}\right)^{1/\mathrm{n}}1\right]\right\}$$
(9)
where the effective radius, $`r_e`$, encloses half the total extrapolated luminosity444For a pure exponential disk, $`r_e=1.678h`$.. $`I_e`$ is the intensity at this radius and $`b_n`$ is chosen to ensure that

$$_0^{\mathrm{}}I_b(r)\mathrm{\hspace{0.17em}2}\pi r𝑑r=2_0^{r_e}I_b(r)\mathrm{\hspace{0.17em}2}\pi r𝑑r.$$
(10)
In magnitudes Eq. 9 translates to

$$\mu _b(r)=\mu _e+2.5\mathrm{log}(e)b_n\left[\left(\frac{r}{r_e}\right)^{1/n}1\right]$$
(11)
where $`\mu _e`$ is the effective surface brightness.

It is trivial to convert from Eq. 2 to Eq. 11 by noting that

$`r_e`$ $`=`$ $`(b_n)^nr_0`$ (12)
$`\mu _e`$ $`=`$ $`\mu _0+2.5\mathrm{log}(e)b_n.`$ (13)
Eq. 10 implies that

$$\mathrm{\Gamma }(2n)=2\gamma (2n,b_n)$$
(14)
where $`\mathrm{\Gamma }(a)`$ is the gamma function and $`\gamma (a,x)`$ is the incomplete gamma function. Unfortunately, Eq. 14 cannot be solved analytically for $`b_n`$. Various numerical approximations have been been given in the literature (Caon et al. 1993; Graham & Prieto 1999; Ciotti & Bertin 1999; Khosroshahi, Wadadekar, & Kembhavi 2000; Möllenhoff & Heidt 2001). One often encounters the approximation $`b_n2n0.32`$, valid supposedly for all values of $`n`$ (sic). Khosroshahi et al. (2000) contend that this approximation is accurate to one part in $`10^5`$, with a range of validity on $`n`$ unspecified. However, because the gamma function diverges near the origin, most utilized approximations are inaccurate for values of the Sérsic exponent $`n1.`$ Differences between numerical solutions for $`b_n`$ (Eq. 14) and commonly adopted approximations can yield brightness differences greater than 0.1 mag arcsec-2 for $`n1`$. As we wish to test for bulges with Sérsic $`n`$ parameter smaller than 1, we have adopted the asymptotic expansion of Ciotti & Bertin (1999) to O$`(n^5)`$ for $`n>0.36`$. For $`n0.36`$ this solution diverges and instead we use a polynomial expression ($`4^{th}`$ order) accurate to one part in $`10^3`$. We compare different numerical solutions for $`b_n`$ (Fig. 23) and present our adopted functional form in Appendix A.

An illustration of profile shapes for different values of the Sérsic $`n`$ parameter is shown in Fig. 2. The top panel shows profiles with $`\mu _e=21`$ mag arcsec-2 and $`r_e`$ = 3$`\stackrel{}{\mathrm{.}}`$5 for values of $`n`$ in the range $`0.2<n<4`$. The bottom panel shows the same profiles but for a constant CSB of $`\mu _0=18`$ mag arcsec-2 . For $`n<1`$ the profiles are shallow at small radii ($`r_e`$) and fall off rapidly with increasing radius. Conversely, profiles with $`n>1`$ are steep at small radii ($`r_e`$), but level off as $`r`$ increases. Given the large differences in the profile shapes above and below $`n=1`$ (exponential case), one might expect different physical mechanisms (formation, transport, dynamics, interactions) to be at work for systems whose light profiles have very different $`n`$ values. Additionally, for the small bulges of late-type galaxies, poor seeing could conceivably smear the image such that an intrinsically $`n>1`$ bulge could be mistaken for an $`n<1`$ structure.

Of potential relevance to the study of galaxy structure is the relative light fraction contributed by the bulge and disk. This is expressed in terms of a bulge-to-disk luminosity ratio, $`B/D`$, derived by integrating the bulge and disk luminosity profiles to infinity. For a face-on Sérsic profile the total extrapolated luminosity is given by

$$L_b=_0^{\mathrm{}}I_b(r)2\pi r𝑑r=\frac{2\pi I_er_e^2e^{b_n}n\mathrm{\Gamma }(2n)}{b_n^{2n}}$$
(15)
and for a face-on exponential disk

$$L_d=_0^{\mathrm{}}I_d(r)2\pi r𝑑r=2\pi I_0h^2$$
(16)
giving a bulge-to-disk light ratio of

$$B/D=\frac{e^{b_n}n\mathrm{\Gamma }(2n)}{b_n^{2n}}\left(\frac{r_e}{h}\right)^2\left(\frac{I_e}{I_o}\right).$$
(17)
Eqs. 15 & 16 should be multiplied by the factor $`(b/a)`$ when considering projections on the plane of the sky. One may use Eq. 17 in a general sense, independent of projection, under the assumption that the bulge and disk density distributions have similar axes ratio (nearly true for late-type galaxies). A weakness of $`B/D`$ ratios for systematic comparisons of galaxy light profiles is its model dependence and the potential covariances between some of the model parameters. Consider Fig. 2 (top) for the relative light fractions contributed by profiles of different $`n`$ values, normalized to $`n=1`$. The integrated bulge light increases steadily as a function of $`n`$, for given values of $`r_e`$ and $`\mu _e`$. Thus, the adopted $`n`$ value in a bulge-to-disk decomposition has a strong influence on the computed $`B/D`$ ratio. Additionally, since larger $`n`$ profiles contribute light out to large $`r`$, the combination of a high $`n`$ and a low $`\mu _e`$ (bright $`r^{1/4}`$ bulge) could take away light from the outer disk and artificially boost the $`B/D`$ ratio. A discussion on non-parametric statistics, such as concentration indices (Kent, 1985; Courteau, 1996a; Graham, 2001) which alleviate model dependences, is presented in Paper II.

We model the total galaxy luminosity profile as a sum of bulge + disk components:

$$I_{tot}(r)=I_b(r)+I_d(r).$$
(18)
Profile smearing by atmospheric turbulence is accounted for in B/D decompositions by convolving Eq. 18 with a Gaussian PSF of the form of Eq. 6.

Similar B/D analyses have also considered additional terms for a Gaussian bar (de Jong, 1996a), a lens or ring (Prieto et al., 2001), spiral arms, and stellar disks with inner and/or outer truncations (Kormendy 1997; Baggett, Baggett, & Anderson 1998). We restrict our choice of fitting functions to a Sérsic bulge and a non-truncated exponential disk for a number of reasons. We find no prominent bars in our sample and most our disk profiles are fairly linear (in magnitude space). Azimuthal averaging for 1D profiles smoothes out spiral arm features (to a different extent depending on whether the position angle was fixed or allowed to vary in the profile extraction. Removal of spiral arm signatures from the light profiles or images would require more time and effort than is warranted by our analysis at this stage.) We do not consider a sharp inner disk truncation for a number of reasons: (i) unsharp masking techniques reveal spiral structure from the inner disk into the galaxy center (Courteau 1992, 1996b; Elmegreen, Elmegreen, & Eberwein 2001); (ii) using HST images of inner disks, Carollo (1999) also finds evidence for inner spiral structure and nuclear star clusters in the centers of early- to intermediate-type spiral galaxies; (iii) all components of the Galaxy have their peak surface brightnesses in the center (e.g. Wyse 1999). Thus, at least some evidence suggests that spiral disks reach in all the way to the center of late-type systems. A lowering of the disk central surface density may occur as stars get heated up into a bulge by the action of a bar-like instability. An exponential profile with a core may thus be a reasonable description of Type-II SB profiles. We do not consider this approach here, but point out that resolved B/D kinematics of nearby galaxies would provide a clear indication whether stellar populations have been strongly depleted and/or systematically scattered vertically into a bulge. Some of our galaxies show outer disk truncation (e.g. Fig. 13), but see § 4.3.2.

The best-fit parameters of the (data$``$model) comparison are those which minimize the reduced chi-square merit function, described in intensity units as

$$\chi _\nu ^2=\frac{1}{\mathrm{N}\mathrm{M}}\underset{i=1}{\overset{\mathrm{N}}{}}\left[\frac{I_{gal}(r_i)I_s(r_i;h,I_0,r_e,I_e,n)}{\sigma _i}\right]^2$$
(19)
where N is the number of data points used, M is the number free parameters (i.e. N - M = $`\nu `$ Degrees of Freedom), and $`\sigma _i`$ is the statistical intensity error at each pixel (2D) or surface brightness level (1D). From here on the $`\nu `$ subscript will be omitted and the $`\chi ^2`$ variable refers to a $`\chi ^2`$ per degree of freedom (unless otherwise specified).

The global $`\chi ^2`$ of intensities is clearly dominated by the contribution from the disk, virtually irrespective of the fitted bulge. This effect would be accentuated in galaxies with prominent features, such as spiral arms, rings, or lenses, which are not accounted for in our pure exponential disk models. Cases are found where B/D decompositions with significantly different bulge exponent $`n`$ values for a given profile have nearly the same global $`\chi ^2`$ value (see Figs. 8 & 9). Thus, in order to refine our parameter search for the best-fit bulge and disk model, we compute a separate, inner, $`\chi ^2`$ statistic out to twice the radius where the bulge and disk contribute equally to the total luminosity of the galaxy ($`r_{b=d}2r(I_b=I_d)`$). We label this statistic as $`\chi _{in}^2`$ (see Graham 2001 for a similar formulation). For cases where the bulges are so small that they never truly dominate the light profile (i.e. $`r_{b=d}`$ is undefined), we compute the $`\chi _{in}^2`$ out to the radius at which $`\nu =1`$.

### 3.3 Reliability of the Decomposition Results

This section describes extensive testing of our bulge-to-disk decomposition programs. Artificial SB profiles and images were created with a wide range of bulge profile shapes and exponential disks including realistic noise and seeing effects. Real galaxies are clearly more complicated than the sum of two idealized mathematical functions, but these tests provide a reasonable base for a global understanding of the reliability and limitations of B/D decomposition algorithms. The mock catalog of SB profiles and images will be used to address the following questions:

* How reliable and meaningful are the bulge-to-disk decompositions and fitted parameters?
* How crucial are initial estimates? Are model fits always converging to the lowest $`\chi ^2`$ minimum?
* How do seeing effects and sky subtraction errors affect the decompositions, and can they be properly accounted for?
* Are the small bulges in late-type disk galaxies sufficiently resolved to permit a reliable solution of the Sérsic $`n`$ parameter as a free parameter?

The literature abounds with investigations of profile fitting algorithms based on artificial data, such as Schombert & Bothun (1987; hereafter SB87) who performed double-blind experiments where one of the authors created mock luminosity profiles and the other independently fitted the data. The SB profiles combined a de Vaucouleurs bulge and an exponential disk. Photon noise, at a level matching typical blue CCD performances, and a systematic 0.5–3.0% error of the sky background were added to the profiles. SB87 found that the simultaneous fitting of disk and bulge using standard NLLS techniques could reproduce the input parameters to within 10–20% in cases where galaxy profiles can be decomposed perfectly as the sum of a bulge and disk (which fails for Type-II profiles.) SB87 claim that a sky estimate uncertainty of up to 3% does not affect their derived parameters, but we find that sky errors as small as 1% can have a significant effect on the shape of the outer disk profile and the derived bulge and disk parameters (see § 3.3.5 below). SB87 did not consider other fitting functions but recognized that bulges may not be adequately described by the de Vaucouleurs $`r^{1/4}`$ function. Andredakis & Sanders (1994) later examined the inadequacy of the $`r^{1/4}`$ functional form for the bulge (1D) profile, and first established the double-exponential nature of late-type spirals.

2D B/D decomposition techniques, which exploit the full galaxy image, were also developed and tested in similar fashion in the mid-nineties (Byun & Freeman, 1995; de Jong, 1996a). De Jong performed extensive tests with mock galaxies modeled as pure exponential bulges and disks, exploring the effects of errors in the measured observables including the seeing FWHM, sky background level, minor over major axis ratio, $`b/a`$, and position angle, PA. These observables were used as fixed input parameters to the fitting routine, and de Jong calibrated the effect of measurement error on the determined parameters by decomposing the artificial galaxies using erroneous values for each observable. He concluded that: errors in $`\mu _0`$ are predominantly caused by sky subtraction errors and can be as large as 0.1 mag arcsec-2 ; errors in $`h`$ can reach 10% and are dominated by sky background and ellipticity measurement errors; bulge parameter errors, of order 20%, are controlled by the B/D size and brightness ratios. Bright bulges are most affected by seeing errors, and fainter bulges can also be affected by sky background errors.

Our own investigation reaches similar conclusions and further extends de Jong’s simulations. We test for the robustness of the fitting procedure and accuracy of the derived parameters with various values of the fit initial estimates, seeing FWHM, sky value, and their errors, and – unlike de Jong (1996a) – we model the bulge with a generalized Sérsic profile. These simulations were initiated by Broeils & Courteau (1997) but are extended here in much greater detail, especially with respect to the determination of the bulge shape parameter $`n`$.

#### 3.3.1 Simulated Profiles and Images

Our tests use a large set of artificial SB profiles and images which span a wide range of the bulge, disk, and seeing parameters. The mathematical forms of the bulge and disk components are those discussed in § 3.1. Noise was added to the model profiles and images from a Gaussian distribution with deviation representative of the standard brightness errors of our luminosity profiles at a given surface brightness level (see Courteau 1996a, Fig. 9; Paper II) .

One hundred SB profiles and fourty images with realistic noise were created for each bulge, disk, and seeing combination. Most of the simulated profiles and images had the same disk parameters,

$`\mu _0=20`$ mag arcsec-2 ,

$`h=12`$″,

which are representative of a typical galaxy in our sample555Different values of $`h`$ were also tested but the results do not change significantly.. For exponential ($`n=1`$) bulges, the structural parameters were selected from

$`\mu _e`$ = 16, 17,…, 22 mag arcsec-2 ,

$`r_e`$ = 0.1, 0.2,…, 3$`\stackrel{}{\mathrm{.}}`$0

corresponding to $`B/D`$ ratios ranging from 0 to 5 and $`B/T`$ from 0 to 0.8 (see Eq. 17).

All profiles and images were convolved with a seeing disk of

$`\mathrm{FWHM}`$ = 1.0, 1.5,…, 3$`\stackrel{}{\mathrm{.}}`$0.

Our models span the full range of parameters typically found in late-type bulges (e.g. de Jong (1996b); Courteau (1996a); Broeils & Courteau (1997)) and the seeing values match the expected range at Lowell Observatory and KPNO (see § 2).

We also explored the following range of Sérsic $`n`$ values

$`n`$ = 0.2, 0.4,…, 4.0

for the set of combinations with $`r_e`$ = 0.8, 1.5, and 2$`\stackrel{}{\mathrm{.}}`$5, $`\mu _e`$ = 18, 20, and 22 mag arcsec-2 and for seeing FWHMs of 1.5, 2.0 and 2$`\stackrel{}{\mathrm{.}}`$5. Here the corresponding $`B/D`$ ratios range from about 0 to 1, and $`B/T`$ from 0 to 0.5. We also simulated the full range of $`r_e`$, FWHM, and $`\mu _e`$ = 18, 20, 22, and 24 mag arcsec-2 for $`n=0.2`$ and $`n=4.0`$ for a test regarding the initial estimates (see § 4.2).

A total of about 223,000 artificial SB profiles and 16,200 images were created and modeled. The mock profiles were all sampled at 0$`\stackrel{}{\mathrm{.}}`$5/pixel to match the data.

The 1D and 2D decompositions of the artificial profiles and images were deemed satisfactory if they met the following (rather liberal) criteria:

* solution found within 100 iterations
* $`0<r_{e(fit)}<500`$
* $`0<\mu _{e(fit)}<30`$ mag arcsec-2
* $`0.05<n_{(fit)}<10`$
* $`0<h_{(fit)}<100`$
* $`0<\mu _{0(fit)}<30`$ mag arcsec-2

For each set of parameters, we require that at least two-thirds of the 100(40) simulated profile (image) decompositions pass these criteria to be included in the analysis. The 2D tests were not developed as fully due to prohibitive computing times. Thus our tests rely more heavily on the 1D technique, but we have confirmed that the results from both techniques corroborate each other.

#### 3.3.2 Disk Initial Estimates

NLLS algorithms require initial estimates as input parameters, and we must verify whether our final solutions are sensitive to our initial guesses. Different initial estimates may yield different solutions with comparable $`\chi ^2`$ values especially if the topology of the $`\chi ^2`$ distribution is non-trivial or shallow (e.g. Schombert & Bothun 1987; de Jong 1996a). Our simulations confirm the robustness of our algorithms to a wide range of initial disk parameter estimates. The results are slightly more sensitive for large values of $`n`$, but in general the disk parameters were perfectly recovered independent of the initial guesses. One must still caution that if the bulge is fit with the wrong $`n`$, the fitted disk parameters will differ from their intrinsic values, even if the initial estimates were good. Fig. 3 shows the relative fit error for $`h`$, $`\mathrm{\Delta }h`$, where

$$\mathrm{\Delta }h\frac{h_{\mathrm{fit}}(\mathrm{mean})h_{\mathrm{model}}}{h_{\mathrm{model}}}$$
(20)
versus the fitted $`n`$ value (held as a fixed parameter in the decomposition) for bulges with $`n`$ = 0.2, 1.0, and 4.0. (i.e. $`\mathrm{\Delta }h=0`$ indicates a perfect recovery of the model parameter $`h`$).

In these decompositions the correct initial values were given for the bulge and disk parameters and the seeing FWHM was assumed to be known. Different seeing values are represented by three point types: circles, triangles, and squares for seeing values of 1.5, 2.0, and 2$`\stackrel{}{\mathrm{.}}`$5 respectively. We see that the fitted disk parameters can have errors as large as 10% when the bulge is modeled with an incorrect shape parameter. The reason for this is obvious if one considers the different shapes for the different values of $`n`$ shown above in Fig. 2, where the different shapes contribute differently to the outer profile. The bulge $`r_e`$ is even more sensitive to the fitted $`n`$ (relative errors in excess of 300% for the worst cases, the figure is not shown).

#### 3.3.3 Bulge Initial Estimates

We now investigate the importance of the bulge initial estimates. Given that the fit baseline for late-type bulges is much smaller than that of the disk and the apparent size of these bulges is comparable to the seeing disk, the bulge model is likely to be much more sensitive to the input parameters. The sensitivity to the initial estimates of bulge parameters $`r_e`$ and $`\mu _e`$ was tested using offsets of $`\pm `$ 25% and $`\pm `$ 50% from the model, as well as the correct model values as initial estimates. This range of offsets matches typical estimate excersions in our data modeling (see § 4.2). The Sérsic $`n`$ was presumed known and the tests were performed for three different values of $`n=`$0.2, 1, and 4.

The tests for the $`n=1`$ case demonstrated that bulge parameter initial estimates are not important in both the 1D and 2D decompositions, as long as $`r_e`$ 0.3$``$FWHM. Below this limit initial estimates are more important; errors on $`r_e`$ can exceed 50% and errors of up to $`\mathrm{\Delta }\mu _e\pm 0.3`$ mag arcsec-2 can occur. In the $`n=0.2`$ case the parameter recovery is largely independent of the seeing FWHM, as expected for profiles that are flat in the center. Given incorrect initial guesses, recovered parameters for profiles with $`r_e`$ 1$`\stackrel{}{\mathrm{.}}`$0 can still be trusted, but for smaller bulges the algorithm is trapped in a local minimum and the output value is nearly the same as the input, e.g. a $`\pm 25(50)`$% input error yields a $`\pm `$25(50)% output error. For $`n=4.0`$ profiles, the parameter recovery is strongly dependent on the seeing FWHM such that decomposition results for profiles with $`r_e`$ 0.7$``$FWHM cannot be trusted. Here again input and output errors are approximately equal. Moreover, even with correct initial estimates, the model parameters are not perfectly recovered for profiles with $`r_e`$ 1$`\stackrel{}{\mathrm{.}}`$0.

We have interpolated these results for different values of $`n`$ and define a parameter space for which our solutions are not affected by the choice of initial estimates:

$$r_e(0.3)^{1/n}\text{FWHM}$$

and

$$r_e\{\begin{array}{cc}0.75n+1.15\hfill & \text{ for }n1.0\hfill \\ \mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.2}n+0.2\hfill & \text{ for }n1.0\hfill \end{array}$$

The corresponding results for the 2D decomposition algorithm closely match those from the 1D tests.

#### 3.3.4 Seeing Effects

The effect of an uncertainty in the seeing FWHM measurement on the model parameters ($`n`$, $`r_e`$, $`\mu _e`$, $`h`$, and $`\mu _0`$) must be accounted for in our simulations. We fit each model using not only the nominal seeing value, but we also varied the seeing FWHM by typical seeing measurement errors (1 $`\sigma `$; $`1520`$% for optical, $`35`$% for infrared; see Paper II). Figs. 4 & 5 show the effect of an incorrect seeing estimate on the fitted $`r_e`$ for the 1D and 2D algorithms respectively. Plotted are the relative fit error on $`r_e`$, $`\mathrm{\Delta }r_e`$, where

$$\mathrm{\Delta }r_e\frac{r_{\mathrm{e},\mathrm{fit}}(\mathrm{mean})r_{\mathrm{e},\mathrm{model}}}{r_{\mathrm{e},\mathrm{model}}}.$$
(21)
The 1D and 2D simulations agree very well. Figs. 4 & 5 show that when the correct seeing is used as input, the bulge and disk parameters are recovered perfectly for the full range of bulge parameters and for all values of the seeing FWHM tested. However, even moderate seeing uncertainties can severely affect bulge parameters depending on the size of the bulge relative to the seeing FWHM. If the seeing width is under(over)-estimated $`r_e`$ is systematically over(under)-estimated, worsening for smaller and fainter bulges and larger seeing values. Similar trends are seen for $`\mu _e`$. Our tests show that the fit errors can be significantly larger if the seeing FWHM is over-estimated than if it is under-estimated. As a rough rule of thumb, for $`r_e\mathrm{FWHM}`$ and a seeing measurement uncertainty at the 35% (15%) level, the bulge $`r_e`$ can be trusted to within 10–25% (0–10%), and $`\mu _e`$ to within $`\pm `$ 0.1–0.4 (0–0.2) mag arcsec-2 , the lower end of the range applying to the brightest bulges and increasing towards the upper end for the fainter bulges. For $`r_e\mathrm{FWHM}+1`$ the errors improve to within 0–15% (0–10%) for $`r_e`$, and $`\pm `$ 0.0–0.2 (0.0–0.05) mag arcsec-2 for $`\mu _e`$.

There is no appreciable effect due to seeing on the disk parameters (less than 1%) except for the worst case of a FWHM of 3$`\stackrel{}{\mathrm{.}}`$0 and a 35% seeing over-estimate. In all other cases, the disk parameters are virtually unaffected by seeing, as the size of the disk is much larger than the seeing profile. However, it is of paramount importance to use accurate seeing estimates and realistic seeing errors in order to sample the true range of bulge parameters in B/D decompositions.

#### 3.3.5 Sky Uncertainty Effects

We now test for the effects of an improper sky subtraction on the decompositions. The tremendous sensitivity of B/D decomposition and scale length determinations to sky errors has been highlighted before (Courteau, 1992; de Jong, 1996a). Here we aim to provide a firm quantitative assessment of such errors. We re-model the same simulated profiles as in the previous section but using sky values that are $`\pm `$1% of the nominal sky level (typical error in the optical passbands), and using a typical optical sky brightness of 21 mag arcsec-2 666Note that the percent sky error in the H-band is much smaller than in the optical, but the sky at H is much brighter, so the effect should be comparable.. Since bulge brightnesses are typically greater than the sky level, at least at optical wavelengths, one might expect bulge parameters to be somewhat insensitive to sky subtraction errors. However, the outer disk is very sensitive to sky subtraction errors and a modified disk ultimately affects bulge structure due to their coupling. Quantitatively, our tests show that if the sky is over- or under-estimated by 1%, the error on the disk scale length, $`\mathrm{\Delta }h`$, will be of order 5–15% and the disk CSB, $`\mathrm{\Delta }\mu _0`$, will be $`\pm `$ 0.1–0.25 mag arcsec-2 . These dispersions hold for the full range of bulge brightnesses except the two faintest bulges which are one and two magnitudes fainter than the sky; these had errors in excess of 50%. The error ranges are controlled by the relative sizes of the bulge and disk such that the disk parameter errors increase slightly from smaller to larger bulges. This is simply because a larger bulge weakens the importance of the disk in the central parts, thus giving more weight to the sky-sensitive outer disk.

The errors on the bulge parameters are negligible for bulges with $`(\mu _{sky}\mu _e)>1`$ mag arcsec-2 for the entire range of $`r_e`$ and seeing FWHMs, but increase up to $`\mathrm{\Delta }r_e`$ 15% and $`|\mathrm{\Delta }\mu _e|0.1`$ mag arcsec-2 (increasing as the bulge gets smaller and as seeing conditions degrade) for bulges with $`(\mu _{sky}\mu _e)<1`$ mag arcsec-2 . In other words, if the bulge effective surface brightness is less than one magnitude greater than the sky brightness, the bulge parameters will be strongly affected by sky subtraction errors. This effect is often neglected in studies of bulge/disk structure.

The bulge and disk parameters are most affected for the case of an under-subtracted sky. This is largely due to a magnitude threshold of 26.5 mag arcsec-2 in our decomposition algorithm. The data are too noisy below this value (Paper II) and we exclude them from the fits. This threshold provides some protection against over-subtracted skies in the measurement of the disk scale length. Similar tests were performed with our 2D decomposition algorithm which confirm, once again, the results above.

#### 3.3.6 Sérsic n Tests

A number of recent studies have described the variation of bulge shapes as a function of Hubble type (Andredakis et al., 1995; Moriondo, Giovanardi, & Hunt, 1998; Khosroshahi et al., 2000; Möllenhoff & Heidt, 2001), going so far as suggesting that precise values of $`n`$ (i.e. $`\pm 0.1`$) could be determined (Graham, 2001). To our knowledge, no study to date has tested the reliability of the recovery of the Sérsic $`n`$ parameter. In order to test the sensitivity of the decomposition to the full range of bulge profile shapes we use mock luminosity profiles with values of the Sérsic $`n`$ parameter ranging from $`n=0.2`$ to $`n=4.0`$. The suite of profiles used all combinations of $`r_e`$ = 0.8, 1.5, and 2$`\stackrel{}{\mathrm{.}}`$5, $`\mu _e=18`$, 20, and 22 mag arcsec-2 and seeing FWHMs of 1.5, 2.0, and 2$`\stackrel{}{\mathrm{.}}`$5. The profile fits used initial estimates of $`n=0.4`$, 1, 2, and 4, and correct initial estimates for $`r_e`$, $`\mu _e`$, and the disk parameters. The seeing FWHM was fixed to the correct model value. The results for the $`n=1`$ initial estimate are presented in Fig. 6, where we plot the average relative fit error on $`n`$ (for the 100 profiles with the same simulated parameters) where

$$\mathrm{\Delta }n\frac{n_{\text{fit}}(\mathrm{mean})n_{\text{model}}}{n_{\text{model}}}$$
(22)
versus the model $`n`$ (the dashed line at $`\mathrm{\Delta }n=0`$ indicates a perfect recovery of the model $`n`$ parameter.) Each panel shows one particular combination of $`\mu _e`$ and $`r_e`$, and the panels are arranged such that the $`B/D`$ ratio, for a given value of $`n`$, decreases from top to bottom and right to left. The different seeing FWHM values are represented by three point types: circles, triangles, and squares for seeing values of 1$`\stackrel{}{\mathrm{.}}`$5, 2$`\stackrel{}{\mathrm{.}}`$0, and 2$`\stackrel{}{\mathrm{.}}`$5 respectively.

Fig. 6 reinforces that the bulges of even nearby late-type spirals are small and not sampled at high enough spatial resolution to yield a stable, robust solution for $`n`$ as a floating parameter. Given the correct value of $`n`$ as an initial estimate (along with the correct initial estimates for the other four parameters), the algorithm normally finds the correct value of the model $`n`$, but any departure from the model value even by a small amount, yields significantly different solutions for $`n`$, or the fit may simply fail (as indicated by the vertical lines in the figures). For most of the parameter combinations, an offset of $``$50% in the initial estimate of $`n`$ yields a $``$50% error on its determined value.

Similar tests using the 2D algorithm show a slightly more robust recovery of the model $`n`$ parameter based on incorrect initial estimates, but the recovery efficiency is still poor and results based on a floating initial estimate of $`n`$ are questionable. We are here faced with an under-determined optimization with too few independent data points for too many model parameters (at least three for the bulge). The strong covariances between $`n`$, $`\mu _e`$, and $`r_e`$ ($`\sigma _{n,\mu _e},\sigma _{n,r_e}`$) prevent a unique determination of $`n`$ with this NLLS code. We actually find the best-fit $`n`$ by grid search, holding $`n`$ as a fixed parameter, solving for a range of values, and using the $`\chi _{in}^2`$ as defined in § 3.1 to determine the best fit. Further simulations showed this technique to be fully reliable for the bulges considered here.

It is difficult to estimate the error on $`n`$. Either we use a grid search and $`n`$ is fixed, or $`n`$ is kept as a floating parameter and varies widely given wrong initial estimates. Based on a test with floating $`n`$ but correct initial estimates for bulge parameters, typical seeing and sky errors modify $`n`$ by no more than 20%.

Based on 2D B/D decompositions of simulated images and JHK images of 40 bright spirals, Möllenhoff & Heidt (2001) estimate that recovery errors for all the standard fit parameters ($`I_d,h,I_e,r_e,`$ and $`n`$) are less than 15%, comparable to our findings. They also tested for variable estimates of the sky level and seeing width, though no clear description of their technique is presented.

### 3.4 Summary of the Simulations

The tests performed in § 3.3, for idealized galaxies, allow us to define a set of guidelines for the reliability and limitations of our 1D/2D decompositions:

* Initial estimates for bulge and disk parameters are unimportant provided that

$$r_e(0.3)^{1/n}\text{FWHM}$$

and

$$r_e\{\begin{array}{cc}0.75n+1.15\hfill & \text{ for }n1.0\hfill \\ \mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.2}n+0.2\hfill & \text{ for }n1.0\hfill \end{array}$$
* Seeing errors must be accounted for in all bulge parameter studies. For $`r_e\mathrm{FWHM}`$ and a seeing measurement uncertainty at the 35% (15%) level, the bulge $`r_e`$ can be trusted to within 10–25% (0–10%), and $`\mu _e`$ to within $`\pm `$ 0.1–0.4 (0–0.2) mag arcsec-2 . For $`r_e\mathrm{FWHM}+1`$ the errors improve to within 0–15% (0–10%) for $`r_e`$, and $`\pm `$ 0.0–0.2 (0.0–0.05) mag arcsec-2 for $`\mu _e`$. There is no appreciable effect due to seeing on the disk parameters (less than 1%).
* Sky subtraction errors dominate disk parameter errors ($``$5–15%) and are non-negligible (up to 25%) for bulges whose effective surface brightnesses are less than one magnitude brighter than the sky brightness (i.e. for $`\{\mu _{sky}\mu _e\}1`$).
* The sampling of late-type nearby bulges may not be high enough to constrain the Sérsic $`n`$ exponent uniquely as a free parameter. Iterative model fitting schemes should be tested for this. Our approach uses a grid search.
* Typical seeing and sky errors modify $`n`$ by no more than 20%.
* The 2D decomposition technique does not provide a significant improvement over the 1D method for the recovery of axisymmetric structural parameters to warrant the extra computational effort.

Armed with these basic guidelines we can now turn our attention to real data decompositions using the 1D technique.

## 4 Bulge-to-Disk Decompositions

### 4.1 Outline

Our study of structural properties and the variation of galaxian parameters as a function of wavelength uses the multi-band (BVRH) data set of late-type spiral galaxies of Courteau, Holtzman, & MacArthur (§ 2; Paper II). Most galaxies have at least one set of BVRH images, and we use multiple observations for 54 galaxies to estimate systematic errors. Other B/D decomposition analyses have used larger samples (e.g. Baggett et al. (1998)) but lack the crucial multi-wavelength information.

We aim to develop a stable and versatile prescription to characterize structural evolution of the bulges and disks of galaxies. However, just as any morphological description of galaxies (e.g. Hubble types) depends on the wave band, intrinsic structural parameters are also expected to vary with wavelength due to stellar population and dust extinction effects. Thus, multi-wavelength information is required for any accurate description of galaxian structural parameters.

Physical differences in the shape and size of bulges among galaxies are also expected depending on how they were formed. Formation by accretion processes (e.g. major/minor mergers) can account for steeply rising de Vaucouleurs light profiles in the central parts of galaxies (e.g. van Albada 1982), while secular evolution would yield exponential distributions, with or without a core, of the central light. The formation of small bulges is indeed largely attributed to secular processes and redistribution of disk material (see § 5.3). The present study is a natural extension of de Jong’s (1996a) structural analysis of 86 face-on spirals with BVRIHK imaging, and Graham’s (2001) re-investigation of de Jong’s data.

De Jong’s 1D and 2D B/D decompositions established significant parametric variations at different wavelengths. Given the intrinsic limitations of the data modeling (i.e. over-determination of the parameter space), his B/D fits also used a fixed Sérsic $`n`$ parameter (see § 3.3.6), but limited to values of $`n=1`$, 2, and 4 bulges. De Jong’s analysis, and that of Courteau et al. (1996) who performed 1D profile decompositions for 290 r-band luminosity profiles, supported the notion of exponential bulges and disks and a tight correlation of B/D scale parameters in late-type spirals. Evidence for this correlation was challenged by Graham & Prieto (1999) but later validated by Graham (2001) who re-modeled de Jong’s thesis sample with a 1D B/D decomposition technique777Graham’s B/D analysis uses an unconstrained (floating) Sérsic shape parameter (see § 4.5.1).. His results support a range in the Sérsic shape parameter from large ($`n23`$) to small ($`n0.5`$) values for early- to late-type spirals. Aware of the inadequacy of basic B/D decompositions in fitting the bulge shape parameter due to poor data resolution and strong covariances with other bulge parameters (de Jong, 1996b; Broeils & Courteau, 1997, see § 3.3.6), we were compeled to revisit this issue with our own well-tested technique and a more extensive data base.

Our approach involves B/D decompositions with fixed $`n`$ values that sample the full parameter space of spiral bulges, from $`n=0.1,0.2,\mathrm{},4.0`$. The fit solutions are filtered out on basis of relative $`\chi ^2`$ and criteria based on our simulations (§ 3.4). We are concerned below with the derivation of robust B/D parameters for each galaxy profile. We compare our results with Graham (2001) and others, and test for any B/D parameter correlations in § 5.

The following is based exclusively on results from 1D B/D decompositions. The fact that we do not model non-axisymmetric shapes (bars, rings, oval distortions) lessens the need for more computationally intensive 2D B/D decompositions, as our simulations showed no improvements using the 2D over the 1D decomposition method for axisymmetric stucture.

### 4.2 B/D Initial Estimates

In order to determine the range of best fitted bulge and disk parameters, we need to assist the minimization program in finding the lowest possible (data$``$model) $`\chi ^2`$. From analysis of our mock images and profiles, we have found that any reasonable initial estimates for the disk parameters yields a robust solution. We base our initial estimates for the disk parameters $`h`$ and $`\mu _0`$ on the “marking the disk” technique, where the linear portion of a luminosity profile is “marked” and the selected range is fit using standard least squares techniques to determine its slope. Clearly, the resulting fits are very sensitive to the adopted baseline. We tested various choices for the fit start and end points for our galaxy profiles including: full profile fit, starting points of 0.2$`r_{max}`$ and 0.4$`r_{max}`$ out to $`r_{max}`$, and a fixed baseline shifted along the length of the profile and tracking 8 different locations. Additionally, we also tested the “moments method” of Willick (1999). The discrepancies between the different fits are large; ($``$10% on average and up to $``$100% for the worst cases), but we found that the 0.2$`r_{max}`$ to $`r_{max}`$ baseline yielded the most reliable fits (as judged by eye). The inner boundary is chosen to exclude the major contribution of a putative bulge or Type-II dip and $`r_{max}`$ is the radius at which the surface brightness error has systematically reached values greater than 0.12 mag arcsec-2 (beyond which the data become too noisy to be trusted). The fits using the 0.2$`r_{max}`$ to $`r_{max}`$ baseline provided fits that were more than adequate as initial estimates for the disk parameters in the decompositions.

Flexibility in the choice of bulge initial parameters is, however, only afforded outside a certain range of bulge sizes relative to the seeing disk. Moreover, observed galaxy profiles show significantly more variety than the idealized profiles from which these conclusions were drawn (e.g. we did not model Type-II galaxies, or the presence of strong spiral features). Accordingly, we explore three different sets of initial bulge parameter estimates to protect against local minima in the parameter space. Initial bulge effective parameters were determined from:

* Subtraction of the disk fit (based on the “marking the disk” technique) from the original profile leaving only the bulge light. $`r_e`$ is then computed non-parametrically from the data by summing up the light up to the radius which encloses half the total light of the bulge. Thus $`\mu _e=\mu (r_e)`$.
* $`r_e=0.15h`$ and $`\mu _e=\mu _0`$, where $`h`$ and $`\mu _0`$ are determined from the “marking the disk” technique.
* $`r_e=(b_n/\mathrm{log}(e))0.15h`$ and $`\mu _e=(b_nb_{n=1})+\mu _0`$.

The second set of initial estimates was motivated by the bulge/disk structural correlation found by Courteau et al. (1996) We added the third set of initial estimates which attempt to scale $`r_e`$ and $`\mu _e`$ more appropriately to the different values of $`n`$. No specific set of initial estimates worked better for all cases, though the $`3^{rd}`$ method may be the least attractive. It failed to provide reliable solutions (i.e. the fit failed or the $`\chi ^2`$ values were large) in most cases, but in a few cases it also yielded the only viable solution.

#### 4.2.1 Seeing and Sky Treatment

“Bulges” of late-type spirals are small and their luminosity profiles can be severely affected by atmospheric blur. In principle, if the blurring from the atmosphere (seeing) can be measured accurately, it can also be corrected by Fourier deconvolution. In practice, however, deconvolution amplifies noise, and the seeing FWHM is subject to measurement errors. In § 3 we used extensive simulations with a wide range of input parameters and various values of $`n`$ to derive a space of recoverable parameters under specific seeing conditions, accounting for the typical measurement errors of our data. Seeing is accounted for by convolving the model profiles with a Gaussian PSF (Eq. 6) whose dispersion is measured from field stars. In order to account for seeing measurement errors, each profile is modeled with three different values of the seeing FWHM: the nominal measured value and $`\pm `$15% of that value. A mean seeing uncertainty of $`\pm `$15% was used rather than the individual errors per measurement as these fluctuate greatly due, in large part, to the different number of stars in each measurement.

Sky subtraction errors, of order $`0.5`$–1.0% in the optical and $`0.02`$% in the H-band, were also examined carefully (§ 3.3.5). The sky brightness measurement error is accounted for in B/D decompositions by using three different sky levels: the nominal measured value and $`\pm `$0.5% (optical) or $`\pm `$0.01% (H-band), of that value.

Each profile is thus reduced 3 times for each different combination of $`r_e`$ and $`\mu _e`$ initial estimates, times 3 seeing FWHM values, times 3 sky values, and times 40 different fixed values of $`n`$, for a total of 1080 decompositions per profile.

### 4.3 Data Filtering

The 1080 decompositions for each profile are first vetted on the basis of structural criteria determined from our simulations (§ 3.3). A decomposition is deemed acceptable if it meets the following criteria:

* $`r_e(0.3)^{1/n}\text{FWHM}`$ and $`r_e\{\begin{array}{cc}0.75n+1.15\hfill & \text{ for }n1.0\hfill \\ \mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.2}n+0.2\hfill & \text{ for }n>1.0\hfill \end{array}`$
* $`B/D<5`$
* $`h<15`$ kpc ; $`r_e<50`$ kpc
* $`r_e/h<1`$

The first constraint is derived from our simulations and effectively eliminates small bulges whose sizes are comparable to, or smaller than, the seeing disk. The remaining constraints are based on physical considerations and help eliminate solutions with small $`\chi ^2`$ values but unrealistic parameters for late-type galaxies. Note, however, that these physical constraints are rather generous and do not contribute any subjective bias.

The successful decompositions are then ranked on the basis of two indicators: (a) a global $`\chi ^2`$, $`\chi _{gl}^2`$, computed for the full SB profile from $`r=0`$ to $`r_{max}`$; and (b) an inner $`\chi ^2`$, $`\chi _{in}^2`$, which includes only the central regions of the galaxy from $`r=0`$ to twice the radius, $`r_{b=d}`$, where the intensities of the fitted bulge and disk are equal (see § 3.2). $`\chi _{in}^2`$ was adopted to increase the sensitivity of the goodness-of-fit indicator to the bulge area888Note that our algorithm minimizes the $`\chi _{gl}^2`$ only. The $`\chi _{in}^2`$ is calculated and used as a discriminator only after the algorithm has converged.. The radius $`r_{b=d}`$ is clearly a function of the bulge shape and may change from small to large $`n`$ (see e.g. Fig. 2). Thus we use a $`\chi ^2`$ per degree of freedom to remove any dependence of the normal $`\chi ^2`$ to a changing $`r_e`$. Because of the presence of spiral arms and other non-axisymmetric features which we do not attempt to model, the reduced $`\chi ^2`$ is always, in principle, greater than unity. However, some of our solutions may have $`\chi ^2`$ values less than unity indicative of an over-determined system (correlated parameters), or over-estimated errors.

We first rank the solutions according to their $`\chi _{gl}^2`$ and preserve only the better half. The reduced set is then ranked according to $`\chi _{in}^2`$ values and the bottom half of the distribution is discarded. This process is iterated at least twice, or until we reach 50 or fewer solutions. Solutions with $`\chi _{gl}^2`$ greater than 50 in this final subset are discarded.

Ideally, the minima for the distributions of $`\chi _{gl}^2`$ and $`\chi _{in}^2`$ values should agree to a common value of $`n`$, but differences may exist. We search the final $`50`$ solutions for a common solution, starting at the minima of each $`\chi ^2`$ distribution. If the $`n`$ values corresponding to the two $`\chi ^2`$ minima do not agree, the $`n`$ values for the next smallest $`\chi ^2`$ values are compared (with the lower $`\chi ^2`$ values and with each other), and this process is iterated up to three times until a match is found. If this process did not converge, i.e. there is no true minimum in the (($`\chi _{gl}^2`$$`\chi _{in}^2`$) - $`n`$ space), a final solution is chosen corresponding to the minimum value of ($`\chi _{gl}^2`$/min($`\chi _{gl}^2`$) + $`\chi _{in}^2`$/min($`\chi _{in}^2`$)).

Figs. 7 & 8 show examples of the distributions of ($`\chi _{gl}^2`$$`\chi _{in}^2`$) versus $`n`$ where $`\chi _{global}^2\chi _{gl}^2`$/min$`(\chi _{gl,filt}^2)`$ and $`\chi _{inner}^2\chi _{in}^2`$/min$`(\chi _{in,filt}^2)`$, where min$`(\chi _{filt}^2)`$ is the minimum $`\chi ^2`$ value from the set of ($`50`$) filtered solutions. Note that these minima do not necessarily correspond to the lowest value of the respective distributions from all 1080 solutions, as the initial absolute minima may have been filtered out (i.e. a poor combination of $`\chi _{gl}^2`$$`\chi _{in}^2`$ for a given solution). Thus, the normalized $`\chi ^2`$’s may be less than one (as is easily seen in the leftmost plot of Fig. 8). In these plots, the left panels show the $`\chi ^2`$ distributions for all 1080 decompositions, while the right panels display only the $`50`$ solutions remaining after the iterative filtering scheme described above.

Fig. 7 highlights the sensitivity of our technique for two V-band observations of UGC 929 taken under different seeing/sky conditions. The left figures show a fairly well-behaved solution favoring $`n=0.6`$ and the figures on the right plot show a rather messy solution favoring $`n=0.8`$. The seeing conditions were worse and the sky was much brighter for the observation shown on the right which could explain the noisy distributions of both the $`\chi _{gl}^2`$ and $`\chi _{in}^2`$.

Fig. 8 shows two different behaviors of $`\chi _{gl}^2`$ for profiles with very well-behaved $`\chi _{in}^2`$. The plot on the left for our UGC 784 B-band profile illustrates the need for an additional, more discriminating statistic for the bulge region. Decompositions based solely on the $`\chi _{gl}^2`$ goodness-of-fit indicator may result in fits, like the one shown on the right side of Fig. 9 (dashed-dotted blue line). However, Fig. 8 a) clearly shows that the fit with $`n=0.6`$ is a far superior match to the bulge shape, as indicated by the $`\chi _{in}^2`$ behaviour.

The final step of our filtering procedure entails a visual inspection of the final decompositions. The criteria for user examination include information from multiple exposures and multi-band reductions for a given galaxy. Profiles and/or solutions with the following pathologies were eliminated from the final sample:

* disk profiles that are too short for proper fitting
* no obvious, extended, underlying exponential structure for the disk (occurs predominantly in Type-II profiles)
* unphysically large fitted bulge
* unrealistic disk fit for Type-II profiles. The fit is tipped below the true disk to account for the Type-II dip near the bulge-disk transition region leading to scale lengths that are biased high (e.g. see Fig. 12 for UGC 12527 for examples of “bad” fits which were eliminated from the final sample)
* large deviations between solutions for multiple observations of a given galaxy.

Not surprisingly, most of the eliminated profiles are Type-II systems. We caution that even the Type-II profile decompositions that survived the full sorting process may not provide the ideal description of their complex surface brightness distributions. Cleary these Type-II profiles cannot be properly modeled with just a Sérsic bulge and exponential disk. Out of 523 images/profiles, a total of 341 passed our acceptance criteria.

#### 4.3.1 Preferred Sky and Seeing

Histograms of the preferred seeing FWHM and sky offsets for all decompositions in all four bands are shown in Fig. 10. Typically, a lower sky brightness level is preferred by our algorithm. In some cases, this could be explained by an over-estimated sky level, but it may also be due in part to profiles with truncated outer disks as in Fig. 13. Our program prefers a slightly larger seeing FWHM than measured. This could be the result of an under-estimated FWHM, or perhaps differences between the idealized Gaussian model and the real seeing PSF. Solutions with variable sky/seeing estimates were retained in the final solution set for assessment of parameter errors.

#### 4.3.2 Effect of $`r_{max}`$

Of significance to the fit results is the maximum radius used in the decompositions. We have used the full observed profile out to radii where the surface brightness errors reached above 0.12 mag arcsec-2 . There is no absolute definition to the edge, $`r_{max}`$, of a disk and a different selection could yield different results. To test the sensitivity of our parameter determinations to the chosen value of $`r_{max}`$, we re-decomposed the profiles as described above, but with a fit baseline extending only to $`0.75\times r_{max}`$. A comparison of the results from the two techniques, prior to eyeball filtering, shows good agreement and we chose to keep the larger baseline to avoid discarding good data.

### 4.4 Decomposition Examples

There is no room for a full display of our catalog of final decompositions, but a few examples are shown in Figs. 1114. The full catalog of decomposition plots is available upon request from the authors.

In these figures, the solid black circles are the data points, the black dots show the sky error envelope (from the measured sky error), the dashed and dashed-dotted lines show the disk and bulge fits respectively, and the solid line is the total (bulge+disk) fit. The fits are all seeing-convolved using the best selected seeing values. The bottom panel shows the fit residuals where $`\mathrm{\Delta }\mu (r)`$ represents the (data$``$model). Fig. 11 shows an example of the quintessential Type-I profile at all wavelengths. Fig. 12 shows a Type-II/Transition galaxy whose Type-II signature significantly weakens from the optical to the infrared. Fig. 13 shows a Type-I profile with an outer truncated disk. Such decompositions will presumably favor an under-subtracted sky in attempt to align the inner and outer parts of the disk. Here is an example where our procedure with an over/under-estimation of the sky and an infinite exponential disk model may not be adequate since the outer disk truncation appears real (as detected in all four bands). We also show an example of a nearly bulgeless system in Fig. 14. Our sample is divided into 52 Type-I, 53 Type-II and 16 transition systems, of which 18 truncated and 7 bulgeless disks are identified999The fraction of a given galaxy class should be interpreted with care since our sample is not volume-limited..

### 4.5 Distribution of the Sérsic $`n`$ parameter

Fig. 15 shows histograms of the Sérsic $`n`$ parameter for all the final fits (left) and good fits only (right) after user-examination as described above. The distribution of $`n`$ has a definite range, implying that not all late-type bulges are best described by an exponential profile, but the mean value is very close to one. This result agrees with Graham (2001) and recent N-body simulations of galaxy evolution (§ 5.3).

#### 4.5.1 Floating Sérsic $`n`$

In § 3.3.6 we showed that resolution limitations prevented stable fitting of the Sérsic shape parameter $`n`$ as a free parameter. To illustrate the effect a floating $`n`$ can have on fitted parameters we re-decomposed all of our galaxy profiles leaving $`n`$ as a free parameter (e.g. akin to Graham 2001). The results are shown in Fig. 16 for three different initial guesses for $`n`$ (0.2, 1.0, and 4.0).

The histograms of the resulting distributions of $`n`$ reveal a strong bias towards the adopted initial estimate. All 3 distributions show a large peak at $`n=0.1`$, indicative of poor bulge fits. The histogram for the $`n=1.0`$ initial estimate looks somewhat similar to our own constrained solution (Fig. 15), but this is somewhat fortuitous given the closely-exponential nature of spiral bulges. Note also the non-Gaussian tail in Fig. 15 is not reproduced in Fig. 16 for the $`n=1`$ initial estimate case. Fig. 17 shows a comparison of the $`\chi ^2`$ values from the floated versus fixed $`n`$ solutions. (When no suitable fit was found, all parameters were set to 0 as indicated by the points lying on the axes; note the large number of fit failures in the floated $`n`$ case.) Note the large discrepancies in $`\chi _{in}^2`$ and $`\chi _{gl}^2`$ between the two methods. Thus, while the final distributions for the $`n=1`$ initial estimate and our constrained $`n`$ procedure look similar, significant differences may exist between individual decompositions.

### 4.6 Error of a Single Measurement

An important feature of any decomposition technique is the stability of the final results for repeat observations of a given system. Our sample has 50 profiles for which multiple (two to four) observations exist, allowing for a direct measure of the reliability of our decompositions. Table 1 gives the mean and mean standard deviation of the five model parameters from repeat observations with,

$$\overline{\sigma _x}=\frac{\underset{N}{}\left\{\left[\frac{\underset{n}{}(x_n\overline{x})^2}{n1}\right]_N^{1/2}\right\}}{N}$$
(23)
where $`x`$ is the fit parameter, $`n`$ is the number of observations for a given profile, and $`N`$ is the number of profiles with repeat observations. The average errors from repeat observations of Type-I profiles are $`\pm \mathrm{\hspace{0.17em}14}`$% for $`n`$, $`\pm \mathrm{\hspace{0.17em}0.2}`$ mag arcsec-2 for $`\mu _e`$, $`\pm \mathrm{\hspace{0.17em}13}`$% for $`r_e`$, $`\pm \mathrm{\hspace{0.17em}0.05}`$ mag arcsec-2 for $`\mu _0`$, and $`\pm \mathrm{\hspace{0.17em}3}`$% for $`h`$. Clearly, determinations of disk parameters are much more stable than for bulges. Error terms quoted below correspond to the 1-$`\sigma `$ deviation, unless otherwise noted.

### 4.7 Comparison with Other Authors

The overlap between our sample and de Jong’s thesis catalog (de Jong & van der Kruit, 1994) amounts to only 3 galaxies. Direct comparison of our SB profiles shows excellent zero-point and overall shape agreement (Paper II); however, our B/D decompositions differ somewhat, as shown in Table 2 (note that de Jong uses fixed $`n=1`$). Also shown in that table are decomposition parameters for the same galaxies by Graham (2001) (same data as de Jong, but using a range of $`n`$). We find scale length differences at the 10% level with de Jong and Graham, consistent with, or slightly better than, typical variations between different authors (Knapen & van der Kruit, 1991). A comparable dispersion is measured between the scale lengths of Graham and de Jong based on 82 R-band profile decompositions. Graham’s scale lengths are, on average, smaller for small galaxies and larger for big galaxies (apparent size) than de Jong’s. We find systematically larger disk scale lengths than de Jong (based on only 7 comparisons.) We verified that sky under/over-estimates cannot account for any difference with de Jong. De Jong’s algorithm gives more weight to the outer part of the disk, possibly explaining the shorter disk scale lengths. For profiles with outer truncated disks or Type-II decrements, greater weight in the outer parts favors the outer disk curvature and thus steeper disk fits.

Bulge parameters between us and Graham match reasonably well for the first two galaxies but differ substantially for UGC 3140. It is however difficult to establish trends based on just 3 comparisons. We can, instead, broadly compare our respective distributions of Sérsic $`n`$ with morphological type. This is done in Fig. 18 for comparison with Fig. 10 of Graham (2001). The general features are similar, but we find a wider range of $`n`$ values for the later-types possibly due to the larger number of Scd/Sd galaxies in our sample. Another favorable comparison of bulge parameters with Graham is shown in Fig. 20 (see § 5.2.1).

We also find an overlap of two galaxies, NGC 3512 and NGC 7782, with the sample of Baggett et al. (1998). As with de Jong and Graham, disk parameters agree within 10%. Bulge parameters from Baggett et al. are missing for NGC 3512, and those listed for NGC 7782 (both V band) differ quite substantially from ours. These authors find $`\mu _e=10.88`$ and $`r_e=0.2`$ for a de Vaucouleurs bulge and we have $`\mu _e=20.1`$ and $`r_e=3.2`$ for n=1 (best fit) or $`\mu _e=25.2`$ and $`r_e=83.8`$ for n=4 (very bad fit with high reduced $`\chi ^2`$). Note that seeing estimates were comparable. Surprisingly, their $`\mu _e`$ is nearly 10 magnitudes brighter than their $`\mu _0=20.3`$ (we also find $`\mu _0=20.3`$)! We find this pathology in a number of their bulge decompositions (see e.g. their Fig. 2) where the models often overshoot the data at the center.

We conclude this section by noting that disk scale lengths between us and other authors differ at the 10% level. Our bulge parameters are also qualitatively consistent with those of Graham.

## 5 Discussion

Simulations of galaxy profiles and images (§ 3) and careful B/D decompositions (§ 4) have led to a final set of structural parameters for late-type spiral galaxies (Table 4 in Appendix B). These data can now be examined for intrinsic structural variations and sensitivity to dust and stellar population effects. The outline of this section is as follows: First, we verify in § 5.1 that our solutions are not affected by projection effects. We then discuss in § 5.2 B/D parameter variations both in the context of profile type differences and wavelength dependence. In light of existing limitations in our modeling of Type-II profiles, our conclusions will be based mostly on properties derived from Type-I profiles. These will enable us to examine the viability of secular evolution models for disk galaxies (see § 5.3).

### 5.1 Inclination Dependence

In order to test for projection effects, we plot the distributions of $`\mu _e`$ and $`r_e`$, as well as disk $`\mu _0`$ and $`h`$ as a function of ellipticity, $`\epsilon =1b/a`$, in Fig. 19. The surface brightnesses are only corrected for Galactic extinction and cosmological dimming (as in § 2); thus the $`\mu _0`$ and $`\mu _e`$ values should be considered as upper limits (i.e. effective brightnesses are too low). No trends with ellipticity are seen, including the Sérsic $`n`$ parameter and ratio of disk scale lengths (not plotted). Furthermore, Types I, II and Transition are not confined to any particular inclination range showing that the Type-II phenomenon is not an accentuated feature due to line-of-sight extinction (e.g. Type-II galaxies are not preferentially inclined with the plane of the sky).

### 5.2 Bulge/Disk Parameters

Table 3 shows the range of fitted parameters at BVRH wavelengths for all galaxy profile types (Type-I, II, and Transition). The number of Type-II and Transition systems included in this table (e.g. only 4 decompositions for Transition galaxies in the B-band) is drastically reduced from our original distribution as many of them did not pass our validity criteria (§ 4.3). Note that the parameters for Transition profiles at H-band broadly match those of Type-I’s at that wavelength.

The Sérsic shape parameter $`n`$ for Type-I galaxies is near unity, within the errors, for all wavelengths. Thus, we advocate that the natural, intrinsic distribution of the Sérsic $`n`$ parameter for late-type spirals has a mean near 1.0 (with $`\sigma _n0.4`$; see Fig. 15). By all accounts, bulges of late-type spirals are well-approximated, on average, by a pure exponential (luminosity/mass) density distribution.

The distributions of Sérsic $`\mu _e`$ and $`r_e`$ are broad, indicative of the range of bulge types in our sample. Effective radii are typically less than 1 kpc. Those of Type-II profiles are even smaller and seemingly better determined than Type-I’s but this is predominently an artifact of our limited 2-component modeling. Examination of Type-II profile fits shows that the model disk is typically shallower (than the “true” disk), as the fit accounts for the fainter bulge/disk transition dip, and bulge effective radii are naturally confined to a smaller range.

The distributions of disk scale lengths and their ratios show a clear decreasing trend as a function of wavelength (as noted by de Jong (1996b)). This statistically significant effect, detected for all profile types, can be explained either by a high concentration of older stars and/or dust in the central regions of the disk. Absorption by dust alone can account for the scale length ratios that we measure (see e.g. Evans (1994), Fig. 5). Evans’ models do not consider scattering but for the nearly face-on galaxies considered here, its effects are negligible (Byun et al. 1994; de Jong (1996c)). The color gradient analysis of de Jong (1996c) using stellar population and dust extinction models suggests, however, that dust and metallicity play a minor role but that age is be the dominant factor. Preliminary analysis of our photometric data with the latest stellar evolutionary and dust models (MacArthur et al., in prep.; hereafter Paper III) suggests a combination of effects. The interpretation of color gradients is non-trivial and may ultimately require a full spectroscopic investigation to convincingly disentangle the effects of age, dust, and metallicity.

#### 5.2.1 B/D Scale Ratios

In Fig. 20, we plot $`r_e`$ vs. $`h`$ for our Type-I decompositions (solid symbols) and those of Graham (2001) for de Jong’s BRK Type-I SB profiles (open symbols). This figure provides the basis for a renewed discussion of the suggestion by Courteau et al. (1996) of structural coupling between the bulge and disk of late-type galaxies. The large dispersions in the $`r_e^\lambda `$ and $`h^\lambda `$ (see Table 3) nearly cancel out to yield tighter $`r_e/h`$ correlations. For Type-I profiles we find $`r_e/h0.22\pm 0.09`$ at all wavelengths, corresponding to $`h_{\mathrm{bulge}}/h_{\mathrm{disk}}=0.13\pm 0.06`$ for $`n=1`$. This result is also borne out in the H-band Transition profiles (see Table 3). For comparison, Courteau et al. (1996) found $`h_{\mathrm{bulge}}/h_{\mathrm{disk}}0.10\pm 0.05`$ (or $`r_e/h=0.15\pm 0.08`$)101010The study of Courteau et al. (1996) combined the r-band decompositions of Broeils & Courteau (1997) and the K-band decompositions de Jong’s (1996a) thesis study but no distinctions were made between Type-I and Type-II profiles.. The latter result is also in agreement with Graham (2001) who finds $`r_e/h=0.2`$ (no quoted dispersion, but it is somewhat larger than ours judging from Figs. 20 and 21) for early and late-type spirals in the K-band. This is consistent with a scenario where bulges of late-type spiral galaxies are more deeply embedded in their host disk, than earlier-type bulges. In such an “iceberg” scenario (e.g. Graham 2001), bulges and disks can preserve a nearly constant $`r_e/h`$ but show a great range of $`\mu _e`$ for any given $`r_e`$.

In Fig. 21, we show $`r_e/h`$ as a function of morphological type from our (solid symbols) and Graham’s (open symbols) decompositions. A mild trend with Hubble type is seen with $`r_e/h=0.200.013(T5)`$ ($`1\sigma =0.09`$), ranging from $`r_e/h0.20`$ for late-type spirals to $`r_e/h0.24`$ for earlier types. Comparison of our and Graham’s decomposition parameters in Table 2 shows that large deviations may exist, thus only our data points were included in the fit of $`r_e/h`$ vs $`T`$ above. Data for different bands scatter evenly about the mean line. More data at earlier and later types would be needed to firm up this trend. It is nonetheless remarkable that early and late-type systems are described by very similar scaling relations, thus suggesting comparable formation and/or evolution scenarios.

### 5.3 Test of Secular Evolution in Late-Type Spirals

This work has provided confirmation of two important structural signatures of spiral galaxies which must be addressed by models of structure formation:

* The underlying surface brightness distribution of late-type spirals has a range for the Sérsic $`n`$ parameter from 0.1–2, but is best described, on average, by a double-exponential model of bulge and disk, such as found in Type-I profile galaxies.
* Bulges and disks of late-type spirals are coupled, with $`r_e/h=0.22\pm 0.09`$, or $`h_{\mathrm{bulge}}/h_{\mathrm{disk}}`$ $`=0.13\pm 0.06`$, at all wavelengths. A mild trend with Hubble type is also detected with a range $`r_e/h`$ 0.20–0.24, from late to early-type spirals.

The first result describes the large-scale appearance of bulges. Analyses of HST images have shown that a significant fraction of bulge nuclei have power-law profiles ($`r<500`$ pc; e.g. Phillips et al. 1996; Balcells 2001) and host a central compact source (Carollo 1999). The extent of these nuclear sources ($`<`$ 0$`\stackrel{}{\mathrm{.}}`$3 for $`cz<2500`$ km s-1) is smaller than our images’ pixel size and smoothed out by seeing. We thus ignore their effects on the bulge light profile in this analysis, but caution that our bulge parameters are to be considered upper limits if a significant nuclear component is present.

A natural interpretation of the near constancy of B/D size ratios in late-type spirals is that their bulges formed via secular evolution of the disk. This scenario is possible if disks are bar-unstable, which can be triggered by the global dynamical instability of a rotationally supported disk or induced by interactions with a satellite and if significant angular momentum transport is feasible (e.g. Martinet 1995; Combes 2000; see the collection of papers in Carollo, Ferguson & Wyse 1999 for comprehensive reviews). For bar-unstable disks, in particular to vertical deformations, the inner disk material is heated up to 1–2 kpc above the plane into a “bulge” via resonant scattering of the stellar orbits by the bar-forming instability. This in turn, catalyzes funneling of disk material into the central regions and generates outward transport of disk material in the outer parts. Gas flows must also be invoked to explain the higher spatial densities of bulges compared to the inner disk. Such a model is expected to produce correlated scale lengths and colors between the disk and its central regions, as observed (e.g. Terndrup et al. 1994; Peletier & Balcells 1996; Courteau 1996b). A “bulge-like” component with a nearly exponential profile is expected from non-axisymmetric disturbances that induce inward radial flow of disk material (Pfenniger & Friedli 1991; Zhang & Wyse 2000, and references therein). The longer the disk-bar heating interaction, the greater the extent of the disk exponential profile (Valenzuela & Klypin 2002). The evolving exponential attractor is an empirical result well established in simulations, but it lacks, at present, a theoretical explanation (Pfenniger 1999).

Although a bar can grow spontaneously ($`\genfrac{}{}{0pt}{}{_<}{^{}}`$20 Myr) from small scale fluctuations in the inner disk, an external finite perturbation can catalyse its growth. However, collisionless mergers seem unsuited to growing the exponential bulges of present-day late-type spirals, though they may contribute to the increase in Sérsic $`n`$ parameter seen toward earlier types in proportion to the accreted satellite mass (Barnes 1988, Aguerri et al. 2001). The spontaneous or triggered formation of bars also suggests that the Hubble type of galaxies can change well after the formation of the disk (Pfenniger 1999). All of our bulges are smaller than a disk scale length and could be created by purely bar-related processes. Instead, accretion of galaxy satellites is required to make bigger bulges, either before or after formation of the host disk.

Secular evolution models of stellar and gaseous disks, especially through cosmologically-motivated three-dimensional N-body simulations, have seen significant developments in the last decade. For example, the cold dark matter (CDM) hierarchical hydrodynamical simulations by Sáiz et al. (2001) and Scannapieco & Tissera, in prep.  show that secular processes can occur naturally during the formation of spiral disks and play an important role in the regulation of star formation and the determination of the dynamical and structural properties of these systems. On average, the simulated disk systems are shown to be characterized by a double exponential profile which naturally emerges within the hierarchical clustering scenario. These results are based on a stellar formation process implemented in such a way that it succeeds in forming compact bulges that stabilize disk-like structure allowing the conservation of an important fraction of their angular momentum during the violent phases of their assembly. Fig. 22 shows the distribution of final Sérsic $`n`$ parameters for relaxed present-day late-type disks by Scannapieco & Tissera, in prep. ; their models reproduce our results (Fig. 15) very nicely. The double-exponential structure of bulge and disk may not always be the final relaxed state of an object, but whenever $`n1`$, the B/D scale ratio $`h_{\mathrm{bulge}}/h_{\mathrm{disk}}`$ takes its nominal value of 0.15. These models do not have bulges with $`n<0.7`$, possibly due to limited resolution and/or excessive angular momentum transfer that supernova feedback could help prevent.

Simulations by Pfenniger (2002; private comm.) of self-gravitating disks forming bars which may later dissolve into a bulge-like component also show a nearly universal ratio $`r_e/h`$, in agreement with observed values, which is related to the stellar dynamics of the barred system (i.e. relative position of the vertical to horizontal resonances). The bar length is related to the initial rising part of the rotation curve (yielding a scale), and the corotation of bars is proportional to their length. The corotation fixes the positions of the other resonances, which in turn fix the maximum extension of bulges made from resonant heating, as indeed the vertical resonances are strong only within the bar. This mechanism would thus set the upper limit for the allowed range in $`r_e/h`$.

The $`N`$-body simulations of Aguerri et al. (2001), which consider the growth of galactic bulges by mergers, also suggest that the final B/D scale ratio $`r_e/h`$ does not scale with the B/D luminosity ratio. These authors show that the disk scale length $`h`$ can increase from 15% (low mass retrograde satellite) to 65% (high mass direct satellite) while $`r_e/h`$ would decrease from 0.21 to 0.14, or 33%, in the most extreme case. One can thus infer $`r_e/h=0.17\pm 0.03`$, independent of B/D luminosity ratio, in good agreement with our findings. Their simulations are however limited to a small range of initial $`r_e`$ and a more complete investigation with a broad range of $`r_e`$ and $`h`$ values is needed to establish the fundamental nature of the B/D scale ratio.

### 5.4 Type-II Profiles

The above scenarios for secular evolution naturally produce the double-exponential character of the bulge and disk radial luminosity profiles for late-type systems. However, over half our sample of 121 late-type spiral galaxies show strong deviations from this simple two-component description. Other authors (Kormendy, 1977; Baggett et al., 1998) have considered inner disk truncation (plus de Vaucouleurs bulges) as an alternative to modeling Type-II light profiles, with

$$I_{disk}(r)=I_{}\mathrm{exp}\left\{\left[r/r_{}+(r_{hole}/r)^n\right]\right\}$$
(24)
where $`r_{hole}`$ is the truncation radius and $`n3`$. As discussed in § 3.2, we do not consider this approach at the present, but its potential merits should not be overlooked.

N-body simulations (e.g. Norman, Sellwood, & Hasan 1996; Valenzuela & Klypin 2002) reproduce Type-II surface density profiles as a result of the redistribution of central stars into a ring by a bar-like perturbation. Approaching the centers where the component called bulge and the component called exponential disk overlap, one cannot tell, in these simulations, if a star or particle belongs to which component. Galaxy centers may recurrently move from a barred to an unbarred phase and undergo continuing bulge building as the bars dissolve111111Only progressively larger bars in the centers of exponential bulges would be allowed to form in a recurring scenario due to the disrupting dynamical effect of a growing nucleus (Rix 1998, as reported in Carollo 1999). (Norman et al., 1996). Thus, the paucity of barred galaxies in our sample does not preclude bar-induced effects as a possible explanation for Type-II profiles (e.g. Gadotti & dos Anjos (2001)). Pre-existing bars may simply have dissolved. For example, out of 8 barred-classified galaxies in our sample, 6 have Type-II profiles thus lending some credence to the bar-lens scenario. On the other hand, the most strongly barred galaxies in the Shellflow sample of $`300`$ bright late-type galaxies (Courteau et al. 2000) have mostly Type-I profiles, indistinguishable in shape and global properties from the profiles of unbarred Type-I galaxies (Courteau et al., in prep.). ¿From an inhomogeneous sample of 167 spiral galaxies, Baggett et al. (1996) find only a weak tendency for barred galaxies to have a higher occurence of Type-II profiles. The link between Type-II profiles and barred galaxies is thus unsecured at present.

Type-II profiles may also be explained by extinction effects in the disk. Increased opacity towards the central disk can cause a depression in the luminosity profile, especially at shorter wavelengths. Realistic Type-II profiles (in shape and colors) have been produced with exponential distributions of stars and dust and variable layering parameters (Evans, 1994). If dust extinction causes the inner disk profile dip, Transition galaxies would just be a case of lesser dust content, whereas bona fide Type-II systems remain optically thick, even at H-band. Using far-infrared (FIR) to B-band flux ratios, and radiation transfer models for the dust (Gordon et al., 2001), we have tested for the origin of Type-II signature as being due to extinction. The FIR/B flux ratio should be higher for the dustier systems. Unfortunately, our measured total FIR/B flux ratios are statistically identical (with large scatter) for Type-I, Type-II, and Transition galaxies (Paper III), thus thwarting any clear interpretation. The IRAS 60 and 100 $`\mathrm{\mu m}`$ fluxes have too low resolution and too large errors to separate the inner disk dust emission from the whole galaxy.

If stellar population effects are relevant (Prieto et al., 1992), age/metallicity gradients should be detected at the bulge/disk transition in Transition systems. We will further investigate the dust and/or stellar population origin of the Type-II dip in Paper III.

## 6 Summary and Concluding Remarks

This study has focused on the development of rigorous B/D decomposition techniques using a new, comprehensive, multi-band survey of late-type spiral galaxies. We examine three types of SB profiles, Freeman Type-I and Type-II, and a third “Transition” class for galaxies whose profiles change from Type-II in the optical to Type-I in the infrared. This distinction is important since Type-II and Transition profiles cannot be adequately modeled by a simple two-component model of the bulge and disk. Thus, our main results are based on Type-I profiles.

Based on extensive simulations, careful treatment of sky and seeing measurement errors, and repeat observations we are confident that systematic errors are $`20`$% for the bulge components, including the Sérsic shape parameter, and $`5`$% for disk components.

The main conclusions from our simulations and final profile decompositions are as follows:

* Simulations to determine the range of acceptable solutions for any B/D decomposition program are crucial. The reliability of bulge model parameters is limited by the relative size of the bulge and seeing disk, seeing errors, the intrinsic bulge shape, sky brightness and errors. Disk parameters are fairly robust to systematic errors, with the exception of improper bulge shapes and sky errors which can have dramatic effects on both modeled disk and bulge components.
* The Sérsic bulge shape parameter for nearby late-type galaxies shows a range between $`n=0.12`$, but, on average, their underlying surface brightness distribution is best described by a double-exponential model of bulge and disk.
* Disk scale lengths decrease at longer wavelengths, indicative of a higher concentration of older stars and/or dust in the central regions relative to the outer disk.
* We confirm and reinforce the result of Courteau et al. (1996) of a structural coupling between the bulge and disk of late-type spirals. We find $`r_e/h=0.22\pm 0.09`$, or $`h_{\mathrm{bulge}}/h_{\mathrm{disk}}=0.13\pm 0.06`$, independent of wavelength. A mild trend with Hubble type is observed with $`r_e/h=0.200.013(T5)`$ ($`1\sigma =0.09`$), ranging from $`r_e/h0.20`$ for late-type spirals to $`r_e/h0.24`$ for earlier types. These results are consistent with scenarios of bulge formation in which bulges of late-type spiral galaxies are more deeply embedded in their host disk than earlier-type bulges. Under this “iceberg” scenario, bulges and disks can thus preserve a nearly constant $`r_e/h`$ but show a great range of $`\mu _e`$ for any given $`r_e`$. The observed scale ratio is consistent with numerical simulations of self-gravitating disks and probably related to the stellar dynamics of an actual or pre-existing barred system.
* The inner brightness profile signatures of Type-II galaxies are likely explained by a combination of dust extinction and stellar population effects and perhaps linked to the occurence of a bar, but no decisive conclusion can be derived at present.

We are grateful to Marc Balcells, Eric Bell, Roelof de Jong, and Daniel Pfenniger for their comments on earlier versions of this manuscript. Cecilia Scannapieco and Patricia Tissera are also thanked for sharing their material (Fig. 22) in advance of publication. Alister Graham kindly provided tables of his profile decompositions of de Jong’s SB profiles for comparison with our and de Jong’s similar results. We also wish to thank the anonymous referee for suggestions and comments that helped improved the presentation and content of the paper. This research has made use of the NASA/IPAC extragalactic database (NED) which is operated by by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. LM and SC acknowledge financial support from the National Science and Engineering Council of Canada.

## Appendix A Functional form for the Sérsic $`b_n`$ parameter

Eq. 14 cannot be solved in explicit closed form for $`b_n`$. Many of the numerical and analytical solutions found in the literature agree well for $`n>1`$ but differ significantly for smaller values of $`n`$. Fig. 23 shows a comparison of the two most commonly used approximations (short- and long-dashed curves) with the exact solution for $`b_n`$, computed to a numerical precision of one part in $`10^7`$ for all $`n10`$ (see also Fig. 2 in Graham 1999).

As we wish to test for spiral bulges with Sérsic $`n`$’s as small as 0.1, we have adopted a formalism that is valid for all $`n`$. To maintain computational simplicity, and ensure a suitably accurate solution we found it practical to divide the curve into two segments. For all $`n>0.36`$ we use the asymptotic expansion of Ciotti & Bertin (1999) up to O$`(n^5)`$ (their Eq. 18),

$$b_n2n\frac{1}{3}+\frac{4}{405n}+\frac{46}{25515n^2}+\frac{131}{1148175n^3}\frac{2194697}{30690717750n^4}+O(n^5)$$
(A1)
which is good to better than one part in $`10^4`$ in that range. However, for $`n0.36`$ this solution diverges. Due to the rapidly changing curvature in the gamma function, and thus $`b_n`$, at small $`n`$, it would be necessary to use an unrealistic number of terms in the asymptotic expansion to achieve the desired accuracy. For $`n0.36`$, we find the best fitting polynomial of the form $`b_n=_{i=0}^ma_in^i`$ where $`m`$ is the order of the polynomial and the $`a_i`$ are the coefficients of the fit given by,

$$\begin{array}{ccccc}a_0=0.01945\hfill & a_1=0.8902\hfill & a_2=10.95\hfill & a_3=19.67\hfill & a_4=13.43.\hfill \end{array}$$
(A2)
This fit is accurate to better than two parts in $`10^3`$. The wiggles in the dotted curve in Fig. 23 result from the polynomial nature of the fit and limited numerical precision where the gamma function approaches infinity.

## Appendix B Decomposition Results for the Type I Profiles

Table 4 gives relevant photometric information and 1D B/D decomposition results for the final set of Type I galaxy profiles. Decomposition results for Type-II and Transition galaxies are available from the authors upon request, with the caution that parameters for these profile types should be interpreted with care. The entries are arranged as follows:

Column (1): (UGC number) (observation number) (passband) for each profile;

Column (2): Ellipticity, $`\epsilon (1b/a)`$. The final ellipticity (and position angle) estimates correspond to an average of those values from the five contours surrounding the best isophotal fit in the outer disk, as determined by eye. This estimate is clearly sensitive to the presence of spiral arms. The typical inclination error is $`3\mathrm{deg}`$, independent of ellipticity;

Column (3): Sky brightness in mag arcsec-2 , measured from 4 sky boxes located between the detector edges and a fair distance away from the galaxy. Typical rms sky errors, computed from the deviations of the mean sky counts in those sky boxes, are $`0.51.0`$% in the optical and 0.05% in the IR. The subscripts indicate the sky offset preferred by our selection process as described in § 4.2.1 and § 4.3 (and see Fig. 10), where “$`+`$” and “$``$” indicate 0.5% for optical and 0.01% for H-band over- and under-subtracted skies respectively. No subscript indicates that the measured sky was preferred;

Column (4): Seeing FWHM values, computed as the average of the FWHMs of all non-saturated stars measured automatically on each image frame; typically 10 to 30 measurements per image were used for each FWHM estimate. The accuracy of the seeing estimate per image is roughly 20% for the optical bands and 30% for the H-band. The subscripts indicate the seeing offset preferred by our selection process as described in § 4.2.1 and § 4.3, (and see Fig. 10) where “$`+`$” and “$``$” indicate 15% over- and under-estimated seeing FWHM respectively. No subscript indicates that the measured seeing FWHM was preferred;

The upper and lower boundaries in the remaining columns correspond to the maximum and minimum values of the $`50`$ (out of 1080 total) solutions left after filtering (see § 4.3);

Column (5): Best fit Sérsic $`n`$ bulge shape parameter;

Column (6): Bulge effective surface brightness, $`\mu _e`$, in mag arcsec-2 , corrected for Galactic extinction and cosmological redshift dimming, as described in § 2.2;

Column (7): Bulge effective radius, $`r_e`$, in arcseconds;

Column (8): Bulge effective radius, $`r_e`$, in kpc. Converted to a physical scale using the Local Standard of Rest velocity, $`V_{LG}`$ (see Paper II);

Column (9): Exponential disk central surface brightness, $`\mu _0`$ in mag arcsec-2 , corrected for Galactic extinction and cosmological redshift dimming as described in § 2.2;

Column (10): Exponential disk scale length $`h`$, in arcseconds;

Column (11): Exponential disk scale length $`h`$, in kpc. Converted to a physical scale using the Local Standard of Rest velocity, $`V_{LG}`$ (see Paper II);

Column (12): Bulge-to-disk luminosity ratio, $`B/D`$, calculated using Eq. 17 in § 3.1.