© ESO 2022

1 Introduction

All information that we can gather from galaxies to constrain models of galaxy formation and evolution is encoded in the light that they emit. From this single source, we need to infer the physical processes that form these objects and distinguish different evolutionary mechanisms (e.g. different quenching mechanisms, continuing star formation, merger history, etc.). Squeezing as much information as possible from this limited resource is hence vital for understanding how galaxies form and evolve.

Modern surveys and analysis have come a long way, partly through the technical and hardware advances of instruments and telescopes and partly through advances in the analysis and software packages used. One such advance in instrumentation is the gathering of integral field unit (IFU) data, which provide spectral information for each part of the sky or galaxy that can be explored in detail in a sophisticated manner. However, gathering these iFU data is very expensive as it requires vast amounts of telescope time and delivers data only for individual or small samples of targeted galaxies or – worse – only the centres of those galaxies due to the limited field of view (FoV) of these instruments at the present time. Only recently have IFUs with a larger FoV become available (e.g. MUSE, Bacon et al. 2010, 1 arcmin²), but even using those, it is – for now – impossible to gather statistically useful samples of galaxies in a cosmological volume. The largest current IFU surveys have observed ~600, ~3000, and ~210 000 galaxies (Califa, Sánchez et al. 2012; SAMI, Croom et al. 2021; and MANGA, Bundy et al. 2015, respectively). Additionally, these surveys are limited to low red-shifts of z ~ 0.03 (Califa) and z ~ 0.1–0.15 (MANGA and SAMI) and still lack a comparison sample at higher redshifts, vital for formation or evolution studies that span a significant fraction of the Hubble time. This redshift issue might be solved somewhat using different wavelengths in different instruments, for example KMOS (Sharples et al. 2013), but this does not solve the issues of sample size and the observing time required for statistically meaningful studies.

On the other hand, imaging of large survey areas can be – and routinely is – done relatively cheaply, for example by SDSS (York et al. 2000), 2MASS (Skrutskie et al. 2006), Cosmos (Scoville et al. 2007), Candels (Grogin et al. 2011; Koekemoer et al. 2011), GOODS (Dickinson et al. 2003; Giavalisco et al. 2004), GEMS (Rix et al. 2004), VIDEO (Jarvis et al. 2013), Combo-17 (Wolf et al. 2003), and countless others. This cheap access to vast amounts of data and large galaxy samples, however, comes at the price of significantly poorer spectral resolution. For example, the Cosmos field is now covered in ~30 different filters with a sky coverage of ~2 deg², Alhambra (Moles et al. 2008) observed ~4 deg² in 20 filters, and JPAS (Benitez et al. 2014) is planning to observe ~8000 deg² in 56 medium-band (and non-overlapping) filters, effectively creating high-quality spectral energy distributions (SEDs) at each image pixel. However, while these datasets do deliver a huge amount of information, only part of it is generally exploited in present-day analyses, for example via aperture photometry. However, even these simple measurements contain valuable information (e.g. one can turn the measured colours of a galaxy into estimates of global stellar masses, star-formation rates, etc.), which allows the development of a general picture of the merger and star-formation history of a given galaxy. Typically, however, a variety of star-formation histories can produce similar integrated properties, making this simple approach sensitive to measurement uncertainties, which are increased when measurements at different wavelengths are carried out independently. This is further complicated by the fact that many of the tools used rely on aperture photometry, which – while summing up all the light from an object – entirely ignore the light distribution within the aperture radius. Different galaxy shapes or the internal distribution of light (e.g. in different components) are lost in such a simple measurement.

A different way to measure galaxy parameters is to measure the distribution of light using non-parametric codes, for example via Gini-M20 (Lotz et al. 2004) or CAS (Conselice 2003). To date, these techniques – again – only work independently in each band, which makes it difficult or impossible to make consistent measurements on images at different wavelengths. Additionally, these and other non-parametric techniques generally do not take point spread function (PSF) effects into account and can hence bias the results since different fractions of the true light of a galaxy are measured at each wavelength due to the different sizes and shapes seen in the PSFS at different wavelengths.

One widely used technique that makes it possible to take the image PSF into account is parametric light profile fitting. Several codes for this work exist – for example Budda (de Souza et al. 2004), GIM2D (Simard 1998), GALFIT (Peng et al. 2002, 2010), IMFIT (Erwin 2015), PROFIT (Robotham et al. 2017), 2DPHOT (La Barbera et al. 2008), and Profiler (Ciambur 2016) – and are more or less routinely used to measure galaxy parameters both on nearby galaxies and in an automated fashion on large samples of galaxies in large-scale galaxy surveys. All these codes, however, suffer from the same problem, in that they – once more – only allow profile fitting in a single observed band¹.

To improve on this situation, several authors have recently developed techniques in which a light profile is fit in one band and then applied to the other bands to derive magnitudes (Bruce et al. 2014a,b; Head et al. 2014; Simard et al. 2011; Mendel et al. 2014). While this does paint a somewhat more consistent picture, it still loses information present in the imaging data themselves since it is not clear whether the profiles are actually the same at all wavelengths, as has to be assumed for these analyses. For example, it is already known that a typical disk galaxy today generally contains an extended blue disk and a compact red bulge. The same is also often found in S0 galaxies (e.g. Bothun & Gregg 1990; Peletier & Balcells 1996; Head et al. 2014; Hudson et al. 2010), but exceptions to this general rule are known for individual S0 galaxies (Johnston et al. 2021) and dwarf galaxies. A single profile fit to such a galaxy at different wavelengths would hence tend to fit something larger with a low Sérsic index (Sérsic 1968) if one of the bluer bands is used in the fit, and something smaller and with a higher Sérsic index if one of the redder bands is used, making it hard to come to a consensus regarding a consistent measurement for different science cases (e.g. at low and high redshift). In fact, several authors (e.g. Kelvin et al. 2012; Vulcani et al. 2014) have reported such variations in galaxy parameters with wavelength, which, using this approach, would stay undetected by design. At least some of these studies have indeed attributed much of the observed change to the different mixing of bulges and disks at different wavelengths. There is, however, a niche market for these kinds of analyses, where the assumption of ‘no wavelength dependence’ is valid, for example, in first order for bulge–disk (hereafter: B/D) decomposition of high signal-to-noise (S/N) galaxy images, where stellar population gradients within the individual galaxy components can be ignored.

For a full multi-wavelength analysis (or what was called ‘oligochromatic’ by De Geyter et al. 2013), it would however be preferred if all images were taken into account simultaneously in an equal manner (i.e. giving the same weight to each image while leaving any possible wavelength dependence to be accommodated for). As part of the MEGAMORPH project (Measuring Galaxy Morphology), we developed and tested such a technique (presented in Häußler et al. 2013; Vika et al. 2013, 2014), which allows the combination of some of the above techniques while avoiding some of the obvious shortcomings. The software developed by the MEGAMORPH team – Galapagos-2 and Galfitm – are based on Galapagos and GALFIT, two publicly available and well-tested single-band codes². Galfitm, however, allows the simultaneous fitting of multiple images at different wavelengths, and Galapagos-2 allows the automatic exploitation of this code on large-scale surveys. Not only does this enable more accurate measurements on simple, single-component profiles (e.g. when fitted values do not actually change with wavelength), as it uses more data and effectively increases the S/N, it also allows structural parameters of galaxies to vary smoothly with wavelength, making additional measurements possible in a physically meaningful way. This has been shown and exploited in Häußler et al. (2013), Vulcani et al. (2014), Kennedy et al. (2015, 2016a,b), Huertas-Company et al. (2016), Dimauro et al. (2018), Mosenkov et al. (2019), Psychogyios et al. (2020), and Nedkova et al. (2021). While these single-component fits provide a great deal of useful information – and are less challenging to perform – there is a further advantage of this multi-band approach when fitting multiple profiles to a galaxy in order to decompose the galaxy into its constituents. The main (first-order) reason why galaxies look different at different wavelengths is that different galaxy components are mixed at different strengths at different wavelengths. In blue bands, much more of the disk light is generally observed, while red bands generally observe more light from the bulge component, explaining at least some of the trends found by Vulcani et al. (2014), Kennedy et al. (2015, 2016a,b), and others. Valuable information needed for a successful separation of these components is hence stored in the colour information of each observed pixel, making the colour information vital in the decomposition process itself. Although conceptually simple, B/D decomposition remains a challenging task due to the variety of structures that galaxies display, not to mention the usual observational limitations of resolution and S/N. Any improvement and increase in parameter reliability is a big step forward for successfully separating bulges from disks so that the individual galaxy components can be studied independently.

This successful decomposition of a galaxy light profile – or better: mass profile – into its individual constituents is the holy grail of galaxy profile fitting, at least at redshifts outside the local universe, and is a vital step towards understanding galaxy formation and evolution. However, this task is incredibly tricky to achieve, especially in an automated manner on a large sample of not-very-well-resolved galaxies spanning a wide range of redshifts, as are currently routinely imaged by large-scale surveys. Previous approaches to this problem mostly worked on bright, well-resolved galaxies, where the profile shapes themselves contain enough information for a successful separation of the individual components, even on single-band data. Some papers have presented large catalogues of galaxies that were successfully separated into bulge and disk components (e.g. Simard et al. 2011; Huertas-Company et al. 2016; Meert et al. 2013), allowing for further analysis, for example of their masses and stellar populations (Mendel et al. 2014; Dimauro et al. 2018). However, the sample used by Simard et al. (2011, hereafter S11) and Mendel et al. (2014) is limited to m_petro,r,corr ≤ 17.77 in SDSS data, Meert et al. (2013, 2015) presented a similar analysis on a sample using the same magnitude limit.

In this work we present B/D decompositions of galaxies that are nearly 2 magnitudes fainter than the ones discussed in S11 and Meert et al. (2013) in order to show how well a multi-wavelength approach can separate different galaxy components and to quantify the advantages provided by multi-band fitting in the decomposition of distant galaxies in typical present-day surveys. We follow the approach taken by Häußler et al. (2013), using both simulated data to measure deviations from ‘true’ values directly and real data to carry out several consistency checks. We complement and expand on the work of Vika et al. (2014) by demonstrating the application of our technique to large surveys in an automated fashion and with greater statistical power. We target a sample of real galaxies with spectroscopic information in GAMA (Driver et al. 2011; Liske et al. 2015; Baldry et al. 2018) with a limiting magnitude of m_r ≤ 19.8 mag, using -ed (Bertin et al. 2002) data from SDSS and UKIDSS LAS archival data (Lawrence et al. 2007). We thus demonstrate how (and why) using multi-band fitting has advantages in terms of stability, improved accuracy, and increased sample sizes, especially for the low S/N bands of a survey.

This paper is part of a series that investigates the benefits of this multi-wavelength approach to measuring structural galaxy properties. In Häußler et al. (2013, hereafter H13) we presented an overview of Galfitm, with further details available on the project website³. We demonstrated our approach by performing single-component fits on a large dataset from the GAMA survey, automating both the preparation of the data and the fitting process itself. The resulting measurements – in particular, the variation in structural parameters with wavelength – are studied further in Vulcani et al. (2014, hereafter V14a) and Kennedy et al. (2016b). In Vika et al. (2013, hereafter V13) we tested this new method by fitting single-Sérsic models (e.g. Sérsic 1968; Graham & Driver 2005) to original and artificially redshifted images of 163 nearby galaxies. In Kennedy et al. (2015), we used the results discussed in V14a to analyse the colour gradients seen in galaxies versus their colour, structure, and luminosity.

Vika et al. (2014, hereafter V14b), using the same galaxy sample as used in V13, showed that multi-band analysis on nearby galaxies allows a more accurate and more stable separation of two components. Johnston et al. (2017) expanded the use of Galfitm to IFU data to cleanly separate spectra of different components in individual galaxies via the Buddi software. A more complete approach to fitting ~1800 objects of the MANGA survey is currently underway (Johnston et al., in prep.). This work will be based on data from SDSS Data Release 16 (Ahumada et al. 2020), which contains ~5000 galaxies, and will use a sample cut based on both the used fibre-bundle size and the pre-existing fits on SDSS data made by Fischer et al. (2019) using PYMORPH (Vikram et al. 2010).

The objective of the current paper is to present the ability of Galfitm to perform B/D decomposition on galaxy images with a wide range of resolution and S/N, in a fashion similar to what was used by H13 for single-component fits. Kennedy et al. (2016a) already used the technique, code, catalogue, and data sample discussed in this work to explain the wavelength dependences seen by V14a as mostly effects of mixing light from different galaxy components. As such, the tests and results presented here are directly transferable to that paper. Other works have already exploited Galapagos-2 on other datasets, for example to examine Subaru SuprimeCam (Kuchner et al. 2017) or HST Candels data (Dimauro et al. 2018; Nedkova et al. 2021). Their results and code behaviours tie in well with the findings in this paper.

This paper is structured in the following way: Sect. 2 explains necessary changes to Galapagos in comparison to H13, which are essential for the performance in this work. This section also explains the setup of both codes – Galapagos-2 and Galfitm– used throughout this paper (Sect. 2.1), gives a description of the starting values used (Sect. 2.2), and presents the constraints used during the fits (Sect. 2.3). For further details about the technique itself, we refer readers to H13. Section 3 discusses the galaxy sample selection for the remaining part of the paper. Section 4 shows tests when this software was applied to simulated data (i.e. galaxies whose true intrinsic values are known). This comparison, while not containing any physical meaning about galaxy populations, allows us to show the improvement obtained by using multi-band fitting in more detail. In Sect. 4.4 we specifically highlight the dangers of using Sérsic indices as a proxy for the bulge-to-total ratio (hereafter B/T). In Sect. 5 we show similar tests with this software applied to real GAMA data, comparing their values to show how much multi-band fitting improves the fitting results both on individual galaxies and on the galaxy population as a whole. We also show a comparison to other works in Sect. 6, including S11, Robotham et al. (2017), and Casura et al. (2019). After briefly discussing the effect of dust in real galaxies on the fitting parameters in Sect. 7, we discuss our results and conclusions in Sect. 8.

We present the technical details of how the simulated images were created in Appendix C. Appendix D presents a brief analysis of the impact of varying starting parameters in Galapagos-2 and Galfitm. In Appendix E we discuss further improvements of the code in the current Galapagos-2 version, compared to the version used in this work, and discuss possible further improvements. We conclude with a description of the GAMA catalogues released with this paper in Appendix F, based on newer and deeper imaging data, which use re-SWARPed data from the Kilo-Degree Survey (KIDS; de Jong et al. 2013; Kuijken et al. 2019) and the VISTA Kilo-Degree Infrared Galaxy Survey (VIKING; Edge et al. 2013).

We should note at this stage that fitting additional components beyond a simple B/D model might be advisable depending on the science case (see e.g. Kruk et al. 2017, 2018; Davis et al. 2019; Sahu et al. 2019). Specifically, as the spatial resolution of images improves, extra components often need to be added to account for the detailed profile of a galaxy. However, it should also be noted that overfitting galaxy profiles is often possible if the spatial resolution of the data does not allow it, which is as bad as under-fitting the data. Care must hence be taken when selecting the appropriate number of components fit to a dataset. As automating fitting of more and more components and subsequently selecting the ‘best-fit’ model becomes more and more challenging, and in order to do one step at a time, we restrict ourselves to the analysis of B/D models in this paper. In Appendix E we show how this paper and the software can be used as a setup for a more in-depth analysis of additional components in selected objects.

2 Changes to the previous software and setup

Galfitm-v1.4.4⁴ has been used in this work, to incorporate some of the improvements of the code in terms of delivering robust measurements in special cases in comparison to H13. Most improvements concern either Galfitm features not used in this analysis (e.g. the implementation of non-parametric components) or bug fixes that do not affect the results in this analysis due to the way the code is used by Galapagos-2 and/or the way the data are set up for analysis, so improvements can be considered minor⁵. This means, however, that the results in this work are not directly comparable to the work presented in H13, where Galfitm-0.1.2.1 has been used. As that work presented single-Sérsic fits only, no such comparison is carried out here.

Equally, we used a more recent version of Galapagos-v2.2.7 throughout this work. This differs from the version used in Sect. 7 of H13 (v2.0.3) mostly by being able to run B/D fits and a few changes to optimise usages of CPU and disk space. Galapagos-2 has since been developed further, several required and convenient additional features (e.g. as desired when running Galapagos on other datasets) have been introduced in newer versions of the code. However, none of the changes should have any impact on the parameters of the fits themselves. optimise For example, this more recent version of Galapagos-2 can deal with surveys with different footprints in each band. For the data used in this analysis – both simulated and real – all images in all bands show complete sky coverage in the areas used, so such a feature does not have any impact on the fitting results discussed in this work. A more complete summary of these new features is given in Appendix E. The newest version of Galapagos-2 can be derived as both a full GitHub repository or a zipped package⁶, including example setup files and an extensive help file that includes practical help on how to properly set up both code and data. The code is distributed under MIT license and can be used and edited by users. We welcome any changes back into the repository if they provide new features, and if they are well tested.

The main purpose of this section is to present and discuss the setups for Galapagos-2 and Galfitm used in this work. Equivalently to H13, we use PSFS from the GAMA survey as created by SIGMA and used in Kelvin et al. (2012). For the simulated images, a mean of these PSFs is used for both images simulation and Galapagos fitting in each band, effectively removing any PSF uncertainty effects from this part of the analysis. Images, object detection and analysis were identical to those used in H13, with the addition of a B/D decomposition being carried out on all objects possible. In fact, for practical reasons, the B/D fits presented in this work are entirely based on the single-Sérsic fits analysed in H13, including the postage stamps for each galaxy, the decision on masking or de-blending, and so forth. We hence assume the reader to be familiar with this paper.

2.1 General code behaviour

While most of the code setup and behaviour is analogous to the single-Sérsic setup, some additional issues have to be considered when carrying out B/D decompositions. Most importantly, we had to decide on the degrees of freedom (DOF) and hence the order of the polynomial used in the B/D decomposition and the starting values for the fit itself. In this section we discuss the general setup of the code and its logical flow in Sect. 2.1.1, and will discuss possible improvements to the code in Sect. 2.1.2 in this context. In Sect. 2.2, we discuss how the starting values for the B/D fits are derived in the code, while Sect. 2.3 we explain the parameter constraints used. As Galfitm and Galapagos-2 are closely linked in this work, separating them is difficult, and we present the general behaviour of Galfitm and, where relevant, the specific implication for Galapagos-2. In Sect. 2.4, finally, we explain the choices made for the analysis in this work in particular, for example the degrees of the polynomials chosen throughout this paper. The main values from these sections are summarised in Table 1.

2.1.1 Code setup and features – Logical flow

Both codes, Galfitm and Galapagos-2, allow great flexibility with respect to the parameters used during the fit, most importantly, the degree of freedom allowed. A user can trivially switch on wavelength variations in the setup files and define the degree of the polynomial used in these variations. They can equally trivially hold the Sérsic index fixed at n == 1 for galaxy disks and n == 4 for bulges if desired (see Appendix E).

The B/D decomposition of a galaxy in Galapagos-2 is set up by directly using the single-Sérsic fit result. Galapagos-2 creates a new setup file for Galfitm, starting with the information available from the single-Sérsic fit and the parameters specified in the Galapagos-2 setup. In this second fit, most of the settings are kept the same. Specifically, the same postage stamps are used for all galaxies as well as the same masks and PSFs. The wavelengths to be given to Galfitm⁷ are identical. As the sky values are directly set to the same values used of single-Sérsic fits, they are not re-measured, but kept at the same values as the initial single-Sérsic fits. Finally, the same decisions regarding fitting (secondaries), masking (tertiaries) and de-blending neighbouring galaxies are used as decided for the single-Sérsic fits. Single-Sérsic and two-component fits are hence directly comparable regarding all effects of neighbouring objects.

However, there is one important difference regarding the last point. Instead of fitting neighbouring galaxies (secondaries and tertiaries; see Barden et al. 2012) again, they are modelled with fixed values (i.e. they are merely subtracted from the image as their parameters have already been optimised in the single-Sérsic fit on the same primary target). It is a small, but important, point to note that the values used are not the true best-fit values from when said neighbouring galaxy was the primary target, but from the associated single-Sérsic fits of the same primary target. The profile of the primary object is ‘simply’ replaced by two profiles, a ‘bulge’ and a ‘disk’ (both characterised by Sérsic profiles). Starting values for the primary object are derived from the single-Sérsic best fit in a way discussed in Sect. 2.2. As Galapagos-2 automatically runs the single-Sérsic fits as a first step, a user does not have to take care of such a ‘pre-analysis’.

This scheme has the big advantage that the single-Sérsic fit and the B/D fits are directly comparable, all values other than the ones of the primary object are identical, which makes any further analysis and a possible selection of the ‘better’ fit easier. There are several additional practical advantages of this setup scheme that starts from the single-Sérsic fits.

Firstly, due to its independence of any additional parameters, for example from neighbouring objects, all B/D fits can be carried out independently of each other. No complicated queuing mechanism is needed to avoid simultaneous fitting of galaxies that might influence each other as was the case when fitting single-Sérsic profiles. As a result, fits can simply be carried out one by one in any order, minimising computing time by Galapagos-2 itself. Galapagos-2 strictly goes from brightest to faintest objects in this step.

Secondly, and with much higher impact in crowded fields: as parameters of neighbouring galaxies are fixed, the Galfitm fit is carried out with fewer free parameters, reducing the time needed to carry out the fit itself (although the primary source itself of course shows a higher number of free parameters). In the dataset used in this analysis this advantage is, however, minimal, due to the relatively low density of objects. Galaxies typically have no or only one neighbour in the fit (mean value: 1.37 neighbouring objects per fit, median: 1), but this is expected to make a bigger difference in densely populated areas, reducing the fitting time by a factor of a few compared to the single-Sérsic fits. However, B/D fits introduce one additional profile to each primary object, which slows the fit down somewhat. In the dataset used here, with a median of one neighbouring galaxy, the DoF of a B/D fit is comparable to that in the single-Sérsic fits, leading to a similar total execution time. The exact values, however, depend on the DOF allowed for bulge and disk parameters.

Additionally, there is a second reason why the overall speed of the B/D fits can be faster than the single-Sérsic fits as a whole. In both steps, it is possible to submit a target list to Galapagos - 2 that includes objects that become ‘primary objects’. In the single-Sérsic fits, it is advisable to not only target the galaxies of interest themselves, but also objects within a certain radius around them, but with higher luminosity. The reason for including these additional objects is that those objects would, in a full Galapagos run (without a target list), be treated first. As such, when they become neighbouring objects for a target of interest, their true best fit values are already known and they can be optimally ‘removed’ from the single-Sérsic fit. Not fitting these neighbouring objects first has the potential to bias the results of the objects one is interested in. In the B/D decomposition fits, these objects are not of interest anymore, as their best fit single-Sérsic values are already known. As they are not targets of interests themselves, running a B/D fit on them can be avoided by removing them from the target list. Not fitting those galaxies in the B/D decomposition step hence speeds up the overall fitting process of that step significantly.

2.1.2 Possible improvements

Several possible improvements are obvious to this fitting scheme. As this scheme depends on the existence of a single-Sérsic output file, it means that if such a file does not exist for some reason – for example, the fit might have crashed or the fit has not been attempted as the object was not in the target list – the B/D fit cannot be carried out, Galapagos-2 will simply continue with the next object. However, we have shown in H13 that these cases of crashed single-Sérsic fits are rare as most fits do produce a result.

There is also the danger of biasing the results against certain objects, for example if there is a particular kind of object for which said single-Sérsic fit would often fail or crash. Such a bias, however, is unknown to us, although we cannot exclude that it exists to a small degree. The alternative would be that these objects with failed single-Sérsic fits could run a B/D fit using different starting values, for example directly from running a second single-Sérsic fit with different starting values or from using SEXTRACTOR (Bertin & Arnouts 1996) values for the setup directly. Due to the small number of objects that fail in the single-Sérsic fit and the risks and complications involved in such a scheme, we decided not to implement B/D fitting from SEXTRACTOR values directly. It is hence down to the user to make sure that any target lists provided to Galapagos-2 are compatible: the B/D targets should be equal to or a subset of the single-Sérsic targets to ensure that a B/D fit is possible.

However, Galapagos-2 is friendly in nature. Individual, important cases where the fit crashed can be dealt with by hand. All input files exist, they are easy to identify in the output catalogue, and it is straightforward to manually ensure that a fit for an object does produce a ‘good’ result. once this is done on a single-Sérsic fit, the missing B/D fit on the same objects can be ‘filled in’ by running the B/D part of Galapagos-2 again. Equally, individual B/D fits can be fixed in a similar manual manner. After fixing such crashed fits manually where deemed necessary, Galapagos-2 can simply be re-run to read out the all fit results as usual.

A second possible improvement would be if all neighbouring galaxies were always taken into account with their ‘real’ fitting values (i.e. from the fit in which they were the primary target themselves), as those are the optimal values for them. This is already true for the neighbours that are brighter than the primary target, as this is how they are dealt with in the single-Sérsic fit in the first place. For the galaxies fainter than the primary target, the improvement of such a scheme would be minimal and would make the single-Sérsic and the B/D fits less comparable, so we did not implement such a system in Galapagos-2. Users are reminded that neighbouring galaxies are simply fit to improve the fit of the primary target, their fit values are not written into the output catalogue.

Similarly, one could obviously fit the secondary objects as multi-component systems to improve their fits. However, this is expected to have very minor effect on the fitting results of the primary source as the important part in the fitting of neighbours – the successful subtraction of the profile wings – can be achieved by a simple single-Sérsic fit. The possible improvements are therefore expected to be small, while introducing a large amount of book-keeping. We hence decided against using this more complex scheme.

2.2 Starting values

Deriving good starting values for all fitting parameters in Galfitm is both difficult and important. This is not very critical in single-Sérsic fits, where we have already shown that even rough estimates, for example on galaxy magnitudes (via a typical galaxy SED), are sufficient for a successful fit, but it becomes important when carrying out B/D decompositions, due to the degeneracies that can be found between different profile parameters (e.g. Lange et al. 2016). The more accurate the starting values are, the more likely one would expect the fit to converge to a physical solution rather than settling in a local minimum. However, getting accurate starting values is difficult in a fully automated code, as it should take care of all eventualities and possible galaxy types. Galapagos-2 has to employ a system that works well in most cases. Below, we explain the behaviour of the code in case of the different parameters.

Before we discuss the general approach of Galapagos-2 to derive starting values for the B/D fits, we need to mention a specific case of how Galfitm handles the starting values it is given, as it explains some of our choices in Galapagos-2. The degrees of variation in the parameters for each image/wavelength are defined by the user as Chebyshev polynomials, and the input values can be provided as the parameter values for each band (hereafter ‘band’ values) or as the underlying Chebyshev parameters. While it is irrelevant which input version is chosen, we find band values more practical for interaction and readability purposes and Galapagos-2 makes use of this notation. Galfitm behaves as naively expected when a degree of freedom DOF = 0 is used for the fit in any parameter, in that the parameter values are simply fixed at the input values at each band, no matter their values or shape of the chosen polynomial. However, depending on the way the input values are being chosen, this is not necessarily true when another DOF is being used. In case of DOF ≥ 1 and giving the parameter values as band values (rather than the Chebyshev parameters themselves), the parameter values given might not resemble the desired polynomial degree or shape (i.e. it is not possible to fit the e.g. nine points with polynomials of a lower DOF). When DOF ≥ 2, Galfitm first fits a polynomial to the input values internally and starts the fit from this polynomial, which delivers the desired behaviour.

The behaviour, however, is different for the specific case of DOF == 1. Here, the fit is technically a fit that allows the parameter to vary, but the values should be the same at all wavelengths. As such, it would be possible that Galfitm simply calculates an average of the input values and uses this as a starting value at each wavelength, which would resemble the behaviour at higher DOF. However, we chose that in this case, while Galfitm is allowed to vary the parameter as a whole, the offsets between the values for each image (i.e. the shape of the polynomial as defined in the input parameters) are kept frozen as specified in the input values. In other words, when DOF == 1, the shape of the final polynomial will be the same as that defined in the input parameters, but with a systematic offset to higher or lower values as derived by Galfitm. The reason for allowing this behaviour (in Galfitm, not Galapagos-2) is that it makes it possible to assign a certain SED to an object – for example, a supernova or active galactic nuclei where the SED are often reasonably well known – but fitting its overall brightness only. A second case where this behaviour could be useful is to assign positions according to a known offset between images (i.e. all images from one telescope are slightly offset to all images from another telescope), so a consistent profile can be fit to all images simultaneously, while leaving its central position variable during the fit. Ideally, such a case should be taken care of by shifting mismatched images to the same pixel-grid as all other images, and Galapagos-2 strictly requires this, but Galfitm technically allows this.

This decision for the behaviour of Galfitm in case of DOF == 1, however, means that if one does desire to fit a value to be the same at all wavelengths, one has to make sure that the input values are indeed identical for all bands. Galapagos-2 takes care of this issue by homogenising the input values in this case if the DOF of any parameter – as given in the setup file – makes this necessary (specifically DOF == 1).

The following list gives a brief overview over both the behaviour of Galapagos-2 in general and the consequences for this work in particular. Table 1 summarises the starting values used by Galapagos-2, as well as the constraints used (see Sect. 2.3) and the DOF used in this work, which we discuss in Sect. 2.4. Values indexed with ‘B’ are bulge values, values indicated with ‘D’ resemble disk values, with ‘SS’ values of the single-Sérsic fit.

Position (x, y). Galapagos-2 uses the output values of the single-Sérsic fit as the starting values for the position of both components. Both components are started at (and in fact constrained to) the same values (see below). If offsets between images are allowed in the single-Sérsic fit, these will be reflected in the starting values for the B/D fit. Only if DOF = 1, the starting values are homogenised by taking the median value of the values at all bands from the single-Sérsic output. Generally, it is strongly advised to work with micro-registered images, so such a varying position can be avoided and DOF = 0 can be used.

Magnitude (m). The magnitude starting values are also taken straight from the single-Sérsic fits. The code simply divides the flux equally into the two components (m_D = m_B = m_SS + 0.75⁸.). This means, critically, that both components are started at both the same brightness and the same SED. While it is known that generally speaking bulges are red and disks are blue, we did not want to bias the fits in this way as it might cover up interesting objects (e.g. objects with blue bulges and red disks). Any findings of a colour difference and SEDs discussed in this work as found by Galapagos-2 are hence a pure results of the fit itself. Again, only if DOF = 1, the starting values are homogenised by taking the median value of the values at all bands from the single-Sérsic output. We would advise code users to leave the SED of objects as freely variable as possible.

Size (half-light radius; r_e). We decided on a slightly different scheme for half-light radii. Generally speaking bulges are smaller than disks. From the fitting of very nearby galaxies carried out in V14b, we have found the fits to occasionally converge slightly better when the starting values reflect this (i.e. when the starting disk profile is bigger than the bulge profile). Using subsets of relatively bright objects that we re-fit with ~ten different sets of starting values, we have found that – at least at the typical resolution and redshift of the objects used in this work – these do not significantly change the fitting outcome. Qualitatively, the fitting results showed the same trends and –on average – same values, although values for individual galaxies can vary (see Appendix D).

Nonetheless, in this work we decided to start the sizes at r_e,D = 1.2 * median(r_e,SS) and r_e,B = 0.3 * median(r_e,SS), with lower limits on the sizes of 1 and 0.5 pixels, respectively. The general behaviour of Galapagos-2 is such that the fit – independent of the DOF of the individual components and the single-Sérsic fit – will always start at a value that is constant with wavelength (i.e. it is the same at all wavelengths) and can vary as a constant value during the fit. This is for simplicity and reflects the fact that at this stage in the program any changes of the sizes within the individual components would be unknown, as the single-Sérsic fit would not be able to recover or predict those. This automatically takes care of any variations that might be frozen when DOF = 1 is used (see discussion above). If a user chooses DOF ≥ 2, such variations will be allowed during the fit, and any different r_e in the fitting results at different wavelengths would indeed be a result of the fitting, not a bias from starting values.

Sérsic index (n). The Sérsic index n is a somewhat special parameter in that it is possible that a user would want to hold it fixed during the fit at a specific value; for example, several authors have fit disks with n == 1 and bulges with n == 4, instead of using free values. As we did not want to restrict the use of the code to either, the code allows both free and fixed (to specific values) Sérsic indices, depending on the setup specified by the user. If the Sérsic index is not fixed, the code always starts the fit at values constant with wavelength, as was the case with half-light radii. It then uses n_D = median(n_SS) and n_B = median(n_SS), but imposes an upper limit 1.5 for the disk Sérsic index starting value, and a lower limit of 1.5 for the bulge Sérsic index starting value, respectively. The limits used here reflect the general case that bulges show higher Sérsic indices than disks. During the fit, however, these limits are not imposed.

Additionally – and in contrast to any other parameter – it is possible to fix the Sérsic index to n_D = 1 and n_B = 4 for disks and bulges, respectively, if preferred, by simply setting the DOF = −1 of the Sérsic index for the respective component. This allows fitting a true exponential disk and/or a de Vaucouleurs bulge as seen in the classic case. In those cases, the starting value of the parameter becomes irrelevant, and is of course set to the fixed value.

Axis ratio (q). For both components, the values q_SS are taken as a starting value directly (including any possible variations with wavelength, if allowed in the single-Sérsic fit). As usual, if DOF = 1, the starting values are homogenised by taking the median of the values at all bands from the single-Sérsic output. Additionally, a minimum starting value for q_B > 0.6 is introduced to reflect that galaxy bulges are generally round-ish. However, the fit itself is allowed to go below this value.

Position angle (θ). For both components, the value θ_SS is taken as a starting value directly, including any possible variations if allowed in the single-Sérsic fit. The case of DOF = 1 is treated in the usual fashion.

This means that, apart from r_e, n, and possibly q_B (depending on q_SS), the bulge and disk in B/D fits in Galapagos-2 are always started at identical parameter values; specifically, the SEDs of the bulge and disk are identical at the start of the B/D fit. Any deviations in the final magnitudes (and hence the SED of the components) are a result of the MEGAMORPH method to use multiple images at different wavelengths to fit one consistent galaxy profile across all wavelengths, while allowing the magnitudes of each component and at each wavelength to vary freely. The colour information embedded in these images at different wavelengths makes it easier for the fit to distinguish different components, as we show in Sects. 4.3 and 5.4.

A radically different approach to two-component fits has been presented by Lange et al. (2016). Instead of starting all B/D fits at commonly derived starting values, they fit each object multiple times, starting on a fixed grid of starting values (i.e. different starting values for r_e and n are used). They then analyse the convergence of these multiple fits to derive average and more reliable ‘best-fit’ parameters. This method allows one to get a better handle on the robustness of fitting parameters but requires every galaxy to be fit multiple times: 40 fits are carried out for each object. For obvious reasons on sample size and required CPU time, we avoided such an approach and instead, using simulated data, show that our approach allows good results.

2.3 Constraints

Generally speaking, the same constraints are being used in the B/D fit that have already been used during the single-Sérsic fits. These have been discussed in detail in H13 and are only quoted here for completeness with additional comments.

Position (x, y). Positions are constrained to lie within a box of size 0.5 * r_e,PB,SS around the object centre as defined by the single-Sérsic fit⁹. Additionally, the position of the disk and the bulge are constrained to be the same. While in nature, the bulge and the disk can in principle be slightly offset, especially in peculiar objects (e.g. post mergers that have not yet entirely relaxed), this offset should be much smaller than 1 pixel, given the resolution and the distance of the galaxies typically fit with Galapagos-2. On larger, nearby galaxies with visible dust lanes that could affect the disk and bulge differently, such an offset might make sense. However, these galaxies are not ideally dealt with using Galapagos-2, so a more flexible implementation seemed too complicated and not useful. The advantage of constraining the B/D positions to be the same, by avoiding many issues with clumpy galaxies or neighbouring objects, outweighs the disadvantages and limitations introduces by such a choice, as discussed in Sect. 3.

Magnitude (m). The constraint on magnitudes is a user-specified value, we use 5 ≤ m_D/B,fit – m_D/B,input ≤ 5 in this work for each band. Such a wide ± 5 magnitude offset has to be allowed as the real brightness of bulge and disk, respectively, are unknown. This limit generally resembles a limit well beyond the brightness or faintness of a component that one would still trust in a fit, so does not impose any significant effect on the output values that should be used in an analysis.

Additionally, we use the same basic constraint of 0 ≤ m ≤ 40 that has been used in the single-Sérsic fitting already. Other than most constraints discussed here, this is currently a limit hard coded into Galapagos-2 (although trivial to change), as it easily covers all current galaxy surveys. Given the above constraint, however, this will only be violated in un-physical fits and catastrophically failed fits, and is only mentioned here for completeness.

Size (half-light radius; r_e). The upper value is a user-specified value. We use r_e,D/B ≤ 400 px throughout this work (i.e. the same as in single-Sérsic fitting) mainly to prevent the fit from returning unphysical results. A lower constraint of 0.3 px ≤ r_e,D/B is currently hard coded into Galapagos-2.

Sérsic index (n). Both upper and lower limits are user defined values in Galapagos-2. In this work, we use 0.2 ≤ n ≤ 8 (i.e. the same as in single-Sérsic fitting) but for each component individually. We note that n_D == 1 in this work (see Sect. 2.4), so this constraint only has an effect on n_B during the B/D fit examined here.

Axis ratio (q). 0.0001 ≤ q ≤ 1, for the reasons given in H13.

Position angle (θ). −180deg < θ < 180deg, in order to prevent numerically different but otherwise identical fits with unnecessarily large values for θ. No other constraints are set: the bulge and disk are allowed to have different θ values.

Neighbouring objects are used identically to the single-Sérsic fits and with fixed values, so no constraints have to be employed for these objects.

2.4 Choices of setup and degrees of freedom

In single-Sérsic fitting, we allowed variation in certain parameters – size (r_e) and Sérsic index (n) – with wavelength and have exploited this in detail in H13, V13, and V14b. This decision effectively allowed colour gradients within galaxies. These gradients are largely a result of different mixing of the different stellar populations in the galaxy bulge (generally older and redder) and galaxy disk (generally younger and bluer) at different radii (Kennedy et al. 2016a).

While, in principle, such gradients could also exist within individual galaxy components, the colour differences are expected to be small compared to the colour differences seen between the two components as a whole, at least in the vast majority of cases, making them very hard to measure. As a result, we decided to not allow any wavelength variation in any component parameters other than magnitude in this work (magnitude is a free parameter in all bands, as in our previous work). We are aware that this is only a first step towards understanding galaxy components in detail, but we find it necessary to understand this step in detail before taking the next step, for example allowing colour gradients within galaxy components in more suitable datasets with brighter and better resolved galaxies (e.g. the ones analysed in V14b).

However, both codes, Galfitm and Galapagos-2, are written in a way that the user can trivially switch on wavelength variations in the setup files, if they so wished. In order to support our choice, we ran B/D decomposition with linear variation allowed for r_e of disks and bulges for a subset of 4000 galaxies, and have found that their values in fact do change slightly with wavelength, but by much smaller factors than reported by Vulcani et al. (2014) for single-Sérsic fits and, in fact, within the measurement uncertainty in most cases. We have found that r_e,D and r_e,B change by ~14 and ~23%, respectively, from the u band to the K band. Similarly, we fitted the 163 nearby galaxies analysed in the V13 sample with linear variation in disk and bulge size (but not Sérsic index) and have found little variation from the u band to the z band covered by these data, with r_e,D changing by ~−10% and r_e,B by ~+10%, the exact value also depending on galaxy type. We conclude from this that keeping these parameters constant with wavelength is a justifiable assumption in case of B/D decomposition of GAMA galaxies. In the following, we discuss the choices for each parameter.

Position. As images at different observed bands are accurately registered (although see Kelvin et al. 2012 and Sect. 4.1), we kept the position constant with wavelength in the B/D fits.

Magnitude. As the polynomial shape of the SED of each component is unknown, we kept the magnitudes with full degree of freedom for both bulge and disk. Users should be reminded of Runge’s phenomenon, whereby a polynomial function can oscillate excessively between data points, particularly at the edges of the considered interval if the degree of freedom is similar to the number of points fit. As such, it is in theory dangerous to use the Chebyshev polynomial of the fit directly to derive values at intermediate wavelengths. While Nedkova et al. (in prep.) have found that such oscillations are rare in their fits on Candels data, great care has to be taken when interpolating the magnitudes.

Sizes. For both component, bulge and disk, we used a size that is constant with wavelength, for the reasons given above.

Sérsic index. We used n_D = 1 (fixed) and a n_B that is constant as a function of wavelength (but variable during the fit). This is also what we used for the simulated galaxies in Sect. 4.1 (see Appendix C), so no biases are introduced in these tests.

Axis ratio. Axis ratios were chosen to be wavelength independent for each component but can vary during the fit.

Position angle. Position angles were also chosen to be wavelength independent (i.e. no rotation is allowed between images at different wavelengths). The value can vary during the fit.

List of target objects. As was mentioned before, the list of target objects for these two-component fits should be equal or a subset of the targets used in the single-Sérsic step of the code. For our runs, we made sure that this was the case. In case of the simulated images, we simply fit all objects at both stages. When fitting real galaxies in Sect. 5, we target all GAMA galaxies at this stage. In the single-Sérsic fits, we had also targeted bright neighbours to these objects in order to be able to take these properly into account when they become neighbouring objects for a fit. Given how these neighbours are being dealt with in the B/D fits, it is no longer necessary to deal with these objects.

3 Data and rejection of bad fits

As usual, before any analysis the catalogues created by Galapagos-2 have to be cleaned of bad fits. In this chapter, we explain only the principles used to define the samples presented in the following chapters. As catalogues used are by design very different (simulated vs. real data), we leave the discussion of the exact samples and object numbers to the respective sections, Sects. 4.2 and 5.3.

To identify ‘good’ fits, we use the same parameter limits used in H13. More precisely, we identify and discard those components with one or more parameters lying on (or very close to) a fitting constraint, as described in Sect. 2.3. Such a fit is unlikely to have found a global minimum in χ² space, but is rather constraint by the shape of χ² space along the boundary box of the allowed parameter space, and is a good indication that the parameter values given cannot be trusted. While GALFIT and Galfitm have ways to try to avoid such local minima in the final fit solution, these measures do not work in extreme cases of very deep local minima. The absolute χ² (or reduced χ², χ²/ν) value itself is also not a good indicator on whether a fit is a good representation of the true galaxy profile, as it depends strongly on the image properties, especially the precise way the pixel-to-pixel noise is correlated. It is also artificially increased in the case that the object of interest cannot precisely be modelled with the model used in the fit. Real galaxies often have features (bars, spiral arms, faint components, merger features) that prevent a perfect fit, and the contribution of neighbouring objects also adds to the absolute χ² value, bright stars being a particularly bad example, so such a mismatch is obvious in real galaxies. We hence avoid using absolute χ² values in this analysis as a measure to define good fits. Instead, we avoid any object where we already know that the fit was not a free fit.

For the analysis throughout this paper, we only keep objects that fulfil the following criteria:

− m_input − 5 < m < m_input + 5,

− abs(x_pos − x_pos,SS) < 0.5 * re_SS

− abs(y_pos − y_pos,SS) < 0.5 * re_SS

− 0 < m < 40,

− 0.205 < n < 7.95,

− 0.305 [pix] < r_e < 395.0 [pix],

− 0.001 < q ≤ 1.0,

where the magnitude input values m_input are derived from the single-Sérsic fit result as described in Sect. 2.2. For the individual single-band fits, this is obvious, but for the multi-band fits, these constraints are equally checked on each band, removing the component from the subsequent analysis if even a single one of these parameters fails this test in any band. The limits on position make sure that the fit stayed on the primary target and did not run away to fit a neighbouring object. As position of bulge and disk are by design fixed to be the same, this limit is in practice never violated (e.g. by the disk fitting one target and the bulge fitting a neighbour). The criteria for n and r_e are slightly more restrictive than the fitting constraints used. When we use limits more tightly around the constraint values (0.201 < n < 7.99 and 0.301 [pix] < r_e < 399.0 [pix]), the number of successful fits in Table 2 are naturally higher as fewer galaxies are excluded (i.e. both bulges and disks in multi-band fitting are 98% ‘successful’). While this might seem like a preferable cutoff to use, it merely changes the definition of ‘success’ and would, in fact, include many galaxies where the fit was in practice not a free fit. This indicates that even values close to a fitting constraint (rather than on the limit precisely) are potentially not the result of a free fit. The numbers in Table 3, however, would not change much, which shows that especially faint objects are removed by these cutoffs, which seems obvious as these objects will be hardest to fit. These are hence the most critical limits we use to clean the galaxy samples and a balance between rigorously avoiding bad fits while allowing large samples must be found.

The obvious difference of this work compared to H13 is that these limits need to be checked for each component individually, instead of each galaxy as a whole. The reason for this is simple: If a faint galaxy component sits within a bright component (e.g. faint bulge within a bright disk or vice versa), it is irrelevant whether the fainter component is fit well in order to decide whether we would believe the values of the brighter component. We might believe the fit values of a disk even if we do not believe the fit result of the bulge in the same galaxy. This can easily be seen in the numbers given in Sects. 4.2 and 5.3, where the bulge and disk samples contain different numbers of objects and do not fully overlap. For such a scheme, it is also wise to include a further brightness limit below which one would not believe the parameters of the fainter component, despite it being ‘well’ fit given the above limits. During this work, we chose to use the fainter component only if it is within 1.5 magnitudes of the brighter component in the r band, which was chosen as the main band in this analysis and within Galapagos-2 in our setup. This value was somewhat empirically established and will be discussed in Sect. 4.2. Realistically, such a choice could additionally be a function of magnitude of the galaxy as a whole (i.e. in brighter objects a larger magnitude difference between the components might be acceptable). Such a scheme is hard to define universally, however, and is beyond the scope of this paper.

4 Application to simulated imaging

The advantages and disadvantages of using simulated versus real data have already been discussed in H13. In order to follow the same approach, we decided to, once again, use both simulated and real data in our analysis. Simulated data have the advantage that one knows the input galaxy parameters exactly, allowing for a direct comparison of recovered galaxy parameters to their true values while assuming a perfect profile match (i.e. in that galaxy components actually are well described by Sérsic profiles). For these tests, the choice of simulated parameters is obviously an important one. If the simulated galaxies do not resemble typical B/D systems, not much can be learned from this analysis. If, on the other hand, parameters are wisely chosen and resemble real objects (and potentially include interesting objects one would like to identify in the real data, e.g. galaxies in which the bulge is larger than the disk or with inverted colours), one can conclude that the parameters found in real galaxies are likely to be well measured and a good description of the galaxies found in nature. As in H13 and similar work (e.g. Meert et al. 2015), simulations are used as an idealised case here and, while allowing a good comparison of single-band to multi-band fitting, they can only give a lower limit on the fitting uncertainties for real data.

In this section we describe the galaxy simulations and the choice of parameters in detail. We then present and discuss our finding when using the galaxy images to test the performance of Galapagos-2 and Galfitm in Sect. 4.3.

4.1 Creating the simulations

The simulations used in this work are created using the same codes/methods used and described in previous papers, so we refer the reader to Sect. 5.1 of H13 and Sect. 2 of (Häussler et al. 2007, hereafter H07) for technical details on how the images are created in detail.

In H07 we presented analysis and testing for single-component, single-band datasets. In H13 we extended this analysis to multi-band fits but still restricted ourselves to single-component objects. In this work, instead of creating single-Sérsic objects, we aim to simulate galaxies with several components. To achieve this, we simply put two objects at the same position, using the same scripts but with different parameters, one generally representing a bulge and one generally representing a disk. As the details of the simulated data might be important to some readers, we give extensive details in Appendix C and refer to that section for the details in these simulated data. The three most important points in the resulting objects are as follows.

First, we simulated fewer faint objects compared to H13, where we used galaxies up to 4 magnitudes below the peak in the magnitude histogram in a size bin. Their main purpose in H13 was to push the single-Sérsic fits to their limits on galaxies along the detection limit. As those objects are extremely unlikely to be decomposed by any code, we restricted the objects simulated here to 1 magnitude below the peak in the magnitude distribution. In the following sections we further restrict our analysis to objects at m_r,B+D < 19.5 only, which would avoid analysing fainter objects. This imposed magnitude limit somewhat matches the magnitude limit of the GAMA survey, which provides spectroscopic redshifts down to a magnitude of m_r < 19.8. While results presented here still hold at this fainter limit qualitatively, it is at significantly lower significance, which is why we chose a somewhat brighter magnitude limit to present our results. Objects at m_r ~ 19.8 are very much on the edge of what could possibly be decomposed in these SDSS/UKIDSS data and the samples derived from overlap with single-band fits are very small.

Second, all disks and all bulges show the same, somewhat extreme, SED, with disks being bluer and bulges being redder. This choice was taken in order to make it easier to analyse and present the results as they are likely colour dependent, disk and bulge with similar SEDs being much harder to decompose. However, initial tests with more general data showed the same trends, but with less statistical significance. Some – 0.1 mag – noise is added to these SEDs, in order to simulate some variation found in real galaxies, however.

Finally, bulge Sérsic indices (n_B) show a variety of values, centred around 4 but with a wide spread. This allows bulges in the simulated data to explicitly not be classical de Vaucouleurs bulges.

The resulting images are realistic looking at all wavelengths with thousands of galaxies that span a large range of parameters, but – overall – show similar parameter distributions to real galaxies with realistic noise properties, and for which we know the true parameter values, including – and especially – their subcomponents.

It should be noted that these simulated images do not contain any stars. While this removes a potential source of error from the analysis, this source has already been tested in previous work (e.g. H07) and should not influence the B/D decompositions more strongly than single-Sérsic fits. Overlapping galaxies are naturally included by our simulation method, and so the effects of blending several objects are still included in our results.

The resulting images are then fed through Galapagos-2 and Galfitm as described in Sect. 2.1, that is, the same pipeline is used for both the simulated and the real data (discussed in Sect. 5.4).

4.2 Samples

In total, 95 143 galaxies have been simulated in the survey area analysed in this work. 75 998 of these objects were simulated as two-component galaxies, the remaining ~20% were one-component systems. While a detailed analysis of these objects, and an attempt to identify them in an automated fashion within our dataset, are beyond the scope as this paper, we had a quick look at the two-component fits of these one-component objects.

These objects seem to often fall into one of three qualitative categories. The objects in the first category are objects in which the fits for one of the components runs into a fitting constraint. For these objects, this is very often the Sérsic index of the bulge, which hits the upper limit of n_B == 8, often in combination with fitting a very small size of the bulge. As fitting constraints are violated, these objects would be removed from the analysis when cleaning the catalogues for a two-component analysis. As the violating component can also be faint, there is a big overlap of this category to the next.

The second category contains objects for which Galfitm returns a large magnitude difference between the components. This behaviour is what one would naively expect in that one component fits the galaxy profile well (which component this is depends on the overall Sérsic index of the galaxy), while the other tries to fit a small correction to this profile, for example even by trying to fit a group of high-flux pixels in the noise pattern. These objects can be identified by a large magnitude difference between the two components in the B/D fit, and often a very small size in the fainter component. This is why we introduce a limit on the magnitude difference in our analysis. We only use the fainter component if it is not more than 1.5 mag fainter than the brighter component, as measured in the r band.

The third category seems to contain objects for which the fit behaves such that both profiles mimic each other and the overall profile. For many of these objects, we find magnitudes (and hence SEDs) of both components to be very similar, as well as similar sizes of the components. As the B/D fits are started at r_e,D = 1.2 * median(r_e,SS)(>1 px) and r_e,B = 0.3 * median(r_e,SS)(>0.5 px), this means that the sizes in the fit actually converged to be the same value as one would expect if they fit the same profile. Effectively, the flux of the galaxy is simply divided into two components, without any physically meaningful separation, making the SEDs of these components largely unconstrained, so it comes as somewhat as a surprise that the flux often seems to be shared equally by the two components. Unfortunately, this does not always seem to coincide with n_B ~ 1 as one would naively expect as these are the galaxies where the profiles would mimic each other best. The galaxies in this category are the hardest to find in large datasets, but we discuss such an effect in Sect. 4.3.

There seems to be a big overlap between all these categories, especially between the first and second. It should be noted that we only mention these general categories here for completeness. While these categories can be found in our analysis, they are not pronounced enough to actually use them to separate out the single-Sérsic objects from the two-component galaxies as they are, and additional development and testing would be required for this purpose. Such an attempt, however, will have to use a more sophisticated method to separate the object classes and is beyond the scope of this paper. For example, it has been found that the classification scheme presented by Allen et al. (2006) works well on Candels data (e.g. Nedkova et al. 2021; Nedkova, in prep.) but it is unclear whether the same scheme would work equally well on the GAMA data used here.

In the following, we restrict our analysis to the two-component galaxies. of these, not all galaxies are recovered by SEXTRACTOR, some are too faint to be detected. Given our analysis limit of m_r < 19.5, these missed galaxies are unlikely to be presented in this work. Depending on the band used for object detection, this fraction is very different, the detection and fitting numbers can be found in Table 2. At best – when running detection in a band-combined/multi-colour stacked image −~70% (53 448/75 998) of the objects can be recovered. This detection completeness is not part of the analysis in this work and will hence not be discussed here. The important part for this work is that we can analyse galaxies all the way down to the detection limit. In this work, we merely analyse what fractions we can successfully fit.

The object numbers given in Tables 2 and 3 show how much multi-band fitting improves the sample size for scientific studies, especially once measurements from several bands are required. Each row shows the numbers of objects in each – single-band and multi-band – Galapagos-2 run as well as the percentages of objects that deliver a successful fit for single-Sérsic (on the left) and B/D fits (on the right) in the event the B/D fits split for bulges and disks. Additionally, we show the object numbers and success rates when combining several single-band fits (in black), for example in case a science case requires values from more than one band (ugrizYJHK and griYHK, respectively), which drastically reduces the available sample size that we would recommend to use. These numbers also include the cut as the fainter component being within 1.5 magnitudes of the brighter component, based on the fitting values. In light grey, we give the same values when using simulated values to make this decision. As one can see, multi-band fits are far more likely to produce a valid fit result. For all galaxies at m_r < 19.8 (GAMA spectroscopy limit, Table 3), only 2369 objects (14% of objects of the 16941 objects for which a B/D decomposition was attempted) have good single-band fit in all griYHK bands in both components (366, 2.2%, in case one adds uzJ bands), immediately reducing the sample size for any analysis that requires these parameters, for example magnitudes in several bands when an SED of a component is required. In comparison, 14 732 (85.5%) of objects have good multi-band fits, increasing a possible science sample by more than a factor of 6. In the following sections, in order to estimate the effect of multi- versus single-band fits, we try to use the largest sample possible in most cases to maximise the statistical significance of our findings and allow a fair comparison of the two different fit performances (i.e. on the same objects). Unfortunately that drastically reduces the sample size available in most plots.

In simulated data, it is possible to define galaxy samples using two different methods: using measured values or using simulated values, for example when identifying the B/T ratio of a galaxy (recovered B/T vs. true B/T). The former would allow a more direct comparison to real, observed objects, as the selection could be entirely identical. The latter defines slightly different samples, but allows a cleaner comparison to simulated values. In this paper, we decided to define galaxy samples using the simulated (and true) values where possible, as the comparison to those is the main purpose of this paper. However, we carried out the same analysis using the measured (and recovered) values instead, which qualitatively leads to the same results. However, plots are somewhat harder to read using this approach as effects can be less pronounced and different effects are harder to separate from each other. The conclusions in this paper have not been significantly altered by this choice. In Tables 2 and 3, we give numbers for both definitions. From these numbers one can see that the sample size increases when using observed values. This indicates that the fainter component in a galaxy often ends up accounting for some flux from the brighter component, boosting its magnitude and hence pushing it into an observed sample.

In Sect. 3 we mention that we only believe the fainter of the two galaxy components if it is less than 1.5 magnitudes fainter than the brighter component (δm_r < 1.5). The reason for this becomes obvious in Figs. 1–5. Figure 1 shows the deviation of the disk r-band magnitude as a function of the difference between disk and bulge magnitude (i.e. B/T ratio) for the ~11 000 galaxies that have successful fits in both multi-band and r-band fits, with a running mean and scatter over-plotted as blue lines. Somewhat arbitrarily, we define ± 0.5 mag deviation from the true value as acceptable (red horizontal lines). A clear trend is visible such that faint disks are badly fit, showing systematically brighter fit values. This is understandable as they are embedded in a much brighter component and the fit compensates fitting residuals of this component. The ± 0.5 mag is reached when the disk is more than 1.5 magnitudes fainter than the bulge within a galaxy (blue median line crosses the red vertical line), slightly more for brighter galaxies (orange median line). The red numbers in the top-left corner indicate the numbers and fraction of objects with deviations larger than 0.5 mag above (top number) and below (lower number). However, it should be noted that these numbers indicating the number of outliers should be taken with a grain of salt, as it is obvious that they are dominated by fits that we already know are bad and refer to relatively small deviations, which explains the high fraction of objects. For comparison, in the bottom panel of this figure, we show the same plot for single-band fitting in the r band. There are a somewhat larger number of outliers (as indicated also by the numbers), and systematic deviations are typically a little larger in this case. The 0.5 magnitude deviation is reached at about 1 magnitude difference between bulge and disk. This exercise can obviously be made in any of the bands used, and it is important to point out that the r band is the deepest and best behaved of the single-band fits, some of the other band looking significantly different, especially the bands with shallower images, for example uzJ. We chose the r band here, as it serves as our main band throughout this work.

Figure 2 shows the deviation in the fits to the simulated bulge magnitudes in the same way as Fig. 1. The fraction of outliers is higher here as would be expected for bulges, and we find the same value of 1.5 magnitude difference between bulges and disks, above which their values are reliably fit. This limit restricts the objects for which we believe the fitting values for both components and targets to be 0.2 < B/T < 0.8, which incidentally are values also found in the literature (i.e. objects with B/T > 0.8 are often called ‘pure spheroids’). There are other trends visible in these plots, which need to be explained. Firstly, it can be seen that very bright (and hence ‘pure’) components, both in the case of the disk and the bulge, are somewhat under-fit (i.e. the fit returns a fainter magnitude value). This can be somewhat understood with the same argument used above. If the fainter component takes away some of the flux of the brighter component, the brighter component will need to compensate and hence return a fainter magnitude. Secondly, there are a large number of outliers in this plot. Points and blue lines show all these galaxies with m_r < 19.8 to match the science sample discussed in Sect. 5. When restricting the sample to m_r < 18 (mean value and scatter shown as orange lines in all plots), both under-fitting and extreme outliers appear less frequently and the plots appear cleaner. For these bright objects, a fainter limit of δm_r ~ 2 mag could in principle be used in multi-band fits, but it is important to point out that several of the brightness bins contain only a few objects. However, for multi-band fits, this improvement is larger than for single-band fits. Additionally, as bulges have been simulated with a range of Sérsic indices, a trend with bulge Sérsic index is to be expected. High-n objects are generally harder to fit (see H07 and H13), so one would assume the same to be true for bulges embedded in a different object. Indeed, when restricting the plot to high-n objects, trends are emphasised in all bulge plots, in that offsets are generally larger for high-n bulges.

Figures 3 and 4 show the deviations in the fits to the disk and bulge sizes, respectively, and show that neither multi-band nor single-band fits display large systematic trends for galaxy disks. Deviations are somewhat more pronounced in single-band fitting of disks, but the bigger difference can again be seen in the scatter of the distribution, especially in bright objects. For bulges in Fig. 4, the same effects can be seen, although both fits seem to fit ‘pure’ spheroids too small, which is true even in bright object. Even a limit of more than ~ 1.5 mag would work for both components when large, bright samples are being used.

For the bulge Sérsic indices, however, this is not true, as can be seen in Fig. 5. Both multi-band and single-band fits have trouble recovering the Sérsic index if the bulge is too faint, as one would also naively expect. In multi-band fits, ~20% accuracy is reached at most δm_r values and about constant until δm_r· = 0 (i.e. B/T = 0.5). However, at δm_r· = 0 mag the single-band Sérsic indices are already badly underestimated, by about 50%. In general, fainter bulges become more and more disk-like, as it becomes more and more likely that the bulge profile fits small residuals of the galaxy disk instead of the actual bulge. There is also an effect at very bright bulges in that Sérsic indices are already ~20% underestimated. This effect is enhanced to ~30% in single-band fits. We confirm these trends in Sect. 4.3.6. Again, these trends weaken significantly when the samples are restricted to brighter galaxy samples (orange lines), where the S/N in the individual components is also higher and multi-band fits can still recover n_B at ≲ 10% at δm_r = 1.5 (~20% for single-band fits). As was expected, and as has been reported by other authors before, parameters of high-n objects indeed seem harder to recover accurately. The Sérsic index is the hardest parameter to fit, but its known degeneracies with galaxy sizes and magnitudes mean that those parameters are also less accurate than in the case of galaxy disks.

Overall, we can establish from these plots that a magnitude limit for the fainter component has to be used and values are generally recovered within 20% (and 0.5mag) when limiting the samples to those components with δm_r < 1.5 mag.

4.3 Results

In this section we discuss the results from the simulations and establish how well we can recover the B/D parameters that we put into the simulated images using both single-band and multi-band fitting. Throughout this section, unless stated otherwise or trivially visible from the axis range, we restrict ourselves to show results of galaxies at m_r < 19.5 in order to match the sample of real galaxies (see Sect. 5.4) for which GAMA provides spectroscopic redshifts at m_r < 19.8, and to avoid very faint galaxies where analysis suffers badly from low-number statistics of successful fits in several single-band fits. Where a different sample cut is used, this will be indicated. Many of our conclusions still hold qualitatively for fainter galaxies: the advantages of multi-band over single-band in recovering the profile parameters and of the sample size being available for a scientific analysis (see numbers in Tables 2 and 3) are expected to become even more pronounced. However, we do not present them in this work in detail. We highly recommend to any user of Galapagos-2 and/or Galfitm to run similar tests on their own data in order to establish the accuracies and issues, especially on faint galaxies.

We also only examine the ‘pure’ fitting modes and what we called Mode_M (multi-band fitting with object detection in a co-added image) and Mode_Sl (single-band fitting with object detection in the according single-band image) in H13. We do not discuss Mode_S2 fitting (single-band fitting with object detection in a co-added image) as this was already shown to be counterproductive since many galaxies are detected in the co-added (and hence deeper) image but are too faint to be sensibly fit in single-band images, confusing the single-band fits even of the visible targets.

4.3.1 Component magnitudes

In Fig. 6 we show how well magnitudes can be recovered using single-band and multi-band fitting for bulges (top) and disks (bottom), respectively. The left plots show the analysis of bright galaxies at m_r,B+D,sim < 17.5 mag from their simulated magnitudes, the right column for faint objects at 18.5 < m_r,B+D,sim < 19.5 mag. This forms the most basic magnitude comparison for each of the components and is equivalent to the way results were presented in H13, where error bars represented symmetric standard deviations. In this work, we use 16 and 84 percentiles in most plots to indicate not only the errors of the fit, but simultaneously allow an estimation of how symmetric the distributions are.

Several effects are visible in this plot. Firstly, for most bands, it becomes apparent that multi-band indeed does improve the recovery of magnitude values, as expected. The user should be reminded that the magnitudes in our multi-band fitting are not constrained by any assumption on a polynomial; the improvement seen here is a result of constraining the profile shape itself, that is, the other fit parameters follow certain functional forms. We already reported in H13 that the improvements of fitting performance on magnitudes is not as strong as on other parameters for this reason. Secondly, it becomes apparent how asymmetric some of the distributions are. Both single-band and multi-band fits tend to fit disks too faint, and bulges too bright, but this effect is much more pronounced in the single-band fits, as is visible from the larger, and asymmetric dispersion of the objects. While neither of the methods show large systematic offsets for these bright galaxies (left) and recover bulge and disk magnitudes that are – on average – nearly perfectly (offsets generally below 0.05 mag), the multi-band fits recover the magnitudes of bulges and disks more accurately at all wavelengths, showing smaller error bars, highlighting the advantage of multi-band fits. This is especially true for the bulge component, which is known to be harder to fit.

In the right column of Fig. 6, we show the same plot, but for fainter galaxies with 18.5 < m_r,B+D,sim < 19.5 mag. As expected, at these magnitudes fitting uncertainties are in general much larger, indicated by the larger error bars. It can also be seen that especially single-band fitting returns values with significant systematic offsets. Systematic offsets in multi-band fits are much smaller where present, with typical offsets of the bulge or disk magnitudes below ~0.1 mag, but larger dispersion. There is, however, a significant effect in the u band of the bulge magnitudes (top-right panel), where a substantial systematic offset can be measured for both multi-band and single-band fits. We emphasise here that the very small number of bulges with good fits in galaxies at these magnitudes and number statistics are a problem. This is largely due to the u-band single-band fit itself being unsuccessful in providing a fit value, as the image is much shallower compared to the images in other bands (compare detection numbers in Table 2). Given the SEDs used in the simulations, the B/T ratio in the u band is also significantly smaller than the B/T ratio in the r band, which we use to define out samples, increasing the difficulty to derive a good value and resulting in many of the bulges of these galaxies indeed being so faint that neither single-band nor multi-band fits can recover their magnitudes accurately. The opposite effect, albeit smaller, can be seen in the disk magnitudes (bottom-right panel). As the offsets seen in the single-band fits go in different directions in different bands, they become especially critical when using colours in a scientific analysis, as colours derived from single-band fitting would be especially affected. Disks are generally still fit too faint, bulges too bright (exception: u band), and this effect is strongly enhanced for these faint objects. This has a catastrophic effect when B/T ratios in a specific band are required (See section 4.3.3 for further discussion).

We can develop a more detailed picture of recovering the disk and bulge magnitudes in Fig. 7, where we show the deviations of the respective values from the simulated values as a function of (B/T)_sim and m_r,sim for multi-band fits, in order to understand their performance in more detail. In each bin of this 2D plane, we show how well the component magnitude is recovered on average via a 2 sigma clipped robust median, with larger squares reflecting larger differences between the magnitudes derived from the fits and those used to simulate the galaxies. We note that in order to ensure reasonable number statistics in each bin, we only plot bins that contain more than 20 objects. We show the scatter in the distributions in the respective bottom panel. In each panel, the hashed areas indicate the area in which we would not trust the fit of the fainter components, and we consider objects pure disks and pure spheroids. In these areas, we would suggest the use of the single-Sérsic fit values instead, as we argue in the following sections. We show below why even the fit to the brighter component is not ideal in the B/D fits in these areas.

The first effect visible in these plots is that, as expected, the fainter component in the extreme areas can indeed not be recovered well. Disks (left panels) are recovered with large systematic offsets in objects with high B/T ratios, bulges show especially large systematic offsets in objects with low B/T ratios. In both cases, the fainter component is recovered too bright by more than 1 magnitude, 3 magnitudes in extreme cases. The same effect can be seen in the lower panels that show the scatter of the recovered values, the values are highest in the area where the component is faintest. This shows that indeed the fits of these fainter components are not accurate and reliable. Similar to the analysis in Sect. 4.2, these plots suggest that the choice of ± 1.5 mag between the components is a good choice, corresponding to B/T < 0.2 and B/T > 0.8. In between these values, both bulge and disks magnitude seem to be recovered relatively well. There is also a (weak) trend with galaxy brightness, in that fainter galaxies are harder to separate, as one would expect. This is especially true in the case of galaxy disks and will be more pronounced in other parameters.

However, there is a second effect, which is somewhat unexpected at first glance. Disks in pure disk objects and bulges in pure spheroidal objects are also fit less well, and show systematic offsets to fainter magnitudes. This might come as somewhat of a surprise, as one would naively assume that, for example, a disk in a pure disk galaxy should be easy to fit. However, this effect can be explained by the fitting process itself and is, in fact, a consequence of the first effect discussed above. We are forcing a B/D fit on an object that is effectively a one-component model with a smooth Sérsic profile and no internal structure. With the constraints on object positions and other factors, this makes it likely that the fainter component fits some flux of the brighter component. This is indeed the effect described and seen above in that disks in pure spheroidal galaxies are overestimated in brightness. However, as a result, in the bulge, this flux is missing such that the bulge brightness would be underestimated, which is the effect seen here. This forced two-component fit also affects the galaxy size and Sérsic index, which we discuss further in Figs. 16 and 20.

4.3.2 Total magnitudes

While in these objects with extreme B/T ratios the flux of the brighter component is underestimated and that of the fainter component overestimated, these two effects cancel each other out nearly completely, as we analyse in Fig. 8, where we show how well both single Sérsic profiles and B/D fits in both single- and multi-band fits can recover the total magnitude of a galaxy. For each band, we show a group of four values, which are most visible in the reddest bands as they are more spread out. The inner two points (orange and blue) show the recovery of single- Sérsic fits of the total magnitude of a B/D object in each band using multi- and single-band fits, respectively. The outer two (red and green) data points show the same values for the B/D fits in multi- and single-band fits, respectively. For bright objects (upper panel) all four methods measure the total magnitude of a galaxy accurately (we note the different ranges on the y axis in comparison to the previous plots) and with small dispersion. This behaviour can be understood by the nature of fitting techniques, in that any residual from one component would be minimised by the other component, balancing out any error in the estimation of the total flux. However, while offsets in recovering the total magnitude are in general small, they are minimised in multi-band fitting, especially in the low S/N bands, for example u. The difference in the scatter – multi-band showing asymmetric error bars towards brighter magnitudes, while single-band show more symmetric errors – is interesting, but has minimal effect on the median values. For faint galaxies (bottom panel) this effect is enhanced and one can see that singleband fits show larger systematic offsets and larger dispersion, especially in case of B/D fits (we note the low-number statistics). Multi-band fits recover total magnitudes more accurately, although small systematic offsets can still be seen.

It is also visible from the error bars shown in this plot that if one only requires a total magnitude of faint galaxies – instead of bulge and disk magnitudes separately – single-Sérsic (multi-band) fits seemingly provide equally good measurements, although the top panel suggests that for bright galaxies, a B/D fit actually recovers the total flux of these objects better. This is possibly an effect of a single-Sérsic profile not being able to perfectly describe the profile of these two-component objects.

However, in order to get a better insight into the systematic effects at work, we can look at the recovery of the total magnitude as a function of B/T ratio and galaxy magnitude, as shown in Fig. 9. This plot includes galaxies at m_r,sim > 19.8 (i.e. our entire dataset). Both fits, single-Sérsic (left) and B/D (right), recover the total galaxy magnitudes well at all magnitudes and B/T. While the scatter/uncertainty (bottom panels) clearly increases for both fainter and more ‘bulge-dominated’ objects in both fits, as one would naively expect, most values for systematic offsets (top panels) are well below 0.1 magnitudes.

From Fig. 9, it is somewhat unclear whether faint galaxies are indeed better fit with single-Sérsic models, as was suggested by Fig. 8. It seems that the situation is somewhat more complicated than that, and B/D fits in fact give – on average – a more accurate total magnitude (less scatter, smaller boxes in bottom panels), depending in detail on the galaxies investigated. This seems to be specifically true for galaxies fainter than m_r ~ 19.5 for which all fits produce bad results much more frequently.

At this point, a brief comparison with Mendel et al. (2014, Appendix B), who show equivalent plots for their catalogues, highlights the strength of the Galapagos-2/Galfitm approach. On Sdss data, they reported typical errors of the order of ~0.1 mag when recovering the total magnitude of an object for both single-Sérsic models and B/D fits at m_r ~ 17.5 mag. In this work we use the same Sdss data in combination with additional Ukidss data, and we derive total magnitudes with better accuracy at much fainter magnitudes, for example ~0.05 mag at around m_r ~ 21 mag, than with either of our methods, B/D or single-Sérsic fits.

For bulge and disk magnitudes, Mendel et al. (2014) see similar effects to our work in that the fainter component is systematically fit badly once it becomes too faint. They report typical disk uncertainties of >1 magnitudes for all galaxies at B/T > 0.9. However, it should be pointed out that for the other extreme (i.e. the brighter component), their method seems to do a better job than what we see in this work, with barely any systematic offset.

Unfortunately, Mendel et al. (2014) only analyse objects at m_r < 17.7 mag, resulting in only minimal overlap with this work. For the objects analysed by both works, we attempt a more direct comparison in Sect. 6.1. From the brief comparison here and the accuracies cited above, it seems that with our analysis we can indeed push B/D fits several magnitudes fainter at the same accuracy and scatter.

4.3.3 Bulge-to-total ratios

Finally, there are two more important parameters in context of measuring galaxy magnitudes other than the component magnitudes themselves: Colours and B/T ratios, both of which are widely used in the community to ‘classify’ galaxies. B/T ratios are probably the value that is most desired when running B/D decompositions of galaxies. Unfortunately, defined as a ratio, they are very sensitive to the fluxes in the two components and hence hard to measure accurately. For example, in a B/T = 0.5 galaxy, distributing the flux between the components wrongly by even 0.2 magnitudes¹⁰ changes the B/T ratio by 0.1. We illustrate this in Fig. 10, where we show the input (B/T)_r,sim and output (B/T)_r,fit ratios of the simulated galaxies. For bright objects (left panels), both multi-band (top) and single-band fitting (bottom) do a similar job at recovering B/T ratios, although the one-to- one correlation is more pronounced in the multi-band fits that also return a result more often, leading to a larger samples size, as indicated in the top left of the panels. For each bin in (B/T)_r,fit of 0.2 width, as one would derive them on real data, we calculate the median (B/T)_r,sim and the 30/70 percentiles as an indication of the width of the distribution for both multi-band fits. While both methods show large scatter, neither method shows any strong systematic offsets. However, this is different when looking at fainter galaxies (right panels, I8.5mag < m_r,B+D.sim < 19mag), where separation of two components is naturally harder. The scatter in both single-band and multi-band fits increases dramatically, measured B/T ratios become very uncertain. While multi-band mean values are not very accurate on a individual basis, their average values are less biased than the values derived from single-band fitting and actually recover the one-to-one line reasonably well. The average values of the single-band fits are closer to the centre of the plot, indicating a random distribution of B/T values. For fainter galaxies, the recovery of the B/T ratio becomes even more challenging and the values less reliable.

We look at the recovered B/T as a function of (B/T)_r,sim and m_r,B+D,sim in Fig. 11, where we plot the difference (B/T)_r,flt – (B/T)_r,sim. A clear trend is visible, in that objects with low (B/T)_r,sim get overestimated in their (B/T)_r,fit, objects with high (B/T)_r,sim get underestimated. This trend is, of course, by design, as it is impossible to overestimate the (B/T) of a galaxy with (B/T)_r,sim ~ 1. However, as we can see at intermediate (B/T)_r,sim values, the recovered B/T ratios are in general measured within 0.2 of the true values over a wide range of magnitudes and 0.2 < B/T < 0.8. Measured B/T values in more extreme galaxies are indeed less well recovered, which we try to explain below. We note that the lack of objects at m_r,B+D,sim < 18 (while Fig. 10 specifically uses objects at m_r,B+D,sim < 17.5) is artificial as we only plot bins with more than 20 objects and entirely ignore objects with m_r,B+D,sim < 17.

In Figs. 3 and 4 in V14b, we have already shown how well galaxy parameters can be measured in real galaxies that are artificially redshifted to larger redshifts (i.e. becoming smaller and fainter) and show a small shift in observed bands, which plays a minor role in this context). We have found that multi-band fits are able to recover B/T ratios more consistently over the redshift range tested. In Fig. 9 of the same paper we have further shown that the measured B/T ratios of different galaxy types are harder and harder to measure with increasing (artificial) redshift, as expected. This effect was visible in both multi-band and single-band fits, but was more pronounced and less smooth (with redshift) in single-band fits. The results presented in this work underline these findings, but provide a more statistical analysis as the galaxy samples used are significantly larger and span a wider range of parameters. Further, as simulated galaxies have been used here, we can compare input and output values directly, rather than relying on the ‘smoothness’ of trends on fitting and artificially redshifting real galaxies.

4.3.4 Colours and SEDs

The second important parameter – and the most critical for many science cases – that users might want to derive from B/D fitting are colours – or SEDs – of the individual components. When trying to measure, for example, stellar populations of bulges and disks, recovering the SEDs of the components – at least on average – becomes vital.

The result of such a test can be seen in Fig. 12, where we show the SEDs recovered by the different fitting methods for galaxy samples of different brightness. For this test, it is important to highlight once more that all bulges and all disks – apart from small variations on the simulated values – show the same SED, respectively, in the simulated data, making this test feasible. These input SEDs are shown in these plots in solid colours, green for disks, orange for bulges. The average SEDs of each component recovered by single-band fits are shown as dotted lines (blue for disks, red for bulges), the average SEDs recovered by multi-band fits are shown as dashed lines. For reasons of readability, we do not show the scatter in the distributions in all panels, but in the second panel only.

For bright galaxies, both fitting methods recover the actual component SEDs very well (top panel, m_r,B+D,sim < 17.5) on average, even in the u band. An investigation into the scatter, however, shows that multi-band fitting recovers the input values accurately for each object, while single-band fits show significantly larger scatter. Depending on the band and component, the scatter in the single-band results is 1.5–2 times larger.

However, in galaxies with 17.5 < m_r,B+D,sim < 18.5 mag (middle panel) average single-band SEDs already deviate significantly from the input SEDs, wrongly suggesting that disks and bulges have similar SEDs. In contrast, multi-band fits recover the average component SEDs much more accurately. The scatter in the SEDs recovered is significantly worse in the case of singleband fitting, as indicated with the thin error bars in the middle panel, thick error bars show the scatter in multi-band fitting. Especially in i, z, and Y band, the scatter in single-band fitting is up to three times worse than in multi-band fitting. Single-band fits make the bulge and disk SEDs statistically indistinguishable at all wavelengths, while multi-band fits statistically allow the SEDs of disks and bulges to be separated much more cleanly. The exception in all cases is the u-band data, in which the bulges of these objects in particular are so faint that their u-band magnitudes cannot be recovered well in either of the setups; the recovered fluxes basically serve as an upper limit.

At 18.5 < m_r,B+D,sim < 19.5 mag, small deviations from the input SEDs can be seen for multi-band fits; however, on average, the bulge and disk SEDs can still be recovered relatively accurately and are clearly different. Single-band SEDs on the other hand, are mostly identical, with very large scatter.

Figure 12, above all others, demonstrates the importance and the significant improvements (besides object numbers) achieved by using multi-band fitting, compared to single-band fitting, as it allows – at least on average – the analysis of component SEDs. In the framework of Galapagos-2, this allows the understanding of stellar populations of large numbers of galaxies and their components, present in typical data derived by present-day surveys and inaccessible by single-band fitting, at least in the brightness regime examined in this work. It is this capability that allows the limits of an investigation and science case to be pushed once multi-band fitting is used.

4.3.5 Galaxy sizes

In Figs. 13–16 we show how well bulge and disk sizes can be recovered by comparing sizes measured to the simulated sizes of the components. Similar to the previous section, we plot the objects with good fits in all ɡriYHK bands in both single- and multi-band fits, to allow easier comparison between bands, and the largest possible subsets of those samples in the uzJ bands, respectively. Although this is not important in most places as we show ratios of sizes, galaxy sizes are given in units of pixel¹¹ throughout this paper, unless indicated otherwise specifically.

In Fig. 13, we show a direct comparison of simulated and fitted bulge sizes (left) for multi-band fits (top) and single-band fits (bottom panel) in the r band, before looking at this parameter in a more statistical approach. This plot shows all objects with fit results in all ɡriYHK bands and multi-band and at m_r,B+D,sim < 19.5 mag. As can be seen, and has been expected, the recovery of the bulge sizes is very noisy, especially when using such a faint limit. However, a correlation between simulated and fitted values is apparent, although the scatter is large. The single-band results additionally show a more significant systematic offset at sizes at r_e,r,sim > 5pix, in that the fits on average return smaller sizes, as can be seen from the rolling median lines that are over-plotted in each panel. We confirm these findings in Fig. 14.

It should be pointed out here that – as the sizes of each component are held constant with wavelength during both the multi – band fit and in the simulations - the multi-band plot would look identical when showing the results of any other band. This fact also explains the identical average and error bars of the multi-band sizes at all wavelengths in several of the successive plots. Where differences between bands can be seen, they are a result of the sample selection in this particular band (i.e. a subset of the other bands).

In the right column of Fig. 13, we show the equivalent plot for disk sizes. Both single- and multi-band fitting show a much narrower relation than for bulges, and outliers favour larger fit sizes, although with only small effects on the median values. The relation of the multi-band results, however, is significantly tighter. No systematic deviations from a one-to-one line can be found, however.

In Fig. 14 we show the ratio of input and output values for samples of bright and faint galaxies. Values >1 indicate that the fit recovers the galaxy component too large. For bright galaxies (left column), it can be seen that both single- and multi-band fits recover the component size without systematic offset, but the scatter in the multi-band fits is significantly smaller. There is also a interesting trend in the single-band fits that in the redder bands, the distribution of bulge values is asymmetrical towards smaller sizes, while in blue bands, it is much more symmetrical, whereas disk sizes (bottom panel) generally seem to be overestimated more often than underestimated. Both effects cannot be seen in multi-band fits, as by design the sizes of both bulge and disk are respectively constant with wavelength. At first glance, it is interesting that the multi-band fits show smaller error bars in the u band than in the other bands, given that in many previous figures, the u band showed larger scatter. This, however, is by design of the multi-band fits. Firstly, other than magnitudes, the size measurements discussed here are not independent from the other bands, as they are constrained with Chebyshev polynomials to follow the trend in the other bands smoothly. As such, they are less dependent on the image depth of the u band itself. Secondly, the u band (as well as z and Y) show only a subset of the objects, presumably the brighter ones that are easier to fit. In fact, given that the error bars returned by Galfitm are, by design, identical in all bands in the case of multi-band fits, these smaller error bars in the u band is precisely an effect of the sample selection.

For faint galaxies (right column), both fits behave similarly in that they slightly underestimate the bulges sizes and overestimate the disk sizes. However, the systematic offsets are smaller in the case of multi-band fits (<5% for disks, approximately −20% for bulge sizes), and the offsets in single-band fits are significantly larger (we notice the logarithmic scale on the y axis). In faint galaxies, single-band fitting cannot reliably measure component sizes, with bulge sizes being badly underestimated and disk sizes badly overestimated, by a factor of up to a few in individual galaxies. Multi-band fits, while showing increased scatter compared to brighter objects, return much more reliable and consistent results. We again note that the different offset of the multi-band values in the u band is a result of using a different sample and reduced sample size in this band, the numbers of fits for each band are shown above the data points. By design sizes in the u band are the same as in all the other bands in the fits. As u-band single-band fits return good fits results in only a few galaxies, a significantly smaller sample can be analysed.

In the top two panels of Fig. 15 we look at the ratio of the measured and simulated sizes of the bulges and discs, respectively, using multi-band fits. The plots follow the same layout as Fig. 7, displaying the ratio (offset to unity) of the sizes (top) and the standard deviation (middle) of the measurements by the size of the squares, plotted as a function of B/T ratio and total galaxy magnitudes. These plots confirm that, in general, the bulge sizes are underestimated by ~10–30% in faint objects, disk sizes are typically overestimated by 5–15%, and that some magnitude trends are visible in that deviations are in general larger in fainter galaxies, as one would expect. The scatter in these distributions are especially high for the fainter component in ‘pure’ galaxies (i.e. bulges in galaxies with B/T < 0.2 and disks in galaxies with B/T > 0.8). In particular, the faint disks embedded in galaxies with high B/T show very large scatter, as one would expect.

As was the case for magnitudes, it is also apparent that both components are hard to fit in objects with B/T < 0.2 and B/T > 0.8, and the largest systematic trends for both components exist in those areas. In Sect. 4.3.1, we indicated that this is an effect of the fainter component fitting some of the flux of the brighter component, and indeed, we can see similar effects in these plots as well. Bulge sizes are overestimated (green, filled squares) more frequently in low B/T objects, which is an effect of galaxy bulges being typically smaller than galaxy disks, although no such restriction was technically imposed on our simulated values. This result fits the scenario in which the fainter component fits some of the flux – and the shape – of the brighter component.

In fact, when we plot the ratio between the bulge and disk halflight radii r_r,B,fit/r_r,D,fit in a similar plot (see Fig. 16), we can see that this ratio is between 0.3 and 0.5 on average over most of the B/T and m_r,sim range covered, reflecting the distributions we used in the simulated data. However, in galaxies with B/T < 0.3 this ratio changes dramatically, reaching values >0.8, at small variation, with the bulge in these galaxies fitting the disk profile instead. The small scatter, despite apparent offsets, suggests that this offset is indeed systematic. This indicates that the bulge profile in the fit accounts for some part of the disk profile, and we will see the same effect again in Sérsic indices in the next section. Interestingly, in galaxies with B/T > 0.7 the opposite can be seen, with the disk becoming very large in comparison to the bulge. We attribute this to the fixed Sérsic index of the disk, which cannot mimic the n_b = 4 shape of the bulge profile, and instead fits a background or neighbouring structure.

In the bottom panels of Fig. 15, we show the relation of the single-Sérsic derived size r_e,r,ss,fit to the input values for bulges and disks, respectively. For the areas with B/T < 0.2 and B/T > 0.8, we can see much smaller offsets for disk and bulges, respectively, than in the B/D fits. This suggests that in these galaxies, the single-Sérsic fits provide a better fit to the profile of the brighter component. In this case, one would want to fall back onto using the single-Sérsic fits values instead of the results from B/D fits and truly count these objects as pure systems. We had already seen previously that the total magnitudes of these objects are recovered well by the single-Sérsic fits. Given the low flux in the fainter component, this total magnitude will largely reflect the flux in the brighter component.

4.3.6 Sérsic indices

As the last of the three important profile parameters, we need to examine how well we can recover the bulge Sérsic index from our simulated data. In general (e.g. H13 and others), Sérsic indices are known to be the parameter hardest to recover, especially for high-n objects. This should be especially true given that the Sér- sic profile of the bulge is now overlaid with an exponential profile of the galaxy disk that will both hide the profile itself and further confuse the Galfitm fits, making n_B even harder to recover. The reader should be reminded at this stage that the disks in our simulated objects are both created and fit using a n_D == 1 profile, so we do not need to examine how well those can be recovered.

Instead, we show the direct comparison of input and output bulge Sérsic indices n_B values in Fig. 17, colour coded by simulated object r-band magnitude. It is clear that n_B is indeed the parameter hardest to fit, with both multi-band and single-band fits showing significant scatter, the correlation for bright objects being more visible for multi-band fits. Additionally, single-band fits show a significantly larger fraction at high n_B;sim values that are fitted with lower n_B,fit values, also visible in the over-plotted lines that indicate rolling median and 16/84 percentiles. From the individual point in the background, it is visible that in both multi- and single-band fits the scatter increases for fainter objects (lighter colour), as one would expect. For galaxies at m_r,B+D,sim < 17.5 mag, this correlation becomes relatively tight, with a scatter of ±1 as can be confirmed in the top panel of Fig. 18. However, the sample becomes small, which is why we include fainter galaxies in Fig. 17. As one would expect, the magnitude of the host galaxy is not the only important parameter in this context. Colour-coding the plots instead by B/T (not shown) shows the expected trends in that Sérsic index values n_B are more easily recovered in prominent bulges. However, while multi-band fits show a very tight correlation for objects with B/T > 0.8, the single-band fits recover these values much less accurately, and a trend with decreasing B/T, while visible, is quite weak.

In Fig. 18, we examine the recovery of the bulge Sérsic index n_B more statistically for bright and faint galaxies. Similar to what we have seen for magnitudes and sizes, both single-band and multi-band fits on average recover the Sérsic index well for bright objects, systematic offsets are small. However, the deviation in single-band fits is already on the level of 40% for individual objects, in multi-band fits on a somewhat smaller level. For faint galaxies, both fitting methods show considerable offset and scatter, both codes underestimate nB significantly. Interestingly, this is contrary to what we found in H13 where we reported that n is generally somewhat overestimated in single-profile galaxies. In all bands, the single-band fits recover values considerably worse than the multi-band fits with offsets of up to ~50% even when excluding the noisy u band. This systematic underestimation of n_B is likely an effect of mixing the bulge profile with an underlying n_D == 1 disk profile, which pulls the fit towards lower n_B values. The scatter of the recovered values is large in both cases, and significantly increased in the single-band fits, as should be expected at this point. The large offsets and scatter make it impossible in single-band fits to distinguish between a classical de Vaucouleurs bulge with n_B == 4 and a pseudo bulge, which generally show lower Sérsic index values of n_B ~ 1.5.

In Figs. 19 and 20 we again take a closer look at the Sérsic indices recovered by the multi-band fitting. In Fig. 19 we show how well n_B can be recovered by comparing the fit value to the simulated values. As noted above, the bulge Sérsic indices are underestimated by 10%, even in multi-band fits of bright objects. The average offset is a function of the object brightness such that n_B is underestimated more in fainter galaxies. However, in faint and disk-dominated galaxies, systematic offsets become very large. In these galaxies, the average n_B is badly underestimated, with very small scatter, indicating that this is a systematic effect.

In fact, when we look at the recovered n_B values directly in Fig. 20, we see that the average n_B recovered by Galfitm is indeed ~3–4 for a large fraction of the galaxies, 4 being the average input value. However, low n_B values are recovered when trying to fit a bulge in faint galaxies (top) or within a bright disk (left). The n_B in these galaxies converges to a value of ~1, which – together with magnitudes and sizes – again shows that the bulge profile of the B/D fits in these objects tends to fit part of the disk light. Magnitudes, sizes, and Sérsic indices all tend towards the values simulated in the disk profile of these objects.

Open red symbols in this plot indicate average value of n_B < 2. A value of 2 is often used in the literature to separate classical bulges from pseudo bulges. In our analysis here, we show that these measurements such a distinction is only reliable in bulge- dominated objects as bulge Sérsic indices are hard to recover in other galaxies.

We have a specific look at simulated pseudo-bulges with 1 < n_B,sim < 2 in Fig. 21. In order to avoid a low number of objects, we increased the bin size in this plot. Equivalently to Fig. 20 we show the average recovered Sérsic index in the top panel. On average, Sérsic indices are recovered with values of 1 < n_B,fit < 2, as this was how they were selected in the simulated values. However, as this might be an effect of Sérsic indices being drawn closer to 1 due to the disk, as seen when looking at all bulges, we try to have a closer look at the recovered values. In the middle panel, we show how well the n_B values can be recovered compared to the input values. Systematic deviations are relatively small. Only for galaxies at B/T < 0.2, larger effects can be seen, which can be easily understood as in these objects the bulge would be too faint to be fit reliably, as discussed before. In these galaxies, the fit values are indeed drawn closer to 1. Comparison with Fig. 19 reveals similar behaviour, at least qualitatively. While on average the values can be recovered, we find scatter values of σn_B,fit/n_B,sim ~ 0.3–0.4 throughout the entire parameter space (not shown). Finally, in the bottom panel, we show how well we can recover r_e,B for these pseudo-bulges, in comparison with Fig. 15. Again, we find similar behaviour for pseudo- as for classical bulges.

From this result, we conclude that the recovery of pseudobulges in these data is similar to the recovery of classical bulges. However, we would like to stress that, while on average the values can be recovered, it is dangerous on an individual galaxy to draw the conclusion whether the bulge is a classical bulge at n_B = 4 or a pseudo-bulge with n_B = 1, the uncertainties on an individual galaxy are too large to allow this conclusion (see Fig. 18). Galapagos-2 offers the option to fix n_B,flt == 4 during the fit, to allow only for classical bulges. However, we did not do so during this analysis, and this test is hence beyond the scope of this work. Nedkova (in prep), found those fits to generally behave more stable in Candels data. We would also like to stress that these results only apply to the faint and ‘distant’ objects analysed in this dataset. In nearby galaxies, even an automated approach might warrant some good results, as results could be more reliable due to higher spatial resolution and generally higher S/N. Certainly, for well resolved galaxies, more care can be taken on individual objects.

4.4 Danger of using Sérsic index to classify galaxies

At this stage, it is important to highlight the danger of using a Sérsic index cut to ‘classify’ galaxies or select specific galaxy samples, as is often done in the literature (e.g. Ravindranath et al. 2004; Shen et al. 2003, and many others) and which has been reported by other authors (e.g. Graham 2019). Figure 22 shows the simulated (B/T)_r,sim values as a function of the measured n_SS values as derived from our single-Sérsic multi-band fits. In order to present a cleaner plot, we restrict ourselves to bright galaxies with m_r,Sim < 18.5. Naively, in such a plot, one would expect a relatively tight correlation between Sérsic index n_ss and B/T ratio, as a more prominent bulge would increase the measured Sérsic index of the galaxy as a whole, as their light peak in the centre of the galaxy and their generally smaller sizes drive the n_ss to higher values. And indeed, the Sérsic index is often used in this way in the literature, using values 1.5 < n_ss < 2.5. While some correlation can be seen, the trend is very weak and not able to select ‘bulge-dominated’ galaxies from a sample, at least without further knowledge. In order to check whether this is an effect of simulating galaxy bulges at a range of Sérsic indices instead of using classical bulges only, we highlight the galaxies with classical bulges (n_B,Sim > 3.5) in red. The correlation is, in fact, very similar, and even for these galaxies, using a cut at n = 2 to separate bulge- from disk-dominated objects should be considered unreliable. Of the 1005 ‘spheroid-dominated’ galaxies with n_r,ss > 2 in the r band, which would be classified as bulge-dominated in such a selection, 350 (~35%) were in fact simulated with (B/T)_r,sim < 0.5, of the 944 galaxies with n_r,ss < 2, 175 (~18.5%) have (B/T)_r,sim > 0.5, so a somewhat cleaner sample can be selected.

This effect is visible in all bands, although it is somewhat less pronounced in the H-band and in fact the opposite for ‘disk-dominated’ objects. The H-band results still show a contamination of ~15% (of 1469 galaxies with n_H,SS > 2 225 have (B/T)_H,sim < 0.5) for ‘spheroid-dominated’ objects, and ~48% (of 480 galaxies with n_H,SS < 2 234 have (B/T)_H,sim > 0.5) for ‘disk-dominated’ objects, respectively. For completeness, it should be stated that we created similar plots for other wavelengths and mixed bands (e.g. checking for correlation of the n_H,ss with (B/T)_r), but no strong correlation can be found; in fact trends look even weaker when mixing wavelengths. From this analysis, we discourage anyone from using such a simple classification approach.

We also checked for other correlations that could potentially allow such an automatic identification of such a class of objects, but we were unable to identify a different, more reliable method.

5 Application to real imaging

In Sect. 4 we present a detailed analysis of code performance on simulated data and have discussed the advantages of such an approach. However, using such simulated light profiles assumes that real galaxies actually do precisely follow Sérsic profiles - although it is already known that they do not (e.g. Graham et al. 2003; Trujillo et al. 2004, and numerous others). Any deviation from Sérsic profiles, for example by additional components, galaxy-galaxy-interactions, dust content (especially dust lanes in edge-on galaxies; e.g. Pastrav et al. 20l3a,b), and/or simply substructure like spiral arms within the disk, will not be included the simulations used in this work, and hence their influence on the test results cannot be tested. We discuss the effect of dust specifically in Sect. 7.

In this section we look at comparisons of single- to multiband fits on real data from the Gama survey, the very data that we tried to replicate in the previous sections. The obvious advantage of this is that all effects apparent in real galaxy surveys can be tested. However, as the true parameter values of objects are unknown, a detailed analysis as in the previous section is impossible. We can merely run sanity checks and see whether already known effects can be recovered better with one method than the other. We can further compare our fit results to results from other works and with fits on alternative (i.e. deeper) data, which we briefly attempt in Sects. 6.1, 6.2, and 6.3. Again, however, these checks can only serve as sanity checks, a decision on which code is ‘better’ is not possible as true parameter values are unknown.

5.1 Data

The data used in this analysis have been well described in Driver et al. (2011), H13, and others, so we refer the reader to those publications. In order to be consistent with H13, we carry out the analysis between single- and multi-band fits on part of the Gama-09 survey region using the same data as H13, although newer and deeper data do exist. Additionally, the simulated data used above have been made to these specifications, so a comparison of results to those data is easier. The data used here have been provided by Gama and use large mosaics, made from re–ing (Bertin et al. 2002) data from Sdss and Ukidss LAS archival data (Lawrence et al. 2007) onto a common pixel-grid. Instead, deeper data now uses deeper input data from KiDS (de Jong et al. 2013; Kuijken et al. 2019) and VIKING (Edge et al. 2013). We use these deeper data in Sects. 6.2 and 6.3 and in Appendix F.

5.2 Setup

The code version and setup used in this section is largely identical to the one used on simulated data. Especially, the same versions of Galfitm and Galapagos-2 are employed. All DOF used for the individual parameters are the same, with the difference that it is – contrary to the simulated data – unknown whether these choices really do fit the real data used here.

The only difference in this setup is that we made use of a newer feature in Galapagos-2, which allows specific galaxies to be targeted, in order to save CPU time. Whereas in the simulated images it was a sensible approach to fit and analyse all objects, we restricted our analysis here to the galaxies that have redshifts provided by the Gama survey, by providing a list of right ascension (RA) and declination (Dec) coordinates to Galapagos-2. However, as Galapagos-2 works through objects in order of brightness, starting with the brightest objects, this approach potentially impacts on the fitting performance – both in speed and accuracy – if bright neighbouring galaxies are not included in the target list. Fitting bright objects several times as ‘neighbours’ costs more CPU time than fitting them once and ‘subtracting’ them from consecutive fits, as Galapagos-2 does. Instead of targeting only the Gama galaxies themselves, we hence additionally target nearby objects (within ~60′′, i.e. 200 pixels at 0.331′′ pix⁻¹) that are brighter than said Gama target. In tests on smaller datasets, it has been established that this minimises the impact of targeting specific galaxies on the fitting parameters (see H07 and H13).

5.3 Samples

The selection of good fits in real galaxies presented in this analysis is somewhat different to the one used in Sect. 4. There, we were able to use simulation values to select galaxies, for example by brightness, size, or B/T ratio, in order to get as clean a sample selection as possible. This is obviously not possible in real data, so a selection based on fit parameters has to be used instead. Additionally, the simulated data did not include stars, which in real data have to be taken care of and filtered out from any analysed galaxy sample.

The sample used in this section is identical to the sample used in H13, and we therefore refer the reader interested in the details to that work. In fact, as Galapagos-2 requires a first step of fitting single-Sérsic profiles to all galaxies, the sample discussed in H13 was directly used as the first step of the fitting process on which these subsequent B/D fits are based. As neighbouring objects are treated as single-Sérsic objects, with parameters kept fixed at their best-fit values in B/D fits in Galapagos-2, the fits of the Gama objects are independent of any B/D fit on neighbouring galaxies. This means it was sufficient to only target the Gama galaxies themselves in this step, saving further CPU time. While we updated the redshifts and masses known for these objects where newer data have since become available, this plays only a minor role in this section as fit results and apparent magnitudes and sizes are compared directly. However, given that we use a redshift selection to separate galaxies from stars, this potentially slightly changes the sample, as these classifications might have changed since. No update of the Gama object sample selection, however, has been carried out; subsequently, we use the same sample used previously.

Equivalently to the simulated data, we excluded any object – or component – where the fit parameters were on – or very close to – one of the constraint values used during the fit, using the same selection presented in Sect. 3. This removes bad fits from our analysis in the same way as in H13.

5.4 Results

The most basic possible test on B/D fits is to compare the total magnitude of objects between different fits, one- and two-component. Figure 23 shows the comparison between the total r-band magnitudes derived from B/D fits and magnitudes from single-Sérsic fits, both for multi-band and single-band fits, using the same object sample. For bright objects, both multi-band fits (top panel, all point) and single-band fits (bottom panel, all points) agree well between the two fits. However, for fainter objects it is clear by visual inspection that single-band fitting produces catastrophic outliers – in which single-Sérsic and total two-component magnitudes do not agree – far more frequently, and much more severely. While in neither case we find any significant systematic offsets (cf. Fig. 8), the number of outliers seen in the multi-band fits is significantly smaller. However, it should be noted that, given the large number of objects in this sample, this effect is small and does not have a significant impact on the median or 16/84 percentiles shown as dashed lines, nicely following the one-to-one line.

The fact that the total magnitude of a B/D fit agrees well with a single-Sérsic fit, however, does not prove that the B/D decomposition itself is sensible and physically meaningful; this could simply be an effect of the minimisation process, which tries to take all flux in the image into account. As the neighbouring galaxies are taken into account in identical fashion to the single-Sérsic fits, the total magnitudes recovered on the primary target should be expected to be relatively similar and simply reflects the ‘remaining’ flux in the image, with single-band fits showing increased scatter as they use less data. The larger amounts of extreme outliers in the single-band plot, however, is both interesting and worrying. As Galapagos-2 constrains the position of the two components in the B/D fits to be the same, and this position to be within 0.5 * re_ss, this cannot be an effect of one or both components trying to fit a different object entirely or part thereof by ‘wandering off’. Similarly, objects with extreme axis ratio, as would be the case if one component tries to include any residual flux of a neighbouring objects, should be ignored in the analysis, and in visual inspection are not found to happen very often. Any extreme such cases would be automatically classified as ‘bad fits’ due to our sample selection, which flags and excludes such extreme objects.

However, by visually checking extreme outliers, we have found that many of these objects are indeed galaxies in which one component shows an extreme size. We have highlighted a number of objects in green in which one of the galaxy components (disk or bulge) shows a size of r_e > 200 pixels (i.e. possibly unrealistically large) but is not filtered out by our sample selection, which selects at 395 pixels. This value has been empirically derived to be ‘large’ compared to typical single-Sérsic sizes, and only serves the purpose to select these extreme galaxies in this figure. Especially users of Galapagos-2 should take care of this limit and carefully select a value for their own dataset, if they wished to take this into consideration in the sample selection.

As can be seen, these galaxies are in general rare in multiband fits, but in single-band fits, these objects account for the vast majority of outliers. Visual inspection revealed that the fitting of the residuals of one or more neighbouring objects or other residual flux in the images is the main reason for these extreme component sizes, but not with such an extreme axis ratio that they would be filtered out by our sample selection above. Such a large, relatively face-on component, despite being quite surface brightness faint, can include a lot flux due to its large size, leading for the brightness of an object to be overestimated badly. The sensitivity of the single-band fits to such effects causes the large number of outliers in this plot, while the robustness of the multi-band fits to avoid fitting such structures prevents this from happening very often.

In order to evaluate how well the B/D fit truly separates the bulge and disk in a galaxy, we need to look not at the total magnitude, but their fractions at each wavelength, or – better yet – the SEDs of the individual components. This is difficult on an individual basis, especially when the true values are not known and each component (and each galaxy) shows very different stellar populations. Figure 24 shows the average SEDs and colours of bulges and disks – normalised to the r band – as derived from the fits for three different magnitude (r-band) ranges and their scatter (16 and 84 percentiles) of the values in the second panel. For this plot we use the largest possible common sample for each component, that is, all multi- and single-band fitting produces a good result. However, we exclude the u and z single-band fits from this definition, as they would decrease the sample size too drastically (i.e. we require griYJHK single-band fits to return a good result). For these two bands only, we additionally limit the sample size in that band (i.e. ugriYJHK and grizYJHK, respectively). While this selection limits the sample size in this plot and in each band, it makes the average SEDs directly comparable, as the same sample is being used at most wavelengths. In all panels, more disks than bulges are used, as disks can generally be recovered better, as we have seen in Sect. 4.3, although the difference in numbers is somewhat more extreme here.

It becomes clear from this plot that, at least statistically, single-band fitting cannot separate the bulge and the disk in galaxies (blue and red dashed lines, respectively), as both components on average are recovered with the same SEDs, even in the brightest galaxies. As such, also any B/T ratios derived using single-band fitting in each band would be highly inaccurate, as we have also shown in Fig. 10. In the middle panel, we show the width of the recovered colour distributions as error bars, thin lines indicating single-band results. It becomes clear that the colours recovered not only are very similar on average, but also their distributions are nearly identical, indicating that the singleband fits are unable to separate the flux of the two components. In contrast, the multi-band fits (blue and red solid lines) successfully recover very different SEDs for both components, and there seems little effect of the galaxy brightness on these SEDs. The error bars, thick lines showing multi-band results, also show much less significant overlap in the distributions. The exceptions here are the u and z bands in which the reduced sample shows significantly wider spread, especially in the galaxy disks. We also notice a dip in the J-band disk SED (blue solid line), which also seems present in many of the objects at all brightnesses. We are unsure what causes this dip, but as it is not present in the average bulge SED (red solid line), we conclude that this is not a data (e.g. zeropoint) issue.

For comparison, we also over-plot the SEDs (using values from NED¹²) of a typical elliptical galaxy (M 89, in the Virgo cluster, orange dashed line) and a known galaxy without a prominent bulge (NGC 0337, classification SB(s)d, green dashed line). Especially in case of the NGC 0337, multi-band fitting can recover this typical SED quite well as the average SED of galaxy disks. The average bulge SED also agrees well with the SED of a quiescent galaxies, but we note the offset in the redder bands. As M 89 is a relatively massive galaxy and in a cluster, this difference might be at least partially attributed to differences in metallicity or in the age of the stellar populations. Indeed, the difference is found to be smallest for the brightest sample, where the most massive galaxies are more likely to be.

As a further sanity check, we took a look at the component r-band sizes in Fig. 25 in comparison to the single-Sérsic sizes for all galaxies with m_r < 18.5 (samples using different magnitude cuts look similar). The samples shown for disks and bulges, respectively, are the same in both single- and multi-band fits, that is, both fits produce a good fit result for the respective components, so the histograms for disks (or bulges) should be directly comparable. Again, one can easily see from the scale on the y axis that Galfitm produces good fits for disks more often, the disk sample being significantly larger. It is well known in the literature that bulges are in general smaller than the disks they reside in Lange et al. (2016, and many others) and one would expect this to be visible in such a plot. And indeed, in both cases – single- and multi-band fits – the bulges of these galaxies are much smaller than the measured single-Sérsic sizes, and the disks are significantly larger on average. For these real objects, it is impossible to tell which of these methods produces the better fit, but it ties in well with our finding in Sect. 4.3 that the two methods produce different results, although some trends seem to be reversed here. For the simulated galaxies, single-band fits were found to provide larger sizes for disks and smaller sizes for bulges (see Fig. 14). In the real galaxies presented here, the behaviour is in fact the opposite with single-band fits producing larger bulges than multi-band fits. It is somewhat unclear where this difference originates. However, we can speculate that it might be connected to the fact that in the simulated data both disk and bulge were simulated (and fit) with a constant size across all wavelength, while this might not true in real galaxies (e.g. Xilouris et al. 1999). However, real galaxies are still being fit with a constant size across all wavelengths. These fits can hence be influenced by the colour gradients present within the individual components and by the effect of dust, neither of which are present in the simulated data.

At this point, we would like to remind the reader that the starting values for bulge and disk sizes in these fits are defined as r_e,D = 1.2 * r_e,ss and r_e,B = 0.3 * r_e,ss (see Sect. 2.2). Hence, it seems suspicious that especially the bulge sizes on average are found to somewhat centre around this value. While in principle this can be induced by the starting values, we verified that it is not the case by re-fitting a large number of objects by using different starting values, e.g r_e,B = 1.0 * r_e,SS, r_e,B = 0.5 * r_e,SS and r_e,D = 1.0 * r_e,SS. The resulting histograms show very similar results, at least for bright galaxies, and we briefly discuss the different starting values in Appendix D. The values found in these tests were in fact one of the reasons why we chose the starting values the way we did.

6 Comparison to other work

In the following sections we present several comparisons to fits carried out with other light profile fitting codes and/or data, and discuss the shortcomings and uncertainties of some of these codes, including our own. These sections can serve as a warning against taking the results from any of these codes and methods as the unbiased truth, especially in certain parameters.

6.1 Comparison to SDSS from Simard, 2011

The most prominent publicly available catalogue of a large sample bulges and disks has been presented by S11. In this work, they use Gim2d to fit 1.2 million galaxies in Sdss data at m_ptro,r,corr ≤ 17.77. By comparison with Fig. 23, it is clear that these galaxies compare to the brighter end of our galaxy sample, so supposedly the ‘easiest’ galaxies to fit. In the previous sections of this paper we showed that we generally reach objects ~2 magnitude fainter than this limit. Since the optical bands of the Gama data that we use throughout this work were created by resampling the same Sdss data used in S11 (i.e. the data used by us and by S11 are nearly identical), these catalogues are well set up for direct comparison. However, due to the magnitude limit imposed in S11, and the much smaller survey area of Gama compared to Sdss, a direct overlap of the sample is small, so we analyse here ~5000 galaxies present in both datasets.

Some caution has to be taken with the results in this section, however, as these galaxies are also the largest and most detail-rich galaxies in our sample. Experience with Galapagos-2 shows that these galaxies are somewhat harder to fit using a fully automated method. Generally speaking, smoother profiles of more distant galaxies are somewhat easier to fit, as no spiral arms or other internal features influence the fit as heavily as in some of these cases. Nearby objects are also more likely to be split up into several detections when using a fully automated approach, which we try to minimise by our ‘hot’ and ‘cold’ setup approach for SExtractor (please see H13 for details).

In Fig. 26, we compare the main parameters of the profile fits in the r band, as this is the deepest image in the dataset and has been used as the ‘main’ band in both these works. Each panel shows our values from multi-band fitting along the x axis, and the values from S11 on the y axis, each plot is labelled with the parameter it shows. In this figure, we use the largest possible sample for bulges and disks, respectively (3038 bulges, 3530 disks). As we have seen before, disks are generally fit well more often, which is why the sample in the panels showing disk parameters is somewhat larger, although in these bright objects, this is less of a problem than in Sect. 5. Where both parameters are required (e.g. for B/T), the sample is restricted to those objects for which both components returns good results. The red points in each panel indicate the additional objects with B/T > 0.8 (109 objects) in case of bulge parameters and B/T < 0.2 (1385 objects) in case of disk parameters. For these objects, we have shown in Sect. 4.3 that the parameters of both components are not reliably recovered, as the fainter component tends to fit part of the brighter component. Instead, we plot the results from the single-Sérsic fit, which we have shown reflects the parameters of the brighter component better in the respective panel.

The top two panels show the agreement between the total galaxy magnitude on the left (parameter RG2D for Gim2d, m_D + m_B for Galfitm), and the single-Sérsic half-light radius (parameter RHLR for Gim2d, converted to pixel values) on the right. To guide the eye, the green dashed line marks the one-to-one line, indicating perfect agreement. For these global parameters, both codes agree very well, although some outliers (e.g. through mis-identified objects, or galaxies split up during the detection process) are visible.

The next two panels compare the bulge magnitudes on the left, and the disk magnitude on the right. The Gim2d values in these panels are derived by using the total magnitude RG2D and the bulge-to-total ratio __B_T_R as given by S11, while the values for Galfitm are a direct output of the code. In general, the agreement between these values is good, but the scatter, even on these bright objects, is substantial. Differences of 1 magnitude between the two codes are common, although only a systematic deviation for the average bulge magnitudes of fainter objects is visible. There also is a small offset in that Gim2d systematically fits disks a little fainter, but this is a small effect. The red points, indicating the single-Sérsic results of pure bulges and disks align well with the findings of the other objects, indicating that our decision to use single-Sérsic results in those cases helps to recover the true component parameters.

The central panels show the size of bulges (left, parameter RE in S11) and disks (right, parameter RD), respectively. As these values are given in kpc in the S11 catalogue, we transferred them back into arcseconds using the redshifts provided by S11, and further into Gama pixels, as this is the parameter we worked with throughout most of this paper. We further take into account a factor of r_e,D = 1.67 * scalelength for disks, as S11 provides those values instead of half-light radii. Disk sizes are generally very consistent between our fits and the ones presented in S11, although some scatter is visible for large galaxies where Galapagos-2 seems to measure somewhat larger sizes. This consistency is not true for the bulge sizes, the scatter between the measurements is very large, and the correlation between the sizes measured is quite weak. Interestingly, the pure bulges in our sample are consistently measured larger in our analysis than in the S11 measurements. We are unsure why the agreement between the codes is so weak, but it is worth noting that after analysing the same plot for bright bulges only (no plot shown), this effect is not dominated by faint bulges, but seems to be apparent over all bulge magnitudes. There seems to be no trend with bulge brightness beyond the general effect that fainter bulges seem to be fit smaller by both codes, which naively makes sense in the context of fainter galaxies likely being further away. One could speculate here that this is again an effect of object detection and might be dominated by split up galaxies, but we cannot show this from the work done here.

A similarly weak correlation between the parameters recovered by the two different codes can be seen in the panel showing the bulge Sérsic index n_B (parameter NB in S11), fourth row, left panel. As is increasingly common in the literature, we show the Sérsic index on a logarithmic axis. While both codes generally seem to recover relatively high values of n_B, the values from Gim2d seem to be somewhat higher than one would expect for classical bulges (green dashed lines). Indeed, we see in Fig. 18 that faint bulges seem to be underestimated by Galapagos-2 in their Sérsic index to some extent. Upon inspection (no plot shown), we can indeed confirm that this effect can be seen somewhat in this plot, galaxies with low n_B values being more likely to be faint objects. We note that S11 seems to impose a lower limit of n_B ≥ 0.5 during their fits, while our values are limited to n_B ≥ 0.2. As these small values are only found in a few objects, either choice seems somewhat arbitrary. Both codes use an upper limit of n_B ≤ 8.

As both we and S11 fit galaxy disks with n_D = 1, it makes no sense to show a comparison of this value. Instead, in the right panel, we show the comparison of the r-band B/T ratio as derived by the codes, where S11 provides this value directly (parameter __B_T_R). Whereas both bulge and disk magnitudes showed clear correlations, this correlation is washed out in this panel, as B/T is defined as a fraction of parameters, and hence more sensitive to small changes. Although Gim2d systematically returns higher B/T values (worse to at higher B/T galaxies), a very noisy correlation can clearly be seen. This difficulty to recover the same B/T values is of course connected to the discrepancies in the bulge sizes and n_b values that are recovered, as these three parameters are known to correlate in Sérsic fits. Upon inspection (no plot shown), this correlation is somewhat tighter for brighter galaxies, as one would expect. Especially catastrophic outliers appear to happen less regularly in objects at m_r ≤ 17, in agreement with our findings shown in Fig. 10. The objects with B/T < 0.2 (left) and B/T > 0.8 (right), which we call ‘pure disks’ and ‘pure bulges’, can be seen in this plot directly. In fact, this is how they are selected. For the pure disks, the recovered B/T values by S11 are indeed predominantly low, while most pure bulges are recovered at high B/T values by S11.

In the last row, we show the axis ratio q_B (parameter E in S11 of the bulge on the left, and the relation between the axis ratio q_D for Galfitm and the disk inclination (parameter I in S11) for Gim2d on the right. In the case of the disk axis ratio, we assumed an intrinsic thickness of the galaxy disk¹³ of 0.18, following Pizagno et al. (2005), to transfer the axis ratio measured by Galfitm into an expected inclination angle to define a one-to-one comparison line (light blue dashed line). As one can see, both parameters agree well between the two codes, despite individual outliers.

In conclusion, Galfitm (as used by Galapagos-2) and Gim2d agree well for the disk parameters in the common sample analysed here. However, the bulge parameters, especially n_B and r_e,B differ significantly between the two codes. As has been mentioned above, this does not necessarily mean that one code performs better than the other. However, in Sect. 4.3 we analysed the accuracy of the Galfitm fits themselves, and find better agreement with the input values than the one shown here.

6.2 Comparison to KiDS/VIKING fits

So far, we have used the same data – based on Sdss and Ukidss data – used in H13 throughout this paper to enable direct comparison. This was also useful for our comparison to S11 in Sect. 6.1, as the data they used was very similar. However, as has been mentioned in Sect. 5.1, newer data have since become available in Gama, which is substantially deeper, as it uses KiDS and VIKING data. As such, these data allow a comparison on how the code performs on these different – and, crucially, entirely independent – datasets.

In Fig. B.1 we show this comparison for the galaxies fit in both datasets. The two leftmost columns of plots show the entire common sample, the right two columns restrict the sample to m_r,SEx ≤ 18, using the SExtractor measurements in the deeper KiDS/VIKING data as independent reference for comparison and discussion. While in Sect. 4.3 we show that such a simple magnitude cut-off does not fully ensure a clean sample – as the fit quality is more likely a function of a mixture of morphological parameters, including B/T – it serves here as a simple cut to restrict the analysis to the objects more likely to produce accurate fit results.

The panels shown in this plot are analogous to the ones shown in Fig. 26, with the difference that we can compare the axis ratios of disks directly, rather than comparing an axis ratio to an inclination angle, assuming some intrinsic disk thickness, as had to be done before. We note that in this plot, we show the results from the deeper KiDS/VIKING data on the x axis, and the results from the Sdss/Ukidss data on the y axis. In both datasets, the selection of good sources was carried out using the same definitions as given above. An additional selection of extended objects, using the CLASS_STAR < 0.8 values provided by SExtractor on the shallow data, has also been applied to remove likely stars within the field. We comment on the effect of this additional selection below.

In the first panels of Fig. B.1, we again compare the single-Sérsic fits. Both fits recover very similar magnitudes, visible as the good agreement between the measurements, with a small offset of 0.07 mag over all objects with brighter objects recovered on KiDS/VIKING data. With ~6000 objects shown in this panel, there are only a few outliers visible. However, for galaxy sizes, it becomes clear that the shallower Sdss/Ukidss data fits somewhat larger sizes to most objects, on a level of 10–15% (median: 12.7%). This is even true for the ~1700 brightest galaxies, shown in the two rightmost panels, where a similar offset is visible. The difference between the recovered magnitudes is reduced to 0.05 mag in this bright sample. From our data, we cannot trace the exact reason for these offsets. However, it seems likely that the reason for the size offset is connected to one of two issues. First, the native pixel scales of the images that went into the -ed images used in the two datasets are different. Sdss has a native pixel scale of 0.396′′/pix and Ukidss has a pixel scale of 0.4′′/pix (using WFCAM), while the data used for the newer/deeper images shows a native pixel scale of 0.21′′/pix for KiDS (OmegaCAM) and 0.339′′/pix for VIKING (VIRCAM). These smaller initial pixels should in principle allow smaller objects to be resolved, even if the data are resampled to a common 0.339′′/pix scale in all cases. Second, the typical seeing reported in these surveys is different, with a reported median PSF full-width half maximum (FWHM) of ~1.4′′ in the r band for Sdss (Stoughton et al. 2002) and <1.2′′ for Ukidss (Lawrence et al. 2007), compared to <0.7′′ for KiDS (e.g. Venemans et al. 2015) and <1.0′′ for VIKING. While Galfitm in principle takes the PSF into account and should be able to measure galaxy sizes significantly smaller than the PSF FWHM if the PSF is well known (see e.g. H07 and H13), it is not unfeasible to assume that this works better in data with a smaller PSF. This could especially be true if the data additionally show a smaller pixel scale (leading to better spatial sampling of the PSF) and if the data are significantly deeper, as is the case here. Furthermore, the quality of the PSF itself will be different between the two datasets even though the same software, PSFEx (Bertin 2011), was used, albeit with slightly different setups. Unfortunately, neither of these are effects that we can easily identify in the simulated data, as the input PSF was by design identical to the one used by Galfitm. Both effects should become visible more strongly for small objects, which indeed seems to be the case, as we discuss below.

Bulge and disk magnitudes (second row of panels) seem to be recovered equally well as in the comparison between Gim2d and Galfitm presented in Fig. 26, especially noting the different magnitude ranges shown. Galaxies in this deep-shallow comparison are up to 1–1.5 magnitudes fainter than the ones examined in the previous section. Overall, no systematic offsets but significant scatter can be seen. Visible trends are very similar when looking at the bright objects only. These plots also confirm our decision to call objects with B/T > 0.8 ‘pure bulges’ (305 objects) and objects with B/T < 0.2 ‘pure disks’ (3402 objects) and use the single-Sérsic results instead, as shown in red. The distribution of these points generally seems to follow the overall distribution of other objects – 1644 bulges and 2569 disks, respectively – but seems to be asymmetrical in the magnitude comparisons. This offset to fainter magnitudes in the shallower Sdss/Ukidss fits of the pure objects is a result of the way we define these pure objects in the deep data, but base the decision of which parameter to plot (B/D or single-Sérsic) on each dataset individually, in order to be consistent with a scientific analysis if only one dataset existed. It is hence possible that a ‘deep’ single-Sérsic parameter is plotted against a ‘shallow’ B/D component parameter, which would, by design, favour fainter magnitudes for the Sdss/Ukidss data. Instead of over-plotting those objects for which both the shallow and the deep fits measure B/T > 0.8 or B/T < 0.2 (i.e. both fits show single-Sérsic parameters), we refer to the top-left panel, where we can see the much better agreement between the single-Sérsic magnitudes from both fits.

Bulge sizes agree relatively well within the scatter between the two datasets, with the exception of very small galaxies at r_e,B < 3[pix], which are fit larger in the shallower Sdss/Ukidss data. As one would expect, small bulges, embedded in a galaxy disk are the most challenging components to fit, and uncertainties are expected to be large. Unfortunately, in our simulated dataset very few bulges showed r_e,B,sim < 2[pix] as those simulations were created to mimic the results from the shallower Sdss/Ukidss data, where those small objects are missing, so we are unable to conclude whether the fits from the Sdss/Ukidss or the KIDS/VIKING data are more accurate in this size range. Figure 13 revealed no such systematic offset between input and output sizes in the few galaxies present, which seems to suggest that the fits on Sdss/Ukidss data indeed did behave as expected on simulated galaxies. However, uncertainties in the PSF shape present in real data, but absent in the simulated dataset, further prevents us from determining which of the datasets and fits reveal the true value. The arguments discussed on image resolution and additional image depth would/should suggest that fits on KIDS/VIKING data are better.

No such size-dependent effect can be seen in the panel showing disk sizes r_e,D. Apart from the general offset discussed above, no systematic trends are seen and the measured disk sizes generally agree well. This might be an effect of no small disks being present in either dataset, so the issue at small radii is simply masked out. Indeed, there is a hint of a similar flattening of the relation just before the cutoff at small values. It is, however, interesting that all pure disks (as defined by the deeper KiDS/VIKING data) seem to recover larger sizes in Sdss/Ukidss data.

The next panel, as in Fig. 26, shows the agreement between the codes in the bulge Sérsic index n_B. There seems to be a trend in that the fits on Sdss/Ukidss data recover somewhat smaller n_B values especially at high n_B values, something that we did not see in Fig. 17. Again, this trend is still visible when only taking the bright objects into account. However, profiles with high n_B values rely particularly strongly on the innermost pixels for a good fit. The PSF and resolution arguments made above hence would also apply here and could explain the trend seen.

The right side panel once again compares the B/T values recovered. As before, this plot shows only weak agreement as ratios of noisy magnitude measurements are inherently sensitive to small changes. The red points indicate galaxies with B/T_deep < 0.2 and B/T_deep > 0.8 as measured on the deep data. It is interesting to see that, while low B/T values seem to agree relatively well between the two fits, high B/T values are not very consistent. On the contrary, objects with B/T_deep > 0.8 show a pretty uniform distribution in B/T_shallow, although number statistics are small. This effect is smaller when only bright objects are considered, but they still do not agree perfectly. This is possibly an indirect effect that the different data quality and depth has on the fitting parameters.

Finally, in the bottom two panels, we show a comparison of the axis ratios of the recovered galaxy components. There is a effect that q values are higher – objects are rounder – when derived on Sdss/Ukidss data. This effect is much more pronounced in bulges, disks seem to agree relatively well apart from a small and constant (~0.05) offset, which decreases towards larger values of q. Once again, this behaviour could be explained as an effect of data resolution and/or PSF uncertainties. Especially in small objects, the axis ratio measurements become very sensitive to data resolution and the accurate knowledge of the PSF becomes critical. Especially if the PSF model used is somewhat too elongated for any reason, measuring an elongated underlying objects becomes very challenging and the code tries to overcompensate for the effects in the PSF.

Most of the effects present in faint galaxies are also visible in the bright objects, presented in the right two columns, which is why we do not discuss those plots in detail. However, we can use them to show that all these features are less pronounced in bright objects, as one would expect. Especially B/T values agree much better and a clearer correlation becomes visible.

However, there are a few features visible in these plots that we have ignored so far. We have already discussed the offset in r_e,B at small sizes, but there are further prominent features at n_B,deep ~ 8 and q_B,deep ~ 0.1, both visible as vertical features. These features are more pronounced for faint objects, and are less visible in the right panels that show bright objects only. However, these clusters of objects are not exclusively connected (i.e. they do not comprise the same objects). While many of the objects at low q_B,deep values indeed also show large n_B,deep values, and while ‘many’ of these objects indeed show small r_e,B sizes, none of these features vanishes by excluding objects in the other. However, the features seem to be connected to stars or stellar like sources, which is why we introduced the cut at CLASS_STAR < 0.8 in this analysis. When further removing faint objects, these features become even weaker. This CLASS_STAR cut removes the vast majority of these objects, but such a hard cutoff value does not work perfectly. We tried cutoffs at different values in the two different datasets and have found that we can achieve a similar effect by using a lower cutoff CLASS_STAR < 0.7 in the KIDS/VIKING dataset.

6.3 Comparison to ProFit

Finally, we can compare our fits with fits derived using PROFIT (Robotham et al. 2017), a new profile fitting code that uses Markov chain Monte Carlo minimisation to derive the best fit profile, potentially outperforming the Levenberg-Marquardt approach used by Galfit and Galfitm.

PROFIT results galaxies on KIDS/VIKING data have been presented by Casura et al. (2019, and in prep.) for the lowest red-shift ~13 000 galaxies in the Gama equatorial survey regions, out of which ~3100 unique galaxies fall into the correct GAMA-9h survey region analysed in this work. Importantly, PROFIT is not run on the -ed mosaics that we use for Galapagos-2, but the original KiDS images (single-band r band only¹⁴), which provide a significantly better spatial resolution than our images. Casura et al. (2019) provide the parameters for three models: single-Sérsic, Sérsic plus exponential, and a point source plus exponential. In this section we only analyse the first two of these models, where appropriate, and select the single Sérsic parameters for all galaxies that were flagged as ‘good’ fits by Casura et al. (2019), and the bulge and disk parameters only for those objects where the flag provided by Casura et al. (2019) indicated that the two-component model is indeed the best fit model to the galaxy.

This selection leaves a sample of 1903 objects for a single-Sérsic comparison (top two rows, right columns in Fig. B.2) and, after selecting only the objects with good fits from Galfitm, 427 two-component objects with good disks (and 0.2 < B/T < 0.8) and an additional 168 pure disk objects, as well as 227 good bulges (and 0.2 < B/T < 0.8) and an additional nine pure bulge objects (remaining panels), for which we plot the single Sérsic parameters in red. Numbers for the overlapping sample of the Sdss/Ukidss fits can be found in the respective panels.

There are a few important differences between how ProFit is run compared to our own fits. Most importantly, contrary to Galfitm, Casura et al. (2019) do not provide values integrated to infinite radii (although PROFIT itself does). Instead, they provide a ‘segment radius’ RAD_SEG within which the model is valid, missing out additional flux at larger radii, and which they suggest is more reliable and should be used. This choice was made to produce better fit results in the postage stamp region and less sensitivity to potential problems in the outskirts of the galaxies, but makes a direct comparison somewhat challenging. In order to allow a fairer comparison, it was necessary to convert the Galfitm values to match this definition for all profiles independently, single Sérsic, bulges and disks. In order to do this, we integrated our derived profiles out to this radius (assuming a semi-major axis, though this is not strictly true as PROFIT allows segments to have any arbitrary shape), correcting the magnitude values, and then we derived a new half-light radius that contains half of this flux. This correction indeed improves the comparison of Galfitm and ProFit values drastically. As Sérsic indices are untouched by this correction, comparing them could be challenging, as often the Sérsic index of a profile is dominated by not only the innermost pixels, but the outskirts of a profile as well; profiles with high n values contain more light in their outskirts.

Figure B.2 shows the comparison of the fit parameters derived by ProFit with our Sdss/Ukidss fits on the left, and our KIDS/VIKING fits on the right. As Casura et al. (2019) use deeper KIDS data, one would expect the comparison to our KIDS/VIKING values to be more accurate. While this seems true in some parameters, for example m_r,B, which shows a smaller number of outliers, it does not seem to be true in all parameters. Particularly, r_e,r,B and q_r,B agree significantly better with the values from the Sdss/Ukidss fits. Specifically, we see the same effects on both parameters that we saw in Sect. 6.2 in that the KIDS/VIKING fits produce smaller and rounder objects than either ProFit or the Sdss/Ukidss fits, while the agreement between those two codes seems to be better. Given the assumed data quality, this seems surprising, and it will require additional investigation to determine why this is the case; however, PSF effects are one of the main suspects for this behaviour. The agreement in the other parameters is generally good, especially for the disk parameters. Also, at least for these bright objects, there is a clear, albeit noisy, correlation in the B/T values, which we deem encouraging.

6.4 Discussion

As we have seen in this section, B/D decomposition of real galaxies is a challenging task, and different methods and codes are prone to give different answers. Even the same code on different data – for example, Galapagos-2 on Sdss/Ukidss and KIDS/VIKING data – can seemingly return different results for certain galaxies. It might hence be worth commenting on some of the issues at hand when fitting galaxies and their components.

It is somewhat speculative what exactly leads to the differences seen, because we cannot prove this with the data at hand in this work. However, apart from the obvious differences in the imaging data used and the advantages of the KiDS/VIKING data regarding image resolution and depth discussed above, the only major difference between our runs on SDSS/UKIDSS and KIDS/VIKING data, was the handling of the PSFs used for the fits. To our understanding, many of the effects we see in the comparison between these fits – and in fact, with the fits carried out using PROFIT– can be explained by PSF effects. In this section we briefly discuss the effects as we might expect them, and compare them to the figures shown above. While investigations into these PSF effects have not revealed any issues with the PSFs in any of the datasets, we discuss here the theoretical implications, as an attempt to explain the differences we see between the different fits.

In the comparison with other codes and data, three different methods for deriving PSFs for the fits were used. We do not include the Gim2d fits in this discussion as the details of deriving the PSFs are unknown to us, although they are briefly discussed in their paper.

First, in our Sdss/Ukidss fits we used PSFs provided by Kelvin et al. (2012). These PSFs were derived on-the-fly, object by object, on a subsection of the images centred on each galaxy, using PSFEX, in a mode that derives a PSF for the central galaxy/position only. PSFEX uses basis vectors to model the images of stars, and then interpolates their parameters to an intermediate position in order to create an artificial PSF image at that position. While Sigma ensures that at least ten stars are available for the task, using an area of up to 1501 × 1501 pixels in size, this method uses the limited number of stars available in this subsection of the image. On average, ~24 stars have been used to create a PSF of 25 × 25 pixels in size (see Kelvin et al. 2012 for details). While it was ensured that no major effects were visible in the resulting PSFs, the limited number of stars used can create some issues, as we will discuss. The advantage of this method is that the PSF in principle was created for the exact position of the galaxy, especially as in the process it was ensured that only stars from the same SDSS and UKIDSS images were being used, avoiding possible issues with galaxies close to regions of overlapping input images.

Second, the PSFs used by ProFit were created following a similar idea, using PROFOUND (Robotham et al. 2Õ18) to identify up to eight good stars around the object of interest. These stars were then fit using Moffat functions, and the median value of each parameter was used to create a (perfectly centred) PSF model, which is then used for the fit. As ProFit fits were run on the KiDS images directly, no issues with overlapping tiles are expected in this dataset. However, the limited number of stars used by this method can again affect the resulting PSF in the same manner discussed above. Bad centring of individual stars could have a significant effect on the PSF shape, for example when a star falls directly between two pixels. In such a case, the resulting fit could be too elongated and too large along its major axis, impacting on the resulting PSF. While quality checks on the PSFs used have been carried out, it is unlikely that they work perfectly in all situations. Specific investigations by Casura et al. (2019) looking for such effects found no indication of elongated PSFs.

Third, for our fits on KiDS/VIKING data, we instead followed the advice of the PSFEx author¹⁵ and ran the code on much larger tiles, ~8500 × 8500 pixels, creating a grid of 10 × 10 or 20 × 20 PSFs (depending on the band, i.e. fewer PSFs in the noisy u band), in order to use the full strength of the PSFEx interpolation over large numbers of stars, using third- and fourth-order polynomials. The advantage of this approach is that it led to significantly higher S/N PSFs as they were typically derived from ~2200–2500 stars in all cases. This method ensured that there was a PSF model within ~300 pixels (~100′′) of the galaxy, but the PSF has not been created for the precise position of the galaxy. While this would have been possible, it would have been unfeasibly CPU intensive to run this method for each object. For the fit, we simply chose the closest 75 × 75 pixel PSF to the galaxy of interest.

As the ‘tiles’ used in Galapagos-2 are not identical to the native KIDS or VIKING tiles, but were arbitrarily sized cuts of the -ed images provided by Gama, edge effects would be expected in those areas where multiple KIDS or VIKING tiles overlap. For example, in an area where several KIDS tiles overlap, a sudden jump of the PSF parameters could be expected, which would be difficult to model using the smooth functions used by PSFEx. This could, for example, lead to PSFs that are systematically too round and uniform, as they smooth over areas that show elongated or drop-shaped PSFs (i.e. the corners of KIDS images), making them non-ideal for the galaxies in those areas of the images. However, due to the number of stars used, centring effects or bad fits in individual stars should play much less of a role. Visual inspection confirmed that the residuals between the derived PSFs and the individual stars used are minimal in the large majority of cases, residual level are typically less than 1.5% in the central 20 × 20 pixels of the PSFs (2.5% in the u band). In the images provided by PSFEx to check the PSF quality, no parameter jumps can be detected in either the PSF parameters themselves or – crucially – in the strengths of the residuals along the overlap areas of images, indicating that these effects, if present, are not very strong.

The potential problems discussed in these cases can lead to different effects on the fitting parameters. Unfortunately, as has been mentioned before, we cannot test our hypothesis directly by using our fits on the simulated data, as by design the PSF used to create the images and the PSF used for the fit were identical. A full analysis of these PSF effects would require very sophisticated simulated images that mimic the overlap of different tiles in the -ing process and is well beyond the scope of this work. However, we can find indications in the plots where we compare the fitting results between different codes.

If the KIDS/VIKING PSFs are generally too round and uniform compared to the true PSF shape present in a galaxy image, one would expect the measured galaxy axis ratios to be systematically underestimated, as some of the ‘roundness’ of a galaxy is already presented in the PSF. This effect should be pronounced in small objects, as the PSF plays a more important role when fitting those. Indeed, such an effect is visible when comparing the KIDS/VIKING fits to Sdss/Ukidss fits in B.1 and ProFit fits in Fig. B.2, especially in the bulge axis ratios. Fits on KiDS/VIKING data show systematically lower axis ratios. Bulges are generally smaller, so such an effect would be more pronounced. However, the effect seen in the – generally larger – disk galaxies, is much weaker. While it is possible that this effect plays some role in our fits, no such effect is seen when we compare our fits to the Gim2d fits by S11, where the axis ratios of even the bulges show no systematic offset and generally agree well.

On the other hand, if the PSF shows some elongation – especially at its centre – due to centring issues or low S/N from a small number of stars going into its creation, one could expect the opposite effect. Fits would be overestimating the axis ratios of small objects, as the PSF already takes care of some of the elongation of a galaxy, at least when the elongations somewhat align. This, again, matches what we see in the bulge axis ratio q_B in Figs. B.1 and B.2. In order to make up for the residual flux perpendicular to the elongation of the PSF, the resulting size of the object would be overestimated to minimise the residuals. Indeed, this is what we see in the measured bulge sizes in both our Sdss/Ukidss and the ProFit fits and additionally in the single-Sérsic sizes in the PROFIT fits. Interestingly, the bulge sizes measured by both Galapagos-2 on Sdss/Ukidss data and by ProFit on KiDS/VIKING seem to show a hard lower limit. Nearly no objects are fit at sizes r_e < 2][pix] or ~0.6′′. There is an argument to be made that PSF and galaxy elongation do not usually align, so this effect should balance out for co-and misaligned (specifically at 90 degrees) galaxies. However, whether this is in practice true, and whether this is a linear effect that could in fact balance out on average, cannot be solved here without significant follow-up studies. Casura et al. (2019) also exclude very small objects from their catalogue during the model selection of two-component objects. While this could explain the lack of objects with very small r_e in PROFIT fits, the remaining objects suggesting a systematic upturn in the comparison with Galapagos-2 fits on KiDS/VIKING data, it cannot explain the differences between those fits and those on Sdss/Ukidss, where such a cleaning of the sample is done in a consistent manner.

We cannot from our data conclusively decide which of the datasets provides the better fit and in this paper we concentrated on comparing multi-band to single-band techniques for B/D decompositions. However, this section highlights that deriving a good PSF model for the fit is paramount to a good fit, no matter which fitting technique is used, especially for small objects, both in single-Sérsic fits and B/D decompositions. If the PSF is perfectly known, in principle all techniques should be able to separate bulges and disks well, as long as the data have sufficiently high S/N.

Besides the PSF effects, there are other issues that should be taken into account when fitting galaxies. For example, we have seen even in the simulated data, which removes some sources of errors, that the Sérsic indices of the bulges are very hard to measure. Multi-band fitting helped deriving this parameter, but even with this technique and on relatively bright galaxies (see Fig. 17), measuring n_B is challenging. Arguably, the uncertainty is too large to separate classical from pseudo-bulges, at least in the kinds of data used in this work. While this might be possible in well-resolved, nearby galaxies, where our multi-band approach using Galfitm would certainly improve the fit quality, this does not seem possible in large surveys of distant galaxies. This could indeed suggest that a fit using a classical n_B == 4 bulge could be considered superior as it removes one large source of error and minimises the confusion between bulge and disk profile when both try to fit the same galaxy component. However, we analysed the fits with Galapagos-2 of a sub-sample of galaxies with fixed n_B values on the Sdss/Ukidss and KiDS/VIKING data. While we do not show plots for this test, we can report that we did not see a significant improvement in the comparison of these fits. In fact, both r_e,b and q_B seem to agree less well between the fits, while the other parameters do not seem to change significantly.

Similarly, it is clear from this work that measurements of B/T flux ratios are very unreliable on a individual basis, at least for faint galaxies. If a science case requires such measurements, care has to be taken in the selection of a bright, reliable and/or sufficiently large galaxy sample such that individual bad measurements do not influence the science results.

Our comparison between codes also highlights the necessity for any user of any light profile fitting code to verify that the galaxy parameters are actually well measured, ideally by comparing different codes on a small sample of objects and/or carrying out tests on simulated images, similar to the ones used in this work. More and more tools for creating such simulated images are becoming available, and the detailed description in this paper can serve as a guideline on how to create such data.

7 Effects of dust

As mentioned previously, we have not incorporated effects of dust attenuation in this work, since these are notoriously difficult to take into account and can only be properly quantified with radiative transfer calculations (Xilouris et al. 1999; Popescu et al. 2000; De Geyter & Baes 2014). This is well beyond the scope of this paper. Nonetheless, it is well known that the effects of dust attenuation are very strong in the ultraviolet, and, while in the optical bands they become very small for face-on disks, they can still produce a significant impact for non face-on galaxies (e.g. Driver et al. 2007). Dust attenuation does not only affect the spatially integrated SEDs (Tuffs et al. 2004; Natale et al. 2022), but also the surface brightness distribution of the direct stellar light (Möllenhoff et al. 2006; Gadotti et al. 2010; Pastrav et al. 2013a,b; Thirlwall et al. 2020). Because of this, the derived photometric parameters, like exponential scale-length, effective Sérsic radius, Sérsic index, and total luminosity, would be different from the intrinsic ones as obtained in the absence of interstellar dust (Pastrav et al. 2013a). In addition, the bulge-to-disk decomposition itself would suffer an extra effect due to dust, meaning that the decomposed disks and bulges in the presence of dust would be different from those derived if the galaxy would only have an attenuated disk or an attenuated bulge (Pastrav et al. 2013b). All these effects would depend on wavelength, inclination angle, dust opacity and bulge-to-disk ratio.

In order to account for dust effects one would need to use realistic, dust attenuated simulations of attenuated disks and bulges, instead of simple exponential or Sérsic functions, since no analytic functions exist to describe the complex modifications induced by dust. Nonetheless, such modified surface brightness distributions can be calculated with radiative transfer codes, and indeed, such simulated images already exist in the literature (Tuffs et al. 2004; Popescu et al. 2011). So the question then arises of why we would not attempt to use these simulated images to test our Galapagos-2/Galfitm codes. The problem is that then, instead of fitting only one or two analytic functions with a few free parameters, as done in this work, we would need to find the best-fit distribution from a large dataset of simulations corresponding to all combinations of parameters describing dust effects. In this paper we already showed that even simple function fitting is computationally a difficult task when dealing with large samples of galaxies. It becomes immediately apparent that complex distribution-fitting, is at present computationally impractical.

There is, however, a different solution to this problem, which has already been addressed and solved outside the development of a photometric code, by Pastrav et al. (2013a,b). The idea in those works was to provide corrections for dust effects on parametric models. Pastrav et al. (2013a,b) followed as closely as possible the procedures and algorithms used in photometric codes like the present one. It is just that instead of using observations of galaxies or the simple simulated images used here, they used radiative transfer simulations for which the input parameters describing the distributions of stellar emissivity and dust were known. By comparing the input values of the parameters that describe the simulations with the values of measured parameters that describe simplified distributions (exponential and Sérsic functions), such as those used in this work, they were able to quantify the degree to which the photometric codes underestimate or overestimate the intrinsic parameters of galaxies. These corrections are listed in Pastrav et al. (2013a,b). They are provided in the form of coefficients of polynomial fits to the corrections as a function of inclination, for different wavelengths and a large range of dust opacity.

We thus recommend using the dust corrections from Pastrav et al. (2013a,b)¹⁶, to adjust the ‘apparent’ (dust attenuated) photometric parameters derived from our code to ‘intrinsic’ (dust de-attenuated) parameters. The corrections could either be considered for a typical dust opacity (e.g. the Milky Way opacity from Natale et al. 2022) or could be derived on an individual basis from other methods (SED fitting, Balmer decrement, etc.).

8 Summary and Conclusions

In this paper we have presented and released Galapagos-2, a multi-band fitting code to be run on large surveys that automatically carries out profile fitting (both single-Sérsic and B/D fits) to all objects in the survey. We have presented a detailed analysis of the capabilities of Galapagos-2 and Galfitm, as well as a comparison with other works. As implications for single-Sérsic fitting have already been discussed in H13, we have focused here on presenting the advantages of multi-band fitting in the context of bulge and disk parameters, by comparing the code performance to the single-band fitting as usually carried out using other codes.

Throughout the paper, and using both simulated and real data, we have shown that multi-band fitting presents a significant and important improvement over single-band fits, not only in terms of parameter accuracy itself but also in terms of returning actual good results. Multi-band analysis allows a significantly larger sample to be used for scientific studies. Improvements on measured parameters are notably also made for the magnitudes and colours and SEDs of the different components (see Figs. 12 and 24). While in single-band fitting it is very difficult to derive a reliable SED for the individual components, multi-band fitting allows – at least on average – the SEDs of bulges and disks to be measured down to fainter magnitudes than is possible with other techniques, both on simulated and real data. Examining parameter accuracy as a function of B/T and galaxy magnitude, we were able to identify which galaxies can be reliably fit with two-component models. In this context, we were able to show that for objects with extreme B/T values (B/T < 0.2 or B/T > 0.8), neither component can be reliably fit; instead, the single-Sérsic fit provides a better match to the profile of a galaxy and the brighter component. As this requires measuring the B/T itself, this is somewhat of a circular argument, however, and cannot seemingly be avoided. In such scenarios, conservative values should be used.

We have further compared the outcome of these B/D fits to the results from other codes, where possible. We have found that, while they generally agree well, there were several specific differences, for example bulge sizes recovered by Gim2d. A detailed comparison to the literature is presented in Sects. 6.1–6.3.

Lastly, we have shown that using a Sérsic index of a single-Sérsic fit as a measure for the prominence of a bulge in a galaxy may not select high B/T galaxies as intended, as the correlation between the Sérsic index and the B/T ratio is at best weak.

We emphasise that this work is based on one-component (Häußler et al. 2013) and two-component profiles (this publication), and the inclusion of additional components beyond a simple B/D model (e.g. bar, halo, central light cusp, etc.) might be advisable depending on the science case, data at hand, spatial resolution, and data quality and depth. However, such profile fits often overfit the data and should hence, for now, only be carried out on an individual basis as they are by nature very hard to automate. To the best of our knowledge, no systematic attempt has been carried out in the literature, but the Galapagos-2 output files might serve as a good starting point for such an approach. Machine learning techniques by themselves, or in combination with profile fitting, might open the way to such an analysis, but they so far have not been proven to solve this issue entirely.

With this paper, we release the code¹⁷ to the community for wider use. In principle, this code is ready to be used on any large dataset. However, due to an IDL limitation of dealing with a maximum 16 galaxies in parallel, it would be advisable to develop a faster wrapper script for application to the large next-generation surveys, such as Euclid, LSST, or WFIRST. This could be done in a one-to-one translation of the IDL code, or by adapting the same principles, and consecutive testing. We invite any developer to further improve the IDL version of Galapagos-2 and submit a pull request on GitHub so that new features can potentially be implemented in the code. We are also happy to provide guidance for such development.

Finally, we are publishing two catalogues of fit results for all Gama objects, using the KiDS/VIKING imaging data (see Appendix F). These catalogues provide both single-Sérsic and B/D fits for 234 239 objects.

Acknowledgements

This publication was made possible by NPRP grant #08-643-1-112 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors. B.H. and M.V. were supported by this NPRP grant. S.P.B. is supported by an STFC Advanced Fellowship. We thank Carnegie Mellon University in Qatar and The University of Nottingham for their hospitality. E.J.J. acknowledges support from FONDECYT de Iniciación en Investigación Project 11200263. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/. This work is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 179.A-2004 and 177.A-3016. The authors thank Kalina Nedkova for publicly providing a python script to interpolate the Chebyshev polynomials to arbitrary wavelengths, and/or for restframe correction.

Appendix A Introduction

The following appendices serve as a resource for additional information regarding the findings above. They also provide a overview of the software and catalogues released with this paper. In Sect. C we describe in detail the creation of the simulated data used throughout Sect. 4. Sect. D will briefly discuss the effect of using different starting values for the B/D fits. Sect. E will introduce and explain features of Galapagos-2 introduced into more recent versions of the code, some motivated by the results presented above, including a brief section on possible further improvements in Sect. E.1 that might be implemented into the code in the future. Sect. F will describe the content of the catalogue of KIDS/VIKING fits - both single-Sérsic and B/D fits -released with this paper.

Appendix B Additional figures

Appendix C Image and galaxy profile simulations

As mentioned in Sect. 4.1, the simulated data used in this work were created using the same scripts described in Sect. 5.1 of H13 and Sect. 2 of H07. Naturally, however, given that two-component objects are the main target in this work, a different input catalogue for the component parameters had to be used. In this section we explain the process of deriving the simulated parameters in detail for those readers interested.

We started from a catalogue of results from fitting profiles to real galaxies (please see Sect. 5), in order to get a realistic distribution of galaxy parameters. In a first step, we created a catalogue of single-Sérsic objects equivalently to the way we did in H13. This ensured that the overall parameters of objects (e.g. the distribution of sizes and total magnitudes) stayed comparable to the analysis carried out in the real data, making sure that our findings on simulated data apply to real data. However, instead of passing this catalogue on to be simulated directly as was previously done, we introduced a few additional steps in order to transfer these objects into B/D systems to be decomposed and to optimise the simulations to their purpose in this work.

While many of the objects in H13 were too faint to be detected by Galapagos-2/SEXTRACTOR, objects just above the detection limit were required to realistically represent faint neighbouring sources. As the purpose of these simulations is somewhat less basic, we decided to include fewer very faint objects, in order to minimise the time needed for the simulations and fits. We have already shown in H13 and H07 that these faint objects do not have a significant effect on the fitting results of the brighter primary sources. As deriving a detection completeness was not part of this analysis, reducing the number of these objects around the detection threshold should not influence our results significantly. However, we tried to keep the object density in the images roughly constant, in order to not introduce additional challenges to the code by increased object overlap/confusion or – vice versa – make it too easy for the code to fit galaxies as no neighbouring objects are present.

While we ran Galapagos-2 on all detected galaxies – the object magnitude of simulated sources peaks around m_r = 20 – objects with a brightness close to the sky brightness are unlikely to be successfully decomposed by any method. As such, these faint galaxies were of somewhat less importance in this analysis, as only a few of them were likely to be fit well, especially in single-band fits, and indeed only presented results for galaxies at m_r,B+D < 19.5, well above the detection limit of these data.

We found that bright galaxies were very rare in this simulated catalogue, making an analysis difficult with ~ten objects in certain bins in certain plots (e.g. galaxies with m_r < 17.5 and successful fits in most/all single-band fits). We hence artificially increased the number of bright objects at the expense of very faint objects, to keep the overall number of galaxies constant, by forcing 1% of objects to have magnitudes in the range 11.5 < maɡ_r < 18. As we already established in H13 that the object density and the number of neighbours are not critical in these data (i.e. the objects do not strongly influence each other’s fits given a typical object density in GAMA), this adaptation is unlikely to change our conclusions, but provides better statistics for the objects we are most interested in.

We also decided to exclude very small objects r_e < 1 px for this work as these objects were not expected to be decomposed into components successfully and the performance of the single-Sérsic fitting has already been tested in H13.

Of all simulated galaxies, we kept 20% as single-Sérsic objects using the parameters chosen at this stage, in order to be able to test any single- versus multi-component classifiers from our simulations. While this is a very important issue, it has not yet been solved satisfactorily. Using these simulations, we are in general able to examine this problem at least in the case of idealised, symmetrical, smooth galaxy profiles. As this is beyond the scope of this paper, we leave this analysis to a later paper and instead only focus on separating the two components of the two-component galaxies themselves.

In an additional step, we transferred the remaining 80% of objects into B/D systems. This is done following these rules for each parameter:

Position: The x and y centre coordinates of the bulge and the disk are kept to be identical in the simulations, we simply kept the randomly chosen position from step 1. While in real galaxies, a slight offset between the two components could in principle be seen, either because of an intrinsic offset or due to dust in the galaxy affecting the different components in different ways, we decided to avoid this additional confusion. This is hence directly in line with the setup used when fitting these images with Galfitm, where the position of bulge and disk are constrained to be identical.

Half-light radius: For the purposes of this test and analysis, we decided to ignore gradients in colour and stellar populations within the galaxy components, mainly for simplicity, but also as in first order (only) this is a good approximation of reality, as has been discussed above. This means that both r_e and Sérsic index n are to be kept constant with wavelength in the simulated galaxies in each component. The sizes of the galaxy components were hence simply calculated from the overall half-light radius as derived in step 1 as the median of the values at the different wavelengths. The distributions used for this purpose followed the distributions found by fitting B/D profiles to the real galaxies (see Sect. 5 for discussion of these fits). We used the following distributions: For disks, we used: $$r_{\text{e,D}}=median\left(r_{\text{e,SS}}\right)*\left(ran{d}_{\text{D}}*0.7+1.2\right),$$(C.1)

where rand_D is a random number (Gaussian, FWHM of 1, centred on 0). This effectively creates most (68%) of disks within 0.5 * median(r_e,SS) < r_e,D < 1.9 * median(r_e,SS), which is typically found in our real galaxies. For bulges, we used $$r_{\text{e,B}}=median\left(r_{\text{e,SS}}\right)*\left(ran{d}_{\text{B}}*0.5+0.85\right),$$(C.2)

with rand_B being a different random number compared to the disk value (i.e. bulge and disk size are determined independently). This effectively creates most (68%) of bulges within 0.35 * median(r_e,SS) < r_eD < 1.35 * median(r_e,SS). We note that we explicitly do not restrict component sizes to r_{e B} < r_e,D as we wanted to test the code for all cases, even though galaxies with r_e,B > r_e,D are rarely found in nature. A similar constraint is not used within Galapagos-2 for the fitting process either.

We do, however, impose a less strict constraint at r_{e B} < 2.5 * r_e,D, in order to avoid extreme values. We further impose our usual minimum and maximum sizes to 1 [px] < r_e,B/D < 400[px]. If any of these limits are violated, we simply derived r_{e B} and r_{e D} again from r_e,SS, using different randomisation.

Magnitudes and Colours: The magnitudes in the simulated objects were by far the trickiest of the parameters to decide on, as they could potentially have strong effects on both the results from the test and its implications on real galaxies. MEGAMORPH will perform best if the two components have very different SEDs, making a separation easier than if both components showed the same colours (i.e. identical SEDs). We created several different sets of simulations, with tens of thousands of simulated galaxies each, before we settled on a simulated dataset to be used in this paper. Previous sets basically lead to the same conclusions, but with lower statistical significance and in a way that was more complicated to present in simple plots (e.g. due to the lack of bright objects, or each disk and bulge showing an individual colour). For all simulations, we used the distribution of the total magnitudes m_SS from step 1 and used a random B/T ratio (uniformly between 0 and 1) to divide the flux of each galaxy in the z band (in the middle of the covered wavelength range) between the two galaxy components. As this split was carried out in the z band, we gave both components a boost of 0.2 mag to account for the typical magnitude difference between the r band used in step 1 and the z band used here. As we wanted to concentrate this analysis on generally brighter galaxies, this worked in our favour.

After defining the z-band brightness this way, we needed to define an SED for disks and bulges, to derive their respective magnitudes in the other bands. The result of this approach is that all disks and all bulges in the simulations actually do show the same SEDs, respectively. This makes a comparison of parameters – especially the magnitudes – much cleaner and easier to present. While leading to similar conclusions in our internal analysis, they show results at higher statistical significance than previous sets of simulations, which is why we present them here.

We settled on the approach to use a median colour of the bulges and disks (i.e. ‘typical’ bulges and disks) from the Galapagos-2 fits on our real galaxies (see Sect. 5), respectively, restricting the sample to a tight subset of bright, large galaxies with high-n bulges, where a B/D decomposition is most likely to be successful. In the end, for disks we used (colours compared to the r band)

• u-r = 2.195

• g-r = 0.910

• i-r = -0.429

• z-r = -0.683

• Y-r = -0.705

• J-r = -0.655

• H-r = -1.047

• K-r = -0.871

and for bulges

• u-r = 3.384

• ɡ-r = 1.442

• i-r = -0.577

• z-r = -0.917

• Y-r = -1.123

• J-r = -1.543

• H-r = -1.820

• K-r = -2.047.

These values are somewhat extreme in their colours, especially in the low S/N bands, for example u, but produced images that looked overall more similar to the real GAMA images than previously tried datasets. Finally, the magnitudes of each component are modified by adding Gaussian noise with a 0.1 standard deviation for each component, in order to simulate observational noise and re-introduce some variety into the component colours.

Sérsic indices: Although for testing purposes this might have been a valid approach, we wanted to avoid the oversimplified case that all bulges have n_B == 4, making a successful decomposition more (and unrealistically) likely. Instead, we wanted bulges to show a variety of n_B values as this was a more realistic approach and tested the codes for a wider variety of objects. We wanted to test whether our codes can successfully decompose galaxies even if they contained low-n bulges (i.e. pseudo-bulges), with a Sérsic index similar to galaxy disks, as long as the two components showed different colours. For simplicity, however, we decided against randomly choosing a Sérsic index from ‘real’ bulges and instead chose to use a more simple (although artificial) distribution of Sérsic indices: $$n_{\text{B}}=rand*2+4,$$(C.3)

with rand again being a random number (Gaussian, FWHM of 1, centred on 0). This results in a n-distribution centred around 4, with a FWHM of 2. We restrict the n_B value as 0.5 < n_B < 8 as values outside of this range are rarely found in nature. We simply re-calculate the value if this constraint is violated.

Equivalently to the r_e values, by deciding to ignore gradients in colour and stellar populations within the galaxy components, Sérsic indices have to be constant with wavelength (i.e. n_B is the same at all wavelengths).

For disks, we simply use n_D == 1 at all wavelengths, as it is well established that the Sérsic index of disk components is very close to this value, while the bulges in real galaxies show a much wider spread.

Axis ratio: For the axis ratio (q) of each component, we also used an artificial distribution for the two galaxy components. We chose a random orientation, Θ, of the galaxy plane (randomly uniform in cos(Θ)) and assumed an intrinsic thickness of 0.18 for disks and 0.5 for bulges, which we used to calculate the intrinsic axis ratio of each component given Θ. These values are assumed to be constant with wavelength using $$q_{\text{D}}=\sqrt{{0.18}^2+\left(1-{0.18}^2\right)*\cos {\left(\Theta \right)}^2}$$(C.4)

for the disk and $$q_{\text{B}}=\sqrt{{0.5}^2+\left(1-{0.5}^2\right)*\cos {\left(\Theta \right)}^2}$$(C.5)

for the bulge.

Position angle: The position angle of the galaxy is chosen randomly and assumed to be the same for both bulges and disks. While there are physical reasons why the position angle could be (at least somewhat) different for the different galaxy components, we chose to ignore these special cases and instead stick to the much more common case that bulge and disk semi-major axes are aligned with each other. Technically, this is the more challenging case for the fits as two elongated components with very different position angles should be easier to separate. However, we note that no such constraint is being imposed during the fit.

The result of this procedure is a catalogue of objects that can be turned into ‘observed’ images using the same scripts as used in H13. These images contain a mixture of single-Sérsic (20%) and two-component objects (80%), which leaves us the possibility of testing any one- versus two-component classifiers in future work.

Following the scripts used in H13, the resulting images used in this work include Poisson noise, and the galaxies are added into an image made up of empty sky patches from real GAMA imaging. This leads to more realistic noise properties in the image and allows a realistic estimation of the recovery of fitting parameters.

Appendix D Varying starting parameters

During the development of Galapagos-2, and specifically during the implementation of B/D fitting, we experimented with different starting parameters for the B/D fits, especially in r_e and n, as the other parameters in our setup are usually free in each band (magnitudes) or held constant at all wavelengths, and usually easier to fit (q and PA). These choices should be presented in the simulations, and in the choice of starting parameters for the fit.

Throughout most of this work, we have used a GALAPAGOA-2 version that uses r_e,D = 1.2 * r_e,SS (lower limit 1 pixel) and r_e,B = 0.3 * r_e,SS (lower limit 0.5 pixels) for component sizes and n_D == 1 and n_B = median(n_SS) (lower limit 1.5) for Sérsic indices for the disk and bulge, respectively.

We, however, tried different combinations of these and the following starting values:

• r_e,B = 1.0 * r_e,SS

• r_e,B = 0.5 * r_e,SS

• r_e,D = 1.0 * r_e,SS

• n_D = 1 (as starting value, not fixed value)

• n_B = 4 (as starting value, not fixed value)

• n_B = 1 (as starting value, not fixed value).

Overall, we found limited impact on the fitting results derived in multi-band fitting derived once extreme values were avoided, at least in bright objects. In single-band fitting, however, there was indeed a difference between the individual test runs, particularly in the number of objects that return a good result at all, as separating bulges and disks becomes increasingly difficult in those data and convergence of B/D fits becomes challenging. Maximising this number became one of the goals of an ideal solution, although the impact on multi-band data was much smaller. However, we saw no reason why multi-band fits should not also benefit from such a choice of starting parameters, even if it might only be for the number of iterations required in the fit, hence speeding up the analysis.

In this section we show examples of the differences between three of these runs on the real GAMA galaxies. Throughout, we label the three different runs as ‘Fit1’, ‘Fit2’, and ‘Fit3 as follows.

Fit1 uses as starting parameters of the B/D fit: r_e,D = r_e,SS (lower limit 1 pixel) and r_e,B = r_e,SS (lower limit 0.5 pixels). Fit2 uses r_e,D = r_e,SS (lower limit 1 pixel) and r_e,B = 0.5 * r_e,SS (lower limit 0.5 pixels). Fit3 (used throughout most parts of this paper) uses r_e,D = 1.2 * r_e,SS (lower limit 1 pixel) and r_e,B = 0.3 * r_e,SS (lower limit 0.5 pixels).

In Fig. D.1 we present the comparison of these fits on the example of the bulge size in the r-band r_e,r,B on the left and the disk size r_e,r,D on the right side. The labels in each panel indicate which values are compared, and the black number shows the number of objects in each panel (different sample sizes largely reflect the fact that different samples of galaxies were fit in these tests). In each panel, we show the ratio of the derived parameters (cat1/cat2) as a function of the single-Sérsic magnitude in cat1.

From these plots, we can see that Fit2 and Fit3 largely agree on the bulge and disk sizes, relatively few galaxies show sizes that disagree by more than 10% (15% of bulges, 12% of disks) between the different runs (bottom panels). However, when comparing either of these runs to the fits from Fit1 (top and middle panels), it is obvious that the disagreement between the runs is pronounced and systematic. Bulges are systematically fit larger in Fit1, disks smaller. We attribute this behaviour to the starting values for the sizes of bulge and disk being identical in this run, indeed making the separation of the components harder. In this case, besides the obvious difficulty in separating the two components, it more likely happens that bulge and disk ‘flip’ position, which explains the larger bulges and smaller disks.

We find similar and consistent behaviour in other fitting parameters like Sérsic index (no plot shown). However, once the bulge is started at a smaller size than the disk, the exact starting value seems to become less critical and fitting values become more consistent, although some offsets can still be seen for faint galaxies. This exercise can serve as a indication for the accuracy of the fitting values (i.e. to what brightness galaxies can be reliably decomposed) and as a warning to not over-interpret the finding of disk and bulge parameters found in faint galaxies.

We further tried using different constraints on the fitting values, including, and most importantly, the position of the object. This was allowed to vary either over the entire image or within 0.5 * r_e,SS, but this choice had a very minor impact on the fitting results. A larger impact was found when we did not constrain the disk and bulge centres to be identical. However, these cases usually turned out to be bad fits in which one of the components tried to fit the residual of a neighbouring object, so we found it important to implement the constraints to keep disk and bulge centres identical in Galapagos-2, which prevents the vast majority of these cases. It also makes these cases easier to detect, as in this setup these objects require extreme axis ratios in order to minimise the residual of neighbouring objects (see Fig. 23).

Additionally, one could individually test the impact of using different DOF for the polynomials used for each parameter. As this is a user choice in GALFIT and Galapagos-2 and far beyond the scope of this work, in which a truly complete analysis of all parameters is impossible, we only present the results for the limited models that we have used in this paper.

Appendix E Additional features in the most recent version of Galapagos-2

As this paper serves as a code-release paper for Galapagos-2, we deem it important to summarise further improvements of the code. in case of further questions, or in case of requests for additional features, please do not hesitate to contact us. As was previously mentioned, we invite users to adapt Galapagos themselves by cloning the GitHub repository¹⁸ and feeding any new features back into the main code via a pull request. in this chapter, we present a brief list of the 12 main improvements of the code since version 2.2.7, which was used throughout this work.

First, a potentially critical bug was fixed in which, depending on why galaxies were fit as neighbours, they were put on a fixed position which was off by 1 pixel in both x and y. This had minor impact on the data that we have seen, but was potentially critical in dense fields with many stars, as the artificial residuals created would impact the fitting parameters of the main target.

Second, a new parameter, primary−χ², was introduced, which measures the χ² value in a area around the main target only, by masking out all other objects (no actual fit is done in the process). This parameter could potentially be used as a measure of the ‘goodness’ of the fit and decide which of the models -single-Sérsic or B/D – better fits the data. However, we have not fully followed this approach through in this work. This selection of the better model is a largely unsolved problem in astronomy, for which several approaches have been attempted, including statistical methods (e.g. Lackner & Gunn 2012; Simard et al. 2011), visual morphologies (e.g. Lintott et al. 2008), machine learning (e.g. Huertas-Company et al. 2015). None of these methods, however, has yet been found to work in a universal fashion, so it can be easily employed to any dataset. We briefly looked at cross-validation to measure how good the fit (in some areas of the images) fits the data in other areas of the images. As this method is very CPU intensive, requiring several hundred fits to the same galaxy, it is unfeasible for large samples. However, it is a promising approach for small samples or individual – well resolved – galaxies and seems to work well on those. our approach here of presenting a primary −χ² value attempts to give a good starting point for statistical approaches, for example BIC (Bayesian information criterion), and AiC (Akaike information Criterion), methods, and has been employed by Nedkova et al. (in prep.) to analyse bulges and disks of galaxies in HFF (Lotz et al. 2017) and Candels data.

Third, several parts of the code have been significantly sped up. Fourth, several bugs have been fixed to avoid crashes of the code in very specific circumstances. For example, the code can now deal with surveys where the footprint of each image is different. Earlier versions of the code would crash if a postage stamp for an object in one band contains only masked pixels.

Fifth, several checks have been implemented at the start of the code, for example whether all data exists. previously, it was possible due of a typo (e.g. in the name of the PSF used) to produce a fit without using a PSF, which was hard to notice after-wards. This is no longer possible, as Galapagos-2 checks that all PSF files exist.

Sixth, the output of the code into the command line has been changed, so it is easier to monitor the progress. seventh, several parameters for object counts now use LONG instead of INTEGER in order to deal with individual images with more than 32,767 detected objects. Eighth, the option \galfitoutput was added, which saves the GALFIT terminal output for each object in case a user wants to examine them at a later stage. Ninth, the option \sex_skip was added to avoid re-running SEXTRACTOR on all the individual tiles, for example in cases where the code previously crashed during the stage of combining those catalogues.

Tenth, the lower constraint for the axis ratio of objects was changed from 0.0001 to 0.01, to avoid very extreme cases that, due to sampling issues, can create a line of dots across the image, rather than a smooth profile. As both values are very extreme, no significant impact on the fits has been expected, other than avoiding these specific issues in very isolated cases.

Eleventh, while Galapagos-2 includes a experimental feature to use super-computing facilities, the limitations of IDL (16 objects in parallel) do not allow this to be used efficiently, and to our experience actually slows the fits down while creating additional work for the user. While this feature has not technically been removed from the code, we stopped any further development of this feature.

Finally, several new utility scripts have been created to make using Galapagos-2 and its data products more efficient. These are somewhat experimental and we do not guarantee that they work perfectly, but they can be used for a variety of useful tasks. (please check the description */src/utilities/Utility_README.md for details.)

E.1 Prospects, plans, and further improvements

It should be stressed that Galapagos-2, in addition to the obvious fits it carries out, generally also works well as a setup for further analysis by readily creating all necessary data and GAL-FITM setup files. Adding another component (point source, bar, etc.) and enabling Fourier or rotational modes becomes a trivial undertaking for specific objects and science cases. As implementing and testing such features is a major task and not easily done in such a general, fully automated code, we do not currently have plans to implement these features. Instead, we encourage users to use the output files from Galapagos-2 to create their own add-on routines for further analysis of specific objects.

However, several improvements of Galapagos-2 are, of course, possible, and some have already been mentioned in Sect. 2.1.2. In this appendix we discuss four of the more obvious improvements. We do not currently have plans to implement any of these, but if users request them, we are happy to work together with them to improve the code in this direction and implement such new features.

First, one obvious possible improvement can be made by using different starting values for the magnitudes. The starting values for magnitudes in different bands are currently defined by using the SEXTRACTOR on the detection image, and applying a typical offset to the other bands. such an approach is obviously not ideal and instead one could run SEXTRACTOR in dual image mode on all images, in order to derive better starting values for the single-Sérsic fits in all bands. It is unclear how this idea could be used for the B/D fits, but it is possible that tweaking the starting values for those fits (currently starting from the same SED for disk and bulge that was derived in the single-Sérsic fit) would improve the robustness and accuracy of the B/D decomposition (i.e. by assuming redder SEDs for bulges and bluer SEDs for disks). This is the most likely new feature to be implemented, but as it would add significant CPU time and a benefit can not currently be estimated, we have no immediate plans to do so.

Second, better starting values for other parameters are similarly possible, for example by using some measure of light concentration to estimate the Sérsic index of the single-component fit. Such an attempt would require extensive testing, however, as these parameters have been found to not correlate well in the past, so it is currently not planned.

Third, depending on the fit result, one could imagine a second fitting step. For example one could re-start the fit of ‘failed’ fit objects using different starting values, even multiple times (e.g. as employed in Lange et al. 2016). Alternatively, one could employ a scheme in which one uses the single-Sérsic axis ratio to create azimuthally averaged 1D profiles and fit those in order (especially small bulges might be more easily detected in such a scheme) to derive better starting values for another 2-D fit using Galfitm, possibly even with a reduced set of free parameters. Such a system, however, would require intensive development and testing.

Finally, a new method for the determination of the sky values in data was presented and discussed by Ji et al. (2018). While we have no current plans to implement this new method as the current method seems to be working well, this is, of course, possible and might further improve the fits.

The ultimate improvement of Galapagos-2 and any other, similar code would be to combine them with their individual strengths and different features into one code and skipping the development of several, largely duplicate codes altogether as a community. A new modular fitting code wrapper – written in a more flexible (and free) coding language than IDL (e.g. Python, Julia, R, or others) – could be implemented in a modular fashion, in which a user could choose to use a multitude of source detection codes (e.g. SEXTRACTOR or PROFOUND), different fitting codes (GALFIT, Galfitm, PROFIT, GIM2D, IMFIT, or others), different setup versions, and so forth. Such a code would be the ultimate tool for galaxy profile fitting, but can only be developed by the community as a whole.

Appendix F Catalogue release and fits of KIDS/VIKING data

The work presented in this paper has been carried out with SDSS/UKIDSS data, in order to be consistent with the test carried out in H13. However, more recent KIDS/VIKING data are much deeper, as discussed in Sect. 6.2, and have recently been fit using Galapagos-2. With this paper, we also release two catalogues of Galapagos-2 fits of these deeper data, as we think this would be highly beneficial for the larger community.¹⁹ In this section we briefly describe the catalogues released.

The catalogues released contain both single-Sérsic and two-component B/D model fits to the 2D surface brightness distributions of 234,239 objects in re–ed KIDS/VIKING data – as provided by GAMA as large mosaics – using Galapagos version 2.3.2. and Galfitm version 1.2.1. This represents the full sample of objects presented in Baldry et al. (2018), with a few exceptions described below.

Re-SWARP-ed imaging data from the GAMA G09, G12, and G15 survey regions are being used. These data include four-band optical (ugri) imaging from KIDS (de Jong et al. 2013; Kuijken et al. 2019) plus five-band near-infrared (ZYJHK) imaging from the VIKING survey.

A number of bright stars and nearby areas are masked out of the data by visually identifying them and setting the values in the weight images in an area around those stars to a value of 0, so both SEXTRACTOR and Galfitm ignore those areas during the object detection and fit, respectively. Similarly, areas where the data were saturated have been masked out generously as to not influence the fits of nearby objects. This removes some objects from the sample that would be difficult to be fit due to their proximity to very bright stars. other than these adaptations, we only cut the large mosaics provided by GAMA into images of ~8500x8500 pixels with overlap (termed ‘tile’), in order to make the file sizes more manageable for Galapagos-2, which further cuts the images into postage stamps centred on each object.

We use SEXTRACTOR on the r-band image to detect sources. This image is deemed the deepest data available amongst the available data. While objects with extremely red colour could potentially be missed, it is sufficient to easily detect any object with GAMA spectroscopy. As was found to be ideal for Galapagos-2, this is done in a two-stage (hot and cold) process process. ‘Cold’ detects the brightest objects without splitting them up, ‘hot’ detects the faint objects. Both catalogues are then combined as outlined in the Galapagos manual (Barden et al. 2012).

The PSFs used for these fits are obtained by using PSFEX to create a large number of PSFs for each input tile, choosing the closest PSF in RA and Dec. to the targeted object for each fit. These PSFs were constructed on a grid of 10x10 (u band) and 20x20 (ɡrizYJHK bands) PSFs per tile of ~8500x8500 pixels in size as outlined in Sect. 6.4. This means that PSFs are available on a regular grid with a density of one PSF every ~400 or ~800 (u-band) pixels. Detailed setup values can be provided upon request.

We fit single-Sérsic profiles to all objects that have a redshift in the GAMA spectra catalogue ‘SpecAll27’ as provided by the GAMA survey, at a search radius of 5″(any object within this radius is fit). Additionally, we fit all objects within 100 arcsec-onds of any of these objects, but brighter than said object. This is ideal for Galapagos-2 as objects are dealt with in brightness order. By leaving these objects out, the sample selection could potentially change the fitting results. In the fitting process itself, fainter neighbours are taken into account as normal, but these will never become primary objects. Successive two-component fits are only carried out on the GAMA objects themselves, which is sufficient (see discussion above).

The DOF used for the individual parameters are as discussed in Sect. 2.4 and in H13. Single-Sérsic fits use free magnitudes, second-order polynomials (three free parameters) for half-light radii, and the Sérsic index. All other parameters are constant with wavelength. Two-component fits use full freedom for magnitudes and values constant with wavelength for all other parameters (where n_D == 1 for disks). The second provided catalogue contains identical single-Sérsic fits but additionally uses n_B == 4 for bulge Sérsic indices. This catalogue is released, but has not been tested in detail in this paper. Results from Nedkova et al. (in prep.) to HFF and Candels data suggest that these bulge parameters are actually more reliable.

We make no attempt to quantify the quality of an individual fit, nor do we indicate whether the single or double component fit is to be preferred for a given galaxy. When using the data from these catalogues, specifically the cleaning of the sample, the guidelines discussed in Sect. 3 and in Vulcani et al. (2014) and Kennedy et al. (2016a) should be followed.

In Fig. F.1, we show a comparison of the main fit parameters for the two catalogues for a bright sample (SEXTRACTOR MAG_BEST< 18.5) and CLASS_STAR < 0.8. The plot looks very similar for fainter objects, but there are too many objects, so that features become difficult to see. As the same single-Sérsic fits are used for both fits, we do not show their comparison as we did in previous figures, as they would exactly show a one-to-one line. In general, the agreement of the B/D parameters is good in all parameters, even the B/T ratio shows a very clear and relatively tight correlation. There is, however, a concentration of objects with B/T ~ 0 in the fits with fixed n_B (on the left side of the B/T panel). While we did not investigate those objects further, we can speculate from the findings above that these are at least in part objects in which the bulge with free Sérsic index fits the disk profile, and the fit with n_B = 4 correctly finds only a very weak bulge. In the panel comparing the n_B values (although for one of the fits this was fixed at n_B = 4), it is interesting to see that most of the pure bulges plotted in red indeed show Sérsic indices in the range of 2.5 < n_ss < 5 as one would expect for such bulge-dominated galaxies. The user should be reminded that this selection, as defined in Sect. 6, is done on the fitted B/T ratio and no Sérsic index selection is used.

There is also a prominent bulk of objects at low q_B values in both fits. This has already been discussed in Sect. 6.2 and can be further reduced by a stricter cut on CLASS_STAR, which also somewhat reduces the scatter in the other panels, indicating that most outliers are indeed very small objects.

In Table F.1, we present a short summary and explanation of the parameters in the released catalogues. As a multitude of parameters are very similar, and a full description would be very long, we only give examples. However, following the nomenclature, other parameters should easily be understood. Any user of either of these catalogues should be careful to follow the advice on deriving clean object samples given in Sect. 3. In summary, the catalogues contain: the GAMA CATAID; several useful SEXTRACTOR parameters; several additional Galapagos-2 parameters; general parameters (single numbers) and Galfitm parameters (nine-element arrays); single-Sérsic parameters indicated as _galfit, BD parameters as _BD, bulge parameters as _B, and disk parameters as _D; parameters (e.g. MAG_GALFIT_BAND) and errors as derived by Galfitm (e.g. MAGERR_GALFIT_BAND); and sizes (half-light radii) given in arcseconds.

Finally, the catalogues also contain two sets of fit parameters, one containing the fit values at the effective wavelengths of the filters used _BAND²⁰ and the other containing the fit parameters for the Chebyshev polynomials _CHEB themselves, where they make sense (i.e. not if values are constant with wavelength). The latter can be used to easily interpolate between values to other wavelengths, but a user should be aware of Sect. 2.4 where interpolation of polynomials with high DOF are being discussed. Where _BAND values are given, they are in the order of r,u,ɡ,i,z,Y,J,H,K band, due to the r band serving as the main band in the GALAPAGOA-2 run²¹.

The catalogues provide fit values even for those objects where the fit crashed. These objects can easily be identified by FLAG_GALFIT and FLAG_GALFIT_BD for single-Sérsic and B/D fits, respectively (see the details in the table). Their parameter values and errors are set to obvious standard values: 0 for positions, −999. for magnitudes, and −99. for sizes and Sérsic indices, generally 99999. for error values. It should be noted that there can be objects with error values of 99999. despite having fit values provided. For these few objects, the error estimation within Galfitm has failed to provide a value, while the fit itself did converge on a final solution. Visually, we have not been able to identify anything special about such objects.

References

---

All Tables

All Figures