DETAILED DECOMPOSITION OF GALAXY IMAGES. II. BEYOND AXISYMMETRIC MODELS

Analyzing astronomical images is challenging because of the diversity in object sizes and shapes, and nowhere is it more difficult than for galaxies. Since the early era of photographic plates, one of the key methods for studying the light distribution of galaxies is to model it by using analytic functions—a technique known as parametric fitting. This technique was first applied to galaxies by de Vaucouleurs (1948) who showed that the light distribution of elliptical galaxies tended to follow a power-law form of exp(−r^1/4). Subsequently, one of the breakthroughs in our understanding of galaxy structure and evolution came when Freeman (1970) showed that dynamically "hot" stars in galaxies make up spheroidal bulges having a de Vaucouleurs light profile, whereas "cold" stellar components make up the more flattened, rotationally supported, exponential disk region.

From that simple beginning, parametric fitting has been the mainstay for galaxy imaging studies, and expanding into many applications whenever the science calls for detailed and rigorous analysis. Among some of the examples, past investigations delved into the structural parameters of disk galaxies (e.g., de Jong 1996), the Tully–Fisher relation (e.g., Tully & Fisher 1977; Hinz et al. 2003; Bedregal et al. 2006), the evolution of disky galaxies (Simard et al. 2002; Ravindranath et al. 2004; Barden et al. 2005), the cosmic evolution of galaxy morphology (e.g., Lilly et al. 1998; Marleau & Simard 1998; Hathi et al. 2009) in both ground-based surveys and Hubble (Ultra-)Deep Fields (Williams et al. 1996; Beckwith et al. 2006), the morphological transformation of galaxies in cluster environments (e.g., Dressler 1980), the fundamental plane of spheroids (Djorgovski & Davis 1987; Dressler et al. 1987; Bender et al. 1992), the red sequence of galaxies (Bell et al. 2004a, 2004b; Faber et al. 2007), morphological dissimilarities between spheroidal galaxies and ellipticals (Kormendy 1985, 1987), the central structure of early-type galaxies (Kormendy 1985; Lauer et al. 1995, 2007; Faber et al. 1997; Ferrarese et al. 2006a, 2006b) and implications for the formation of massive black holes (Ravindranath et al. 2002; Kormendy & Bender 2009), black hole versus bulge relations (Kormendy & Richstone 1995) and their evolution (Rix et al. 2001; Peng et al. 2006a, 2006b), the "extra light" due to gas dissipation in galaxy centers (Kormendy 1999; Cote et al. 2006, 2007; Kormendy et al. 2009; Hopkins et al. 2008, 2009), quasar host galaxies (e.g., Hutchings et al. 1984; McLeod & Rieke 1994; McLure et al. 2000; Jahnke et al. 2004; Sánchez et al. 2004; Kim et al. 2008b), gravitational lensing of quasar host galaxies (Rix et al. 2001; Peng et al. 2006b), and the clustering of dark matter through weak lensing (e.g., Heymans et al. 2006, 2008).

Since the original development of galaxy fitting nearly 70 years ago, where the analysis was performed on one-dimensional (1D) surface brightness profiles (see also Kormendy 1977; Burstein 1979; Boroson 1981; Kent 1985; Andredakis & Sanders 1994; MacArthur et al. 2003), newer techniques have emerged to directly analyze two-dimensional (2D) images (e.g., Shaw & Gilmore 1989; McLeod & Rieke 1994; Byun & Freeman 1995; de Jong 1996; Moriondo et al. 1998; Simard 1998; Ratnatunga et al. 1999; Wadadekar et al. 1999; Simard et al. 2002; de Souza et al. 2004; Laurikainen et al. 2004; Gadotti 2008). The benefit of performing 2D image analysis is to potentially make full use of all spatial information and to properly account for image smearing by the point-spread function (PSF).

Even though 2D analysis can be quite sophisticated, there are legitimate questions about whether it is more beneficial than 1D for profile analysis. Proponents of the 1D technique are skeptical that perfect ellipsoid models are suitable to use for galaxies that show isophotal twists, or that are non-elliptical in shape. They note that not only is 1D analysis more appealing because it is more straightforward to implement, but also the surface brightness profiles serve as visual confirmation about the reality of fitting multiple-component models.

However, beneath the apparent simplicity there are a number of important subtleties to weigh. For instance, the decision about how to extract 1D profiles is often not unique, nor are there strong reasons to prefer major or minor axis profiles, or a profile along some arc traced by spiral arms or isophote twists that result from the superposition of multiple components oriented at different angles. When symmetry is broken, it is also unclear that there is an optimal or unique way to extract a 1D profile, such as in irregular galaxies, overlapping galaxies, and galaxies with double nuclei. Another factor to consider is that the process of extracting 1D profiles reduces spatial information content: in many situations, a bulge, disk, and bar can all have different axis ratios, position angles (P.A.s), and profiles that help to break model degeneracies, but this information is lost when the data are collapsed into 1D. Lastly, for compact galaxies, 1D profile fitting cannot properly correct for image smearing by the PSF because 1D profile convolution is not mathematically equivalent to convolution in a 2D image. While some of the above concerns also affect 2D analysis (i.e., irregular galaxies), most others benefit from treatment using 2D techniques. When it comes to judging which models are more plausible, there are few diagnostics more discerning than a moment's glance at 2D models and residual images; a good fit in 2D always means that 1D profiles are necessarily a good fit. Proponents of 2D analysis therefore believe that the benefits outweigh the drawbacks. Moreover, many drawbacks can be mitigated by breaking free from axisymmetry in 2D analysis, which we aim to show in this study.

In the box of tools for morphology analysis, a complimentary approach is non-parametric analysis. While we do not use non-parametric methods in this study, it is useful to understand the conceptual differences between the two approaches. We thus provide a brief overview. In contrast to function fitting, the non-parametric approach does not involve deciding what functional form to use or how many. One method is to decompose an image into "shapelets" or "wavelets" (e.g., Refregier 2003; Massey et al. 2004), which is analogous to taking a 2D Fourier transform of an image using mathematically orthogonal basis functions. The main conceptual difference with parametric fitting is that the shapelet basis functions do not represent physical subcomponents of a galaxy. Moreover, the power spectrum of the basis functions is quite useful for diagnosing the degree of galaxy distortions. There are also other non-parametric techniques (e.g., Abraham et al. 1994; Rudnick & Rix 1998; Conselice et al. 2000; Lotz et al. 2004). To measure concentration non-parametrically, one way is to compare fluxes within apertures of different radii; whereas to measure asymmetry one can rotate an image by 180° and subtract it from the original image and measure the residuals (e.g., Abraham et al. 1994; Conselice et al. 2000). Toward the same goals, two studies, Abraham et al. (2003) and Lotz et al. (2004), introduce the Gini index to measure the concentration of a galaxy by comparing the relative distribution of pixel flux values within a certain area. Lotz et al. (2004) also introduce a method for measuring asymmetry through the M₂₀ parameter, which is the second-order moment of the brightest 20% of a galaxy's flux.

The application of non-parametric analysis has mostly been to quantify galaxy mergers (e.g., Conselice et al. 2003; Lotz et al. 2008). These techniques are generally much simpler to implement than parametric fitting and have a strong virtue that no assumptions are made about the galaxy profiles and shapes. The trade-off is that the techniques often do not deal with image smearing by the PSF and different sensitivity thresholds between different surveys. Consequently, one has to take particular care to compare compact with extended objects, measured in different apertures, or measured from images of different surface brightness depths (Lisker 2008). One also should guard against contamination by intervening galaxies or stars because the techniques do not have a rigorous way to separate overlapping objects. For separating objects, extracting structural components of a galaxy, and extrapolating galaxy wings well into the background noise, there are few, if any, alternatives to parametric analysis that are more rigorous.

When comparing the merits of non-parametric and parametric analysis, the idea of ellipsoid models in parametric analysis is sometimes considered to be a weakness, because galaxies, after all, are not perfect ellipses in projection. However, it is worth pointing out that the notion of there being a global average size inherently implies comparison against some kind of approximate shape. Indeed, even in non-parametric techniques, to measure a size in a 2D image, one assumes a basic shape either explicitly (through using aperture photometry) or implicitly (through calculating flux moments, which requires a center to be defined). An ellipsoid is one of the simplest and most natural low-order shapes against which all galaxies can be compared, especially for measuring an average size. This notion is useful: deviations from the basic ellipsoid shape can then be considered as higher order modifications, even for highly irregular galaxies.

Nevertheless, there are many situations where it is desirable to use models that deviate from ellipsoid shapes. Contrary to the common practice, there are numerous ways to break from axisymmetry. However, the harder challenge is to devise a scheme that is intuitive to grasp and well motivated. Breaking free from axisymmetry allows for other interesting science applications, including a promising new way to quantify asymmetry.

In this study, we present, as a proof of concept, new capabilities in 2D image fitting that progress beyond the limitations of traditional parametric fitting models. One key aspect of our approach is to first identify a minimum basis set of features that spans the range of galaxy morphologies and shapes. From experience, we determine those four "basic elements" to be bending, Fourier, coordinate rotation, and truncation modes. Second, one of the main reasons why parametric fitting is useful is that the profile parameters are intuitive to grasp (e.g., concentration index, effective radius, total luminosity, etc.). Therefore, another key requirement is that the traditional profile parameters must retain their original, intuitive, meaning even under detailed shape refinements, and even under such extreme cases as irregular galaxies. This can be accomplished if the basic premise starts with the traditional ellipsoid function, on top of which one can add perturbations, rotations, irregularities, and curvature. This is possible because of the fact that simple ellipsoid fits are a reasonable way to quantify global average properties, and other details can be considered to be higher order perturbations that may be of other practical interest.

As we attempt to demonstrate, combining just the four basic morphology elements can quickly yield a dizzying array of possibilities for fitting galaxies. The end result can look highly "realistic." Indeed, it is now possible to fit many spiral galaxies, asymmetric tidal features, irregular galaxies, ring galaxies, dust lanes, truncated galaxies, arcs, among others (though, certainly, there are limitations). However, it is important to realize that "being possible" often does not mean "being necessary" or "being practical." Necessity ought to be judged in the scientific context of whether it is worth the extra effort to obtain diminishing returns. For instance, to measure total galaxy luminosity, it is often unnecessary to fit high-order Fourier modes or spiral rotations. For many science studies interested in global parameterization, often a single ellipsoid component would suffice. It is therefore important to always let the science determine what kind of analysis is required, rather than to use the new capabilities in the absence of a clearly defined goal. Having provided some foregoing disclaimers, some of the key scientific reasons motivating the new capabilities are to: (1) quantify global asymmetry or substructure asymmetry; (2) quantify bending modes for weak-lensing applications, or fit arcs in the image plane for strong gravitational lenses; (3) obtain more accurate substructure decomposition in the presence of bars, spirals, rings, etc.; (4) obtain more accurate global photometry; (5) quantify profile deviation in inner or outer regions of a galaxy, such as disk truncations, deviations from a Sérsic function, etc.; (6) Extract parametric information to the limits imposed by resolution, signal-to-noise ratio (S/N), and other small-scale fluctuations; and (7) quantify model-dependent errors in the decomposition.

We thus begin by giving an overview of the Galfit software (Section 2). Then, we introduce the radial profile functions that one can use (Section 3), and illustrate how symmetric and asymmetric shapes can be generated by modifying the coordinate system in Section 4. Next, we introduce a new capability that allows for radial profile truncation (Section 5). Enabling all the capabilities may result in extremely complex galaxy shapes, the interpretation of which may give concerns to those new to the analysis. Therefore, we discuss the interpretations and model degeneracies of the parameters in Section 6. We then apply these new features to real galaxies in Section 7, followed by concluding remarks (Section 8).

Galfit is a nonlinear least-squares fitting algorithm that uses the Levenberg–Marquardt technique to find the optimum solution to a fit. The Levenberg–Marquardt algorithm is currently the most efficient one for searching large parameter spaces, allowing for the possibility of fitting complex images with multiple components and a large number of parameters. Galfit determines the goodness of fit by calculating χ² and computing how to adjust the parameters for the next step. It continues to iterate until the χ² no longer decreases appreciably. The indicator of goodness of fit is the normalized or reduced χ², χ²_ν:

$$\begin{equation} \chi ^2_\nu = \frac{1}{N_{\rm dof}}\sum _{x=1}^{nx}\sum _{y=1}^{ny} \frac{\left(f_{\rm data}(x,y) - f_{\rm model}(x,y)\right)^2}{{\sigma (x,y)}^2}, \end{equation} \tag{ 1 }$$

where

$$\begin{equation} f_{\rm model} (x,y) = \sum _{\nu =1}^{m} f_{\nu }(x, y; \alpha _1,\ldots,\alpha _n). \end{equation} \tag{ 2 }$$

Here N_dof is the number of degrees of freedom in the fit; nx and ny are the x and y image dimensions; and f_data(x, y) is the image flux at pixel (x, y). The f_model(x, y) is the sum of m functions of f_ν(x, y; α₁, ..., α_n), with n free parameters (α₁, ..., α_n) in the 2D model. The uncertainty as a function of pixel position, σ(x, y), is the Poisson error at each pixel, which can be provided as an input image. If no σ-image is given, one is generated based on the gain and read-noise parameters contained in the image header. Pixels in the image marked as being bad do not enter into the calculation of χ².

In the Levenberg–Marquardt algorithm, the minimization process involves computing a Hessian matrix, which is closely related to the covariance matrix of the parameters (e.g., see Press et al. 1992). The covariance matrix is then directly related to the formal uncertainty in the fitting parameters that Galfit reports. However, the usefulness of the formal uncertainty is limited to ideal situations where the fluctuations in the residual image are only due to Poisson noise after removing the model. This situation is mostly realized in idealized situations, such as image simulations. In real images, the residuals are due to structures like stars and galaxies that are not fitted, flat-fielding errors, and imperfect functional match to the data. These factors cause formal uncertainties reported in numerical fits to be only lower estimates. In image fitting, more realistic uncertainties are necessarily obtained by other processes, such as comparing fit results based on different assumptions about the model rather than through a formal covariance matrix.

In summary, the three images Galfit takes as input to calculate least squares are the data, a σ-image, and an optional bad pixel mask. To account for image smearing by the PSF, Galfit will also require a PSF image.

The wavefront of light from distant sources is always perturbed by the act of producing an image, distortions due to imperfect optics, and sometimes by Earth's atmosphere, resulting in some blurring. To accurately compare the intrinsic shape of an object with a model, image blurring must be taken into account. In image fitting, this is often done by convolving a model image with the input PSF before comparing with the data. The process of performing convolution is mathematically rigorous, but the actual implementation has several subtleties.

One consideration is the computation speed, as the process of convolving a model is frequently the most time-consuming part of parametric fitting. The trade-off is that the smaller the convolution region, the faster the computation time, but also the less accurate. To achieve a compromise, Galfit allows the user to decide on the size of the convolution region. This gives the flexibility for one to hone in on a solution quickly before trying to obtain higher accuracy in the final step.

Another important issue to consider is whether to convolve each component separately or all of them together just once in the final image. This is an important consideration because even though the model functions are analytic, they are resampled by discrete pixel grids, resulting in a "pixellated" profile instead of one that is infinitely smooth. If an intrinsic model is sufficiently sharp, the curvature may not be critically subsampled by the pixels prior to convolution, regardless of whether the recorded data are Nyquist sampled. The resulting profile after image convolution therefore can depend very sensitively on how the model is centered on a pixel. If such a model is created off-center, pixellation effectively broadens out the model ever so slightly more than normal once convolution is applied, but the effect is noticeable in high-contrast imaging studies. Therefore, the better way to deal with "pixellation broadening" is to convolve each model component individually rather than the entire image at once. To do so, Galfit creates every model on a pixel center; the pixel fluxes near the center of the models are integrated over the pixel area adaptively. Then to effect an off-centered model, Galfit makes use of the convolution theorem by shifting the PSF by the required amount before convolving it with the model. This process circumvents the problem of artificial pixellation broadening because whereas the model core region may not be sufficiently resolved, the PSF ought to be.⁶

Shifting the PSF, however, can be quite problematic when it is marginally Nyquist sampled, or if the diffraction patterns are not critically sampled. Accurate shifting of the PSF is of basic importance in high-contrast imaging studies. For instance, in the case of studying active galaxies with a strong central point source, issues of contrast, resolution, and sampling all conspire to make the PSF fitting crucial to deriving a reliable host model. In such situations, the standard interpolation techniques (e.g., linear or spline) tend to broaden out the PSF core, so they are only accurate in the extended outskirts where the gradient is shallow. One alternative is to interpolate using the sinc kernel, which is theoretically the perfect interpolation kernel for critically sampled images, and preserves the intrinsic width of the data. However, significant "ringing" appears around sharp features (i.e., PSF or galaxy core). This effect can be nearly as bad on a fit as pixellation broadening. An improved solution is to taper the wing of the sinc kernel using a windowing function (e.g., Lanczos), but the ringing often may still be quite large beyond the PSF core, which must be further suppressed.

Galfit seeks a compromise by using a hybrid scheme where the interpolation in the PSF core is done by using a sinc kernel with a Kaiser window function so as to faithfully preserve the width, but a bicubic spline interpolation is used in the wings. The result of this scheme is that for a Gaussian having an FWHM of 2 pixels, the interpolation is accurate to 0.1% in the center, and 0.03% at the distance of the FWHM relative to the peak (or 1% relative to the local flux). For oversampled PSFs, the interpolation is even more accurate. Compared to bicubic spline interpolation, our scheme is about 20 times more accurate. From a more practical standpoint, the mismatch in the PSFs between data taken using the Hubble Space Telescope (HST) imaging cameras and synchronously observed PSFs is rarely better than 3% in the core. For all practical purposes, our interpolation scheme therefore will more than suffice for the most demanding high-contrast studies of quasar host galaxies at high redshift.

When the data are undersampled, convolution of the model can still be done correctly if the convolution PSF provided to Galfit is either critically sampled or oversampled. In this situation, Galfit will generate a model on a finer grid, convolve it with the PSF, then bin the result down to the resolution of the data for comparison. One way for users to obtain an oversampled PSF compared to the data is to dither the PSF observations by fractional pixels. Another way is to numerically reconstruct a more oversampled PSF star by extracting multiple stars from the data image itself (e.g., via DAOphot, Stetson 1987).

However, lastly, we note that when the data and the convolution PSF are both undersampled (i.e., with PSF FWHM <2 pixels), convolution cannot be done accurately. In such a situation, for the purpose of image fitting, it is often better to broaden out the data and the PSF to critical sampling than to perform the analysis in the original resolution (Kim et al. 2008a).

Using the new features, each model can take on a shape that is completely unrecognizable from a traditional ellipsoid shape. It is therefore necessary to clarify what constitutes a single model component. In Galfit, each model component is referred to by the name of the surface brightness profile, just as it is standard practice to call something a Sérsic, Gaussian, or exponential component in traditional models. As implied by this notion, no matter how complex the shape, the flux declines monotonically (unless modified by a truncation function, Section 5) from a peak in every radial direction in a non-rotating frame, or along an arc in a rotating frame, strictly following the functional form specified by the user. The radial profile parameters are mathematically decoupled from the azimuthal shape because the radial profile functions are self-similar in the expression of the radius parameter, i.e., with powers of (r/r_e), whereas the complex azimuthal shapes are obtained by simply stretching the coordinate metric into more exotic grids than the standard Cartesian grid. This idea is in fact implicit in all 2D image-fitting algorithms, where the axis ratio parameter, q, turns a circular profile into an ellipse by compressing the coordinate axis along one direction, even though the functional form of the profile remains the same in every direction. In the same manner, the definition of a scale or effective radius in a component, no matter how complex the shape, corresponds closely to that of the best-fitting ellipse in the direction of the semimajor axis.

Figure 1 demonstrates how Sérsic profiles can be modified by bending modes, Fourier modes, and a spiral rotation function in Galfit—the results of which are shown in Figure 2. In the example, there are only two Sérsic model components, despite the appearance of numerous parameters: both the "radial" and "power" functions are modifications to the Sérsic profiles. Furthermore, the Fourier and bending modes can modify the Sérsic profiles, or modify the modifiers to the light profiles. Each radial surface brightness profile has a single peak, and the flux decline is monotonic radially (unless truncated by a truncation function called "radial" in Figure 1) or in a rotating coordinate system (called "power" in Figure 1). Therefore, for each component, it is still meaningful to talk about, for example, an "average" light profile (e.g., Sérsic), with an average Sérsic concentration index n—no matter what the galaxy may look like azimuthally. In this manner, even irregular galaxies can be parameterized in terms of their average light profile. When the average peak of an irregular galaxy is not located at the geometric center, it has a high-amplitude m = 1 Fourier mode (i.e., lopsidedness).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In such a way, no matter how complex the azimuthal shape, interpreting the surface brightness profile parameters is just as straightforward as the traditional ellipsoid.

In Galfit, each truncation function can modify one or more light profile models. Also, any number of light profiles can share the same truncation function. The truncation function in Galfit is a hyperbolic tangent function (see Equation (B2) in Appendix B). Schematically, a truncated component is created by multiplying a radial light profile function, f_0,i(x, y; ...), by one or more truncation functions, P_m or 1 − P_n (depending on whether the type is an inner or an outer truncation), in the following way:

$$\begin{eqnarray} &&f_i (x,y;\ldots) = f_{0,i} (x, y; x_{c,i}, y_{c,i,}\ldots,\ q_i, \theta _{{\rm P.A.}, i})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \nonumber\\ &&\times \prod ^m P_m (x, y; x_{c, m}, y_{c, m}, r_{{\rm break},m}, \Delta r_{{\rm soft},m}, q_m, \theta _{{\rm P.A.},m}) \nonumber\\ &&\times \prod ^n [1-P_n (x, y; x_{c, n}, y_{c, n}, r_{{\rm break},n}, \Delta r_{{\rm soft},n}, q_n, \theta _{{\rm P.A.},n})].\nonumber\\ \end{eqnarray} \tag{ 30 }$$

The break radius, r_break, is defined to be the location where the profile is 99% of the original (i.e., untruncated) model flux at that radius. The parameter Δr_soft is the softening length, so that r = r_break ± Δr_soft is where the flux drops to 1% of that of an untruncated model at the same radius (the ± sign depends on whether the truncation is inner or outer). The inner truncation function (P_m) tapers a light profile in the inner regions of a light profile, whereas the outer truncation function (1 − P_n) tapers a light profile in the wings.

The behavior of the hyperbolic tangent function is ideal for truncation because it asymptotes to 1 at the break radius r ≳ r_break and 0 at the softening radius r < r_soft, and vice versa for the complement function. Thus, when multiplied to a light profile f(r), the functional behavior exterior to the break radius has intuitively obvious meanings. For example, as shown in Figure 16(a), if a Sérsic function with n = 4 is truncated in the wings (shown in red), the core has exactly an n = 4 profile interior to r_break (marked with a vertical dashed line), which is a free parameter to fit. Likewise, an n = 4 profile truncated in the core (green) has exactly an n = 4 profile exterior to the outer break radius. Thus, when one sums two functions of different Sérsic indices n (Figure 16(b)) the asymptotic profiles of the wing and core retain their original meaning, and there is very little cross talk outside of the truncation region (denoted by vertical dashed lines in Figure 16).

Zoom In Zoom Out Reset image size

Use of the truncation functions is highly flexible. There can be an unrestricted number of inner and outer truncation functions for each light profile model. Furthermore, multiple light profile models can share in the same truncation functions. This is useful, for instance, when trying to fit a dust lane (inner truncation) in a fairly edge-on galaxy that may affect both the bulge and the disk components. Just as with light profile models, the truncation functions can be modified by Fourier modes, bending modes, etc., independent of the higher order modes for the light profile they are modifying.

Truncation models appear in many physical contexts, such as dust lanes, rings, spirals that do not reach the center, joining a spiral with a bar, or cutoff of the outer disk. To allow the truncation parameters to be more intuitive to understand given situations at hand, Galfit offers several variations. In addition to inner and outer truncations, truncation functions can share in the same parameters as the parent light profile. There are radial and length/height truncations, softening radius versus softening length (default versus type 2), inclined versus non-inclined (default versus type b) truncations, and, lastly, four different ways to normalize the flux—the most sensible choice depends on how a profile is truncated. We now discuss each of these variations in more detail.

Parameter sharing. In the most general form, each truncation function has its own set of free parameters: x₀, y₀, r_break, Δr_soft, q, and θ_P.A.. However, by default, the parameters x₀, y₀, q, and θ_P.A. are tied to the light profile model, and are activated only when the user explicitly specifies a value for them.

Radial ("radial") versus length ("length") or height ("height") truncations. The most useful type of truncation is one that has radial symmetry to first order, i.e., where it has a center, an ellipticity, and an axis ratio. However, in the case of a perfectly edge-on disk galaxy ("edgedisk" model), an additional type is allowed that truncates linearly in length or in height. For instance, a dust lane running through the length of the galaxy has an inner height truncation. For the "edgedisk" profile, Galfit also allows for a radial truncation, as with all other functions. The one drawback to height and length truncations is that they cannot be modified by Fourier and higher order modes like the radial truncations.

Softening length ("radial") versus softening radius ("radial2"). Sometimes, instead of softening length (Δr_soft), it is more useful for the fit parameter to be a softening radius (r_soft), especially when one desires to hold the parameter fixed. That is also allowed in Galfit as a type 2 truncation function, designated, for example, as "radial2." The default option does not have a numerical suffix.

Inclined (default, "radial") versus non-inclined ("radial-b") truncations.) A spiral rotation function is an infinitesimally thin, planar structure. Nevertheless, it should be thought of as a 3D structure in the sense that the plane of the spiral can be rotated through three Euler angles, not just in P.A. on the sky. When a truncation function is modifying a spiral model, it is therefore sometimes useful to think about the truncation in the plane of the spiral model. When Fourier modes and radial truncations are modifying a spiral structure, the default ("radial") is for the modification to take place in the plane of the spiral structure. However, there are some instances when that may not be ideal (e.g., a face-on spiral may actually be ellipsoidal). In those situations, one can choose "radial-b," which would allow a truncation function to modify the spiral structure in the plane of the sky, even though the spiral structure can tip and tilt as needed.

Lastly, the truncation function can be type 2b (i.e., "radial2-b") as well.

Flux normalization. The most intuitive flux normalization for a truncated profile is the total luminosity. Unfortunately, both the total luminosity and the derivative of the free parameters with respect to the total luminosity are especially time consuming to work out computationally and slow down the iteration process. There are generally no closed form analytic solutions to the problem. Therefore, the alternative is to allow for different ways to normalize a component flux. The user may choose whichever one is more sensible, given the situation and the science task at hand. The default depends on the truncation type.

• 1.  
  Inner truncation: the flux is normalized at the break radius. This is most sensible for a ring model because this radius roughly corresponds to the peak flux of the ring.
• 2.  
  Outer truncation: flux normalized at the center.
• 3.  
  Both inner and outer truncations: same as the case for inner truncation.

However, there are many situations when the default is not desirable. Instead, the user can choose the radius where the flux is normalized. To be pedagogical, we explicitly show here the normalization for just the Sérsic function:

• 1.  
  function: default (e.g., "sersic," "nuker," "king," etc.). See the details for individual functions.
• 2.  
  function1: flux normalized at the center r = 0 (i.e., Σ₀). A function that is given originally by f_orig(r) is now defined as $f_{\rm mod}(r) = \Sigma _0 \frac{f_{\rm orig}(r)}{f_{\rm orig}(0)}$. For the Sérsic profile (i.e., called "sersic1"), the profile function is redefined in the following way, written explicitly as

  $$\begin{equation} f_{\rm mod}(r)=\Sigma _0\frac{\exp {\left[-\kappa \left(\left({\frac{r}{r_e}}\right)^{1/n} - 1\right)\right]}}{\exp {\left[\kappa \right]}}. \end{equation} \tag{ 31 }$$

  For the Ferrer and King profiles, this normalization is the same as the default normalization.
• 3.  
  function2: flux parameter is the surface brightness at a model's native size parameter (parameter 4 of the light profile model). For a Sérsic profile, called "sersic2," this means the effective radius r_e. So, $f_{\rm mod}(r) = \Sigma _e\frac{f_{\rm orig}(r)}{f_{\rm orig}(r_e)}$. For example, a Sérsic profile now has the following explicit form:

  $$\begin{equation} f_{\rm mod}(r)=\Sigma _e \exp {\left[-\kappa \left(\left({\frac{r}{r_e}}\right)^{1/n} - 1\right)\right]}. \end{equation} \tag{ 32 }$$

  For the Nuker profile this normalization is the same as the default normalization.
• 4.  
  function3: flux parameter is the surface brightness (Σ_break) at the break radius (r_break). This is the most useful situation when a truncation results in a large-scale galaxy ring, so that the surface brightness parameter corresponds closely to the peak of the light profile model. When the truncation is not concentric with the light profile model, this kind of normalization is not very intuitive. For "radial" truncation, r_break is parameter 4, whereas for "radial2," r_break is parameter 4 for outer truncation and parameter 5 for inner truncation. When the "sersic3" option is chosen, the r_break parameter comes automatically from the first truncation component with which a certain light profile model is associated.In our example of the Sérsic profile, $f_{\rm mod}(r) = \Sigma _{\rm break}\frac{f_{\rm orig}(r)}{f_{\rm orig}(r_{\rm break})}$. For example, a Sérsic profile now has the following explicit form:

  $$\begin{equation} f_{\rm mod}(r)=\Sigma _{\rm break} \frac{\exp {\left[-\kappa \left(\left({\frac{r}{r_e}}\right)^{1/n} - 1\right)\right]}}{\exp {\left[-\kappa \left(\left({\frac{r_{\rm break}}{r_e}}\right)^{1/n} - 1\right)\right]}}. \end{equation} \tag{ 33 }$$

Figure 17 demonstrates just some of the possibilities allowed when fitting truncations. In addition to the regular ellipsoid shape, the higher order modes like diskiness/boxiness parameters, bending modes, and Fourier modes can also modify the shape of the truncation functions. One can also use the truncation function on a spiral model, on models with Fourier and bending modes, and diskiness/boxiness models, some of which are shown in Figures 17(d), 17(e), 16(f), and 17(i). When a truncation function acts on a spiral component, it can do so either in the plane of the disk ("type a") or in the plane of the sky ("type b;" e.g., "radial-b"). While the default is in the plane of the disk, the parameters are more intuitive in type b cases when the disk is tilted and rotated.

Zoom In Zoom Out Reset image size

The use of truncation functions should be carefully supervised because unexpected things can happen, such as the size or the concentration index of a component can grow without bound. This behavior is due to the fact that there are degeneracies between the sharpness of truncation and the steepness/size of the galaxy. Therefore, truncation functions should only be used on objects that clearly have truncations.

When two functions are joined by using a truncation function, the cross-talk region is located in between the two truncation radii: it is worth bearing in mind the definition that at the break and softening radii, the fluxes are 99% and 1% that of the same model without truncation, respectively. In other words, the larger the truncation length, the larger the cross-talk region. Therefore, when one (or more) of the parameters r_break, r_break + Δr_soft, or r_soft is either too small (≲ few pixels) or larger than the image size, it probably indicates that profile truncation parameters are not meaningful. Rather, it more likely reveals the fact that there is a mismatch between the light profile model and the actual galaxy profile.

There are well-known situations when different parameters in one or more functions are capable of modeling the same profile behavior. This scenario is the one most commonly referred to in generic discussions about model degeneracies. For instance, very large Sérsic index values (n ≳ 4) have highly extended wings, the presence of which is non-unique with the sky parameter. A high n, caused by profile mismatch or poor model prior, can often suppress the sky estimate. It is therefore advisable to estimate the sky independent of the fit, and to hold it fixed to the best estimate. As a second example, in the Nuker profile (Equation (17)), there are three parameters (α, β, and γ) that control the inner/outer slopes and sharpness of the bending (Figure 7). When the break radius r_b of a Nuker profile is sufficiently small and profile mismatch sufficiently large, model discrimination relies entirely on the power law $\frac{\gamma -\beta }{\alpha }$. Because there are numerous ways to yield a specific value for $\frac{\gamma -\beta }{\alpha }$ in the model, it leads to a degenerate situation involving three parameters.

As another example, a low-amplitude second Fourier mode and the first bending mode (shear) can both be degenerate with the axis ratio q parameter of an ellipse; therefore, they should not be used together except in obvious situations where doing so is useful. Lastly, in the spiral rotation function, the periodicity of the rotation function can sometimes be a source of "degeneracy." Multiple windings can approximate a smooth continuous model, whether or not there is a spiral structure present. For instance, a classical Sérsic ellipsoid can be simulated by a spiral model with a very large θ_out. While the fit is not good and easy to diagnose by an end user, it is nevertheless a numerically allowed solution.

The above situations are not meant to be a complete laundry list, but they are the most common situations. In complex analysis, one always needs to be circumspect about the potential hazards of mixing and matching different functions whose parameters produce similar profile behaviors. Even though Galfit allows for a great deal of flexibility in the analysis, it is ultimately up to the user to decide on what to allow, based on the goals of the science, and to understand when potentially degenerate parameters may be used effectively.

The above discussion may also seem to imply degeneracies or non-uniqueness are too numerous for complex analysis to be practical or reliable. That notion is only true when it is not possible to verify the results of a fit and to try out other solutions. Such a scenario is more common for large-scale galaxy surveys, in which automated, detailed, analysis is admittedly quite difficult to conduct sensibly. However, even in those scenarios, there are many situations where mutually coupled parameters do not affect the other main parameters of scientific interest: degeneracies in the Fourier modes often do not have any bearing on the total luminosity or size of a component. Moreover, when an analysis is done manually, it is reassuring that the problems are almost always easy to recognize and remedy when they do happen, even by simple inspection of the model and residual images.

Occasionally, what appears to be numerical degeneracy problems may be caused by someone using a code outside the algorithm's physical capabilities. As such, it is a pseudo-degeneracy. Different algorithms have different limitations that affect convergence, whether the code is gradient descent, Metropolis, or otherwise. This situation may appear like parameter degeneracy because restarting the fit does indeed yield a different solution, but in fact the code may be hamstrung in its convergence. For example, gradient descent algorithms require the calculation of a gradient, and thus can run into problems when the gradient cannot be calculated properly. In simulated annealing algorithms (e.g., Press et al. 1992), parameter boundaries and annealing speed control the algorithmic behavior: anneal too quickly, the solution may settle into a local minimum. To search larger parameter spaces requires longer annealing times.

While all algorithms have conditions under which they perform poorly, pseudo-degeneracies can always be recognized and mitigated. Galfit is based on a Levenberg–Marquardt subroutine that performs the least-squares minimization. In part a gradient descent algorithm, the convergence behavior is affected by the calculation of gradient images that determines the direction of steepest χ² descent. When the gradient images are problematic, they affect the convergence to a proper solution. There are three main problematic situations. The first, and most common, occurs when a model becomes extremely compact (FWHM ≲0.5 pixels), so that the profile gradient cannot be resolved: all the gradient information in the model fits into a single pixel. This situation mostly arises when working with high-contrast imaging data, such as quasar host galaxy decomposition, when one of the subcomponents may be used to reduce the strong residuals caused by a PSF mismatch. A similar situation arises when a model is very thin (axis ratio q ≲ 0.05) and the object is compact; here, the gradient does not exist along one spatial direction because of the lack of pixel resolution. Another rare example occurs when the inclination angle of a spiral rotation component is close to perfectly face-on (θ_incl → 0 in Equations (28) and (29)), when the derivative image for the inclination parameter approaches zero.

Another abnormal numerical behavior may occur when one places parameter constraints on a model to prevent some parameters from wandering too far from their initial values. Doing so may cause poor convergence by forcing the solution into a tight "corner," when the best solution is somewhere beyond it. A typical situation is where there are other sources in the image that are not masked or fitted by models, but that are sufficiently luminous to influence the fit of the target of interest. In this situation, no amount of effort will produce a sensible solution, because the best solution is outside of the parameter boundaries, even though the desired solution may be within. Pseudo-degeneracies occur in this situation both because there is an abnormal condition imposed and because the input prior for the model is poor.

While technical issues with code operation add a layer of complexity to image analysis, in practice the majority of situations one encounters are straightforward to recognize by observing when the parameters take on extremely large or small values. However, clearly recognizing the problem as being pseudo-degeneracies is the key. Once diagnosed, these situations are easy to guard against, by holding those parameters fixed when they go below certain values. In practice, technical issues are not problematic even when Galfit is used for automated analysis (Häussler et al. 2007).⁷

One of the most common causes of degeneracy problems in galaxy fitting analysis comes from using input priors that are not well suited to the data. The most common "input priors" involve the choice of the type or the number of components in a model.⁸ Input priors are ideal when the number of components of a model used in a fit matches the number of luminous components in a galaxy. However, often times one may choose to use either fewer or more components than needed by the data.

A common example where the input prior is bad is when one uses fewer model components than called for by the data. Two of the main reasons for doing so are to reduce the number of components/free parameters, or to allow automated analysis, where it is not yet possible to tailor fits to individual galaxies. This approach is often an intentional course of action taken by many studies, especially when it comes to automating two-component analysis; the goal, ostensibly, is to decompose a galaxy into bulge and disk components. Seemingly reasonable and justifiable on the notion of reducing the potential for degeneracy, the approach is generally regarded by most people to be a positive attribute, rather than a source of problem itself. Yet, that intuitive notion conflicts with the basic principle of how least-squares algorithms work, and leads to perhaps the most common causes of (pseudo-)degeneracy problems cautioned by the literature.

To understand why using fewer components than necessary is bad, it is important to appreciate that galaxy fitting analysis is fundamentally flux weighted. Thus, when a luminous structure is not accounted for, other subcomponents try to compensate, however imperfectly, for its presence. For instance, one may use a two-component model fit to a galaxy that has a bulge, disk, and bar. Doing so may have several different outcomes. One solution is where one component is a sum of (disk+bar) while the other is the bulge. Another can be (bulge+bar) and disk, or perhaps a compromise (e.g., bulge + 0.7 bar; disk + 0.3 bar). Which scenario occurs depends on the relative contrast (i.e., flux weighting) of the bar to the bulge and disk, and potentially on the initial parameters of the three components; small perturbations may "bump" the solution out from one minimum into another. It is quite possible for there to be a "global minimum" solution to this problem. However, when the most meaningful solution, physically, is simply disallowed by the input prior, a globally minimum χ² cannot lend much credence to the reality of the model components.

An input model prior might also be bad if the model involves using more subcomponents than inherently present in a galaxy. In this situation, the results depend strongly on the degree of profile mismatch between the model function and the data. If there is a significant mismatch, all the components cooperate to reduce the residuals. For instance, it is always possible to fit multiple exponential models to a single-component de Vaucouleurs profile. If the goal is to obtain the total flux, the sum would do a better job than using a single exponential. However, individually, the structural parameters may not have much physical meaning.

Another example involving model prior is in the area of high-contrast imaging, where the goal is to deblend a central, unresolved, point source from a diffuse underlying object (e.g., quasar and host galaxy). To do so reliably requires an accurate PSF model for the unresolved source, or else the residuals may overwhelm the extended object, causing unreliable fits. Here, the prior is the PSF model. Quantifying how the prior affects the fitting results involves trying out different PSFs, or to include extra components to account for the PSF residuals, depending on the science goal.

These examples illustrate some of the most common situations where the reliability of a fit depends less on the number of free parameters and more on having a proper model to describe the data. Beyond a single-component analysis, the need to make such a decision means that it will be difficult to automate highly detailed decompositions of galaxies. However, while multi-subcomponent fitting is difficult to automate, it is reassuring that making a wise decision, interactively, is often not particularly difficult when a science goal is clearly defined. Moreover, for galaxy surveys, where the goal is to fit single-component profiles to galaxies, multi-object decomposition is quite feasible to automate (e.g., M. Barden et al. 2010, in preparation; Häussler et al. 2007).

In summary, pseudo-degeneracy conditions exist because least-squares fitting fundamentally involves flux weighting: when luminous flux distributions are present in an image, the models are attracted toward them to reduce the residuals. Therefore, when all components are not properly modeled, the result may be tricky to interpret not because of potential for model degeneracies, but that the solution may have no physical meaning even if there is a global minimum. The solution is to increase the complexity of the analysis progressively until all luminous components are properly accounted. This process does not imply, however, that it is necessary to account for every component inside a galaxy for the solution to have any meaning, only that components of similar flux ratios ought to be simultaneously accounted in detailed analysis; with a few exceptions (e.g., locally dominant features like nuclear star cluster, nuclear ring), components with low fluxes generally do not significantly disturb the parameters of the much more luminous subcomponents.

One of the most common notions regarding fitting degeneracy is that the more free parameters there are, the greater is the potential for degeneracy problems. However, the sheer number of parameters is often not itself an indication of a potential for parameter cross talk. Consider, for instance, that it is equally robust to fit thousands of well-separated stars as it is to fit an isolated one. The same is true for galaxies, even though they are considerably more extended and may overlap: in large-scale image simulations, Häussler et al. (2007) studied automated batch analysis of galaxies using one Sérsic profile per galaxy. They find that simultaneously fitting overlapping or neighboring objects using multiple components (often 3–10 Sérsic models at a time) recovers the input simulated parameters more accurately than fitting a galaxy singly while masking out the neighbors.

Indeed, spatially well-localized sources, like a bar, ring, or off-nuclear star clusters, are virtually free from degeneracies caused by cross-talk with other components. A galaxy bar is well determined because it is more elongated, has a flatter radial profile, and is more sharply defined than the surrounding bulge and disk components, despite being embedded within. Compact objects that are off-centered may also be well determined if the rest of the galaxy can be modeled accurately. Contrary to notions that more model components lead to greater degeneracy, it is important to consider the qualitative aspects of those components: not accounting for strong features explicitly can yield a less reliable and less physically meaningful fit because the solution is a compromise between the different subcomponents.

The issue of parameter degeneracies closely ties into the topic of measurement uncertainties, especially when the result of the analysis may depend on the input model in fitting galaxies. When the model fits the data perfectly (i.e., the residuals are only due to Poisson noise) it is possible to infer parameter uncertainties from the covariance matrix of free parameters, which is produced during least-squares minimization by the Levenberg–Marquardt algorithm. In galaxy fitting, ideal situations are often not realized because the differences between the data and the model profile involve not only random (e.g., Poisson) sources, but also systematics from non-stochastic (e.g., profile function or shape mismatch, neighboring galaxies, etc.), and stochastic factors (overall smoothness, for instance, due to star clusters). The one exception is under low S/N situations, when Poisson noise exceeds model imperfection. In most other situations, non-random factors dominate the residuals, causing uncertainties inferred from covariance matrices to be underestimated. Therefore, it is frequently not very meaningful in galaxy fitting to cite measurement uncertainties for the fitting parameters in the traditional sense.

One way to quantify uncertainties, possible in large galaxy surveys, is to allow the scatter of the data points in physical relations (e.g., radius versus luminosity, luminosity versus metallicity, etc.) to articulate the overall uncertainty of the measurements, even if individual errors could not be easily obtained. Such a scatter inherently involves a convolution of several error sources: the intrinsic scatter present in a physical relation, Poisson measurement error, stochastic and non-stochastic systematic errors due to model imperfections. Intrinsic scatter, being a fact of nature, remains present in physical relationships even should the data have infinite S/N, and even if the models are perfect fits to the data. Intrinsic scatter is often a scientifically interesting quantity, but it is difficult to differentiate from scatter caused by systematic and stochastic errors, which do not vanish given infinite S/N.

In the absence of large galaxy surveys, it is then important to quantify stochastic and non-stochastic systematic errors for individual objects. Some example situations include the black hole mass versus galaxy relation studies (Kormendy & Richstone 1995; Gebhardt et al. 2000; Ferrarese & Merritt 2000) and the fundamental plane (Djorgovski & Davis 1987).

In general, it is very difficult to pinpoint all the causes of non-stochastic systematic errors in an analysis, and to quantify their magnitude. However, one common cause is profile model mismatch: to the extent that one does not know the intrinsic model of a galaxy a priori, the uncertainty in measuring the parameters is wedded to one's assumptions about the model. Therefore, the process of quantifying systematic, model-dependent errors involves exploring the degree of parameter coupling, by trying out different models and seeing how the parameters of key scientific interest change. Another source of systematic error is due to comparing results from different algorithms. In this scenario, the most sensible practice is therefore to only compare parameters that are derived using the same fitting technique (rather than 1D versus 2D), and using the same pixel and flux weighting scheme (instead of Poisson versus non-Poisson) during analysis.

In contrast, stochastic errors arising from general non-smoothness of a galaxy profile are caused by, for example, star-forming patches, dust lanes, etc. Existing on small scales and widely dispersed, non-smoothness cannot be easily identified and modeled in a practical manner using multiple components. Even if it is possible to do so, whether they ought to be fitted explicitly, masked, or not treated at all, falls under the purview of the science goal. Stochastic fluctuations often influence the analysis in a manner analogous to having large correlated noise in the data. If the fluctuations can be quantified, one possible solution is to include them in the fit as a variance term of χ² (Equation (1)). To estimate the fluctuations requires first obtaining a smooth underlying model, which is not always easy to do if galaxies have steep and/or irregular profiles.

While it is generally difficult to disentangle stochastic from non-stochastic sources of systematic errors, there also do not seem to be obvious benefits for doing so from a scientific standpoint. For most applications, one should only be interested in the overall magnitude of the systematic errors in a collective sense. One way forward is therefore to understand which parameters are most strongly coupled, then compare the results of different solutions judging by which ones are physically plausible. For instance, one common interest in bulge-to-disk (B/D) decomposition is to quantify the uncertainty of the Sérsic index n. We know that the Sérsic index n takes on a large value when a profile has both a steep core and an extended wing (see Figure 3). Therefore, quantifying systematic errors in measuring the Sérsic n might involve masking or fitting nuclear sources/neighboring contamination, trying out different PSFs, or fitting the disk by allowing for different disk Sérsic index values. Properly judging the causes of systematic errors and accounting for them often would lead to more natural fits and more sensible parameter values, without the need to hold certain parameters fixed to preconceived answers.

In exploring the parameter space as described, there is often a concern that parameter degeneracies are too numerous or problematic to understand, which brings the discussion back full circle. As discussed in previous sections, when the cause of parameter degeneracy is properly identified, and when the model priors are well conceived, our experience has been that spatial information in 2D can often effectively break many potential degeneracies between the components. Even when the size, luminosity, and central concentration of the different components correlate, they often interact in fairly superficial ways, and do not dramatically change what the model components represent physically. However, in situations where cross talk is significant and there is no reason to prefer one solution over another (when the input prior is befitting), then differences in the answer speak to the degree of the parameter uncertainty that is of key interest to quantify, rather than to avoid, because ultimately the models are empirically motivated.

In summary, to the extent that the results may depend on model assumptions, parameter exploration is the only viable way to quantify true measurement errors in the fit parameters. Thus, when used properly, parameter coupling/degeneracy, rather than complicating the interpretation, offers a deeper insight into the reliability of the overall analysis. We illustrate the above ideas more explicitly in the following examples.

IC 4710 is an SB(s)m galaxy, which has a bar-like feature in the middle of a roundish outer structure, as shown in the R-band image of Figure 18, which comes from the Carnegie–Irvine Nearby Galaxy Survey (CINGS) project.⁹ Prior to the analysis, we masked out the stars using the SExtractor software (Bertin & Arnouts 1996). We analyze this galaxy using, for comparison, both one- and two-component regular and higher order models with Fourier modes, shown in Figures 18(b)–(i). The best-fit parameters are given in Table 1, which illustrates three different sets of analysis parameters: best fit using two components (Figure 18(c)), a model using just the traditional ellipsoid component (Figure 18(g)), and the same single-component model with Fourier modes added (Figure 18(b)). Figure 18(i) shows the radial surface brightness profile of the data and the individual subcomponents of the best model.

Zoom In Zoom Out Reset image size

There are several points to understand from comparing detailed and simple analyses. The best-fitting ellipsoid model (Figure 18(g)) is oriented more parallel to the bar-like, higher surface brightness component than the lower surface brightness body; however, the ellipsoid model is much broader than the bar (Figure 18(e)). This happens because a single-component fit is a compromise between the various subcomponents of a galaxy, and, as such, it reflects neither one perfectly. Allowing the azimuthal shape to change by adding nine Fourier modes results in a shape shown in Figure 18(b). Note that because the profile is restricted to having a Sérsic functional form in every direction radially from the peak, the shape does not have complete freedom to take on any shape, as opposed to a shapelet or wavelet-type Fourier inversion: it is merely a higher order perturbation of the best-fitting ellipse. Indeed, in comparing single-component fit parameters in Table 1 for the two models, the main Sérsic structural parameters hardly budged, despite the Fourier model having 18 more free parameters. Therefore, the marginal returns in using more free parameters are negligible in this situation when it comes to the main Sérsic structural parameters. However, if the scientific interest is to quantify the global symmetry, then the higher order modes are of interest.

Another point of interest is how higher order models affect the accuracy of the global photometry. It is natural to expect when a model is unrealistic for a galaxy that the photometry is also unreliable. In Figure 18(c), it is evident that a two-component model is more appropriate than the single-component fits of Figures 18(b) and (g). However, when the flux of the two-component model is summed, one finds that the difference with the single-component fits is only 0.03 mag. This and subsequent examples illustrate empirically that the process of least-squares minimization using even naïve ellipsoids is often capable of providing accurate photometry to within 0.1–0.2 mag, even if the galaxy shape departs from ellipsoid models quite drastically.

Lastly, Figures 18(e) and (f) demonstrate that it is quite feasible to unambiguously disentangle embedded components that have different shapes, using higher order Fourier modes. Despite there being a large number of parameters, it is visually clear based on Figures 18(e) and (f) that parameter degeneracy is not an issue, because the shapes of the components are quite different. In part, this is possible because of the way Fourier modes are implemented in Galfit: the profile function is preserved in every direction radially from the peak, even in situations where the shape is irregular, as in Figure 18(e).

This edge-on galaxy (Figure 19, Table 2) comes from the Galaxy Evolution from Morphology and SED (GEMS; Rix et al. 2004) project, which is an HST imaging survey of the Chandra Deep Field-South. Belying a benign morphological appearance is a dust lane (Figure 19(e)) that courses through the center, complicating the traditional ellipsoid fitting technique.

Zoom In Zoom Out Reset image size

The analysis of even this simple object can be quite involved. The best-fitting model involves three components: a fairly compact bulge, an edge-on disk component, and a puffy stellar halo enveloping both. Since the halo component is more luminous than the bulge component, a two-component model fit would naturally ascribe the halo component to the bulge, despite there being a distinctly rounder component at the center. Like the previous example, each of the three components (Figures 19(f)–(h)) is modified by Fourier modes. Furthermore, the best fit includes an actual model for the dust lane (component 4, Table 2). The dust lane is modeled by an inner truncation function as discussed in Section 5.

A truncation model is shown as a model "component" in the fit; it is unique because it is not a light profile model, and one cannot generate an image to see what it looks like. Instead, its influence is to be seen on all the light profile models (i.e., components 1–3; Figures 19(f)–(h)), where it reduces the light by the same fraction for all components. In every other way, the truncation function behaves exactly like a light profile model: it can have its own centroid (or not), and it can be modified by Fourier modes, as shown in Table 2. The benefit of using a single truncation model for all three light profile models is not only to reduce the degrees of freedom, but it is also physically motivated because foreground dust attenuates all background light sources by an equal fractional amount. Nevertheless, if desired, it is also possible to allow each component to be attenuated differently.

This example also demonstrates how the result of the analysis depends on the input prior of the model. In the fit using traditional ellipsoid parameters, a mask is used to minimize the effect of the dust on the analysis. Yet, the effects cannot be completely removed. As shown in Table 2, the inclusion of the truncation model can significantly affect the structural parameters: the surface brightnesses can differ by 0.8 mag arcsec⁻², and the sizes by 10%–20%, even in this seemingly uncomplicated situation. Moreover, the differences far surpass the statistical uncertainties shown in Table 2. To the extent that it is not possible to judge which model is more physically correct, both measurements ought to be treated as equally valid. In that situation, the uncertainties, due entirely to model assumptions, are roughly ∼ 0.4 mag in surface brightness and ∼ 10% in size.

The HST/F814W image of the field Arp 147 contains two ring galaxies (Figure 20, Table 3), one of which has a bulge-like component with a tidally disturbed outer region (Galaxy 1), and the other is a pure ring (Galaxy 2). The best-fitting model for Galaxy 2 is a single-component ring, modified by Fourier modes, as seen in Figure 20(b), whereas Galaxy 1 requires two ring components, a bulge, and an inner fine-structure component (Figures 20(e)–(h)). The fine-structure component of Galaxy 1 can easily be seen in the surface brightness profile as an upturn within r = 02 of Figure 20(i, left). In addition, the tidal component is slightly bent, which is modeled elegantly using the bending modes of Equation (27). As in the case of a dust lane, the ring model comes about by truncating the inner region of a pure Sérsic profile (see Section 5. The only difference here is that the truncation radii are a larger fraction of the galaxy size. Whereas for the edge-on galaxy, it makes more sense to normalize the flux at the effective radius (Equation (32)), for ring galaxies, normalizing the flux at the break radius (Equation (33)) is more intuitive, because it is closer to the peak of the profile model. In fact, the peak of the ring is about half-way between r_break and r_break + Δr_soft, but the exact location depends on the profile type.

Zoom In Zoom Out Reset image size

It is again instructive to compare a traditional fit using simple Sérsic ellipsoid models (Table 3, bottom) with more sophisticated analysis (Table 3, top). In terms of the total flux for Galaxy 1, the magnitude of the most sophisticated model is m = 14.18, compared to m = 14.32 for a model based on classical ellipsoids. Interestingly, a single-component fit (not shown) to Galaxy 1 yields a magnitude of m = 14.15. For Galaxy 2, we know from the outset that classical ellipsoid models are entirely inappropriate to use. Yet, despite every reason to believe that the photometry would be inaccurate, we find that the total flux of the traditional ellipsoid fit is only 0.2 mag different from the most realistic ring model. These two examples show once again that a single-component Sérsic ellipsoid fit to complicated galaxies can produce quite accurate measurement of the total flux.

It is sometimes desirable to conduct B/D decompositions, and Galaxy 1 is an ideal candidate to conduct a comparison. In the traditional ellipsoid model (Table 3, bottom), the B/D ratio is 0.65. The more sophisticated model (Table 3, top) requires summing the ring+tidal feature components to obtain the disk component, which yields 14.65 mag, thus a B/D ratio of 0.54. In this situation, most of the differences arise from measuring the disk component, which differs by 0.2 mag, whereas the bulge component is quite robust, with a difference of only 0.01 mag.

It is also of interest to understand how the structural parameters are affected by different model choices, in particular for the ring Galaxy 2. Whereas the effective radius for the ring model is only 078, for the ellipsoid model it is 828. This is understandable, bearing in mind that the ring has a radius of nearly 8''. To a classical Sérsic profile, the galaxy appears to have a very flat (in fact, a deficit) core, which leads to a low Sérsic index of n = 0.12. As most of the flux is at 8'', beyond which the ring flux quickly fades, the ring radius is closely related to the effective radius for a classical Sérsic model. For the inner-truncated ring model (component 7), however, the physical size of the ring is captured by the break radius r_break parameter, whereas the r_e term no longer has the classical meaning of the effective radius (i.e., half the light is within r_e). Instead, r_e for component 7 is essentially an exponentially declining scale-length parameter, given by Equation (33). As the flux dies away quickly beyond the peak, as shown in Figure 20(i), the scale length r_e for the ring model must therefore be quite small. The differences in the r_e parameter between the traditional model and the truncated model are therefore only due to definitions, and not due to systematic or random measurement uncertainties.

The classical Whirlpool galaxy is a beautiful system where a grand-design spiral, M51A, is interacting with another spiral, M51B (Figure 21, Table 4). In addition to there being obvious spiral structures for both galaxies, there are large tidal disturbances that emerge from M51B, as seen in the Sloan Digital Sky Survey (SDSS) r-band image provided by D. Finkbeiner. Because they are closely overlapping, a desirable goal is to deblend M51A and B, as well as to model the spiral and tidal structures, simultaneously.

Zoom In Zoom Out Reset image size

As with previous examples, we fit this galaxy using both the most sophisticated analysis (Table 4, top) in our toolbox, and comparing the results to the traditional axisymmetric ellipsoids (Table 4, bottom) analysis. The traditional analysis requires two components each, in order to decompose a galaxy ostensibly into a bulge and a disk. The reduction in χ²_ν between the two methods is modest, because most of the residuals come from high-frequency star-forming regions that are not removed by models which are fundamentally smooth, despite being modified by radial Fourier modes and spiral rotations.

In the most detailed analysis of M51A, we use two spiral arm components and two components for the bulge. There is actually not a strong need to use two components for the bulge except to better capture the detailed profile shape, which has an inflection at r ≈ 04, as seen in Figure 21(i). On the other hand, the use of two spiral components is necessary because the spiral arm has a "kink" in the rotation that cannot be created by using a single smooth rotation function. The spiral structures are modified by Fourier modes to create both a slight lopsidedness and other subtle features. Because there are more degrees of freedom in a two-arm spiral, the higher order Fourier modes also can "see" detailed structures, like the reverse flaring of the spiral structure in Figure 21(f).

For M51B, we employ three components in the fit, a bulge (component 5 in Table 4, top), a tidal feature component (component 6), and a spiral function (component 7), which model the three most visually striking components. The tidal feature is mostly obtained by using the second and third bending modes of Equation (27), as illustrated in Figure 10. However, bending modes 2 and 3 are symmetric functions, so the high degree of asymmetry comes about because of combined action with the Fourier modes, which is shown to have a high amplitude of 0.23 for the m = 1 mode, as well as moderate values for other modes. Incidentally, despite the complexity of the higher order structures, all the parameter values are determined automatically by Galfit without the need for an user to provide initial guesses (i.e., initially all 0 values) and without tweaking at any point in the analysis (which is hardly feasible anyhow).

For even those who are experienced with detailed parametric fitting, one of the alarming facts about this analysis is that it employs 103 free parameters in the best-fit model. So there are natural concerns about parameter degeneracies. However, as we have discussed in Section 6.5, parameter degeneracies do not arise purely based on the number of free parameters, but rather on the types of parameters involved. The availability of spatial information in 2D provides one of the most important ways to break parameter degeneracies. We see this explicitly in Figures 21(e)–(h), where there is little evidence that the subcomponents for M51A are strongly influenced by M51B, and vice versa. Furthermore, within each galaxy, the subcomponents are so different in shape, both qualitatively and quantitatively, that the amount of cross talk between them is also not significant. Therefore, despite the extreme complexity of this system, and the use of 103 free parameters, we find that degeneracies between the parameters are not an issue. Or, if they exist, they do so at a low enough level that they do not significantly affect the main parameters of interest, like the luminosity of the subcomponents, or the profile shapes and sizes.

There are, however, seemingly degenerate conditions that have little to do with parameter coupling. Instead, these are attributed to the fact that M51 has many non-smooth structural features, caused by dust lanes, star-forming regions, tidal disturbances, and so forth. Such spatially localized features, if strong enough, can influence Galfit to "lock" on to them if the initial conditions happened to be sufficiently close. The consequences appear as degeneracies when, in fact, there are many small local minima solutions. This graininess in the χ² terrain introduces slight perturbations to the models, and may even cause fairly large shape differences in the final solutions. However, to a large extent, it rarely affects the main parameters of interest, such as the luminosity of a particular component or its size, which are determined by much more global features than the nuisances of local fluctuations to which higher order parameters are more sensitive.

To gain some intuitive insight into the effects of complex analysis, it is instructive to compare simple and complex methods with regard to global and subcomponent properties. In terms of the total luminosity, here we find excellent agreement between sophisticated and traditional analysis, respectively, of m_r = 8.24 versus m_r = 8.25 for M51A, and m_r = 8.80 versus m_r = 8.73 for M51B. While this level of agreement may at first seem surprising, it is expected given the basic premise of least-squares minimization. In fact, even a single-component fit to M51A yields m_r = 8.0, and for M51B m_r = 9.0, which are both quite close to the overall best-fit models, despite the complications in the image. The main reason for the discrepancy here is the uncertainty in the sky, due to there being a large gradient. This fundamentally sets the limit on the accuracy of the photometry to perhaps no better than 0.1 to 0.2 mag, independent of the analysis method.

The most sensitive benchmark for understanding differences in the analysis is in detailed decompositions. Here we compare the B/D decomposition results. In the traditional ellipsoid analysis, we find a B/D ratio of 0.23 for M51A and 4.9 for M51B. The large B/D ratio for M51B is clearly unphysical, and is driven by the large Sérsic index (n = 8.0) of the bulge component, which is increased to accommodate the flux in the outskirts due to tidal features. In the most detailed analysis, the B/D ratio for M51A is 0.16, whereas for M51B it is merely 0.17. Examining the bulge of M51A more closely, we find that the detailed analysis yields a total flux of 10.38 mag, whereas the traditional analysis extracts a brighter bulge of 10.05 mag. The differences come from the fact that the light of the inner spiral is in part driving up the Sérsic index of the bulge when it is not properly accounted. It is probably safe to conclude that a magnitude of 10.05 is a firm upper limit to the bulge luminosity.

Finally, it is worthwhile to compare how the disk parameters differ between the analyses to gain an understanding for how coordinate rotation affects the interpretation of the parameters for the spiral models. From Table 4, we find that the Sérsic indices of the simple and complex models are essentially identical for M51A, at n ≈ 0.33. The interpretation for M51B is more complicated, because the "disk" in an ellipsoidal model is not qualitatively the same structure as the spiral analysis. In fact, it is necessary to hold the Sérsic index of the tidal component 6 fixed in the analysis. Nevertheless, there are clearly quantitative differences in that the simple analysis is larger by 55% in n. With regard to the effective radius, the traditional analysis of M51A finds the disk size to be about 22, which compares favorably with the spiral model size of 28, or a 25% difference. Furthermore, the disk magnitudes for M51A differ only by 0.03 mag between simple and complex.

These comparisons therefore demonstrate that despite the complex analysis being much more realistic looking, fundamentally the meaning of the structural parameters (size, luminosity, concentration index) is unchanged from the original definition, even in the situation of spiral components. This is a useful fact because our prior intuitions, honed on fitting ellipsoidal models, continue to be applicable. We note that the generally good agreement between detailed and simplistic analysis witnessed here and in previous examples is not entirely coincidental. It so happens because all shapes are fundamentally perturbations of the best-fitting ellipsoidal model, even if the result bears no resemblance to the original ellipse.

NGC 289 is an SAB(rs)bc galaxy, with a weak bar and a complex inner spiral system (Figure 22, Table 5) that resembles a ring. Upon closer examination, the ring appearance comes about because there exists a bifurcation in the spiral structure that connects up with the opposing spiral arm. Furthermore, the bar is also multi-component, with the inner component oriented at an angle nearly 45° from the strong outer bar.

Zoom In Zoom Out Reset image size

The best-fit analysis involves three spiral components, an inner and an outer bar component, and a bulge (Table 5, top). All except for the bulge component are modified by five Fourier modes, and are shown in Figures 22(e)–(h). The requirement of components 3 and 4 (Figures 22(f) and (g)) is clear, because they are what form the most striking and intricate patterns in the center, while the requirement of component 5 (Figure 22(h)) is only evident in the residuals, and makes up some of the diffuse light within the inner 60'' region. Although it does not seem like an essential component, the inner bar structure (Figure 22(e)) qualitatively affects the detailed residual pattern at the center, and is therefore included. When all the detailed inner structures are properly accounted for, it is straightforward to infer the bulge component, and assess the uncertainties by varying different parameters of the bulge. Doing so does not affect the inner fine structures because they are sharp and well localized.

Conducting the same decomposition using traditional ellipsoid models (Table 5, bottom), we opted to fit three components, ostensibly to model a bulge, disk, and a bar. The result produces residuals seen in Figure 22(d), revealing the intricate details of the inner spiral system. From the fit, even though the disk and the bar component are sensible, the bulge component is actually fitting a diffuse disk component, which, in retrospect, is that shown in Figure 22(h). Because that inner spiral component is quite luminous, and because there exists a bulge component superposed on top of it, this quasi-bulge model is almost 0.7 mag brighter than that inferred through the detailed modeling above. Adding a fourth ellipsoid model is not possible, because the central spiral residuals are so great that they completely suppress the addition of another component, causing the flux to go to zero.

Once again, comparing the total luminosity between the best-fit model with the ellipsoid fit, we find an excellent agreement of m = 10.37 mag versus m = 10.42, respectively. For a single-component ellipsoid fit, there is also an excellent agreement of m = 10.46, despite the main structural details not being unaccounted.