A Kinematic Analysis of Ionised Extraplanar Gas in the Spiral Galaxies NGC 3982 and NGC 4152

We present a kinematic study of ionised extraplanar gas in two low-inclination late-type galaxies (NGC 3982 and NGC 4152) using integral field spectroscopy data from the DiskMass H$\alpha$ sample. We first isolate the extraplanar gas emission by masking the H$\alpha$ flux from the regularly rotating disc. The extraplanar gas emission is then modelled in the three-dimensional position-velocity domain using a parametric model described by three structural and four kinematic parameters. Best-fit values for the model are determined via a Bayesian MCMC approach. The reliability and accuracy of our modelling method are carefully determined via tests using mock data. We detect ionised extraplanar gas in both galaxies, with scale heights $0.83^{+0.27}_{-0.40}\,\mathrm{kpc}$ (NGC 3982) and $1.87^{+0.43}_{-0.56}\,\mathrm{kpc}$ (NGC 4152) and flux fraction between the extraplanar gas and the regularly rotating gas within the disc of 27% and 15% respectively, consistent with previous determinations in other systems. We find lagging rotation of the ionized extraplanar gas in both galaxies, with vertical rotational gradients $-22.24^{+6.60}_{-13.13} \,\mathrm{km\,s^{-1}\,kpc^{-1}}$ and $-11.18^{+3.49}_{-4.06}\,\mathrm{km\,s^{-1}\,kpc^{-1}}$, respectively, and weak evidence for vertical and radial inflow in both galaxies. The above results are similar to the kinematics of the neutral extraplanar gas found in several galaxies, though this is the first time that 3D kinematic modelling of ionised extraplanar gas has been carried out. Our results are broadly consistent with a galactic fountain origin combined with gas accretion. However, a dynamical model is required to better understand the formation of ionised extraplanar gas.


INTRODUCTION
Spiral galaxies are surrounded by multiphase gas layers extending up to several kpc from their disc planes (see Reynolds et al. 1973;Lockman 2002;Marasco & Fraternali 2011 for observations of the Milky Way's gas layers and Dettmar 1990;Rand et al. 1990;Swaters et al. 1997;Oosterloo et al. 2007 for observations of other galaxies), which are detected both in emission (H and H ) and in low-and high-ions absorption against bright background sources (Lehner & Howk 2011;Zheng et al. 2017). These gas layers, known as extraplanar gas (EPG), are located at the disc-halo interface where gas accretion and stellar feedback take place. Studying their properties gives us fundamental insight into the gas cycle between galactic discs and the circumgalactic medium (CGM), which in turn is key to understanding galaxy formation and evolution mechanisms.
Neutral hydrogen gas (H ) is a primary tracer of EPG, and extensive observations of neutral EPG have been conducted. Studies find that neutral EPG is ubiquitous in late-type galaxies, with H mass fractions relative to the underlying neutral gas mass in the disc varying from 10% to 30% (e.g. Oosterloo et al. 2007;Marasco et al. 2019). Observations of nearby spiral galaxies with different inclination angles (nearly face-on: Boomsma et al. 2008;intermediate inclination: Fraternali et al. 2001;Hess et al. 2009; ★ E-mail: li@astro.rug.nl 2019; and edge-on: Oosterloo et al. 2007;Rand & Benjamin 2008;Zschaechner et al. 2011) provide a 3D view of the structure and kinematics of neutral EPG in disc galaxies. Neutral EPG has a scale height around 1-3 kpc, and its velocity is characterised by differential rotation similar to the disc, but lagging behind the disc gas with a vertical velocity gradient of typically −10 to −20 km s −1 kpc −1 (Fraternali et al. 2001(Fraternali et al. , 2005Oosterloo et al. 2007;Zschaechner et al. 2011;Marasco et al. 2019). Non-circular motions like radial inflow and vertical flows are also detected in the EPG of some galaxies (e.g. Fraternali et al. 2002a; Barbieri et al. 2005;Marasco et al. 2019).
The presence of EPG in the Milky Way has been known for several decades. It was first detected in the form of extended H complexes with velocities that deviate strongly from the Galactic rotation (Muller et al. 1963;Oort 1966). Clouds with deviation velocities larger than 90 km s −1 are called high-velocity clouds (HVCs) (Wakker & van Woerden 1997;Westmeier et al. 2007). HVCs usually have sub-solar metallicities and are located at several kpc from the Galactic plane (Wakker 2001). Although first detected in HVCs, most extraplanar H clouds in the Milky Way are at intermediate velocities (referred to as IVCs), which have near-solar metallicities and are within 1-2 kpc of the Sun (Wakker 2001). A kinematic modelling study of the global EPG emission in the Milky Way shows that the EPG in the Milky Way has a velocity gradient around −15 km s −1 kpc −1 and inflow in the radial and vertical directions (Marasco & Fraternali 2011). This model shows that the "classical" IVCs are simply the local (Solar Neighbourhood) manifestation of a widespread (across the whole Galaxy) EPG component.
The similarity between neutral and ionised EPG kinematics suggests that these two phases might have the same origin. Rossa & Dettmar (2003a) show that galaxies with larger star-formation rates tend to possess more prominent ionised EPG (also seen in Jones et al. 2017). Previous studies find that the distribution of the EPG is often correlated with the locations of the H regions in the disc (Miller & Veilleux 2003a;Boomsma et al. 2008). The ionisation source for ionised EPG is most likely to be photons leaking from star-forming regions (Levy et al. 2019). Thus, a possible scenario is that the ionised EPG and neutral EPG are counterparts: part of neutral EPG is photoionised by star formation activity and turns into ionised EPG.
The origin of EPG is a long-standing issue, and two different mechanisms have been proposed. One is an internal mechanism: a galactic fountain powered by star formation activities like supernova explosions or stellar winds ejects gas above the galactic disc to form the EPG (e.g. Shapiro & Field 1976;Bregman 1980). The other is an external mechanism: accretion of gas from the CGM (Binney 2005;Kaufmann et al. 2006). A more plausible explanation of the origin of EPG involves both mechanisms (discussed below).
A galactic fountain is the best-studied mechanism for EPG formation and maintenance. In general, H gas in the disc is partly ionised by photons from star-forming regions and pushed above the disc by supernova explosions or stellar winds. The gas then travels through the halo within the gravitational potential of the galaxy, cools down and eventually falls back to the disc. A strong prediction of such a model is that the EPG kinematics should follow the galaxy rotation, which is confirmed by the observations. Simple ballistic models for the galactic fountains reproduce many of the observed properties of the EPGs (e.g. Fraternali & Binney 2006). Nevertheless, an external origin of EPG cannot be ruled out, as a pure galactic fountain model tends to predict a radial outflow and underestimate the velocity gradient of the EPG (Heald et al. 2006;Fraternali & Binney 2006). To solve these issues, an external mechanism that can alter the angular momentum of the EPG is required. A possible scenario is fountain-driven accretion, where interaction with the fountain's neutral clouds decreases the cooling time of the hot coronal gas to a time scale shorter than a cloud's orbit time, and thus part of the hot gas is condensed and accreted onto to the disc (Fraternali 2017). Models based on this scenario explain the EPG kinematics observed in external galaxies (Fraternali & Binney 2008) and the phase-space distribution of the cold (neutral) and warm (ionised) EPG in the Milky Way (Marasco et al. 2012Fraternali et al. 2013).
A comprehensive 3D view of the multi-phase EPG kinematics is needed to disentangle the origin of the EPG. Marasco et al. (2019) studied neutral EPG in a sample of 15 galaxies at intermediate inclinations using the HALOGAS survey . They modelled the EPG and characterised its structural and kinematic properties, focusing in particular on the spatial distribution, vertical and radial velocities, rotational gradient, and velocity dispersion, providing the first systematic and homogeneous study of neutral EPG in nearby disc galaxies.
As for ionised EPG, however, nearly all samples are composed of edge-on galaxies, which limits or even prevents analysis of noncircular (especially vertical) motions of the EPG. We therefore have chosen to study the ionised EPG of a pair of intermediate-inclination galaxies in this paper. In Sec. 2 we introduce the DiskMass sample from which we select our galaxies. The strategy used for modelling ionised EPG, adapted from Marasco et al. (2019), is described in Sec. 3 and tested with mock galaxies in Sec. 4. We present our analysis and results for ionised EPG in a pair of DiskMass sample galaxies in Sec. 5 and discuss the results in Sec. 6. We summarise our results in Sec. 7.

The H sample of the DiskMass Survey
In this paper, we use the DiskMass survey (DMS) H sample (Bershady et al. 2010;Swaters et al., in prep.) to study the existence and kinematic properties of ionised extraplanar gas (EPG). This sample contains integral field spectroscopy of 138 late-type galaxies in the DMS parent sample (Bershady et al. 2010) observed with the SparsePak integral field unit on the WIYN telescope . SparsePak is a sparsely distributed fibre optic bundle containing 82 fibres, including seven sky fibres. Each fibre has a diameter of 4. 7 and is separated by 5. 6, and the fibres are sparsely distributed on the focal plane. All 138 galaxies have low or intermediate inclination angles and diameters comparable to the field-of-view of SparsePak (1 < 25 < 3. 5).
Each galaxy in the DMS H sample has been observed with three pointings of different offsets. One pointing is centred at the galaxy centre, a second pointing is 5. 6 to the south, and a third pointing is 4. 9 to the west and 2. 8 to the north. Since the distribution of fibres follows a hexagonal pattern , the SparsePak field-of-view can be almost fully sampled with this 3-pointing strategy.
The SparsePak signal is fed into the Bench Spectrograph, a spectrograph optimised for use with fibre optic bundles. For the DMS H sample, the 316 lines mm −1 echelle grating in order 8 was used, which covered a wavelength range of 6475-6880 Å, with a dispersion of 0.195 pix −1 and an instrumental FWHM of 30 km s −1 .
The spectral data we use in this paper have been reduced in the following manner (Swaters et al, in prep). The raw data have been over-scan corrected and trimmed, and a bias level has been subtracted. Cosmic rays have been removed using the method introduced in Andersen et al. (2006). The IRAF 1 task dohydra has been used to flat-field, wavelength-calibrate, and extract the spectra. Sky and (stellar) continuum have then been removed. The data product for each pointing is a FITS file in row-stacked spectra (RSS) format that contains 82 emission-line spectra with a wavelength range of 6475-6880 Å and a wavelength sampling of 0.2 pix −1 , which corresponds to 9.14 km s −1 around the H line.

The subsample in this paper
In this work we focus on two galaxies from the DMS sample, NGC 3982 and NGC 4152. We select these two galaxies because they have sufficient S/N to reveal the relatively weak EPG emission, and their structures are well-resolved under the SparsePak spatial resolution. Besides, their position-velocity (pv) diagrams along their major axes have emission features called 'beard' patterns (Schaap et al. 2000;Fraternali et al. 2001) that deviate from the bulk rotation and extended towards the systemic velocity. These two galaxies are the best cases in the DMS to study the EPG. This does not mean that other galaxies in the DMS sample do not have EPG, but rather they do not have high-enough S/N or spatial resolution to allow us comment on the incidence of EPG.
NGC3982 is a late-type Seyfert 2 galaxy (Cardamone et al. 2007) with a relatively low inclination angle (26. • 2) and a distance of 21.6 Mpc (Tully et al. 2013). It is a high surface brightness member of the Ursa Major cluster, with on-going star-forming activity (Martinsson et al. 2013), making it an ideal target for the kinematic study of EPG. NGC 4152 is a late-type galaxy with on-going star formation. It lies at a distance of 30.6 Mpc (Tully et al. 2009), with an inclination angle of 38. • 1. Figure 1 shows the -band images of NGC 3982 and NGC 4152. Their physical properties are listed in Table 1.
To get a first glimpse of the EPG in the data, we start by stacking spectra of NGC 3982 before implementing the 3D-modelling analysis. We fit a Gaussian function to every H velocity profile in NGC 3982 that is within ±30 • of the major axis and then use the fitted intensities to normalise each profile. All normalised profiles are shifted using the fitted central velocity of the Gaussian (spectra on the receding side are flipped). In this way all velocity profiles are centred at zero velocity while positive (negative) values represent rotation faster (slower) than the disc. We then stack those profiles together and get a stacked H spectrum for NGC 3982. The wavelength range of our spectra data also covers the [N ] doublet (6583 Å and 6548 Å), although the S/N of the [N ] doublet does not allow us to detect EPG in most single velocity profiles. Using a similar method, we also derive a stacked [N ] 6583 spectrum. The stacked spectra overlaid by Gaussian fits (based on the upper 40% of the profiles) are shown in the top panel of Fig. 2, also shown are the residuals of the Gaussian fit. It is clear that both H and [N ] 6583 profiles show wings that cannot be fitted by a single Gaussian, which are very likely due to EPG. Note that NGC 3982 is a Seyfert galaxy, the stacked spectra might be contaminated by active galactic nucleus (AGN) feedback. We have therefore also tried masking the central 1 kpc region and verified that the pattern of the wings does not change with respect to that shown in Fig 2. As mentioned, the kinematics of the EPG are characterised both by a lag in rotation with respect to the disc gas and non-circular (vertical and radial) motions. The fact that the wing at the slower-rotation side is stronger is what we would expect from EPG, since the lag only contributes to the slower-rotation side while vertical motion results in wings on both sides. We also look into the [N ] 6583/H ratio of the stacked spectra, which is shown in the bottom panel of Fig. 2 and normalised at zero velocity. The ratio increases for gas at velocities that depart from the disc rotation (peak of the Gaussian) as expected, since EPG usually shows a higher [N ] 6583/H ratio compared to gas in the disc (Hoopes et al. 1999), similar to what Fraternali et al. (2004) found for ionised EPG in NGC 2403.
The stacked spectra along the major axis of NGC 3982 do not give us information about radial motions, which are actually often detected in EPG, as mentioned in Sec. 1. Radial motions are visible in spectra along the minor axis, where inflow (outflow) will cause a blue-shifted (red-shifted) wing at the far side of the galaxy and a red-shifted (blue-shifted) wing at the near side of the galaxy. We therefore also derive stacked H spectra within ±30 • of the minor axis using similar method as described above, but stacking the far side and the near side separately, as shown in Fig. 3. The blue-shifted wing dominants in the far-side stacked spectrum while the red-shifted wing is relatively more prominent in the near-side stacked spectrum, suggesting a radial inflow for EPG in NGC 3982, consistent with studies of other galaxies (see Sec. 1).
The above analysis on the stacked spectra strongly suggests that ionised EPG is present in NGC 3982. In the following we further study the kinematics of ionised EPG by 3D modelling.

Building H datacubes
The 3D modelling code we use to study EPG is implemented on datacubes (Marasco et al. 2019), but the DMS H data is in RSS format, as described above. We take the following steps to reshape the spectra of NGC 3982 and NGC 4152 into datacubes.
First we build a model FITS file for each galaxy with the Groningen imaging processing system (GIPSY, VAN Hulst et al. 1992) task create. This FITS file is a 300 × 300 pixel 2D image with a pixel size of 0. 5. As mentioned above in Sec. 2.1, the distribution of SparsePak fibres follows a hexagonal pattern. We need an over-sampled FITS file to allocate fibres with high accuracy, and therefore the grid size is much smaller than the 4. 7 fibre diameter of SparsePak. The FITS file is centred at the galaxy centre.
Next we convert the wavelengths of the spectra into optical velocities with respect to the H line and use the IRAF task rvcorrect to derive the heliocentric velocity (all velocities in this paper are heliocentric velocities or velocities relative to the galaxy centre). We have spatial offsets for each of the 75 fibres (ignoring seven sky fibres) in each pointing, as well as the spectrum at each pointing. For every point along this velocity grid, we then use the GIPSY task gds2text to fill the model FITS image with intensities at that particular velocity of all 75 × 3 spectra from the three pointings. Finally, we have an image for each velocity along the velocity axis, which is equivalent to a slice of a datacube along the velocity axis. Then we use GIPSY task create again to build a 3D model FITS file with a size of 300 × 300 × 81 pixels. Its spatial grid size (0. 5 pix −1 ) is that of the above-mentioned 2D model FITS file while the third axis is velocity, centred at the systematic velocity of the galaxy with a velocity interval of 9.14 km s −1 (equivalent at H to the underlying spectrum's wavelength interval of 0.2 pix −1 ) and a velocity resolution of 30 km s −1 , corresponding to spectral resolution of 10049. With 81 channels we cover a velocity range of ±365.44 km s −1 around the systemic velocity of the galaxy, sufficient to capture the kinemat-   ics of the disc gas and EPG. We then pile the image slices from the above step into this 3D model FITS file.
Finally, we have a datacube for each galaxy in the subsample that contains spectra from all three pointings. Note that because our field-of-view was not fully sampled, there are blank regions in the datacube. We are aware that the above method is different from the usual way of building IFU datacubes where interpolation is done between adjacent pixels and blank regions are filled. We choose this method under the consideration that the study of EPG kinematics prefers a small beam-smearing effect. However, the beam-smearing effect with our relatively large fibre size (4. 7, corresponding to ∼ 0.6 kpc in our sample) is already non-negligible and would increase further if we interpolated between different fibres. Thus we build the datacube as appropriate to the original observation, at the cost of blank regions in the resulting datacube where no data were available.

MODELLING THE EXTRAPLANAR GAS
We use the EPG modelling code written by Marasco et al. (2019) to study the EPG in NGC 3982 and NGC 4152, which was originally written with the aim of studying the neutral H EPG in the HALO-GAS survey (Marasco et al. 2019) but in principle is also capable of characterising structure and kinematics of generic line-emitting medium outside the disc.
The spatial resolution in the HALOGAS sample is 30 , corresponding to 2.6 kpc for NGC 4414, which is the furthest galaxy in the HALOGAS sample for which the code manages to model the EPG (Marasco et al. 2019). As a comparison, the spatial resolution in our datacube is 4. 7, corresponding to 0.41 kpc for NGC 3982 and 0.60 kpc for NGC 4152. The typical velocity resolution in the HALO-GAS sample is FWHM ∼ 8 km s −1 Marasco et al. 2019) while in our datacube the FWHM is around 30 km s −1 .
Although not comparable with the velocity resolution in the HALO-GAS sample, we argue that a FWHM of 30 km s −1 is sufficient for kinematic study of ionised gas, and mock galaxy tests in Sec. 4 verify  that the code can model the EPG successfully under such velocity resolution. This code (with some customisation) is therefore suitable for our datacubes.
In general, the code first separates EPG from disc emission and then fits the EPG emission with an axisymmetric, smoothly distributed EPG model, which is characterised by three structural and four kinematic parameters. An MCMC algorithm is used to find the best-fit parameters for the data. In the following parts, we briefly introduce the algorithm of the code and the modifications we made for H data; more details about the code can be found in Marasco et al. (2019).

Isolating the EPG signal
Separation of EPG emission from the underlying gas disc and noise in the datacube is done through two masks: an external mask, which filters out noise voxels without emission from the galaxy; and an internal mask, which filters out the disc emission. The external mask is created by spatially smoothing the datacube by a factor of 5 (using a 2D Gaussian kernel), calculating a smoothed rms noise level, then sigma-clipping at S/N = 4, which removes the majority of noise voxels and leaves only voxels with emission.
The internal mask, which separates disc emission from the EPG emission, is created using the method introduced by Fraternali et al. (2002a), under the assumption that the EPG velocity deviates from the disc motion because of EPG's lagging rotation, vertical and radial motion, as well as larger velocity dispersion (see Sec. 2.2 in Marasco et al. 2019 for details). For the velocity profile at a certain spaxel, the disc component is assumed to be described by a Gaussian profile. The lagging rotation of the EPG adds a wing to the profile at the slower rotation side, and the radial and vertical motion of the EPG contribute wings at both sides (a sketch of this can be found in Fig. 1 of Marasco et al. 2019). Although the disc and EPG profiles are blended together, the ratio of the EPG's mass to that of the disc is relatively small, and thus it is reasonable to neglect the contribution of the EPG around the peak of the emission line profile and use this region to fit disc emission. We perform a Gaussian fit for the disc component using only the upper 40% (in intensity) of the line profile. We implement this method on all pixels that have passed through the external masking. By considering all the fitted Gaussian profiles together, we build a non-parametric model cube for the H emission from the disc. Finally, the internal mask is generated from this synthetic-disc datacube: all voxels with disc emission higher than times the rms-noise are masked out. The mask threshold should be low enough to reject the majority of disc emission and high enough to leave as much EPG emission in the data as possible. We choose = 2, which was decided empirically and has been tested to be eligible for separating disc and EPG (see Sec. 2.2 of Marasco et al. 2019). Moreover, in Sec. 5.1 we show that for our data the results of MCMC fit are not significantly sensitive to the mask threshold. Figure 4 shows two example velocity profiles illustrating how the internal mask works.

Extraplanar gas model
After applying external and internal masks to the datacube, we fit the remaining emission with a synthetic datacube built from a parametric model of the EPG. The density distribution of this EPG model follows the profiles introduced in Oosterloo et al. (2007) cylindrical geometry. The surface density profile is given by where Σ 0 is the central surface density, is an exponent regulating the density decline towards the centre, and is the scale length of the EPG. The density distribution along the -direction is given by where ℎ is the scale height of the EPG. , , and ℎ together control the spatial distribution of the EPG model. Four other parameters characterise the EPG kinematics: the rotation velocity gradient along the -direction that describes the lagging rotation of the EPG, / ; the radial velocity, ; the vertical velocity, ; and the velocity dispersion, . These seven parameters ( , , ℎ, / , , , ) together define an EPG model.
The realisation of a parametric model into a datacube is described in detail in Section 2.4 of Marasco et al. (2019). However, we note that while for the H -emitting gas the line intensity is simply proportional to the gas column density at given velocity 2 , this is not the case for the H emission, thus modifications to the previous implementation have to be made. These are described in details in Sec. 3.4. When building the model datacube, the distance, position angle (PA), inclination angle (INCL), and rotation curve of the galaxy are also necessary ingredients. We use distances from the literature ( Table 1). The rotation curve is derived from the datacube before internal masking, using the 3D-tilted-ring modelling code 3D Barolo (Di Teodoro & Fraternali 2015). When running 3D Barolo, we fix the kinematic centre at the optical centre and use PA and INCL derived from GIPSY task rotcur as input values, which are then slightly adjusted (by less than 1 • ) by 3D Barolo and later implemented in the MCMC fitting described below. The inferred rotation curves for 2 Assuming an optically-thin regime, appropriate for low-density gas in the halo

MCMC fitting
We now estimate the values of the seven EPG model parameters using a Bayesian fitting procedure that uses MCMC for sampling the posterior. According to Bayes' theorem, for a given data D, the posterior probability P for parameter vector x is given by where P (D|x) is the likelihood function and P (x) is the prior. We do not make any prior assumptions about the EPG parameters, and thus the prior for each parameter is uniform within a reasonable and wide range (see Table 3 of Marasco et al. 2019). The likelihood function is given by where M represents the model datacube built from parameter x, R is the sum of the absolute residuals between the model and the data, and is the uncertainty of the data. Since the EPG model is axisymmetric, we should take the asymmetry of the data into consideration while calculating the uncertainty. Thus we define an "effective uncertainty" , given by = × , where is an estimator for the deviation of the data from pure axisymmetry, and is the number of voxels per resolution element (see Eq. 6 of Marasco et al. 2019 and their Sec 2.5).
We sample the posterior using , a Python implementation of the Markov Chain Monte Carlo (MCMC) method (Foreman-Mackey et al. 2013). During the fitting, we use 100 walkers and a chain length that varies between 1000 to 2500 steps. The walkers are initialised at a minimum, which is found by a downhill-simplex minimisation routine (Nelder & Mead 1965).

Customisation for H data
As mentioned above, the EPG modelling code was written to analyse neutral EPG in the HALOGAS sample (Marasco et al. 2019). The application of the code to the DMS H data requires customisation, which is related to the difference in beam shapes, the conversion from gas-density to line-intensity, and the treatment of dust extinction.
The synthetic EPG datacube that we build from a parametric model to fit the EPG emission should have the same beam shape as the data. This is achieved by spatially smoothing each channel of the synthetic datacube with the beam of the data. As described in Sec. 2.3, our datacube is built directly from the SparsePak fibre data without interpolation, and thus the data beam is equivalent to the SparsePak fibre's top-hat beam shape, i.e., uniform within each fibre's radius (2. 35) and zero outside fibres, which is contrary to the nearly Gaussian beam shape of H data (as it is in HALOGAS survey). We have therefore changed the smoothing treatment for our model datacubes accordingly.
While H emission intensity is linearly proportional to the density of the emitting gas, H emission depends on the square of the electron density of the extraplanar plasma (Spitzer 1978;Collins et al. 2002). We account for this effect by simply squaring the model datacube derived by assuming linear proportionality between emission intensity and column density per velocity bin. As in Marasco et al. (2019), the cube normalisation is then determined a posteriori, by comparison with the total observed flux. We stress that this squared proportionality affects the shape of the line profiles, and as a consequence, producing a factor √ 2 difference between the density-weighted and flux-weighted velocity dispersions.
The model datacube is built in an optically-thin regime, which is a reasonable assumption for H gas. The ionised EPG gas, however, suffers from dust extinction. Due to the lack of dust extinction data for our sample, we do not introduce extinction factor into our model when fitting the EPG data. Nevertheless, mock galaxy tests show that dust extinction does not influence our analysis on ionised EPG significantly (see Sec. 4 and discussion in Sec. 6.1).
Except for the above mentioned differences, our procedure is the same as that in Marasco et al. (2019).

TESTING EXTRAPLANAR GAS PARAMETER ESTIMATION WITH MOCK GALAXIES
The robustness of the code can be checked via mock data tests, where we build mock galaxies with certain parameters and see if the code can recover these parameters. Mock data tests have been conducted extensively in Marasco et al. (2019) and have verified the reliability of the code. However, because of the above-mentioned differences between H data and H data and the corresponding changes we make to the code, we have also conducted mock data tests.

Building a mock galaxy
The algorithm to build a mock galaxy is described in Appendix B of Marasco et al. (2019). We begin by choosing a distance, PA, INCL, and rotation curve for the mock galaxy. We then construct a regularly rotating thin disc model and choose an EPG model. We use NGC 3982 from the DMS H sample as a template for building mock galaxies. All mock galaxies (except MG7, see below) have the same distance, PA, INCL, and rotation curve as NGC 3982, using parameters inferred using the method described in Sec. 3.2. The disc model assumes a Gaussian profile for both vertical and surface density distribution. We set the disc scale length d to be 1.6 kpc (which is comparable to NGC 3982) and the disc scale height ℎ d to be 0.2 kpc. For the disc kinematics, we assume a velocity dispersion of d = 27 km s −1 ; we note that this is a density-weighted dispersion due to the way we build our datacube (see Sec 3.2), which corresponds to a flux-weighted velocity dispersion of 19.1 km s −1 , a typical value for disc ionised gas (e.g. Martinsson et al. 2013). In Sec. 3.2 above we have discussed how to construct an EPG model.
We then build two mock datacubes from these ingredients: a mock disc datacube and a mock EPG datacube. To produce a mock galaxy datacube, we must first decide on the ratio of the EPG flux to the disc flux. To explore the performance of the code on different flux fractions of EPG and match these with the observations -the EPG flux fractions for NGC 3982 and NGC 4152 are 0.27 and 0.15 (Section 5) -we choose flux fractions of 0.15, 0.20, and 0.30, covering the upper and lower fractions and an intermediate value. Once we have chosen the flux fraction, we scale the EPG and disc mock datacubes by the appropriate flux fractions and then add them to produce the mock galaxy datacube.
Finally, we inject noise into the mock data, assuming a constant rms noise level throughout the datacube, and quantify the corresponding "data quality" using two parameters: [S/N] MED and N ind . [S/N] MED is defined as the median S/N of the voxels in the EPG cube with intensities above twice the assumed rms noise while N ind is defined as the total number of voxels in EPG residual with intensity above twice the rms-noise divided by the number of voxels per resolution element (see Eq. 6 of Marasco et al. 2019). As [S/N] MED and N ind will also influence code performance, we explore models with different values by changing the rms noise.

Mock galaxy parameters
We first build an extreme case, MG1, where the EPG signal is strong (high S/N and flux fraction). MG2 is created by doubling the rms noise of MG1 and decreasing velocity gradient / and velocity dispersion . MG3 is the same as MG2 except that we double the rms noise. Then we move back to the rms noise of MG2, lower the flux ratio of EPG to 0.15, flip the values for and , and build MG4.
The above tests explore different S/N or EPG fractions, but it is noteworthy that in these mock galaxies the EPG has kinematics which significantly deviates from the disc (large / and ℎ). To further explore the parameter space, we test models with moderate S/N and EPG fraction but have EPG kinematics which deviates less from the disc. We therefore build MG5 with smaller / and ℎ. From MG5, we further reduce ℎ, slightly increase / and build MG6. MG7 is similar to MG5, except that we set a higher INCL in order to explore the influence of INCL.
H emission can suffer from dust attenuation. While dust can be present in the inner halo (which itself supports an internal origin for the EPG, Howk & Savage 2000), we expect that the most severe attenuation is produced by dust in the midplane, which would cause a drop in the H intensity produced by portion of EPG located "behind" the disc with respect to the observer. To explore this effect, we use a extinction fraction to quantify the percentage of flux we observe from EPG located in the background of the midplane. Based on the above simplified scenario, we build four more mock galaxies MG8-MG11, where = 0.7, 0.5, 0.25, 0.0, respectively; i.e., in MG8 we can still observe 70% of the flux from EPG behind the disc while in MG11 we cannot see any emission from this EPG. Table 2 summarises the parameters we choose for the mock galaxies.

Mock galaxy test results
We now analyse the mock galaxy datacubes built using the parameters listed in Table 2 with the method described in Sec. 3. We summarise the fitting results in Table 3. Figure 6 shows the MCMC fitting result for MG1. For most parameters, the fit looks robust, except for the strong degeneracy between and , which also occurs in the H data analysis (Marasco et al. 2019), and we will not comment on results of these two parameters for the rest of this section. We notice that the distribution of and have peaks consistent with input values but also show stronger tails towards one side that however do not influence the recovery of these parameters. The values of ℎ, / , , , and (marked with blue lines in Fig. 6) always fall within the 1 range of the best-fit values, showing that the EPG parameters are recovered successfully. The position-velocity (pv) slices of the best-fit model are also consistent with the mock data. We find a value of 0.32 for the flux fraction of EPG in MG1, which is derived as the total flux ratio between the EPG best-fit model and the datacube, overestimating the input value 0.3 by only 5%. In general, we believe the code can reliably recover the EPG parameters in a case where the EPG emission is strong.
The results for other mock galaxies are shown in Appendix A. In MG2-MG7 the code works well for five EPG parameters (ℎ, / , , , and ), whose input values always fall within (or very close to) 1 range of the best-fit values, showing the reliability of our code for datacubes with even very low S/N or low EPG fractions. We also notice that the code tends to overestimate the flux fraction of EPG, which is not surprising given the poor estimation we have for and . For MG8-MG10, the five EPG parameters also fall within (or close to) 1 range, except for / in MG10, which is more than 2 away from the input value. The parameters (especially / , , and ) are not well recovered in MG11, which is the case where we assume no EPG behind the disc is visible. This is an extreme assumption and is not likely to be the case for our galaxies (see the discussion in Sec 6.1). The above results suggest that dust extinction will not influence our analysis on EPG significantly, even in an extreme case where we can observe only 25% emission of the EPG behind the disc.
Our mock galaxy tests are limited to this relatively small sample because the fitting is time consuming. The parameter space we explore may not cover the entire range of EPG properties in real galaxies. Nevertheless, the tests give us confidence in the following statements. (1) Implementing the EPG modelling code on H data is possible, our modifications are reliable, and even moderate dust extinction does not influence our results significantly. (2) The code is still reliable for very low S/N of EPG (e.g. MG3) or when EPG flux fraction is as low as 15%. (3) The code works well for five EPG parameters (ℎ, / , , , and ), whose input values always fall within (or close to) 1 range of the best-fit values. (4) Because of the intrinsic degeneracy between and , the code cannot in general reproduce these two parameters correctly. As a consequence, flux fraction of the EPG is not fully recovered either. But this does not affect the kinematic analysis, which is our prime goal here.

MODELLING THE EXTRAPLANAR GAS OF THE DISKMASS H SAMPLE
We now apply the method introduced in Sec. 3 to determine the properties of the EPG in NGC 3982 and NGC 4152.

NGC 3982
Following the method described in Sec. 3.1, we extract EPG emission from the NGC 3982 data by implementing external and internal masks. These masks, as shown in Fig. 7, seem to effectively separate EPG emission from the disc emission and pure noise. The internal mask shown in Fig. 7 is computed using an empirical mask threshold of 2 × rms noise (see Sec. 3.1). To explore the influence of the mask threshold on the fitting result, we also generate two additional internal masks using 1.5× and 2.5 × rms noise as thresholds. Once the rotation curve (derived from 3D Barolo, see Sec. 3.2 and Fig. 5) and masks are ready, we fit the residual EPG signal. We show our results in Fig. 8. The pv slices from the data and the best-fit model are very consistent. There are some local features, however, which cannot be recovered by our smooth and axisymmetric EPG model intrinsically, for example the bright filament in the pv-slice −0. 1 offset from the minor axis (indicated by a black arrow). We also note the high-velocity part in the central region (pointed out by red arrows), which is possibly due to outflows from the AGN or non-circular motions caused by the nuclear bar (see Knapen et al. 2002). Modelling such features is beyond the scope of this paper, although we discuss their influence on our fitting results in Sec. 6.1. The corner plots look robust, except for the expected degeneracy between and . The best-fit parameters are shown in Table 4. Our results indicate that there is abundant ionised EPG in NGC 3982, with lagging rotation, vertical and radial inflow, and large velocity dispersion 3 . We note the large uncertainties on the radial and vertical velocities and will discuss these below in Sec 6.2. Since our H data   is not flux calibrated (see above Sec 2.1), the absolute flux of EPG is not available; consequently we can not conduct a flux-to-mass conversion and estimate the total mass of the EPG in NGC 3982.
For the other two internal masks, with thresholds 1.5 × rms and 2.5×rms, we also implement the MCMC fitting method. We find that using these thresholds leads to best-fit parameters which are perfectly compatible with those reported in Table 4 (Sec. 6.1).  (2)

NGC 4152
We apply the method discussed in Sec. 5.1 to the NGC 4152 H datacube. The fitting results for NGC 4152 are shown in Fig. 9. We note that, similarly to NGC 3982, there are also central regions (indicated by red arrows) with high-velocity that cannot be reproduced by the EPG model. Another irregular pattern is shown in the pv-slice offset −0. 1 from the major axis, which seems to be a local feature than cannot be reproduced by our EPG model intrinsically (indicated by a black arrow). Best-fit values are shown in Table 4. The best-fit values for NGC 4152 suggest that this galaxy also has ionised EPG with lagging rotation, radial and vertical inflow, and large velocity dispersion.

Limitations of our EPG model
The EPG modelling code introduced in Marasco et al. (2019) allows us to study ionised EPG with 3D modelling for the first time. Limitations of the code are discussed in detail in Sec. 4.1 of Marasco et al. (2019). Below we discuss some limitations and their influences on our analysis of ionised EPG.
The code assumes a smooth and symmetric (both axisymmetric and across the galactic plane) EPG model which is not the case for all galaxies. Indeed, both NGC 3982 and NGC 4152 show local features that cannot be fitted by our EPG model (see Sec. 5). We note a slight trend that the model tends to 'overshoot' the normal EPG emission component, indicating that the model attempts to recover local features. Nevertheless, those local features do not influence the fitting results significantly. Apart from those local features, our data seem to be symmetric and smooth. We have also attempted to exclude some of the local features by hand (i.e., a manual mask) and the fitting results are compatible.
Another assumption made by the model is constant PA and INCL for the galaxy. Galaxies are typically regular within their effective radius (e.g. approximately within the size of their star-forming disc), while warps and other asymmetries are more significant at larger radii. In the DMS sample, only galaxies with regular morphology and kinematics of H data are selected (Bershady et al. 2010). We therefore conclude that constant PA and INCL is a well motivated assumption.
Ionised gas suffers from dust extinction, of course (Osterbrock & Ferland 2006), which, however, is not included in our model. Nevertheless, mock galaxy tests (MG8-MG10) show that our code can recover EPG parameters well for mock datacubes where ionised EPG located behind the disc is only partly (down to = 25%) visible due to extinction. But in an extreme case (MG11), where gas behind the disc is completely invisible, EPG parameters are not well recovered. Luminous star-forming galaxies like NGC 3982 and NGC 4152 typically have extinction < 1.0 (Calzetti 2001) in most regions, which corresponding to ∼ 40%. We therefore conclude that our results are reliable for these two galaxies. When spatially resolved dust extinction is available (e.g. from H /H line ratios) in future studies, they could be included in our model. An important step in our method is to separate the EPG and the disc gas, for which we use an empirical threshold (2 × rms noise) to mask out disc emission. In principle mask threshold could influence the EPG residuals and thus the fitting results. Our experiment in Sec. 5.1 shows that the fitting results is not sensitive to the mask threshold.

Reliability of the fitting results
To characterise the EPG in our galaxies we fit a model with seven independent parameters to the data. Given the number of parameters and the quality of our data, one may question the reliability of our results. The mock galaxy tests in Sec. 4 provide some insights into this question. Comparing the S/N, EPG parameters and flux fractions of NGC 3982 and NGC 4152 with mock galaxies, we find that they have comparable S/N, flux fraction, and EPG parameters. We conclude that these mock galaxies MG1-MG10 can well represent NGC 3982 and NGC 4152. We do not include MG11, since it assumes an extreme case of (complete) dust extinction and is therefore not comparable with our sample.
From the mock galaxy tests we find that our method does not recover and well due to the intrinsic degeneracy between these two parameters. We therefore should not trust the and values  inferred for NGC 3982 and NGC 4152. In the following we will therefore not discuss these two parameters.
The reliability of the estimated EPG scale height ℎ is influenced by the inclination (INCL) of the galaxy. For example, it is easier to measure the scale height of EPG in an edge-on galaxy than in a galaxy with intermediate inclination. The estimated INCL of NGC 3982 is 26. • 2, which raises the question if our estimate for scale height ℎ is still accurate at such low inclination. With the same INCL (except MG7) as NGC 3982, the scale heights ℎ of MG1-MG10 are recovered correctly (within 1 range of best-fit value in MG1-MG9, and only 1.3 away in MG10), suggesting that our scale height ℎ values for NGC 3982 (even with this low inclination angle) and NGC 4152 are reliable.
Our method also recovers the velocity gradient / well. In most cases the differences between best-fit values and input values of / are less than (or close to) 1 . We note that the code overestimates the velocity gradient by more than 1.5 in MG9 and MG10, where we introduce dust extinction. We conclude that NGC 3982 and NGC 4152 undoubtedly show a lagging rotation, but perhaps with a smaller velocity gradient than suggested by the best-fit values if dust extinction in their discs are significant.
As for the radial velocity and vertical velocity , our method manages to recover them correctly (within or very close to 1 ). Although these are successfully recovered in mock galaxy tests, we note that the error bars for and in NGC 3982 and NGC 4152 are relatively large, even drifting towards zero velocity. From the corner plots in Figs. 8 and 9, we find that the peaks of and are clearly located at the negative value side, and it only is the small tail to the positive side that results in the large error-bar. Such tails also appear in mock galaxy tests (see, e.g., Fig. 6). We conclude that the code indeed prefers radial and vertical inflows for NGC 3982 and NGC 4152, but with these large error-bars, we also cannot completely exclude the possibility of zero radial and/or vertical EPG velocities.

Comparison with previous work
We now compare our ionised EPG parameters to those previously measured in the literature. Rossa & Dettmar (2003a) detected ionised EPG in 30 edge-on galaxies with photometric observations and found a typical scale height ℎ of 1-2 kpc. Another photometric study conducted by Miller & Veilleux (2003a) reported scale heights of ionised EPG varying from 0.4 kpc to 17.9 kpc for 16 edge-on galaxies, with a mean value of 4.3 kpc. This relatively high value of scale height is likely due to the deep flux limit in their data. Using a MaNGA sample of 65 edge-on galaxies, Bizyaev et al. (2017) derived a median scale height of 1.2 ± 0.5 kpc for ionised EPG. The EDGE-CALIFA Survey used a sample of 25 edge-on galaxies and measured a median scale height of 0.8 +0.7 −0.4 kpc (Levy et al. 2019). They also summarised a list of previous studies on ionised EPG and found a median scale height of 1.0 ± 2.2 kpc. The scale heights we derive for NGC 3982 (0.83 +0.27 −0.40 kpc) and NGC 4152 (1.87 +0.43 −0.56 kpc) are consistent with these studies.
The velocity gradients of ionised EPG was also studied in several spectroscopic studies. For example, Heald et al. (2005) reported a velocity gradient of −15 km s −1 kpc −1 for ionised EPG in NGC 891 and −8 km s −1 kpc −1 for NGC 5775. And the EDGE-CALIFA survey found a velocity gradient of −21 km s −1 kpc −1 (Levy et al. 2019). The velocity gradients of NGC 3982 ( / = −22.24 +6.60 −13.13 km s −1 kpc −1 ) and NGC 4152 ( / = −11.18 +3.49 −4.06 km s −1 kpc −1 ) are also consistent with previous studies. Those previous measurements of EPG properties are from edgeon galaxies where the scale height can be directly measured and velocity gradients are calculated by analysing velocity profiles. It is noteworthy that our results for scale height and velocity gradient are very consistent with these studies although we use a completely different method, which is 3D modelling of the datacube. This method allows us to study a much broader sample of galaxies, as geometrical constraints become significantly less important.
Our analysis for NGC 3982 and NGC 4152 suggests that ionised EPG shows radial and vertical inflow, which has also been found in previous studies on other galaxies. The ionised EPG in NGC 2403 shows evidence of localised outflows, but the inflow seems to be diffused over the scales of the entire star-forming disc (Fraternali et al. 2004). Zheng et al. (2017) found disc-wide ionised gas inflow in M33 via absorption line study against bright UV sources in the disc. Warm gas in the Milky Way traced by absorption lines has more detections at negative than positive velocities (Lehner et al. 2012; see also left panel of Fig.4 in Marasco et al. 2013), indicating that inflow is more prominent than outflow for ionised EPG in the Milky Way. Other previous studies were limited to edge-on samples as mentioned above and lack simultaneous vertical and radial velocity analysis. But the inflow of H EPG has been found in many studies (e.g. Fraternali et al. 2001;Marasco et al. 2019).

Is the EPG produced by a galactic fountain?
In the classic galactic fountain scenario, the EPG originates from the disc, is pushed up by star-formation activity (stellar wind and supernova explosions), and eventually falls back to the disc (e.g. Bregman 1980). Therefore the gas is outflowing in the first stage and inflowing in the second stage. However, H EPG is usually found to be preferentially inflowing in observations. A plausible explanation is that the outflowing gas might be highly ionised and only cools down near the apocentre, thus H is only partly visible during the ascending phase and fully visible during the descending phase (e.g. Fraternali & Binney 2006Marasco et al. 2012). In the above mentioned scenario, we would expect to see outflows to dominate the ionised gas kinematics, contrary to what we infer for the ionised EPG in NGC 3982 and NGC 4152.
In what follows we discuss the possible explanations for this discrepancy, either within the galactic fountain framework or outside of it. One possible mechanism to explain this discrepancy is that the outflowing gas, being highly ionised, is at a temperature of around 10 5 K or more, in which case the gas is too hot to be seen in H but may be visible through the absorption of O or even X-ray emission (e.g. Fraternali et al. 2002b). Later the hot gas cools down and turns into warm gas and neutral gas. Furthermore, it is likely that part of the cooling H gas is ionised by photons leaking from the disc again during the descending phase, and emits in H , which can explain why we see similar kinematic properties for the ionised (H ) and neutral (H ) EPG.
In Sec. 1 we mention that a model of pure galactic fountain cannot fully explain the properties of EPG, in particular its slowly rotating and inflowing kinematics require fountain-driven condensation of the hot corona (Fraternali 2017). Hydrodynamic simulations are required to fully capture the physics of this mechanism (Marinacci et al. 2011;Armillotta et al. 2016;Grønnow et al. 2018). When fountain gas moves through the hot corona, it mixes with the hot gas. The mixture increases the metallicity and decreases the temperature of the hot corona and therefore decreases its cooling time, making part of the hot corona condense into warm or cold gas. A sketch of this process can be seen in Fig. 1 of Fraternali et al. (2013). It is possible that this accretion mechanism can also explain the ionised EPG inflow found in our data, since the amount of warm EPG increases due to condensation in the descending part (also see Fig. 1 of Fraternali et al. 2013 and Fig. 1 of Marasco et al. 2013).
While our model assigns a unique value to the inflow velocity of the EPG, the galactic fountain kinematics can be very complex. We stress that a preference for inflow does not mean that there is no outflow. It is likely that the inflow of ionised EPG due to condensation of the corona together with the H EPG being ionised in the descending phase leads to our models to prefer inflow for H EPG. Moreover, as mentioned in Sec. 6.3, outflowing gas is likely to be more localised (e.g. Fraternali et al. 2004;Cicone et al. 2016;Boettcher et al. 2017) since it comes from H regions, which are clumpy and preferentially confined within the spiral arms. As fountain clouds pass through the halo they fragment and expand (Marinacci et al. 2011), eventually producing a smooth and diffused rain of multi-phase gas onto the disc. In the above scenario, our model, which is based on the largescale smooth emission, would also prefer the inflow. It is, however, not possible to test the above hypothesis with our current kinematic model. To better resolve this issue, a dynamical model for ionised EPG data is required and it will be subject of future investigations.
Another possible explanation for the observed inflow is that of an external origin for the ionised EPG. In this scenario the H -emitting material would not originate from the circulation of fountain clouds, but it would rather trace gas accretion from the CGM (Kaufmann et al. 2006;Zheng et al. 2017) produced either by the spontaneous cooling of the hot corona (e.g. Voit et al. 2015;Sormani & Sobacchi 2019) or by cold filaments penetrating the hot gas down to the disc (e.g. Fernández et al. 2012;Mandelker et al. 2019). However, we stress that the dominating motion in our EPG is not the inflow but the rotation, which has remarkable proximity to the disc with only a mild vertical gradient. This indicates that the accreting gas close to the disc must have a relatively high angular momentum (Pezzulli & Fraternali 2016), and/or be efficiently accelerated by the galactic fountain (Marinacci et al. 2011).

CONCLUSION
In this paper we have studied the kinematic properties of ionised EPG in two galaxies from the DiskMass survey, NGC 3982 and NGC 4152. We have modelled the H datacubes of these galaxies using the 3D kinematic modelling code from Marasco et al. (2019), which first separates the EPG emission from the disc emission and then fits the EPG residual with a model controlled by seven parameters via a Bayesian MCMC approach. The seven parameters characterising EPG properties are the scale length , an exponential parameter (see Eq. 1), the scale height ℎ, the velocity gradient / , the radial velocity , the vertical velocity , and the velocity dispersion . We have tested the reliability of our method with mock galaxies, which we have built assuming ionised gas properties compatible with those of NGC 3982 and NGC 4152.
Our results for NGC 3982 and NGC 4152 can be summarised as follows.
(i) We have detected ionised EPG in both galaxies, with flux fraction between the EPG and the regularly rotating gas within the disc of 27% (NGC 3982) and 15% (NGC 4152).
(iii) The ionised EPG in both galaxies show lagging rotation, with rotational gradients of −22.24 +6.60 −13.13 km s −1 kpc −1 (NGC 3982) and −11.18 +3.49 −4.06 km s −1 kpc −1 (NGC 4152). (iv) Our best-fit model suggests radial and vertical inflow for both galaxies. But we cannot exclude the possibility of zero radial and/or vertical velocity as the uncertainties for these two parameters are large.
The kinematic signature that we have derived is consistent with the EPG layer being mostly produced by galactic fountains, but gas accretion, either induced or spontaneous, may also play an important role. A dynamical model of EPG is needed to better understand the origin of the ionised EPG. A large sample of galaxies with better spatial resolution is also essential for further studies of ionised EPG. In the near future, WEAVE will provide IFU data for a large sample of late-type galaxies at a comparable spectral resolution, a wider wavelength range, and a higher angular resolution, which will enable a more comprehensive study of the ionised EPG of disc galaxies.