Properties of extragalactic thick discs recovered from ultra-deep Stripe82 imaging

TL;DR

This study investigates the structure of thick galactic discs using ultra-deep imaging, emphasizing the importance of PSF effects and providing detailed measurements of disc parameters in five edge-on galaxies.

Contribution

It presents a method to accurately model thick discs in ultra-deep images by accounting for PSF effects, extending previous surface brightness studies to deeper levels.

Findings

PSF effects significantly influence low surface brightness profiles.

Thick disc components are necessary for accurate galaxy modelling.

Deeper imaging reveals more complex disc structures and better constrains formation theories.

Abstract

We study the structure of thick discs in extraordinary depth by reaching a surface brightness limit of $\mu_{r_{\rm deep}} \sim$28.5$-$29~mag arcsec$^{-2}$ with combined $g$,$r$,$i$ band images from the IAC Stripe 82 Legacy Project. We present the characterisation of the thick discs in a sample of five edge-on galaxies. A study of the radial and vertical surface brightness profiles is presented by comparing our data with point spread function (PSF) deconvolved models. Our method begins with an analysis of the background and masking processes. Then we consider the effects of the PSF through galaxy modelling. The galaxy disc components are fitted considering that the thin and thick discs are two stellar fluids that are gravitationally coupled in hydrostatic equilibrium. We find that effects due to the PSF are significant when low surface brightness is reached, especially in the vertical profiles, but it can be accounted for by careful modelling. The galaxy outskirts are strongly affected by the faint wings of the PSF, mainly by PSF-redistributed light from the thin disc. This is a central problem for ultra-deep imaging. The thick-disc component is required to reach satisfactory fit results in the more complex galaxies in our sample, although it is not required for all galaxies. When the PSF is ignored, the brightness of these structures may be overestimated by up to a factor of $\sim$4. In general, our results are in good agreement with those of previous works, although we reach deeper surface brightness levels, so that the PSF effects are stronger. We obtain scale heights and mass ratios of thin and thick discs ($z_{{\rm t}}$, $z_{{\rm T}}$, and $M_{{\rm T}}/M_{{\rm t}}$), which provide excellent agreement with previous studies. Our small initial sample provides evidence for aspects of a wide variety of formation theories for the thick discs in disc galaxies.

Tables

Table 1. Table 1: List of the 51 rejected galaxies from the shortlist of 56 galaxies after performing the first selection using the HyperLeda database. The name and the primary reason why we rejected the galaxy from the final data sample are indicated in each column.

Galaxy ID	Rejection Reason
2MASX J2231+0026	No distance
IC 0330	Not edge-on
IC 1756	Foreground objects
NGC 0585	Dust lane
NGC 0955	Dust lane
NGC 7648	No disc galaxy
PGC 011579	Not edge-on
PGC 011777	No distance
PGC 012721	No distance
PGC 066048	Foreground objects
PGC 072794	Not edge-on
PGC 091784	No distance
PGC 092868	Foreground objects
PGC 1149944	No distance
PGC 1157026	Not edge-on
PGC 1171976	Interaction
PGC 173375	No distance

Table 2. Table 2: Final data sample properties and physical parameters.

Galaxy	Morph.	D	Incl.	Mabs	$v_{rot}$	$M_{\mathrm{dyn}}$	M⋆	$d_{25}$	${\mu }_{r_{\mathrm{deep}}}$
ID	type	[Mpc]	[deg.]	[B-mag]	[km/s]	[M⊙]	[M⊙]	[arcmin]	mag arcsec-2
NGC 0429	S0-b?	93.5	90.0	-20.29 $\pm$ 0.353	-	-	4.4-8.3$\times$1010	1.51 $\pm$ 0.11	28.5 $\pm$ 0.1
NGC 1032	S0-sp / E(d)7	38.7	88.1	-20.75 $\pm$ 0.43	261.0	3.3$\times$1011	-	3.55 $\pm$ 0.11	29.7 $\pm$ 0.3
UGC 00931	Sc	31.9	78.1	-17.46 $\pm$ 0.678	50.7 $\pm$ 2.1	3.1$\times$109	-	1.02 $\pm$ 0.11	29.3 $\pm$ 0.2
UGC 01040	S0-a?	44.1	90.0	-19.62 $\pm$ 0.37	-	-	1.3-2.8$\times$1010	1.17 $\pm$ 0.12	29.0 $\pm$ 0.2
UGC 01839	Sdm:sp	24.7	90.0	-17.40 $\pm$ 0.517	65.1	5.9$\times$109	-	1.20 $\pm$ 0.12	28.0 $\pm$ 0.1

Table 3. Table 3: Physical parameters of the thin and thick discs of each galaxy in our sample.

Phys. Param.	NGC   0429		NGC   01040		UGC   00931		UGC   01839		UGC   1032
$z_t$ [kpc]	1.14	1.38	0.31	0.36	0.31	0.35	0.31	0.28	1.74	1.8
$z_T$ [kpc]	5.15	5.75	1.13	0.91	0.61	0.53	0.6	0.47	8.82	8.81
$z_T/z_t$	4.51	4.17	3.64	2.53	1.96	1.51	1.93	1.68	5.07	4.89
${\rho }_{0T}/{\rho }_{0t}$	0.09	0.04	0.16	0.09	0.24	0.33	0.29	1.31	0.11	0.13
${\Sigma }_T/{\Sigma }_t$	0.34	0.15	0.5	0.22	0.45	0.46	0.51	2.19	0.42	0.49
${{\rm Y}}_T/{{\rm Y}}_t$	-	0.99	-	1.07	-	1.97	-	1.46	-	0.56

Equations

$$\frac{d^2{\rho }_i}{d{z}^2}=\frac{{\rho }_i}{{\sigma }_i^2}\left[-4\pi G({\rho }_{\mathrm{t}}+{\rho }_{\mathrm{T}}+{\rho }_{\mathrm{g}})+\frac{d{K}_{\mathrm{DM}}}{dz}\right]+\frac{1}{{\rho }_i}{\left(\frac{d{\rho }_i}{dz}\right)}^2,$$

$$\frac{\mathrm{erro}{\mathrm{r}}_{\mathrm{bin}}}{\#count{\mathrm{s}}_{\mathrm{bin}}}>0.2.$$

$${\log }_{10}(M/L)=a_{\lambda }+b_{\lambda }\times \mathrm{colour},$$

Full text

Aims. We study the structure of thick discs in extraordinary depth by reaching a surface brightness limit of $\mu_{r_{\rm deep}}\sim$28.5$-$29 mag arcsec-2 with combined $g$,$r$,$i$ band images from the IAC Stripe 82 Legacy Project.

*Methods. *We present the characterisation of the thick discs in a sample of five edge-on galaxies. A study of the radial and vertical surface brightness profiles is presented by comparing our data with point spread function (PSF) deconvolved models. Our method begins with an analysis of the background and masking processes. Then we consider the effects of the PSF through galaxy modelling. The galaxy disc components are fitted considering that the thin and thick discs are two stellar fluids that are gravitationally coupled in hydrostatic equilibrium.

*Results. *We find that effects due to the PSF are significant when low surface brightness is reached, especially in the vertical profiles, but it can be accounted for by careful modelling. The galaxy outskirts are strongly affected by the faint wings of the PSF, mainly by PSF-redistributed light from the thin disc. This is a central problem for ultra-deep imaging. The thick-disc component is required to reach satisfactory fit results in the more complex galaxies in our sample, although it is not required for all galaxies. When the PSF is ignored, the brightness of these structures may be overestimated by up to a factor of $\sim$4.

*Conclusions. *In general, our results are in good agreement with those of previous works, although we reach deeper surface brightness levels, so that the PSF effects are stronger. We obtain scale heights and mass ratios of thin and thick discs ($z_{{\rm t}}$, $z_{{\rm T}}$, and $M_{{\rm T}}/M_{{\rm t}}$), which provide excellent agreement with previous studies. Our small initial sample provides evidence for aspects of a wide variety of formation theories for the thick discs in disc galaxies.

1 Introduction

Most if not all disc galaxies have a thin (or classical) disc and a thick disc (Burstein, 1979; Tsikoudi, 1979; Yoshii, 1982; Fuhrmann, 1998; Dalcanton & Bernstein, 2002; Yoachim & Dalcanton, 2006; Fuhrmann, 2011; Snaith et al., 2014; Comerón et al., 2011, 2012, 2014, 2018). In most models thick discs are thought to be a necessary consequence of the disc formation process and/or of the evolution of the galaxy. However, because of their faint surface brightness ($\mu_{g}\sim$26.0 mag arcsec*-2*), the detailed study of the physical properties of thick discs has thus far been very difficult. The goal of this paper is to shed light on this problem by characterising the thick-disc component in an initial sample of five edge-on nearby spiral galaxies in order to compare the properties of their thin and thick discs.

The origin of thick discs is still a matter of debate (Yoachim & Dalcanton, 2008, and references therein). Although there are many different approaches, three main families of models have been proposed (e.g. Comerón et al., 2011). The first family of models is secular. Here, the original idea was that the thick disc was caused by vertical heating and/or radial migration of stars (e.g. Schönrich & Binney, 2009; Loebman et al., 2011), but radial migration was later found to be less likely as a mechanism (e.g. Minchev et al., 2012; Vera-Ciro et al., 2014) and a revised scenario has been proposed where thick discs form from the nested flares of mono-age populations (Minchev et al., 2015; Martig et al., 2016). In the second family, the disc was always thick from its beginning, as a result of a high-velocity dispersion of the interstellar medium at high redshift; the thin disc is accreted with a lower velocity dispersion (e.g. Elmegreen & Elmegreen, 2006; Bournaud et al., 2009). In the third class of models, interactions with satellite galaxies either due to dynamical heating (e.g. Quinn et al., 1993) or accretion of stars (e.g. Abadi et al., 2003) are responsible for the formation scenario of the thick disc. The three families of models suggest distinct outcomes depending on the physical properties of the thick disc. Thus, a proper analysis of their parameters can provide constraints on disc formation models.

Thick discs can give invaluable information about the evolution of galaxy structures. Our goal is to study their structure in unprecedented detail. However, the study of low surface brightness structures is a technical challenge from an observational point of view. In the past several years significant advances have been made in detector technology, observational strategies, and data reduction techniques in order to fulfil the requirements of new ultra-deep imaging science. These improvements allow astrophysicists to reach surface brightnesses fainter than 28 mag arcsec*-2* in the $r$ band. These depths have been achieved on many scales, from large imaging surveys to small galaxy samples, or even in the study of individual objects. These observational depths may be achieved either through very long exposures on small telescopes or using shorter exposures on the largest telescopes. The most remarkable works have achieved a surface brightness of $\mu_{r}\sim 29$ mag arcsec*-2* (3$\sigma$, 10$\times$10 arcsec boxes; e.g. Martínez-Delgado et al., 2010; Ferrarese et al., 2012; Merritt et al., 2014; Duc et al., 2014; Koda et al., 2015; Capaccioli et al., 2015; Peters et al., 2017; Martínez-Lombilla et al., 2019), which is equivalent to a surface brightness $\sim$1500 times fainter than the darkest sky on Earth. Only a few studies have passed beyond the $\mu_{r}=30$ mag arcsec*-2* frontier (3$\sigma$, 10$\times$10 arcsec boxes; e.g. Trujillo & Fliri, 2016; Borlaff et al., 2019). Further progress in the field of low surface brightness science is difficult not because of the collection capability of the telescopes and their detectors, which provide deep imaging and good photon statistics, but rather because of systematic errors. It is therefore crucial to perform very careful data reduction and treatment of both sky and object images (for a more detailed description, see the review by Knapen & Trujillo, 2017). For this purpose, we have used very deep data that reach surface brightness limits of 28.5$-$29 mag arcsec*-2* (3$\sigma$, 10$\times$10 arcsec boxes) from the combined $g$-, $r$-, and $i$ -band images from the IAC Stripe 82 Legacy Project (Fliri & Trujillo, 2016).

The light of a source is affected by the telescope, the instrument, and the detector, as well as by the atmosphere. This results, among other effects, in an increase of the scattered light. The point spread function (PSF) characterises how the light of a source is affected and measures the extent of that scattered light. How to account for the PSF effect on astronomical objects and the importance of this process is well documented (e.g. Michard, 2002; Slater et al., 2009; Sandin, 2015; Trujillo & Fliri, 2016). However, the correction process is still challenging. The shape of the PSF may vary widely depending on the atmospheric conditions and the instrument used, and is generally time-variable on a timescale that can be as short as minutes. The projected surface brightness structure of a source is convolved with the PSF to make up the observed form. If not corrected for, the scattered light therefore adds a systematic component to the observed intensities of our low surface brightness objects. This effect extents to large angular radii, so that although the PSF rapidly becomes faint with increasing radius, the integrated amount of light in its faint extended wings can contribute significantly to observed structures. In the case of extended objects such as galaxies, Sandin (2014, 2015) found that the scattered light of such sources (if not accounted for) can in principle be misinterpreted as a bright stellar halo or a thick disc. Other works have shown that the properties of stellar haloes and the faint outskirts of galactic discs can indeed be significantly influenced by the wings of the stellar PSFs (e.g. Zibetti et al., 2004; de Jong, 2008; Sandin, 2014; Trujillo & Fliri, 2016; Peters et al., 2017; Comerón et al., 2018). In general, the effect becomes stronger the sharper the light distribution of the source (Trujillo et al., 2001a, b). Thus, a proper consideration of the scattered light distribution is crucial, especially in light of upcoming all-sky deep surveys such as the LSST (Large Synoptic Survey Telescope; Gressler et al., 2014).

In this work, we characterise the thick discs in five edge-on galaxies with a wide range of mass. A study of the radial and vertical surface brightness profiles is presented by comparing our data with model predictions. We find that a thick-disc component is required to reach satisfactory fit results, especially in intermediate- to high-mass galaxies, but its brightness is overestimated if the PSF is ignored. Thus, effects due to the PSF are significant when a lower surface brightness is reached, especially in vertical profiles, but this can be accounted for by careful modelling. We show that galaxy outskirts are strongly affected by the faint wings of the PSF, mainly by light from the thin disc that is redistributed by the PSF.

2 IAC Stripe 82 Legacy Project

Our data are taken from the publicly available IAC Stripe 82 Legacy Project111Public data release: http://www.iac.es/proyecto/stripe82 (-50*∘$\leq$ RA $\leq$ +60∘* and -1.25*∘* $\leq$ Dec. $\leq$ +1.25*∘; Fliri & Trujillo, 2016). This project produced optimal co-adds of the SDSS Stripe 82 dataset after reprocessing of all images. The main goal of the IAC Stripe 82 Legacy Project is to maintain the low surface brightness features on all spatial and intensity scales using a non-aggressive sky-subtraction strategy. Thus, the IAC team in particular preserves the characteristics of the background during the reduction, that is, the sky and the diffuse light. The IAC Stripe 82 Legacy Project reduction allows us to reach a surface brightness limit of $\sim$28.5 mag arcsec-2* in the r band (3 $\sigma$, 10$\times$10 arcsec2), which improves by 1.7-2.0 mag upon the single-epoch SDSS images. Thus, it is possible to trace surface brightness profiles down to $\sim$ 30 r-mag arcsec*-2* (Peters et al., 2017), a depth that before has been reached by only a few studies (e.g. Zibetti et al., 2004; Zibetti & Ferguson, 2004; Jablonka et al., 2010; Trujillo & Bakos, 2013; Trujillo & Fliri, 2016). Our images are in the deepest band, a deep r band co-add ($r_{\rm deep}$), which is the linear average of the g, r, and i co-adds.

In order to determine whether the faint wings of the PSF can add artificial light to faint components in the outskirts of galaxies, we need to construct our own PSF to model the galaxies. To estimate this effect in the Stripe 82 data, the IAC Stripe 82 Legacy Project provides ultra-deep and very extended PSFs that extend to 8 arcmin in radius and to over 20 mag in dynamic range for any particular field (details in Fliri & Trujillo, 2016). We use the ultra-deep extended PSFs in the $r_{\rm deep}$ band of our regions of interest (FWHM${}_{r_{\rm deep}}=1.08$ arcsec), obtained by combining the wings of very bright saturated stars from the whole Stripe 82 area (i.e., the outer part of the extended PSF) and PSFEx-generated PSFs with high signal-to-noise ratio (S/N) of a particular field (i.e. the inner part of the extended PSF).

3 Selection criteria and final sample

Thick discs are low surface brightness structures that have been detected by different observational approaches. In particular, their study using integrated light of edge-on galaxies has been performed by Yoachim & Dalcanton (2006) and Comerón et al. (2011, 2012, 2018), for instance. Thus, we need to search for highly inclined disc galaxies with available deep imaging in a multitude of bands and with a spatial resolution that is high enough to resolve these structures. We now describe the selection criteria required to obtain a reliable and unbiased sample of nearby highly inclined galaxies that are imaged in the SDSS g, r, and i bands.

In the first step we used the HyperLeda333http://leda.univ-lyon1.fr/ online catalogue (Makarov et al., 2014) to select all highly inclined galaxies (inclination $\geq 70^{\circ}$ or log${}_{\rm{r25}}\geq 0.477$) that fall within the IAC Stripe 82 Legacy Project area ($\sim$ 275 deg2, $-$50${}^{\circ}\leq$ RA $\leq$ $+$60*∘* and $-$1.25*∘* $\leq$ Dec. $\leq$ $+$1.25*∘*) and have large apparent sizes ($d_{\rm{d25}}\geq$1 arcmin or log${}_{\rm{d25}}\geq 1$) in order to have enough spatial resolution in the discs. At this stage, we had a pre-selection of 56 galaxies. The next step was to make a second-order selection by visual inspection of the SDSS Data Release 12 colour images. Table 1 shows the list of 51 rejected galaxies from that pre-selected sample as well as the primary reason for rejection for all of them.

Morphological features such as warps or possible interacting satellites are undesirable. We also rejected galaxies with very bright foreground objects in the field because they can affect the measurement results. Crowded areas are not recommended, especially if the area of the thick disc is affected. To illustrate this, Fig. 1 shows the bright stars around IC 1756 in detail. It is clear that the extremely crowded region adds undesired light to the low surface brightness galaxy structures.

In order to have discs that are sufficiently visible to derive reliable radial and vertical surface brightness profiles, we rejected galaxies with insufficiently high inclination or non-disc galaxies, as well as galaxies with very strong dust lanes. The latter is particularly noteworthy when the dust lanes are irregular (i.e. not just a straight line) and/or when they are thick and cover most of the visible disc. In the particular case of NGC 1032 (see Fig. 2), the dust lane is extremely thin and lies over the very central part of the disc, so that we are able to study the dust-free regions of the thin- and thick-disc structures. This galaxy is especially interesting as it is the only one in our sample that was studied in previous papers on thick discs using S4G data (Comerón et al., 2011, 2012, 2014, 2018).

As mentioned earlier, we require galaxies of a large apparent size to resolve the disc structure. For the same reason, we need stellar discs with a width that allows us to distinguish thick-disc structures.

As our study is focused on the analysis of disc components, we need to know the distance to the object in order to obtain the physical size and the extent of both the thin and thick disc. Thus, the final selection criterion is the availability of distance measurements in the on-line databases NASA$/$IPAC Extragalactic Database444http://ned.ipac.caltech.edu/ (NED: Helou et al., 1991) and SIMBAD555http://simbad.u-strasbg.fr/simbad/ (Wenger et al., 2000).

The final sample of galaxies consists of three massive galaxies with complex structures (NGC 0429, NGC 1032, and UGC 01040) and two low-mass galaxies with a simple morphology (UGC 01839 and UGC 00931) from a total initial sample of 56 galaxies. Table 2 summarises the general information and physical parameters of the five galaxies in the final sample. The Hubble type was taken from the de Vaucouleurs system from the HyperLeda database (for NGC 0429, UGC 01040 and UGC 00931; Makarov et al., 2014) or from the Spitzer Survey of Stellar Structure in Galaxies (S4G) morphology classification of Buta et al. (2015) in the case of NGC 1032 and UGC 01839. The distance was obtained from the average of the most recent values (since 2010) in NED (Helou et al., 1991). The other properties were all taken from the HyperLeda database (Makarov et al., 2014). When the maximum rotational velocity was not available in Hyperleda, we took the values of the stellar mass from Chang et al. (2015) instead. Fig. 2 shows the SDSS optical colour images of the final galaxy sample, shown overlaid on the $r_{\rm deep}$ IAC Stripe 82 Legacy Project images in order to illustrate the data quality.

4 Method

4.1 Masking process

Ultra-deep imaging requires an accurate background treatment to remove the sky and foreground objects (see, e.g. Trujillo & Fliri, 2016). The goal of the masking process is to avoid any undesirable flux from fore- and background objects over the target in order to only keep the flux from the source we wish to analyse. However, this is difficult because we reach surface brightness levels of $\sim$ 29-30 r-mag arcsec*-2*. At these depths, even the light from bright stars far outside the frame can contribute to the galaxy flux. To address this problem, we used a semi-automatic code developed by Martínez-Lombilla et al. (2019, see details in their Section 3.1). It was optimised for low surface brightness data, allowing us to build a mask of each image covering every emission source except for the galaxy to be analysed.

4.2 Sky subtraction and background treatment

The IAC Stripe 82 Legacy Project dataset has been obtained under generally photometric conditions with good seeing and low sky background. In this work, we use their sky-rectified versions of the co-added images in which the IAC team have used a non-aggressive sky subtraction strategy (details on the reduction process are given in Trujillo & Fliri, 2016).

We performed a second-order (non-aggressive) background subtraction in our region of interest. To do this, we measured the background level following the procedures in Pohlen & Trujillo (2006) and in Martínez-Lombilla et al. (2019, see details in their section 3.2), and we obtained the sky value from our region of interest through very extended radial surface brightness profiles. This permits us to clearly identify where the galaxy surface brightness profile reaches the sky limit (with the current depth of our images), and consequently, where the profile starts to be dominated by the sky background component (see an example in Fig. 3). It is important to note that our sample of galaxies is free of large-scale gradients in the sky background emission so that the sky level does not depend on the location of the outskirts of a galaxy (Peters et al., 2017).

4.3 Profile extraction

A galaxy surface brightness profile in any given direction can be extracted by calculating the flux values through a slit along the corresponding axis (see Fig. 4 for a graphical overview of the following description). As we illustrate in Fig. 4, the slits are double-wedge-shaped and their widths are determined in terms of the physical distance (kpc) in the galaxy, which depends on its apparent size and distance. In our study, the width of the slit increases in the direction perpendicular to the profile as we move away from the centre of the galaxy. The slit is split into logarithmic scale bins that are larger towards the outskirts. The method we followed to obtain the profiles is the same for all SDSS bands.

After the galaxy image was masked (Sect. 4.1) and sky-subtracted (Sect. 4.2), we rotated it such that the plane of the galaxy disc was aligned horizontally. We began profile extraction by obtaining a central radial surface brightness profile along the galactic plane. We calculated the right- and left-side surface brightness profiles separately from the centre of the galaxy to the outer regions. We binned logarithmically in the direction of the profile and perpendicularly to it, which allowed us to increase the S/N in the outskirts of the galaxies (this causes the double-wedge shape of the slits). The final surface brightness profile was obtained by averaging the left and right radial logarithmic bins at the same distance from the centre of the galaxy. Despite small asymmetries in the radial direction, we combined the two sides to explore the surface brightness profiles down to very faint regimes (below $\mu\sim 27.0$ mag arcsec*-2*).

After producing the central radial profile, we obtained a shifted radial surface brightness profile at a given height above and below the galaxy mid-plane. We proceeded similarly to the central profile, but averaged the right and left radial logarithmic bins of the profile above and below the galaxy mid-plane (i.e. the average of four bins for each profile point).

The vertical surface brightness profiles were obtained in the same manner as the radial ones, but we placed a slit along the minor axis of the galaxy. We divided the vertical slit into vertical logarithmic scale bins with a variable radial width (increasing with height above or below the galaxy mid-plane), and then average the vertical logarithmic bins of the upper and the lower profiles at the same distance from the centre of the source. We extracted two vertical surface brightness profiles, one along the minor (vertical) axis of the galaxy, and a second shifted profile at some radial distance from the galaxy centre.

The position of the shifted surface brightness profiles was determined for each galaxy depending on its size. As our main scientific goal is to study the thick discs of our targets, the shifted radial surface brightness profile is located above or below the mid-plane but we tried to avoid (as estimated by eye) the central and inner structures such as the thin disc, the bar, the dust lane, or the bulge. At the same time, any given radial profile cannot be too far above the mid-plane, or we might start to lose S/N. We considered the S/N to be too low when the mean size of the error bars exceeded the intrinsic oscillations of the profile due to the galaxy structure (i.e. spiral arms, regions of star formation, etc). The shifted vertical profile follows a similar procedure to avoid the inner structures as much as possible (i.e. bulge and/or bar). The uncertainties in each bin of the profiles for all cases are defined as in Martínez-Lombilla et al. (2019, in their Section 3.3).

5 Fittings

5.1 Galaxy image model fitting

As shown by Zibetti et al. (2004), Sandin (2014), Trujillo & Fliri (2016), Peters et al. (2017), or Comerón et al. (2018), the properties of the outskirts of a galactic disc can be strongly influenced by scattered light (see Sect. 1). Thus, for a given galaxy, we require a single PSF that characterises the scattered light in the region and is extended enough to cover the whole source (1.5–2 times larger than the galaxy; Sandin 2014). In the case of the IAC Stripe 82 Legacy Project data, we already have an ultra-deep extended PSF for each field ($R\sim 8$ arcmin; see details in Sect. LABEL:chap3/sect2sub1:PSFextendedStripe_82). The next step is to fit the galaxy in order to obtain a model without the PSF effect, and consequently, to study the true structure of the low surface brightness regions.

There are two principal ways of reproducing the PSF effect on the surface brightness distribution of a given object. The first approach is to directly deconvolve the original image by building a very well characterised PSF (e.g. Karabal et al., 2017). The main problem of this method arises in the regions where the noise level of the image is similar to the surface brightness of the PSF, as the results in these areas can be quite uncertain. The second option consists of modelling the intrinsic light of the target and convolving it with the image PSF. The PSF-convolved model of the object is then fitted to the observed data. This process allows an observer to build a new image of the source from the combination (sum) of the PSF-deconvolved model of the target, plus the PSF-convolved fitting residuals (e.g. de Jong, 2008; Trujillo & Fliri, 2016; Peters et al., 2017; Martínez-Lombilla et al., 2019). These residuals include all the non-symmetric features, such as the spiral arms, that the model cannot fit.

Our work is based on the second method listed here, as it requires fewer assumptions about the background sky and is therefore less prone to error in the low surface brightness regime. Consequently, the uncertainty in the determination of the scattered light in the outer regions of the sources is lower than in the first method.

In order to build reliable 2D PSF-deconvolved models of our targets, we used imfit777Precompiled binaries, documentation, and full source code (released under the GNU Public License) are available at the following web site: http://www.mpe.mpg.de/erwin/code/imfit/ . (Erwin, 2015). We modelled the intrinsic light of our observed galaxy sample using a set of 2D image functions that reproduce the profiles of the galaxy components. These functions were fitted to the data by a non-linear minimisation of the total chi-squared statistic, $\chi^{2}$, using a Levenberg-Marquardt optimisation algorithm (Levenberg, 1944; Marquardt, 1963; Moré, 1978). The surface brightness profiles of edge-on galaxies are principally divided into five components: a central bulge, a bar, a thin disc, a thick disc, and (if required) a halo. Depending on the particular case, this could be reduced to only two components (bulge + disc) or even to just a single disc. It is important to point out that imfit provides analytic functions that perform line-of-sight integration through 3D luminosity density models of discs seen at arbitrary inclinations. This is a very useful tool for modelling the disc structure of our edge-on galaxies. In Sect. 6 we explain in detail the particular set of functions we used for each galaxy according to their morphological classification (Table 3).

After performing a reliable and comprehensive 2D fitting for each galaxy and its components, imfit produces the PSF-convolved and -deconvolved models. In order to quantify the effect of the PSF on the surface brightness profiles, we extracted the surface brightness profiles from the deconvolved images following exactly the same slit configurations as in the original data set. The primary benefit is that this process allows us to build a new 2D image of the observed galaxy by the combination (sum) of the PSF-deconvolved model of the target, plus the PSF-convolved fitting residuals (hereafter, deconvolved model). The deconvolved models of the galaxies are used in all subsequent steps of the analysis as the actual image of the galaxy because they are free of any background or scattered-light influence.

Although imfit uses a multi-threaded version of the fftw (fastest Fourier transform in the West888http://www.fftw.org/) library, a relatively speedy function, adding PSF convolution to the image-fitting process does act to slow things down considerably. In addition, the high number of free parameters needed to accurately fit the functions to the observed data (around three or four, depending on the function) also makes model minimisation a time-consuming process. As a result, a large amount of time is needed for a single galaxy model fit. In the case of two or more functions, the fitting to the images was first performed for the innermost components (i.e. bulge and/or bar) and then, with the parameters for these functions fixed, the parameters of the outer components (i.e. thin and/or thick disc) were left free. The outputs for a reasonable fit from these initial runs were used to re-run the code by allowing all (or most of) the model parameters to be free. However, this issue is related to the multiple possible solutions for the modelling process. Consequently, two or more different solutions can constitute a good fit for a given galaxy. To resolve this degeneracy problem, imfit allows bootstrap resampling, which is a method for estimating the variability of our results by making a more detailed analysis of the parameter distributions, including potential correlations between free parameters. The combined set of bootstrapped parameter values can be used to estimate confidence intervals, and consequently, to discern between the multiple good solutions. However, the bootstrapping process is slow because it essentially repeats the fitting process $n$ times (although a somewhat faster version). For a reasonable and reliable set of confidence intervals at least 200 iterations are required, although 1000 or more is preferable. It is therefore better to postpone the bootstrap resampling until the final (and best) fit has been obtained, and then use it to estimate the uniqueness and the uncertainties of the fit.

In the following sections we describe the details of an additional analytical profile-fitting method (Comerón et al., 2011, 2012, 2014, 2018) we used to obtain a more accurate estimate of the physical parameters of the galaxy discs from the deconvolved models. Our main goal is to evaluate the effect of the PSF upon the galaxy structure and to elucidate how the thick discs were formed and how they have evolved.

5.2 Vertical luminosity profiles fitting

After we generated the deconvolved models of each galaxy in the sample, we fitted their vertical luminosity profiles in a similar way as described by Comerón et al. (2011, 2012, 2014). In a later paper (Comerón et al., 2018), the authors explain the overall procedure in enough detail to allow reproduction or further development. We replicated their codes in Python-based routines. The main goal is to decompose a disc into the thin- and thick-disc components in order to obtain their physical parameters from the PSF-free images of the galaxies, that is, from the galaxy luminosity profiles of the deconvolved models. This procedure uses the solutions of the equations of equilibrium for two stellar and one gaseous coupled isothermal discs and therefore has the advantage of being physically motivated. There are two main reasons to make this particular vertical profile fit instead of just taking a surface brightness profile of the fitted image functions that are returned by imfit. The first is that the method of Comerón et al. of fitting vertical luminosity profiles provides an independent solution for the thin- and thick-disc decomposition and their parameters. In the case of similar solutions, it can therefore give a measure of robustness to the results. The second reason is that the physical parameters can be more accurately determined because the decomposition takes into account the gravitational interactions between the thin- and thick-disc components.

In essence, we derived the vertical height of each disc from the luminosity profiles in the vertical direction (i.e. perpendicular to the galaxy mid-planes). To do this, we fitted the vertical profiles of the data (in this case, the deconvolved model of each galaxy provided by imfit) by comparing them with a grid of models that are the result of a system of equations that describe three gravitationally coupled vertically isothermal discs (thin, thick, and gas discs). We thus assumed that the stellar discs are composed of two vertically isothermal discs (at a certain distance from the galaxy centre) in hydrostatic equilibrium. Both discs can interact gravitationally and have a unique vertical velocity dispersion at all heights. The third is a non-stellar (gas) disc with a given mass but no optical luminosity. The system of equations that describes the vertical density distribution for a point at a fixed height $z$ above and below the galaxy mid-plane can be written as follows (adopting the formalism from Narayan & Jog, 2002):

[TABLE]

where $\rho_{i}$ denotes the density volume of either the thin (subscript $\rm t$), thick ($\rm T$), or gas ($\rm g$) disc; $\sigma$ is the velocity dispersion; and $K_{{\rm DM}}$ is a term related to the behaviour of the dark matter halo. However, we did not take into account the last term ($dK_{{\rm DM}}/dz=0$) in our fitting as it only introduces a small bias (in the worst case, $M_{\rm T}/M_{\rm t}$ could be overestimated by $\sim$10 % Comerón et al., 2011, 2012). This allowed us to assume cylindrical symmetry for the galaxy, which considerably simplifies the numerical treatment. As we are mainly interested in the study of the luminosity profiles, the gas component of the disc was not considered in this work in order to simplify the process by reducing the number of free parameters, so $\rho_{rmg}=0$. The galaxy disc components were then fitted by considering that the thin and thick discs are two stellar fluids in hydrostatic equilibrium. However, a downside of this method is that the fitted function is not analytical and needs to be solved by numerical integration.

The luminosity profiles have a certain width given by the radial bins of $0.2\,r_{25}<|R|<0.5\,r_{25}$ and $0.5\,r_{25}<|R|<0.8\,r_{25}$ on each side of the bulge, where $r_{25}$ is the radius of the 25 mag arcsec*-2* level in the $B$ band, taken from the Hyperleda database (Makarov et al., 2014). To determine where to stop fitting in height above and below the galaxy mid-planes, we established a 20 % threshold for the maximum contribution of the error of the bin (${\rm error_{bin}}$) to its total amount of light (${\rm\#counts_{bin}}$):

[TABLE]

This method for cutting the fit is different from that adopted in Comerón et al. (2011, 2012, 2014, 2018). We used robust means to obtain the mean value of each bin, therefore we minimised the noise in the fainter regions of the surface brightness profiles. Furthermore, we solved the problem with “bright haloes” because the functional shape of the disc components drops quickly and the very outer regions with excess light were not considered. It is also important to mention that the fit was made over the profile in magnitudes instead of number of counts in order to properly weight differing galaxy structures. The primary difference between our version and the code of Comerón et al. is that we do not need to deconvolve the observed data with the PSF of the images because we have already taken this into account (see details in Sect. 5.1). Finally, we note that if the scale heights of the thin and the thick discs remain roughly constant and the scale lengths of both discs are similar, the line-of-sight integration is not a problem when Eq. 1 is applied.

5.2.1 Mass-to-light ratio

The decomposition of galaxy discs into two gravitationally coupled vertically isothermal discs by solving the system of equations above (Eq. 1) is based on volume densities as a function of height. However, in our data, we measured surface brightness rather than volume density. We required the mass-to-light ratio ($\Upsilon_{i}=M_{i}/L_{i}$) of the components of our targets, that is, the thin and thick discs, to convert densities into light. In addition, structural decomposition provides the ratio of disc masses ($M_{\rm T}/M_{\rm t}$). It is therefore important to be cautious with the assumed value of mass-to-light ratios of these components, $\Upsilon_{\rm T}/\Upsilon_{\rm t}$, because this quantity contributes most to the uncertainty in the fit.

Comerón et al. (2011, 2012, 2014) discussed the advantages and disadvantages of differing star formation histories (SFHs) with the purpose of obtaining the mass-to-light ratio of a given galaxy using SFHs for the thin and thick disc of the Milky Way. They used a fixed SFH with a constant $\Upsilon_{\rm T}/\Upsilon_{\rm t}$ for all the galaxies, independent of their individual properties. To avoid making this assumption, we used the relation between the photometric broad-band colour and the mass-to-light ratio from Bell & de Jong (2001),

[TABLE]

where $a_{\lambda}$ and $b_{\lambda}$ are coefficients taken from Roediger & Courteau (2015), which provide linear least-squares fits in the SDSS bands to the mass-to-light versus colour relations (MLCRs) created on the basis of the spectral energy distributions by Bruzual & Charlot (2003). Equation 3 allows us to obtain the mass-to-light ratio for each disc component $M_{i}/L_{i}$ by taking radial surface brightness profiles along the thin and the thick disc separately and obtaining the corresponding colour profiles.

5.2.2 Physical parameters of thin and thick discs

From the vertical fitting (Sect. 5.2) of the luminosity profiles of the deconvolved model of each galaxy, we estimated a set of physical parameters that allowed us to analyse the thin- and thick-disc properties. These quantities help to discern between the various disc formation scenarios and are described as follows:

• •

Scale heights: the scale height of the thin- and thick-disc components ($z_{\rm 0t}$ and $z_{\rm 0T}$, respectively) is defined as the vertical distance (in arcsec) over which the density of a particular physical variable (the luminosity obtained from the profile in this case) drops by a factor of $e$ ($\sim$2.718) from its centre and it is heavily dependent on the type of object. For each bin the scale height was determined by establishing the distance between the location where the flux is a factor 100 lower than in the galaxy mid-plane and the position where the flux is 100$\times e$ lower than at the mid-plane in the fitted surface brightness profiles for each of the two stellar disc components.

• •

Mid-plane mass-density ratio: the ratio of the thick-disc over the thin-disc mass density in their central regions ($\rho_{\rm 0T}/\rho_{\rm 0t}$). This parameter is directly obtained from the fit, as the solution of Eq. 1 yields the corresponding values of $\rho_{\rm 0T}$ and $\rho_{\rm 0t}$ for each galaxy.

• •

Thick- to thin-disc mass ratio: the ratio of the cumulative mass density (column mass) along the profile of each disc component, that is, a surface-density ratio ($\Sigma\rho_{\rm T}/\Sigma\rho_{\rm t}=\Sigma_{\rm T}/\Sigma_{\rm t}$).

6 Results

In this section, we present the main results of our fits for each galaxy individually in the $r_{\rm deep}$ band of the IAC Stripe 82 Legacy Project (i.e. the SDSS g, r, and i combined band). The aim is to compare the observed data and the corresponding PSF-deconvolved model to evaluate the PSF influence using surface brightness profiles in both the radial and vertical directions (as explained in Sect. 4.3). It is crucial to perform the PSF-deconvolution in order to work with what is assumed to be the actual flux from these galaxies. We therefore analysed the radial and vertical surface brightness profiles for the observed data, the PSF-convolved model, and the PSF-deconvolved model. For each model, the corresponding galaxy components were also studied. We performed a bootstrap resampling (described in Appendix A) for each 2D galaxy model using the astronomical software imfit (Erwin, 2015). The main reason for a bootstrap resampling is that there is a strong degeneracy for the disc parameters with one another, so that multiple possible solutions may appear for a given galaxy in the modelling process.

After we obtained the PSF-deconvolved image (plus the residuals of the PSF-deconvolved fit), we fitted its vertical light distribution by solving a system of equations that describe two gravitationally coupled vertically isothermal discs (thin and thick discs), as described in Comerón et al. (2011, 2012, 2014, 2018), and in Sect.  5.1. The main goal was to decompose the disc of the PSF-deconvolved model into the thin- and thick-disc components considering their interactions, in order to obtain their main physical parameters such as the scale heights ($z_{{\rm t}}$ and $z_{{\rm T}}$, respectively), the initial mass density ratio ($\rho_{\rm 0T}/\rho_{\rm 0t}$), or the thick- to thin-disc mass ratio ($\Sigma_{{\rm T}}/\Sigma_{{\rm t}}$). The weighted mean values of all these physical parameters of the thin and thick discs of each galaxy are shown in Table 3. All fitted surface brightness profiles for the whole galaxy sample are shown in Martínez-Lombilla (2018).

6.1 Intermediate- to high-mass galaxies: NGC 0429 and UGC 01040

NGC 0429 and UGC 01040 (see optical colour images from the SDSS DR12 in Fig. 2) are lenticular galaxies of intermediate to high stellar mass (4.4-8.3$\times$1010 $M_{\odot}$ and 1.3-2.8$\times$1010 $M_{\odot}$, respectively). NGC 0429 is the most distant galaxy in our sample at more than twice the distance to UGC 01040, the next galaxy (see Table 2). Table 2 also shows that it is unclear whether UGC 01040 has a bar, but our models return very good fits when a bar is included.

In Fig. 5 we show an example of the radial and vertical surface brightness profiles for the observed data of UGC 01040, the PSF-convolved (and -deconvolved) model, and their components. In this particular case, it was necessary to use four galaxy components to find a reliable fit: a bulge using a 2D Sérsic function, a bar represented by an elliptical 2D Sérsic function using generalised ellipses (“box” to “disc” shapes) for the isophotes instead of pure ellipses (Athanassoula et al., 1990), and finally, a thin and a thick disc using a 3D broken exponential disc function. For NGC  0429, we required three components: a bulge based on a 2D Sérsic function, and both a thin and a thick disc, but this time using a 3D pure exponential disc function (for more details about the analytical functions, see Erwin, 2015). This type of 2D model fit was performed for all the galaxies in the sample.

In the case of NGC 0429, the PSF effect in the radial surface brightness profiles becomes relevant below $\sim$28 mag arcsec*-2*. Thus, if we do not account for the PSF along the radial profile, we overestimate the outer part of the galaxy but by a similar amount in all the components. For the vertical profiles, this effect is much more dramatic as it starts to be significant beyond $\sim$25 mag arcsec*-2*. Consequently, the galaxy flux and thus its mass (for a given $M/L$) is overestimated by a factor of $\sim$10 in the outskirts ($z\sim$16 kpc) because of the way in which the PSF affects the bulge and the thin disc, rather than the thick disc. Consequently, we infer that the PSF-deconvolved components of the galaxy model contribute less to the total flux when the slope of the profile is steeper (although this also depends on the central intensity).

For UGC 01040, the most complex galaxy of the sample with four structural components for the 2D model fit, the profiles do not show any remarkable discrepancy between the observed data and the PSF-convolved model (see the residual plots in the lower part of each panel in Fig. 5), except for a slight excess of light in the data in comparison with the PSF-convolved model at $z\sim$2-4 kpc in the shifted vertical surface brightness profile, probably due to some asymmetry. The PSF effect in the radial direction (top panel Fig. 5) becomes relevant below $\sim$28 mag arcsec*-2*, producing a slight overestimate in the outer part of the galaxy if not accounted for. Once again, the PSF has a dramatic effect on the vertical light distribution in the galaxy beyond $\sim$25 mag arcsec*-2* (bottom panel Fig. 5). The galaxy UGC 01040 is most affected by the PSF in the sample. The main reason is the very steep slope of the vertical surface brightness profiles. However, in this case, the dominant component in the outer regions of the PSF-convolved model is the thick disc, although the thin disc and the bar also contribute significantly, mainly in the case of the minor axis profile. After the PSF is removed, the vertical light distribution reaches a level of $\sim$40 % of the observed data, which corresponds to a decrease in flux (i.e. mass for a given $M/L$ value) of a factor of at least 40 at $z\sim$4 kpc. Thus, it is critical to account for the PSF effect for these intermediate- to high-mass galaxies as the scattered light produces an overestimation in all galaxy components.

We proceed to analyse the results of the vertical fits to the light distribution through the resolution of a system of equations that describe gravitationally coupled vertically isothermal discs (thin and thick discs), as described in Comerón et al. (2011, 2012, 2014, 2018, and Sect. 5.2 of this work). Figure 6 shows the thin- and thick-disc decomposition for the fit of the PSF-deconvolved model at a given radial distance of UGC 01040. In comparison with previous 2D models, these fits just reach surface brightness values above $\sim$27-27.5 mag arcsec*-2*, depending on the bin, and seem to have problems when the profile is more noisy or when it presents asymmetries. These fits were not considered in our results.

We note that after we added the residuals of the PSF-convolved model to the PSF-deconvolved galaxy image (i.e. green point profiles in Fig. 6), the dramatic effect of the PSF diminished along the vertical surface brightness distribution in comparison with the bottom panel of Fig. 5, but it still adds artificial light to the galaxy outskirts. The difference between the observed data and the PSF-deconvolved model is now $\sim$1 mag arcsec*-2* at $z\gtrsim$10 kpc for NGC 0429, and $\sim$1.5 mag arcsec*-2* at $z\sim$4 kpc for UGC 01040, which corresponds to an excess of galaxy flux, and therefore mass, of a factor of $\sim$2.5 and $\sim$4, respectively. This is a critical result that reveals just how crucial careful PSF characterisation is.

For each galaxy we obtained a set of physical parameters that are shown in Table 3. We did not derive uncertainties because they are very difficult to obtain due to the complicated fitting process (we solve a system of equations that describe two gravitationally coupled discs). The reliability of the fits is well studied in Comerón et al. (2011), therefore we use the same number of significant digits. In our case, the deviations between the different values of the physical parameters in the panels of each of those figures may be representative of the uncertainties in the method.

It is important to note that in the case of NGC 0429, against the expectations, the vertical scale height of the thick disc ($z_{{\rm T}}$) is larger after the PSF is accounted for, so that the thick disc becomes thicker. The ratio between the observed data and the PSF-deconvolved thick- to thin-disc mass ratios (i.e. $(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm obs}/(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm PSFdec}$) is equal to 2.3, which is in relatively good agreement with the previously estimated excess of mass (a factor of $\sim$2.5). In UGC 01040, the vertical scale of the thick disc ($z_{{\rm T}}$) is lower when the PSF is accounted for, so the vertical distribution of the thick disc becomes steeper, in full agreement with expectations. The ratio between the observed data and the PSF-deconvolved thick- to thin-disc mass ratios (i.e. $(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm obs}/(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm PSFdec}$) is equal to 2.3, which is lower than the previously estimated excess of mass (a factor of $\sim$4). In any case, both values are high enough to claim that a proper PSF treatment of the observed data is absolutely indispensable. The thin- and thick-disc components may have different shapes in the vertical luminosity profile fits (e.g. Fig. 6) than in the 2D models (e.g. bottom panel in Fig. 5) because in the first physical interactions between the thin and thick discs are considered.

6.2 Low-mass galaxies: UGC 00931 and UGC 01839

UGC 00931 (optical colour image from the SDSS DR12 in Fig. 2) is the galaxy with the lowest inclination in the sample and also the least massive (3.1$\times$109 $M_{\odot}$), while UGC 01839 (colour image in Fig. 2) is a low-mass (5.9$\times$109 $M_{\odot}$) late-type spiral and the closest galaxy in the sample (see details in Table 2). Visually, the disc of UGC 00931 is very diffuse and asymmetrical, which leads to enhanced residuals in the 2D model. UGC 01839 also has an asymmetric morphology with diffuse and broken arms, and no central bulge. However, as indicated in Table 2 with the letters “sp” after the morphological type, its classification has considerable uncertainty because its internal structures are not easily distinguished due to the high inclination of the galaxy.

We only required a single 3D broken exponential disc for a satisfactory fit of these two low-mass edge-on galaxies. Thus, no thick (or thin) disc component was necessary to model them. To illustrate this, the vertical surface brightness profile of UGC 01839 and its 2D model are shown in Fig. 7. This vertical surface brightness profile is affected by the PSF when values below $\sim$26 mag arcsec*-2* ($z\sim$2 kpc) are reached. At $z\sim$3 kpc, the mass is overestimated by a factor of $\sim$3. UGC 00931 is affected in an equivalent way.

The following results are based on fits along the vertical light distribution, considering the physical interactions between two gravitationally coupled vertically isothermal discs (details in Comerón et al., 2011, 2012, 2014, 2018, and in Sect. 5.2 of this work). In Fig. 6 we illustrate the thin- and thick-disc decomposition for the PSF-deconvolved model of UGC 01839, down to a surface brightness of $\sim$27.5-28 mag arcsec*-2*.

By comparing these results with the previously presented 2D models (e.g. Fig. 7), a discrepancy becomes clear: in this case, the fit can be made using two discs instead of just one. Comerón et al. (2018) found only 17 galaxies out of 141 that showed no clear traces of a thick disc. These 17 galaxies are on average less massive than those that show at least two discs. These fits based on gravitational interactions are limited to a two-disc decomposition, so that we always obtained two components. Thus, “no clear traces of a thick disc” means that the result is similar to ours in the performed fits for these two low-mass galaxies, in the sense that either one disc clearly dominates the other (e.g. Fig. 6), or the two discs are very similar and therefore confirm a single disc-like profile. The vertical scale height ratios of the thin and thick fitted discs are quite similar in both cases (in comparison with other galaxies in the sample, less than a factor of 2). In UGC 00931, the disc mass excess in the galaxy outskirts due to the PSF effect (taking into account all the non-symmetric features such as the spiral arms, which the model cannot fit, and the physical interactions between both discs) is now a factor of $\sim$1.5, while for UGC 01839, the disc mass is overestimated a factor of $\sim$2. Thus, as for the two intermediate- to high-mass galaxies, the low surface brightness structures of UGC 00931 and UGC 01839 are significantly overestimated when the PSF is not considered because of the light that is scattered from the inner to the outer regions.

The weighted mean values of the physical parameters obtained for our two low-mass galaxies are shown in Table 3. The scale heights of both the thin and thick disc ($z_{{\rm T}}$ and $z_{{\rm t}}$) are smaller when the PSF is accounted for, as expected. This is because after removing the PSF effect, the fit uses a dominant component (the thick disc) that fits the PSF-deconvolved model points, therefore its scale height should be smaller, plus a faint thin disc. These results are consistent with an increase in the ratio of thick- to thin-disc masses.

6.3 A particular case: NGC 1032

NGC 1032 (see an optical colour image from the SDSS DR12 in Fig. 2) is a very interesting galaxy because it is the only one in common with the S4G-based sample of Comerón et al. (2011, 2012, 2014, 2018). It has a high dynamical mass (3.3$\times$1011 $M_{\odot}$, see Table 2), which is isolated (Tully, 2015), and it has a straight and thin dust lane along its mid-plane, out to $\sim$7.5 kpc.

The radial and vertical surface brightness profiles of NGC 1032, as for NGC 0429, required three components: a bulge based on a 2D Sérsic function, and both a thin and a thick disc, but this time using a 3D pure exponential disc function (see Erwin, 2015). The most interesting aspect in this galaxy is that the PSF does not play an important role, not even in the outermost parts of the galaxy. This is mainly because the slope of the profiles is smoother and less steep than in the other galaxies of the sample. By analysing the structural components individually, it is clear that the PSF effect is stronger in the bulge and thin-disc components than in the thick disc, and along the vertical direction, although the thick disc remains the dominant component.

In Fig. 6 we show the thin- and thick-disc decomposition for the fit of the PSF-deconvolved model. The surface brightness of the fit reaches deeper values than for the rest of the sample ($\gtrsim$28 mag arcsec*-2*) because the profiles are less noisy and better sampled (the spatial resolution of NGC 1032 is better because of its larger apparent size). Fig. 6 shows that the PSF effect is almost negligible, in agreement with the 2D fitting, because the PSF-deconvolved model is systematically immediately below the observed data and has the same shape.

The weighted mean values for the physical parameters are shown in Table 3. In this case, the scale heights of the thick discs are the same (within errors) in the observed data as in the PSF-deconvolved model. Thus, the mass ratios are also very similar (the small difference comes from the change in the scale height of the thin disc).

7 Discussion

We have presented a new method as applied to a fairly diverse sample of galaxies (from early to late type ) to reach two main goals: (1) properly account for the PSF effect in different types of galaxies, and (2) parametrise the galaxy disc components by considering their physical interactions. To do this, we combined two techniques that have already been partly used in previous works (see details and references in Sect. 5). Our resulting method allowed us to reach a surface brightness limit of $\mu_{r_{\rm deep}}\sim$28.5$-$29 mag arcsec*-2* (3 $\sigma$, 10$\times$10 arcsec2). In this section we discuss the advantages and disadvantages of the different methods when the PSF effect in our sample of galaxies is analysed. We also consider possible formation and evolution scenarios of the thick disc. Throughout this discussion we use the weighted mean values of all the physical parameters of the thin and thick discs of each galaxy, as in Table 3.

7.1 PSF effect

Sandin (2014, 2015) found that the scattered light of extended sources (if not accounted for) might be misinterpreted as a bright stellar halo or a thick disc. Other works have shown that the properties of stellar haloes and the faint outskirts of a galactic discs can indeed be influenced significantly by the wings of the stellar PSFs (e.g. Zibetti et al., 2004; de Jong, 2008; Sandin, 2014; Trujillo & Fliri, 2016; Peters et al., 2017; Comerón et al., 2018). In general, the effect becomes stronger when the light distribution of the source is sharper (Trujillo et al., 2001a, b). Here, we have used a very extended PSF for each galaxy (Sect. LABEL:chap3/sect2sub1:PSFextendedStripe_82) in order to properly account for time, angle, and field dependencies, and with a diameter of at least twice that of the source (Sandin, 2014).

We used imfit to perform reliable 2D fitting of each galaxy and its components, producing PSF-convolved and -deconvolved models. From the surface brightness profiles in both radial and vertical directions of all the data and the models, we can quantify the PSF effects. We find that the observed radial surface brightness profile deviates less from the PSF-deconvolved profile than the vertical one, as expected. UGC 01040 is a clear example of a galaxy with a very steep distribution of light, while NGC 1032, in contrast, has very smooth surface brightness profiles along both axes.

In the radial direction, the PSF influence becomes significant at $\mu_{r_{\rm deep}}\sim$ 28 mag arsec*-2* (i.e. the radial surface brightness profile of the observed data starts to deviate from the PSF-deconvolved profile). At fainter levels, the scattered light contribution is up to $\sim$1.5$-$2.5 times that of the intrinsic stellar disc brightness, which is slightly lower than the values reported by Trujillo & Fliri (2016, they found mass excesses up to three times higher, probably because they went much deeper). However, it is important to note that these values were extracted from the PSF-deconvolved models without considering the residuals (i.e. all non-symmetric features such as the spiral arms that the model cannot fit) and the model was obtained from an analytical function rather than by taking the physical interactions of the discs into account.

In the vertical direction, we achieved satisfactory fits to the PSF-deconvolved models, but in this case, considering the residuals and the physical interactions of the discs (gravitationally coupled vertically isothermal thin and thick discs in hydrostatic equilibrium). In our galaxy sample, the PSF effect starts to be relevant at $\mu_{r_{\rm deep}}\sim$ 27 mag arsec*-2* for the less massive galaxies, UGC 00931 and UGC 01839, but at $\mu_{r_{\rm deep}}\sim$ 25 mag arsec*-2* for the galaxies of intermediate to high mass, NGC 0429 and UGC 01040 (in the high-mass galaxy NGC 1032, the PSF is negligible because its light distribution is very smooth). At lower values, the PSF (if not accounted for) adds artificial light to the outskirts of the galaxies, which corresponds to an excess in galaxy flux, and therefore in mass (for a given $M/L$), of a factor of $\sim$1.5$-$2 times for our two low-mass sources, and $\sim$2.5$-$4 times for our two galaxies of intermediate to high mass. This is a critical result that reveals that a careful PSF treatment of observed data is absolutely indispensable for deep imaging of extended objects.

We also analysed how the wings of the PSF impact each galaxy component. The PSF mainly spreads the light from the brightest regions (usually located in the centre) to the faintest (the outer parts). Thus, a galaxy component is more affected by the PSF wings if the light distribution is sharp. In particular, the intrinsic stellar disc brightness of NGC 0429 is mainly composed of the PSF wings deriving from its prominent bulge and also of the thin disc, rather than the thick disc. If not corrected for, this produces an incorrect quantification of the physical properties of the galaxy structures. In the case of NGC 1032, the PSF-influenced thick disc is the dominant component, although the thin disc and the bulge contribute more when high altitudes above or below the galaxy mid-plane are reached, beyond $z\sim$35 kpc (see all the plots and more details of NGC 0429 and NGC 1032 in Martínez-Lombilla, 2018). Finally, the low surface brightness regions of UGC 01040 (top panel in Fig. 5) are dominated by the effect of the PSF wings in the thick disc, closely followed by those in the bar and the thick disc. This is important because properly accounting for the PSF is not only necessary for a better estimation of the actual sizes and masses of the galaxies, but also for an accurate determination of the shapes, sizes, and masses of their structural components.

7.2 Possible contribution of a halo component

Our work reaches surface brightnesses of $\mu_{r_{\rm deep}}\sim$28.5$-$29 mag arcsec*-2* (3 $\sigma$, 10$\times$10 arcsec2) such that the stellar halo may start to become visible (Martín-Navarro et al., 2014, and references therein). Thus, it is worthwhile to study the cases of NGC 0429 and NGC 1032 because in both galaxies the fit reaches regions at $R>$30 kpc. In NGC 0429, the presence of a stellar halo is not clear, as the model fits all the data points up to the very outskirts of the galaxy. However, in NGC 1032, it is not possible to obtain a good fit when all the data points in the outer parts are included. Perhaps a stellar halo component plus a thick disc with a shorter scale height (and probably also a shorter scale length) might constitute a more realistic model of the galaxy. The physical parameters of a stellar halo component are different from those of the thick disc, but it is not practical to account for a halo in this work. A logical next step for future improvement of our method therefore is to add this component into our fits. This implies accounting for physical interactions between the halo and the disc components, so that more computer time and capability will be needed as the number of free parameters will increase. Ideally, deeper imaging should also be used.

7.3 Origin of thick discs

Most if not all disc galaxies have a thin (classical) disc and a thick disc. In most models, thick discs are thought to be a necessary consequence of the disc formation process and/or of the evolution of the galaxy. However, because of their faint surface brightness ($\mu_{rmg}\sim$26.0 mag arcsec*-2*), the detailed study of the physical properties of thick discs has been very difficult. We here studied their structure at significant depth by reaching a surface brightness limit of 28.5$-$29 mag arcsec*-2* with the combined $g$, $r,$ and $i$ band images from The IAC Stripe82 Legacy Project.

Yoachim & Dalcanton (2006) and Comerón et al. (2011, 2012, 2018) established that thick discs are more dominant in low-mass galaxies than they are in massive galaxies. The ratio between the thick- and thin-disc masses is typically in the range 0.2 $<M_{\rm T}/M_{\rm t}<$ 1 for galaxies with $v_{c}>$ 120 km/s, but low-mass galaxies can have much higher values, up to $M_{\rm T}/M_{\rm t}\sim$3. In Fig. 8 we show the thick- to thin-disc mass ratio (i.e. $\Sigma_{{\rm T}}/\Sigma_{{\rm t}}$) that results from our PSF-deconvolved model fits, plotted as a function of the circular velocity ($v_{c}$) for the galaxies in our sample with circular (i.e. rotational) velocity information (see Table 2). Our results are in good agreement with those of Comerón et al. (2018, the overploted filled dots in Fig. 8), and we confirm that thick discs are comparable in mass to their thin-disc counterparts.

It is common to think of a thin disc as being dominant, and embedded in a diffuse and extended thick-disc component. The low-mass galaxy characteristics point to a thick disc as the dominant component with a thinner component inside, however, which is compatible with Fig. 8, which shows a clear tendency towards higher $M_{\rm T}/M_{\rm t}$ when $v_{c}$ decreases (i.e. the mass of the galaxy). This in turn agrees very well with our results for the low-mass galaxies, UGC 00931 and UGC 01839, where the 2D fits suggest a single disc component (see also detailed plots in Chapter 3 of Martínez-Lombilla, 2018). In the fits to the PSF-deconvolved models of these galaxies (e.g. Fig.  6), we found two discs in both cases. However, in UGC 00931 the discs are quite similar, consistent with a single disc, whereas for UGC 01839 the thick disc clearly dominates the thin disc. This is a natural consequence in a model where thick discs were first formed from a turbulent gas disc with intense star formation (Elmegreen & Elmegreen, 2006; Bournaud et al., 2009; Comerón et al., 2014, 2018). This corresponds to a Universe that evolves in line with the philosophy of down-sizing, where low-mass galaxies take longer to evolve (Cowie et al., 1996). In this picture, thin discs in low-mass galaxies are still in their infancy, whereas those in high-mass galaxies are already mature. An alternative interpretation may be that the star formation in the thin discs of low-mass galaxies is not efficient enough, and the thin disc cannot be distinguished from the thick disc (or could even be considered as the same structure). Thus, thick (or thin) discs might not be ubiquitous for all types of disc galaxies (contrary to the result of Comerón et al., 2018).

Figure 9 shows the scale heights of the thin ($z_{{\rm t}}$) and the thick ($z_{{\rm T}}$) discs as a function of the circular velocity ($v_{c}$) for the three galaxies in the sample with available circular velocity values (equivalent to Figure 17 from Comerón et al., 2018). Compared to the results of Comerón et al. (2018), UGC 00931 and UGC 01839 fall in the expected area when we consider the uncertainties in the distance measurements and all the uncertainties associated with our rather complicated fits. The thin- and thick-disc scale heights of these galaxies correlate with the circular velocities of the galaxies, in good agreement with Comerón et al. (2018). In particular, the positions of UGC 00931 and UGC 01839 in Fig. 9 seem to fit the gas-poor galaxy sample better.

NGC 1032 is the only galaxy in common with the sample of Comerón et al. (2018), so that here we show an optical (and deeper) counterpart to their work. When we return to Fig. 9, the data point for NGC 1032 is even more of an outlier than in Comerón et al. (2018) for the thin- and thick-disc scale heights. This galaxy is a gas-poor object and has one of the two largest scale heights for a thin disc in the results of Comerón et al. (2018). However, our results suggest a higher discrepancy in the scale heights of both discs ($z_{{\rm t,C18}}=0.8$ kpc versus $z_{{\rm t,TW}}=$1.8 kpc, and $z_{{\rm T,C18}}=$3.3 kpc versus $z_{{\rm t,TW}}=$8.8 kpc; where we used C18 for Comerón et al. 2018, and TW for results of PSF-deconvolved models in this work), even though the scatter in the relation is larger for gas-poor galaxies than for gas-rich ones. As NGC 1032 is isolated (Tully, 2015), the disc heating may have been caused by the merger with a companion.

The scale heights of the thin and thick discs of NGC 1032 are very different from those reported by Comerón et al. (2018), which might be due to various reasons. The main differences between this work and that of Comerón et al. (2018) are the different wavelength ranges, the PSF treatment, and the depth of the data. The latter is a key factor as it may allow us to see the stellar halo, so that our thin- and thick-disc scale heights could be larger that they otherwise may be. In the case of NGC 1032, we already saw that the PSF does not cause significant effects because the light distribution of the galaxy is very smooth. On the other hand, each wavelength gives different information about the underlying stellar population properties. Because near-infrared light traces older populations than the optical, our larger scale heights could point to a young stellar population component in the outer parts of both discs. In agreement with this, the thick- to thin-disc scale height ratios of this galaxy are quite similar in both works: $(z_{{\rm T}}/z_{{\rm t}})_{\rm C18}=$4.1, and $(z_{{\rm T}}/z_{{\rm t}})_{\rm TW}=$4.9. Another fact in favour of this hypothesis is that the thick- to thin-disc mass ratios are also similar: $(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm C18}=$0.4, and $(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm TW}=$0.5. When we take the thick- to thin-disc mass ratio of our observed data ($(\Sigma_{{\rm T}}/\Sigma_{{\rm t}})_{\rm TWobs}=$0.4), the result is even more similar, which may suggest that the scattered light does not affect the ratio much.

Comparing all our fits with those in Comerón et al. (2018), we reach deeper surface brightness levels, which highlights the PSF effect over galaxy structures beyond $\mu_{r_{\rm deep}}\sim$27 mag arcsec*-2*. While Comerón et al. (2018) did not find a significant excess of mass due to scattered light, we obtain overestimates in the mass of the thick disc of up to a factor of four in the PSF-deconvolved model than in the data at $\sim$28 mag arcsec*-2*. However, Comerón et al. (2018) found that even though the thin disc is less extended vertically than the thick disc, it dominates the emission at large heights because of the light scattered from the mid-plane to large heights, in agreement with our results. This effect would remain unnoticed if we had not accounted for the extended wings of the PSF.

We finish by analysing the mass-to-light ratios between the thin and thick discs ($(M/L)_{\rm T}/(M/L)_{\rm t}=\Upsilon_{\rm T}/\Upsilon_{\rm t}$, obtained as described in Sect. 5.2.1), which gives an idea of the stellar population age of each component. For the 2D PSF-deconvolved model of our sample of galaxies, we obtain the following values: 0.99, 1.07, 1.97, 1.46, and 0.56 (for NGC 0429, UGC 01040, UGC 00931, UGC 01839, and NGC 1032, respectively, with $\sim$15 % uncertainty). In the cases of NGC 0429 and UGC 01040, $(M/L)_{\rm T}/(M/L)_{\rm t}$ is almost unity, so that both disc components seem to have similar stellar populations, suggesting a common origin for the thin- and thick-disc components. NGC 1032 shows younger stellar populations in the thick disc than those in the thin disc ($(M/L)_{\rm T}/(M/L)_{\rm t}=$0.56). In the low-mass galaxies UGC 00931 and UGC 01839, the thick disc looks older than the thin disc, in agreement with the hypothesis that their thick discs formed first from a turbulent gas disc with intense star formation (Elmegreen & Elmegreen, 2006; Bournaud et al., 2009; Comerón et al., 2014, 2018).

In general, our results are in good agreement with those in Comerón et al. (2018), although there are some discrepancies because we analysed data in a different wavelength and down to a deeper surface brightness limit. This suggests that our method is able to properly reproduce the disc components beyond $\mu_{r_{\rm deep}}\sim$27 mag arcsec*-2*, bringing light to the formation and evolution scenario of low surface brightness galaxy structures.

Although our sample of galaxies is small, we find that the three formation scenarios proposed to explain the origin of the thick discs are likely. It seems that the formation of these structures depends on the galaxy mass and morphology, and on the mass assembly history of a given galaxy. In low-mass galaxies, UGC 00931 and UGC 01939, only a single disc component is required, in agreement with a model where thick discs were first formed from a turbulent gas disc with intense star formation. In the intermediate-high mass galaxies NGC 0429 and UGC 01040, a common origin is suggested for both the thin- and thick-disc components. Finally, NGC 1032 seems to have formed a thick disc as a result of merger events.

8 Conclusions

We characterised the thick discs of five edge-on galaxies down to a surface brightness limit of $\mu_{r{\rm deep}}\sim 28.5-29$ mag arcsec*-2* (3 $\sigma$, 10$\times$10 arcsec2) with combined $g$-, $r$-, and $i$ -band images from the IAC Stripe 82 Legacy Project. The method accounts for a careful analysis of the background, masking, and a meticulous PSF treatment through 2D galaxy modelling. Our results are based on a study of radial and vertical surface brightness profiles, and we compared our data with models. The galaxy disc components were fitted by considering that the thin and thick discs are two stellar fluids in hydrostatic equilibrium. The main conclusions of this work are listed below.

We developed and tested a new method and applied it to a fairly diverse sample of five galaxies, from early- (i.e. S0) to late-type galaxies (Sdm). The main goals were to 1) properly account for the PSF effect in very deep imaging of extended objects, and 2) reproduce their disc components in unprecedented detail. 2. 2.

We used this technique to reproduce 2D PSF-deconvolved galaxy models (including low surface brightness structures) and to obtain thin- and thick-disc decompositions by considering physical interactions. The method is highly time-consuming, therefore large galaxy samples cannot be analysed at present. 3. 3.

A careful PSF treatment of observed data is absolutely indispensable for deep imaging of extended objects. The PSF effect (if not accounted for) produces an overestimate of the mass in the galaxy outskirts at $\sim$28 mag arcsec*-2* by a factor of 1.5–2 in our two low-mass sources and of 2.5–4 in our two intermediate- to high-mass galaxies (previous works such as Comerón et al. 2018 did not find this mass excess because they reached $\sim 2$ mag arcsec*-2* brighter than we did). The wings of the PSF affect the structures (not only the thick discs) by adding artificial light to the outer parts, which hinders the determination of the real shapes and sizes of galaxy components. 4. 4.

A logical next step to improve our method is to add a stellar halo component to our fits. This would account for physical interactions between the halo and the disc components, but more computer time is required because the number of free parameters will increase. 5. 5.

Our sample of five galaxies provides evidence for aspects of a wide variety of the formation scenarios of thick discs in disc galaxies. The origin of these structures depends on the galaxy mass and morphology, and on the mass assembly history of a given galaxy:

1. (a)

In 2D fittings of the low-mass galaxies UGC 00931 and UGC 01839, only a single disc component is required. The mass-to-light ratios between the thin and thick discs is $(M/L)_{\rm T}/(M/L)_{\rm t}>$1 for these galaxies, suggesting a redder and older thick disc. These two arguments are in agreement with a model where thick discs were first formed from a turbulent gas disc with intense star formation (Elmegreen & Elmegreen, 2006; Bournaud et al., 2009; Comerón et al., 2014, 2018). An alternative interpretation is a thin disc with low star formation efficiency, which cannot be distinguished from the thick disc (it could even be considered as the same structure). Thus, thick (or thin) discs might not be ubiquitous for all types of disc galaxies (contrary to the result of Comerón et al., 2018). 2. (b)

In general, our results on scale heights and mass ratios of thin and thick discs ($z_{{\rm t}}$, $z_{{\rm T}}$, $M_{{\rm T}}/M_{{\rm t}}$, respectively) confirm the results of Comerón et al. (2018) (see Figs 8 and 9). In the case of the one galaxy we have in common, NGC 1032, the main discrepancy is in the values of $z_{{\rm t}}$ and $z_{{\rm T}}$, which are perhaps due to a young stellar population component that the near-infrared data cannot trace. 3. (c)

For NGC 0429 and UGC 01839, $(M/L)_{\rm T}/(M/L)_{\rm t}$ is almost unity, so that both disc components seem to have similar stellar populations. This suggests a common origin for the thin- and thick-disc components.

Acknowledgements.

We thank Ignacio Trujillo, Sebastién Comerón, Victor Debattista, and Lee Kelvin for many useful suggestions and comments. We acknowledge support from the Spanish Ministry of Economy and Competitiveness (MINECO) under grant numbers AYA2013-41243-P and AYA2016-76219-P. J.H.K. acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under Marie Sk$\l$odowska-Curie grant agreement No. 721463 to the SUNDIAL ITN network, from the Fundación BBVA under its 2017 programme of assistance to scientific research groups, for the project “Using machine-learning techniques to drag galaxies from the noise in deep imaging”, and from the Leverhulme Trust for the award of a Visiting Professorship at Liverpool John Moores University. CML thanks the Astrophysics Research Institute of Liverpool John Moores University for their hospitality, and the Spanish Ministry of Science, Innovation, and Universities for financial support of her visit there, through grant number EST2019-013110-I. This research has made use of different Python packages and of NASA’s Astrophysics Data (NED) System Bibliographic Services which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautic and Space Administration.

Appendix A Bootstrap resampling of the 2D galaxy models

Each 2D model of the galaxies in our sample was obtained using the astronomical software imfit (Erwin 2015). However, as described in Sect. 5.1, each galaxy is made of components that need an analytical function that describes its light distribution. To model just one component, it is necessary to leave at least three parameters free. In the case of two or more functions, the image fits were first performed for the innermost components (i.e. bulge and/or bar) and then, with the parameters for those functions fixed, the parameters of the components in the outer part (i.e. thin and/or thick disc) were left free. The output for a reasonable fit from such an initial run was then used to re-run the code by allowing all (or most) model parameters to be free. However, there is a strong degeneracy of the disc parameters, which can result in multiple possible solutions for the modelling process: two or more different solutions can constitute a similarly good fit for a given galaxy. In order to solve this degeneracy problem, imfit allows bootstrap resampling. Bootstrap resampling is a method for estimating the variability of our results by making a more detailed analysis of the parameter distributions, including potential correlations between free parameters. The combined set of bootstrapped parameter values can be used to estimate confidence intervals, and in consequence, to discern between the multiple good solutions. With this procedure, we were able to ensure that the 2D PSF-convolved model fits in our work are not degenerate.

In the following figures we show a set of bootstrapped parameter values and their confidence intervals obtained from the 2D model fitting of each edge-on galaxy in our sample. The bootstrap resampling performed 1500 iterations of each fitting process. In general, the parameters for the discs were the central intensity ($J_{0}$), the value of which for each pixel comes from line-of-sight integration through a 3D luminosity-density model, generating a projected 2D model image given input specifications of the orientation and inclination to the line of sight (i.e. it is not the counts in the image, therefore is not straightforward to convert into surface brightness units); the inner and outer exponential scale lengths ($h_{1}$, $h_{2}$) in kpc for broken exponential functions; the exponential scale lengths ($h$) in kpc for exponential functions; and the vertical scale heights ($z$) also in kpc. When the disc has distinct thin and thick disc components, the corresponding values are indicated with the subscripts $t$ and $T$, respectively. If there is just one disc component, it is shown with the subscript $d$.

Figures 10, 11, 14, and 15 show the results for NGC 0429, UGC 01040, UGC 01839, and NGC 1032, respectively. For all these galaxies, one single set of parameters is obtained as the only (and not degenerate) possible solution. However, for UGC 00931 (Fig. 12), we obtain two different solutions, mainly due to a high discrepancy in the central intensity ($J_{0,d}$) value. We chose the set of parameters that is most physically meaningful and better fits the surface brightness profiles (see Fig. 13). This is the set with lower central intensity, $J_{0,d}=9.53^{+2.33}_{-1.81}$ values/pixel, $h_{1,d}=3.19^{+0.21}_{-0.22}$ kpc, and $z_{d}=0.90^{+0.03}_{-0.03}$ kpc.