Structural properties of disk galaxies

# Structural properties of disk galaxies

I. The intrinsic equatorial ellipticity of bulges
J. Méndez-Abreu INAF-Osservatorio Astronomico di Padova, vicolo dell’Osservatorio 5, I-35122 Padova, Italy
Universidad de La Laguna, Av. Astrofísico Francisco Sánchez s/n, E-38206 La Laguna, Spain Dipartimento di Astronomia, Università di Padova, vicolo dell’Osservatorio 3, I-35122 Padova, Italy
J. A. L Aguerri Instituto Astrofísico de Canarias, Calle Vía Láctea s/n, E-38200 La Laguna, Spain
E. M. Corsini Dipartimento di Astronomia, Università di Padova, vicolo dell’Osservatorio 3, I-35122 Padova, Italy
E. Simonneau Institut d’Astrophysique de Paris, C.N.R.S., 98bis Boul. Arago, F-75014 Paris, France
July 5, 2019
###### Key Words.:
galaxies: bulges – galaxies: formation – galaxies: fundamental parameters – galaxies: photometry – galaxies: statistics – galaxies: structure
offprints: J. Méndez-Abreu
###### Abstract

Context:A variety of formation scenarios have been proposed to explain the diversity of properties observed in bulges. Studying their intrinsic shape can help to constraint the dominant mechanisms at the epochs of their assembly.

Aims:The structural parameters of a magnitude-limited sample of 148 unbarred S0–Sb galaxies were derived in order to study the correlations between bulges and disks, as well as the probability distribution function of the intrinsic equatorial ellipticity of bulges.

Methods:We present a new fitting algorithm (GASP2D) to perform two-dimensional photometric decomposition of the galaxy surface-brightness distribution. This was assumed to be the sum of the contribution of a bulge and disk component characterized by elliptical and concentric isophotes with constant (but possibly different) ellipticity and position angles. Bulge and disk parameters of the sample galaxies were derived from the band images, which were available in the Two Micron All Sky Survey. The probability distribution function of the equatorial ellipticity of the bulges was derived from the distribution of the observed ellipticities of bulges and misalignments between bulges and disks.

Results:Strong correlations between the bulge and disk parameters were found. About of bulges in unbarred lenticular and early-to-intermediate spiral galaxies are not oblate but triaxial ellipsoids. Their mean axial ratio in the equatorial plane is . Their probability distribution function is not significantly dependent on morphology, light concentration or luminosity. The possible presence of nuclear bars does not influence our results.

Conclusions:The interplay between bulge and disk parameters favors scenarios in which bulges have assembled from mergers and/or have grown over long times through disk secular evolution. However, all these mechanisms have to be tested against the derived distribution of bulge intrinsic ellipticities.

## 1 Introduction

The relative prominence of galactic bulges with respect to their disks is important in the definition of galaxy types. Therefore, understanding the formation of bulges is key to understanding the origin of the Hubble sequence.

Bulges are diverse and heterogeneous (see the reviews by kormendy93, ; wyse97, ; kormendy04, ). The big bulges of lenticulars and early-type spirals are similar to low-luminosity elliptical galaxies. Their surface-brightness radial profiles generally follow a De Vaucouleurs law (andredakis95, ; carollo98, ; mollenhoff01, , hereafter MH01). The majority of these bulges appear rounder than their associated disks (kent85, ). Their kinematical properties are well described by dynamical models of rotationally flattened oblate spheroids with little or no anisotropy (kormendy82, ; davies83, ; cappellari06, ). They have photometrical and kinematical properties, which satisfy the fundamental plane (FP) correlation (bender92, ; bender93, ; burstein97, ; aguerri05a, ). On the contrary, the small bulges of late-type spiral galaxies seems to be reminiscent of disks. Their surface-brightness radial profiles have an almost exponential falloff (andredakis94, ; dejong96, ; macarthur03, ). In some cases they have apparent flattenings that are similar or even larger than their associated disks (fathi03, ) and rotate as fast as disks (kormendy93, ; kormendy02, ). Late-type bulges deviate from the FP (carollo99, ).

Different formation mechanisms (or at least a variety of dominant mechanisms at the epochs of star formation and mass assembly) were proposed to explain the variety of properties observed in bulges. Some of these formation processes are rapid. They include early formation from the dissipative collapse of protogalactic gas clouds (eggen62, ; sandage90, ; gilmore98, ; merlinchiosi06, ) or later assembly from mergers between pre-existing disks (kauffmann96, ; baugh96, ; cole00, ). In both scenarios the disk forms after the bulge as a consequence of either a long star-formation time compared to the collapse time or a re-accretion around the newly formed bulge.

Bulges can also grow over long timescales through the disk secular evolution driven by bars and/or environmental effects. Bars are present in more than half of disk galaxies in the local universe (eskridge00, ; menendez-delmestre07, ) and out to (elmegreen04, ; jogee04, ). They are efficient mechanisms for driving gas inward to the galactic center and feed the galactic supermassive black hole (see corsini03a, , and references therein). In addition, bar dissolution due to the growth of a central mass (pfenninger90, ), scattering of disk stars at vertical resonances (combes90, ), and coherent bending of the bar perpendicular to the disk plane (raha91, ; debattista04, ; athanassoula05, ; martinez-valpuesta06, ) are efficient mechanisms in building central bulge-like structures, the so-called boxy/peanut bulges. Moreover, the growth of the bulge out of disk material may also be externally triggered by satellite accretion during minor merging events (searlezinn78, ; aguerri01, ; eliche-moral06, ) and gas infall (thakarryden98, ).

Traditionally, the study of the relations between the structural parameters of the galaxies have been used to understand the bulge formation processes, e.g., the correlation between the bulge effective radius and the scale length of the disk in many galaxy samples has always been interpreted as an indication that bulges were formed by secular evolution of their disks (see macarthur03, ). However, one piece lost in this study is the three-dimensional shape of the bulges. By studying this, one might be able to provide the relative importance of rapid and slow processes in assembling the dense central components of disk galaxies. A statistical study can provide a crucial piece of information for testing the results of numerical simulations of bulge formation for different galaxy type along the morphological sequence.

In this paper, we analyze a sample of unbarred early-type disk galaxies to derive the intrinsic ellipticity of their bulges in the galactic plane. The twisting of bulge isophotes (lindblad56, ; zaritsky86, ) and misalignment between the major axes of the bulge and disk (bertola91, ) are not possible if the bulge and disk are both oblate. Therefore, they were interpreted as a signature of bulge triaxiality. This idea is supported by the presence of non-circular gas motions (e.g., gerhardvietri86, ; bertola89, ; gerhard89, ; berman01, ) and a velocity gradient along the galaxy minor axis (corsini03b, ; coccato04, ; coccato05, ). We improve the previous works in several aspects. First, we use near-infrared images to map the distribution of the mass-carrying evolved stars and avoid contamination of dust and bright young stars. Second, we retrieve the structural parameters of the bulge and disk by applying a new algorithm for two-dimensional photometric decomposition of the observed surface-brightness distribution. Finally, we obtain the probability distribution function (PDF) of the intrinsic equatorial ellipticity of bulges by using a new mathematical treatment of the equations describing their three-dimensional shape.

The paper is organized as follow. The selection criteria of our sample galaxies and the analysis of their near-infrared images are described in Sect. 2. Our new photometric decomposition method for deriving the structural parameters of the bulge and disk by analyzing the two-dimensional surface brightness distribution of galaxies is presented in Sect. 3. The correlations between the structural parameters of the sample galaxies are discussed in Sect. 4. The PDF of intrinsic equatorial ellipticity of the studied bulges is derived in Sect. 5. Our conclusions and a summary of the results are given in Sect. 6.

## 2 Sample selection and data acquisition

Our objective was to select a well-defined complete sample of nearby unbarred disk galaxies to study in a systematic way the photometric properties of their structural components. Since these properties are strongly dependent on the dominating stellar population at the observed wavelength, it is preferable to consider near-infrared images to map the mass-carrying evolved stars and avoid contamination due to dust and bright young stars (e.g., MH01). The complete sample is drawn from the Extended Source Catalogue (XSC) (Jarrett00, ) of the Two Micron All Sky Survey (2MASS) (skrutskie06, ). Our sample consists of galaxies that meet the following requirements: (1) Hubble type classification from S0 to Sb (; de Vaucouleurs et al. 1991, hereafter RC3) to ensure that bulges are fully resolved in 2MASS images; (2) unbarred classification in RC3; (3) total band magnitude mag (2MASS/XSC); (4) inclination (RC3) to measure the misalignment between the position angle of the bulge and disk; (5) Galactic latitude (RC3) to minimize both Galactic extinction and contamination due to Galactic foreground stars.

We ended up with a sample of 184 bona-fide unbarred galaxies. We retrieved their 2MASS band images from the NASA/IPAC Infrared Science Archive. The galaxy images were reduced and flux calibrated with the standard 2MASS extended source processor GALWORKS (Jarrett00, ). Images have a typical field of view of few arc-minutes and a spatial scale of pixel. They were obtained with an average seeing 3 as measured by fitting a circular two-dimensional Gaussian to the field stars.

After a visual inspection of the images, we realized that some of the sample galaxies were not suitable for our study. We rejected paired and interacting objects as well as those galaxies that resulted in being barred after performing the photometric decomposition (see Sect. 3). Therefore, the final sample presented in this paper contains 148 galaxies (90 lenticular and 58 early-type spiral galaxies). A compilation of their main properties is given in Table 3. Figure 1 shows the distribution of the sample galaxies over the Hubble types. The lenticular galaxies are predominant over the spirals due to our magnitude selection, which favors red galaxies. Moreover, we show the distribution of radial velocities of the sample galaxies with respect to the 3K background. The mean radial velocity is 2000 km s (corresponding to a distance of 27 Mpc by assuming km s Mpc), but we include galaxies as far as 8500 km s (113 Mpc) because the sample is magnitude limited.

Tonry et al. (2001) derived the distance of 30 galaxies of our sample from the measurement of their surface brightness fluctuations. The difference between the distances obtained from radial velocities and those derived from surface brightness fluctuations was calculated for all these galaxies. The standard deviation of the distance differences is 5 Mpc and it was assumed as being a typical distance error. For the 4 common galaxies in the Virgo cluster it is 2 Mpc.

## 3 Two-dimensional bulge-disk parametric decomposition

Conventional bulge-disk decompositions based on elliptically averaged surface-brightness profiles usually do not take into account the intrinsic shapes (e.g., prieto01, ) or the position angle (e.g., trujillo01c, ) of the bulge and disk components, which can produce systematic errors in the results (e.g., byun95, ).

For this reason a number of two-dimensional parametric decomposition techniques have been developed in the last years. As examples we may point out the algorithms developed by Simard (1998, GIM2D), Peng et al. (2002, GALFIT), de Souza et al. (2004, BUDDA) and Pignatelli et al. (2006, GASPHOT). These methods were developed to solve different problems of galaxy decomposition when fitting the two-dimensional galaxy surface-brightness distribution. They use different functions to parametrize the galaxy component and different minimizations routines to perform the fit.

In this paper we present our new decomposition algorithm named GASP2D (GAlaxy Surface Photometry 2 Dimensional Decomposition). The code works like GIM2D and GASPHOT in minimizing the interaction with the user. It works in an automatical way to be more efficient when dealing with a large amount of galaxies. However, like GALFIT and BUDDA it also adopts a Leverberg-Marquard algorithm to fit the two-dimensional surface-brightness distribution of the galaxy. This reduces the amount of computational time needed to obtain a robust and reliable estimate of the galaxy structural parameters.

In the present work, we show the first version of the code. We assume that the galaxy can be modeled with only two components, the bulge and the disk. In a forthcoming paper, we will show an improved version of GASP2D with the possibility to fit other galaxy components, like bars.

### 3.1 Photometric model

We assumed the galaxy surface-brightness distribution to be the sum of the contribution of a bulge and disk component. Both of them are characterized by elliptical and concentric isophotes with constant (but possibly different) ellipticity and position angles.

Let , , be the Cartesian coordinates with the origin in the galaxy center, the axis parallel to the direction of right ascension and pointing westward, the axis parallel to the direction of declination and pointing northward, and the axis along the line-of-sight and pointing toward the observer. The plane of the sky is confined to the plane.

We adopted the Sérsic law (sersic68, ) to describe the surface brightness of the bulge component. The Sérsic law has been extensively used in the literature for modeling the surface brightness profiles of galaxies. For instance, it has been used to model the surface brightness of elliptical galaxies (grahamguzman03, ), bulges of early and late type galaxies (andredakis95, ; prieto01, ; aguerri04, ; mollenhoff04, ), the low surface brightness host galaxy of blue compact galaxies (caon05, ; amorin07, ) and dwarf elliptical galaxies (binggeli98, ; aguerri05b, ; grahamguzman03, ). It is given by

 Ib(ξ,η)=Ie10−bn[(rb/re)1/n−1], (1)

where , , and are the effective (or half-light) radius, the surface brightness at , and a shape parameter describing the curvature of the surface-brightness profile, respectively. The value of is coupled to so that half of the total luminosity of the bulge is within and can be approximated as (caon93, ). Bulge isophotes are ellipses centered on with constant position angle PA and constant ellipticity . The radius is given by

 rb = [(−(ξ−ξ0)sinPAb+(η−η0)cosPAb)2+ (2) −((ξ−ξ0)cosPAb+(η−η0)sinPAb)2/q2b]1/2.

We adopted the exponential law (freeman70, ) to describe the surface brightness of the disk component

 Id(ξ,η)=I0e−rd/h, (3)

where and are the central surface brightness and scalelength of the disk, respectively. Disk isophotes are ellipses centered on with constant position angle PA and constant ellipticity . Disk inclination is . The radius is given by

To derive the coordinates of the galaxy center and the photometric parameters of the bulge (, , , PA, and ) and disk (, , PA, and ) we fitted iteratively a model of the surface brightness =+ to the observations using a non-linear least-squares minimization method. It was based on the robust Levenberg-Marquardt method (e.g., press96, ) implemented by More et al. (1980). The actual computation was done using the MPFIT algorithm implemented by C. B. Markwardt under the IDL environment111The updated version of this code is available on http://cow.physics.wisc.edu/ craigm/idl/idl.html. MPFIT allows the user to keep constant or impose boundary constraints on any parameter during the fitting process.

For each pixel , the observed galaxy photon counts are compared with those predicted from the model . Each pixel is weighted according to the variance of its total observed photon counts due to the contribution of both galaxy and sky, and determined assuming photon noise limitation by taking into account the detector readout noise (RON). Therefore, the to be minimized can be written as

 χ2=N∑ξ=1M∑η=1[Im(ξ,η)−Ig(ξ,η)]2Ig(ξ,η)+Is(ξ,η)+RON2, (5)

with and ranging over the whole pixel image.

An important point to consider here is the weight function used to calculate the ; some authors claim that is better to assign to each pixel an uncertainty given by the Poissonian noise (e.g., peng02, ) while others adopt constant weights to obtain better results (e.g., MH01; de Souza et al. 2004). Both possibilities were implemented in the fitting algorithm. We adopted the poissonian weights after extensive testing with artificial galaxies.

### 3.2 Seeing effects

The ground-based images are affected by seeing, which scatters the light of the objects and produces a loss of spatial resolution. This is particularly critical in the central regions of galaxies, where the slope of the radial surface brightness profile is steeper. Since the bulge contribution dominates the surface brightness distribution at small radii, seeing mostly affects bulge structural parameters. Seeing effects on the scale parameters of a Sérsic surface brightness profile have been extensively discussed by Trujillo et al. (2001a,b).

During each iteration of the fitting algorithm, the seeing effects were taken into account by convolving the model image with a circular two-dimensional Gaussian point spread function (PSF). The PSF FWHM was chosen to match the one measured from the foreground stars in the field of the 2MASS galaxy image. The convolution was performed in the Fourier domain using the fast Fourier transform (FFT) algorithm (press96, ) before the calculation. Our code allows us to introduce also a Moffat or a star image to reproduce the PSF. We tested that the results are not improved by adopting a circular two-dimensional Moffat PSF or by computing the convolution integrals.

### 3.3 Technical procedure of the fit

Since the fitting algorithm is based on a minimization, it is important to adopt initial trials for free parameters as close as possible to their actual values. This would ensure that the iteration procedure does not just stop on a local minimum of the distribution. To this aim we proceeded through different steps.

First, the photometric package SExtractor (see bertinarnouts96, , for details) was used to measure position, magnitude and ellipticity of the sources (e.g., foreground stars, background and companion galaxies, as well as bad pixels) in the images.

We then derived the elliptically-averaged radial profiles of the surface brightness, ellipticity and position angle of the galaxy. We fitted ellipses to the galaxy isophotes with the ELLIPSE task in IRAF222IRAF is distributed by NOAO, which is operated by AURA Inc., under contract with the National Science Foundation.. After masking the spurious sources using the parameters provided by SExtractor, we fitted ellipses centered on the position of the galaxy center. This could be estimated by either a visual inspection of the image or with SExtractor. The coordinates of the galaxy center were adopted as initial trials for in the two-dimensional fit.

In the third step we derived some of the trial values by performing a standard one dimensional decomposition technique similar to that adopted by several authors (e.g., kormendy77, ; prieto01, ). We began by fitting an exponential law to the radial surface-brightness profile at large radii, where the light distribution of the galaxy is expected to be dominated by the disk contribution. The central surface brightness and scalelength of the fitted exponential profile were adopted as initial trials for and , respectively. The fitted profile was extrapolated at small radii and then subtracted from the observed radial surface-brightness profile. The residual radial surface-brightness profile was assumed to be a first estimate of the light distribution of the bulge. We fitted it with a Sérsic law by assuming the bulge shape parameter to be . The bulge effective radius, effective surface brightness, and shape parameter that (together with the disk parameters) gave the best fit to the radial surface-brightness profile were adopted as initial trials for , , and , respectively.

We also obtained the initial trials for ellipticity and position angle of the disk and bulge, respectively. The values for and PA were found by averaging the values in the outermost portion of the radial profiles of ellipticity and position angle. The initial trials for and PA were estimated by interpolating at the radial profiles of the ellipticity and position angle, respectively.

Finally, the initial guesses were adopted to initialize the non-linear least-squares fit to galaxy image, where all the parameters, including , were allowed to vary. A convergent model was reached when the had a minimum and the relative change of the between the iterations was less than . A model of the galaxy surface-brightness distribution was built using the fitted parameters. It was convolved with the adopted circular two-dimensional Gaussian PSF and subtracted from the observed image to obtain a residual image. In order to ensure the minimum in the -space found in this first iteration, we perform two more iterations. In these iterations, all the pixels and/or regions of the residual image with values greater or less than a fixed threshold, controlled by the user, were rejected. Those regions were masked out and the fit was repeated assuming, as initial trials for the free parameters, the values obtained in the previous iteration. These kind of masks are usually useful when galaxies have other prominent structures, different from the bulge and disk (e.g. spiral arms and dust lanes), which can affect the fitted parameters, solving the problem in an automatic way. We found that after three iterations the algorithm converges and the parameters do not change.

### 3.4 Test on simulated galaxies

To test the reliability and accuracy of our two-dimensional technique for bulge-disk decomposition, we carried out extensive simulations on a large set of artificial disk galaxies. We generated 1000 images of galaxies with a Sérsic bulge and an exponential disk. The central surface-brightness, scalelength, and apparent axial ratios of the bulge and disk of the artificial galaxies were randomly chosen in the range of values observed in the band by MH01 for a sample of 40 bright spiral galaxies. The adopted ranges were

 1≤re≤3 kpc0.4≤qb≤0.90.5≤n≤6 (6)

for the bulge parameters, and

 2≤h≤5 kpc0.4≤qd≤0.9 (7)

for the disk parameters. The parameters of the artificial galaxies also have to satisfy the following conditions

 qd

All the simulated galaxies were assumed to be at a distance of Mpc. This is the average distance of the galaxies of our sample and it corresponds to a scale of pc arcsec. The pixel scale used was 1 arcsec pixel, and the CCD gain and RON were 8 e ADU and 40 e, respectively, in order to mimic the instrumental setup of the 2MASS data. Finally, a background level and photon noise were added to the resulting images to yield a signal-to-noise ratio similar to that of the available 2MASS images.

The two-dimensional parametric decomposition was applied to analyze the images of the artificial galaxies as if they were real. Errors on the fitted parameters were estimated by comparing the input and output values. Relative errors () were assumed to be normally distributed, with mean and standard deviation corresponding to the systematic and typical error on the relevant parameter, respectively.

The results of the simulations are given in Table 1. In Run 1 we built the artificial galaxies by assuming the correct values of PSF FWHM and sky level, so only errors due to the Poisson noise are studied. The mean relative errors on the fitted parameters are smaller than (absolute value) and their standard deviations are smaller than (absolute value) for all galaxies with mag, proving the reliability of our derived structural parameters. Relative errors increase for fainter galaxies, which are not included in our sample.

Systematic errors given by a wrong estimation of PSF and sky level are the most significant contributors to the error budget. To understand how a typical error in the measurement of the PSF FWHM affects our results, we analyzed the artificial galaxies by adopting the correct sky level and a PSF FWHM that was larger (Run 2) or smaller (Run 3) than the actual one. As expected (Sect. 3.2), the parameters of the surface-brightness profiles show larger errors for the bulge than for the disk. We recovered larger values for the Sérsic parameters (r, ) when the PSF FWHM is overestimated, and lower values when it is underestimated. Relative errors are correlated with the values of effective radius and shape parameter of the bulge but not with the magnitude of the galaxy. In Run 4 we built the artificial galaxies by adopting the correct PSF FWHM and a sky level that was larger than the actual one. For brighter galaxies an improper sky subtraction mostly affects the parameters of the disk surface-brightness profile. For fainter galaxies, the large relative errors on the bulge parameters are due to their coupling with the disk parameters. This is consistent with the results of similar tests performed by Byun & Freeman (1995).

The structural parameters to be measured to derive the intrinsic shape of bulges are the ellipticity of the bulge and position angles of the bulge and disk (Sect. 5.2). In all the runs, the relative errors are smaller than (absolute value) for galaxies with mag. Larger errors ( up to about 0.1) were found for fainter galaxies after an improper subtraction of the sky level from the image. However, this is not the case for our sample galaxies.

### 3.5 Results and comparison with previous studies

The parameters derived for the structural components of the sample galaxies are collected in Table 3. All the listed values are corrected for seeing smearing and galaxy inclination. Surface brightnesses were calibrated by adopting, for the 2MASS images, the flux zero point given in the image headers (Jarrett00, ).

The comparison of the structural parameters obtained for the same galaxy by different authors is often not straightforward on account of possible differences in the observed bandpass, parameterization of the surface-brightness distribution, and fitting method.

MH01 already studied 11 of our sample galaxies (NGC 772, NGC 2775, NGC 2841, NGC 2985, NGC 3169, NGC 3626, NGC 3675, NGC 3898, NGC 4450, NGC 4501 and NGC 4826). They performed a two-dimensional parametric decomposition of the band surface brightness distribution. They considered ellipticities and position angles of both the bulge and disk as free parameters. Therefore, we considered their results as the most suitable to be compared with ours. The structural parameters we measured are consistent with those given by MH01, within for all the common galaxies but NGC 4826. We argue that they strongly underestimated the scale length of its disk (and consequently obtained a wrong estimate of the other parameters) because of the small field of view () of their image. In Fig. 2 we show the comparison between our axial ratios and position angles of the bulge and disk and those measured by MH01.

## 4 Correlations between structural parameters

The study of correlations between the structural parameters of bulges and disks of our sample galaxies will help us to both cross check our results with those available in literature and identify and rule out peculiar bulges from any further analysis.

### 4.1 Bulge parameters

We did not find any correlations between the bulge parameters and Hubble type. Neither the effective radius (Fig. 3A), effective surface brightness (Fig. 3B) nor the shape parameter (Fig. 3C) show a statistically significant Pearson correlation coefficient ().

From near-infrared observations of spiral galaxies, Andredakis et al. (1995) found that bulges of early-type spirals are characterized by (i.e., they have a de Vaucouleurs radial surface brightness profile), while the bulges of late-type spirals are characterized by (i.e., they have an exponential radial surface brightness profile). This early result was confirmed in various studies (e.g., de Jong et al. 1996; Khosroshahi et al. 2000; MacArthur et al. 2003; Möllenhoff & Heidt 2001; Möllenhoff 2004; Hunt et al. 2004). We argue that our data does not show such a correlation due to the smaller range of Hubble types covered by our sample (S0–Sb) with respect to the cited works, where it is mostly evident for Hubble types later than Sb.

The shape parameter increases with effective radius. Larger bulges have a surface-brightness radial profile that which is more centrally peaked than that of the smaller bulges (Fig. 3D). We obtained

 logn=0.38(±0.02)+0.18(±0.05)logre (r=0.28). (9)

The effective surface brightness is dependent on the effective radius. Larger bulges have a lower effective surface brightness (Fig. 3E). We found a linear regression

 logμe=17.74(±0.07)+1.7(±0.2)logre (r=0.55). (10)

This is in agreement, within the errors, with the correlation by MH01. If we use the mean surface brightness inside one effective radius instead of the effective surface brightness this relation becomes the so-called Kormendy relation, already known for bulges and elliptical galaxies (kormendy77, ).

Finally, the absolute luminosity of the bulge is correlated with the effective radius. Larger bulges are more luminous (Fig. 3F). This result in

 Mb=−21.93(±0.06)−3.4(±0.2)logre (r=−0.80), (11)

where is the band magnitude of the bulge. A similar result was obtained by MH01 for a smaller sample of disk galaxies spanning a larger range of Hubble types.

Bulges and elliptical galaxies follow a tight relation, the FP, defined by the effective radius, mean surface brightness within effective radius, and central velocity dispersion (djorgovskidavies87, ; dressler87, ). Therefore, we derived the FP for the bulges of our sample galaxies. The measurements of the central stellar velocity dispersion for a subsample of 98 galaxies were available in literature and were retrieved from the on line HyperLeda catalog (paturel03, ). The velocity dispersions given by the catalogue are corrected to a circular aperture of radius of 0.595 h kpc, which is equivalent to an angular diameter of at the distance of Coma, following the prescription by Jorgensen et al. (1995). The aperture-corrected velocity dispersions are given in Table 3. The coefficients describing the FP

 logre=1.08(±0.09)logσ0+0.25(±0.02)⟨μe⟩−6.61(±0.40), (12)

were derived by minimizing the square root of the residuals along the axis. Errors given for every coefficient were calculated by performing a bootstrap analysis with 1000 iterations. No statistically significant difference was observed when only bulges of lenticular or early-to-intermediate spiral galaxies were considered. The dispersion around this relation is dex and was measured as the rms scatter in the residuals of . The observational error on the FP is 0.066 dex and includes the measurement errors in (0.055 dex), (0.029 dex), and (0.021 mag arcsec). Errors in and are not independent (Kormendy 1977). Compared with the dispersion around the relation, this gives an intrinsic scatter of 0.088 dex. Figure 4 shows an edge-on view of the FP. Our coefficients and those by Falcón-Barroso et al. (2002) are consistent within the errors, although they analyzed band data. Unfortunately, we have not found in the literature the coefficient of the FP in the band for a direct comparison (see bernardi03, ).

One of the projections of the FP is the so-called Faber-Jackson relation (FJ), which relates the luminosity of elliptical galaxies and bulges to their central velocity dispersion (faberjackson76, ). We derived the band FJ relation for the bulge subsample obtaining

 logσ0=0.1(±0.2)−0.095(±0.009)Mb (r=−0.71). (13)

This result also holds when we consider only galaxies with errors on the central velocity dispersion lower than 10 km s ( dex). In fact, we derived

 logσ0=0.4(±0.3)−0.085(±0.013)Mb (r=−0.65), (14)

which is consistent within errors with Eq. 13. From Eq. 13 we derived , which is very close to the virial relation and indicates that our bulges share important characteristics with bright elliptical galaxies matkovicguzman05 (). On the other hand, Balcells et al. (2007) found (close to faint ellipticals) observing a sample of bulges with the Hubble Space Telescope in the band. This discrepancy is not due to the adopted fitting method and it is not observed if only bright bulges ( mag) are considered. Indeed, we found for our sample and for the sample by Balcells et al. (2007). The different behaviour of faint bulges ( mag) requires further investigation to be explained.

Khosroshahi et al. (2000) noticed that the shape parameter, effective radius and central surface brightness of elliptical galaxies and bulges are correlated. This relation was termed photometric plane (PP). Figure 6 shows an edge-on view of the PP of our bulge sample

 logn=0.17(±0.02)logre−0.088(±0.004)μ0+1.48(±0.05). (15)

The coefficients were derived by minimizing the square root of the residuals along the axis. Errors given for every coefficient were calculated by performing a bootstrap analysis with 1000 iterations. No statistically significant difference was observed when only bulges of lenticular or early-to-intermediate spiral galaxies were considered. The dispersion around this relation is and was measured as the rms scatter in the residuals of . Our coefficients and those by MH01 are consistent within the errors, although their sample is dominated by bulges of late-type spirals. The presence of the bulges of lenticular and spirals on the same PP hints towards a common formation scenario.

### 4.2 Disk parameters

Regarding the disk parameters, we found no correlation between the scale length and Hubble type (, Fig. 7A). The same is true for the central surface brightness. In fact, it shows a large scatter also with Hubble type (, Fig. 7B). This is consistent with the results of de Jong et al. (1996) and MH01.

On the other hand, the central surface brightness and the luminosity of the disks are dependent on the scale length. Larger disks have a lower central surface brightness (Fig. 7D). We found a linear regression

 logμ0=17.36(±0.1)+1.4(±0.2)logh (r=0.49), (16)

and brighter disks show larger scale lengths (Fig. 7C)

 Md=−21.21(±0.09)−3.5(±0.2)logh (r=−0.80), (17)

where is the band magnitude of the disk.

The coefficients are in agreement within the errors with those given by MH01.

### 4.3 Bulge and disk interplay

We have found that the disk scale length increases with central velocity dispersion. Since central velocity dispersion correlates with the virial mass of the bulge ( with , Cappellari et al. 2006), we conclude that larger disks are located in galaxies with more massive bulges (Fig. 8A). For the subsample of 98 early-to-intermediate bulges with a measured velocity dispersion we found

 logσ0=2.13(±0.03)+0.27(±0.06)logh (r=0.42). (18)

We also found a strong correlation between the bulge effective radius and the disk scale length. Larger bulges reside in larger disks (Fig. 8B). This relation was already observed by Courteau et al. (1996) and later observed in NIR by MH01, MacArthur et al. (2003). We obtained a linear regression

 logre=−0.45(±0.03)+0.91(±0.07)logh (r=0.74), (19)

which is in agreement within error bars with the correlation found by MH01. All these correlations between bulge and disk parameters indicate a link between the bulge and disk formation and evolution history. This connection was interpreted as an indication of the formation of late-type bulges via secular evolution of the disks (courteau96, ). However, our measurements of the scale lengths of the bulge and disk and (Fig. 8C) are also fully consistent with the predictions of the numerical simulations by Scannapieco & Tissera (2003) and Tissera et al. (2006). They studied the effects of mergers on the mass distribution of bulges and disks of galaxies formed in hierarchical clustering scenarios. Our mean value is also in good agreement with found by Naab & Trujillo (2006) for a series of major mergers’ remnants. These results indicate that these relations are not enough to distinguish between bulges formed by mergers or by secular evolution of the disk, even if a strong crosstalk between both components is present.

## 5 The equatorial intrinsic ellipticity of bulges

In the previous section, we realized that the knowledge of correlations between structural parameters of the bulge and disk is not sufficient to distinguish between the different scenarios that were proposed to explain the formation of bulges. Therefore, we decided to study the intrinsic shape of bulges in order to give a further constraint to these scenarios.

Independently of its internal structure, we can consider a bulge of a spiral galaxy as an ellipsoidal stellar system located in the center of the galaxy,which stands out against the disk in the photometric observations. We assume that both the bulge and disk share the same center, which coincides with the galactic one, and they have the same polar axis (i.e., the equatorial plane of disk coincides with that of bulge).

### 5.1 Geometrical formalism

Let be Cartesian coordinates with the origin in the galaxy center, the axis and axis corresponding to the principal axes of the bulge equatorial ellipse, and the axis corresponding to the polar axis. As the equatorial plane of the bulge coincides with the equatorial plane of the disk, the axis is also the polar axis of the disk. If , , and are lengths of the ellipsoid semi-axes, the corresponding equation of the bulge in its own reference system is given by

 x2A2+y2B2+z2C2=1. (20)

Let be now the Cartesian coordinates of the observer system. It has its origin in the galaxy center, the polar axis along the line-of-sight (LOS) and pointing toward the galaxy, and the plane of the sky lies on the plane.

The projection of the disk onto the sky plane is an ellipse whose major axis is the line of nodes (LON), i.e., the intersection between the galactic and the sky planes. The angle between the axis and axis corresponds to the inclination of the bulge ellipsoid; it can be derived as from the length and of the two semi-axes of the projected ellipse of the disk. We defined () as the angle between the -axis and the LON on the equatorial plane of the bulge . Finally, we also defined () as the angle between the -axis and the LON on the sky plane . The three angles , , and are the usual Euler angles and relate the reference system of the ellipsoid with that of the observer by means of three rotations. Indeed, since the location of the LON is known, we can choose the axis along it, and consequently it holds that . By applying these two rotations to Eq. 20 it is possible to derive the equation of the ellipsoidal bulge in the reference system of the observer, as well as the equation of the ellipse corresponding to its projection on the sky plane (simonneau98, ). Now, if we identify this ellipse with the ellipse that forms the observed ellipsoidal bulge, we can determine the corresponding axes of symmetry x and y. The first one, on which we measured the semi-axis , forms an angle with the LON (the x’-axis in the observed plane); the semi-axis is taken over the y-axis. We always choose , so it is possible that either be the major or the minor semi-axis, and vice versa for .

We have

 a2b2 = A2C2sin2θcos2ϕ+B2C2sin2θsin2ϕ+A2B2cos2θ. (21) a2+b2 = A2(cos2ϕ+cos2θsin2ϕ)+B2(sin2ϕ+cos2θcos2ϕ)+C2sin2θ. (22) tan2δ = (B2−A2)cosθsin2ϕA2(cos2θsin2ϕ−cos2ϕ)+B2(cos2θcos2ϕ−sin2)+C2sin2θ. (23)

If the ellipsoidal bulge is triaxial () then it is possible to observe a twist (; see Eq. 23) between the axes of the projected ellipses of the bulge and the disk.

### 5.2 Inverse problem or deprojection

We will now focus our attention on the inverse problem, i.e., deprojection. Following Simonneau et al. (1998), from Eqs. 21, 22, and 23, we are able to express the length of the bulge semi-axes (i.e. , , and ) as a function of the length of the semi-axes (i.e. , ) of the projected ellipse and the position angle ().

 A2 = a2+b22[1+e(cos2δ+sin2δsinϕcosϕ1cosθ)] (24) B2 = a2+b22[1+e(cos2δ−sin2δcosϕsinϕ1cosθ)] (25) C2 = a2+b22[1+e(2sin2δcosθsin2θcos2ϕsin2ϕ−cos2δ1+cos2θsin2θ)], (26)

where is, in some way, a measure of the ellipticity. It will be .

Notice that and are all observed variables. Unfortunately, the relation between the intrinsic and projected variables also depends on the spatial position of the bulge (i.e., on the angle), which is not directly accessible to observations. For this reason, only a statistical determination can be performed to assess the intrinsic shape of bulges.

As and are the semi-axis of the equatorial ellipse of the bulge, we have to distinguish between two cases, according to Eqs. 24 and 25. If (or equivalently ) then . Otherwise, if (or equivalently ) then . Thus, if the equatorial plane of the bulge ellipsoid is not circular and the bulge ellipsoid is triaxial.

From Eqs. 24, 25 and 26 we can write the axial ratios and as explicit functions of . Moreover, we assume that the angle is random and independent of the length of ellipsoid semi-axes. Thus, the normalized probability distribution of getting a given value of in + is

 P(ϕ)=2/π  ;∫π/20P(ϕ)dϕ=1. (27)

According to a fundamental theorem of statistics, the probability of obtaining a given value of any function (e.g., one of the axial ratios) will be equal to the probability of getting the corresponding value of , provided that the ratio between the corresponding differential elements is taken into account.

In this work, we will focus our attention on the intrinsic equatorial ellipticity of bulges. We define it as -)/(+ with . In a forthcoming paper, we will study the intrinsic flattening of the bulge ellipsoids defined as the ratio between the length of polar semi-axis and the mean length of the equatorial semi-axes.

From Eqs. 24 and 25, it is straightforward to derive a relation among the intrinsic variables (equatorial ellipticity and position angle ), and the measured ( i.e. , , and ), which is

 Esin(2ϕ)1+Ecos(2ϕ)=1cosθesin2δ1+ecos2δ≡Q. (28)

The second member of the equality in Eq. 28 allows us to define the observable in terms of the measured variables , , and . It must be stressed that for each specific bulge the relation between the equatorial ellipticity and the unknown parameter embraces the whole of the measured variables through the single variable .

On the one hand, Eq. 28 will yield the conditional probability that a given bulge, with a measured value of , takes on any particular value of (individual statistic); on the other hand, this equation will give the probability associated to each value of for a bulge with intrinsic equatorial ellipticity . This latter probability will be the kernel of an integral equation that relates the observed statistical distribution , corresponding to a sample of galaxies, with the statistical distribution of the equatorial ellipticity for the same sample.

### 5.3 Individual statistics

For a given galaxy, we can measure the values of , and , and then derive the value of through Eq. 28. We want to determine the probability that such a galaxy (i.e. with such a value of ) will take on a value of in the range +. The subindex specifies this galaxy. All the galaxies with the same value of shall partake the same probability distribution .

Once the value of is prescribed, for some values of there are not values of () that satisfy Eq. 28. Hence, it shall hold that . Only for those values of such that

 E2≥Q21+Q2≡T2  ;  0≤T≤1, (29)

will there will exist two values of that fulfill Eq. 28.

Then for any value of the probability will be given by

 PQ(E)=∑j=1,2(P(ϕ)∣∣∣δϕδE∣∣∣Q)ϕj . (30)

By calculating the partial derivative of Eq. 30 and normalizing we obtain

 PT(E)=TE1√E2−T2arccosT , (31)

where the subindex plays the same role as the subindex ; both are related through Eq. 29.

This means that the possible values of are very concentrated and slightly larger than . To get an idea of how is peaked near the value of , we calculated the value for which the total probability that is equal to the probability that . For every bulge is a sort of mean value of , and is given by

 E1/2=√21+TT. (32)

In Fig.9 we show, as an example, the probability for one galaxy of our sample.

### 5.4 Global statistics

Likewise, we can define the probability associated to each value of for a given bulge of intrinsic equatorial ellipticity .

For a prescribed value of , it will only be possible to get the values of that satisfy Eq. 28 when . The probability for the two values of is given by Eq. 27. The probability is equal to the sum of the two probabilities , weighted with the ratio of the differential elements:

 PE(Q)=∑j=1,2(P(ϕ)∣∣∣δϕδQ∣∣∣E)ϕj. (33)

Once the partial derivative from Eq. 28 is computed, we obtain

 PE(T)=2π1√E2−T2  ; E≥T . (34)

Our purpose here is to determine the probability distribution , that is the number of galaxies with a value of inside the interval , starting from a sample that is sufficiently representative. We have measured for such a sample the distribution , namely the number of bulges whose values of are within and . We must write now as the integral over all the values of of the product of the conditional probability by the so far unknown distribution :

 P(T)=∫1−1P(E)PE(T)dE. (35)

Then, by making use of Eq. 34, we obtain an Abel-like integral equation

 P(T)=2π∫1TP(E)√E2−T2dE, (36)

which will allow us to derive from the observed distribution .

However, as usual, the data of our statistical problem takes the form of histograms, hence, the relevant equations must be formulated accordingly.

Let with (, ) be a set of discrete ordinates that defines the histogram of the observed function (Fig. 10). The -element of this histogram is defined by

 Pk(T)=1Tk−Tk−1∫TkTk−1P(T)dT. (37)

We must now seek the integral equation that relates the variables with the probability distribution .

We notice that the integral of in Eq. 37 is equivalent to the difference between the two quadratures of over the intervals (,1) and (,1). In both of them, we will replace by its integral form, given by Eq. 36.

Then, by inverting the order of integration, we can easily rewrite the two resulting integrals to obtain

 Pk(T)=1Tk−Tk−1∫TkTk−1P(E)dE
 −2/πTk−Tk−1(M(Tk−1)−M(Tk)), (38)

where we defined

 M(Tk)=∫1TkP(E)arcsin(Tk/E)dE, (39)

and in a similar way.

Equation 38 is the integral equation that will allow us to obtain the values of . Since it is consistent with the numerical structure of the data, we are confident that we have eliminated most of the numerical problems that arise from the direct inversion of Eq. 36, which constitutes a typical ill-posed problem.

At this point, we require for the same histogram representation that we have already introduced for the data . We introduce a similar set of discrete ordinates for the variable (i.e. ), and an analogous definition to obtain the elements of the histogram

 Pj(E)=1Ej−Ej−1∫EjEj−1P(E)dE. (40)

Thus Eq. 38 can be rewritten into the form

 Pk(T)=Pk(E)−2/πTk−Tk−1(M(Tk−1)−M(Tk)). (41)

In order to express the integrals and as linear functions of the so far unknown values of (where ), we consider that is constant and equal to the unknown values of over each interval (, ), according to its histogram representation. Therefore, it is straightforward to compute the coefficients , defined by the linear relation

 M(Tk)=N∑j=k+1CM(k,j)Pj(E). (42)

Thus Eq. 41 becomes a simple linear algebraical equation that relates the terms of the two histograms and through a triangular matrix. Once we have the integral Eq. 36 into a suitable matrix form, according to the histogram representation of the data and results, a simple matrix inversion could be, in principle, enough to obtain the resulting . However, such a procedure may add to the intrinsic difficulties, due to the lack of precision typical of the observational data, which naturally arise in the matrix inversion process; a catastrophic mixture when dealing with an inverse problem.

### 5.5 Inversion methods for the integral equation

We have considered two different approaches to tackle the numerical problem. The first one is suggested by the method that leads us to the integral Eq. 41 for the set of discrete elements . We notice that there are two different terms in its right-hand side. The first one is the identity operator. The second one is the difference between two integrals of , multiplied by a kernel that is quickly decreasing. We can consider the former as the leading term, and treat the latter as a corrective term in a iterative perturbation method. Consequently, we can write Eq. 41 as

 Pk(E)=Pk(T)+2/πTk−Tk−1(M(Tk−1)−M(Tk)), (43)

where we can determine P at the -iteration making use of the form of P from the -iteration to compute the correction term . As an initial guess, we consider equal to at a zero-order approximation. This iterative process is repeated until convergence is achieved.

From the data of the histogram shown in Fig. 10, we have obtained a satisfactory solution with a few number (5-10) of iterations. We will discuss later the physical quality of this solution, but we must recognize here the stability of the method. Actually, we have always achieved the same solution with a few number of iterations, starting from any trial initial distribution (namely, any zero-order approximation for ). Moreover, the greater advantage of solving the integral equation by means of an iterative perturbation method is the possibility to recover, according to Eq. 41, the approximate diagram that corresponds to any iterative solution , and this yields a double check on the evolution and quality of results.

However, when dealing with this kind of inverse problem, it may often happen in the practice that the results are correct mathematically, but not from the physical standpoint. In view of this difficulty, and in spite of the excellent quality of the foregoing iterative method, we wished to develop an alternative method of inversion for the integral Eq. 36, in order to double check the results.

The other way to numerically treat the integral Eq. 36 comes from the analytical inversion

 ∫1T′P(E)=∫1T′T√T2−T′2P(T)dT′. (44)

Once we have this analytical form for the required solution, we rewrite it for the histogram representation of , by defining as in Eq. 40. Again the difference between the two integrals of show up. Now, we impose the histogram model of this to derive analytically the matrix elements that relate any to all the elements . The corresponding solutions obtained with the two methods are the same, allowing for the small differences due to round-off errors of the two different numerical algorithms employed.

Now that we are confident of the reliability of both methods of inversion of the integral equation, we can came back to the aforesaid difficulties.

When applying either method to the observed distribution , in form of a histogram with bins as shown in Fig. 10, we may obtain non-physical results, these are negative values for some bins . This occurs due to the fact that in the frame of the adopted histogram representation for and and the associated numerical algorithm chosen to represent the matrix operator for the integral Eq. 36, the measured values cannot be the integral transform of any physical distribution . However, another set of values that are slightly different from the original, and consequently compatible with the observations, might satisfy the above requirement; i.e., it can be the integral transform of some physical and its inverse transform will solve our problem.

These considerations claim a statistical regularization process, which can be achieved by considering the histogram to be the statistical mean of 1000 histograms, all of them compatible with the observations according to Poisson statistics. For each one of the 1000 possible realizations for , we have obtained the corresponding by means of the inversion of the integral equation. The non-physical histograms , i.e., those with some negative bins, were rejected. From the physical solutions, we have obtained the mean histogram and the corresponding error bars, as shown in Fig. 11. The statistical regularization process also allows us to estimate errors due to the possible lack of statistics in the sample.

### 5.6 The probability distribution function of intrinsic ellipticities

In Fig. 11 we present the PDF of the bulge intrinsic ellipticities . It was obtained by applying the procedure described in Sect. 5.5 using the PDF shown in Fig. 10. The values for each galaxy were calculated by means of Eqs. 28 and 29 from the measured values of , and .

The PDF is characterized by a significant decrease of probability for (or equivalently ), suggesting that the shape of bulge ellipsoids in their equatorial plane is most probably elliptical rather than circular. Such a decrease is caused neither by the lack of statistics (because in the regularization method we took into account the Poisson noise) nor by the width of the bins (because we tried different bin widths).

We have calculated the average value weighted with the PDF through

 ⟨E⟩=∑iPDF(Ei)∗Ei∑iPDF(Ei), (45)

obtaining a value of (). This is fully consistent with previous findings by Bertola et al. (1991) and Fathi & Peletier (2003), based on the analysis of smaller samples of bulges. For the sake of comparison, the value of was derived from their data using Eq. 45. It is for the bulges studied by Bertola et al. (1991). They adopted a different approach to derive the PDF for intrinsic axial ratio of bulges from the misalignment of the major axes of bulges and disks and the apparent ellipticity of bulges. It is for the early-type bulges of the sample studied by Fathi & Peletier (2003). They measured the bulge equatorial ellipticity by analyzing the deprojected ellipticity of the ellipses fitting the galaxy isophotes within the bulge radius.

A further important result derived from is that there are not bulges with (). This is also in good agreement with Bertola et al. (1991) and Fathi Peletier (2003). They found a minimum axial ratio and , respectively.

We also studied the possible differences in the shape of bulges depending on their observational characteristics (i.e., morphology, light concentration, and luminosity, see Fig. 12). First of all, we subdivided our bulges according to the morphological classification and of their host galaxies to look for differences between lenticular and early-type spiral galaxies. A second test was done by subdividing the sample bulge between those with a Sérsic index and those with to investigate possible correlations of bulge shape with light concentration. Finally, we subdivided our bulge into faint () and bright () in order to search for differences of bulge shape with the band total luminosity. We did not find any significant difference between the studied subsamples. They are characterized by the same distribution of , as confirmed at a high confidence level by a Kolmogorov-Smirnov test.

Several authors discussed the problem of the intrinsic shape of elliptical galaxies by means of observations and/or numerical simulations. Ryden (1992), Lambas et al. (1992), and Bak & Statler (2000) agree that the observed distribution of ellipticities cannot be reproduced by any distribution of either prolate or oblate spheroidal systems. Any acceptable distribution of triaxial systems is dominated by nearly-oblate spheroidal rather than nearly-prolate spheroidal systems. The formation of triaxial elliptical galaxies via simulation of merging events in the framework of a hierarchical clustering assembly was studied by Barnes Hernquist (1996), Naab Burkert (2003) and Gonzalez-Garcia Balcells (2005). On the other hand, in the monolithic scenario where the galaxy formation occurs at high redshift after a rapid collapse, we may expect that the final galaxy shape would be nearly spherical or axisymmetric, as recently found in numerical experiments by Merlin & Chiosi (2006). But there is no extensive testing of the predictions of numerical simulations against the derived distribution of bulge intrinsic ellipticities.