text:
# Detailed Structural Decomposition of Galaxy ImagesBased on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.

## 1 INTRODUCTION

Galaxies span a wide range of morphology and luminosity, and a very useful way to quantify them is to fit their light distribution with parametric functions. The de Vaucouleurs $`R^{1/4}`$ and exponential disk functions became standard functions to use after de Vaucouleurs (1948) showed many elliptical galaxies to have $`R^{1/4}`$ light distributions, while Freeman (1970) found later-type galaxies to be well described by a de Vaucouleurs bulge plus an exponential disk. Since then, the empirical techniques of galaxy fitting and decomposition have led to a number of notable advances in understanding galaxy formation and evolution. These include investigations into the Tully-Fisher relationship (Tully & Fischer 1977), the fundamental plane of spheroids (Faber et al. 1987; Dressler et al. 1987; Djorgovski & Davis 1987; Bender, Burstein, & Faber 1992), the morphological transformation of galaxies in cluster environments (e.g., Dressler 1980; van Dokkum, & Franx 2001), the bimodality of galaxy nuclear cusps (Lauer et al. 1995; Faber et al. 1997) and its implications for the formation of massive black holes (Ravindranath, Ho, & Filippenko 2002), and the cosmic evolution of galaxy morphology (e.g., Lilly et al. 1998; Marleau & Simard 1998).

There are two general types of galaxy fitting: one-dimensional (1-D) fitting of surface brightness profiles (e.g., Kormendy 1977; Burstein 1979; Boroson 1981; Kent 1985; Baggett, Baggett, & Anderson 1998), and 2-D fitting of galaxy images (e.g., Shaw & Gilmore 1989; Byun & Freeman 1995; de Jong 1996; Simard 1998; Wadadekar, Robbason, & Kembhavi 1999; Khosroshahi, Wadadekar, & Kembhavi 2000), with each its own tradeoffs and benefits.

In 1-D, an important consideration is how to first obtain a radial surface brightness profile from a 2-D image, for which there is no universally agreed upon procedure. A common practice is to use isophote fitting, which is a powerful technique when performed on well-resolved galaxies because it averages over elliptical annuli to increase the signal-to-noise ratio (S/N) at a given radius. However, as many galaxies have isophote twists and changing ellipticity as a function of radius, the galaxy profile is extracted along a radial arc that is ill-defined. An alternative approach is to use a direct 1-D slice across an image. Burstein (1979) argues that only cuts along the major axis should be used in bulge-to-disk (B/D) decompositions. Meanwhile, Ferrarese et al. (1994) point out that galaxies with power-law central profiles may have different profiles along the major and minor axis.

Fitting profiles in 1-D is frequently used because it suffices for certain goals, and is simple to implement. But many studies now resort to 2-D techniques. For B/D decompositions, a number of authors (e.g., Byun & Freeman 1995; Wadadekar et al. 1999) have used idealized simulations to show that 2-D modeling can better recover the true parameter values. In 1-D, while the galaxy bulge and disk may appear to merge smoothly, which causes non-uniqueness in the decompositions, in 2-D, isophote twists and ellipticity changes provide additional constraints to break those degeneracies.

While exponential and de Vaucouleurs functions can model a wide range of global galaxy profiles, galaxies are generally more complex. Well resolved, nearby galaxies are often poorly fit by the standard models in detail, especially in 2-D. Yet, few studies have gone beyond a global bulge and disk decomposition, except in clear-cut cases where there are nuclear bars (e.g., de Jong 1996) or nuclear point sources (e.g., Wadadekar et al. 1999; Portal 1999; Quillen et al. 2001; and many recent studies of quasar host galaxies).

More recently, a growing number of theories and Hubble Space Telescope (HST) observations show that many clues of galaxy formation lay hidden in the fine details of galaxy structure. New HST images reveal striking correlations between nuclear cusps of galaxies with their mass, stellar velocity dispersion, radius, and other large-scale properties (Jaffe et al. 1994; Lauer et al. 1995; Faber et al. 1997; Carollo et al. 1997; Rest et al. 2001; Ravindranath et al. 2001; Gebhardt et al. 2000; Ferrarese & Merritt 2000). Those observations have led to on-going debates about whether the correlations can be explained by an adiabatic growth of black holes in isothermal cores and binary black hole mergers (e.g. van der Marel 1999, and Ravindranath et al. 2002). Other fossil remnants of galaxy formation manifest as photometrically distinct nuclei, nuclear disks of stars and gas, dust lanes, and nuclear spiral patterns (e.g., Phillips et al. 1996; Carollo, Stiavelli, & Mack 1998; Tomita et al. 2000; Tran et al. 2001). If elliptical galaxies form from 3:1 and 1:1 mass mergers, then Naab & Burkert (2001) merger simulations predict that all ellipticals should contain a significant embedded stellar disk component in order to explain the kinematic line profiles. Is this seen, and if so, how often? High resolution and high S/N images now permit this to be tested. Active galactic nuclei (AGNs) are intimately related to formation of bulges and black holes; their fuel-starved remnants are common tenants in nearby galaxies (Ho, Filippenko, & Sargent 1997a). Yet, quantitative investigations of these central sources have often been hampered by the difficulties encountered in separating the faint nucleus from the surrounding bright bulge (see, e.g., Ho 1999; Ho & Peng 2001; Sarzi et al. 2001). In order to understand the physical implications of all these clues, one must be able to extract accurate, quantitative morphological information, via detailed galaxy decomposition.

Motivated by these and other possible applications, we developed a technique (GALFIT) to accurately model galaxy profiles. Section 2 discusses the mechanics of our program and its computational requirements. Section 3 provides an overview of our fitting procedures. Section 4 presents a number of case studies to illustrate the versatility of the program. We then discuss 2-D decomposition to extract nuclear sources in § 5, using simulations to evaluate the recovery of point sources in 1-D compared to 2-D. Then we apply our method to some nearby galaxies that contain weak nuclear sources . Conclusions are given in § 6. We make the source code readily available to the public upon request.

## 2 TWO-DIMENSIONAL FITTING PROGRAM

We created a 2-D image decomposition program called GALFIT, written in the C language. In order to model galaxy profiles with a maximum degree of flexibility, GALFIT uses a number of functions and can combine an arbitrary number of them simultaneously. To use the program the user provides a simple input file, as shown in Figure 1. In the example we specify a fit using five functions simultaneously: a Nuker function, a Sérsic, an exponential disk, a Gaussian, and a uniform sky pedestal (see below for details). Additional components can be added by extending the list without limit, except by the computer memory and speed available to the user.

### 2.1 Accounting for Telescope and Atmospheric Seeing

All recorded images carry an imprint of the observing set-up, due either to imperfections or diffraction in the telescope optics, and sometimes by the atmosphere. To measure the intrinsic profile of an object, the seeing must be taken into account either by convolution or deconvolution techniques. Even with HST WFPC2, whose optics are noted for near perfectly sharp, diffraction-limited, point-source images, $`10\%`$ of stellar light falls outside $`0\stackrel{}{\mathrm{.}}5`$ of its core. The point-spread function (PSF) of NICMOS on HST has even more pronounced structure; the diffraction wings are extended, and $`10\%`$ of the flux lies beyond $``$$`1\mathrm{}`$. One way to remove the seeing is to deconvolve the image. Lauer et al. (1995) show that galaxy images with high S/N at the center can be accurately deconvolved using $`4080`$ iterations of Lucy-Richardson algorithm, if an accurate PSF is available. An alternative approach uses convolution, where one constructs a model image and convolves it with a PSF before comparing the result to data. While both techniques rely on knowing an accurate PSF, deconvolution has benefits when the S/N is high — the observed image needs only to be deconvolved once and one does not have to assume a particular model for the intrinsic image. However, deconvolution will not work well for low-S/N or undersampled images because it may amplify Poisson and pixellation noise. In contrast, the convolution scheme works on images with both high or relatively low S/N. The drawback is that fitting an image proceeds more slowly because convolution is done repeatedly, which is computationally intensive. Another disadvantage is that one needs to have reasonable a priori knowledge of the input model.

GALFIT uses convolution, but it can be turned off if not needed. Convolution is done by using the convolution theorem: we multiply the Fourier transforms of the PSF and the models, and then inverse transform them. Details of the convolution and fitting processes are elaborated below.

### 2.2 Computing the $`\chi _\nu ^2`$

The four input images for GALFIT are the CCD image of the galaxy, a noise array, a PSF, and an optional dust (or bad pixel) mask — all in FITS image file format. Pixels in the dust mask are rejected from the fit.

During the fit, $`\chi _\nu ^2`$ (the reduced $`\chi ^2`$) is minimized, defined in the standard way as

$$\chi _\nu ^2=\frac{1}{N_{\mathrm{dof}}}\underset{x=1}{\overset{nx}{}}\underset{y=1}{\overset{ny}{}}\frac{\left(\text{flux}_{x,y}\text{model}_{x,y}\right)^2}{\sigma _{x,y}^{}{}_{}{}^{2}},$$
(1)
where

$$\text{model}_{x,y}=\underset{\nu =1}{\overset{nf}{}}f_{\nu ,x,y}(\alpha _1\mathrm{}\alpha _n).$$
(2)
$`N_{\mathrm{dof}}`$ is the number of degrees of freedom in the fit; $`nx`$ and $`ny`$ are the $`x`$ and $`y`$ image dimensions; and $`\text{flux}_{x,y}`$ is the image flux at pixel $`(x,y)`$. The $`\text{model}_{x,y}`$ is the sum of the $`nf`$ functions $`f_{\nu ,x,y}(\alpha _1\mathrm{}\alpha _n)`$ employed, where $`\{\alpha _1\mathrm{}\alpha _n\}`$ are the 2-D model parameters. The uncertainty as a function of pixel position, $`\sigma _{x,y}`$, is the Poisson error at each pixel, which can be provided as an input image. If no noise image is given, one is generated based on the gain and read-noise parameters contained in the image header.

### 2.3 Galaxy and Compact Source Profiles

#### 2.3.1 Azimuthal Shape and Pixel Sampling

To fit galaxy and compact source profiles we assume the profiles all are axially symmetric, generalized ellipses. The form of the generalized ellipse is motivated by Athanassoula et al. (1990) who originally use it to describe galaxy bar shapes. When the principle axes of the ellipse are aligned with the coordinate axes, the radial pixel coordinate is given by

$$r=\left(\left|x\right|^{c+2}+\left|\frac{y}{q}\right|^{c+2}\right)^{\frac{1}{c+2}}.$$
(3)
An offset of 2 in the exponents is used so that a pure ellipse has $`c=0`$, a boxy shape $`c>0`$, and a disky shape $`c<0`$. This serves exactly the same purpose as the cos $`4\theta `$ Fourier coefficient frequently associated with isophote fitting (Jedrzejewski 1987). The parameter $`c`$ differs from the cos $`4\theta `$ Fourier coefficient in that the latter applies to individual isophotes locally while the former applies to an entire component. The parameter $`q`$ is the ratio of the minor to major axis of an ellipse. Both $`q`$ and $`c`$ for each model are constant as a function of radius, although they are free parameters. Figure 2a and 2b show two examples of azimuthal shapes as a function of $`c`$ for $`q=1.0`$ and $`q=0.5`$. The PA of a component is defined with respect to the image pixel coordinate system such that a major axis positioned vertically is $`0^{}`$ (nominally North if rotated to the standard orientation) and increases counter-clockwise (nominally toward the East).

When creating a model image, GALFIT decides whether or not to oversample the pixels based on their distance from the centroid of a model component. For pixels far away from the center, the gradient is shallow; it suffices to sample only at the center of each pixel. But pixels near the center of a component must be over-sampled. For many cases dividing the pixels up by square grids then summing the flux:

$$\text{model}_{\overline{x},\overline{y}}=\underset{pixel}{}f_{model}(x\overline{x},y\overline{y})\text{ d}x\text{ d}y$$
(4)
is sufficient. However, for the Nuker profile (see below), when $`\gamma `$ is large the function changes very quickly around $`r=0`$, and the simple method quickly becomes inaccurate. We devise a technique to integrate the pixels around the center ($`r<3`$) by using an elliptical polar grid that adapts to the ellipticity of the galaxy. Figure 3a and 3b demonstrate the gridding for axis ratios of 0.9 and 0.3. The angular spacing is a 1-degree interval in a circular polar coordinate, modulo a tilt by the inclination of the component. The radial spacing of the grids increases geometrically toward the center following

$$r_n=\frac{r_{n1}}{(1+\text{astep})^n},$$
(5)
where $`n`$ increases closer to the center. Empirically, astep = 0.01 works well for even the steepest Nuker profile, and we keep it as a constant in the fit.

#### 2.3.2 Radial Profiles

The Sérsic Profile. The Sérsic (1968) profile has the following form:

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_e\text{ e}^{\kappa \left[(\frac{r}{r_e})^{1/n}1\right]}$$
(6)
where $`r_e`$ is the effective radius of the galaxy, $`\mathrm{\Sigma }_e`$ is the surface brightness at $`r_e`$, $`n`$ is the power-law index, and $`\kappa `$ is coupled to $`n`$ such that half of the total flux is always within $`r_e`$. For $`n2`$, $`\kappa 2n0.331`$; at low $`n`$, $`\kappa (n)`$ flattens out towards 0 and is obtained by interpolation. The original de Vaucouleurs (1948) profile is a special case with $`n=4`$ and $`\kappa =7.67`$. While the de Vaucouleurs profile is well suited for “classical” bulges, some bulges may be better represented by exponential profiles (e.g., Kormendy & Bruzual 1978; Shaw & Gilmore 1989; Kent, Dame, & Fazio 1991; Andrekakis & Sanders 1994). The elegance of the Sérsic profile is that it forms a continuous sequence from a Gaussian ($`n=0.5`$) to an exponential ($`n=1`$) to a de Vaucouleurs ($`n=4`$) profile simply by varying the exponent. It is very useful for modeling bars and flat disks; the smaller the index $`n`$ is, the faster the core flattens within $`r<r_e`$, and the steeper the intensity drop beyond $`r>r_e`$. The flux, integrated over all radii for an elliptical Sérsic profile with an axis ratio $`q`$ is

$$F_{\mathrm{tot}}=2\pi r_e^2\mathrm{\Sigma }_ee^\kappa n\kappa ^{2n}\mathrm{\Gamma }(2n)q/R(c),$$
(7)
where $`\mathrm{\Gamma }(2n)`$ is the Gamma function. $`F_{\mathrm{tot}}`$ is converted into a magnitude by GALFIT using the standard FITS exposure time parameter (EXPTIME) in the image header. $`R(c)`$ is a function that accounts for the area ratio between a perfect ellipse and a generalized ellipse of diskiness/boxiness parameter $`c`$, given by

$$R(c)=\frac{\pi c}{4\beta (1/c,1+1/c)},$$
(8)
where $`\beta (1/c,1+1/c)`$ is the “Beta” function with two arguments. In the 2-D implementation, the Sérsic model has eight free parameters: $`x_{cent},y_{cent},F_{\mathrm{tot}},r_e,n,c,q`$, PA. We note that in place of fitting $`\mathrm{\Sigma }_e`$, $`F_{\mathrm{tot}}`$ is fitted instead, which is more often a useful parameter. This is also done for all other models below, except for the Nuker function.

The Exponential Disk Profile. The exponential profile and the total flux are simply

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_0\text{ e}^{\frac{r}{r_s}}$$
(9)
and

$$F_{\mathrm{tot}}=2\pi r_s^2\mathrm{\Sigma }_0q/R(c),$$
(10)
where $`\mathrm{\Sigma }_0`$ is the central surface brightness and $`r_s`$ is the disk scale length. The relationship between the half-light radius, $`r_e`$, and the scale length, $`r_s`$ is $`r_s=1.678r_e`$ for this profile. Most disky galaxies are not composed of a single exponential disk, but also have either a central Sérsic or de Vaucouleurs component. They may also have either a flat core or a truncated disk, which deviates from a simple exponential (e.g., van der Kruit 1979; Pohlen, Dettmar, & Lütticke 2000).

The Nuker Law. The “modified Nuker” law was initially proposed by Lauer et al. (1995) to fit the diverse inner 1-D profiles of galaxies observed with HST. Its high degree of flexibility makes it an excellent model for fitting most 1-D galaxy profiles. The functional form of the Nuker law, which can be thought of as a double power law mediated by a smooth transition, is the following:

$$I(r)=I_b2^{\frac{\beta \gamma }{\alpha }}\left(\frac{r}{r_b}\right)^\gamma \left[1+\left(\frac{r}{r_b}\right)^\alpha \right]^{\frac{\gamma \beta }{\alpha }}$$
(11)
This function has five adjustable parameters: $`I_b,r_b,\alpha ,\beta ,`$ and $`\gamma `$. Taken to the limits of large and small radii, the parameter $`\gamma `$ is the slope of the inner power law, and $`\beta `$ is the slope of the outer power law. The break radius $`r_b`$ is the location where the profile changes slope, $`I_b`$ is the surface brightness at $`r_b`$, and $`\alpha `$ describes how sharply the two power laws connect. The more positive $`\alpha `$ is, the sharper is the break at $`r_b`$. In GALFIT, we fit the surface brightness magnitude $`\mu _b`$ instead of $`I_b`$. $`\mu _b`$ is obtained from $`I_b`$ via the use of the EXPTIME image header parameter and the pixel size specified in item K of Fig. 1. Although the profile appears singular at the center when sampled at $`r=0`$, the integrated flux is finite for $`\gamma 2`$. The 2-D Nuker profile has a total of 10 free parameters: $`x_{cent},y_{cent},\mu _b,r_b,\alpha ,\beta ,\gamma ,c,q`$, PA.

Although the parameters in Equation 1 correlate, the larger $`\alpha `$ is (i.e. the sharper the break), the less coupled are $`\alpha `$, $`\beta `$, and $`\gamma `$, when $`r_b`$ is well resolved (but much smaller than the image radius). The parameter coupling between $`\beta `$ and $`\gamma `$ is small for $`\alpha `$ fixed at a value $`1`$, which is a useful reference point. In this case, varying $`\beta `$ across its whole range of typical values ($`0<\beta <2.5`$) would affect the slope of $`\gamma `$ measured within $`r0.1r_b`$ by $`0.04`$, for any given value of $`\gamma `$.

The Gaussian Profile. The Gaussian function and its total flux are

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_0\text{ e}^{\frac{r^2}{2\sigma ^2}}$$
(12)
and

$$F_{\mathrm{tot}}=2\pi \sigma ^2\mathrm{\Sigma }_0q/R(c).$$
(13)
The full width at half maximum (FWHM) is 2.355$`\sigma `$. The 2-D model has seven free parameters: $`x_{cent},y_{cent},F_{\mathrm{tot}},\sigma ,c,q`$, PA.

The Moffat/Lorentzian Profile. The generalized Moffat function has the following form,

$$\mathrm{\Sigma }(r)=\frac{\mathrm{\Sigma }_0}{\left[1+(r/r_d)^2\right]^n},$$
(14)
with the total flux given by

$$F_{\mathrm{tot}}=\frac{\mathrm{\Sigma }_0\pi r_d^2q}{(n1)R(c)},$$
(15)
where $`r_d`$ is the dispersion radius and $`n`$ is the power-law index. The Moffat profile with $`n`$ of 1.5 or 2.5 is empirically similar in shape to the observed WFPC2 PSF, while $`n=1`$ corresponds to a Lorentzian function. It has FWHM $`=2r_d\sqrt{2^{1/n}1}`$. The Moffat model has eight free parameters in 2-D.

#### 2.3.3 Creating Point Sources

To create a true point source, one normally creates a $`\delta `$-function that is then convolved with the input PSF. But in GALFIT we approximate a $`\delta `$-function by using a Gaussian (or Moffat) function having a small width, usually FWHM $`0.3`$ pixel. The benefit of using a functional representation is that the algorithm can smoothly transit between fitting a true point source or a compact source without a loss of generality.

#### 2.3.4 Creating Bars

Spiral galaxies can have embedded bars, and to model them we prefer to use the Sérsic profile with initially a relatively flat inner and a steep outer profile ($`n<1`$), and a boxy shape ($`c>0`$). The true light distribution of a bar also has a bulbous component at the center, distinct from the bulge. Therefore one should in principle use two components. However, the round component may be partially degenerate with the bulge itself unless the two profiles are significantly different.

### 2.4 Computational Considerations

One design consideration of GALFIT is to efficiently fit well-resolved galaxies in large-format images, using an arbitrary number of model components. In our examples, the nearby galaxies extend beyond the Planetary Camera (PC) chip of 800$`\times `$800 pixels. As will be shown, they can be fit accurately with typically three or more components, involving 20 or more free parameters. Fitting uncropped images is computationally intensive, and convolving entire images can often usurp $`>99\%`$ of the computing time and memory. However, in nearly all cases the seeing has significant impact only near the central few arcseconds of a galaxy. To make fitting large-format images manageable, we convolve only the area most affected by the seeing (Fig. 1, parameters H and I); the process is described below. But the region can be enlarged by the user as required. We find that a convolution radius of at least 20 to 30 seeing disks away from the center suffices for a wide range of nuclear profile shapes. For the HST WFPC2, this can correspond to a region as small as $`40\times 40`$ pixels.

In our implementation, computation time depends on the convolution size, the number of parameters to fit, and the size of the fitted image. Table 1 shows the approximate resources used in running GALFIT for Sérsic and exponential disk fits where all parameters are free, as well as more complicated cases. The total computation time required to reach convergence is estimated for an Intel Pentium III 450 MHz computer with 128 megabytes of memory.

### 2.5 GALFIT Implementation

We give an overview of the inner working of GALFIT to clarify the use of the parameters shown in the top section of Figure 1. Although we will go into significant detail, we emphasize that the the actual fitting process is completely transparent to the user, who needs only to prepare an input template shown in Fig. 1.

GALFIT iterates steps 3-8 until convergence, with further explanations to follow.

1. Normalize and prepare the PSF for convolution (item D in Fig. 1).
2. “Cut out” a section of the image centered on the object to fit from the original data image (item G in Fig. 1).
3. Create model images and derivative images based on new or initial input parameters.
4. “Cut out” the convolution region (items H & I) from the model and derivative images in the previous step, and pad them around the edges with values of the models.
5. Convolve the convolution regions (both model and derivative images) in previous step with the PSF using a Fast Fourier Transform (FFT) technique.
6. Copy the convolution region back into the model/derivative images of step 3.
7. Compare with data image. Minimization is done using the Levenberg-Marquardt downhill-gradient method/parabolic expansion (Press et al. 1997).
8. Iterate from 3 until convergence is achieved.
9. Output images and generate final parameter files.

The Input PSF. For convolution, GALFIT normalizes the input PSF image (item D in Fig. 1) and rearranges it into a wrap-around order (by splitting a PSF into four quadrants, then transposing them to diagonal corners). The PSF is then Fourier transformed using FFT. In GALFIT, all convolutions are done with FFT, a technique found in many numerical computation texts (e.g., Press et al. 1997). FFT reduces the computation time by a factor proportional to $`N^2/(N\text{ log}_2N)`$ compared to brute-force numerical integrations, where $`N`$ is the number of pixels.

For HST observations, we obtain the PSF using either the Tiny Tim (Krist & Hook 1999) software or through the HST data archive. For WFPC2, the plate scale and PSFs are stable, and a number of studies find that Tiny Tim can model the core structures of the PSFs, although the fine details of the diffraction structures may be harder to reproduce. The encircled energy diagrams of the real and observed PSFs are very similar. On the other hand, data from the first year of NICMOS suffered from plate-scale breathing and significant focus changes. The small field of view (FOV) of NICMOS ($`20\mathrm{}40\mathrm{}`$) makes it difficult to find a PSF calibrator in the same image for nearby galaxies. However, Tiny Tim also allows the modeling of time-varying plate-scale and breathing changes and can reproduce the core structure of the PSF. As with WFPC2, the diffraction rings and structures are more difficult to reproduce. For convolution and deconvolution, the accuracy mostly depends on the encircled energy curves and general structures of the diffraction pattern, and to a lesser extent on the fine details. Hence, the use of Tiny Tim PSFs should be adequate for (de)convolution. The user should supply a PSF large enough that the flux amplitude around the image edges is negligible compared to the peak ($`1\%`$); any sky pedestal must also be removed in the PSF image.

Extracting a Sub-Image. GALFIT can either fit the entire image or a sub-region (item G in Fig. 1), allowing reduced computation time for a small object. This is also useful when one might want to fit the nucleus of a galaxy accurately, before enlarging the fit to outer regions.

Creating Model and Derivative Images. To optimize the galaxy profile parameters, GALFIT uses a down-hill gradient method. Unlike most other methods which create only model images, the down-hill gradient technique requires both the model and flux derivative images with respect to all free parameters at each pixel. Because they are stored in memory, the memory usage scales as $`N_{\mathrm{free}}\times N_{\mathrm{pix}}+1`$, where $`N_{\mathrm{free}}`$ is the number of free parameters in the fit and $`N_{\mathrm{pix}}`$ is the number of pixels in an image.

Extracting the Convolution Region. Once the model and derivative images are created, they are convolved with the PSF to account for the seeing. In order to minimize the computation time and memory in the convolution, GALFIT removes a small region (specified in items H and I in Fig. 1) from the model and extends it with additional padding. Padding is needed because of edge effects resulting from convolution (see Press et al. 1997). GALFIT performs two padding operations, first by half the size of the PSF all around the user-specified convolution region. That region is filled with real  model values rather than with zeros — padding with zeros would also corrupt the model around the edges into the convolution region. The image is then further extended with zero padding up to the next $`2^N`$ ($`N`$ integer) number of pixels in both width and length. This is the actual image size used in the convolution. This process is performed for all the models, as well as the derivative images.

Convolving with the PSF. The convolution regions removed from the model and derivative images are convolved with the PSF using FFT. If a PSF is not specified, as when an image has already been deconvolved, the convolution step is bypassed.

Copying the Convolution Region back to the Model and Derivative Images. Once convolution is complete, the model/derivative images are copied back into the original model/derivative images. The two surrounding layers of (corrupted) padding regions are discarded.

Minimizing Residuals. To minimize residuals, we choose the Levenberg-Marquardt method (Press et al. 1997) as the engine instead of the Simulated Annealing algorithm. We discuss our decision in § 3.1. The process of minimization repeats until convergence is achieved, which happens, artificially, when the $`\chi ^2`$ does not change more than 5 parts in $`10^4`$ for 5 iterations.

In complicated fits as the examples we will perform, there is a non-zero chance that GALFIT will fail to produce a good solution. One kind of failure is when the program quits (“crashes”) without a solution because the solution matrix is singular, either caused by poor parameter values, or trying to fit more components into the image than appropriate. GALFIT is generally highly forgiving about poor initial values and frequently converges to good fits ($`\chi _\nu ^21`$-2), down from $`\chi _\nu ^2`$ as high as $`10^5`$-$`10^7`$. When GALFIT does crash, more often than not it is caused by numerical overflows when the the power law indices ($`\alpha `$ and $`\gamma `$ for Nuker, and $`n`$ for Sérsic) become too big ($`>10`$) or too small ($`0.01`$). Another kind of failure is when the solution settles into a local $`\chi ^2`$ minimum and does not get out, but this will not cause GALFIT to crash. For example, in B/D decompositions, there are correlations between the scale-length and luminosity parameters (see, e.g., Byun & Freeman 1995), while ellipticities, centers, and position angles are generally more decoupled and better constrained by the data. We discuss the issue of degeneracy in more detail in § 3.3.

Output Images. Upon completion of the minimization, a FITS image block is created, consisting of three images: the original, the model, and the residual. The software provides the option to not subtract a fitted profile, thereby leaving it in the residual image (item Z).

### 2.6 Parameter Uncertainty Estimation

Once a fit has been optimized, GALFIT estimates the uncertainties analytically. We assume that the surface of constant $`\mathrm{\Delta }\chi ^2`$, as a function of n free parameters, away from the minimum $`\chi ^2`$ can be approximated by an n-dimensional ellipsoid. This allows the uncertainties to be obtained from the covariance matrix, defined as the local curvature of the $`\chi ^2`$ surface at $`\chi _{min}^2`$ with respect to each parameter such that $`\text{C}_{ij}=\text{d}^2\chi ^2/(\text{d}a_i\text{d}a_j)`$, where $`a_i`$ and $`a_j`$ are the fitted parameters. The standard 68% bound on the confidence interval of the fitted parameters, individually, is within a boundary marked by a shell at $`\mathrm{\Delta }\chi ^2=1`$ (e.g. Press et al. 1997). Figure 4 shows an example for a 2-parameter ($`a_1`$ and $`a_2`$) fit.

When the parameters are uncorrelated, the off-diagonal entries of the covariance matrix are zero, and the ellipsoid has major axes running parallel to the unit coordinate axes $`\widehat{a}`$. Formally, the parameter uncertainties are related to the eigenvalues of the covariance matrix, and the standard 68%, 1-$`\sigma `$ uncertainty is $`\sigma a_i=\sqrt{2/\text{C}_{ii}}`$.

However, when some of the parameters are correlated, and the $`\mathrm{\Delta }\chi ^2`$ ellipsoid axes are no longer aligned with the parameter axes. Figure 4 shows an example for a 2-parameter ($`a_1`$ and $`a_2`$) fit, where $`a_1`$ and $`a_2`$ are correlated. Generalizing to higher dimensions, to estimate uncertainties in one parameter we first obtain the eigenvectors and eigenvalues of the n-dimensional hyper-ellipsoid. The semi-major axis vector $`i`$ is given by:

$$\stackrel{}{v}_i=\sqrt{\frac{\mathrm{\Delta }\chi ^2}{\lambda _i}}\widehat{v}_i,$$
(16)
where $`\lambda _i`$ is the eigenvalue corresponding to the eigenvector $`\widehat{v}_i`$. We define the uncertainty for parameter $`a_i`$ by vector summing all the semi-major axis vectors in such a way that the projection onto a given axis $`\widehat{a}_i`$ is largest. In Figure 4, these are vectors represented by the arrow-tipped dotted lines. We then multiply this quantity by a factor of $`\sqrt{\chi _\nu ^2}`$, i.e. we renormalize the errors, to avoid underestimating the errors when $`\chi _\nu ^2`$ is not ideally 1. Written explicitly,

$$\sigma a_i=\sqrt{\chi _\nu ^2}\underset{j=1}{\overset{n}{}}|\stackrel{}{v}_j\widehat{a}_i|.$$
(17)
Strictly speaking, our definition of the uncertainty in Equation 17 is equal to, or slightly larger than, the true uncertainties (which are $`\sigma a_1\mathrm{}`$ and $`\sigma a_2\mathrm{}`$ in Figure 4). We settle for this possible overestimate, but do not consider it to be serious, because the fundamental limitation in determining precise errors is that we do not know how the $`\chi ^2`$ “valley” may twist near the best fit $`\chi ^2`$ or how well it is represented by a hyper-ellipsoid. Furthermore, the parameters we fit may not be all normally distributed. We verified using artificial data with Poisson noise that our method produces reasonable uncertainty estimates for simple models consisting of one and two input components.

Uncertainty estimates quoted in the text and table below are obtained using the method described unless otherwise stated.

## 3 GALAXY FITTING: GENERAL CONSIDERATIONS

### 3.1 The Choice of the Minimization Algorithm

The uncertainty estimate just described is only valid when the solution has reached a global minimum. However, the $`\chi ^2`$ topology for galaxy fitting is complex with many local minima (degeneracies) because the parameter space is large (e.g., $`16`$ parameters for two or more Sérsic components). While no algorithm can guarantee convergence on the global minimum $`\chi ^2`$, sophisticated Simulated Annealing (also known as Metropolis or Annealing for short) algorithms are often the engine of choice for automating B/D decompositions, where they are reputed to be robust. In contrast, though vastly more efficient, the gradient method may head downhill blindly, regardless of whether it is going toward a global $`\chi ^2`$ minimum. Despite this shortcoming, we use the downhill-gradient/parabolic expansion method of Levenberg-Marquardt (Press et al. 1997) as our minimization engine, for reasons given below.

An algorithm that uses Simulated Annealing is GIM2D, developed by Simard (1998) for two-component B/D decomposition. In this method, the annealing (i.e., cooling) “temperature” step size controls the rate of convergence on a solution. At each iteration, a set of parameters is perturbed randomly by some amount; high temperatures correspond to large perturbations. Then, the probability $`P_1`$ that the parameter set is the true one is calculated. If this new probability is higher than the previous, $`P_0`$, the new parameter set is adopted. However if the new set is less likely, there is nonetheless a finite chance $`P_1/P_0`$ that it is adopted. At some point the temperature is decreased and the iteration continues until convergence. The application of the software is demonstrated in several studies (e.g., Marleau & Simard 1998; Simard 1998).

Fitting more than two components dramatically increases the complexity of the $`\chi ^2`$ topology, and extending the Annealing technique to fit more components would seem to be sensical. However, in practice it is tricky to find an optimum cooling rate that can generally accommodate the different number, and type, of components a galaxy might require. Anneal too quickly, the solution settles into a local $`\chi ^2`$ minimum; too slowly, the program converges even more inefficiently. Efficiency is an issue because increasing the number of parameters makes it exponentially more cumbersome to adequately sample the parameter space for calculating probabilities. As an extreme example, to fit the double nucleus of M31 accurately, one needs 6 components, having 41 to 50 free parameters (Peng 2002). Moreover, for highly resolved nearby galaxies, one generally has to create dust masks and test out the components iteratively. This occasional need for interaction makes both the speed and the straightforward implementation of the downhill-gradient method attractive.

The downhill-gradient method may not be quite as well optimized for B/D decomposition as the Simulated Annealing method. However, a direct comparison between the two methods is needed to evaluate their respective strengths and weaknesses. We discuss below (§ 3.3) how to test for degeneracies. It may also be possible to tailor a hybrid gradient-search and Annealing method, which tries to climb out from a local minimum by briefly using Monte-Carlo sampling after the solution has been optimized by the gradient search. This is being considered for future updates to GALFIT.

### 3.2 Fitting Procedures

Most bulges and late-type galaxy nuclei that are well resolved by HST appear as smooth and simple spheroids on large scales when one disregards the occasional dust. But, as we will show below, simple one or two component fits frequently produce large residuals that have bipolar or quadrupolar symmetry, which can be reduced with additional components. This subtlety makes estimating the required number of components a trial-and-error process, and Occam’s razor suggests using the fewest number necessary, based on the goal of the experiment. We determine the number and types of components iteratively, often by starting out with a Sérsic or a Nuker profile for ellipticals, and a Sérsic plus an exponential for spirals. Based on $`\chi ^2`$ and the pattern of residuals, we determine if we need to replace them with more flexible models or add in more components.

To use GALFIT in a typical manner such as for B/D decomposition, we point out that there is little to no human interaction beyond preparing an initial parameter file. However, our task of accurately fitting nearby large galaxies is considerably more challenging, thus we outline the steps we normally take:

1. Create a dust/bad pixel mask by hand. We create a mask by outlining the affected region with a polygon, then feed the list of pixels into GALFIT either as a text file or a FITS image.
2. Estimate by eye the number of components needed.
3. Although not needed, it sometimes helps to initially hold the diskiness/boxiness parameter $`c`$ fixed until the algorithm has found a plausible solution. The radicands in the radius Equation  have absolute value quantities; hence, the derivative images contain regions of discontinuity that cause the downhill-gradient method to react sensitively unless near a minimum.
4. If convolution is needed, leaving it off initially saves time until a reasonable solution has been attained.
5. If a galaxy is difficult to fit near the center, we restrict fitting to that region until a good solution is found, then enlarge the fit.
6. Optimize all the parameters in a $`\chi ^2`$ sense.
7. Examine the residuals to decide whether more components need to be added, make better pixel masks, free parameter values that have first been held fixed, or turn on convolution.
8. Repeat previous steps as necessary.

### 3.3 Degeneracy and the Significance of the Components

In this paper and in other related studies (Peng 2002; Ho et al. 2002), we find that galaxy light profiles can generally be modeled accurately (in a $`\chi ^2`$ sense) using three to five components. The larger-than-usual number of components immediately prompts the question of how unique is the decomposition.

Degeneracy is a common problem in galaxy decompositions, mostly because the model functions used are not selected based on physical criteria. One can, in principle, decompose a galaxy with as many “basis functions” as one might contrive. Even when a set is well defined, such as in two-component B/D decompositions, the large number of parameters involved makes degeneracy sometimes an issue, as shown in the simulations of Byun & Freeman (1995) and Wadadekar et al. (1999). Extending to more components might, at first sight, make the situation unmanageable. One way to reinforce our confidence in parameter uniqueness is to use Monte-Carlo simulations to search the parameter space, but these say nothing about whether a given combination of models is the appropriate or unique set. Moreover, searching randomly may be unfeasible because the dimensionality is over 20, and occasionally over 40. However, while there are potentially many local $`\chi ^2`$ minima solutions, in practice the situation is not hopelessly degenerate. Imposing the demand that models fit galaxies with high fidelity ($`\chi _\nu ^21`$) severely reduces the possible solution space for a given set of models; most solutions can be easily dismissed based on $`\chi ^2`$ analysis. Such demand produces components that have significantly different scale lengths, axis ratios, etc..

While there is no recipe to guarantee uniqueness in the fitting, there are ways to probe the parameter spaces exhaustively and to satisfy to one’s own confidence that a solution is stable and insensitive to initial conditions. One way is to vary two parameters at a time on a Cartesian grid of values in order to trace out the local $`\chi ^2`$ contours. This technique is very useful when the number of parameters is low. However, the same technique quickly becomes intractable with increasing number of fitting parameters. To more effectively explore a wider $`\chi ^2`$ topology, we propose the following technique of combining Monte-Carlo simulations and minimization.

1. Randomly select a full set of parameters, possibly drawing from distributions of parameter values centered on the best-fit values.
2. Using these as initial values, minimize the $`\chi ^2`$.
3. Repeat step 1 for many different sets of random initial guesses to see if they return the same optimized solution or to other equally plausible ones.

The degeneracy situation, together with the use of more than three components, makes interpretation difficult, leading one to ask: “what do all the components mean?” or “why use so many components?” The answer depends on the prior goal: Perhaps the most mundane answer is simply that a bulge is triaxial, or otherwise not well represented by our limited models. In this case two or more functions may be needed to describe the same entity that has meaning only when summed, but not individually. A different, but perhaps more interesting tactic is to explicitly seek out components motivated by theory. A case in point is the prediction by Naab & Burkert (2001) that there should be disks embedded in elliptical galaxies formed in galaxy mergers. The use of extra components may also be motivated by external data, such as evidence for a nuclear point source based on AGN spectra, or decoupled cores based on kinematic data. To look for weak AGN point sources, it is especially important to deblend a large portion of an image accurately, rather than biasing the fit to regions that one can fit well. Lastly, from a purely empirical standpoint, one can decompose a galaxy simply to look for structures too subtle to be seen in full light, such as weak nuclear bars, stellar disks, and nuclear point sources. Again, accurate decomposition here is crucial so that the residuals do not undermine the detection or believability of weak sources.

The last scenario brings the discussion back in full circle to the issue of degeneracy. Despite possible degeneracies in decomposing large-scale components, smaller structures are better localized, and better defined in size and shape. Even if there is doubt as to what the exact profile and parameters are, often different models can reproduce a component with similar overall characteristics (e.g., profile, scale-length, shape, and orientation), even if they are not exactly identical. That a similar component is repeatedly borne out through different profile assumptions is in principle the keystone of its reality, since its inclusion is essential to achieving a reasonable model fit in that region of the image. We show a few such cases below.

## 4 MODEL FITTING: A CASE BY CASE STUDY

In the discussion to follow, we present GALFIT decomposition of galaxies having a wide variety of properties and shapes in order to test how well nearby galaxies can be modeled in 2-D by parameterized functions. We apply GALFIT to images of 7 galaxies that have a wide range of interesting properties: different shapes, morphologies, isophote twists, and internal disks. Our WFPC2 data come from the HST archive. We remove the edges which are affected by CCD readouts, and fit models to $`721\times 721`$ pixel images. The galaxy nuclei in these images have very high S/N and are centrally placed in the PC. They often extend well beyond the WFPC2 FOV, and the exposures are short (several hundred seconds, see Table 1), making the sky background and image bias level hard to estimate from the Wide Field chips. Given the typical exposure times of our images, the 26 mag arcsec<sup>-2</sup> isophote is about 0.01 – 0.05 ADU. We set the sky to zero in our fitting, unless specifically mentioned otherwise.

Figures $`510`$ and Figure $`14`$ show the original images as well as the residuals in positive greyscale. We show two panels of residuals. The first (panel b) illustrates the traditional decomposition technique using only a one-component bulge and disk, occasionally adding a bar and a point source when obviously needed. We use a Sérsic or Nuker function instead of de Vaucouleurs model for the bulge. The second residual panel (panel c) shows a more accurate decomposition using the components listed in Table 3.

Beneath the images are plots of the surface brightness profile, PA, ellipticity, and cos $`4\theta `$ obtained by running the IRAF<sup>1</sup><sup>1</sup>1IRAF (Image Reduction and Analysis Facility) is distributed by the National Optical Astronomy Observatories, which are operated by AURA, Inc., under cooperative agreement with the National Science Foundation. task ellipse. The surface brightness plot shows the observed profile (i.e., without deconvolution) as solid circles with error bars. The short-dashed line going through those points is a 1-D Nuker law that best fits the observed profile; this is merely to illustrate the flexibility of the function, as well as to guide the eye for deviations from smoothness or power-law behavior. The long-dashed lines show the intrinsic profile of each component listed in Table 3 from 2-D decomposition. Finally, the solid line is the net sum of all the components that make up the intrinsic profile. Therefore, it will rise above the observed points at the nucleus in cases when there is no point source. The dotted vertical line in these plots shows the region most affected by PSF smearing, but in some cases the effect can be seen as far out as $`0\stackrel{}{\mathrm{.}}6`$ (e.g., NGC 221).

The ellipticity is defined as $`ϵ=1q`$, where $`q`$ is the axis ratio, and the cos $`4\theta `$ plot shows whether the isophotes are disky ($`>0`$) or boxy ($`<0`$), similar in definition to the $`c`$ parameter in GALFIT, but opposite in sign.

Table 3 summarizes all the fitting parameters; the magnitudes have not been corrected for Galactic extinction. The magnitudes are based on the synthetic zeropoint calibrations by Holtzman et al. (1995), which transform data units onto the Landolt UBVRI magnitude system. Table 3, column 15, lists the $`\chi _\nu ^2`$ for the best fit shown in Figures c, and column 16 lists the increase in $`\chi _\nu ^2`$ for fits shown in Figures b. We do not quote uncertainties for all the parameters individually. Instead, in the last row of Table 3, we quote representative uncertainties derived in the fit (see § 2.6) as an ensemble. For most fits, “representative” means a conservative upper bound on the fitting parameters.

### 4.1 A Caveat on Galaxy Morphology Comparisons

Comparing galaxy morphology/decompositions using inhomogeneous studies is difficult if they are based on images of limited FOV and S/N. Galaxies in low-S/N images have wings that are barely detectable above the sky noise, often biasing scale-length measurements toward lower values, if the noise is not accurately accounted or the suitable function is unknown. While low S/N in the wings is offset by a larger number of pixels, there is a practical tradeoff between areal increase and the rate of profile decline in the presence of sky and detector noise. In decomposing nearby galaxies there is an additional concern if they extend well beyond the FOV of the CCD. Since galaxy profiles tend to change shape, PA, and ellipticity with radius, in addition to the degeneracies discussed above, it is generally unsafe to extrapolate profiles beyond the FOV.

The galaxies in our examples below were all observed with WFPC2. And although they have very high S/N in the core, they have short exposure times (several hundred seconds). Moreover, they show only the central regions of the galaxies. While they are perfectly suited for studying galaxy nuclei, profile comparisons with larger FOV images are not straightforward. When we fit components whose scale-lengths are comparable to, or larger than, the FOV, we do not necessarily believe the numbers to be accurate. However, the reality of such a component is often suggested by analysis of the residuals.

We leave further interpretations to future studies (Ho et al. 2002).

### 4.2 A “Simple” Elliptical, NGC 221 (M32)

This E2 galaxy has an apparent axis ratio of 0.72. Michard and Nieto (1991) find that M32 has pointed isophotes between 1$`\stackrel{}{\mathrm{.}}`$2 and 5$`\mathrm{}`$. Color images obtained with HST and presented in Lauer et al. (1998) show strong surface brightness fluctuations, which were studied by Ajhar et al. (1997). The surface photometry of M32 is shown in Figure 5, uncorrected for PSF smearing. Lauer et al. (1992, 1998), through deconvolution of WFPC and WFPC2 data, find that the ellipticity is constant throughout the inner portions of the galaxy. They see no signs of significantly boxy or disky isophotes, nor evidence for dust, disks, or other structure down to a few percent in local surface brightness.

In 1-D, the Nuker function fits the surface brightness profile remarkably well out to $`r\mathrm{\hspace{0.17em}20}\mathrm{}`$. However, a 2-D Nuker fit leaves behind large residuals near the center (Fig. 5b). We need two components (see Table 3) to produce a good fit (Fig. 5c). The requirement of two components is not due to ellipticity and PA changes; in fact, in agreement with Lauer et al. (1998), the PA, ellipticity, and diskiness/boxiness parameters all appear to be constant with radius in our two components. Rather, two components are needed because of a subtle upturn in the profile in the region $`0\stackrel{}{\mathrm{.}}72\mathrm{}`$. This is nearly indiscernible in the 1-D surface brightness profile of Figure 5, but much more so in Figure 5b where a one-component Nuker fit leaves a “doughnut-shaped” residual, with an increase in $`\chi _\nu ^2`$ by 0.13 from the best fit, Fig. 5c.

Much of the structure can be removed by fitting two components to the image. Figure 5c shows a Sérsic + exponential disk fit. A reasonable fit can also be obtained with a de Vaucouleurs bulge plus a nuclear disk described by Nuker profile. We adopt the first set of parameters because formally it is better. Despite the excellent fit, the residuals show a very low diffuse positive halo out to $`1\stackrel{}{\mathrm{.}}3`$, at a level about $`0.1\%1\%`$ of the local flux. The need for two components point to the intriguing possibility that there might be a nuclear stellar disk with a semi-major axis of $`r_s1.5`$ pc. There are also slight negative residuals around ($`20\mathrm{},18\mathrm{}`$), which may suggest the presence of dust at a level $`3\%8\%`$ below the local mean intensity, which is not seen in the region reflected about the major axis. Hence, this is probably not caused by a profile mismatch because mismatches would produce bipole or quadrupole patterns. We confirm Michard & Nieto’s (1991) finding that M32 has a slight disky component, as evidenced in both Figure 5 and in our $`c`$ parameter ($`c=0.06\pm 0.01`$).

### 4.3 “Disky” Galaxies

#### 4.3.1 NGC 4111

This nearby edge-on S0 galaxy is distinguished by highly disky isophotes, a peanut-shaped nucleus (pinched at the minor axis), and a complex surface brightness profile (Fig. 6) having multiple breaks and inflections at large radii. Burstein (1979) and Tsikoudi (1980) analyze photographic plate scans of NGC 4111 that cover 1–4 in radius. Upon decomposing it in 1-D, they find several distinct photometric components. Tsikoudi (1980) finds that the inner spheroidal component with a classic $`R^{1/4}`$ profile has an effective radius $`r_e=15\mathrm{}`$ along the semi-major axis; it contributes 48$`\%`$ of the total luminosity in the $`B`$ band. The remainder of the galaxy luminosity is contained in the “lens” component of the galaxy and in an exponential disk that extends beyond $`102\mathrm{}`$. There are several profile inflections at large radii, in particular at $`a30\mathrm{}`$ and $`70\mathrm{}`$; Tsikoudi suggests that the hump between $`30\mathrm{}<a<102\mathrm{}`$ may be a ring or weak spiral arm. Cross-cuts perpendicular to the major axis (perpendicular profile) at $`10\mathrm{}`$, $`30\mathrm{}`$, and $`50\mathrm{}`$ from the center show that the exponential envelope itself is composed of three components, where the inner (which extends out to 12$`\mathrm{}`$) and intermediate regions ($`<24\mathrm{}`$) are Gaussians, while the outer region is characterized by an exponential disk with a steep drop-off, extending well beyond $`60\mathrm{}`$ (Tsikoudi 1980). Similarly, Burstein finds that the perpendicular profiles of NGC 4111 can be described by an exponential disk profile, and that the rounded center can be fitted by a bulge component. Moreover, the residuals suggest there is a “thick disk” component, which is also seen in other S0s in his sample.

The innermost 20$`\mathrm{}`$ of NGC 4111, as seen with HST, proves to be just as complicated as the region exterior, with several breaks in the light profile (Fig. 6). 2-D decomposition reveals that much of the dip at $`r=1\mathrm{}`$ is caused by a circumnuclear dust ring. The dust lane bisecting the nucleus is not fully symmetric in shape, and gives rise to the peanut appearance. The two-component B/D decomposition (Nuker and exponential) presented in Figure 6b shows large over- and under-subtractions. There is a prominent disky component with a boxy nucleus in the residuals. Moreover, at the junction between the two, the disk appears to be pinched.

A more accurate decomposition of the nuclear region requires four disky components and a boxy component. The five components give an excellent fit to the nucleus, despite all the complexities, and we claim that at least four are real and photometrically distinct. Of the two dominant components, one is an exponential disk ($`r_s=21\stackrel{}{\mathrm{.}}4`$) with a disky shape ($`c=0.65`$). The other is a Sérsic ($`r_e=7\stackrel{}{\mathrm{.}}85`$) with $`n`$ = 2.51 and a boxy shape ($`c=0.28`$), giving rise to the peanut-bulge appearance. There is also a nuclear stellar disk that stands out very prominently in the decomposition. It consists of two components that have nearly identical half-light radius of $`r_e=4\stackrel{}{\mathrm{.}}5`$: a Gaussian halo ($`n=0.48`$) of $`q=0.23`$, and a thin ($`q=0.14`$), flat disk with a steep fall-off ($`n=0.19`$) in the profile beyond $`r_e`$. Together these two components add up to a total apparent magnitude of $`m_{\mathrm{F547M}}=13.1`$, or roughly 10% of the entire $`V`$-band output of the galaxy. In the disk plane, about $`5\mathrm{}`$ on both sides from the center, there appears to be over-subtractions perpendicular to the major axis of the galaxy. Because they are not fully symmetric across the galaxy plane nor across the semi-minor axis, its suggests that there is another set of dust lanes that gives rise to the pinched appearance in the disk in Figure 6b. Our decomposition shows that the bulge does not have a spherical component with a de Vaucouleurs profile, as suggested by Tsikoudi (1980), but rather is a complicated blend of several disky and boxy components.

#### 4.3.2 NGC 4621

NGC 4621 is an E5 galaxy with a very steep nucleus where $`\gamma =0.8`$. The nuclear cusp slope is defined by $`\gamma \text{dlog}(I)/\text{dlog}(r)`$ within $`r<0\stackrel{}{\mathrm{.}}1`$. The isophotes are disky but become gradually less so at large radii; the axis ratio decreases outward from 0.88 to 0.64 in the same regime, and the PA remains essentially constant with radius. There are a number of globular clusters in the WFPC2 F555W image, but there is no evidence for nuclear dust. Mizuno & Oikawa (1996) perform 2-D decomposition of a ground-based image of this galaxy and find that it is made up of three distinct photometric components: an exponential disk, a de Vaucouleurs bulge, and an outer envelope. The exponential disk and de Vaucouleurs bulge, which account for $`16\%`$ and $`62\%`$ of the total light, respectively, are both disky in shape, while the outer envelope is spherical and contains $`22\%`$ of the galaxy luminosity. Scorza & Bender (1995) and Simien & Michard (1990) use two-component models, but they arrive at different conclusions regarding the galaxy structure. Simien & Michard (1990) find the disk scale-length to be $`28\mathrm{}`$, while Scorza & Bender (1995) obtain $`7\stackrel{}{\mathrm{.}}5`$. The discrepancy between the two studies may be due in part to the use of different size images in the decomposition.

Figure 7b illustrates a fit using a Nuker bulge and an exponential disk ($`r_s=4\stackrel{}{\mathrm{.}}7`$). The fit produces residuals that have high symmetry. Their diffuseness and symmetry, suggest there is an intrinsic profile and shape mismatch rather than there being additional components, in contrast to the sharp features seen in Fig. 6b of NGC 4111.

Despite the changing galaxy shape and inflections in the surface brightness profile, the entire image of NGC 4621 on the PC chip (inner $`25\mathrm{}`$ radius) can be subtracted remarkably well by employing two Sérsic components and a Nuker component (Fig. 7c). The two innermost components have similar scale-lengths, even though one is disky ($`c=0.35`$) and steep ($`n=5.98`$) at the center and the other is boxy ($`c=0.78`$) and flat ($`n=0.89`$). Both are required to fit the “eye-shaped” nucleus well, but our combination may not be unique. Hence, it is possible that the two may add to be a single photometric component. The outer envelope is characterized by a disky ($`c=0.16`$) Nuker profile with $`q=0.70`$. This component is not the spherical halo component found by Mizuno & Oikawa (1996) because our FOV is much smaller. In conjunction with their finding of a round galaxy halo, we support the notion that NGC 4621 has $``$3, possibly 4, photometric components.

Upon casual visual inspection, the nucleus of NGC 4621 may appear to be pointlike because the isophotes are round and the profile is steep ($`\gamma =0.8`$). But the three-component model fits the nucleus cleanly without the need to add in an unresolved point source.

### 4.4 A “Boxy” Galaxy, NGC 5982

The nuclear region (within $`30\mathrm{}`$) of this E3 galaxy has boxy isophotes that become increasingly round at $`r3\mathrm{}`$, but is otherwise uniform with no sign of dust lanes (Figure 8). The kinematics of this galaxy have been studied in detail because the stellar velocities reverse direction $`2\mathrm{}3\mathrm{}`$ from the center, appearing as though the galaxy has a kinematically distinct core with minor-axis rotation (Wagner, Bender, Möllenhoff 1988). One scenario to explain kinematically distinct cores invokes mergers of galaxies that preserve their initial angular momentum vectors. An alternative is proposed by Statler (1991), who shows that dynamically smooth, circulating orbits in a triaxial potential can, in projection, cause the appearance of kinematic twists such as those observed in NGC 5982. To distinguish between the orbital projection model from the merger model, Statler (1991) predicts that the former should produce a minimum in the velocity field along a direction somewhere between the minor and major axis. Moreover, there should be a photometric core that is similar in size to the kinematic core. Observations by Oosterloo, Balcells, & Carter (1994) indeed suggest that there is such a trough present in the velocity field at PA = $`78^{}`$. But the ground-based surface photometry of this galaxy published by Michard & Simien (1988) and Bender, Döbereiner, & Möllenhoff (1988) do not reveal a photometric core down to at least $`3\mathrm{}`$. Although the kinematic peculiarity may be due to projection, Oosterloo et al. (1994) suggest that the two effects: orbital projection and the merger of a small clump, may not be mutually exclusive events.

With HST images, we find that the radial profile nearly flattens off into a core interior to $`0\stackrel{}{\mathrm{.}}6`$ (Fig. 8). There also appears to be a subtle break at $`r3\mathrm{}`$, where the nucleus becomes significantly rounder. The 3″ radius also corresponds closely to the peak of the rotation curve. The finding that the photometric core is at least an order of magnitude smaller than 7″ rules out Statler (1991) models (T. S. Statler, 2002, private communication). However, while the kinematically decoupled core may not be best explained by those models, the orbital projection effects proposed by Statler (1991) may still play a role.

Fig. 8b shows a one-component Nuker fit that produces large quadrupole residuals at the center, but a relatively good fit farther out. The decomposition of this galaxy is complex and requires four components because of the large ellipticity gradient that tends toward zero near the center. Nonetheless, the galaxy can be accurately decomposed (Fig. 8c); the slight quadrupole residual pattern (with RMS fluctuations 1%-3% of the galaxy peak). The ellipticity gradient in Figure 8 might be due to triaxiality. However, purely triaxial galaxies should also have confocal isophotes, but we find that the first component is significantly displaced from the rest by ($`\mathrm{\Delta }`$RA, $`\mathrm{\Delta }`$DEC) = ($`0\stackrel{}{\mathrm{.}}38,0\stackrel{}{\mathrm{.}}16`$), suggesting there may have been a merger event proposed by Oosterloo et al. (1994). Forcing all the components to have the same position produces asymmetric residuals near the center, and increases $`\chi _\nu ^2`$ by $`3\times 10^3`$ with $`N_{\mathrm{dof}}=5.2\times 10^5`$.

About $`13\stackrel{}{\mathrm{.}}9`$ (2.6 kpc in projected separation) from the center of NGC 5982, there is a galaxy that appears to have a double nucleus, although there is a remote chance it may be caused by dust. The dwarf galaxy’s orientation points nearly directly toward the center of NGC 5982. If this galaxy is indeed within the halo of NGC 5982, not just due to projection, one is tempted to speculate that the dwarf companion is in the process of being tidally disrupted by the larger galaxy.

To summarize, the profile of NGC 5982 flattens out into a core at $`0\stackrel{}{\mathrm{.}}6`$, which is $`10`$ times smaller than that predicted by Statler (1991). Although this finding rules out Statler (1991) models, similar projection effects may still play a role in explaining the kinematically decoupled core. Our finding of an off-centered, luminous, component at the core is probably evidence that this galaxy has had a recent merger episode.

### 4.5 Spiral Galaxies

While early-type galaxies can be modeled well by ellipsoids, the obvious complication when fitting spiral galaxies is the presence of spiral arms and large-scale bars that are embedded in at least 50%–60% of the objects (e.g., Eskridge et al. 2000). Spiral arms cannot be easily parameterized with our models, but often there is a diffuse disk of older stars underneath the arms that can be well described by an exponential disk (e.g., de Jong & van der Kruit 1994; Prieto et al. 2001), whose contrast with the arms is smaller in the near-infrared. The situation is different for galactic bars, as they can be fitted with a superposition of one or more analytic functions. We show a decomposition of an ordinary spiral and a barred galaxy.

#### 4.5.1 NGC 4450

NGC 4450 (Fig. 9) is an Sab galaxy viewed at a low inclination angle; its LINER nucleus exhibits broad emission lines (Stauffer 1982; Ho et al. 1997b), which appear double-peaked in HST spectra (Ho et al. 2000). We will later revisit NGC 4450 in the context of point-source extraction (§ 5.2.3), but for now we focus only on the global galaxy morphology. The 1-D decomposition by Baggett et al. (1998), assuming a spherical bulge and a circular disk, gives a B/D ratio of $`0.3`$.

The large scale disk of the galaxy extends well beyond the PC chip, thus we first fit the entire WFPC2 mosaic image to get a better handle on the disk parameters. To do so, we mask out the echelon region missing in the WFPC2 FOV, and fit the remaining three quadrants plus the PC chip by assuming the components are individually axisymmetric. Then, holding the parameters of the disk fixed, we fit only the PC chip while allowing the nuclear components to optimize. Decomposing the galaxy using the classic two-component model (de Vaucouleurs and exponential disk) gives B/D = 0.23, with $`r_e=24\stackrel{}{\mathrm{.}}0`$ and $`r_s=72\stackrel{}{\mathrm{.}}2`$, but the fit is poor. In comparison to Baggett et al. (1998), our scale-length measurement is about twice as large. Part of this large difference may be attributed to our inability to determine an accurate background sky value because of the limited FOV of WFPC2; we set the sky level to 26 mag arcsec<sup>-2</sup> (0.07 ADU).

An initial two-component decomposition of the WFPC2 mosaic reveals what appears to be a bar. This identification is reinforced in a three-component fit (Sérsic bulge, exponential disk, and a bar) that shows the bar to have a flat profile within $`r_e`$ but steep fall off beyond it ($`n=0.11`$). However, the identification may not be secure because of the missing quadrant in the mosaic. For the purpose of decomposing the PC chip, this is not crucial, so we assume that such a component exists. After optimizing the parameters for the disk and bar, we hold them fixed when fitting the bulge on the PC chip. Through trial and error, we find that the bulge itself is better fit by summing two Sérsic components; neither one alone should be regarded with high significance because the bulge is severely affected by dust. If one were to fit the bulge with a single component, Figure 9b shows that the northern and southern regions of the nucleus would be oversubtracted. The central dark patches are probably oversubtractions, rather than dust obscuration, because the round and bipolar shape of the residuals contrasts sharply in appearance with the irregularity of the surrounding dust lanes. With some confidence in the global components, we fit the remaining point source on top of the nucleus, and the result of the fit is listed in Table 3.

To summarize, we use five components for NGC 4450: a bar, a disk, a nuclear point source, and a bulge made up of two Sérsic functions. In this decomposition, the bulge is somewhat fainter than the standard B/D analysis, in part because of the contribution of the bar, but also because the nucleus is steeper than a de Vaucouleurs profile. We find B/D $``$ 0.09 if the bar is left out, or B/D $``$ 0.13 if it is included with the bulge.

#### 4.5.2 A Barred Galaxy, NGC 7421

NGC 7421 is an SBbc galaxy seen at low inclination. The bulge and the bar are displaced from the center of the disky envelope, while the western bar truncates in a “bow-shock”-like spiral arm. H I column-density maps show a “wake” that suggests that NGC 7421 is plowing through, and interacting with, the intracluster medium in which it is embedded (Ryder et al. 2000).

Decomposition of this galaxy in 2-D requires at least four components in order to fit a bulge, a disk, and a bar (Fig. 10). As with other spiral galaxies, obtaining a plausible disk contribution requires using a FOV beyond that of the PC. So we proceed in a similar fashion to NGC 4450, by first fitting the mosaic image, followed by optimizing over the PC region separately while holding the large-scale disk component fixed. A sky pedestal of 21.6 mag arcsec<sup>-2</sup> (1.1 ADU) was determined from the edges of the mosaic. Removing the model of the galaxy reveals a faint dust lane that runs along the bar through the nucleus. Moreover, there is an unresolved compact nuclear component that can be characterized by a Gaussian; this component has also been noted by Carollo et al. (1998). The “bulge” is composed both of an intrinsic galaxy bulge plus the round component of the bar (components 2 and 3); the two are difficult to distinguish. Without the fainter of the two bulge components, Figure 10b shows doughnut-shaped residuals. The luminosity of the summed components of the “bulge” rivals that of the bar (component 2), while the single most luminous component is the exponential disk. The bulge-to-disk ratio, excluding the bar component, is B/D $``$ 0.06; including the bar with the bulge, B/D $``$ 0.15.

## 5 EXTRACTION OF NUCLEAR POINT SOURCES

### 5.1 General Considerations

A common application of galaxy decomposition is to extract reliable measurements of nuclear point sources embedded in a host galaxy. While luminous quasars can overpower the light of the host galaxy, the opposite extreme is the case for the more garden-variety AGNs in nearby galaxies, where the nucleus is overwhelmed by the glare of the galaxy bulge. Photometric decomposition is absolutely crucial to properly disentangle the nuclear component. To do so, a central compact source or point source can be added on top of the normal galaxy profile to represent the AGN (e.g., Kotilainen et al. 1992; Portal 1999). Implicit in this technique is the assumption that the galaxy profile can be extrapolated all the way into the center by fitting the asymptotic galaxy wings using standard functions such as the de Vaucouleurs plus exponential (Portal 1999) or Gaussian profiles (Carollo et al. 1997).

When fitting faint nuclear point sources, the adopted form of the galaxy profile is crucial because all functions differ widely in steepness near $`r=0`$. This is compounded by the fact it is difficult to fit profiles accurately on large and small scales simultaneously; most of the weight is dominated by outer regions where there are more pixels, at the expense of few pixels in the center. To overcome these complications, it is tempting to simply exclude the outer regions. But in so doing the intrinsic steepness of the galaxy profile is less well constrained. As galaxy nuclei continue to rise as a power law even in the absence of central sources (Lauer et al. 1995), this practice easily produces fake compact source detections – a power law profile can always be artificially decomposed into a steeply rising line plus a compact source. Finally, if even the slowly varying region surrounding the bulge cannot be well fit, the harder task of fitting the steeply rising inner region should be even more dubious. Therefore, the desired goal should be to fit the entire nuclear region as accurately as possible when trying to interpret and extract faint compact sources at the central few pixels.

We first compare the merits of 1-D and 2-D point-source extractions. Some arguments made in favor of 2-D techniques in the context of B/D decomposition are directly applicable here. The 1-D profile extracted from images is affected by several factors: the seeing, the different ellipticities, profiles, and shapes of the various galaxy components, a possible point source, and other nuclear features such as dust. There are few pixels near the nucleus, so the inner isophotes are sensitive to small-scale features such as dust, which can affect the centroiding and the shape of the nuclear profile. Miscentering can lead to uncertainties in the central source profile in images where the PSF may be undersampled, as is the case with WFPC2. As an example of the ambiguity that can be encountered in extracting a weak nuclear point source in 1-D, Figure 11 shows fits of the azimuthally averaged profiles (without deconvolution) of NGC 4143 and NGC 4450. Different assumptions can also lead to conflicting evidence for a compact source. In 1-D, for example, it is impossible to tell whether the inflection at $`r0\stackrel{}{\mathrm{.}}2`$ is caused by circumnuclear dust or is a true sign of a photometrically distinct compact component. It is far easier to gauge the reality of compact sources in 2-D based on the residuals of the model-subtracted image.

To compare the relative effectiveness of 1-D versus 2-D recovery techniques, we added point sources to three real galaxies used in this study. A point source has only three free parameters ($`x`$ and $`y`$ position and brightness), thus is easier to simulate and recover uniquely than compact sources, providing the PSF is accurate. Our simulations give only qualitative guidelines for comparing the relative merits of 1-D and 2-D decompositions; we do not cover all possible scenarios. Other simulations of point-source recovery have been done by Wadadekar et al. (1999), who use idealized galaxies with exponential disks and Sérsic profiles to explore a much wider range of parameter space.

We simulate faint nuclear sources by placing a PSF on top of NGC 221, NGC 4621, and NGC 5982. Figures 5a, 7a, and 8a show the original 2-D images of the galaxies with contours overlaid; Figures 5c, 7c, and 8c give the 2-D model residuals. As shown in Fig. 8d NGC 5982 has a global axis ratio of $`q`$ = 0.68, a slightly boxy isophote shape, and a flat nuclear cusp of $`\gamma =0.04`$. For NGC 221, $`q=0.72`$ and $`\gamma =0.5`$. For NGC 4621, $`q=0.70`$ at large radius and is disky at $`r8\mathrm{}`$; the nucleus is rounded by the seeing ($`q=0.9`$) and has a steep power-law cusp of $`\gamma =0.8`$. These examples have no apparent compact source at the nucleus, no dust features, and span a wide range of nuclear slopes. NGC 4621 appears to have a pointlike nucleus only because it has a very steep cusp and because the nucleus is round. We have shown in § 4.3.2 that NGC 4621 can be fit accurately without requiring a nuclear point source. We can model all three galaxies quite well in both 1-D (using a single Nuker profile) and 2-D.

We add Tiny Tim PSFs with brightnesses 2 to 7 magnitudes fainter than the bulge, in 1 magnitude intervals, then extract them with the same PSF. The bulge magnitude is defined here as the flux within $`r0\stackrel{}{\mathrm{.}}5`$. Figure 12a–15c illustrate the projected profile resulting from adding the point sources. Figure 13 compares the results of the extractions in 1-D versus 2-D, showing deviations between input and recovered magnitudes. The S/N of the images are as high as $`7080`$ at the peak. When the point source is bright and the nucleus is moderately flat, the 1-D technique can recover the point sources well. But the recovery in 1-D becomes poor very quickly for a moderately dim source and a steep profile. For $`\gamma 0.5`$ there is little accuracy even for relatively luminous sources. At the steep end of $`\gamma `$, the technique is completely unreliable.

By comparison, 2-D modeling can recover much fainter nuclear sources, although there is a slight systematic offset. The recovered magnitude in most cases is slightly fainter than the input magnitude, a consequence of the degeneracy between a steep galaxy profile (parameterized using a Nuker function) and a point source. The use of the Nuker function can overfit a steep cusp at the expense of a point source. For faint point sources and steep cusps, the opposite scenario may occur; the recovered magnitude may be brighter than the input at the expense of a shallower cusp.

In AGN studies, it is often informative to know the spatial extent of the nuclear source. Direct non-stellar emission from a central engine is expected to be unresolved, even for the nearest AGNs, whereas reprocessed non-stellar light or a nuclear star cluster may be resolved. To test the source extent, one can fit an analytical PSF function, such as the Moffat profile described in § 2.3.2, with an adjustable width. In our simulations in 1-D, we find that the luminosity and width of the compact source depend sensitively on initial parameters because of degeneracies in the fits. We conclude that 1-D fitting generally gives unreliable measures of compactness for central sources. In our 2-D simulations, the widths converge to that of a point source in every case, except for the faintest test sources.

For the purpose of extracting faint compact sources, we can constrain their parameters by exploring the degeneracies of the local minima in $`\chi ^2`$, and by fitting different nuclear models. In fact, the latter is a necessary exercise because different model assumptions may produce contradictory results because of profile ambiguities. For instance, a compact source on top of an exponential bulge may, in some cases, be mimicked by a single Nuker profile without a compact source. In the next section we give some examples of the 2-D technique of fitting point sources.

### 5.2 Examples

#### 5.2.1 NGC 2787

The SB0 galaxy NGC 2787 contains a LINER nucleus with broad Balmer lines (Ho et al. 1997b). The association of this object with AGN activity is corroborated by the recent HST detection of a central dark mass of $`7\times 10^7M_{}`$, presumed to be a massive black hole (Sarzi et al. 2001). The smooth and otherwise featureless bulge is beautifully adorned by a circumnuclear dust disk with ripples that lie nearly orthogonal to the major axis of the bulge (Fig. 14). Isophotal analysis shows a dip between $`1\mathrm{}6\mathrm{}`$ in the profile cut and shape parameters, caused by the dust disk. The bulge itself has an ellipticity of 0.32 at large radii with a large and fluctuating gradient; the gradient is due in part to dust.

In 2-D, the bulge can best be fit by an exponential and a Nuker component. Although the bulge component can also be well fitted by a superposition of a Sérsic ($`m_{\mathrm{F547M}}=12.43`$ mag, $`n=2.71`$, and $`r_e=7\stackrel{}{\mathrm{.}}43`$) and an exponential ($`m_{\mathrm{F547m}}=11.56`$ mag, $`r_s=11\stackrel{}{\mathrm{.}}74`$), overall the first combination gives formally a better fit. We list their parameters in Table 3. Figure 14b shows a fit where we use a single Nuker component. Although most of the bulge is removed, considerable over- and under-subtractions remain.

Various trials suggest that there is a compact source at the center of the nucleus, but there is considerable ambiguity as to whether it is resolved. While we can force the object to be a true point source in the fit, there is a small but significant pedestal in the residuals. The compact source is resolved in an unconstrained fit. However, we cannot rule it out as being unresolved because of its faintness and because the excess pedestal may be caused by model mismatch near the galaxy core. If a point source is present, we can limit its magnitude to $`m_{\mathrm{F547M}}19.8`$ mag, or $`M_V10.8`$ mag.

Another interesting component that remains in the residuals is a luminous stellar disk at the center of the dust disk. The stellar disk is flat and has sharp, truncated edges of 100 pc ($`1\stackrel{}{\mathrm{.}}6`$) in semi-major axis. The sharpness of the disk truncation, the uniformity of the disk, and the clean subtraction around its immediate surroundings all suggest the disk is likely to be real, rather than an image artifact created by a model mismatch.

#### 5.2.2 NGC 4450

The overall fit of NGC 4450 was discussed in § 4.5.2. Surface photometry of this object (Fig. 9d) shows a double power law with a break at $`r0\stackrel{}{\mathrm{.}}7`$. The apparent power-law behavior of the inner region is only slightly violated by an inflection at $`r=0\stackrel{}{\mathrm{.}}2`$, which hints at the presence of a compact central source. However, Figure 9c shows that the inflection in the 1-D profile may be caused by the extensive dust present at small radii. The situation is much clearer in 2-D. Our modeling finds a nuclear point source with $`m_V=18.9\pm 0.2`$ mag, or $`M_V=12.2\pm 0.2`$ mag.

## 6 SUMMARY

We have presented a general algorithm to decompose a galaxy into components in two dimensions. To illustrate its merits and flexibility, we have applied it to high-resolution optical HST images of 11 galaxies of various shapes and morphological types. As another example, Peng (2002) presents a detailed decomposition of the double nucleus of M31. We show that we can model the central regions of galaxies accurately, even for difficult cases that show large isophotal twists and changes in shape. These morphological complexities are often signatures of distinct galaxy components. After the major components are removed, in a number of systems we discover evidence of galaxy substructures too subtle to be seen in the original image. These include features such as nuclear point sources, low-level dust patterns, stellar disks, stellar bars, and other distinctive large-scale components. Despite the featureless appearance of some giant elliptical galaxies (e.g., NGC 4621, NGC 5982), in detail they reveal unusual shapes and slight misalignment of the sub-components. The physical interpretation of these subtle features is not yet clear. In theory, the amount of distortion in a galaxy — as evidenced, for example, by the need for fitting multiple components, their different shapes and displacements — may give clues to its evolutionary history. Substructure may also signal triaxiality in the bulge potential. A fruitful avenue of future research would be to compare the structural decomposition of real galaxies with a similar analysis applied to high-resolution N-body simulations of galaxy formation. Another is to couple the 2-D structural analysis with integral-field kinematic maps.

We thank Alice Quillen for enlightening discussions about various galaxies in our sample, and Dennis Zaritsky, Daniel Eisenstein, and Adam Burrows for discussions about error analysis. We also thank the referee for comments to improve this work. L. C. H. acknowledges financial support through NASA grants from the Space Telescope Science Institute (operated by AURA, Inc., under NASA contract NAS5-26555). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
```
===============================================================================
# IMAGE PARAMETERS
A) f547m.fits          # Input data image (FITS file)
B) f547m-out.fits      # Name for the output image
C) none                # Noise image name (made from data if blank or "none")
D) f547m-psf.fits      # Input PSF image for convolution (FITS file)
F) dust                # Pixel mask (ASCII file or FITS file with non-0 values)
G) 47 800 57 800       # Image region to fit (xmin xmax ymin ymax)
H) 427 435             # Convolution box center
I) 60   60             # Size of convolution box (x y)
J) 21.689              # Magnitude photometric zeropoint
K) 0.046   0.046       # Plate scale (dx dy).  Relevant only for Nuker model.

# INITIAL FITTING PARAMETERS
#
#   For object type, allowed functions are: sersic, nuker,
#                       expdisk, devauc, moffat, gaussian.
#

# Objtype:      Fit?         Parameters

===============================================================================

 0) nuker              # Object type
 1) 427.17  435.46 0 1 # position x, y
 3) 12.16      1       #    mu(Rb)
 4) 4.45       1       #     Rb
 5) 1.62       1       #    alpha
 6) 1.26       1       #     beta
 7) 0.21       1       #    gamma
 8) 0.72       1       # axis ratio (b/a)
 9) 39.48      1       # position angle (PA)
10) 0.00       0       # diskiness/boxiness
 Z) 0                  # Output image type (0 = residual, 1 = Don’t subtract)

 0) sersic             # Object type
 1) 427.17  435.46 1 1 # position x, y
 3) 15.22      1       # total magnitude
 4) 39.70      1       #     R_e
 5) 3.03       1       # sersic index (deVauc=4)
 8) 0.9        1       # axis ratio (b/a)
 9) -18.75     1       # position angle (PA)
10) -1.70      1       # diskiness/boxiness
 Z) 0                  # Output image type (0 = residual, 1 = Don’t subtract)

 0) expdisk            # Object type
 1) 427.17  435.46 0 0 # position x, y
 3) 14.96      1       # total magnitude
 4) 63.76      1       #     Rs
 8) 0.75       1       # axis ratio (b/a)
 9) -34.03     1       # position angle (PA)
10) 0.30       1       # diskiness/boxiness
 Z) 0                  # Output image type (0 = residual, 1 = Don’t subtract)

 0) gaussian           # Object type
 1) 427.17  435.46 1 1 # position x, y
 3) 20.50      1       # magnitude
 4) 0.50       0       #    FWHM
 8) 1.0        1       # axis ratio (b/a)
 9) -34.03     1       # position angle (PA)
10) 0.00       1       # diskiness/boxiness)
 Z) 0                  # Output image type (0 = residual, 1 = Don’t subtract)

 0) sky
 1) 0.60       1       # sky background
 Z) 0                  # Output image type (0 = residual, 1 = Don’t subtract)

```


format: text

text:
# References

PACS numbers: 36.40.Qv, 71.15.Nc, 71.15.Pd.

Metastable quasi-one-dimensional ensembles of nitrogen clusters N8

N. N. Degtyarenko, V. F. Elesin, L. A. Openov, A. I. Podlivaev

Moscow Engineering Physics Institute (State University),

Kashirskoe sh. 31, Moscow 115409, Russia

> By means of ab initio and tight-binding calculations it is shown that metastable nitrogen clusters N8 (boats) can form quasi-one-dimensional ensembles in which the nearest clusters N8 are bound to each other by covalent bonds. Those ensembles are characterized by rather high energy barriers ($`0.3`$ eV) that prevent the fission of the ensembles into isolated N8 clusters and/or N2 molecules.

A nitrogen molecule N2 is known to be the most stable form of nitrogen at standard conditions. Contrary to the overwhelming majority of other chemical elements, the cohesive energy per atom, $`E_{coh}(n)=E_{tot}(n)E_1`$ (where $`E_{tot}(n)`$ is the total energy of an $`n`$-atom cluster Nn, $`E_1`$ is the energy of an isolated nitrogen atom) decreases with $`n`$, so it is energetically favorable for clusters Nn with $`n>2`$ to decompose into separated N2 molecules and thus to lower their energy. However, recently it was shown theoretically  that several nitrogen clusters, e.g., N4 and N8, can exist in metastable states which correspond to the local minima of the potential energy as a function of atomic coordinates. Although metastable clusters Nn with $`n4`$ have not been synthesized yet, they attract much theoretical attention as potential high energy density materials (HEDM). The energy $`E_{acc}(n)=E_{coh}(2)E_{coh}(n)`$ accumulated in the Nn ($`n4`$) systems approaches the value of $`E_{acc}3`$ eV/atom .

In order to make use of Nn ($`n4`$) clusters as new efficient and environmentally safe energy sources, one should be able to create bulk metastable nitrogen structures composed of Nn ($`n4`$) clusters. On the one hand, in order such structures could exist, the nitrogen clusters should be connected with each other by covalent bonds. On the other hand, those bonds should not be ”too strong” in order the clusters preserved their individuality and the energy accumulated in such structures was not substantially lower than the energy accumulated in an isolated cluster, i.e., the clusters should form ensembles. Besides, there should be a sufficiently high energy barrier preventing the fission of ensembles into isolated nitrogen clusters and/or nitrogen molecules.

An extensive theoretical study of Mailhiot et al.  suggests that nitrogen can form bulk polymeric metastable structures at high pressures. However, the polymeric structures cannot be considered as composed of nitrogen clusters, i.e., as ensembles of clusters. In this paper, we present the preliminary results of numerical simulations of bonding between nitrogen clusters and show that boats N8 can form metastable quasi-one-dimensional ensembles in which the nearest clusters N8 are bound to each other by covalent bonds.

The version of the self-consistent-field theory with 3-21G, 6-31G, and 6-311G∗∗ basis sets in HyperChem was used in the present work. First, we studied the geometry and energetics of isolated clusters N4 (tetrahedron) and N8 (cubane and boat). Those clusters are shown in Figs. 1-3. The accumulated energies are $`E_{acc}(4)=(1.98÷2.21)`$ eV/atom for the cluster N4, $`E_{acc}(8)=(2.44÷2.59)`$ eV/atom for the cubane N8, and $`E_{acc}(8)=(1.53÷1.66)`$ eV/atom for the boat N8, depending on the basis set used. Our results agree with previous works .

We have carefully analyzed a lot of different mutual arrangements of nitrogen clusters. In both cases of N4 and cubane N8 clusters, we have found no metastable configurations. However, in the case of boats N8, we have revealed that two boats can form a ”cluster molecule” (N8)2. It is shown in Fig. 4. The energy accumulated in this molecule, $`E_{acc}(16)=1.79`$ eV/atom (3-21G basis set), exceeds slightly the energy accumulated in an isolated boat N8. We have demonstrated by direct calculations that the molecule (N8)2 is mechanically stable, i.e., all its oscillation frequencies are real, giving evidence that the configuration shown in Fig. 4 corresponds to a local minimum in the total-energy surface as a function of atomic coordinates.

In order to study the formation of large ensembles of boats N8, we employed a version of TBMD (tight binding molecular dynamics) method with transferable interatomic potential  used by us earlier for carbon clusters  and modified for calculations on nitrogen systems. Although the TBMD method is not as rigorous as ab initio approaches, it allows one to study relatively large systems as well as to calculate directly the lifetime of the metastable state and the height of the energy barrier. For N4 and N8 clusters and (N8)2 ”boat dimers”, the results of TBMD simulations agree qualitatively with ab initio calculations. We have also shown that rather long quasi-one-dimensional ensembles of boats can be formed through attaching of new boats N8 to the (N8)2 ”molecule”. The ensembles (N8)m with $`m`$ up to 6 have been shown to be metastable with accumulated energy about 90 percents of the energy $`E_{acc}(8)`$ accumulated in an isolated boat. The heights of the energy barriers preventing the fission of metastable configurations exceed 0.3 eV.

Finally, we note a recent success in recovering the non-molecular phase of nitrogen at ambient pressure at temperatures below 100 K . Structure and physical properties of a new phase are not understood yet. Its characterization remains an important issue. Presumably, it is polymeric. This can be considered as a first step towards synthesis of bulk metastable nitrogen structures, including ”cluster phases” studied in this paper.

Acknowledgments

The work was supported by the CRDF (project B-P/252), the Contract DSWA01-98-C-0001, and the Russian Federal Program ”Integration” (Project AO133).

Figure captions

Fig.1. Tetrahedron N4.

Fig.2. Cubane N8.

Fig.3. Boat N8.

Fig.4. ”Cluster molecule” (N8)2 composed of two boats N8.