Forward modeling of galaxy kinematics in slitless spectroscopy

Slitless spectroscopy has long been considered as a complicated and confused technique. Nonetheless, with the advent of Hubble Space Telescope (HST) instruments characterized by a low sky background level and a high spatial resolution (most notably WFC3), slitless spectroscopy has become an adopted survey tool to study galaxy evolution from space. We investigate its application to single object studies to measure not only redshift and integrated spectral features, but also spatially resolved quantities such as galaxy kinematics. We build a complete forward model to be quantitatively compared to actual slitless observations. This model depends on a simplified thin cold disk galaxy description -- including flux distribution, intrinsic spectrum and kinematic parameters -- and on the instrumental signature. It is used to improve redshifts and constrain basic rotation curve parameters, i.e. plateau velocity $v_{0}$ (in $\rm km.s^{-1}$) and central velocity gradient $w_{0}$ (in $\rm km.s^{-1}.arcsec^{-1}$). The model is tested on selected observations from 3D-HST and GLASS surveys, to estimate redshift and kinematic parameters on several galaxies measured with one or more roll angles. Our forward approach allows to mitigate the self-contamination effect, a primary drawback of slitless spectroscopy, and therefore has the potential to increase precision on redshifts. In a limited sample of well-resolved spiral galaxies from HST surveys, it is possible to significantly constrain galaxy rotation curve parameters. This proof-of-concept work is promising for future large slitless spectroscopic surveys such as EUCLID and WFIRST.


Introduction
Spectroscopy surveys play a fundamental role in the understanding of galaxy formation and evolution with cosmic time and in cosmology. These surveys have been achieved using different techniques. Fiber-fed multi-object spectrographs are commonly used to measure redshift and integrated properties on pre-selected targets. For instance, the Sloan Digital Sky Survey (SDSS) has measured more than a million of local galaxy spectra (York et al. 2000;Strauss et al. 2002), and has revealed how the star formation rates, metallicities, stellar populations vary with environment, mass and redshift (Gómez et al. 2003;Brinchmann et al. 2004;Kauffmann et al. 2004). Similar surveys like the 2dF and 6dF Galaxy Redshift Surveys (Folkes et al. 1999;Jones et al. 2004) have also allowed to constrain cosmological parameters by mapping the distribution of galaxies along cosmic time (Cole et al. 2005). More recently, the VIPERS survey (de la Torre et al. 2013) allowed to map with an unprecedented precision the large-scale distribution of galaxies by measuring more than 100 000 redshifts at 0.5 < z < 1.2. Nevertheless, these fiberbased surveys suffer from drawbacks: galaxy central regions are integrated so two-dimensional internal structure cannot be properly recovered, and objects need to be selected and targeted a priori.
Alternatively, Integral-Field Spectrographs (IFS) are ideal to study resolved objects over their spatial extent, but still require explicit pointing of individual galaxies. In the last years, the development of highly successful IFS surveys (R 1500)such as SAURON (de Zeeuw et al. 2002), ATLAS 3D , SAMI (Croom et al. 2012), CALIFA (Sánchez et al. 2012), and MaNGA (Bundy et al. 2015) -has pushed our understanding of galaxy properties further. Not limited to integrated measurements, these surveys in the nearby Universe (z 0.1) allowed to accurately map the gas, stellar populations and kinematics, and led to a new kinematical classification scheme of early-type galaxies (Emsellem et al. 2007. Contrary to inherently targeted multi-object spectrographs, "panoramic" IFS can be used to carry on untargeted surveys, but only on very limited sky areas because of observation time/cost constraints. E.g., MUSE (Bacon et al. 2010) has probed the evolution of gas kinematics of low-mass galaxies (M ≤ 10 10 M ) up to z = 1.4 on the Hubble Deep Field South (Contini et al. 2016).
In contrast to graceful IFS, slitless spectroscopy has generally been seen as a clumsy technique with well recognized drawbacks. As a matter of fact, the absence of an independent spatial sampling before spectral dispersion induces two major contamination issues: self-contamination -the effective spectral resolution is directly related to the size and shape of the spatially resolved object in the dispersion direction -and crosscontamination -signal pollution from nearby objects; both effects make the data reduction difficult and the redshift measurement less accurate. In addition, slitless spectroscopy is affected by a comparatively high integrated background level (particularly from the ground), limiting the depth or SNR of observations. However, this technique has some pros of its own: the ease of instrumental design and observational use, a large field of view and a very high multiplexing capability, all leading to very large object catalogs (∼ 100 000 galaxies in recent surveys and 15M in future surveys). Furthermore, it has the potential to Thereby, since the advent of the Hubble Space Telescope (HST) grism instruments, slitless spectroscopy has de facto become a tool of choice to study galaxy evolution from space: a low background level and a fine spatial resolution both mitigate the aforementioned shortcomings (Freudling et al. 2008;Kümmel et al. 2011;Dressel 2012). Dedicated HST surveys such as WISP (Atek et al. 2010), 3D-HST Momcheva et al. 2016), GLASS (Schmidt et al. 2014;Treu et al. 2015), FIGS (Pirzkal et al. 2017) have led to consider this technique as appropriate to derive redshift and integrated galaxy properties over large samples, and ready to be used in future missions as EUCLID (Grupp et al. 2012) and WFIRST (Spergel et al. 2015).
Traditional approaches in slitless spectroscopy (e.g. Kümmel et al. 2009) use standard "inverse" data reduction and analysis methods, extracting parameters from observations using successive and dedicated data manipulation steps. Typically, it involves an empirical modeling of the spectral trace, a cross-dispersion summation to estimate the 1D galaxy spectrum, the ad-hoc combination of spectra obtained at different position angles, and any subsequent spectral analyses performed on the averaged spectrum. Not only the proper error propagation is difficult between the different data-reduction steps, but such reverse approach can hardly correct for or quantify the impact of spatially resolved galactic properties, such as internal dynamics or metallicity gradients.
Alternatively, a "forward" approach allows to constrain physical or instrumental parameters directly in the observation space, properly accounting for degeneracy and covariances, and allowing for the inclusion of bayesian-like priors. By constructing a predictive model of the galaxy 2D dispersed image (hereafter coined spectrogram) depending on a set of observationally or physically motivated parameters, we investigate the possibility to measure not only intrinsic mean spectral quantities -e.g., redshift, emission line intensities and widths -independently of self-confusion, but also spatially resolved quantities such as internal kinematics. Thus, by combining forward methods to derive resolved quantities on single object and the large multiplexing power of multi-object spectrographs, slitless spectroscopy surveys offer a unique opportunity to study galaxy properties at an unprecedented scale.
In this paper, we will detail how to forward-model slitless spectrograms from a galaxy model -including flux distribution, intrinsic spectrum and kinematic parameters -and an instrumental signature. Considering our targets are mainly lineemitting disk galaxies, we will use two majors assumptions: an axi-symmetric thin cold disk geometry for the galaxy, and a separability hypothesis under which the intrinsic galaxy spectrum is supposed uniform over its whole extent. Using this approach, we will investigate the application of slitless spectroscopy to single object studies to measure internal kinematic parameters, namely the plateau velocity v 0 (Kalinova et al. 2017;Varidel et al. 2019) and the central velocity gradient (CVG) w 0 (Lelli et al. 2013;Erroz-Ferrer et al. 2016).
The article is organized as follows. In Sec. 2 we present galaxy kinematics in slitless spectroscopy and describe the model parametrization. We test our method on simulated spectrograms in Sec. 3, and apply it on selected galaxies from 3D-HST and GLASS survey in Sec. 4. We discuss the results in Sec. 5, and conclude and open some perspectives in Sec. 6.

Resolved kinematics in slitless spectroscopy
In this section, we investigate the kinematic signature in a slitless resolved galaxy spectrum. The internal velocity induces differential Doppler shifts (in addition to the systemic cosmological redshift), inducing small offsets of observed wavelengths as a function of position, and therefore distorting the overall spectral shape.
The slitless spectrogram I(x, y) can be derived from two key ingredients, first the spectro-spatial flux distribution cube of the galaxy C(r, λ), which contains all the observable information -spatial profile, intrinsic spectrum, velocity field, instrumental transmission and PSF, etc. -and second the 2D dispersion law D(λ) from the spectrograph, relating the wavelength to the (x, y)offset on the detector. More details on how the spectrogram is computed is given in Sec. 2.3.

Pedagogic case
To illustrate the effect of resolved kinematics in slitless spectroscopy, we build a pedagogic simulation mimicking the observation of an Hα-emitting disk galaxy at z ∼ 0.9 with an exponential thin disk density profile, with an inclination of 60 • and a scale length of r d = 6 px, and a dispersion direction perpendicular to the major axis; a uniform intrinsic spectrum made of a constant continuum and an Hα + [N ii] emission line complex at z = 0.9; a typical plateau velocity curve with v → v 0 = 300 km s −1 beyond transition radius r 0 = 10 px (see Sec. 2.3.4).
For illustration purpose, the simulated instrument is similar to an HST-like slitless spectrograph but with an unrealistic spectral sampling of D = 2.5 Å px −1 , ten times better than actual slitless instruments. Simulated slitless spectrograms are shown in Fig. 1. In both cases, one can distinguish the Hα + [N ii] emission line complex from the constant continuum which spreads out on each side. When the velocity field is ignored (upper panel), the emission lines have a shape similar to the flux distribution; on the opposite, when the kinematic effects are included (lower panel), the velocity-induced Doppler offset significantly distorts the emission line spectral shape. As can be seen, the kinematic signature on the slitless spectrogram is somewhat similar to the one traditionally observed in long slit spectroscopy.
Two important lessons can be learned from this simple simulation. In this particular model, the disc scale length r d has been chosen twice smaller than the turnover radius of the velocity curve r 0 . As a consequence, only the inner solid body rotation part of the velocity field has a clear observational signature, the plateau region being too far out the exponential disk extent to have any significant impact on spectrogram shape (see galaxy rotation curves plotted as red dashed lines in the lower panel of Fig. 1). This is discussed more in Sec. 3.
Secondly, it appears that the relative position angle (PA), defined as the angle between the (projected) galaxy major axis and the cross-dispersion direction, plays a critical role in the kinematic signature, PA = 0 (the major axis is perpendicular to the dispersion direction) being the most favorable case. In the central region of the galaxy where the kinematics is dominated by the solid body rotation, we can approximate the effective PA observed in the spectrogram (PA eff ) from the relative PA of the galaxy and the CVG w 0 (see Fig. 2 Fig. 1. Toy simulation of the spectrogram of a typical Hα-emitting disk galaxy at z ∼ 0.9 (intrinsic spectrum in uppermost panel, spatial profile on the leftmost panels) as observed with an HST-like slitless spectrograph but with an improved spectral resolution of R ∼ 2500 (see text). Top: without the kinematic effects; bottom: the signature of the intrinsic velocity field is clearly visible as a distortion of the spectrogram. The red dashed lines correspond to the galaxy Rotation Curves at each emission lines position and are tracing the distortion in the spectrogram.  where s is the spatial sampling of the instrument (in arcsec px −1 ), w 0 is expressed in km s −1 arcsec −1 , and kinematic sampling R kin is defined below in Eq. (2). Since relative PA and CVG w 0 are so correlated, they cannot be constrained independently from the spectrogram alone: the kinematic major axis needs to be set a priori from external photometric observations (more discussions in Sec. 2.3.4).
Overall, due to self confusion effects inherent to slitless spectroscopy, the kinematic parameters are entangled with spectrospatial flux distribution of the galaxy. This is the goal of our analysis to estimate them using an accurate modeling of the slitless spectrogram.

Kinematic sampling
The impact of resolved kinematic on slitless spectra can be roughly quantified using the kinematic sampling R kin , defined as the line-of-sight velocity resolution in km s −1 px −1 : where D is the spectral sampling (in Å px −1 ), λ o = (1 + z)λ e the cosmologically redshifted wavelength of a line at rest-frame wavelength λ e , and R ≡ λ/∆λ ≈ λ o /D the resolving power. As defined, a smaller R kin corresponds to a higher sensitivity to internal kinematics. The kinematic sampling estimated for various current and future slitless surveys is presented in Table 1. Unfortunately, current HST-based surveys only have a kinematic sampling 800 km s −1 px −1 , which might prove barely sufficient to derive precise velocity parameters for a large fraction of the sample. Future surveys however will all reach R kin ∼ 150−300 km s −1 px −1 , more appropriate for detailed kinematic analyses.
In order to properly constrain galaxy internal kinematics from slitless spectroscopy, we now build a predictive model sufficiently realistic to be quantitatively compared to actual observations in a forward approach. As noticed earlier, this model only depends on a galaxy model (Sec. 2.3) -including flux distribution, intrinsic spectrum and kinematic parameters -and on an instrumental signature (Sec. 2.4) -allowing a complete simulation of the observed spectrogram.

Core assumptions
By design, slitless observations essentially target Hα-emitting galaxies at high redshift (0.5 < z < 1.5, see Sec. 4). Additionally, given the limited spatial resolution and kinematic sampling (Table 1), the observations cannot yet constrain elaborate models. Therefore, we are making the physical assumption (yet reasonable for galaxies considered in this analysis) of an axi-symmetric thin cold disk. As a consequence, the galaxy flux distribution, spectrum and velocity field only depends on the internal radius R.
We further adopt the separability assumption: the rest-frame galaxy spectrum is supposed to be uniform over its whole extent, and only modulated by Doppler shift from internal kinematics. As a result, the observer-frame galaxy datacube C(r, λ) -with two spatial dimensions and a spectral one -can easily be reconstructed from the normalized spatial flux distribution F(r), the cosmologically redshifted intrinsic spectrum S (λ) and the lineof-sight velocity field v(r) of the galaxy:

Spatial flux distribution
As mentioned earlier, slitless spectrography is plagued by selfconfusion, mixing different spatial and spectral contributions A&A proofs: manuscript no. aa Notes.
Under the assumption of an axi-symmetric thin cold disc, the galaxy morphology is characterized by its inclination, position angle and intrinsic radial flux profile. Such a model will be used to estimate F(r) in simulations (Sec. 3).
For real observations (Sec. 4.3), however, given the morphological variety of galaxies -e.g. spiral arms or disc warps -we choose to estimate the galaxy flux distribution directly from the thumbnail broadband images B(r) systematically acquired along slitless spectroscopic observations. We note that, even though our model will essentially be constrained by emission lines (see Sec. 2.3.3), the broadband image, acquired in a band covering the spectrograph band pass, is a direct observation of the integrated flux mostly originating from the continuum. For high redshift galaxies (z 1), the spatial Hα emission line profile is shown to be similar to the continuum but more extended on average (Nelson et al. 2012). To allow for this difference, we relate the internal flux distribution of the galaxy F(r) to the peak-normalized broadband image B(r) through a simple power law: where a flux distribution index η < 1 corresponds to a more diffuse distribution than broadband one's. We discuss further the impact of our choices -in particular the fact that the flux distribution might not be the same in the continuum and in the emission line -in Sec. 5.

Intrinsic spectrum
The 3D-HST and GLASS surveys use grisms G141 and G102 in the infrared domain and cover the 7500 to 17 500 Å wavelength range (see Sec. 4.1). As shown in Sec. 2.2, the kinematic impact is expected to be at most subtle in the slitless spectrograms, and only significant for strong emission lines, namely the complex Hαλ6563 + [N ii]λλ6548, 6584 + [S ii]λλ6718, 6732 (for a redshift 0.3 < z < 1.7), or the doublet [O iii]λλ4959, 5007 (for 0.7 < z < 2.5). We do not consider fainter emission lines such as Hβ.
The intrinsic galaxy spectrum is modeled as a sum of individual Gaussian lines on top of a smooth continuum: where A i and λ i are respectively the amplitude and the restframe wavelength for each line i, σ the supposedly constant line width, and C(λ) an ad hoc continuum. For the line doublets, we assume a constant amplitude ratio of Since the spectrogram adjustment is performed on a very restricted range around the modeled emission lines, and the instrumental transmission is assumed to be known (see Sec. 2.4), the continuum C is simply modeled by a constant C 0 which is the same for all lines.
Overall, only a handful of parameters are needed to describe the intrinsic spectrum, namely 6 for an Hα-emitting galaxy: effective redshift z, line amplitudes A Hα , A [N ii]λ6584 and A [S ii]λ6718 , effective dispersion σ and continuum constant C 0 .

Velocity field
In this section, we present how we construct a model for the galaxy velocity field v(r). Even under simplifying hypotheses, a physical modelling of galaxy rotation curves (hereafter RC) requires a detailed description of the contributions from the disc, bulge and halo components to the galaxy dynamics, to be constrained only with high precision morphologic and spectroscopic observations (e.g. Courteau 1997), out of reach to low-dispersion slitless spectrography.
Under the assumption of an axi-symmetric thin cold disc, one can revert to an analytic expansion to reproduce the overall shape of the intrinsic rotation velocity curve v rot (r) with a restricted number of empirical parameters. In our case, we use a simple hyperbolic tangent profile, very similar to the commonly used arc-tangent profile (Stott et al. 2016;Pelliccia et al. 2017): where v 0 is the plateau value of the RC and r 0 is the transition (or turnover) radius. In practice, we use expression (7), since the CVG w 0 ≡ v 0 /r 0 is the dominant term at small radius, leading to less degeneracy between the parameters. We note that v 0 can only be significantly constrained if the observations go beyond r ∼ r 0 . As briefly explained in Sec. 2, there is a competition between the RC turnover radius r 0 and the galaxy disc scale length r d : if r 0 r d , then both v 0 and w 0 can be reasonably constrained; alternatively, if r 0 is significantly larger than r d , only the solid body rotation parameter w 0 can be sensibly measured.
For a thin cold disk, the observed mean velocity field v(r) along the line of sight is straightforwardly given by: where cz is the systemic velocity, v rot (R) is the rotation curve, i the galaxy inclination and θ the azimuthal angle in the plane of the galaxy: with PA the relative position angle, r = (x, y) the Cartesian coordinates in the sky, and (x 0 , y 0 ) the galactic center coordinates. The galaxy being modeled as a cold rotating thin disc, the kinematic and morphologic position angles are assumed to be the same.
Since PA is highly degenerate with w 0 (see Fig. 1), it is crucial to constrain it independently from photometry. To do so, we estimate the projection angle i and relative position angle PA from a Sersic fit (Sersic 1968) to the broadband image used for the flux distribution. However, since only v rot (R) sin i is adjusted, we stress out that the inclinaison i is not needed in the fit per se, but only for post-analysis velocity deprojection if needed.

Instrumental model
The second key ingredient needed to simulate the slitless spectrogram is the dispersion law D(λ) from the spectrograph, as well as its transmission curve T (λ).
Dispersion law. The dispersion law D(λ) gives the (x, y)-offset on the detector (with respect to a reference position) as a function of wavelength. This is mostly an instrumental quantity, derived from dedicated calibration procedures (e.g. Kuntschner et al. 2009a,b).
Even though the forward approach described in Copin (in prep.) would be an appropriate way to calibrate the dispersion law, we rely in this analysis on the WCS solution computed and delivered for each galaxy by standard data reduction. Given the required precision of our model, we observed minor inconsistencies for some spectrograms, in the form of a px-scale offset in the cross-dispersion direction, as a result of a small registration error between broadband and dispersed images. To account for this effect, we introduce a nuisance parameter ∆y (in px). We note that a similar mis-registration along the dispersion axis is corrected to first-order by the effective redshift z adjusted in the procedure.
Transmission. The transmission T (λ) conveys the chromatic evolution of the instrumental response, and is derived from specific flux calibration (Kuntschner et al. 2011). In our model, we simply include the provided transmission into the galaxy datacube C(r, λ) derived from Sec. 2.3.

Spectrogram reconstruction
Once the galaxy datacube C(r, λ) and the dispersion law D(λ) are known, we compute the resulting spectrogram using (e.g. Freudling et al. 2008): In practice, the spatial convolution and wavelength integration are performed in Fourier space (Copin, in prep.).
Since both quantities C(r, λ) and D(λ) depends on various parameters p, we are able to simulate a spectrogram to be compared to the observations. This "forward" approach allows to constrain physical or instrumental parameters directly in the observation space, in opposition to more traditional "inverse" methods extracting fully or partially free parameters from data using ad hoc procedures.

Maximum likelihood
For an observed spectrogram D i j with estimated variance σ 2 i j , we compute a spectrogram model I i j ( p) and constrain the set of free parameters p by minimizing the following χ 2 : The maximum likelihood procedure further provides the full parameter covariance matrix, from which one can compute the 1σ uncertainties. Since only the emission lines are significantly distorted by resolved kinematics, the adjustment is performed on a restricted area about 40 × 60 px around emission lines of interest, equivalent to 5 . 2 × 700 Å (resp. 1400 Å) for galaxy respectively observed with grism G102 (resp. G141).

Fiducial case
We present realistic simulations designed to mimic WFC3-GLASS observations with grism G102 (see Table 1) in a favorable case where the disc scale length r d is similar to the RC turnover radius r 0 . We first simulate spectrograms with the following properties (see Fig. 3): flux distribution: axisymetric exponential profile with a disc scale length r d = 0 . 6 ∼ 5 px, an inclination i = 60 • and a relative position angle PA = 0 • . intrinsic spectrum: from 7500 to 11 500 Å with 5 emission lines to simulate Hα+[N ii]+[S ii] complex at z = 0.6, with a constant line width σ = 5 Å and a constant continuum; velocity field: v 0 sin i = 250 km s −1 and w 0 sin i = 420 km s −1 arcsec −1 corresponding to a turnover radius r 0 = r d ; instrumental model: linear dispersion law D(λ) aligned along the x-axis, and simplified apodized transmission.
A constant gaussian noise component is finally added, so that the Peak Signal-to-Noise Ratio (PSNR) is around 40 (a typical value for high quality slitless HST spectra).
Article number, page 5 of 13  The fit residuals without (resp. with) kinematics is shown in the middle (resp. bottom) panel of Fig. 4. The minimal χ 2 computed on a rectangular region of 40 × 60 px around the Hα and [S ii] lines decreases from 2688 for 40×60−7 = 2392 degrees of freedom (DoF) to 2380 with only two additional kinematic parameters (w 0 , v 0 ). This ∆χ 2 = −106 decrement for a model with only two extra parameters has a one-tailed p-value of ∼ 10 −23 , corresponding to a 10.0σ detection level. As expected, the residual map without kinematics displays a clear signature of the velocity field as a coherent quadrupolar structure around Hα line position; this structure disappears in the residual map with kinematics.
The correlation matrices are presented in Fig. 5 for the fit without and with kinematics. Note that the kinematics parameters v 0 sin i and w 0 sin i are slightly anti-correlated, but are almost uncorrelated to the other adjusted parameters: no other galactic or instrumental parameter can mimic a kinematic signature.
To test the parameter distribution, we perform the fit for 500 different realizations of the gaussian noise with a PSNR of 40 in the same configuration, and we present the marginalized distributions of measured parameters v 0 sin i, w 0 sin i and redshift z in Fig. 4. All distributions are consistent with the input param- eters in the simulation, and reasonably gaussian. As expected from correlation matrix (Fig. 5), there is a slight anti-correlation between kinematic parameters.

Impact of simulation parameters
The following paragraphs are dedicated to different effects which can impact the measure of both the kinematic parameters and redshift.
Impact of Signal-to-Noise Ratio. We consider the same model as described previously, but perform the fit for 500 realizations of a gaussian noise with a lower PSNR of 20. As expected, the kinematic parameter marginalized distributions are broader around the true values, but there is no hint of biases at lower SNR.
Impact of r d /r 0 . We construct simulations with the same attributes as the fiducial case but with r d /r 0 = 1/2 or 2, with fixed r d = 0 . 6 and v 0 sin i = 250 km s −1 . The marginalized distributions of v 0 sin i, w 0 sin i and redshift z are presented in Fig. 7. As previously illustrated in the pedagogic case (Sec. 2.1), the case r 0 /r d = 1/2 (left panels) is favorable: the plateau of the galaxy RC is reached within the photometric extent of the galaxy, and the maximum velocity v 0 sin i is well measured. In the opposite case r 0 /r d = 2 (right panels), the velocity turn-over radius lays outside the extent of the galaxy disc, and velocity v 0 sin i is only marginally constrained. However, since the inner solid body rotation part always lies within the disc extent, the CVG w 0 sin i is always estimated with a similar accuracy.

Application to HST-based observations
In this section, we present the results obtained by applying the forward model to real galaxy spectrograms from both 3D-HST and GLASS surveys.

3D-HST and GLASS surveys
The 3D The GLASS survey is a HST large program (140 orbits) with the goal of measuring grism spectra over the field of ten massive galaxy clusters at redshift z = 0.31 − 0.69 (for Hα emitters) (Schmidt et al. 2014;Treu et al. 2015). The WFC3 is also used in this survey covering wavelength between 0.75 and 1.75 µm using both grism G102 and G141 to observe the cluster cores. The spectra are acquired at two almost orthogonal roll angles (90 ± 10) • to ease cross-decontamination. It results in a catalog of 1272 redshifts down to M AB ≤ 26 (1060 redshifts with M AB ≤ 24). The WFC3 G102 grism has a spectral coverage from 0.75 to 1.15 µm, a dispersion about 24.5 Å px −1 (twice smaller in resampled spectrograms) and a resolution of 210 in the primary spectral order. The analysis was performed using the 3D-HST reduction pipeline .
We would have liked to test other HST-based surveys such as FIGS (Pirzkal et al. 2017), but the reduced spectrograms are not yet publicly available. The on-going FIGS survey is providing more than 10 000 spectrograms of 2000 different sources (each object is observed at five different roll angles). It has a huge potential for our analysis since multi-roll joint adjustments (as presented in Sec. 4.3.2) would improve constraints on the redshift and kinematic parameters.

Kinematic sample preselection
Keeping in mind that the kinematic signature is expected to often be evasive on the HST spectrograms, we do not foresee to detect it for all targets. We therefore applied some criteria on both 3D-HST 1 and GLASS 2 catalogs to select the most promising targets.
From the 33 559 initial Hα and [O iii] emitters in the G102 and G141 wavelength domains, we first selected highly significant emission lines (F line > 25 × 10 −17 erg s −1 cm −2 and > 15σ detection level on flux for 3D-HST and quality factor Q ≥ 3 for GLASS). We further selected well-resolved and bright galaxies -allowing an accurate measurement of relative PA -by applying cuts on their effective radius (r e > 0 . 6 ∼ 5 px at WFC3 scale) and integrated magnitude (M F140W < 22), as well as moderately inclined galaxies (20 • < i < 80 • ): face-on galaxies have a vanishing apparent velocity field, while edge-on galaxies are generally not well approximated by a thin cold disk.
A final visual inspection of the spectrograms of the 386 preselected candidates was performed to discard severely contaminated spectra, and other data issues. Broadband images were also examined to remove galaxies with highly asymmetric flux distributions, on-going mergers or any other complex structures that could not be handled within our model assumptions.
This selection process picked out 87 galaxies (57 with  ) from GLASS survey. This subsample is only a minimal pre-selection, from which we present the most promising cases in terms of kinematic signature. In this proof-of-concept analysis, we do not try to estimate the overall fraction of targets over which the kinematic analysis is prone to provide accurate measurements; this should be the subject of a forthcoming study once the methodology is applied in a systematic way over large simulated and/or observed samples.

Fitting procedure
The adopted method to probe kinematics from the spectrograms is the following.
1. The position angle is estimated from a Sersic fit to the broadband image, and is kept fixed afterwards. This preliminary measurement was necessary since it is a critical parameter of the model. 2. A kinematic-less fit (v 0 = w 0 ≡ 0) is performed as a reference, providing an estimate of spectral and nuisance parameters. 3. The final fit including kinematic parameters v 0 and w 0 is performed using previous estimates as initial guess. 4. The one-tailed p-value is computed from the best-fit χ 2 without and with kinematics to assess the significance of the fit improvement with the addition of two kinematic parameters. The p-value is converted into a kinematic detection z-score, expressed in σ.
In the following sections, we present some particular kinematic detection by computing the objective χ 2 on a rectangular region of 40 × 60 px around the emission line of interest, Hα complex or [O iii] doublet. Table 2 summarizes the main parameters values and uncertainties of the various fits detailed below.

GLASS results
Single line from single roll angle. We present the result of the fit performed on galaxy #1134 from the GLASS survey, illus-   (7) 55.7 ± 2.9 45.9 ± 4.4 -50.7 ± 1.3 A O iii (7) -43.5 ± 3.0 96.9 ± 1.5 σ (8) 2.96 ± 0.15 5.69 ± 0.27 6.38 ± 0.09 9.46 ± 0.49 η (9) 1.014 ± 0.012 0.981 ± 0.011 1.048 ± 0.005 1.008 ± 0.015 z-score (10) 12.7 13.2 25.5 6.2 trating a very significant kinematic detection. This galaxy has a strong emission complex Hα+[N ii]+[S ii] and a disc scale length r d = 0 . 45. The fit residuals without (resp. with) kinematics is shown in the bottom left (resp. right) panel of Fig. 8. We note that the residual map without kinematics displays a coherent structure around Hα line position, which is significantly reduced in the residual map with kinematics; this indicates that this structure was a signature of the velocity field. The fit gives the following results for the kinematic parameters: v 0 sin i = (205 ± 28) km s −1 , w 0 sin i = (242 ± 28) km s −1 arcsec −1 which corresponds to a RC turnover radius r 0 = 0 . 84 2r d ≈ 6 px. The χ 2 improvement between the two models (without and with kinematics) corresponds to a 12.7 σ detection (p ∼ 10 −33 ) of the kinematic signature. total χ 2 = χ 2 (Hα) + χ 2 (O iii), where each contribution is computed from its own region. Both components share the same intrinsic parameters (redshift, v 0 sin i, etc.), with the exception of the continuum level C 0 and cross-dispersion offset ∆y which can be different. We show an example of this approach on galaxy #451 from the GLASS survey. The fit residuals with kinematics on each emission line region is shown in Fig. 9.
Multiple roll angles. As noticed earlier, each galaxy in the GLASS survey was observed with two almost orthogonal satellite roll angles, corresponding to different dispersion directions and therefore relative PAs. Parameters describing its intrinsic spectrum, flux distribution and velocity field are common to both spectrograms: only the transverse offsets ∆y of the dispersion law may be different. As a matter of fact, we also observed inconsistencies between roll angles in the longitudinal component of the dispersion law (i.e. wavelength solution), leading to spuriously different redshifts when estimated from individual spectrograms. To account for this effect, we introduced another nuisance parameter, the wavelength solution offset ∆x (in Å) such that the observed wavelength of an emission line may differ in both roll angles: where z is the joint (effective) redshift and λ 0 is the restframe wavelength.
We are now able to construct a model describing both spectrograms and fit them simultaneously by minimizing the joint χ 2 = χ 2 (PA 1 ) + χ 2 (PA 2 ), where each χ 2 (PA i ) is the objective function defined in Sec. 2.5.2 for an individual spectrogram with given relative PA i .
We apply this simultaneous fit on the [O iii] emission line of galaxy #399 from the cluster MACS0717 of the GLASS survey, measured with a disc scale length r d = 0 . 44. Each roll-angle spectrogram is first adjusted on its own, leading to a 14.7 σ and 21.6 σ internal kinematic detection respectively; the joint fit on the two spectrograms provides a kinematic z-score of 25.5 σ (see Fig. 10). The fit gives the following results for the common adjusted kinematic parameters: v 0 sin i = (287 ± 15) km s −1 , w 0 sin i = (605 ± 32) km s −1 arcsec −1 which corresponds to a RC turnover radius r 0 = 0 . 47 r d . We find a final redshift z = 1.691 75 ± 0.000 06, accounting for the wavelength offset ∆x = (−12.37 ± 0.32) Å; this corresponds to a shift of about ± ∼ 0.5 px between the wavelength solutions. It is therefore crucial to include it in the model to constrain subtle sub-pixel effects such as kinematics from multi-roll angle spectrograms.

3D-HST results
We present the result for galaxy #19843 from the 3D-HST survey, illustrating a weaker kinematic detection than previously shown. This galaxy has a strong Hα+[N ii]+[S ii] emission complex and a disc scale length r d = 0 . 48. As shown in Fig. 11, the kinematic detection is only at the 6.2σ level (p ∼ 2 × 10 −10 ).

Core model assumptions
Axisymetric thin cold disk. Under the cold rotating thin disc hypothesis, morphological and kinematic position angles are assumed to be the same. However, the misalignment between these two position angles can be quite significant (> 30 • ). Contini et al. (2016) and Harrison et al. (2017) estimated that such misalignment can be attributed to either low/high inclinaison, dispersion dominated systems or complex morphological substructures, such as central bars, spiral arms or clumps in a low-surface brightness disc.
As the CVG w 0 is strongly correlated with the assumed kinematic axis (see Eq. (1)), our method is particularly sensitive to the presence of central asymmetric structures, such as bars and other kinematically distinct components (KDC), incompatible with the cold thin disk approximation. The method presented here therefore remains best suited to rotation-dominated isolated late-type spiral galaxies.
Separability assumption. The second strong hypothesis of our model is more of a technical one helping to construct the threedimensional galaxy datacube; assuming the rest-frame galaxy spectrum to be uniform over its whole extent. This assumption is acceptable for homogeneous objects but less valid for KDCs or star-forming galaxies. Nelson et al. (2012) studied a sample of 57 galaxies from the 3D-HST survey and compared Hα and stellar continuum maps: for high-redshift galaxies (z ∼ 1), Hα emission globally follows the rest-frame R-band light but tends to be more extended and clumpier. Comparing Hα and continuum effective radii, they found that r e (Hα)/r e (R) = 1.3 ± 0.1. However, this difference has not been confirmed with IFU observations for galaxies at z ∼ 2 (Förster Schreiber et al. 2011).
For our study, we have selected only reasonably compact objects, and systematically found η ∼ 1 ± 0.1, indicating that Hα and continuum flux distribution are reasonably similar.
Overall, the two core assumptions are well justified for large and bright isolated late-type galaxies. Furthermore, given the available instrumental setup, with limited spatial and kinematic samplings, slitless observations can only constrain simplified models (disk rotation curves, age or metallicity gradients, etc.).

Systematics from the position angle
As mentioned in Sec. 2, galaxy relative PA and CVG w 0 induce similar line shape distortions in the spectrogram, and are therefore strongly correlated. By differentiating Eq. (1) at constant effective PA, one can estimate the systematic error ∆w 0 associated to an error on the relative PA of the galaxy: For the galaxy #1134 from the GLASS survey, we find a systematic error on the CVG ∆w 0 ≈ 38 km s −1 arcsec −1 . For both the galaxies #451 and #399, ∆w 0 ≈ 123 km s −1 arcsec −1 . Finally, for the galaxy #19843 of the 3D-HST survey, we compute ∆w 0 ≈ 108 km s −1 arcsec −1 .
Overall, the systematic uncertainties on the CVG w 0 are larger than the statistical ones (except for galaxy #19843) but remains in a similar range since the PA is very well constrained from the broadband images (see Table 2) and the kinematic sampling is the governing term in Eq. (14) where s = 0.13 arcsec px −1 .

Implications on redshift precision
Self-confusion effect in slitless spectroscopy is generally assumed to induce redshift measurements less accurate than in more traditional spectroscopic observations. However, the forward model presented here significantly attenuate this effect, as the effective spectral resolution does not depend on the object extent anymore, even in the absence of kinematic parameters. In Table 2, we compare the redshift uncertainties from the forward model described in Sec. 4.3 to the total uncertainties (including systematics) estimated by 3D-HST and GLASS analysis. We emphasize that only the statistical error on the redshift measurement is quoted in our case, as we do not have access to the data reduction details for HST-based surveys. Nonetheless, we point out that the forward modeling delivers a redshift accuracy down to ∼ 10 −4 or less. Using a similar forward approach for slitless spectroscopic reduction would probably provide an equivalent precision gain on instrumental calibrations (notably dispersion solution), which would in turn decrease the final systematic errors to a comparable level. Ultimately, our study suggests that a consistently forward analysis of slitless observations -including for data reduction procedures -could allow a significant gain in redshift precision.  Fig. 10. Same as Fig. 8 for a joint fit on [O iii] emission line of galaxy #399 from GLASS survey observed with a roll angle of 52 • (left panels) and 152 • (right panels). We note that for the first spectrogram, the large residual structure which is visible at the bottom results from an incomplete source decontamination: it does not significantly affect the model adjustment result, but only increases the final χ 2 value. The decrease in the joint χ 2 corresponds to a 25.5 σ-level detection.

Kinematic measurements
Notwithstanding its direct degeneracy with the assumed kinematic angle, the CVG w 0 is reasonably constrained by the core 2D shape of the emission lines in the slitless spectrogram. On the other hand, the plateau velocity v 0 remains difficult to measure, since spectroscopic observations are rarely deep enough in available redshift surveys to probe regions well beyond effective radius r d ≈ r 0 where the velocity flattening would be manifest. Furthermore, current instrumental setups, with a limited kinematic sampling (R kin > 650 km s −1 px −1 ), do not favor precise measurements of internal kinematic signatures, inducing sub-px spectral distortions.
We note that in this proof-of-concept analysis, we did not take into account the "beam smearing", i.e. the degradation of kinematic resolution due to limited spatial resolution. Formally, one would need to construct an "infinitely" spatially resolved model cube before applying the spatial PSF, while, in our case, we built the model cube from an already PSF-convolved flux model (Eq. (3)). This would mostly change the estimated value of w 0 for very steep unresolved CVGs.
Last, the natural drawback of the forward approach is that the model is assumed to be a fair representation of the observed data.
If this is not the case (specifically, non cold thin disc-like objects and/or non-radial structures), errors are dominated by systematics and the resulting kinematic parameters have no adequate physical interpretations.

Conclusions and perspectives
In this article, we explored the possibility to probe single object with slitless spectroscopy by measuring not only integrated spectral features but also spatially resolved quantities such as internal kinematics. To achieve that, we presented a forward-model of slitless spectrograms from a galaxy model -including flux distribution, intrinsic spectrum and kinematic parameters -and an instrumental signature. This method relies on two major assumptions: the axi-symmetric thin cold disk approximation and the separability hypothesis. We applied this method on HST-surveys galaxies to measure internal kinematics parameters: the plateau velocity v 0 and the central velocity gradient w 0 .
The kinematic signature on slitless spectrograms is very specific, as a quadrupole structure in the kinematic-less residuals around emission lines. It is therefore difficult to confuse with other effects such as line flux distribution or radial structures. Even if the kinematic signature is subtle, it extends on the full galaxy scale, which makes the statistical detection significant even with limited kinematic sampling.
The mains results of our method can be summarized as follows: -We present the first detection of resolved internal kinematics in galaxies from slitless spectroscopic observations. Using a simple asymptotic velocity curve, we note that the CVG w 0 is more precisely determined than the plateau velocity v 0 , but is directly sensitive to systematic errors on the assumed relative position angle of the galaxy, which must be estimated from external photometry. -The forward model naturally corrects for the self-confusion effect: the spectral resolution does not explicitly depend on the object shape/extent anymore. It implies a more precise redshift measurement (down to 10 −4 ) with a simple model without kinematics and even better with a model including kinematics. Improved redshifts from slitless surveys will have a direct impact on cosmological probes as Baryon Acoustic Oscillations and Redshift Space Distortions measurements, where the expected precision is 10 −3 to 10 −4 . -We observe sub-px inconsistencies in the dispersion law of the spectrograms from 3D-HST and GLASS surveys in the cross-dispersion direction, and use a nuisance parameter to account for it. However, calibration errors along the dispersion axis induce a systematic error in the effective redshift determination, but do not affect the kinematic parameters. It should be kept in mind that this analysis is a proof of concept applied at the data analysis level only: optimal results would be obtained with a full forward model of slitless spectrograms, including in data calibration stages. -We demonstrate the great flexibility of our forward model, where one can simultaneously fit different emission lines, spectrograms observed with different dispersion direction, and even spectrograms from different instruments, allowing improved constraints on redshift and kinematic parameters.
The spectroscopic resolving power R is a key parameter determining the amplitude of the Doppler distortion in the spectrogram and therefore is a deciding quantity to constrains precisely the CVG and the plateau velocity. Furthermore, the spectroscopic observations must be deep enough to probe the velocity curve well beyond the photometric disk radius and reach the asymptotic regime.
We stress out that, building on a finer understanding of the spectrogram formation, this forward model of slitless spectra has the potential to obtain more precise redshift measurements than standard approaches. This will be massively tested on futures slitless surveys such as Euclid which will acquire spectrograms for 30M galaxies in the Wide Field and will also produces deeper spectroscopic observations in the Deep Fields (Laureijs et al. 2011). Moreover, these surveys will be characterized by a finer kinematic sampling R kin ∼ 200 km s −1 px −1 , more suited for kinematic measurements and galaxy scaling relation studies such as the Tully-Fisher relation (e.g. Aquino-Ortíz et al. 2018;Lelli et al. 2019) or morpho-kinematic classification (e.g. Cortese et al. 2014;Kalinova et al. 2017;Schulze et al. 2018).