© ESO 2019

1. Introduction

Thick discs are the dim and vertically extended counterparts of canonical thin discs. They were discovered in external galaxies by Burstein (1979) and Tsikoudi (1979). In the Milky Way, the thick disc was found by Yoshii (1982) – although he identified it as the stellar halo – and Gilmore & Reid (1983). In spite of their low surface brightness at the heights where they dominate the disc luminosity, thick discs are fairly massive; their stellar masses typically come second only to that of the thin disc, as shown by Yoachim & Dalcanton (2006) and Comerón et al. (2011, 2012, 2018) for external galaxies and by Fuhrmann (2008), Fuhrmann et al. (2012), and Snaith et al. (2014) for the Milky Way.

The growing consensus is that the bulk of the Milky Way thick disc has an internal origin. Ruchti et al. (2014) found no evidence for a significant accreted component in the Galaxy disc. By studying the stellar velocity dispersions in the solar neighbourhood, Kroupa (2002) deduced that, in the past, star formation episodes in clusters with embedded masses up to ∼10⁶ M_⊙ have occurred in the Milky Way. Stars in such massive clusters must have a large velocity dispersion and the cluster dissolution would have resulted in a thick component. Lehnert et al. (2014) and Haywood et al. (2015) proposed that the thick disc in the Milky Way was born thick out of a turbulent and well-mixed interstellar medium during the first gigayears of the Galaxy. The Milky Way thick disc might however contain a small fraction of accreted stars in streams such as those detected by Gilmore et al. (2002) and Martin et al. (2004).

While several studies of local galaxies have been interpreted to favour an internal origin for thick discs (Comerón et al. 2014, 2015, 2016), we are aware of the existence of at least one retrograde thick disc (in FGC 227; Yoachim & Dalcanton 2005, 2008), which strongly indicates that it was accreted. The retrograde nature of the thick disc of FGC 227 was found by studying its rotation curve. The possibility of thick discs being, at least partly, accreted was already pointed out by Abadi et al. (2003). A recent study of FCC 170, an S0 galaxy in the Fornax Cluster, shows clear signs of an accreted thick disc component, although the bulk of the thick disc is believed to have formed in situ (Pinna et al. 2019).

The number of published thick-disc rotation curves is small due to the difficulties in obtaining spectra of regions with low surface brightness. Furthermore, it is hard to pinpoint the heights at which the thick disc dominates (Comerón et al. 2018). Hence, several rotation curves thought to characterise the thick disc might actually sample the thin disc. To our knowledge, rotation curves thought to belong to the thick disc are available for seven galaxies in Yoachim & Dalcanton (2005, 2008)¹, two in Comerón et al. (2015, 2016), one in Guérou et al. (2016), and one more in Kasparova et al. (2016). This last paper studies single-slit spectroscopy of three edge-on galaxies. The spectra were obtained for two slits: one was placed at the midplane and the second one was parallel to the first one and was placed at a height of 5″ − 7″, depending on the galaxy. However, according to the decompositions in Comerón et al. (2018), only in NGC 5422 was the slit placed high enough so as to include a significant thick-disc contribution. Based on the above references, we estimate the total number of galaxies with available thick-disc rotation curves to be eleven, of which ten have prograde thick discs. Increasing the number of available rotation curves for thick discs is necessary to check whether a population of FGC 227-like thick discs exists or whether FGC 227 is a peculiar object with an exceptionally troubled youth.

In this paper we study the rotation curves of a small sample of galaxies at different heights above the midplanes. The data were obtained for an ongoing Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) survey of edge-on galaxies aimed at studying the kinematics and the stellar populations of a statistical sample of thick discs. Six out of eight galaxies obtained so far in our survey have a distinct thick disc (according to the decompositions in Comerón et al. 2018). As a consequence, we are increasing the number of available thick-disc rotation curves by ∼50%.

We note that the reported properties of thick discs depend on how thick discs are defined. Indeed, in the Galaxy the thick disc is defined by selecting individual stars based on their kinematics and/or composition (e.g. Carollo et al. 2010; Bensby et al. 2011; Bovy et al. 2012; Cheng et al. 2012) or using star counts (e.g. Reylé & Robin 2001; Jurić et al. 2008; Jayaraman et al. 2013; López-Corredoira & Molgó 2014; Wang et al. 2018). In external galaxies, the thick disc is most often defined based on photometric decompositions. The different definitions might explain why the thick disc of the Milky Way is often – mostly in studies where the thick-disc stars are kinematically- and/or abundance-selected – found to have a shorter radial scale length than that of the thin disc, whereas thick discs characterised in extragalactic objects are typically as extended as the thin discs, or more so. The not-quite-overlapping thick-disc definitions have caused some authors to demand that the expression “thick disc” be used for the geometrically defined features only (Hayden et al. 2017). The different thick-disc definitions are tentatively reconciled in models such as those by Minchev et al. (2015), where what is defined as a thick disc in photometric decompositions is actually the sum of the flares of mono-age populations. In these models, the metal-poor, α-enhanced, centrally concentrated thick disc found in the Galaxy would correspond to the inner portion of the thick discs in external galaxies, whereas the outer portion would be made of much younger flared populations. We note however that these models only account for internal secular evolution and do not include the possible accretion of external material.

This paper is organised as follows. We describe our sample in Sect. 2. In Sect. 3 we explain how the data were collected and processed. In Sect. 4 we show our main observational results and in Sect. 5 we present a method to estimate the fraction of retrograde material in the thick discs. We present our discussion in Sect. 6 and summarise our findings and conclusions in Sect. 7.

2. The sample

The sample of eight galaxies studied in this paper was selected from the sample of 70 galaxies in Comerón et al. (2012) where we studied edge-on galaxies in the Spitzer Survey of Stellar Structure of Galaxies (S⁴G; Sheth et al. 2010). The S⁴G was made using the IRAC camera (Fazio et al. 2004) on-board the Spitzer Space Telescope. The galaxies were decomposed assuming two vertically isothermal discs in hydrostatic equilibrium. The eight galaxies that we study here were selected to be observable from the southern hemisphere and to have an r₂₅ radius smaller than 60″ according to HyperLeda² (Makarov et al. 2014). In that way, it was guaranteed that the region between the centre of a target and its edge would be covered in a single MUSE pointing (this is necessary in order to be able to build a rotation curve out of the observations). Some of the properties of the galaxies in the sample are summarised in Table 1. The sample is made of galaxies less massive than the Milky Way, with maximum circular velocities obtained from H I kinematics in the range 99 km s⁻¹ ≤v_c ≤180 km s⁻¹.

Our original decompositions in Comerón et al. (2012) had an over-simplistic point spread function (PSF) treatment and are superseded by those in Comerón et al. (2018), where the extended wings of the IRAC PSF are accounted for (the “correct” IRAC PSF is described in Hora et al. 2012)³. With the updated PSF at hand we found that two of the galaxies in the sample that had been observed with MUSE – ESO 443-21 and PGC 30591 – have no distinct thin and thick-disc components. Another galaxy – PGC 28308 – has a thick-disc component that is dimmer than the thin disc component at all observed heights once the scattering of the thin-disc light in IRAC is accounted for.

As the IRAC and the MUSE PSFs are different, the range of heights at which the thick-disc light dominates in Comerón et al. (2018) is unlikely to be the same as in the MUSE datacubes⁴. The MUSE PSF has a smaller full width at half maximum (FWHM) than that of IRAC. However, when studying thick discs, the extended PSF wings play a role that is more important than that of the PSF core⁵. These wings are not described by the FWHM and are an unknown for MUSE. This is not critical in our case, as integral field data cover a large region, meaning that different heights of a galaxy can be studied in a single exposure. Hence, even if the range of heights where the thick-disc light dominates in MUSE is different to that in IRAC, we still expect to see a significant thick-disc contribution at some heights. We would only lose the thick-disc kinematic signatures if the discs were vertically unresolved. This is clearly not the case for our observations because the rotation curves change with height, as shown in Sect. 4.

3. Data acquisition and processing

The data were obtained with the Integral Field Unit MUSE at Unit 4 of the VLT. The instrument has a 1′×1′ field of view (FOV), a spaxel size, and a spectral range that stretches between 4750 Åand 9351 Åwith a 1.25 Åspectral pixel size. The spectral resolution is R ≈ 1700 in the blue side and R ≈ 3700 in the red side. This translates into a ∼2.5 Åspectral resolution that varies little – but noticeably – with the wavelength, as discussed below. The sample galaxies were observed in service mode between December 2015 and August 2016. Each of the eight galaxies had four 2624 s on-target exposures, except for IC 217, for which the number of exposures was three. The exposures for a given galaxy were all centred at the same position and taken with different rotation angles separated by 90°. Thus, the total exposure time per galaxy was almost three hours, except for IC 217. Off-target sky frames with a 240 s exposure were also taken, but were not used in the reduction (see below).

The data were reduced using version 1.6.2 of the MUSE pipeline (Weilbacher et al. 2012) within the Reflex environment (Freudling et al. 2013). Because of our relatively long on-target exposures, the sky variations between the on- and the off-target exposures were large and the sky frames were unsuitable for sky correction. Instead, the sky was modelled from regions in the on-target exposures that were distant from the galaxy midplane (skymethod=model in the pipeline). This is possible because the galaxies are edge-on, and thus fill only a portion of the field of view (FOV).

In our first data-reduction attempts the sky model was obtained from selecting the spaxels with the least light in a white-light image, that is, an image obtained from collapsing the full datacube along its spectral direction. We however later found that several of our galaxies have regions with significant off-plane ionised gas emission that would cause genuine galaxy emission to be included in a sky defined from white-light images. We thus ran the pipeline twice for each galaxy. The first time, we did it for a spectral window 20 Å in width centred in the galaxy rest-frame Hα. The pipeline created sky masks based on Hα images obtained from these datacubes. They were built by selecting the 10% least luminous spaxels (skyfr_2 = 0.10). Additionally, the pipeline, by default, excludes the 5% dimmest spaxels as possible artefacts. The sky masks based on Hα were then fed into a second data reduction, this time run over the whole MUSE spectral range.

Because the observations were scheduled as a filler, most of the exposures had relatively poor seeing ( to in the combined datacubes). Hence, the auto-alignment recipe of the pipeline, muse_exp_align, often failed. We manually aligned the cubes for each of the on-target exposures by blinking white-light images built from the datacubes. The small angular shifts found with this method were then fed into the muse_exp_combine recipe to obtain a single combined datacube.

The combined datacubes were cleaned from small sky residuals with version 2.1 of ZAP (Soto et al. 2016). ZAP creates a sky mask similar to that made in the main MUSE pipeline. It then performs a principal component analysis of the sky residuals and corrects the datacube for them. In this case the sky was defined to be the 7% of spaxels with the least emission in a Hα image obtained from the combined datacube⁶. ZAP worked well across most of the FOV, but spurious features appeared for the regions with the largest line emission. Specifically, a significant noise was introduced. To circumvent this problem, we found the spectra where ZAP introduced noise by measuring the standard deviation of the fluxes in the range between spectral pixels 600 and 1000 for both the cleaned and the uncleaned cubes⁷. We found those spaxels for which this standard deviation was increased by ZAP. We then replaced the noisy spectra in the cleaned datacube with those from the uncleaned datacube.

We produced white-light images out of the datacubes to obtain an optical-near-infrared image of the area under study (top-left panels in Figs. A.1–A.8). The images were used to create a mask of foreground stars (or globular clusters around the target) and background galaxies. We also masked the regions where the median signal-to-noise ratio (S/N) is smaller than 0.5 in the spectral range between 5490 Å and 5510 Å.

The unmasked spaxels were Voronoi binned or tessellated to produce spacial bins with a median S/N of 25 in the spectral range between 5490 Å and 5510 Å with version 2.6 of the Voronoi binning code by Cappellari & Copin (2003). This kind of binning is adaptive, meaning that regions with a high surface brightness have smaller spacial bins than those with a low surface brightness. In many cases, individual spaxels close to the midplanes exceeded the S/N requirement and were left unbinned by the code. The number of spacial bins created per galaxy was typically a few thousand. The galaxy with the fewest bins is IC 217, with 847, and that with the most bins is IC 1553, with 7008. A spectrum was created for each bin by coadding all those in all spaxels within a bin. We call this binning based on the stellar emission the “stellar binning”.

For each spacial bin we obtained a velocity – V – and a velocity dispersion – σ – with version 6.7.0 of the python version of the Penalized Pixel-Fitting code (pPXF; Cappellari & Emsellem 2004). pPXF is a full-spectral-fitting code that fits spectra with a combination of stellar population spectral energy distribution (SED) templates. The output is a parametrised line-of-sight velocity distribution (LOSVD). We masked all the regions with possible strong emission lines and all the regions with sky lines with a peak flux larger than 6 × 10⁻¹⁹ W m⁻² Å⁻¹ arcsec⁻² according to Hanuschik (2003). This number was selected so that the three lines in the calcium triplet remain unmasked (we verified that a three-times-smaller masking threshold does not change our kinematical maps significantly). The width of the masking window was 1000 km s⁻¹. We also masked the O₂-A telluric feature at around 7610 Å with a masking window width of 4400 km s⁻¹.

The SED templates that we used for the fits come from the MIUSCAT library (Vazdekis et al. 2012). We used the models with a Kroupa universal initial mass function (Kroupa 2001) and with Padova isochrones (Girardi et al. 2000). They have a FWHM of 2.51 Åwhich is similar to that of MUSE. However, some of the objects under study here have small masses and hence small velocity dispersions whose determination can be significantly affected if the FWHM mismatch between the MUSE and the template is not accounted for. We also tested in one of our galaxies – ESO157-49 – whether fitting different spectral ranges would affect the results. The selected spectral ranges were “blue” (4750 Å–6000 Å), “red” (6000 Å–8400 Å), and “infrared” (8400 Å–8800 Å; this last spectral window covers the calcium triplet)⁸. In the maps in the top row in Fig. 1 we can see that in general σ is larger in blue than it is in infrared. This potentially indicates that the spectral resolution of MUSE at the calcium triplet is better than that in blue⁹. This is a well-know effect. For example Bacon et al. (2017) modelled the FWHM of the line spread function as

Fig. 1. Velocity dispersion maps of ESO 157-49 for different wavelength ranges (“blue” for the left column, “red” for the middle column, and “infrared” for the right column). Maps in the top row were computed from spectra directly obtained from the processed datacubes. Maps in the lower row were made using spectra and SED templates that have been convolved to a common spectral FWHM as explained in the text. The axes are in arcseconds. The north and the east are indicated by the “backwards-L” symbol.

(1)

with λ in Å based on a MUSE survey of the Hubble Deep Field.

To account for the changes in resolution as a function of wavelength and also for the difference in resolution between the spectra and the templates we did the following: for regions where the instrumental resolution is better than that of the template, we convolved the spectra with a Gaussian kernel of varying width so all the wavelengths had the same resolution as the template. For regions with a worse resolution than the templates, the templates were convolved to the instrumental resolution. This was again done with a varying-width Gaussian kernel to account for changes in the resolution with wavelength. We assumed the FWHM of the line spread function to be described by Eq. (1). Applying such a convolution significantly lowers the σ in the blue map (lower-left panel in Fig. 1). The blue FWHM might be slightly overestimated which might result in a slight underestimation of σ. Indeed, whereas for the unconvolved maps the median difference between the blue and infrared σ values is 6 km s⁻¹, for the convolved ones it is −8 km s⁻¹ (the median values are calculated over all bins and not over the surface)¹⁰.

The final pPXF fits were done for the spectral range between 4750 Å and 8800 Å. The continuum was modelled with an eight-order additive Legendre polynomial. The initial values that we used for the fits were the galaxy radial velocity as reported in the NED and a velocity dispersion of 100 km s⁻¹. The final velocity and velocity dispersion maps are shown in the second and third panels in the left column in Figs. A.1–A.8.

We stress that all the velocities and velocity dispersions obtained with the fits are line-of-sight (LOS) integrated quantities and are thus not representative of any particular radius in the galaxy. This point is important for understanding the analysis in Sect. 5.

We also produced emission line velocity maps. The procedure that we followed is that used in the MaNGA survey (M. Cappellari, private communication). We used the stellar binning to produce new pPXF fits. This time we did not mask the emission lines and included gas line templates to account for the gas emission. We kept the stellar velocities fixed and equal to those in the star-only pPXF fit. The line emission templates were chosen to be Gaussian profiles with a FWHM corresponding to that to which the spectra were degraded at the specific wavelength of the line. The line templates were centred at the wavelengths of the most prominent lines (Hβ, [O III] λ4959, [O III] λ5007, [O I] λ6300, [O I] λ6364, [N II] λ6548, Hα, [N II] λ6583, [S II] λ6716, and [S II] λ6731). All the emission lines were assumed to share the same kinematics. The resulting fitted stellar spectrum was then stored.

We then created a Voronoi tessellation of the spectra based on the Hα emission (we refer to this second tessellation as the “line binning”). We started by creating an Hα line image of the galaxy from the datacubes by integrating in a window of 22 Å centred in Hα and we subtracted the continuum emission calculated in a window of the same width 50 Å redwards. We ran the Voronoi binning code with the condition of S/N = 50 per bin in Hα. We excluded the spaxels where the Hα line had S/N <  0.5. The extracted spectra were again fitted with pPXF. For each bin we used as a stellar template the result of the pPXF fits in the previous step (that where gas and stars were fitted using the stellar binning) whose bin centre was the closest to that of the line bin under study. We kept the velocity of this template fixed and we only allowed it to change by a normalisation factor. This is the fit from which we obtained the final gas kinematics. The line velocity maps are shown in the lower-left panels in Figs. A.1–A.8.

We used the line emission velocity maps to estimate the circular velocity of galaxies as well as the position at which the rotation curve starts to flatten, x_f. To do so we took the midplane gas rotation curve and folded it through the galaxy centre. Prior to the folding, the rotation curves were resampled to a regular sampling. Subsequently, as in Yoachim & Dalcanton (2005), we fitted the rotation curve with the following expression (the same as in Courteau 1997, when setting their β = 0)

(2)

where v_c is the circular velocity, χ is the axial coordinate normalised by the position of the flattening x_f (χ ≡ x_f/x), and γ parametrises the sharpness of the transition from a rising to a flat rotation curve. The larger γ, the sharper the transition. The fits are shown in Fig. 2. The resulting v_c and x_f values are reported in Table 2.

Fig. 2. Symmetrised midplane rotation curves extracted from the line emission velocity maps in Figs. A.1–A.8 (black line) and their fit with Eq. (2) (red line). All panels cover the same physical scale (horizontal axis) and the same velocity scale (vertical axis). The lower-right corner of each panel shows the best fit parameters where vc is the circular velocity, xf is the axial coordinate of the flattening of the rotation curve, and γ describes the sharpness of the transition between a rising and the flat sections. The horizontal dashed lines indicate vc and the vertical dashed lines indicate xf. The panels are ordered by descending vc.

We compared our v_c values to those obtained from H I surveys in the literature (Table 1). We find that for some galaxies both values have significant discrepancies. This might be attributed to the fact that we are only sampling the inner parts of rotation curves, which in turn might cause us to poorly constrain v_c. However, this might not be the case for the galaxy with the largest discrepancy, PGC 28308 (v_c = 130 km s⁻¹ for us and v_c = 180 km s⁻¹ in the literature). The Extragalactic Distance Database (EDD; Tully et al. 2009) H I line profile for this galaxy is not symmetric, which might indicate that the galaxy is disturbed in its outskirts or that the line comes from the superposition of two objects. However, within the radius that we have studied with MUSE, our maps indicate undisturbed kinematics (Fig. A.7). At larger radii, the S⁴G images show that the galaxy is asymmetric and has a warp on its Eastern side. PGC 28308 has been identified to form a pair with MCG-02-25-019 in the HIPASS catalogue (Meyer et al. 2004). The angular separation between the two galaxies is , which is equivalent to 59 kpc in projection at an estimated distance of 43.1 Mpc.

To fully analyse the velocity maps, we produced velocity curves at different heights. Initially the curves were made at the midplane and at the heights ±z_c1 (above and below the midplane) at which the thick disc starts to dominate the surface brightness budget in Comerón et al. (2018). We then created additional velocity curves at equidistant heights between the midplane and z_c1. The distance between cuts was set to be about . We created more rotation curves at heights z >  z_c1 with the same vertical spacing as above. Where no well-defined z_c1 height existed we produced rotation curves at heights separated by . The adopted position angles (PA) came from Salo et al. (2015). We did not use the PA values in Comerón et al. (2018) because we found that for ESO 157-49 the midplane dust lane was obviously misaligned with the Comerón et al. (2018) PA by . For the other galaxies the differences are smaller, ranging from 0° to . The rotation curves are shown in the top-right panels in Figs. A.1–A.8.

The panels below those displaying the rotation curves in Figs. A.1–A.8 show the differences in circular velocity between the different vertical cuts and that at the midplane. The curves were produced by measuring the mean velocity at a given height for a given axial range. The size of the axial bins was set to be between and . The exact bin size was set so that an odd number of bins would fit within the region where the rotation is not flat (marked by the vertical dashed lines in the plots).

The third row of panels in the right column in Figs. A.1–A.8 shows the mean velocity of the rotation curve as a function of height at the regions beyond 1.3 times the flattening radius x_f. This has been computed for both the blue- and the red-shifted halves of the galaxies when available. We could have chosen to compute these curves beyond x_f, but for rotation curves where the transition between the rising and the flat sections is not sharp, regions at x ≈ x_f show a substantial slope. Ideally, we would have liked to put the limit at 1.5 x_f or more, but the MUSE FOV did not cover this region in all our galaxies.

In the lower-right panels in Figs. A.1–A.8 we compare the stellar rotation curves at the midplane and at the maximum heights under study with those of the ionised gas.

4. Results

The regions studied in this paper are largely the inner parts of galaxies. The maximum radius at which rotation curves have been obtained is almost always confined within r₂₅. The one exception is ESO 443-21. The fact that for ESO 443-21 we can obtain kinematics out to 1.5 r₂₅ with the same exposure time as for the other galaxies suggests some inaccuracy in the HyperLeda r₂₅ measurement for this galaxy.

Since we are looking at the inner parts of galaxies, we are studying a region where the different definitions of thick disc do not diverge significantly.

For all galaxies in the sample we are able to reach the region where the rotation curve starts to flatten (we can study regions at x ≳ x_f). However, our rotation curves are never traced far beyond the flattening. In Fig. 3 we present the velocity for x >  1.3 x_f as a function of height averaged over the two or four quadrants covered by pointings. The data points are computed at the heights of the rotation curves in the top-right panels in Figs. A.1–A.8. We did so by averaging the velocity values obtained for all the relevant bins (we gave an equal weight to each bin). The grey bands indicate the standard deviation in the velocity measurements for the bins corresponding to a given data point. Fig. 3 also shows the velocity dispersion σ calculated as

Fig. 3. Symmetrised velocity and velocity dispersion as a function of height (red-orange lines) for the regions at x >  1.3 xf. The plots are ordered by descending maximum circular velocity at a height z = 0 as calculated from the parametrisation of the midplane gas rotation curves (Fig. 2 and Table 1). The grey bands correspond to the standard deviations of velocities and velocity dispersions at the bins used to obtain each point in the curves. The vertical solid lines show the height where the thick disc starts to dominate, if applicable to the galaxy. The black curves indicate fits to the velocity and velocity dispersion maps obtained with the formalism presented in Sect. 5.3 for the galaxies that have two discs. For these fits the non-convolved vertical profiles from photometric decompositions of Comerón et al. (2018) were adopted and then convolved with the MUSE FWHM (Eqs. (21) and (22) and Table 3). The horizontal axes are set so they stretch the same vertical span for all galaxies (1.85 kpc). The vertical axis is also the same in all the panels.

(3)

where σ_i are the velocity dispersion values in the individual bins and N is the number of bins. The grey bands indicate the standard deviation in the velocity dispersion fitted values.

One way to detect a retrograde thick-disc component would be to find a large velocity dispersion in a thick disc. If a significant retrograde component were present, we would expect to find velocity dispersions of the order of the galaxy circular velocity. This only occurs in PGC 28308, where high σ values at large heights and small radii can be seen (Fig. A.7). This is probably due to the central spheroid and not to a retrograde thick disc. Indeed, PGC 28308 has a central mass concentration with a mass of 10% of the baryonic galaxy mass, which makes it the largest one in the sample (Comerón et al. 2018). The other galaxies fail to show a coherent cluster of bins with a large velocity dispersion.

Small fractions of retrograde material are hard to detect. The presence of such material would slightly increase the velocity dispersion and would cause some lag in the rotation curves. However, due to asymmetric drift, both the velocity dispersion and the lag increase with height even for a fully prograde disc. For example, in Comerón et al. (2015) we made N-body models of ESO 533-4 and found that lags of a few tens of kilometres per second are to be expected at the height of the thick disc due to asymmetric drift and a non-perfectly-edge-on orientation. ESO 533-4 has a maximum circular velocity of ∼150 km s⁻¹, comparable to that of the most massive galaxies investigated here. In all our galaxies the lag in the rotation curve at x >  1.3 x_f is of a few tens of kilometres per second at the height of the thick disc. We discuss how to disentangle the effects of asymmetric drift from those of retrograde motions in Sect. 5.

It is possible that the vertical velocity gradients shown in Fig. 3 could be affected by dust extinction (see e.g. the effect of different extinction levels in rotation curves in Bosma et al. 1992; Zasov & Khoperskov 2003). We checked this possibility by comparing the rotation curves as obtained from the blue and the infrared spectral ranges (see Sect. 3) at different heights. The results are shown in Fig. 4 where the dark red lines represent the difference between the absolute velocities determined in blue an infrared at the midplane and the yellow lines represent the same at the height where the thick disc starts to dominate the surface brightness in Comerón et al. (2018). For visualisation purposes, the curves were smoothed with an 8″ window. Uncertainties were computed as the error of the mean over the smoothing window. We use the velocity errors as computed by pPXF and assumed that the errors in blue and infrared are not independent. As expected, at large heights dust has little to no effect and the difference between the two velocities is close to zero. At the midplane, this is not always the case. In ESO 157-49, ESO 544-27, and IC 1553 the plot shows two humps around galactocentric projected distance x = 0. These indicate that the rotation curve grows faster in the infrared than in blue, which in turn indicates a large dust extinction in the central parts of the galaxy.

Fig. 4. Difference between the absolute infrared and blue rotation velocities at the midplane (thick dark red line) and the height where the thick-disc light starts to dominate (thick yellow line). Uncertainty intervals are indicated by the thin lines. The position at which x = 1.3 xf is indicated with dashed lines. Only the galaxies for which the thick-disc light dominates a height range are shown. The curves have been smoothed with an 8″ window.

However, our interest here is to look at the behaviour of the rotation curve at x >  x_f. There, large differences between the infrared and the blue velocities (larger than |V_IR| − |V_Blue| = 10 km s⁻¹) are always associated with large error bars that make those velocity differences compatible with |V_IR| − |V_Blue| <  10 km s⁻¹. We conclude that the midplane velocities – and hence the asymmetric drift – are underestimated at a level of ∼10 km s⁻¹ at most in the flat regions of the rotation curves. The effect would however be small, in accordance with previous findings that the outer parts of the discs are optically thin (Bosma et al. 1992).

In areas of the galaxies with radii smaller than that of the flat rotation curve region, some galaxies show large lags. In one, PGC 28308 (Fig. A.7), this is because of the presence of the spheroid (just as can be seen in NGC 3115; Guérou et al. 2016). ESO 469-15 also has a similar peak in the lag curve at x ≈ 10″ (Fig. A.3). Although this could also be attributed to a small spheroid, none were fitted in Comerón et al. (2018). A large lag of up to 60 km s⁻¹ is observed in some regions of IC 1553 at an axial coordinate x ≈ 10″ − 20″ and at a height of z ≈ 8″ (Fig. A.6). This lag appears well above the limit where the thick disc starts to dominate the luminosity. It might correspond to a small fraction of retrograde material, but no significant increase in the velocity dispersion is found in the lagging bins. IC 1553 also shows a significant asymmetry and a prominent warp (Fig. A.6). It is tempting to link this warp to a galaxy covered by the MUSE FOV and found to the southeast of IC 1553. However we found it to be a background galaxy, at a redshift z ≈ 0.05 (vs. z ≈ 0.01 for IC 1553). According to the NED, IC 1553 has no large nearby galaxy with a comparable recession velocity.

5. Estimating the contribution of retrograde thick-disc stars

We assume that the galaxies consist of a thin- and thick-disc component, and denote by f_T the relative contribution of the thick-disc stars at a given height, and by , and σ_t, σ_T, the mean LOS velocities and velocity dispersions of the thin and thick components, respectively. Moreover, we assume that a certain fraction f_r of the thick-disc stars are on retrograde orbits: these stars are assumed to have the same velocity dispersion as the prograde thick-disc component and a mean velocity . This implies that the effective mean value of the thick-disc component is . The combined thin+thick disc LOS velocities V and velocity dispersions σ can then be written as

(4)

(5)

(6)

In what follows we use the observed 2D velocity and velocity dispersions maps to estimate f_r. To do that we approximate the circular velocities with the gas rotation curve and apply Jeans’ equations to account for the velocity dispersion support which reduces the mean tangential velocity needed to balance the central gravity (“asymmetric drift”). The derivation of the asymmetric drift equation is summarised in Sect. 5.1, following Binney & Tremaine (2008). We also need to deal with the difference between the actual quantities and their LOS-integrated observed counterparts: we use a self-consistent N-body model to check how to apply the asymmetric drift equation to LOS velocities (Sect. 5.2). We then devise a method to fit the fraction of retrograde thick-disc stars, and use the N-body data to estimate the different sources of uncertainty, the main contribution coming from the uncertainty in f_T(z). Finally, in Sect. 5.3, we apply this method to MUSE data, using the photometric thin/thick-disc decompositions made in Comerón et al. (2018), taking into account the effect of the PSF on the observed vertical profiles.

5.1. The effect of asymmetric drift and retrograde stars

For an axisymmetric steady-state galaxy the radial component of Jeans’ velocity-moment equation reads (see Eq. (4.222a) in Binney & Tremaine 2008)

(7)

where R, z, and ϕ are the radial, vertical, and azimuthal coordinates, and v_R, v_z, v_ϕ the corresponding velocity components; brackets indicate averaging over a small volume around R, z, ϕ. Here n = n(R, z) stands for the number density of stars and Φ(R, z) denotes the gravitational potential. In a galaxy with thin and thick-disc components, stars of both components separately obey this equation. Following Binney & Tremaine (2008), we assume a zero mean radial velocity so that

(8)

separate the mean azimuthal velocity and the azimuthal velocity dispersion,

(9)

and write the derivative of the potential in terms of circular velocity v_c,

(10)

Eq. (7) can then be recast to the form

(11)

To estimate the radial and vertical gradient terms we assume exponential profiles for the number density and the radial velocity dispersion

(12)

and

(13)

where h_R and h_z are the disc scale length and scale height and h_v is the radial scale length of the velocity dispersion. We then approximate

(14)

where the multiplier k₀ describes how the orientation of the velocity ellipsoid depends on vertical height; k₀ = 0 corresponds to cylindrical rotation, while k₀ = 1 corresponds to a situation where the principal axis of the velocity ellipsoid points to the galaxy centre. Inserting these expressions into Eq. (11) finally leads to

(15)

This deviation of the mean tangential velocity from the circular velocity due to velocity dispersion is the so-called asymmetric drift. The lag is larger for the thick disc, which has a larger velocity dispersion than the thin disc. The thick-disc contribution f_T increases with height, thereby making the combined velocity and velocity dispersion profiles depend on height. However, this is not the only factor that can cause a z-dependence: other contributions may come from the possible dependence of circular velocity and velocity dispersion on height. Moreover, Eq. (15) contains a term with an explicit z-dependence.

Figure 5 shows a simple toy model illustrating how the vertical velocity and velocity dispersion profiles can help to indicate the presence of a retrograde thick-disc component. We assume that the thin- and thick-disc components have the same total flux, and that they follow exponential vertical distributions with scale heights h_zt and h_zT; this specifies the relative contribution f_T(z). If we ignore for a moment the difference between actual and LOS velocities, we approximate Eq. (15) with

Fig. 5. Toy models for the influence of asymmetric drift and retrograde thick-disc stars on the vertical velocity and velocity dispersion gradients at |z| = 3hzt. We use Eqs. (16) and (17) with vc = 1, ki = 1, σt = 0.2, assuming hzT/hzt = 1 − 4 and a fraction of retrograde stars fr = 0 − 0.9. The crosses indicate hzT/hzt = 3. Circles stand for the N-body simulation fits of Fig. 8 below. Relatively small retrograde fractions of fr = 0.1 − 0.2 cause lags of a few tens of per cent if the thin and the thick discs are well differenciated (i.e. hzT/hzt ≤ 3, as is the case for the galaxies in our sample).

(16)

and assume that velocity dispersions depend on scale height as

(17)

Fixing v_c = 1, σ_t = 0.2, and for simplicity using k_i = 1, we use Eqs. (4)–(6) to calculate the expected vertical gradients in combined mean velocity and velocity dispersion, for a set of h_zT/h_zt and f_r values. In Fig. 5 the gradients are measured by the relative change of V and σ between z = 0 and |z| = 3h_zt. In the absence of a retrograde component, the velocity drop in the vertical direction is due to the stronger asymmetric drift of the thick-disc component which becomes dominant at large height. In this case a drop in with height is achieved only in combination with an increase in σ. On the other hand, a larger relative decrease in is expected in this toy model if the reduction of the mean velocity with height is due to the retrograde contribution to the effective mean velocity of the thick disc.

5.2. N-body modelling

We next turn to an N-body model to assess, in a more quantitative manner, how to apply this approach to actual data. Figure 6 displays a simulation model with superposed thin- and thick-disc components, with equal masses and with scale heights differing by a factor of three. As in Salo & Laurikainen (2017), the initial values have been constructed with the GalactICS software (Kuijken & Dubinski 1995). This code was used to create two separate, self-consistent axisymmetric bulge-disc-halo models, with similar disc-to-total mass ratio and total rotation curve but different vertical thickness of the disc. The disc radial velocity dispersions were also different for the two models and were set by the ratio σ_z/σ_R ≈ 0.8; the Toomre Q at two radial scale lengths was about 1.4 and 2.8 for the thinner and the thicker disc models, respectively. The two disc models were then superposed together with the bulge and the halo of one of the single-disc models. The combined system was evolved with GADGET-2 (Springel 2005) for 200 time units (a time unit is equivalent to the dynamical timescale at one exponential scale length, or about 10 Myr) to ensure that the two discs are mutually relaxed and also with respect to the halo and bulge. No significant non-axisymmetric perturbations developed during this time, meaning that the simulation is expected to fulfill the axisymmetric steady-state Jeans’ equations. The parameters of the N-body model were chosen in a manner that should be relevant to the observed, fairly low-mass galaxies. In the simulations the total disc-to-halo mass ratio is approximately 0.25 within five scale lengths, with the thick and thin components contributing equally to the disc mass. The vertical thicknesses are 0.15 and 0.45 times the radial scale length (we use the initial radial scale length as a unit in all simulations). The model also has a small classical bulge component, but since this amounts to only 0.5% of the disc mass its influence is insignificant. Both thin and thick disc are realised with 10⁶ particles. This number is sufficiently large to allow the construction of mean stellar velocities and velocity dispersions, both locally at any point (R, z), and along a line of sight with a given axial distance x at height z. Using the accelerations calculated with GADGET-2, we can also reconstruct the circular velocities at the locations of particles; we use this to create a “gas-particle” population following circular orbits. To reduce the noise, five simulation snapshots are superposed when calculating kinematic quantities.

Fig. 6. Illustration of the N-body model used for checking the method for extracting the fraction of retrograde thick-disc stars. Upper left panel: frame edge-on snapshot of the simulated disc, while the other two frames display the thin and thick components separately (each with 106 particles). Upper right frame: combined LOS-averaged rotation curves: black and grey lines correspond to the circular velocities at z/hR = 0 and 0.9, respectively, while the circles indicate the same for stars. Lower right frame: combined vertical velocity and velocity dispersion profiles, averaged over the radial range 1.8 − 2.2 hR, indicated by the white box on the N-body snapshot. The dashed line stands for the circular velocity, open circles for , and squares for σ. The continuous blue and dash-dot red curves show and σ for the thin and thick discs, respectively.

With this N-body model we can numerically check several factors affecting the application of Jeans’ equations. First of all, the asymmetric drift equation is indeed well satisfied at any point for the two components. In Fig. 7 the thick blue and red curves indicate the vertical v_ϕ profiles obtained directly from velocities of simulation particles, while solid circles indicate Eq. (15), when numerically evaluating all the different terms contained within it. To obtain this agreement, the drop of circular velocity with height is important to include, indicating that the same effect must be accounted for when modelling the vertical profiles of observed LOS velocities. The precise shape of the velocity ellipsoid has a smaller but still noticeable effect (compare diamonds with filled circles; for clarity, only shown for the thick disc), and the same concerns the radial gradient terms (compare open and filled circles). The comparison in Fig. 7 is made at the region where the rotation curve is turning flat: in this region k₀ ≈ 1; however, ignoring this term has little effect.

Fig. 7. Example of application of Eq. (15) to N-body data. Frame (a) displays velocities and velocity dispersions measured in a radial annulus at R = 1.8 − 2.2, while in (b) LOS-integrated quantities at axial distance x = 1.8 − 2.2 are shown. In (a) the dashed line indicates the circular velocity as a function of z, and the thick blue (above) and red curves (below) are the mean tangential velocities for the thin and thick discs, respectively. Thin curves are the radial-velocity dispersions. The symbols refer to variations in the application of Eq. (15): the filled black circles include all terms, evaluated numerically at the radial annulus (hRt = 1.1, hRT = 1.2, hvt = 1.4, hvT = 1.5, k0 = 1; measured σR(z), σϕ(z), and σz(z) values are used). For the thick disc, other approximations are also shown. The crosses correspond to the case where σR(z), σϕ(z), and σz(z) are replaced with their mean values for |z| <  2hzt (0.27, 0.23, and 0.28, respectively; this improves agreement, probably because the values are less noisy). Diamonds indicate σz/σR = 0.7 and σϕ/σR = 0.7. Finally, open circles correspond to hR = 1, hv = 1. Frame (b) shows LOS integrated quantities separately for thin (blue curves; below) and thick discs (red; above). For the thin disc, the solid circles are calculated using the corrected σt(c) (see the text), while open circles use the non-corrected velocity dispersion, dominated by the large range of mean velocities along the LOS. In both cases, Eq. (15) is applied with the hR and hv values measured from the LOS-integrated radial profiles, while σz/σR = 0.7, σϕ/σR = 0.7, and k0 = 1 are assumed.

The right-hand frame of Fig. 7 shows a similar comparison for the LOS-integrated velocities, substituting and σ_i for v_ϕ and σ_R in Eq. (15). Though the trends are similar to those seen when using actual velocities, the amount of predicted asymmetric drift is clearly too large. In particular, for the thin disc it is several times too large. This follows from the fact that since the intrinsic velocity dispersion of the thin disc is small, σ_t is dominated by the large range of mean velocities traversed in the LOS integration; for the thick disc this effect is smaller (compare the intrinsic and LOS velocity dispersions in Fig. 7a and b). To some degree this extra dispersion can be accounted for by using a correction, where the gas-velocity dispersion is quadratically subtracted from the thin-disc velocity dispersion,

(18)

This correction (dashed curve in Fig. 7b) assumes that gas has zero intrinsic dispersion and that its LOS dispersion arises solely from the mixture of mean velocities (this applies exactly to our simulated “gas” particles). Clearly this correction improves the applicability of Eq. (15) to the thin disc.

Our N-body model allows us to mimic a variable degree of retrograde rotation. To achieve this we reverse the velocities of a fraction f_r of randomly selected thick-disc particles. We also checked that the Eqs. (4)–(6) are satisfied by such modified particle realisations, basically indicating that the first two moments are sufficient to describe the LOS velocity distributions with good accuracy.

Armed with this experience we next devise a method for extracting the fraction of stars in retrograde orbits f_r from synthetic N-body data, using only quantities which can be obtained also from real data. We write

(19)

(20)

The difference with the toy model is that we now take into account the vertical drop of v_c, and include σ_t and σ_T as unknowns; a₀ is an unknown scale parameter between the LOS and the actual velocities. The function f(x) collects the velocity ellipsoid and radius-dependent terms of Eq. (15): replacing R with x, and assuming and k₀ = 1 this function becomes f(x)= − 1 + x/h_R + x/h_v. Combining these expressions with Eqs. (4)–(6) we form models for V(x, z) and σ²(x, z), and use a χ² minimisation to estimate the parameters a₀, σ_t, σ_T, f_r which best match the synthetic profiles constructed from N-body realisations with different f_r. The fraction f_T is measured from the simulation, and the scale heights are fitted to the simulated surface density profiles (in practice they are close to the input values used in GalactICS, even after the relaxation of the discs).

Figure 8 shows examples of fitted profiles overplotted on synthetic data, for various amounts of retrograde thick-disc stars, while Fig. 9 collects the fit results to a full range f_r = 0.0, 0.1, …, 0.9, 1.0. In the example shown, the synthetic profile is constructed at axial range x/h_R = 2.0 − 2.5, and the above model is used with f(x)=3.5. Figure 9 shows that altogether the fit recovers retrograde fractions remarkably well: the deviation |f_r(fit)−f_r(data)| stays well below 0.01. Also the fitted values of σ_t and σ_T are practically the same for all values of f_r and within 10% of the input N-body values.

Fig. 8. Examples of fitting the simulated LOS velocity (filled circles) and velocity dispersion profiles (open circles) when 0, 20, and 40% of thick-disc stars are in a retrograde component. The curves are the fitted values using the fT(z) ratio from simulated vertical profiles.

Fig. 9. In black, fitted retrograde fraction for a range of input fractions; boxes indicate the input values. Continuous blue and dash-dot red curves are the fitted “effective” velocity dispersions for the thin- and thick-disc components, respectively.

The values obtained for this N-body model are also marked in Fig. 5, overplotted on the toy-model curves. Although there is a systematic offset compared to the toy model, the qualitative trends are quite close to those in the N-body model. The only exception is the behaviour for f_r close to zero, where the actual vertical gradient in σ is much smaller than what the toy model suggests. This is related to the above-discussed effect of LOS integration on the thin-disc velocity dispersion, which hides the actual intrinsic velocity dispersion. In this respect, observed galaxies can be expected to be close to N-body data.

The above success in retrieving f_r from N-body data is not fully representative of the application to true data, in which many uncertainties are present. To see how robust the fitting is, we experiment by varying the axial range over which the fit is done and also assess the effects of using inappropriate values for the term f(x). Moreover, we introduce an uncertainty into the f_T factor because it is the most uncertain in real data. In Fig. 10 various factors are explored. In the upper left, synthetic data from different axial distances are used, while in the upper right, the details of the fitted model are varied; both have negligible effect on the results. These and similar tests with varying f(x) (e.g. dropping the explicit z dependence, or omitting the a₀ parameters) indicate that the precise form of the model function for is not important for recovering the f_r fraction.

Fig. 10. Various factors affecting the error between the actual input fraction of retrograde thick-disc stars in the N-body model, and the fraction returned by the fit. Upper left panel: fits with Eq. (20), using data from various axial distance intervals as indicated by the labels. Upper right panel: different variants of Eq. (20) are tested, using the data from the axial interval x/hR = 2.8 − 3.2: compared to the nominal case (“model=1”), in “model=2” the latter term in Eq. (20) has been omitted, while for “model=3” the has been replaced with , to correct for the contamination of LOS-velocity dispersions by systematic velocities. The middle panels (notice the change in the y-axis range) show the effect of using an incorrect fraction for thick disc, illustrated in the corresponding panel in the lowermost row. In the left panel frame the shape of fT is retained, but the contribution of the thick disc in the equatorial plane has been multiplied with a factor indicated by the labels. In the right panel frame, the fT profile has been stretched/dilated in the vertical direction by the amount indicated by the labels. The dashed lines indicated the fit obtained using the fT estimated from a photometric thin-thick disc decomposition to the simulation data.

On the other hand, the middle row of Fig. 10 explores uncertainties due to the adopted thin/thick decomposition, quantified by using different f_T profiles, where either the fractional contribution of thick-disc stars at the equatorial plane is multiplied by some factor (left), or the f_T profile is stretched/dilated in the vertical direction by some factor. Also shown is a photometric decomposition made on N-body data, using the same method that was used in Comerón et al. (2018) for observed galaxies (see below). In this photometric decomposition, the central plane contribution of the thick disc was overestimated by a factor of approximately 1.25 (36% instead of 29%) and the f_T profile was additionally stretched by a factor of about 1.1 compared to the actual N-body profile. As expected, the errors in the fitted f_r fraction are considerably larger than those caused by the details of the fitting model. Nevertheless, even assuming relatively large errors, the ratios f_r are retrieved with a better than 10% accuracy.

5.3. Application of the formalism to the MUSE data

We are now ready to estimate retrograde fractions from MUSE datacubes. In the application to MUSE data we use the density profiles from the thin-thick disc decompositions made in Comerón et al. (2018) for the 3.6 μm S⁴G data. These decompositions were based on the assumption that the two components can be represented with isothermal models in mutual hydrodynamical equilibrium. The effective scale heights and the common radial scale lengths from that paper are given in Table 3. We start the fitting the stellar velocity and σ maps at x >  1.3x_f using Eqs. (19) and (20) to obtain V_t and V_T; we use the observed gas velocities as a proxy for circular velocities. To obtain the combined V and σ we use the non-convolved model profiles ρ_t(z) and ρ_T(z) from Comerón et al. (2018), and convolve the first and second velocity moments with the 1D MUSE PSF,

(21)

(22)

yielding . We assumed the MUSE PSF to be a Gaussian with the FWHM determined using IRAF’s IMEXAMINE¹¹. The adopted FWHM values are shown in Table 3.

The application of our formalism to six MUSE galaxies is shown in Fig. 3. The labels in the plot indicate the fitted f_r, together with an upper and a lower bound obtained by assuming that the photometric profiles are stretched/dilated by factors of 1.3 and 0.7 with respect to the real one: this corresponds to roughly twice the error in the similar photometric fit to the N-body data. In all cases the retrograde fractions, allowing the above uncertainty, are consistent with zero.

It is possible that our stellar velocity dispersion determinations are slightly underestimated if the FWHM of the observations has been overestimated. However, we have tested whether a negligible retrograde material function would also be found if we assume an underestimated FWHM. For this experiment we modelled the wavelength-varying FWHM of the instrumental line spread function provided by the MUSE pipeline calibration files. These files describe the profile of the lines in the pixel tables. When pixel tables are processed to build the datacubes this results in a small widening of the FWHM. The FWHM derived from the calibration files is ∼0.2 Ånarrower than that estimated by Eq. (1) in the blue side of the spectrum and ∼0.1 Ånarrower at the calcium triplet. Even when such an underestimated FWHM is used, only IC 217 has a mildly significant positive fraction of retrograde stars ().

We conclude that none of the six galaxies with a clear thick disc exhibit a significant retrograde thick-disc component.

6. Discussion

The kinematic data that we present show that retrograde material is not generally dominant at any height in the discs of the galaxies under study. For a couple of galaxies there is some evidence for strongly lagging – non-rotating or counter-rotating – material close to the centre, where a spheroid could be found. However, such evidence vanishes at large distances from the galaxy centre, where the rotation curves become flat.

With the numbers at hand (from this paper and from previous studies), only one thick disc in seventeen – that in FGC 227 (Yoachim & Dalcanton 2005, 2008) – necessarily contains a large fraction of retrograde material.

Large amounts of retrograde thick-disc material would strongly support an external origin. The opposite, that a lack of such material excludes an accretion origin, is not necessarily true. The reasons for that are related to the way dynamical friction (Chandrasekhar 1943) works, as shown below.

To first order, we can think of galaxies as rotating discs within a non-rotating or slowly rotating, pressure-supported, dark matter halo. At large distances from the galaxy centre, an accreted object is mostly influenced by the friction caused by the non-rotating halo and its behaviour does not depend on whether the accreted object is in a prograde or a retrograde orbit (see, e.g. the mass-loss rates as a function of time in Peñarrubia et al. 2002, when the infalling galaxy is far away from the centre of the halo). Differences between these two kinds of orbit arise when the interaction with the disc increases, as dynamical friction is more efficient for an accreted object with a velocity similar to that of the disc (a situation of near-resonance). Hence, the disc of a galaxy is more efficient at braking a prograde satellite. This has two effects. The most obvious one is that mergers are faster in the prograde case (but only by a factor of at most two; Peñarrubia et al. 2002; Villalobos & Helmi 2008). The second is that the accreted objects are strongly dragged towards the main galaxy disc plane in the case of a prograde merger with a low inclination (i <  20°; Peñarrubia et al. 2002; Read et al. 2008). Therefore, prograde mergers are the most likely to result in stars deposited into a thick-disc component (as opposed to being deposited into the stellar halo). However, a very-low-inclination retrograde merger would also result in a stellar distribution resembling that of a thick disc (Read et al. 2008).

The fundamental question here is: if thick-disc stars were accreted in a small number of merger events, what would be the likelihood for a given galaxy to have a large fraction of retrograde material? One difficulty in answering this question is that many numerical studies on dynamical friction sample the prograde encounter parameter space much better than they do the retrograde one. However, satellite systems are pressure-supported – with small net rotation – so retrograde and prograde orbits are nearly equally likely (Sales et al. 2007; Read et al. 2008).

In the discussion that follows, we only consider relatively major mergers as potential thick-disc progenitors. Small satellites take a long time to decay because dynamical friction scales with the mass of the infalling object to the second power. Furthermore, when small satellites are finally dissolved, they do so into structures that do not resemble a thick disc (Read et al. 2008).

In Peñarrubia et al. (2002) and Villalobos & Helmi (2008), a massive retrograde object would take at most two times longer than a prograde object to dissolve. Since the decay time for the massive satellites in those simulations is smaller than the Hubble-Lemaître time, both prograde and retrograde objects bound to the galaxy at its formation time are expected to have merged by now. However, among objects falling into the main galaxy potential well at a later time, not all have merged yet, and retrograde ones are the most likely not to have merged. Let us make a rough estimate of the number of prograde mergers versus that of retrograde mergers. For the models in Peñarrubia et al. (2002), a massive satellite with an initial apogalacticon of 55 kpc and a moderate eccentricity e = 0.45 would dissolve in ∼5 Gyr if in a prograde orbit and in ∼7 Gyr if in a polar orbit (we assume that the decay time of a polar orbit is similar to that of a retrograde orbit, as shown for objects with larger eccentricities in Peñarrubia et al. 2002). If we assume the dark matter halo merger rate per gigayear to go like ∝(1 + z)² (e.g. Fakhouri & Ma 2008) we obtain that the fraction of retrograde mergers is ∼45% for a galaxy of 10 Gyr old. Here we have not accounted for the primary galaxy mass evolution – which is likely to result in a smaller Coulomb logarithm at early times and hence into a longer decay time – or for the possibility of nearly circular orbits with a decay time close to a Hubble-Lemaître time for a polar orbit (on the other hand, satellites with an eccentricity larger than e = 0.45 would decay faster than stated here; Peñarrubia et al. 2002). Therefore, to an order of magnitude, we can estimate the ratio of retrograde-to-prograde mergers in a Hubble-Lemaître time to be about one third.

Among the prograde mergers, only the ones with i <  20° would be plane-dragged and dissolve into a thick disc, whereas the others would disperse into the halo. This can also be seen in Qu et al. (2011) where, using GalMer giant-dwarf merger simulations (Chilingarian et al. 2010) with i = 33° and i = 60° retograde encounters, the authors find a very small resulting fraction of retrograde stars at and below z ∼ 1 kpc (their primary galaxy, however, is more massive than those studied here). It is to be expected that for retrograde encounters a smaller range of inclinations would result in a thick disc. In Read et al. (2008) the only retrograde simulation ends up forming a thick disc and has an inclination of i = 10°. Therefore, for the model galaxies in that set of simulations, this is a lower limit for the inclination of retrograde objects forming thick discs. Using that number, and assuming isotropic accretion, we obtain that a retrograde satellite has a chance that is cos 70 ° /cos 80 ° ≈2 smaller than that of a prograde one to end in an accreted thick disc.

If we combine the two factors – that from the decay time and that from the in-plane drag – we obtain that at the present time approximately one sixth of the objects accreted into a thick disc should be retrograde, and the others prograde. We warn that this number is only an order-of-magnitude estimate and certainly depends heavily on the mass distributions involved and on the orbital elements of the merging objects.

If we assume that all thick discs are accreted in a single merger event, a sixth of thick discs should have retrograde material. Considering a binomial distribution with p = 1/6, there is a two sigma discrepancy between our prediction and the observed fraction of thick discs with retrograde material (1 in 17). If we assume that the thick disc is made in a few (two or three) accretion events, the tension worsens, since having any one of these coming from a retrograde orbit would likely leave an imprint detectable with our formalism, as described in Sect. 5.

The above numbers, as well as the inferred origin for the Milky Way thick disc, are evidence against accretion of stars as the main thick-disc formation mechanism. There have been simulations indicating that only up to ∼50% of the masses of thick discs have been accreted (Abadi et al. 2003; Read et al. 2008), but our spectroscopic data and those in the literature seem to indicate that the fraction is lower.

We note that a non-accretion origin does not necessarily imply that thick discs are not created by interactions. Indeed, merger events – even those for which stars of a disrupted satellite end up in the halo – heat the pre-existing thin disc and can participate in the formation of the thick disc (Kazantzidis et al. 2008).

7. Summary and conclusions

In this paper we present new kinematical information for thick discs. Specifically, we observed eight nearby galaxies with the integral field unit MUSE at the VLT. Six of the galaxies have distinct thick discs according to Comerón et al. (2018) and hence we can study the thick-disc kinematics. This increases the number of galaxies with known thick-disc kinematics to seventeen.

In Sect. 5 we present a formalism able to estimate the fraction of thick-disc retrograde material from the velocity and velocity dispersion maps. None of the thick discs in our sample show signs of retrograde star fractions larger than 10%. Within the uncertainties, the retrograde star fraction is probably compatible with zero for the six galaxies.

Of the 17 thick discs with known kinematics only one shows clear evidence for retrograde material, which is a clear signature for an accretion origin. Using information from the simulations in Peñarrubia et al. (2002) and Read et al. (2008) we estimate that approximately one sixth of the satellites accreted into a thick disc had a retrograde orbit. This is in tension with the observed fraction of galaxies with retrograde thick-disc material as estimated from the rotation curves (1 in 17). We are planning new MUSE observations to confirm this result. Undoubtedly, additional new MUSE data for edge-on galaxies will be obtained by other groups and can also be used to increase the sample size.

Our findings, and the lack of evidence for large fractions of accreted material in the Milky Way thick disc, seem to discard accretion as the main thick-disc formation mechanism for most galaxies. The thick-disc origin is now narrowed down to either an internal origin or a scenario where a pre-existing thin disc is dynamically heated by encounters. Either way, the source of most of the thick-disc material is most likely to be internal.

Acknowledgments

We thank the referee for in-depth comments that helped to quantify the fraction of retrograde material in the discs. We thank Dr. Lodovico Coccato for his help at fine-tuning the data reduction. We thank Dr. Simon Conseil for his help with ZAP and Dr. Peter Weilbacher for his suggestions on the line spread function handling. SC and HS acknowledge support from the Academy of Finland (grant No. 297738). HS, JHK, and RFP acknowledge financial support from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No. 721463 to the SUNDIAL ITN network. JHK acknowledges additional support from the Spanish Ministry of Economy and Competitiveness (MINECO) under grant number AYA2016-76219-P, as well as 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 through the award of a Visiting Professorship at LJMU. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Individual galaxies map and plots

In this appendix we show the maps and plots derived for each of the galaxies in our sample. The plots are structured in two columns.

In the left column of plots we show the following.

1. The white-light image obtained from the MUSE datacube.

2. The stellar velocity map; the systemic velocity of the galaxy has been subtracted.

3. The stellar velocity dispersion map.

4. The gas velocity map as measured from emission lines; the systemic velocity of the galaxy has been subtracted.

In the above four frames the diamond symbols indicate the galaxy centre and the north and the east are indicated by the “backwards-L” symbol at the top-right corner. The labels in the axes are in arcseconds. The “NE” and “SW” (“NW” and “SE”) labels indicate the northeastern and the southwestern (northwestern and southeastern) sides of the galaxy, respectively. The dividing line between the two halves is the midplane which is indicated by a thin continuous line. The two thick continuous line indicate the height above which the thick disc dominates the luminosity according to the fits in Comerón et al. (2018). These three lines, as well as the dashed lines, indicate the heights for which velocity curves have been constructed (see the plots in the right column).

The plots in the right column are structured as follows.

1. Velocity curves for the heights indicated by the thin, the thick, and the dashed lines in the maps in the left column of plots. The curves for the NE (NW) half of the galaxy are shown at left and those for the SW (SE) half are shown at right. The vertical dashed lines indicate 1.3 times the position of the rotation curve flattening, x_f. The colour coding indicates the height above the midplane to which a curve corresponds.

2. Rotation velocity difference at a given height with respect to that of the midplane. The curves for the NE (NW) half of the galaxy are shown at left and those for the SW (SE) half are shown at right. The vertical dashed lines indicate 1.3 times the position of the rotation curve flattening according to our parametrisation in Eq. (2), 1.3 x_f. The colour coding indicates the height above the midplane to which a curve corresponds.

3. Average rotation velocity of the flatter part of the rotation curve (x >  1.3 x_f) as a function of height for the red-shifted side of the galaxy (R) and the blue-shifted side of the galaxy (B). The vertical solid lines show the height where the thick disc starts to dominate, if applicable to the galaxy.

4. Stellar and the gas rotation curves for selected heights. The circular symbols correspond to the stellar rotation curves and the lines to the gas rotation curves. Rotation curves for the midplanes are indicated by dark points and blue lines, and rotation curves for the largest height for which we constructed rotation curves are indicated by light points and red lines. The top panel is for the NE (NW) half of the galaxy and the lower panel is for the SW (SE) half. The vertical dashed lines indicate x = 1.3 x_f. The dot lines indicate the zero axial distance to the centre and the zero velocity lines.

Fig. A.1. Maps and plots corresponding to ESO 157-49.

Fig. A.2. As in Fig. A.1, but for ESO 443-21. Comerón et al. (2018) found no distinct thin- and thick-disc structure for this galaxy.

Fig. A.3. As in Fig. A.1, but for ESO 469-15.

Fig. A.4. As in Fig. A.1, but for ESO 544-27.

Fig. A.5. As in Fig. A.1, but for IC 217.

Fig. A.6. As in Fig. A.1, but for IC 1553. The small galaxy at the lower-left of the top-left panel is unrelated to IC 1553 and is found to have redshift z = 0.05.

Fig. A.7. As in Fig. A.1, but for PGC 28308. In this galaxy there is no height where the thin disc has a lower surface brightness than the thick disc once the PSF is accounted for (Comerón et al. 2018).

Fig. A.8. As in Fig. A.1, but for PGC 30591. Comerón et al. (2018) found no distinct thin- and thick-disc structure for this galaxy.

All Tables

All Figures

Fig. 1. Velocity dispersion maps of ESO 157-49 for different wavelength ranges (“blue” for the left column, “red” for the middle column, and “infrared” for the right column). Maps in the top row were computed from spectra directly obtained from the processed datacubes. Maps in the lower row were made using spectra and SED templates that have been convolved to a common spectral FWHM as explained in the text. The axes are in arcseconds. The north and the east are indicated by the “backwards-L” symbol.

In the text

Fig. 2. Symmetrised midplane rotation curves extracted from the line emission velocity maps in Figs. A.1–A.8 (black line) and their fit with Eq. (2) (red line). All panels cover the same physical scale (horizontal axis) and the same velocity scale (vertical axis). The lower-right corner of each panel shows the best fit parameters where vc is the circular velocity, xf is the axial coordinate of the flattening of the rotation curve, and γ describes the sharpness of the transition between a rising and the flat sections. The horizontal dashed lines indicate vc and the vertical dashed lines indicate xf. The panels are ordered by descending vc.

In the text

Fig. 3. Symmetrised velocity and velocity dispersion as a function of height (red-orange lines) for the regions at x >  1.3 xf. The plots are ordered by descending maximum circular velocity at a height z = 0 as calculated from the parametrisation of the midplane gas rotation curves (Fig. 2 and Table 1). The grey bands correspond to the standard deviations of velocities and velocity dispersions at the bins used to obtain each point in the curves. The vertical solid lines show the height where the thick disc starts to dominate, if applicable to the galaxy. The black curves indicate fits to the velocity and velocity dispersion maps obtained with the formalism presented in Sect. 5.3 for the galaxies that have two discs. For these fits the non-convolved vertical profiles from photometric decompositions of Comerón et al. (2018) were adopted and then convolved with the MUSE FWHM (Eqs. (21) and (22) and Table 3). The horizontal axes are set so they stretch the same vertical span for all galaxies (1.85 kpc). The vertical axis is also the same in all the panels.

In the text

Fig. 4. Difference between the absolute infrared and blue rotation velocities at the midplane (thick dark red line) and the height where the thick-disc light starts to dominate (thick yellow line). Uncertainty intervals are indicated by the thin lines. The position at which x = 1.3 xf is indicated with dashed lines. Only the galaxies for which the thick-disc light dominates a height range are shown. The curves have been smoothed with an 8″ window.

In the text

Fig. 5. Toy models for the influence of asymmetric drift and retrograde thick-disc stars on the vertical velocity and velocity dispersion gradients at |z| = 3hzt. We use Eqs. (16) and (17) with vc = 1, ki = 1, σt = 0.2, assuming hzT/hzt = 1 − 4 and a fraction of retrograde stars fr = 0 − 0.9. The crosses indicate hzT/hzt = 3. Circles stand for the N-body simulation fits of Fig. 8 below. Relatively small retrograde fractions of fr = 0.1 − 0.2 cause lags of a few tens of per cent if the thin and the thick discs are well differenciated (i.e. hzT/hzt ≤ 3, as is the case for the galaxies in our sample).

In the text

Fig. 6. Illustration of the N-body model used for checking the method for extracting the fraction of retrograde thick-disc stars. Upper left panel: frame edge-on snapshot of the simulated disc, while the other two frames display the thin and thick components separately (each with 106 particles). Upper right frame: combined LOS-averaged rotation curves: black and grey lines correspond to the circular velocities at z/hR = 0 and 0.9, respectively, while the circles indicate the same for stars. Lower right frame: combined vertical velocity and velocity dispersion profiles, averaged over the radial range 1.8 − 2.2 hR, indicated by the white box on the N-body snapshot. The dashed line stands for the circular velocity, open circles for , and squares for σ. The continuous blue and dash-dot red curves show and σ for the thin and thick discs, respectively.

In the text

Fig. 7. Example of application of Eq. (15) to N-body data. Frame (a) displays velocities and velocity dispersions measured in a radial annulus at R = 1.8 − 2.2, while in (b) LOS-integrated quantities at axial distance x = 1.8 − 2.2 are shown. In (a) the dashed line indicates the circular velocity as a function of z, and the thick blue (above) and red curves (below) are the mean tangential velocities for the thin and thick discs, respectively. Thin curves are the radial-velocity dispersions. The symbols refer to variations in the application of Eq. (15): the filled black circles include all terms, evaluated numerically at the radial annulus (hRt = 1.1, hRT = 1.2, hvt = 1.4, hvT = 1.5, k0 = 1; measured σR(z), σϕ(z), and σz(z) values are used). For the thick disc, other approximations are also shown. The crosses correspond to the case where σR(z), σϕ(z), and σz(z) are replaced with their mean values for |z| <  2hzt (0.27, 0.23, and 0.28, respectively; this improves agreement, probably because the values are less noisy). Diamonds indicate σz/σR = 0.7 and σϕ/σR = 0.7. Finally, open circles correspond to hR = 1, hv = 1. Frame (b) shows LOS integrated quantities separately for thin (blue curves; below) and thick discs (red; above). For the thin disc, the solid circles are calculated using the corrected σt(c) (see the text), while open circles use the non-corrected velocity dispersion, dominated by the large range of mean velocities along the LOS. In both cases, Eq. (15) is applied with the hR and hv values measured from the LOS-integrated radial profiles, while σz/σR = 0.7, σϕ/σR = 0.7, and k0 = 1 are assumed.

In the text

	Fig. 8. Examples of fitting the simulated LOS velocity (filled circles) and velocity dispersion profiles (open circles) when 0, 20, and 40% of thick-disc stars are in a retrograde component. The curves are the fitted values using the fT(z) ratio from simulated vertical profiles.
In the text

	Fig. 9. In black, fitted retrograde fraction for a range of input fractions; boxes indicate the input values. Continuous blue and dash-dot red curves are the fitted “effective” velocity dispersions for the thin- and thick-disc components, respectively.
In the text

Fig. 10. Various factors affecting the error between the actual input fraction of retrograde thick-disc stars in the N-body model, and the fraction returned by the fit. Upper left panel: fits with Eq. (20), using data from various axial distance intervals as indicated by the labels. Upper right panel: different variants of Eq. (20) are tested, using the data from the axial interval x/hR = 2.8 − 3.2: compared to the nominal case (“model=1”), in “model=2” the latter term in Eq. (20) has been omitted, while for “model=3” the has been replaced with , to correct for the contamination of LOS-velocity dispersions by systematic velocities. The middle panels (notice the change in the y-axis range) show the effect of using an incorrect fraction for thick disc, illustrated in the corresponding panel in the lowermost row. In the left panel frame the shape of fT is retained, but the contribution of the thick disc in the equatorial plane has been multiplied with a factor indicated by the labels. In the right panel frame, the fT profile has been stretched/dilated in the vertical direction by the amount indicated by the labels. The dashed lines indicated the fit obtained using the fT estimated from a photometric thin-thick disc decomposition to the simulation data.

In the text

	Fig. A.1. Maps and plots corresponding to ESO 157-49.
In the text

	Fig. A.2. As in Fig. A.1, but for ESO 443-21. Comerón et al. (2018) found no distinct thin- and thick-disc structure for this galaxy.
In the text

	Fig. A.3. As in Fig. A.1, but for ESO 469-15.
In the text

	Fig. A.4. As in Fig. A.1, but for ESO 544-27.
In the text

	Fig. A.5. As in Fig. A.1, but for IC 217.
In the text

	Fig. A.6. As in Fig. A.1, but for IC 1553. The small galaxy at the lower-left of the top-left panel is unrelated to IC 1553 and is found to have redshift z = 0.05.
In the text

	Fig. A.7. As in Fig. A.1, but for PGC 28308. In this galaxy there is no height where the thin disc has a lower surface brightness than the thick disc once the PSF is accounted for (Comerón et al. 2018).
In the text

	Fig. A.8. As in Fig. A.1, but for PGC 30591. Comerón et al. (2018) found no distinct thin- and thick-disc structure for this galaxy.
In the text