I start by providing an updated summary of the penalized pixelfitting PPXF method that is used to extract the stellar and gas kinematics as well as the stellar population of galaxies via full spectrum fitting. I then focus on the problem of extracting the kinematics when the velocity dispersion  is smaller than the velocity sampling V that is generally by design close to the instrumental dispersion SUBinstSUB. The standard approach consists of convolving templates with a discretized kernel while fitting for its parameters. This is obviously very inaccurate when   V2 due to undersampling. Oversampling can prevent this but it has drawbacks. Here I present a more accurate and efficient alternative. It avoids the evaluation of the undersampled kernel and instead directly computes its wellsampled analytic Fourier transform for use with the convolution theorem. A simple analytic transform exists when the kernel is described by the popular GaussHermite parametrization which includes the Gaussian as special case for the lineofsight velocity distribution. I describe how this idea was implemented in a significant upgrade to the publicly available PPXF software. The key advantage of the new approach is that it provides accurate velocities regardless of . This is important e.g. for spectroscopic surveys targeting galaxies with   SUBinstSUB for galaxy redshift determinations or for measuring lineofsight velocities of individual stars. The proposed method could also be used to fix Gaussian convolution algorithms used in todays popular software packages.
['techniques: radial velocities', 'techniques: spectroscopic', 'galaxies: kinematics and dynamics', 'Astrophysics - Astrophysics of Galaxies']
Improving the full spectrum fitting method accurate convolution with GaussHermite functions
{'Michele Cappellari ?': 'Sub-Department of Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH  \nAccepted 2016 November 18. Received 2016 November 18; in original form 2016 July 28', 'ABSTRACT': "I start by providing an updated summary of the penalized pixel-fitting ( ppxf ) method, which is used to extract the stellar and gas kinematics, as well as the stellar population of galaxies, via full spectrum fitting. I then focus on the problem of extracting the kinematic when the velocity dispersion GLYPH<27> is smaller than the velocity sampling GLYPH<1> V , which is generally, by design, close to the instrumental dispersion GLYPH<27> inst. The standard approach consists of convolving templates with a discretized kernel, while fitting for its parameters. This is obviously very inaccurate when GLYPH<27> < GLYPH<24> GLYPH<1> V = 2, due to undersampling. Oversampling can prevent this, but it has drawbacks. Here I present a more accurate and e GLYPH<14> cient alternative. It avoids the evaluation of the under-sampled kernel, and instead directly computes its well-sampled analytic Fourier transform, for use with the convolution theorem. A simple analytic transform exists when the kernel is described by the popular Gauss-Hermite parametrization (which includes the Gaussian as special case) for the line-of-sight velocity distribution. I describe how this idea was implemented in a significant upgrade to the publicly available ppxf software. The key advantage of the new approach is that it provides accurate velocities regardless of GLYPH<27> . This is important e.g. for spectroscopic surveys targeting galaxies with GLYPH<27> GLYPH<28> GLYPH<27> inst, for galaxy redshift determinations, or for measuring line-of-sight velocities of individual stars. The proposed method could also be used to fix Gaussian convolution algorithms used in today's popular software packages.  \nKey words: techniques: radial velocities - techniques: spectroscopic - galaxies: kinematics and dynamics", '1 INTRODUCTION': "The kinematics of the stars and gas in galaxies is a key ingredient in our understanding of how they form and evolve. Nowadays, this information is routinely extracted from integral-field spectroscopic (IFS) data, which provide a three dimensional view of galaxies. The technology is now available on all major telescopes and has become the standard way of obtaining spectra for galaxy evolution studies (see Cappellari 2016, for a review).  \nThe IFS data provide a 'fossil record' of galaxy formation. Specifically, the IFS kinematics of the stars allow one to distinguish a galaxy assembly dominated by gas poor merging from a growth driven by gas accretion and star formation (e.g. Emsellem et al. 2011; Cappellari et al. 2013b; Naab et al. 2014). Dynamical models based on IFS kinematics allow one to infer galaxy mass distributions to study scaling relations (e.g. Cappellari et al. 2013a; Scott et al. 2015), the stellar and dark matter content (e.g. Cappellari et al. 2012, 2015), or measure black hole masses (e.g. Krajnovi'c et al. 2009; Seth et al. 2014; Walsh et al. 2016). The gas content and kinematics tell us about the role of gas accretion in galaxy assembly (e.g. Sarzi et al. 2006; Davis et al. 2011; Barrera-Ballesteros et al.  \n2015) or the mechanism that regulates star formation in galaxies (e.g. Alatalo et al. 2011; Cheung et al. 2016; Ho et al. 2016).  \nFor more than a decade, this information has been extracted from IFS data by surveys targeting one galaxy at a time (SAURON de Zeeuw et al. 2002; DiskMass Bershady et al. 2010; ATLAS 3D Cappellari et al. 2011; CALIFA Sánchez et al. 2012). But the observational panorama is undergoing a revolution, with the arrival of large multiplexed IFS surveys targeting 10-20 galaxies at a time. Thousands of galaxies have already been observed in this manner (MaNGA Bundy et al. 2015; SAMI Bryant et al. 2015).  \nThese ongoing IFS surveys depend critically on the full spectrum fitting technique to deliver their science. The method is used to extract stellar and gas kinematics, as well as stellar population from the spectral data cubes. But, given the large number of objects and the variety of galaxy morphological types, these surveys are pushing the existing techniques to their limits.  \nThis paper is motivated by the existence of these ongoing large surveys and mostly arises from initial experiences with the analysis of the MaNGA data. In fact, one new characteristic of both the MaNGA and SAMI survey is that they are observing large numbers of galaxies (especially spirals) with stellar velocity dispersion GLYPH<27> well below the instrumental dispersion GLYPH<27> inst of the spectrographs. This situation was until now not very common, as the  \nobservers generally tried to target galaxies with GLYPH<27> inst < GLYPH<24> GLYPH<27> . However this restriction is actually not necessary to obtain useful kinematic information. In fact, although it is true that GLYPH<27> becomes intrinsically di GLYPH<14> cult to measure reliably when GLYPH<27> GLYPH<28> GLYPH<27> inst, the kinematics of lowGLYPH<27> galaxies is dominated by the stellar velocity V (e.g. Cappellari 2016), which remains a well defined observable. I realized that in this observational regime, all the spectrum fitting approaches, which are nearly universally adopted to extract stellar and gas kinematics, su GLYPH<11> er from limitations and can be significantly improved. Here in Section 2 I describe general concepts about kinematic extraction, in Section 3 I give an updated summary of the ppxf method, in Section 4 I discuss problems of the current approach and propose a clean solution, I summarize my paper in Section 5.", '2 GENERAL CONCEPTS': 'For external galaxies, or other stellar systems not resolved into individual stars, the stellar kinematics information is extracted by comparing the observed spectra, at di GLYPH<11> erent spatial locations, with broadened stellar templates. This comparison is nowadays almost invariably performed in pixel space rather than in Fourier space (e.g. Rix & White 1992; van der Marel 1994; Kelson et al. 2000; Cappellari & Emsellem 2004). This is because it makes the masking of bad pixels, or the simultaneous extraction of gas kinematics and stellar population, an easy and accurate process, while also allowing for other observational e GLYPH<11> ects to be easily modelled.  \nIn this section I review some general concepts, useful for understanding how the kinematics is generally extracted from galaxy spectra. I also try to clarify aspects that I discovered can lead to potential mistakes, from more than a decade of feedback I received from users of my publicly-available Penalized Pixel-Fitting software ppxf (Cappellari & Emsellem 2004).', '2.1 Templates observed like the galaxy spectrum': 'More than a decade ago, the stellar templates were generally observed with the same spectrograph used to obtain the galaxy spectrum under analysis (e.g. van der Marel et al. 1994). In this way one could assume that all the di GLYPH<11> erence between the galaxy spectrum and the templates was due to the broadening by the galaxy lineof-sight velocity distribution (LOSVD). In this idealized case one could schematically model the galaxy spectrum G mod as follows  \nG mod( x ) = T ( x ) GLYPH<3> L ( cx ) ; (1)  \nwhere x = ln GLYPH<21> is the natural logarithm of wavelength, T ( x ) is the stellar template, L is the LOSVD, with c the speed of light and GLYPH<3> denoting convolution.  \nIt is worth noting that the use of a logarithmic axis, although universally adopted in the past few decades, is not strictly required for the kinematic extraction. It ensures a constant velocity scale for every x . This means that L ( v ) remains constant along the x coordinate and allows one to use e GLYPH<14> cient Fast Fourier Transform (FFT) methods when performing the convolution.  \nIn principle however, it would also be possible to perform the convolution in linear GLYPH<21> . For this one would need to allow for the L ( v ) to vary with x , but one could still construct a relatively e GLYPH<14> -cient convolution, via direct summation, by making use of the fact that the convolution kernel L ( v ) generally consists of a few orders of magnitude fewer elements than the galaxy and templates spectra.  \nHowever, given that at least one resampling is unavoidable when reducing the spectral data, there are currently no compelling reasons to stop using logarithmic coordinates.', '2.2 Templates from spectral libraries': 'In the past decade, thanks to the availability of high resolution empirical stellar libraries spanning large ranges of atmospheric parameters over a wide wavelength range (e.g. Cenarro et al. 2001; Prugniel & Soubiran 2001; Valdes et al. 2004; Sánchez-Blázquez et al. 2006), libraries of synthetic stellar spectra (e.g. Coelho et al. 2005; Munari et al. 2005; Gustafsson et al. 2008), and stellar population models based on these libraries (e.g. Bruzual & Charlot 2003; Vazdekis et al. 2010; Maraston & Strömbäck 2011; Conroy & van Dokkum 2012), it has become common practice to extract the stellar kinematics by employing as templates linear combinations of a large number (up to a thousand) of template spectra observed with a di GLYPH<11> erent spectrograph than the one used to obtain the galaxy spectra.  \nDue to this di GLYPH<11> erence, before the kinematic extraction, and assuming the templates have a better resolution that the galaxy spectrum, the templates need to be matched to the galaxy spectrum by convolving them with a kernel K ( x ; GLYPH<21> ), which generally varies with wavelength GLYPH<21> , defined by the relation  \nLSFinst( x ; GLYPH<21> ) = LSFtemp( x ; GLYPH<21> ) GLYPH<3> K ( x ; GLYPH<21> ) : (2)  \nHere LSFinst( x ; GLYPH<21> ) is the line-spread function of the spectrograph used to observe the galaxy under study, which describes the observed shape of an intrinsically very narrow spectral line (in principle a GLYPH<14> function) at the wavelength GLYPH<21> , due to purely instrumental broadening, including the e GLYPH<11> ect of integrating over the detector pixels. LSFtemp( x ; GLYPH<21> ) is the LSF of the spectrograph used to observe the templates. After the templates have been convolved as  \n˜ T ( x ) = T ( x ) GLYPH<3> K ( x ; GLYPH<21> ) ; (3)  \nthey become as if they had been observed with the same spectrograph as the galaxy, and one can still model the galaxy spectrum with equation (1). Note that equation (3) is not a standard convolution, due to the GLYPH<21> dependence in the kernel. In signal processing this is called general convolution or general superposition integral.  \nIn principle the kernel K ( x ; GLYPH<21> ) could be measured in the very same way one uses to extract the LOSVD from galaxies, using solar templates to fit spectra of the twilight sky, while allowing for a non-Gaussian shape for K ( x ; GLYPH<21> ). However, in practice, it is quite di GLYPH<14> cult to accurately characterize the spectral, spatial and temporal variation of both the LSFinst and LSFtemp in the data. Moreover, the LSFs turn out to be generally well represented by Gaussian functions and, for this reason, are generally assumed to be Gaussians.  \nIt is well know that the convolution of two normalized Gaussians with dispersion GLYPH<27> a and GLYPH<27> b is another normalized Gaussian with dispersion GLYPH<27> 2 c = GLYPH<27> 2 a + GLYPH<27> 2 b (e.g. Brigham 1974, §2). For these reasons, when the LSF are assumed Gaussian, the kernel in equation (3) is a normalized Gaussian  \nG di GLYPH<11> ( x ; GLYPH<21> ) = exp n GLYPH<0> x 2 = h 2 GLYPH<27> 2 di GLYPH<11> ( GLYPH<21> ) io GLYPH<27> di GLYPH<11> ( GLYPH<21> ) p 2 GLYPH<25> (4)  \nwith dispersion  \nGLYPH<27> 2 di GLYPH<11> ( GLYPH<21> ) = GLYPH<27> 2 inst ( GLYPH<21> ) GLYPH<0> GLYPH<27> 2 temp ( GLYPH<21> ) ; (5)  \nwhere GLYPH<27> inst and GLYPH<27> temp are the Gaussian dispersion of the LSF for the instrument used to observe the galaxy and for the one used to obtain the stellar templates, respectively.  \nGiven the GLYPH<21> dependence of GLYPH<27> di GLYPH<11> , the general convolution in equation (3) needs to be performed via direct summation rather than using FFTs. But this only need to be done once, before fitting the LOSVD and for this reason its computation is not time critical. However, special care needs to be taken to avoid under-sampling during the computation of equation (3), as it can be a GLYPH<11> ected by the same under-sampling issues discussed in Section 4.2.', '2.3 From measured velocity to observed redshift': 'One important aspect that often causes confusion to users of spectral fitting software, is the connection between the measured lineof-sight velocity V and the galaxy redshift z . For this reason I clarify this aspect here.  \nAs described in Section 2.1, spectral fitting codes rebin the spectra logarithmically and measure velocities from the shift GLYPH<1> x required to match the spectral templates to the galaxy spectrum. This shift is precisely related to the redshift as follows  \nGLYPH<1> x = GLYPH<1> ln GLYPH<21> = ln( GLYPH<21> obsv =GLYPH<21> emit) = ln(1 + z ) : (6)  \nHowever, spectral fitting codes also need to define a velocity scale to be used as variable in the LOSVD (e.g. to quantify the galaxy velocity dispersion). For small velocity di GLYPH<11> erences V GLYPH<17> GLYPH<1> V , the Doppler formula is  \nV GLYPH<25> c GLYPH<1> GLYPH<21> GLYPH<21> emit GLYPH<25> c GLYPH<1> ln GLYPH<21>: (7)  \nSo one can define a constant velocity scale per pixel, which reduces to the Doppler formula in the limit of small V , namely on the scale of the LOSVD  \nV GLYPH<17> c GLYPH<1> x = c GLYPH<1> ln GLYPH<21> = c ln(1 + z ) : (8)  \nThis definition is used to convert the pixels di GLYPH<11> erences GLYPH<1> x measured by the codes into velocities, to provide the V and GLYPH<27> in km s GLYPH<0> 1 . Given this definition, one has to use equation (8) to convert back the measured velocities into redshifts, when needed. In particular, the well known approximation z GLYPH<25> V = c should never be used, being already less accurate than the typical measurement uncertainties when z > GLYPH<24> 0 : 01.', '2.4 Separating peculiar velocities and cosmological redshift': "In the common case of IFS data, one is interested in the stellar and gas velocities V bin measured within di GLYPH<11> erent spatial bins with respect to the galaxy barycentre (e.g. to construct dynamical models, or measure gas outflows) and not in the actual redshift z bin of the individual spectra. When the redshift is negligible ( z < GLYPH<24> 0 : 01), the V bin can be computed with su GLYPH<14> cient accuracy by simply subtracting the velocity V syst of the galaxy barycentre from the measured velocities V directly provided by the spectral fitting program. However this simple approach will lead to dramatic errors at larger redshift, where a quite di GLYPH<11> erent approach must be used.  \nIn cosmology, the V bin are called 'peculiar velocities', to distinguish them from the cosmological recession velocity or redshift z cosm. The latter is formally defined as the redshift of the galaxy barycentre. For galaxies at significant redshift, an accurate and simple approach to measure V bin consists of bringing all the galaxy spectra onto the rest frame wavelength GLYPH<21> rest = GLYPH<21> obs = (1 + z 0 cosm ), where GLYPH<21> obs is the observed wavelength and z 0 cosm is an initial estimate of the galaxy redshift. This division is actually equivalent to merely re-defining the origin of the velocity coordinate, while keeping the  \nspectra unchanged. Crucially, the galaxy LSFinst( x ; GLYPH<21> ) must be compressed by the same factor or, when this is assumed Gaussian, its instrumental dispersion in units of wavelength GLYPH<27> rest inst = GLYPH<27> obs inst = (1 + z 0 cosm ). When quantified as R = GLYPH<1> GLYPH<21>=GLYPH<21> , in km s GLYPH<0> 1 , or as GLYPH<1> ln GLYPH<21> , the instrumental resolution is independent of redshift, for every spectral pixel, but the pixels wavelength changes. This implies that, in general, the matching of the templates resolution of equation (5) must still be performed for every di GLYPH<11> erent redshift.  \nOnce these two steps have been performed, the velocities V returned by the spectral fitting program will be the desired V bin. From the extracted rest-frame velocity field, one can accurately measure possible residual o GLYPH<11> sets V bary of the barycentre velocity from zero, by enforcing symmetries in the field (e.g. Krajnovi'c et al. 2006, appendix C). This small o GLYPH<11> set can just be subtracted from all V bin. Finally, an improved estimate of the galaxy redshift can be obtained using the relation (e.g. Cappellari et al. 2009, eq. 2)  \n1 + z cosm = (1 + z 0 cosm ) GLYPH<2> (1 + V bary = c ) : (9)  \nwhich derives form the general expression linking the cosmological redshift z cosm, the bins peculiar velocities V bin and the observed bins redshift z bin (e.g. Hogg 1999, eq. 10)  \n1 + z bin = (1 + z cosm) GLYPH<2> (1 + V bin = c ) : (10)", '3 OVERVIEW OF THE PPXF METHOD': 'One particular technique used to extract stellar kinematics is called the Penalized Pixel-Fitting ( ppxf ) method. It was originally described in Cappellari & Emsellem (2004), but it has significantly evolved over more than a decade of intense usage, to address our own specific needs and requests or feedback from users. A number of the new features in ppxf were not described in detail in the literature, this is true in particular for the use of the covariance matrix (Section 3.3), the implementation of regularization, for stellar population studies (Section 3.5), the fitting of gas emission lines (Section 3.6) and the kinematic bulge-disk decomposition (Section 3.7). For this reason, both the new and old features are summarized in this section. I describe here the current version 6.0 of the Python implementation of the software 1 . But most of the features are also implemented in the IDL version of ppxf .', '3.1 Model for the galaxy spectrum': 'The ppxf method approximates the observed galaxy spectrum via the following quite general parametrization  \nG mod( x ) = N X n = 1 wn 8 > > < > > : [ Tn ( x ) GLYPH<3> L n ( cx )] K X k = 1 ak P k ( x ) 9 > > = > > ; + L X l = 0 bl P l ( x ) + J X j = 1 cjS j ( x ) ; (11)  \nwhere the L n are the LOSVDs, which can be di GLYPH<11> erent for the N templates Tn (e.g. Johnston et al. 2013), the P k and P l are multiplicative or additive orthogonal polynomials of degree k and l respectively (of Legendre type or a truncated Fourier series), and S j are spectra of the sky (e.g. Weijmans et al. 2009). The additive polynomials can significantly minimize template mismatch by changing the strength of individual absorption lines. Mismatch can be non', '4 M. Cappellari': "negligible even when using today's large stellar libraries as templates. Additive polynomials can additional correct for imperfect sky subtraction or scattered light. The multiplicative polynomials can correct for inaccuracies in the spectral calibration, and make the fit insensitive to reddening by dust, without the need to adopt a specific reddening curve. Note that the 0th degree multiplicative polynomial is not included. However, if the calibration is reliable and one is interested in estimating the reddening from the shape of the stellar continuum, the multiplicative polynomials P k ( x ) can be replaced by the expression  \nf ( x ) = 10 GLYPH<0> 0 : 4 E ( B GLYPH<0> V ) k ( x ) (12)  \nwhere k ( x ) is a specific reddening curve (e.g. Cardelli et al. 1989; Calzetti et al. 2000) and E ( B GLYPH<0> V ) quantifies the amount of reddening, which is constrained by ppxf to be positive. Any of the polynomials or sky components is optional and does not need to be actually used during a ppxf fit.  \nThere is considerable freedom in the choice of templates depending on the specific application. The templates can consist of combinations of e.g. (i) stellar population models with specific parameters (e.g. Bruzual & Charlot 2003; Vazdekis et al. 2010; Maraston & Strömbäck 2011; Conroy & van Dokkum 2012); (ii) individual empirical (e.g. Cenarro et al. 2001; Prugniel & Soubiran 2001; Valdes et al. 2004; Sánchez-Blázquez et al. 2006) or synthetic (e.g. Coelho et al. 2005; Munari et al. 2005; Gustafsson et al. 2008) stars; (iii) principal components derived from a library of spectra; (iv) weighted sum of di GLYPH<11> erent stars; (v) gas emission lines; or (vi) sky spectra. The first option is used when ppxf is employed for 'full spectrum' fitting to study the galaxies stellar population (e.g. Onodera et al. 2012; Cappellari et al. 2012; Morelli et al. 2013, 2015; McDermid et al. 2015; Shetty & Cappellari 2015), but it can also be used when extracting the stellar kinematics with ppxf (e.g. Emsellem et al. 2004; Oh et al. 2011; Thomas et al. 2013), while the remaining choices are useful when extracting the stellar or gas kinematics with ppxf , due to the extra freedom which allows one e.g. to fit possible variations in elemental abundances (e.g. Cappellari et al. 2011; Blanc et al. 2013; van de Sande et al. 2016).", '3.2 Parametrization for the LOSVD': "The LOSVD L n ( v ), for both the stellar and gas templates, are parametrized using the Gauss-Hermite parametrization introduced for this purpose by van der Marel & Franx (1993) and Gerhard (1993). Note however that the two papers did not quite define the same parametrization for the LOSVD. The crucial di GLYPH<11> erence is that, while Gerhard (1993) fitted all moments in the Gauss-Hermite expansion, van der Marel & Franx (1993) chose to explicitly set the first three coe GLYPH<14> cients to ( h 0 ; h 1 ; h 2) = (1 ; 0 ; 0), and only fit the higher coe GLYPH<14> cients, in such a way that the LOSVD has the form  \nL ( y ) = exp GLYPH<16> GLYPH<0> y 2 = 2 GLYPH<17> GLYPH<27> p 2 GLYPH<25> 2 6 6 6 6 6 4 1 + M X m = 3 hmHm ( y ) 3 7 7 7 7 7 5 ; (13) y = ( v GLYPH<0> V ) =GLYPH<27>  \nIt was the latter form, of equation (13), which has become the current standard in the field (e.g. Bender et al. 1994) and was also adopted by ppxf . Here, Hm are the Hermite polynomials, standardized (by definition) in such a way that the terms of equation (13) are quantum-mechanical wave functions for the one-dimensional harmonic oscillator (e.g. Schi GLYPH<11> 1968, § 4)  \nHm ( y ) = H Abr m ( y ) p m ! 2 m ; (14)  \nwere the H Abr m are the so-called 'physicists' Hermite polynomials, defined as in Abramowitz & Stegun (1964, eq. 22.2.14). The H Abr m polynomials are the ones most commonly provided by default by today's popular software. They are e.g. the form returned by M athematica 's (Wolfram 2003) function HermiteH , or by N umpy 's (Oliphant 2007) polynomial.hermite.hermval , or by MATLAB's function hermiteH . The first two Hermite polynomials Hm in equation (13) are (e.g. van der Marel & Franx 1993, eq. A5)  \nH 3 = y GLYPH<16> 2 y 2 GLYPH<0> 3 GLYPH<17> p 3 ; H 4 = 4 GLYPH<16> y 2 GLYPH<0> 3 GLYPH<17> y 2 + 3 p 24 : (15)", '3.3 Linear fitting procedure': 'For every choice for the possible non-linear parameters in the model of equation (11)  \n[ V ; GLYPH<27>; h 3 ; : : : ; hM ; a 1 ; : : : ; aK ; E ( B GLYPH<0> V )] ; (16)  \nppxf minimizes the functional (e.g. Press et al. 2007, §19.5)  \nF = GLYPH<31> 2 + GLYPH<21> B (17)  \nwhere the first term measures the agreement between the model spectrum G mod and the observed galaxy spectrum G , while the second is an adjustable term, which can be zero, and which quantifies the smoothness of the weights, in a space spanned by the population parameters (e.g. age, metallicity, GLYPH<11> enhancement or the stellar initial mass function [IMF]). The computation of the first term of equation (17) is described in this section, while the ppxf implementation of the second one is described in Section 3.5.  \nIn the general case in which one knows the covariance matrix GLYPH<6> , with elements cov( xj ; xk ), were xj and xk are a pair of spectral pixels, the agreement between the data and the model, making the standard assumption of Gaussian uncertainties, is quantified by  \nGLYPH<31> 2 = GLYPH<2> A GLYPH<1> x GLYPH<0> y GLYPH<3> T GLYPH<1> GLYPH<6> GLYPH<0> 1 GLYPH<1> GLYPH<2> A GLYPH<1> x GLYPH<0> y GLYPH<3> = GLYPH<13> GLYPH<13> GLYPH<13> L GLYPH<0> 1 GLYPH<1> ( A GLYPH<1> x GLYPH<0> y ) GLYPH<13> GLYPH<13> GLYPH<13> 2 (18) = k r k 2 ;  \nwhere kGLYPH<1>k is the Euclidean norm, the columns of the matrix A consist of the convolved templates (multiplied by the polynomials), the additive polynomials and the sky spectra, y = G is the galaxy spectrum,  \nx = ( w 1 ; : : : ; wN ; b 0 ; : : : ; bL ; c 1 ; : : : ; cJ ) (19)  \np  \nis the solution vector, L = GLYPH<6> is the square-root or the positivedefinite covariance matrix GLYPH<6> = L GLYPH<1> L T , computed via the Cholesky decomposition (e.g. Press et al. 2007, §2.9), and r is the weighted vector of residuals.  \nThe minimum GLYPH<31> 2 of equation (18) is found by solving the following special quadratic programming problem  \nminimize f ( x ) = GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<16> L GLYPH<0> 1 GLYPH<1> A GLYPH<17> GLYPH<1> x GLYPH<0> L GLYPH<0> 1 GLYPH<1> y GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> 2 subject to wn GLYPH<21> 0 ; n = 1 ; : : : ; N (20) cj GLYPH<21> 0 ; j = 1 ; : : : ; J :  \nSpecific, e GLYPH<14> cient and robust algorithms exist to solve this type of problems, with guaranteed convergence in a finite number of steps to the global minimum. Currently the IDL version of ppxf uses the Bounded-Variables Least Squares ( bvls ) algorithm by Lawson & Hanson (1974), while the Python version uses the Non-Negative Least Squares ( nnls ) code by the same authors, made available by the Scipy (Jones et al. 2001) procedure optimize.nnls . In the  \nnnls solution, I employ slack variables to remove the positivity constraints from the additive polynomials coe GLYPH<14> cients ( b 0 ; : : : ; bL ).  \nIn the common situation where spectral covariance is ignored, or unknown, and one only has the error spectrum, then  \nL GLYPH<0> 1 = diag [1 = GLYPH<1> G ( x 1) ; : : : ; 1 = GLYPH<1> G ( xP )] (21)  \nreduces to a diagonal matrix, where GLYPH<1> G ( xp ) is the 1 GLYPH<27> uncertainty of every pixel xp in the galaxy spectrum. In this special, but usual case, the computation can be simplified, as there is no need to decompose and multiply by the inverse covariance matrix, and the residual vector r has elements  \nrp = G mod( xp ) GLYPH<0> G ( xp ) GLYPH<1> G ( xp ) ; p = 1 ; : : : ; P : (22)', '3.4 Non-linear fitting procedure': "A key feature of ppxf , from which the method derives its name, is the fact that it automatically penalizes non-Gaussian solutions, to reduce the noise in the recovered kinematics, when the data do not contain su GLYPH<14> cient information to constrain the full shape of the LOSVD. This is done by minimizing a new objective function  \nGLYPH<31> 2 p = GLYPH<31> 2 + GLYPH<11> P ; (23)  \nwhere the GLYPH<31> 2 is given in equation (18) or equation (22), P is a penalty function, which describes the deviation of the LOSVD from a Gaussian shape, while GLYPH<11> is an adjustable penalty, which depends on the data quality. The penalty is implemented in an e GLYPH<14> cient way, by perturbing the residual vector r in such a way that the sum-ofsquares nature of the problem is preserved. This is described in the original ppxf paper (Cappellari & Emsellem 2004, section 3.3) and will not be repeated here.  \nIn both cases, with or without the covariance matrix, the minimization of the GLYPH<31> 2 p as a function of the non-linear parameters in equation (16) is performed with the specific Levenberg-Marquardt non-linear least-squares optimization algorithm (e.g. Press et al. 2007, §15.5.2). This exploits the sum-of-squares nature of the optimization problem by requiring the user to provide the vector of residuals r , instead of the scalar GLYPH<31> 2 of equation (18), to compute explicitly the Hessian matrix of the GLYPH<31> 2 merit function. This makes the method much more robust and e GLYPH<14> cient than generic optimizers of scalar functions.  \nppxf uses the state-of-the-art minpack implementation of the Levenberg-Marquardt method by Moré et al. (1980), which is based on the robust trust region method. It was converted into IDL and named mpfit by Markwardt (2009), with the important addition of the ability to set bounds or keep any of the variables fixed. mpfit was ported to Python by Mark Rivers and later adapted for use with Numpy by Sergey Koposov.  \nIn ppxf , every individual non-linear parameter of the LOSVD, for every kinematic component separately, can be: (i) left free; (ii) kept fixed at a specified value; or (iii) bound within specified values. This flexibility has numerous applications. As an example, one can fit for the kinematic (e.g. V and GLYPH<27> ), templates and polynomials for the gas emission-line components in individual spaxels, while keeping the stellar kinematics (e.g. V , GLYPH<27> , h 3 and h 4) fixed to the value determined from a larger spatial bin.  \nA possible alternative for the non-linear optimization under Python would be the Scipy function optimize.least\\_squares , introduced in version 0.18 in 2016, which is also designed for leastsquares problems with bounds and is also based on the trust region method. I performed some tests using ppxf with least \\_ squares for  \nthe non-linear optimization, trying to match as closely as possible both the accuracy requirement and the step size for the numerical derivatives, that I currently use in ppxf with mpfit . I extracted the stellar and gas kinematics both from simulated spectra, with a range of GLYPH<27> , as well as from real MaNGA data cubes. I found that ppxf using least \\_ squares requires an execution time and a number of function evaluations which is generally comparable, but never smaller than the public ppxf version using mp -fit . This is true both when using the trust-region reflective algorithm ( method='trf'; Branch et al. 1999 ) and the dogleg algorithm ( method='dogbox' , Voglis & Lagaris 2004, Nocedal & Wright 2006, § 4) in least \\_ squares . It was reassuring to see both least \\_ squares methods invariably converge to the same mpfit solution within the specified accuracy. These tests indicates there is currently no compelling reason to replace mpfit in ppxf .", '3.5 Weights regularization for stellar population': 'To study the star formation history and stellar population of a galaxy, ppxf uses the common approach of modelling its unobscured rest-frame spectrum by discretizing the following integral equation (e.g. Cid Fernandes et al. 2005; Ocvirk et al. 2006; Tojeiro et al. 2007)  \nG mod( GLYPH<21> ) = Z t = T t = 0 SSP GLYPH<21> ( t ; Z ) GLYPH<1> SFR( T GLYPH<0> t ) d t ; (24)  \nwhere SFR is the star formation rate, SSP GLYPH<21> is a Single Stellar Population spectrum per unit mass, with age t and metallicity Z , while T is the age of the Universe at the redshift of the galaxy. This expression is straightforwardly generalized in ppxf to study the distribution of more parameters, like e.g. metallicity, GLYPH<11> enhancement or IMF, in addition to the SFR (e.g. Conroy 2013, §2.3).  \nThe equation (24) is an inhomogeneous Fredholm equation of the first kind, with kernel SSP GLYPH<21> . And the recovery of the SFR( t ) from the observed G mod is a textbook example of ill-conditioned inverse problem (e.g. Press et al. 2007, §19.0). This means that the recovery su GLYPH<11> ers from severe degeneracies and a unique solution cannot be found without further assumptions.  \nA standard way of dealing with this type of problem is by using regularization (e.g. Tikhonov & Arsenin 1977; Hansen 1998), which can be thought of as damping the high-frequency variations in the solution, unless they are actually required to describe the data, or finding a trade-o GLYPH<11> between the quality of the fit and the noise in the solution (e.g. Press et al. 2007, §19.4). In other words, regularization allows one to select the smoothest solution, among the many degenerate solutions that are equally consistent with the data.  \nA conceptual mistake that is often made regarding regularization, is to think that this is equivalent to assuming the solution has to be smooth. This is not correct. In fact, the solution can be as non-smooth as required by the data and e.g. it will allow for multiple sharp bursts of star formation if the data require them. The solution will only be smooth if this is consistent with the data.  \nUnder some simplifying assumptions, the regularized solution has a simple Bayesian interpretation: it represents the most likely solution for the weights, given an adjustable prior on the amplitude of the fluctuations (e.g. Press et al. 2007, §19.4.1). An alternative to regularization consists of explicitly exploring the posterior of the allowed solutions using a Markov Chain Monte Carlo approach (e.g. Gelman et al. 2013, § 11), while using the individual weights as non-linear model variables. Although potentially interesting and  \nconceptually simple, we found the latter approach quite time consuming for general usage and did not yet include it in the public version of ppxf .  \nAssuming, without loss of generality, that one is performing a spectral-fitting study of the stellar population, while varying a single parameter (e.g. age), the smoothness functional which appears in equation (17) is defined as follows (e.g. Press et al. 2007, §19.5)  \nGLYPH<21> B = GLYPH<21> Z w 00 ( t ) 2 d t / N GLYPH<0> 1 X n = 2 ( wn GLYPH<0> 1 GLYPH<0> 2 wn + wn + 1) 2 GLYPH<1> (25)  \nwhere the wn are the spectral weights and the functional is zero only for a linear function of the weights. This expression is straightforwardly extended within ppxf to dimensions large than one, by performing the finite di GLYPH<11> erences separately along every axis. In this case the software requires the stellar population templates to form a full grid in 2-dim or 3-dim, in such a way that one can easily map the elements of the weights vector onto the 2-dim or 3-dim coordinates in population parameters.  \nIn practice, the regularization is implemented in ppxf by augmenting the matrix in equation (20) with one extra row, namely one extra equation  \nwn GLYPH<0> 1 GLYPH<0> 2 wn + wn + 1 GLYPH<1> = 0 ; (26)  \nfor every term in the summation of equation (25), in such a way that the regularized solution of the problem in equation (20) becomes  \nminimize f ( x ) = GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> L GLYPH<0> 1 GLYPH<1> A GLYPH<21> C ! GLYPH<1> x GLYPH<0> L GLYPH<0> 1 GLYPH<1> y 0 ! GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> 2 subject to wn GLYPH<21> 0 ; n = 1 ; : : : ; N (27) cj GLYPH<21> 0 ; j = 1 ; : : : ; J ;  \nwith C the array with the coe GLYPH<14> cients of equation (26). This augmented system is still solved by bvls or nnls as in the nonregularized case. Also here, the computation is specialized and simplified when the covariance matrix is not given. Applications of ppxf with regularization to the study of the star formation history or stellar population in galaxies were presented e.g. in Onodera et al. (2012), Cappellari et al. (2012), Morelli et al. (2013, 2015) and McDermid et al. (2015).  \nA conceptually similar regularization approach was used in the stecmap full-spectrum fitting code (Ocvirk et al. 2006). However the stecmap implementation is very di GLYPH<11> erent. Their approach only regularizes in one dimension. Moreover it does not exploit the quadratic nature of the sub-problem of equation (27), but instead solves it as a generic non-linear problem, using a variablemetric non-linear optimization method. This makes their approach less robust as well as much slower and more complex than the nnls approach used by ppxf . Fast and accurate convergence of the solution of the quadratic sub-problem of equation (27) is essential for ppxf to be able to reliably solve for the non-linear parameters of the LOSVD together with the template weights.', '3.6 Modelling gas emission lines': 'The ability of ppxf to fit di GLYPH<11> erent kinematics, or a di GLYPH<11> erent LOSVD L n , for every individual template also allows one to fit the gas emission lines together with the stellar kinematics and stellar population. This is achieved with ppxf by simply passing a set of gas emission-lines templates, together with the stellar ones. In this way, ppxf can fit for the kinematics and weights of both the stellar and  \ngas components, without the need for the code to know, or make any algorithmic distinction, between the two sets of templates. However, given that the gas kinematics is often dominated by rotation rather than by random motions, the gas emission lines can have quite low dispersion. This implies that an accurate treatment of undersampling, which is the focus of this paper, is critical for unbiased gas measurements.  \nThe generality of the ppxf approach, where emission lines are just normal templates, gives flexibility to the users. One can freely specify the emission lines templates to include in the fit, and each line can be described by an arbitrary number of components if desired. The kinematics of every line can be tied to each other, or they can be fitted independently. And the relative fluxes of emission line doublets can be fixed, by placing two lines in the same template. The LOSVD of the gas emission lines are not restricted to being Gaussians, but they can be described by the same general GaussHermite parametrization of equation (13), when needed, and in this case one can use the same penalty P of equation (23). Application of ppxf to the simultaneous fit of the gas emission lines and the stellar population where presented e.g. in Johnston et al. (2013), Shetty &Cappellari (2015) and Mitzkus et al. (2016).  \nIt is worth mentioning that a simultaneous fit of both the stellar and gas kinematics should not be used as the standard approach. In general, more robust stellar kinematics is obtained by masking the gas emissions and including additive polynomials to reduce template mismatch. And the stellar kinematics is generally extracted from larger spatial bins than those needed for the more clumpy gas. Moreover, for extracting the gas fluxes and kinematics, as well as for the stellar population, one generally only includes multiplicative polynomials, to prevent changes in the line strength of the absorption features in the templates. For these reasons, in general, the gas kinematics is extracted from a separate ppxf fit, at fixed stellar kinematics. The simultaneous fitting of gas and stellar templates is especially useful when studying stellar population, as one does not need to mask, or clean, regions of possible emission. This allows one to make full use e.g. of the spectral region of the Balmer series, where gas emission is common, and which provides strong constraints on the SFR( t ).  \nIn practice, the gas templates should represent the LSFinst( x ; GLYPH<21> ) at the wavelength of a given set of emission lines. Usually the LSF is approximated by a Gaussian, and for maximum accuracy, the gas templates need to be integrated over the spectral pixels, as is automatically the case for the stellar templates. Pixels integration is generally ignored by other software, however it has a small but measurable e GLYPH<11> ect on the recovered gas dispersion as I will show in Section 4.4. For a normalized (unit area) Gaussian LSF, the pixelintegrated gas template, for a single line of wavelength GLYPH<21> line, in the same logarithmic coordinates xp as the stellar templates, is  \nT gas( xp ) = 1 2 ( erf " xp GLYPH<0> x 0 + ( GLYPH<1> x )pixel = 2 GLYPH<27> inst( x 0) p 2 # GLYPH<0> erf " xp GLYPH<0> x 0 GLYPH<0> ( GLYPH<1> x )pixel = 2 GLYPH<27> inst( x 0) p 2 #) ; p = 1 ; : : : ; P (28)  \nwhere x 0 = ln GLYPH<21> line, ( GLYPH<1> x )pixel = ( GLYPH<1> ln GLYPH<21> )pixel is the size of one spectral pixel and erf is the error function (Abramowitz & Stegun 1964, eq. 7.1.1). The gas template is band-limited as it is the convolution of a box function and a Gaussian with dispersion GLYPH<27> inst. When the template is Nyquist sampled at intervals ( GLYPH<1> x )pixel GLYPH<25> GLYPH<27> inst, it is uniquely defined by its samples and will result in accurate convolutions with the LOSVD to near machine accuracy. The Nyquist sampling is guaranteed by the fact that, if the galaxy spectra or tem-  \nplates were sampled at steps larger than GLYPH<27> inst, one would be wasting useful information in the data!  \nThe ppxf approach is closely related, but conceptually di GLYPH<11> erent from the one adopted by the popular and state-of-the-art G an -dalf gas-fitting code (Sarzi et al. 2006), which was in fact derived from an earlier version of ppxf , and first introduced the simultaneous fitting of multiple stellar and gas components, for the accurate subtraction of the stellar continuum. The G andalf code treats the gas emission lines and the stellar templates in a separate manner: the non-convolved Gaussian emission lines are added to the stellar templates, which are convolved with a fixed stellar kinematics. Ultimately both the ppxf and G andalf approaches are expected to give comparable results for the gas emissions, when adopting similar settings. The key di GLYPH<11> erence is that ppxf was designed to extract stellar kinematics, for multiple components, and regularized stellar population, in addition to the gas kinematics.', '3.7 Kinematic bulge disk decomposition': 'A recent addition to ppxf is the possibility of constraining the ratio of the first two kinematic components to a desired value. This can be useful e.g. to perform kinematic bulge-disk decompositions (Tabor et al. 2016). One can force the template spectra for the bulge component to contribute a prescribed fraction f bulge of the total flux in the fitted spectrum  \nf bulge = P w bulge P w bulge + P w disk ; (29)  \nwhere w bulge and w disk are the weights assigned to the set of spectral templates used to fit the bulge and disk respectively.  \nThis constraint is enforced in the same way as the regularization constraints of equation (27), by adding the following extra equation to the matrix of equation (20)  \n( f bulge GLYPH<0> 1) P w bulge + f bulge P w disk GLYPH<10> = 0 ; (30)  \nwith the parameter GLYPH<10> set to a very small number (e.g. GLYPH<10> = 10 GLYPH<0> 9 ), which specifies the relative accuracy at which this equation needs to be satisfied. When both the spectral templates and the galaxy spectrum are normalized to have a mean flux of order unity, the best fitting weights w bulge and w disk are generally smaller than unity and equation (30) is satisfied as an equality constraint to numerical accuracy.', '4 EXTRACTING KINEMATICS AT LOW DISPERSION': 'In this section I discuss the problem of extracting accurate kinematics using full spectrum fitting, in the special case in which the velocity dispersion is smaller than the velocity step ( GLYPH<1> V )pixel. This is generally chosen to be similar to the instrumental dispersion GLYPH<27> inst. In fact, if ( GLYPH<1> V )pixel < GLYPH<27> inst the pixels have significant covariance and some data are redundant, while if ( GLYPH<1> V )pixel > GLYPH<27> inst the spectrum is not Nyquist sampled and one is wasting useful information. In what follows I assume ( GLYPH<1> V )pixel = GLYPH<27> inst, and sometimes use the two interchangeably, when referring to the critical scale where issues appear.  \nI illustrate the problems and solution discussed in this section using numerical experiments. My goal is not to asses the general ability of ppxf to recover galaxy kinematics, or other practical issues, like the biases as a function of the quality of the stellar templates. This has already been covered in the literature (e.g. Westfall  \net al. 2011). Instead, I focus here exclusively on the important numerical problems due to the discretization of the LOSVD.', '4.1 Setting up the numerical experiments': "For the experiments I restrict myself to the wavelength range 360750 nm, which is a spectral region containing a number of useful absorption features for the kinematic extraction. This wavelength range is fully covered e.g. by the MaNGA (Bundy et al. 2015) and CALIFA (Sánchez et al. 2012) integral-field spectroscopic surveys. This is also the region covered by the spectra of the MILES stellar library (Sánchez-Blázquez et al. 2006). However the results I present are general and are quite insensitive to the adopted wavelength region.  \nI use as my 'galaxy' spectrum the solar spectrum 2 by Kurucz (2005), which has a resolution R = GLYPH<1> GLYPH<21>=GLYPH<21> GLYPH<25> 300 000 FWHM and a typical S = N GLYPH<25> 3000 per spectral element. This spectrum can be regarded as having essentially infinite resolution and S = N for all practical purposes of my tests. Also in this case, I have verified that the results of the experiments are general and weakly depend on the specific choice of the input spectrum, as long as it is representative of real high-resolution galaxy spectra. Nearly indistinguishable results were obtained, in all my tests, using e.g. a high-resolution R = 20 000, solar metallicity, 10 Gyr, single stellar population model of Maraston & Strömbäck (2011), based on the MARCS theoretical stellar library (Gustafsson et al. 2008).  \nTo test the recovery of the gas kinematics, on top of the Solar spectrum, I add a realistic set of Gaussian gas emission lines. I include the Balmer series H GLYPH<11> , H GLYPH<12> , H GLYPH<13> , H GLYPH<14> , with relative fluxes following an unreddened Balmer decrement (Storey & Hummer 1995). I additionally include the [O iii ] GLYPH<21>GLYPH<21> 4959,5007 and [N ii ] GLYPH<21>GLYPH<21> 6548,6583 doublets, with fluxes in the standard ratio 1 / 3, and the [S ii ] GLYPH<21>GLYPH<21> 6716,6731 doublet with flux ratio of one. The doublets are scaled with respect to the Balmer lines according to the ratios [O iii ] GLYPH<21> 5007 = H GLYPH<12> = 2, [N ii ] GLYPH<21> 6583 = H GLYPH<11> = 1 = 2 and [S ii ] GLYPH<21>GLYPH<21> 6716 ; 6731 = H GLYPH<11> = 1 = 2. These ratios are nearly irrelevant for the tests, but are chosen to be representative of realistic ionization levels for what Ho et al. (1997) calls 'transition objects' between star-forming regions and LINERS. The Balmer series and the three doublets are all treated as separate kinematic components, in addition to the stellar one. This means that the ppxf fits presented here have five kinematic components.  \nI adopt an idealized spectrograph, with a Gaussian LSF, with a constant instrumental dispersion GLYPH<27> inst = 70 km s GLYPH<0> 1 (i.e. R GLYPH<25> 1800 FWHM). This is the median resolution of the MaNGA spectrograph, but is also similar to the one used by other galaxy surveys: it is similar e.g. to the resolution of the ATLAS 3D survey (Cappellari et al. 2011), to the high-resolution mode used by the CALIFA survey (Sánchez et al. 2012), and to the low-resolution mode used by the SAMI survey (Bryant et al. 2015). I assume the spectra are Nyquist sampled by the detector, and consequently adopt a spectral pixel ( GLYPH<1> x )pixel = GLYPH<27> inst = 70 km s GLYPH<0> 1 . This is the same velocity sampling adopted e.g. by the SDSS data release 13 (SDSS Collaboration 2016).  \nI investigate the case where I use as input stellar template in ppxf the same spectrum used for the galaxy spectrum. I further assume a quite high S = N = 200 per spectral interval ( GLYPH<1> x )pixel. Both choices are made again to provide a clean experiment which isolates the under-sampling problem from e.g. the unrelated issue of  \nFigure 1. Solar spectrum with added gas emission lines. The figure illustrates a typical ppxf fit to the solar spectrum by Kurucz (2005), and to the gas emissions, used for our tests. The black line (mostly hidden by the fit) is the relative flux of the observed spectrum (with noise added at S = N = 200), for an adopted dispersion GLYPH<27> in = 140 km s GLYPH<0> 1 , for both the gas and the stars. The red line is the ppxf fit for the stellar component, while the orange line is a fit to the gas emission lines. The green symbols at the bottom are the fit residuals, while the blue lines is the gas-only best-fitting spectrum. The main absorption and emission features are indicated at the top of the plot. The two inset plots show an enlarged view of the ppxf fits and residuals in the regions including the H GLYPH<11> (top) and H GLYPH<12> (bottom) lines.  \n<!-- image -->  \ntemplate mismatch. I adopt the default degree=4 additive polynomials, but my results are totally insensitive to this choice. A representative ppxf fit to the solar spectrum and to the added emission lines, is shown in Fig. 1.  \nIn my experiments, the solar spectrum was initially logarithmically sampled with velocity scale of 2 km s GLYPH<0> 1 per pixel, which is an integer factor smaller than the final detector pixels. Then the gas emission lines were added. Subsequently, the whole spectrum was accurately convolved with a very well-sampled, discretized LOSVD and with a Gaussian LSF. Then the spectrum was integrated over the 70 km s GLYPH<0> 1 wide pixels, by summing every adjacent set of 35 pixels. Finally noise was added to every pixel. In the tests, the input velocity was chosen randomly for every realization, to prevent the LOSVD from being aligned in a constant manner with respect to the pixels boundaries. The same applies to the starting guess for V and GLYPH<27> in ppxf , which were chosen randomly for every Monte Carlo realization.  \nI emphasize the fact that for the initial LOSVD convolution I did not use the analytic Fourier transform introduced in this paper, nor the analytic pixel integration of the gas emission lines. Instead I performed these steps numerically on the finely-sampled spectrum. This provides a useful debugging of my software implementation, for both gas and stars. In fact a good recovery can only be achieved if the analytic Fourier approach correctly and ac-  \ncurately corresponds to the numeric pixel integration and LOSVD convolution, in the limit of a well-sampled kernel, and when the spectrum contains information on the Gauss-Hermite coe GLYPH<14> cients, namely when GLYPH<27> > GLYPH<24> GLYPH<27> inst.", '4.2 Description of the problem': "The convolution of the templates with the LOSVDs in equation (11) was until now performed by ppxf according to the standard discrete definition  \n( T GLYPH<3> L ) p GLYPH<17> Q = 2 X q = GLYPH<0> Q = 2 Tq + p L q ; (31)  \nwhere Q is the number of elements where the kernel L q GLYPH<17> L ( cxq ) is non-zero and Tq GLYPH<17> T ( xq ). In practice, for computational e GLYPH<14> ciency, this convolution is performed using the standard, mathematically equivalent, Fourier approach (e.g. Press et al. 2007, §13.1)  \nT GLYPH<3> L = F GLYPH<0> 1 [ F ( T ) GLYPH<1> F ( L )] (32)  \nwhere F is the Discrete Fourier Transform (DFT) and F GLYPH<0> 1 is its inverse, which can be computed e GLYPH<14> ciently with a number of operations proportional to P log P using the classic Fast Fourier Transform (FFT) algorithm (Cooley & Tukey 1965). For e GLYPH<14> ciency, ppxf pre-computes the F ( T ) of all the templates, as they do not change  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Stellar kinematics recovery with a Gaussian LOSVD. The three panels illustrate the ppxf recovery of the velocity V and dispersion GLYPH<27> for an adopted Gaussian LOSVD with known velocity V in and dispersion GLYPH<27> in. The vertical dashed line indicates the instrumental dispersion GLYPH<27> inst = 70 km s GLYPH<0> 1 , which coincides with the adopted velocity sampling ( GLYPH<1> x )pix. The grey region with GLYPH<27> < GLYPH<27> inst = 10 is not allowed by the program. The thin lines in the bottom panels illustrate the e GLYPH<11> ect on GLYPH<27> of an error in the instrumental dispersion GLYPH<1> GLYPH<27> inst = GLYPH<6> 3%. Left Panel: recovery without oversampling, with the old ppxf . Note the dramatic increase of the errors below GLYPH<27> in < GLYPH<24> GLYPH<27> inst = 2. Middle Panel: recovery with a well-sampled LOSVD and a template oversampled by a factor of 10, with the old ppxf . This approach is slower and produces a significant drop in the dispersion below GLYPH<27> in < GLYPH<24> GLYPH<27> inst = 2. Right Panel: recovery with the solution proposed in this paper, namely using an analytic Fourier transform of the LOSVD, with the new ppxf . This is at least as fast as the left panel, but here both V and GLYPH<27> are recovered without bias.  \n<!-- image -->  \nduring the GLYPH<31> 2 minimization, and uses the specific FFT for real input, Numpy's fft.rfft , which is based on fftpack (Swarztrauber 1982), to further decrease the FFT computation time by a factor of two.  \nFor the discrete convolution, the LOSVD, which acts as convolution kernel, was until now sampled by ppxf at discrete intervals ( GLYPH<1> x )pixel = ( GLYPH<1> V )pixel = c before the computation of F ( L ). This straightforward approach of discretely sampling the kernel, represents the standard practice, when one needs to perform convolutions of spectra or images with Gaussian kernels. It is used e.g. in the GAUSS\\_SMOOTH function in IDL, in the ndimage.filters.gaussian\\_filter function of the Python package Scipy (Jones et al. 2001), or in the convolution.Gaussian1DKernel function of the Python package A stropy (Astropy Collaboration 2013).  \nHowever, it is clear that, when GLYPH<27> < GLYPH<24> ( GLYPH<1> x )pixel, where GLYPH<27> is the dispersion of the Gaussian, the kernel starts becoming severely under-sampled and cannot be expected to accurately represent the LOSVD any more. As an illustration, when the Gaussian is centred on a pixel and GLYPH<27> = ( GLYPH<1> x )pixel = 2, the kernel has essentially only three non-zero elements L = (0 : 14 ; 1 ; 0 : 14), fully dominated by the middle one. I note here that integrating the LOSVD over the pixels instead of sampling it at the pixel centre, as done for a Gaussian in equation (28), is not a correct solution to the undersampling problem. Instead, it is equivalent to incorrectly performing an extra convolution of the templates with a box function of the size of a pixel.  \nOne should not expect a reliable recovery of the LOSVD at these low dispersion, and in particular, both the velocity and the dispersion should not be expected to be determined to better than about half a pixel. This was the reason why, until now, ppxf included the oversample keyword to analytically calculate a wellsampled LOSVD before convolving it with a template oversampled by interpolation onto an equally-fine grid. This approach however su GLYPH<11> ers from two limitations: (i) it requires a significant increase in the size of the template spectra, resulting in an increase of the  \ncomputation time and (ii) the oversampling of the templates is necessarily smooth below the observed scale, and cannot represent the dense forest of thin absorption lines in real stellar spectra. In fact the smooth oversampling has a similar e GLYPH<11> ect as an extra convolution of true spectrum with a box function of the size of one spectral pixel.  \nThe left panel of Fig. 2 illustrates the dramatic problems which appear when GLYPH<27> in < GLYPH<24> ( GLYPH<1> x )pixel = 2, while adopting a Gaussian input LOSVD: ppxf becomes essentially unable to reliably recover both the velocity and the dispersion of the LOSVD. In fact, although the velocity has an rms error of only 0.3 km s GLYPH<0> 1 , when GLYPH<27> in > ( GLYPH<1> x )pixel, the errors become as high as 20 km s GLYPH<0> 1 at low GLYPH<27> in. This problem is not specific to ppxf but will a GLYPH<11> ect any spectral fitting program that uses the standard equation (31) to define convolution!  \nThe middle panel of Fig. 2 shows that, as expected, by computing a well-sampled LOSVD and convolving it with an oversampled template, the velocity can be well recovered, however the dispersion tends to be under-estimated at low GLYPH<27> in, because the oversampled template is smoother than the real spectrum. In the commonsituation adopted for my tests, where ( GLYPH<1> x )pixel GLYPH<25> GLYPH<27> inst, the oversampling has a similar e GLYPH<11> ect as a 3% over-estimation of GLYPH<27> inst.  \nSimilar problems can be seen when the LOSVD is assumed to be described by a significantly non-Gaussian shape, parametrized by equation (13). I adopted as input some realistic values h 3 = h 4 = 0 : 1 (e.g. Emsellem et al. 2004). Note that, contrary to what we did in Cappellari & Emsellem (2004), here I assume the input LOSVD to be precisely described by the adopted parametrization, because I want to test the accurate recovery of the known input Gauss-Hermite moments.  \nThe left panel of Fig. 3 again illustrates the dramatic problems in the recovery of the velocity and dispersion, when GLYPH<27> in < GLYPH<24> ( GLYPH<1> x )pixel = 2. In this case, with non-Gaussian LOSVD, one can see the expected convergence towards zero of the h 3 and h 4 parameters, due to the penalty term in equation (23). For all the examples in Fig. 3 I adopted the same penalty bias=1 in ppxf . This becomes important for GLYPH<27> in < GLYPH<24> ( GLYPH<1> x )pixel, where the broadening by the instru-  \nmental dispersion makes the deviations from a Gaussian hard to measure and ppxf tries to penalize the LOSVD towards a Gaussian. But at the smallest GLYPH<27> in even the penalty looses it e GLYPH<11> ectiveness, and one can see sharp variations, especially in the recovered h 4, instead of the desired convergence toward zero (i.e. Gaussian shape).  \nThe middle panel of Fig. 3 again shows that by computing a well-sampled LOSVD and oversampling the template one can overcome the main problems with the recovery of the velocity. However the velocity dispersion shows the same behaviour as in the middle panel of Fig. 2, with the recovered GLYPH<27> being underestimated, and hitting the lower boundary, when GLYPH<27> in < GLYPH<24> 20 km s GLYPH<0> 1 .", '4.3 Solution of the problem': "A clean solution to the under-sampling problem comes by noting that the Fourier transform of the kernel happens to be analytic in the special case of interest, where the LOSVD is parametrized by the Gauss-Hermite series of equation (13). This suggests one can use an analytic rather than discrete Fourier transform to evaluate F ( L ) in equation (32).  \nIf I define the continuous Fourier transform of a time-domain function h ( t ), using the modern-physics convention, with ! the angular frequency, as  \nH ( ! ) = 1 p 2 GLYPH<25> Z 1 GLYPH<0>1 h ( t ) e i ! t d t ; (33)  \nthe H ( ! ) of the LOSVD has an exceedingly simple form. This is because the individual terms of L are eigenfunctions of the H ( ! ) operator. This fact was already pointed out by van der Marel & Franx (1993, eq. A7), who quote the result from Gradshteyn & Ryzhik (1980, eq. 7.376.1), and of course it can be verified analytically. It is a textbook result (e.g. Schi GLYPH<11> 1968, §4), related to the fact that the terms of equation (13) are quantum-mechanical wave functions for the harmonic oscillator. These have the same form in either position or momentum-space, which are related via the Fourier transform. Specifically, if one defines  \nL m ( t ) = exp GLYPH<16> GLYPH<0> t 2 = 2 GLYPH<17> GLYPH<27> p 2 GLYPH<25> Hm ( t ) (34)  \nthe following simple transform pair holds  \nh ( t ) = L m ( t ) () H ( ! ) = i m L m ( ! ) : (35)  \nIn other words, the H ( ! ) have the same form as the original L m functions, simply transformed into the frequency domain, while the even / odd terms are alternating real / imaginary respectively. This latter fact comes from two general properties of H ( ! ): (i) if h ( t ) is real and even, then the same is true for H ( ! ); (ii) if h ( t ) is real and odd, then H ( ! ) is imaginary and odd.  \nThe analytic calculation of F ( L ) has two major advantages, with respect to its discrete version. They derive from general properties of the Fourier Transform (e.g. Brigham 1974, §3):  \nGood kernel sampling: When the LOSVD is narrow, ad consequently poorly sampled, in the 'time' or spatial domain, it will become broad in the frequency domain according to the transform pair  \n1 j GLYPH<27> j h GLYPH<18> t GLYPH<27> GLYPH<19> () H ( GLYPH<27>! ) frequency scaling : (36)  \nThis implies that the LOSVD will never be under-sampled in the frequency domain because the LOSVD is necessarily non-zero only for a small fraction of the size of the full template spectrum.  \nSimple velocity shifting: When evaluating the kernel of equation (13) for the discrete convolution of equation (31), one needed to sample a symmetric range of velocities from GLYPH<0> v max to v max, fully enclosing the o GLYPH<11> set LOSVD. This produced many zero elements in the kernel, causing practical issues with large V . In Fourier space this problem disappears as one can use  \nh ( t GLYPH<0> t 0) () H ( ! ) e i ! t 0 time shifting (37)  \nto shift the velocity of the templates by simply changing the phase of all elements of the complex transform.", '4.4 Final formula and tests': 'Combining equation (35) to equation (37), I find that the analytic Fourier transform of the LOSVD of equation (13) is  \nH ( ! ) = exp GLYPH<16> i ! V GLYPH<0> GLYPH<27> 2 ! 2 = 2 GLYPH<17> p 2 GLYPH<25> 2 6 6 6 6 6 4 1 + M X m = 3 i m hmHm ( GLYPH<27>! ) 3 7 7 7 7 7 5 : (38)  \nThis expression is used to compute the DFT F ( L ) which appears in equation (32). Given that L is a real function, it is uniquely defined over half of the angular frequencies  \nh ( t ) is real = ) H ( GLYPH<0> ! ) = [ H ( ! )] GLYPH<3> : (39)  \nThis implies that one only needs to evaluate H ( ! ) of equation (38) in the interval ! = [0 ; GLYPH<25> ] to obtain the DFT of real input. This needs to have the same number of elements as the DFT of real input F ( T ), namely half the size of T .  \nThe right panel of Fig. 2 illustrates the dramatic improvement in the recovery of both the velocity and the dispersion, when using the analytic Fourier transform approach. No biases can be detected in the recovery, and importantly, the scatter in the residuals is nearly independent of GLYPH<27> in. In particular, the velocity recovery is as accurate as in the case of the oversampled templates, down to the smallest GLYPH<27> in. Additionally, even the very small systematic trend in velocity, visible in the middle panel of Fig. 2, now completely disappears. The GLYPH<27> is also recovered without detectable bias, and in particular, the underestimation of GLYPH<27> that a GLYPH<11> ects the oversampled case, also disappears.  \nSimilar conclusions apply to the case with non-Gaussian LOSVDin the right panel of Fig. 3. Also in this case, ppxf using the analytic F ( L ) looks similar to the oversampled case, but now the GLYPH<27> underestimation disappears. Note that the penalty of ppxf is still important to prevent wildly degenerate solutions and consequently large noise in all parameters when GLYPH<27> in < GLYPH<24> GLYPH<27> inst. This important aspect was described in detail in Cappellari & Emsellem (2004). It is still relevant here, but the discussion will not be repeated.  \nIn Fig. 4 I illustrate the accuracy in the recovery of the gas kinematics using ppxf with the analytic Fourier transform, and in particular here for the case of the [O iii ] GLYPH<21>GLYPH<21> 4959,5007 doublet. The results are the same for the other emission lines. The left panel shows the e GLYPH<11> ect of neglecting the pixel integration of the Gaussian lines. In this case the bias in GLYPH<27> is almost equivalent to an underestimation of GLYPH<27> inst of over 3%, due to the missing convolution by a box function of the size of one spectral pixel. As one may have expected, the e GLYPH<11> ect on GLYPH<27> is nearly the opposite as that of oversampling the galaxy spectrum, in the middle panel of Fig. 2. The right panel of Fig. 4 presents the ppxf recovery when including the analytic pixel integral of equation (28). The gas kinematics is unbiased and recovered as accurately as the stellar one, at any GLYPH<27> .  \nAlso shown, with thin lines, in the bottom panels of Fig. 2 and Fig. 4 is the variation in the recovered GLYPH<27> one should expect when the  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Stellar kinematics recovery with LOSVD described by the Gauss-Hermite parametrization. The three panels illustrate the ppxf recovery of the kinematics for a LOSVD with adopted velocity V in, dispersion GLYPH<27> in and realistic h 3 = h 4 = 0 : 1. The vertical dashed line indicates the instrumental dispersion GLYPH<27> inst = 70 km s GLYPH<0> 1 , which coincides with the adopted velocity sampling ( GLYPH<1> x )pix. The grey region with GLYPH<27> < GLYPH<27> inst = 10 is not allowed by the program. Left Panel: recovery without oversampling, with the old ppxf . Note the dramatic increase of the errors below GLYPH<27> in < GLYPH<24> GLYPH<27> inst = 2. Middle Panel: recovery with a well-sampled LOSVD and a template oversampled by a factor of 10, with the old ppxf . This approach is slower and produces a significant drop in the dispersion below GLYPH<27> in < GLYPH<24> GLYPH<27> inst = 2. Right Panel: recovery with the solution proposed in this paper, namely using an analytic Fourier transform of the LOSVD, with the new ppxf . This is at least as fast as the left panel, but here all kinematic parameters are properly recovered. Note that h 3 and h 4 start converging towards zero (i.e. Gaussian LOSVD) when GLYPH<27> in < GLYPH<24> GLYPH<27> inst. This is intentional and unavoidable. It is used to prevent a dramatic increase of the errors in V and GLYPH<27> when the parameters of the LOSVD become degenerate at low GLYPH<27> in. This e GLYPH<11> ect is extensively discussed in Cappellari & Emsellem (2004).  \n<!-- image -->  \n<!-- image -->  \nFigure 4. Gas kinematics with ppxf . The two panels illustrate the recovery of the gas kinematics for the [O iii ] GLYPH<21>GLYPH<21> 4959,5007 doublet shown in Fig. 1. The input lines are assumed Gaussian, with adopted velocity V in and dispersion GLYPH<27> in, and are fitted with Gaussian gas templates. The vertical dashed line indicates the instrumental dispersion GLYPH<27> inst = 70 km s GLYPH<0> 1 , which coincides with the adopted velocity sampling ( GLYPH<1> x )pix. The thin lines in the bottom panels illustrate the e GLYPH<11> ect on GLYPH<27> of an error in the instrumental dispersion GLYPH<1> GLYPH<27> inst = GLYPH<6> 3%. Left Panel: recovery of gas kinematics while ignoring the integration of the Gaussian line over the spectral pixels. This is nearly equivalent to an underestimation of the GLYPH<27> inst. Right Panel: gas recovery with the proper integration of the Gaussian line over the pixels of equation (28). The gas kinematics is now accurately recovered down to the smallest dispersion.  \n<!-- image -->  \nassumed instrumental dispersion is in error by GLYPH<1> GLYPH<27> inst = GLYPH<6> 3%, which represents the best accuracy achievable today (e.g. Yan et al. 2016, § 7.4). The curves illustrate the fact that, at GLYPH<27> < GLYPH<24> GLYPH<27> inst = 2, biases in GLYPH<27> will be dominated by the uncertainty in the LSF, which becomes a fundamental barrier to the accurate recovery of GLYPH<27> in this regime. For this reason, it is important to emphasize that the key advantage of the proposed ppxf upgrade is not the unbiased GLYPH<27> , but rather the fact that one can measure accurate velocities at any GLYPH<27> . The velocity is empirically an accurately defined quantity regardless of GLYPH<27> and actually becomes even better defined at low GLYPH<27> . It is crucial for a kinematic extraction method to be able to keep providing reliable velocities even when GLYPH<27> is too low to be measurable.  \nA complementary way of dealing with the undersampling of the LOSVD consists of using stellar templates observed with higher resolution that the galaxy, while keeping them at their native resolution. For this reason ppxf was modified to allow for templates with a smaller velocity scale than the galaxy. It is also possible to skip the broadening of the templates to the galaxy LSF described in Section 2.2 and instead include the Gaussian convolution in equation (38) for H ( w ), using the associative property of the convolution. An additional convolution by a Gaussian is obtained by multiplying H ( w ) by the Gaussian Fourier transform. This implies that a convolution with the Gaussian of equation (4) is obtained by replacing the term inside the exponential in equation (38) with  \ni ! V GLYPH<0> GLYPH<16> GLYPH<27> 2 di GLYPH<11> + GLYPH<27> 2 GLYPH<17> ! 2 = 2 : (40)  \nTo use this expression GLYPH<27> di GLYPH<11> must be constant in km s GLYPH<0> 1 . This must be achieved by homogenizing the template resolution in such a way that it has a constant GLYPH<27> di GLYPH<11> at all wavelengths.', '4.5 Relation to previous work': 'I have not found previous usage in the literature of an analytic Fourier transform to improve the accuracy of the spectral fitting method. However the general idea is obviously a good one and for this reason it was used before in the literature, in a similar context, for the special case of Gaussian functions.  \nThe first generation of methods to extract the galaxies stellar kinematics generally performed the GLYPH<31> 2 minimization in Fourier space, rather than in pixels space as ppxf does. Given that the data to fit consisted of a Fourier transform, and it is common practice to use analytic fitting functions, it was natural to use as fitting function an analytic Fourier transform of the LOSVD, when available. In fact, this general idea, specifically the use of an analytic Fourier transform of a Gaussian broadening function, was employed by Sargent et al. (1977) and Schechter & Gunn (1979) for their improved version of the Fourier quotient method. The analytic transform was fitted in Fourier space to the ratio of the FFT of the galaxy and template. An analytic Fourier transform of the Gaussian was also used by Tonry & Davis (1979) in their analysis of the errors of the cross-correlation method.  \nI found applications of an analytic Fourier transform for convolution in image processing. For example, the algorithm I describe in Section 4.4 for the Gauss-Hermite functions is a generalization of the algorithm 5 in the survey of Gaussian convolution algorithms by Getreuer (2013). A key di GLYPH<11> erence in this paper is that the templates we want to convolve are bandwidth limited (by the Gaussian LSF) and Nyquist sampled (at steps GLYPH<27> inst). Due to the sampling theorem, this implies the templates are completely determined by their samples. This is the reason we obtain an unbiased recovery of GLYPH<27> down to the smallest values in the right panels of Fig. 2 and Fig. 4,  \nand we are not a GLYPH<11> ected by the convolution bias at small GLYPH<27> described by Getreuer (2013), for generic non band-limited functions.  \nIn another broadly related work in image processing, Berger & Simental (1999) noted the fact that two-dimensional GaussHermite functions can be convolved analytically with Gaussians, thanks to their analytic Fourier transform, and propose to use this fact to accurately deconvolve images, to remove PSF e GLYPH<11> ects.', '5 SUMMARY': 'In the first part of the paper, I provided an overview, or tutorial, of general concepts useful to understand and properly interpret the extraction of kinematics from galaxy spectra. I tried to clarify in particular the questions I received more often, over more than a decade, from users of my publicly available ppxf software.  \nThen I gave an updated overview of the ppxf method. I concentrated especially on the description of features of the method which were included after the publication of the original paper in 2004, some of which had never been properly explained and precisely documented in the literature.  \nSubsequently, I focused on the problem of extracting kinematics via full spectrum fitting, when the velocity dispersion is smaller than the spectral sampling, which is generally chosen to be the same as the instrumental dispersion. I illustrated the obvious but dramatic problems that arise when one completely ignores the issue, as well as the limitations of the previous solution, which consist of oversampling the spectra.  \nFinally I provided a clean solution to the long-standing undersampling issue, which consist of using the analytic Fourier transform of the LOSVD in conjunction with the convolution theorem. This approach completely removes the need for oversampling and makes the full spectrum fitting method suitable for measuring reliable kinematics at any velocity dispersion. This is especially crucial for the mean velocity, which now becomes a well-determined quantity even when the dispersion becomes negligible, and consequently impossible to reliably recover from real data.  \nThe approach described in this paper was implemented in a significant upgrade to the publicly available ppxf code, and is already being used as part of the MaNGA Data Analysis Pipeline (Wesfall et al. in preparation). Further tests on real IFS data will be published elsewhere.  \nThe proposed solution appears quite natural, however, perhaps surprisingly, it is currently not being used by any of the popular software packages. Given the simplicity of our approach for accurate convolutions, we argue it should become standard practice.', 'ACKNOWLEDGEMENTS': 'This paper was motivated by productive discussions within the MaNGA Data Analysis Pipeline team and in particular with Kyle Wesfall and Matt Bershady. I thank Samantha Penny for starting the discussion on this subject. I am also grateful to Eric Emsellem and Daniel Thomas for comments. I acknowledge support from a Royal Society University Research Fellowship. This paper made use of M atplotlib (Hunter 2007) and of the lineid \\_ plot 3 Python program by Prasanth Nair.', '531, 159': 'Voglis C., Lagaris I., 2004, in WSEAS International Conference on', '14 M. Cappellari': 'Applied Mathematics. Corfu, Greece, http: // www.wseas.us / e-library / conferences / corfu2004 / papers / 488-317.pdf  \nWalsh J. L., van den Bosch R. C. E., Gebhardt K., Yıldırım A., Richstone  \nD. O., Gültekin K., Husemann B., 2016, ApJ, 817, 2  \nWeijmans A.-M., et al., 2009, MNRAS, 398, 561  \nWestfall K. B., Bershady M. A., Verheijen M. A. W., 2011, ApJS, 193, 21 Wolfram S., 2003, The Mathematica Book, 5th ed.. Wolfram Media, Cham-  \npaign, IL, http: // reference.wolfram.com /  \nYan R., et al., 2016, preprint, (  \narXiv:1607.08613  \n)  \nde Zeeuw P. T., et al., 2002, MNRAS, 329, 513  \nvan de Sande J., et al., 2016, submitted  \nvan der Marel R. P., 1994, MNRAS, 270, 271  \nvan der Marel R. P., Franx M., 1993, ApJ, 407, 525  \nvan der Marel R. P., Rix H. W., Carter D., Franx M., White S. D. M., de Zeeuw T., 1994, MNRAS, 268, 521  \nThis paper has been typeset from a T E X / L A T E X file prepared by the author.'}