A Virgo Environmental Survey Tracing Ionised Gas Emission (VESTIGE). VIII. Modeling ram pressure stripping of diffuse gas in the Virgo cluster spiral galaxy NGC 4330

NGC 4330 is one of the Virgo galaxies whose UV distribution shows a tail structure. An associated tail structure is also observed in the HI and H$\alpha$ emission distributions. Previous dynamical modeling showed that the galaxy is approaching the cluster center and is therefore undergoing increasing ram pressure stripping. Recent stellar population fitting of deep optical spectra together with multiband photometry lead to the determination of the time when star formation was quenched in the galactic disk. We introduce a new version of the dynamical model that includes the diffuse ionized gas and aim to reproduce the HI, H$\alpha$, UV distributions together with the star formation histories of the outer gas-free parts of the galactic disk. The results of 50 simulations with five different Lorentzian temporal ram-pressure profiles and five different delays between the simulation onset and peak ram pressure are presented. The inclusion of diffuse gas stripping changes significantly the HI, UV, and H$\alpha$ emission distributions. The simulations with diffuse gas stripping naturally lead to vertical low surface density filaments in the downwind region of the galactic disk. These filaments occur less frequently in the simulations without diffuse gas stripping. The simulations with diffuse gas stripping lead to better joint fits to the SEDs and optical spectra. The HI, NUV, and H$\alpha$ morphologies of the model snapshots which best reproduce the SEDs and optical spectra are sufficiently different to permit a selection of best-fit models. We conclude that the inclusion of diffuse gas stripping significantly improves the resemblance between the model and observations. Our preferred model yields a time to peak ram pressure of 140 Myr in the future. The spatial coincidence of the radio continuum and diffuse H$\alpha$ tails suggests that both gas phases are stripped together.


Introduction
Once a galaxy has entered a galaxy cluster, its evolution will change, sometimes drastically. Multiple gravitational interactions together with the influence of the cluster potential can remove the stellar and gaseous content of its outer disk (galaxy harassment; Moore et al. 1996). The rapid motion of the galaxy within the cluster atmosphere, the hot (temperature T ∼ 10 7 K) and tenuous (density n ∼ 10 −4 cm −3 ) intracluster medium (ICM), causes the removal of the outer gas disk via ram pressure. In contrast to galaxy harassment, ram pressure stripping does not affect the stellar content of the galaxy. Spiral galaxies which underwent or undergo ram pressure stripping show a truncated gas disk together with a symmetric stellar disk (e.g., VIVA Chung et al. 2009). If the interaction is ongoing, a gas tail mainly detected in Hi is present (e.g., Chung et al. 2007). The gas truncation 2 Vollmer et al.: Modeling ram pressure stripping of diffuse gas in NGC 4330 radius is set by the galaxy's closest passage to the cluster center via the criterion introduced by Gunn & Gott (1972): where ρ ICM is the ICM density, v gal the galaxy velocity, Σ the surface density of the interstellar medium (ISM), v rot the rotation velocity of the galaxy, and R the stripping radius. If peak ram pressure occurred more than 400-500 Myr ago, the gas tails have disappeared and the truncated gas disk has become symmetric again (Vollmer 2009). Under extreme conditions, ram-pressure stripped galaxies often show extraplanar, one-sided optical and UV emission and important tails of ionized gas (e.g., Yagi et al. 2010Yagi et al. , 2017Poggianti et al. 2017). Because of their optical appearance these objects are called jellyfish-galaxies. Most of the Hα emission in these tails is due to photoionization by massive stars born in situ in the tails (Poggianti et al. 2019). Often extraplanar molecular gas is found in jellyfish galaxies (Jachym et al. 2014(Jachym et al. , 2019Moretti et al. 2018), from which the ionizing stars are formed.
Once the gas has been removed from the outer galactic disk under the action of ram pressure, the corresponding stellar content remains gas and dust-free and gradually stops star formation. Optical spectroscopy together with multiwavelength photometry can be used to determine the stellar populations and thus the star formation history of the ram-pressure stripped parts of the galactic disk (Boselli et al. 2006, Crowl et al. 2008. These star formation histories can then be compared to those of dynamical models. For the Virgo cluster galaxy NGC 4388 an observed quenching timescale of ∼ 200 Myr was determined by Pappalardo et al. (2010) based on a VLT spectrum and optical/UV photometry. An updated dynamical model including ram pressure stripping was able to reproduce the observed quenching timescale together with additional observational characteristics (Vollmer et al. 2018): Hi distribution and velocity field, FUV, Hα, and polarized radio continuum emission distribution. A realistic model has to fulfill all these observational constraints. From the model the ram pressure temporal profile, time to peak ram pressure, and the angle between the galaxy velocity within the cluster and its disk plane can be determined. A model-based time series of ram pressure-stripped galaxies was established by Vollmer (2009). The Virgo spiral galaxy NGC 4330 is part of this time series.
The interstellar medium (ISM) of a spiral galaxy consists of different phases: (i) cool dense molecular gas (T ∼ 10-30 K, n > ∼ 100 cm −1 ), (ii) cold neutral hydrogen (T ∼ 100 K, n ∼ 100 cm −1 ), (iii) warm neutral hydrogen (T ∼ 8000 K, n ∼ 1 cm −1 ), and (iv) warm diffuse ionized hydrogen (T ∼ 8000 K, n ∼ 0.4 cm −1 ). Since the acceleration of a gas cloud by ram pressure is inversely proportional to the gas surface density (Eq. 1), it is expected that the different phases are stripped with different efficiencies. Indeed, the stripped diffuse ionized gas is accelerated to higher velocities and is therefore observed at larger distances from the host galaxy than the dense gas (Vollmer et al. 2009, Moretti et al. 2018). On the other hand, the stripped dense CO-emitting gas seems to be well mixed with the warm Hi (Vollmer et al. 2008Jáchym et al. 2014;Nehlig et al. 2016), which is consistent with a scenario where the neutral turbulent ISM is stripped as an entity. We therefore expect that the distribution of the stripped diffuse ionized gas is more extended than that of the neutral gas. The detection of this gas in deep Hα observations represents an important additional constraint on simulations of ram-pressure stripped galaxies. An additional constraint on the simulations comes from the star formation history of stripped gas-free regions of the stellar disk. Detailed simulations of a ram-pressure stripped galaxy can thus not only give access to the interaction parameters (see Vollmer 2009), but also give insight into the reaction of the different gas phases to ram pressure.
In this article we (i) test a new version of the dynamical model that includes not only the dense neutral gas, but also the diffuse ionized gas and apply it to NGC 4330 and (ii) aim to reproduce the observed Hi, Hα, NUV distributions of NGC 4330 together with the star formation histories of the outer gas-free parts of the galactic disk. This article represents a continuation and update of the work of .

NGC 4330
NGC 4330 is one of the Virgo spiral galaxies with a long Hi tail observed by Chung et al. (2007). It has a maximum rotation velocity of v rot ∼ 180 km s −1 and a total Hi mass of M HI = 4.5 × 10 8 M ). NGC 4330 is located at a projected distance of ∼ 2 • (600 kpc) from the cluster center, i.e. it is relatively close to M 87, and has a radial velocity of 400 km s −1 with respect to the Virgo cluster mean. The Hi deficiency of NGC 4330, which is defined as the logarithm of the ratio of the Hi content of a field galaxy of same morphological type and diameter to the observed Hi mass, is 0.8 (Chung et al. 2007), i.e. the galaxy has lost about 80 % of its atomic hydrogen. The Hi distribution in the galactic disk (upper left panel of Fig. 1) is truncated at about half the optical radius. In addition, NGC 4330 is one of the rare Virgo galaxies showing an extended UV tail (Abramson et al. 2011). The Hi and UV tails show a significant offset, with the Hi tail being downwind of the UV tail. At the leading edge of the interaction, the Hα emission and dust extinction distribution is bent sharply out of the galactic disk. These features are signs of active ram pressure stripping. Roediger & Hensler (2005) elaborated the subsequent stripping phases for face-on ram pressure stripping. NGC 4330 undergoes the instantaneous stripping phase where the outer part of the gas disk is displaced but only partially unbound. This leads to straightly bent but coherent outer gas disk.  presented a detailed dynamical model for NGC 4330. Their best-fit model qualitatively reproduced the observed projected position, the radial velocity of the galaxy, the molecular and atomic gas distribution and velocity field, and the UV distribution in the region where a gas tail is present. In contrast to other Virgo spiral galaxies affected by ram pressure stripping, NGC 4330 does not show an asymmetric ridge of polarized radio continuum emission. This is due to the particular projection of NGC 4330 where the compressed large-scale magnetic fields are oriented along the light-of-sight.  concluded that the ISM is stripped as a whole entity and that only a tiny fraction of dense clouds can be left behind. Extraplanar star formation proceeds in dense gas arms pushed by ram pressure. The collapsing clouds decouple rapidly from the ram pressure wind and the young, UV-emitting stars stay behind the gas. On the basis of dynamical models, the galaxy moves to the north and still approaches the cluster center with the closest approach occurring in ∼ 100 Myr. Despite the success of their model,  noted that the observed red UV color on the windward side was not reproduced by the model. Moreover, their simulations could only produce the model distribution of Hii regions, i.e. compact regions of high ionized gas density and high Hα surface brightness ( Fig. 9 of . Since the diffuse ionized gas is much less dense than the gas in the Hii regions, it has a much lower surface brightness and can only be detected through deep narrow-band observations. Such observations for NGC 4330 are only available now. Therefore, the model of  was upgraded to enable the simulation of diffuse ionized gas stripping. Within the Virgo Environmental Survey Tracing Ionised Gas Emission (VESTIGE; Boselli et al. 2018 detected a low surface brightness 10 kpc Hα tail exhibiting a peculiar filamentary structure. In addition, these authors collected literature photometry in 15 bands from the far-UV to the far-IR and VLT FORS2 deep optical longslit spectroscopy. Using a Monte Carlo code that jointly fits spectroscopy and photometry, they reconstructed the star formation histories in apertures along the major axis of NGC 4330. In this way they found a clear outside-in gradient with radius of the time when the ram pressure induced star formation quenching started: the outermost radii were stripped about 500 Myr ago, while the stripping reached the inner 5 kpc from the center in the last 100 Myr. In addition, Fossati et al. (2018) discovered a low surface brightness 10 kpc tail of ionized gas exhibiting a peculiar filamentary structure (upper right panel of Fig. 1). The ionized gas tail is associated to the Hi tail, but extends further to the south.

The dynamical model
We used the N-body code described in Vollmer et al. (2001), which consists of two components: a non-collisional component that simulates the stellar bulge/disk and the dark halo, and a collisional component that simulates the ISM as an ensemble of gas clouds. A scheme for star formation was implemented, where stars were formed during cloud collisions and then evolved as non-collisional particles.

Halo, stars, and gas
The non-collisional component consists of 81 920 particles, which simulate the galactic halo, bulge, and disk. The characteristics of the different galactic components are shown in Table 1. The resulting rotation velocity is ∼180 km s −1 and the rotation curve becomes flat at a radius of about 5 kpc. We adopted a model where the ISM is simulated as a collisional component, i.e. as discrete particles that possess a mass and a radius and can have inelastic collisions (sticky particles). The advantage of our approach is that ram pressure can be easily included as an additional acceleration on particles that are not protected by other particles (see Vollmer et al. 2001).
The 20 000 particles of the collisional component represent gas cloud complexes that evolve in the gravitational potential of the galaxy. The total assumed gas mass is M tot gas = 4.3 × 10 9 M , which corresponds to the total neutral gas mass before stripping, i.e. to an Hi deficiency of 0. A radius is attributed to each particle which depends on its mass assuming a constant surface density. The normalization of the mass-size relation was taken from Vollmer et al. (2012a). During each cloud-cloud collision the overlapping parts of the clouds are calculated. Let b be impact parameter and r 1 and r 2 the radii of the larger and smaller clouds. If r 1 + r 2 > b > r 1 − r 2 the collision can result into fragmentation (high-speed encounter) or mass exchange. If b < r 1 − r 2 mass exchange or coalescence (low speed encounter) can occur. The outcome of the collision is simplified following Wiegel (1994). If the maximum number of gas particles/cloud (20000) is reached, only coalescent or mass exchanging collisions are allowed. In this way a cloud mass distribution is naturally produced. The cloud-cloud collisions result in an effective gas viscosity in the disk.
As the galaxy moves through the ICM, its clouds are accelerated by ram pressure a = ρ ICM v 2 gal /Σ. In addition, the gas clouds are accelerated by the gradients of the gravitational potential a = −∇φ. The ISM surface density is assumed to be Σ(r) = (1 + exp(−r/2 kpc)) × 10 M pc −2 . The constant surface density at the outer disk is motivated by the observed saturation of the Hi surface density (e.g., Leroy et al. 2008). The radial dependence takes into account the increased surface density due to the presence of molecular gas. Within the galaxy's inertial system, its clouds are exposed to a wind coming from a direction opposite to that of the galaxy's motion through the ICM. The temporal ram pressure profile has the form of a Lorentzian, which is realistic for galaxies on highly eccentric orbits within the Virgo cluster (Vollmer et al. 2001). The effect of ram pressure on the clouds is simulated by an additional force on the clouds in the wind direction. Only clouds that are not protected by other clouds against the wind are affected. This results in a finite penetration length of the intracluster medium into the ISM (see Appendix A). Since the gas cannot develop instabilities, the influence of turbulence on the stripped gas is not included in the model. The mixing of the intracluster medium into the ISM is very crudely approximated by the finite penetration length of the intracluster medium into the ISM, i.e. up to this penetration length the clouds undergo an additional acceleration caused by ram pressure.

Star formation
We assume that the star formation rate is proportional to the cloud collision rate. A mass-size relation based on a constant gas surface density is assumed for the collisional particles. The collisions are explicitly calculated for each timestep (∼ 10 4 yr) via the simplified geometrical recipes of Wiegel (1994). We verified that the star formation rate of an unperturbed galaxy is constant within 1 Gyr.  Vollmer & Beckert (2003) elaborated an analytical model for a galactic gas disk which considers the warm, cold, and molecular phases of the ISM as a single, turbulent gas. This gas is assumed to be in vertical hydrostatic equilibrium, with the midplane pressure balancing the weight of the gas and stellar disk. The gas is taken to be clumpy, so that the local density is enhanced relative to the average density of the disk. Using this local density, the free-fall time of an individual gas clump, the governing timescale for star formation, can be determined. The star formation rate is used to calculate the rate of energy injection by supernovae. This rate is related to the turbulent velocity dispersion and the driving scale of turbulence. These quantities in turn provide estimates of the clumpiness of gas in the disk (i.e., the contrast between local and average den-sity) and the rate at which viscosity moves matter inward. Vollmer & Leroy (2011) applied the model successfully to a sample of local spiral galaxies. Within the framework of the model the star formation rate per unit volume is given where Φ V and t −1 ff, cl are the volume filling factor and freefall time of selfgravitating gas clouds and ρ is the overall gas density. The volume filling factor is Φ V = n cl 4π/3 r 3 cl , where n cl is the internal number density of clouds and r cl is the cloud radius. For a selfgravitating cloud the free-fall time equals the turbulent crossing time t ff,cl = 2 r cl /v turb,cl . With the collision time given by t coll = (n cl π r 2 cl v turb ) −1 , this leads toρ * = 2/3 ρ t −1 coll . Thus, the star formation rates of the two recipes are formally equivalent.
In numerical simulations, the star formation recipe usually involves the gas density ρ and free-fall time t ff = 3 π/(32 Gρ):ρ * ∝ ρ t −1 ff ∝ ρ 1.5 . We verified that our star formation recipe based on cloud-cloud collisions leads to the same exponent of the gas density in a simulation of an isolated spiral galaxy. As a consequence, our code reproduces the observed SFR-total gas surface density, SFR-molecular gas surface density, and SFR-stellar surface density relations .
During the simulations, stars are formed in cloud-cloud collisions. At each collision, a collisionless particle is created, which is added to the ensemble of collisional and collisionless particles. The newly created collisionless particles have zero mass (they are test particles) and the positions and velocities of the colliding clouds after the collision. These particles are then evolved passively with the whole system. Since in our sticky-particle scheme there is mass exchange, coalescence, or fragmentation at the end of a collision, the same clouds do not collide infinitely. The local collision rate traces the cloud density and the velocity dispersion of the collisional component.
The information about the time of creation is attached to each newly created star particle. In this way, the Hα emission distribution caused by Hii regions can be modeled by the distribution of star particles with ages younger than 20 Myr. The UV emission of a star particle in the two GALEX bands is modeled by the UV flux from single stellar population models from STARBURST99 (Leitherer et al. 1999). The age of the stellar population equals the time since the creation of the star particle. The total UV distribution is then the extinction-free distribution of the UV emission of the newly created star particles.

Diffuse gas stripping
Our numerical code is not able to treat diffuse gas in a consistent way. For a realistic treatment 3D hydrodynamical simulations should be adopted. Nevertheless, we can mimic the action of ram pressure on diffuse gas by applying very simple recipes based on the fact that the acceleration by ram pressure is inversely proportional to the gas surface density which in turn depends on the gas density for a gas cloud of constant mass. For simplicity we do not modify the radii of the diffuse gas clouds for the calculation of the cloud-cloud collisions. The diffuse clouds are thus mostly ballistic particles under the influence of a ram-pressure induced acceleration. The recipes only take into account gas cooling of the hot phase.
We assume that the warm (∼ 10 4 -10 5 K) gas clouds become diffuse, i.e. their sizes and volume filling factor increase and their densities and column densities decrease, if it is stripped out of the galactic plane and if its density falls below the critical density of n warm crit = 5 × 10 −3 cm −3 . We further assume that the first condition is fulfilled if the stellar density at the location of the gas particle is lower than ρ crit * = 2.5 × 10 −4 M pc −3 . For our model galaxy this density is reached at a disk height of ∼ 3.5-4 kpc.
Once the ionized gas has a high volume filling factor, i.e. it becomes diffuse, its volume and surface densities decrease. The decrease of the surface density is taken into account by considering gas cloud of constant mass: for ρ * < ρ crit * the cloud size is proportional to n −1/3 , the cloud surface density is proportional to n 2/3 , and the acceleration caused by ram pressure a = ρ ICM v 2 gal /Σ is increased by a factor (0.044 cm −3 /n ISM ) 2/3 . The latter normalization and the critical densities were chosen such that they led to acceptable results for several observed galaxies undergoing ram pressure stripping. The critical stellar density is motivated by the fact that a significant change of the ISM properties only occurs once the ISM has entirely left the galactic disk. It turned out that this condition is necessary to avoid excessive gas stripping. The gas and stellar densities are calculated via the 50 nearest neighbouring particles.
The hot (> 10 6 K) diffuse gas is taken into account in the following way: we assume that once the stripped gas has left the galactic disk, it is mixing with the ambient ICM. The initial temperature of the ISM is T ISM = 10 4 K. If (i) the gas density falls below the critical value of n hot crit = 5×10 −4 cm −3 , (ii) the stellar density is below ρ crit * , and (iii) the ISM temperature is below 9 × 10 6 K, ICM-ISM mixing begins to act. On the other hand, if the stellar density exceeds its critical value, i.e. the gas is located within or close to the galactic disk, the gas is assumed to cool rapidly and its temperature is set to T ISM = 10 5 K 1 . We assume that mixing occurs instantaneously and raises the temperature of the mixed gas clouds to where n ICM is the density of the mixed ICM which is calculated for each particle. If an exceedingly high mixed ICM density (n ICM > 2 × 10 −3 cm −3 ) is needed to heat the gas cloud to 9 × 10 6 K, ICM-ISM heating is not allowed. The temperature of the mixed hot gas clouds T ISM , if not heated, decreases with a cooling timescale of t cool = 1.5 × 10 7 n −1 ISM yr, where the ISM density n ISM is given in cm −3 . Once mixed, the gas cloud is not allowed to cool below a temperature of 10 6 K 2 . This algorithm ensures that stripped gas below the critical density n crit is kept at a temperature of 9 × 10 6 K, whereas denser gas can cool to a temperature of 10 6 K.
For the stripped hot gas clouds (> 10 6 K) we modified the equation of motion by adding the acceleration caused by the gas pressure gradient a = −∇p ISM /ρ ISM using an SPH formalism. For simplicity an isothermal stripped ISM with a sound speed of c s = 235 km s −1 is assumed for this purpose. We thus do not solve an explicit energy equation and the calculation of hydrodynamic effects is approximate. Its main purpose is the avoiding of clumping of the hot diffuse ISM.
For the calculation of the acceleration caused by ram pressure a = ρ ICM v 2 gal /Σ, the acceleration is further increased by a heuristic factor of 10 if the ISM temperature exceeds 10 5 K and the stellar density is below ρ * = 2.5 × 10 −4 M pc −3 . We note that without the inclusion of hydrodynamical effects thin gas filaments cannot be pro-duced by our simple model. The diffuse gas is assumed to emit in the Hα line if its temperature is below 10 6 K.

Parameters of the ram pressure stripping event
Following Vollmer et al. (2001) we used Lorentzian profile for the time evolution of ram pressure stripping: where p max is the maximum ram pressure occurring at the galaxy's closest passage to the cluster center, t HW is the width of the profile, and t peak is the time of peak ram pressure. The simulations were calculated from t = 0 to t = 800 Myr. We set p max = 5000, 6000, 10000 cm −3 (km s −1 ) 2 and t HW = 100, 200, 300 Myr. Following Nehlig et al. (2016) and Vollmer et al. (2018), we investigated the influence of galactic structure, i.e. the position of spiral arms, on the results of ram pressure stripping by varying the time of peak ram pressure between t peak = 570 Myr and 690 Myr. The ram pressure profiles are specified in Table 2   We adopted the  values for the angle between the ram pressure wind and the disk plane (75 • ). The inclination of the model galaxy is i = 90 • to the lineof-sight.

Observations and model fitting
The deep optical spectrum was obtained with the FOcal Reducer and low dispersion Spectrograph (FORS2; Appenzeller et al. 1998) mounted on UT1 of the ESO Very Large Telescope. With a length of 6.8 , the FORS2 slit covered the full extent of the galaxy. A slit aperture of 1.3 was used to optimize the covered area without compromising the spectral resolution of R ∼ 1000. Six exposures of 900 s each using the 600B+22 grism which covers the wavelength range 3300-6210Å, and three exposures of 900 s each using the 600RI+19 grism covering the wavelength range 5120-8450Å. The deep continuum subtracted Hα+[Nii] observations were obtained within the VESTIGE survey ) using MegaCam at the CFHT with the NB filter (MP9603) and the broad-band r filter (MP9602). The ancillary GALEX, CFHT, Calar Alto 2.2m telescope, Spitzer, and Herschel photometric data are described in . In addition, for the comparison with our models we use VIVA Hi and 20 cm continuum data ) and 6 cm radio continuum data from Vollmer et al. (2013). Fossati et al. (2018) reconstructed the quenching histories of NGC 4330 in various regions along the major axis by means of stellar population fitting. They performed a joint fitting of spectra and photometry using Monte Carlo techniques to optimally explore the parameter space. The emission of the galactic disk was divided into 13 apertures along the major axis of NGC 4330, as shown in Fig. 3. Aperture photometry was obtained by summing the calibrated flux values in each pixel contributing to a given aperture minus the background value in the same aperture. To estimate the background each aperture was randomly placed in empty regions of the images and the sum of the counts was taken. This procedure was repeated 1000 times and the median value was adopted as the best estimator for the background value. For the UV, optical, and near infrared bands, the measured fluxes were corrected for the Galactic attenuation. The spectra from the FORS2 reduced 2D spectral images were extracted by summing the flux in the spatial pixels contributing to each region along the slit direction. For the reconstruction of the quenching history of each radial bin a model of its unperturbed radial star formation history (SFH) is needed. For this purpose Fossati et al. (2018) adopted the multizone models for the chemical and spectrophotometric unperturbed disc evolution originally presented in Boissier & Prantzos (2000) and later updated by Muñoz-Mateos et al. (2011). The model was projected to fit the high inclination of the disc of NGC 4330 and the SFH were derived in each of the apertures. Fossati et al. (2018) further modified the unperturbed SFHs by introducing an exponential decline of the SFR to parametrize the quenching event.
Parametric models to the spectro-photometric data were produced using MULTINEST (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2013. A grid of stellar spectra was constructed using the quenched SFHs coupled to Bruzual & Charlot (2003) high-resolution models. Nebular emission lines were added to the stellar template. To do so  first computed the number of Lyman continuum (λ < 912Å) photons Q H0 for a given stellar spectrum. They then assumed that all the ionizing radiation is absorbed by the gas and that the ionizing photons do not contribute to the heating of dust. However, they modeled the uncertainty in the conversion of Q H0 into emission line flux by adding a nuisance parameter (Ly scale ) in their MC-SPF code. This parameter also takes into account differential aperture effects between the 1.3 spectral slit and the larger photometric apertures. A double Calzetti et al. (2000) attenuation law was assumed to include the effects of dust attenuation on the stellar and nebular emission line spectrum. Fossati et al. (2018) coupled the stellar spectrum with dust emission models from the mid-to the far-infrared by means of an energy balance between the stellar flux absorbed by the dust and re-emitted in the infrared. They jointly fitted the photometric data points and the FORS2 spectra at their native resolution. The multidimensional likelihood space was sampled by using PYMULTINEST (Buchner et al. 2014), a python wrapper for the MULTINEST code. For direct comparison with the results obtained from our model SFHs, we reproduce the results of Fossati et al. (2018) in Fig. 10 to 15.
For our purpose we use the modeling of Fossati et al. (2018) to obtain the synthetic photometry and the optical spectra for our model SFHs. The model photometry and spectrum is then compared to observations using MULTINEST and the total evidence is calculated. The goodness is evaluated using the Bayesian evidence. This value is formally defined as the integral of the likelihood of a given model over the entire parameter space. When two models are compared, their goodness in describing the data can be summarized in the ratio of their evidences (or the difference in log space). The Bayesian evidence automatically penalises models with more free parameters thus optimally highlighting the best model even in case the models under considerations have a different number of free parameters. In our case, these free parameters are the temporal ram pressure profile, the chosen timestep, and the inclusion of diffuse gas stripping. The variation of parameters which are used to compare the model to the data (stellar libraries, attenuation model, Ly scale parameter) can change the goodnesses and thus the order of the preferred models. A better model yields a higher goodness.

Results
We calculated the model photometry and spectrum for all timesteps (∆t = 10 Myr) of the 50 simulations presented in Table 2 with (format Nanew) and without (format Na) a diffuse gas component. In addition, we decided to handle the Ly scale nuisance parameter in two different ways: (i) we set it to the fixed values found by Fossati et al. (2018) and (ii) we treated it as a free parameter and kept the value that lead to the best fit to observations. The distributions of the total goodnesses for all model timesteps are shown in Fig. B.1. The distributions show a tail to high values containing about 10-20 timesteps. The resulting total goodnesses are presented in Table 3. The last number of the model names is the simulation timestep N . The time to peak ram pressure is ((N − 1) × 10) − t peak Myr, where t peak is given in Table 2.
We then searched for the 12 timesteps with the highest goodnesses for the downturn, tail, and total (down-turn+tail) regions. These different regions are generally not best fit by a single timestep (Table 4). Fig. B.2 compared to Fig. 12 gives an impression of the meaning of the total goodness g: the model and observed spectra and SEDs of region 3 for two timesteps (models 4new 47 and 1cnew 46) with a goodness difference of ∆g ∼ 200 show that the FUV and NUV fluxes are significantly better reproduced by model 4new 47. We therefore presume that goodness differences in excess of ∆g > ∼ 200 are significant. All of our total model goodnesses are smaller than that of parametric model of Fossati et al. (2018). This is expected because the parametric model had a higher degree of freedom. The models with a diffuse stripped ISM lead to somewhat higher goodnesses than those without a diffuse stripped ISM. Moreover, models with a variable Ly scale parameter systematically yield higher goodnesses thanks to the physical relevance of this parameter in reproducing the observed emission lines in the spectra. This is true for the downturn, tail, and total regions.
Model 4new 47 with diffuse gas stripping has the highest total goodness of all tested models (Table 3). Given that a goodness difference of ∆g > ∼ 200 is significant, models 1cnew 47 (diffuse gas stripping), 4b 67, and 2b 80 are also acceptable amongst the models with a variable Ly scale parameter. For a constant Ly scale parameter, only model 2b 80 is acceptable.
If we look at the goodnesses of the downturn and tail regions separately (Table 4), the results change: in the case of a variable Ly scale parameter all 12 best-fit models of the downturn region are acceptable if diffuse gas stripping is included or not. For the tail region all models without diffuse gas stripping are acceptable. If the stripping of diffuse gas is included in the simulations, only the first three models are acceptable. In the case of a constant Ly scale parameter the four and two best-fit models are acceptable when diffuse gas stripping is included or not. In the tail region only the eight best-fit models without diffuse gas stripping are acceptable.
We conclude that the best-fit models based on the total goodness are pre-peak ram pressure models if diffuse gas stripping is included, otherwise the preferred timesteps are post-peak ram pressure (see Table 2).

The NUV and Hi morphologies
The observed Hi emission stems from warm (∼ 8000 K) gas with densities > ∼ 1 cm −3 . The NUV is emitted by the photospheres of massive hot stars. It traces the star formation on timescales < 300 Myr.
The 48 simulation snapshots of Table 3 yield the highest total goodnesses for each kind of simulation (with or without diffuse gas stripping, with a variable or constant Ly scale parameter). Models with total goodnesses which ex-ceed g = 32000 for a variable Ly scale parameter (4new 47, 1cnew 47, 4b 67, and 2b 80) and g = 31650 for a constant Ly scale parameter (2b 80) are acceptable fits to the spectrophotometric observations. However, the observed Hi, NUV, and Hα morphologies have also to be reproduced by the model snapshots. To assess the quality of the fits, we compared the model and observed morphologies of the downturn and tail regions by eye. The results are presented in Table 3. The key features used for the quality assessment are: Hi: truncated disk, extraplanar gas south of the galactic disk, southwestern tail; NUV: disk emission extending further than the Hi emission, southwestern tail offset from the Hi tail; Hα: truncated disk, southwestern tail offset from the NUV tail, low surface brightness north-south filaments. The adopted signs are: +: reproduced; ∼: approximately reproduced; -: not reproduced.
The model NUV and Hi distributions of the six model snapshots with highest goodnesses from the simulations with a diffuse ISM component and a variable Ly scale scale parameter are presented in Fig. 4. They can be directly compared to the upper left panel of Fig. 1.
The observed truncation radius of the Hi disk of 7.5 kpc as well as its symmetry along the major axis are well reproduced by all model snapshots. The observed strong asymmetry of the Hi distribution along the minor axis with a significant amount of gas pushed by ram pressure downwind of the major axis is also reproduced by all model snapshots. The observed southwestern tail is present in all model snapshots. However, the gas tails of models 4new 47, 1cnew 47, 3dnew 47, and 4anew 47 are somewhat more extended than the observed gas tail. Moreover, the gas distribution of the eastern downturn region of all model snapshots are more extended in the downwind direction than it is observed. The observed gas distribution of the downturn region is best reproduced by model 4new 47.
In the downturn region the observed NUV distribution is more extended along the disk than the Hi distribution. This feature is reproduced by all model snapshots. However, the extent of all six model NUV disks are smaller than the observed one. In addition, the NUV emission of all model snapshots is extended in the downwind direction. Whereas the model NUV downturns filaments are perpendicular to the galactic disk in models 1cnew 47, 4anew 47, 2anew 64, and 2anew 65, they are further bent toward the galactic disk in models 4new 47 and 3dnew 47. In this respect, the latter two model snapshots are most similar to observations.
We observe an interesting feature in the NUV distribution of models 2anew 64, 2anew 65: an extension in the upwind direction on the northwestern side of the galactic disk. This extent is caused by star created in the UV tail which have the time to rotate into the downturn region. Since their angular momentum is that of the tail gas, these stars rotate in a plane which is tilted with respect to the galactic plane leading to a position angle opposite to that of the tail in the downturn region. Because these stars need to have the time to fulfill about half of a galactic rotation, this effect occurs in simulation snapshots more than several 10 Myr after peak ram pressure. These model snapshots can thus be excluded. The observed bending of the southwestern NUV tail is reproduced by all model snapshots. Its observed offset with respect to the Hi tail is only present in model 1cnew 47. We conclude that the observed Hi and NUV morphologies are best reproduced by model 4new 47 followed by model 1cnew 47.
The time to peak ram pressure is ((N − 1) × 10) − t peak Myr, where N is the number of the timestep (last number of the model names). The quality assessment of the reproduction of the observational characteristics is based on the joint fit of the SEDs and optical spectra.   Fig. 3.
To demonstrate the effect of diffuse gas stripping on the model Hi and NUV distribution, we show the same model snapshots as in Fig. 4 for simulation without diffuse gas stripping in Fig. 5. The resulting Hi and NUV distributions of all models snapshots are more extended than those of the snapshots of the corresponding models with diffuse best-fit gas stripping and observations. Indeed, the snapshots with the highest goodnesses for the models without diffuse gas stripping are not the same as those presented in Fig. 4 and 5 (see Table 3).
The snapshots with the highest goodnesses for the models without diffuse gas are presented in Fig. 6. All model with high snapshot numbers (4b 67, 2b 80, 4b 68, and 2d 70) show extensions in the upwind direction on the northwestern side of the galactic disk and therefore are excluded. The Hi and NUV tails of the two other models, 3c 47 and 3 47, are less bent than the "best-fit" models with diffuse gas stripping (4new 47 and 1cnew 47). In addition, extended Hi tails in the downturn region of 3c 47 and 3 47 are not observed. We thus conclude that the models with diffuse gas stripping reproduce the observed Hi and NUV morphologies significantly better than the models without diffuse gas stripping.

The Hα morphology
The model Hα images consist of two components: (i) the Hii regions ionized by young massive stars and (ii) diffuse ionized gas that is ionized by the stellar UV radiation or possibly by strong shocks induced by ram pressure stripping. The first component is modeled by the distribution of stellar particles with ages less than 20 Myr. This is the only component used for the images of the models without diffuse gas stripping. For the second component the emission measure for all gas particles with temperatures smaller than 10 6 K was calculated: EM = ρ 2 d, where ρ is the gas density calculated with SPH algorithm and the cloud diameter is given by d = 2 (3 M cl /(4 πρ)) 1/3 with the cloud mass M cl . In the absence of a model for the ISM ionization fraction we assume a constant ionization for the diffuse and dense gas. We then added the Hii region component with a heuristic normalization which worked best to insure that extraplanar Hii region are well visible within the emission of the diffuse ionized component. In these images the surface brightness of the disk emission with respect to the extraplanar emission is realistic in a qualitative but not in a quantitative sense. Fossati et al. (2018) observed an Hα tail whose high surface brightness part is clumpy and follows the NUV tail (upper right panel of Fig. 1). A diffuse component of low surface brightness extends downwind to the south. In addition, low surface brightness filaments extend further from the tail to the south. Another low surface brightness Hα filament is detected close to the downturn in the northeast. As noted by Fossati et al. (2018) all filaments indicate the direction of the motion of the galaxy.
The resulting model Hα emission maps for the snapshots of highest goodness of the diffuse gas stripping model are presented in Fig. 7. All model snapshots show an Hα tail with the same extent as the NUV tail as it is observed. Most of the models show extraplanar linear filaments in the downwind region of the galactic disk. These filaments are thicker than the observed Hα filaments (see Sect. 3.3). They are labeled for model 4new 47 in the upper left panel of Fig. 7. Their lengths are comparable to those of the observed Hα filaments. All models show an extraplanar linear filament originating from the downturn region which extends out of the disk toward the south. This filament has a north-south direction in models 4new 47, 1cnew 47, and 3dnew 47. It has a position angle of ∼ 135 • , i.e. it extends toward the southeast, in the other models. A second linear filament originating close to the galaxy center and extending into the downwind region is observed in models 4new 47, 1cnew 47, and 4anew 47. A third northsouth filament is observed in the middle of the tail of model 4new 47. A fourth north-south filament starting from the outer edge of the Hα/NUV tail is present in models 4new 47, 1cnew 47, and 4anew 47. The direct comparison between model 4new 47 and the VESTIGE observations (Fig. 1) shows that the latter filament might be interpreted as an extension of the gas tail. Because of the vertical orientation of its outermost part (compare to model 1cnew 47), we rather interpret this feature as a filament. Fossati et al. (2018) identified five extraplanar filaments in the downwind region: one originating close to the downturn region and four filaments in the tail region. Model 4new 47 with its four Hα filaments resembles most the observed Hα morphology. We therefore conclude that model 4new 47 reproduces best the deep VESTIGE Hα observations (lower right panel of Fig. 1).
To demonstrate the effect of diffuse gas stripping on the model Hα emission distribution, we show the same model snapshots as in Fig. 7 for simulation without diffuse gas stripping in Fig. 8. All model snapshots show an Hα tail comparable to that of the simulations with diffuse gas stripping. As expected, the number of extraplanar linear filaments is significantly decreased compared to the simulations with diffuse gas stripping. Only models 4 47, 1c 47, and 3d 47 show such a filament starting from the downturn region. However, its position angle differs significantly from zero (north-south direction). Vertical filaments along the projected direction of ram pressure stripping are only present in the diffuse gas because of its increased stripping efficiency. The ratio between tidal and ram-pressure acceleration is higher for the dense gas which is thus less accelerated along the wind direction. Based on the absence of north-south filaments we conclude that the presence of diffuse gas stripping increases the number of extraplanar Hα filaments in the downwind region and leads to a northsouth alignment of these filaments. The simulations with diffuse gas stripping thus reproduce observations best.
The snapshots with the highest goodnesses for the models without diffuse gas are presented in Fig. 9. The Hα tail is comparable to that of the previously discussed models. The extraplanar filaments starting from the downturn region are still present in all model snapshots. However, their position angle is not zero as it is observed. Only model 2d 70 shows north-south filaments starting from the galactic disk and the tail region. Only the filament originating at the end of the tail deviates from the north-south alignment of the other filaments, in contrast to the corresponding filament 4 in model 4new 47 (upper left panel of Fig. 7). We think that this difference is significant.
We conclude that the existence of north-south Hα filaments is a natural consequence of the presence of diffuse gas stripping. However, in rare cases, such filaments can also be created in models without diffuse gas stripping. Models 4new 47 and 2d 70 reproduce best the deep VESTIGE Hα observations.

The "best-fit" model
The results of the comparison between the model and observed Hi, NUV, and deep Hα emission distributions are Fig. 6: NUV (grayscale), Hi (solid contours), and stellar (dotted contours) images. Models without a diffuse ISM component. For the spectral fitting a variable nuisance parameter Ly scale was used. presented in Table 3. Model 4new 47 reproduces the Hi and Hα emission distribution in the downturn and tail region in a satisfactory way. Moreover, the overall morphology of the observed NUV emission is qualitatively reproduced by the model. However, the NUV disk in the tail region is less extended than observed. This is mainly an effect of NUV surface brightness, which is lower in the model NUV tail compared to observations (see also the UV photometry in Figs. 10 to 12). As a result, the offset between the Hα and NUV tails is much less pronounced than it is observed. Moreover, the model NUV emission of the downturn region shows too much extraplanar emission on the downwind side compared to the GALEX observations. The observed offset between the NUV and the Hi tail is absent in this model. This offset is only present in model 1cnew 47 and the observed NUV emission distribution is best reproduced by this model snapshot. However, model 1cnew 47 shows too much extraplanar Hi emission in the downturn region. Based on the comparison of all observed characteristics with our model snapshots ( Table 3) we conclude that model 4new 47 fits the available observations including the optical spectra best.

The optical spectra
The model optical spectra and the model spectral energy distribution (SED) for the "best-fit" model 4new 47 (including diffuse gas stripping and a variable nuisance parameter Ly scale ) are presented in Figs. 10 to 15 together with the star formation history for the regions labeled in Fig. 3. The optical to NIR SED only depends on the SFH prior to the onset of the simulations and is always well reproduced by the models. The FUV and NUV flux densities are well reproduced by the model in regions 2, 3, and 11. They are somewhat underestimated by the model in downturn regions 12 and 13 and overestimated in the outer tail region 1. The optical spectra of all positions are well reproduced by the model. For a seamless comparison we added the results of Fossati et al. (2018) for each position to Figs. 10 to 15 3 . It becomes clear that the good resemblance between the model and observations for all regions is caused by the good resemblance between the observed and modeled star formation histories. Whereas the SFHs of Fossati et al. (2018) are analytical with exponentially declining SFR after gas stripping, the model SFHs show more variations. If the model SFHs show a late peak of the SFR after the general decline of the SFR (region 2 and 12), the onset of the decline occurs earlier than that derived with a simple analytical function as in Fossati et al. (2018) because this additional young stellar population has to be compensated by a removal of stellar populations created just before the onset of SFR quenching. We conclude that our "best-fit" model 4new 47 reproduces the observed optical spectra and SED in a satisfactory way.

A constant nuisance parameter Ly scale
For the calculation of the model goodness the nuisance parameter Ly scale was varied systematically between 1/3 and 3 to search for the value which lead to the highest goodness. The parameter Ly scale takes into account that either some UV flux from the massive stars located in a given aperture escapes the region without ionizing the gas or that massive star outside the aperture ionize gas located inside the aperture. In addition, ram pressure induced shock can in principle lead to the ionization of the ambient gas. In the first case Ly scale is smaller than one, in the two latter cases it exceeds unity. The resulting Ly scale parameters for Fossati et al. (2018) model SFH and the SFHs derived from the dynamical models with the highest goodnesses are presented in Table 5. The analytical SFHs of Fossati et Its values range from 0.33 to 0.75. The Ly scale parameter found for the dynamical models can exceed unity ranging from 0.33 to 1.83. Especially our "best-fit" model 4new 47 leads to Ly scale > 1 in five out of six regions. Since these values are not exceedingly high, we think that this model is well acceptable. Values exceeding unity could imply that ram pressure induced shock probably play a role in ionizing the ambient ISM. However, the high flexibility of the fitting model and the PSF mismatch between the spectrum and the photometry prevent such a conclusion. All spectra and SEDs were recalculated using the values of the Ly scale parameter found by Fossati et al. (2018). The resulting model snapshots with the highest goodnesses are presented in Table 3 for the simulations with and without a diffuse gas stripping. As stated previously, all goodnesses of the models without diffuse gas stripping are smaller than those of the model with diffuse gas stripping. In addition, a fixed Ly scale parameter for each region lead to highest goodnesses for snapshots at later times than those for the models with a variable Ly scale parameter.
The resulting model Hi, NUV, and Hα images are presented in Fig. C.1 to C.4. Whereas the model Hi distributions are similar to the models with a variable Ly scale parameter, the snapshots from the models including a fixed Ly scale parameter all show extraplanar NUV emission in the upturn region on the upwind side (Figs. C.1 and C.3). This NUV emission, which is not observed, is due to massive stars created in the tail region which had the time to rotate in to the upturn region (see Sect. 5.1). In addition, none of the models with a fixed Ly scale parameter shows the vertical low surface brightness Hα filaments on the downwind side of the NUV/Hα tail. We therefore conclude that the models with a fixed Ly scale parameter should be discarded and this parameter has primarily a nuisance role in the fit rather than a physical meaning.

Discussion
The inclusion of diffuse gas stripping significantly changed the results of our simulations. Since a variable Ly scale parameter leads to the highest goodnesses, we do not consider the models with a constant Ly scale parameter in the following analysis. When the tail region is considered separately, the timesteps of interest are ∼ 50 Myr earlier before peak ram pressure when diffuse gas stripping is included in the simulations (Table 4). This is expected, because the gas which is pushed into the tail becomes diffuse and is stripped with a higher efficiency than the denser gas. Unexpectedly, the opposite trend is found for the downturn region: the timesteps of interest are ∼ 70 Myr earlier before peak ram pressure when diffuse gas stripping is not included. The reason for this behaviour is the following: as expected, the diffuse gas is stripped rapidly out of the galactic disk on the windward side when diffuse gas stripping is included. At the same time, the gas which remains within the disk is compressed to higher densities compared to the simulations without diffuse gas stripping. This results into a more compact gas distribution at t > ∼ 200 Myr after the onset of the simulations, i.e. the surface density contrast of the gas located within the galactic disk is enhanced. Since more time is needed to push the outer dense arms to smaller galactic radii, the timesteps of interest for the downturn (windward) side are delayed when diffuse gas stripping is included. The bulk of the star formation occurs in the dense gas arms.
The situation changes when both, the downturn and tail regions, are fitted together (Table 3). Generally, the Hα and NUV emission distributions of the downturn region are more often reproduced when diffuse gas stripping is not included in the simulation. Overall, the Hi distribution of the downturn region is hard to reproduce (4new 47 and 4 d70). The best reproductions of the observed Hi distribution of the tail are only found in models with diffuse gas stripping. The Hα emission distribution of the tail region is as often reproduced when diffuse gas stripping is not included, but only for peak and post-peak models. Most frequently, postpeak models are not allowed, because they show a NUV horn feature in the leeward side of the downturn region (see Sect. 5.1).
Whereas the two acceptable models with diffuse gas stripping according to the total goodness (4new 47 and 1cnew 47) have ∆t peak = t − t peak ∼ −170 Myr, i.e. they are pre-peak models, the two acceptable models without diffuse gas stripping (4b 67 and 2b 80) are post-peak models (∆t peak > 0 Myr). Based on the comparison between the model and observed Hi, NUV, and Hα distributions the pre-peak models 4new 47 and 1cnew 47 are clearly preferred (Sect. 5.3).

The "best-fit" model
The "best-fit" model is 4new 47 with a maximum ram pressure of p max = 6000 cm −3 (km s −1 ) 2 , a width of the Lorentzian temporal profile of t HW = 200 Myr, and a time to peak ram pressure of ∆t = t − t peak = −140 Myr. The face-on projection of the distribution of the total gas surface density at the beginning of the simulation and at ∆t = −140 Myr is presented in Fig. 16. The initial gas distribution of the galaxy has multiple spiral arms with surface densities between 10 and 20 M pc −2 . During the ram pressure stripping event we observe a highly non-linear interaction between these arms and the induced gas compression. The gas distribution at ∆t = −140 Myr has a compressed windward side and a gas tail with a multi-arm structure on the opposite side.
The model includes the stripping of diffuse gas and a variable Ly scale parameter for the calculation of the model optical spectra and thus the model goodness. The interaction parameters are close to those found by : p max = 5000 cm −3 (km s −1 ) 2 , t HW = 100 Myr, and ∆t = −100 Myr. The  was not able to reproduce the UV color in the stripped gas-free regions, because the width of the Lorentzian profile was too small. A larger width permits an earlier gas stripping which is consistent with the available UV observations and the optical spectra. In addition, the inclusion of diffuse gas stripping is essential. Only if this ingredient is added to the simulations the gas leaves the galactic disk early enough to lead to FUV/NUV flux densities and optical spectra consistent with observations. The diffuse gas stripping leads to the formation of linear filaments of diffuse ionized gas which is present in deep VESTIGE Hα observations. The model filaments are thicker than the observed filaments because of the coarse spatial resolution of the diffuse gas particles. It should be noted that pronounced filamentary structures are naturally produced in simulations whenever magnetic fields are taken into account (Ruszkowski et al. 2014;Tonnesen & Stone 2014). It has to be shown for other galaxies undergoing ram pressure stripping if this effect is ubiquitous.
6.2. Ram pressure stripping of cosmic ray electrons NGC 4330 does not show an asymmetric ridge of polarized emission as other ram pressure stripped galaxies ). This is due to the NGC 4330's projection which places the region of compressed gas along the line of sight. On the other hand, Vollmer et al. (2012) discovered a tail structure bending southward in their 6 cm radio continuum data (upper panel of Fig. 17). This tail structure is also visible in the VIVA 20 cm data (lower panel of Fig. 17). In addition, an outer tail of low surface brightness extending toward the southwest is detected at 20 cm. The tail, which is detected at 6 and 20 cm, bends to the south at the southwestern edge of the starforming disk. It coincides with the northern part of the low surface brightness diffuse Hα emission. The cosmic ray electrons produced in the luminous Hii regions of the galactic disk are most probably pushed to the south by ram pressure. This hypothesis is corroborated by the vertical magnetic field lines detected within the radio continuum tail (Vollmer et al. 2013). The spatial coincidence of the radio continuum and diffuse Hα tails suggests that both gas phases are stripped together.
We do not observe a steepening of the spectral index in the 6 and 20 cm tail as expected for an aging cosmic electron population. The spectral index is given by α = ln(f 6 /f 20 )/ln(ν 6 /ν 20 ), where f is the surface brightness or flux density. For a purely aging population of cosmic ray electrons the evolution of the flux density is given by f = f 0 ×exp(−t/t syn ). The time which is needed for a steepening of the spectral index from an initial value α in to a final value α end can be calculated in the following way: Based on equipartition between the energy densities of the cosmic ray electrons and the magnetic field we estimate the magnetic field strength to be B ∼ 6-7 µG in the tail region assuming an extent along the line-of-sight of l ∼ 4 kpc. The associated synchrotron timescale at 6 cm is then t syn = 4.5 × 10 7 ( B 10 µG ) −3/2 ν −1/2 GHz yr ∼ 40 Myr .
If we assume that α in − α end < ∼ 0.2 in the 6 and 20 cm tail the timescale for the steepening of the spectral index is t SI < ∼ 24 Myr. With a north-south extent of 3.5 kpc the southward bulk velocity is > ∼ 140 km s −1 within the sky plane.
Assuming that the outer radio continuum tail is barely detected at a 2σ level at 6 cm, we estimate a spectral index between 6 and 20 cm of ∼ −1.5. With an initial spectral index of α in = −0.7 to −0.5, a final spectral index α in = −1.5, and a magnetic field strength of B = 6 µG we obtain t SI ∼ 90 to 120 Myr. If the cosmic ray electrons were produced in the normally starforming disk before gas stripping t SI should be comparable to the star formation quenching time derived for that region which corresponds to region 2 in Fossati et al. (2018). Whereas Fossati et al. (2018) determined a quenching time of 175 Myr based on their simple analytical star formation histories, our preferred model shows that a more recent star formation episode (about 50 Myr ago; Fig. 11) is also consistent with the data. It is not clear if such an episode can create enough cosmic ray electrons to make the synchrotron emission detectable. On the other hand, a somewhat lower magnetic field strength of B = 4-5 µG would lead to t SI ∼ 175 Myr. We thus conclude that a scenario in which the cosmic ray electrons were produced in the healthy disk before gas stripping and were then pushed southward together with the diffuse ionized gas is consistent with the available data.
It is surprising that there is no 20 cm emission associated with the upper part of the outer Hα tail where we see compact emission, most probably Hii regions. Either these Hii regions are too young and did not produce cosmic ray electrons yet or the production of cosmic ray electrons by these regions is low. We suggest that the present star formation in the outer tail is sporadic and low level and this explains the absence of a significant amount of cosmic ray electrons there.

How to ionize the stripped diffuse gas
It is surprising that the surface brightness of the Hα filaments does not systematically decrease southward. These filaments therefore need to be confined by the ambient ICM during ∼ 60 Myr without recombining, i.e. there must be a continuous ionization mechanism. In the following we argue that the stripped diffuse gas is ionized by thermal conduction. With a recombination timescale t rec = 0.1 (n e /1 cm −3 ) −1 Myr=60 Myr this yields an electron density of n e ∼ 2 × 10 −3 cm −3 . The mean electron density of the Hα filaments can be estimated by the relation L(Hα) = n e n p α eff (Osterbrock & Ferland 2006) where n p = n e is the proton density, α eff Hα = 1.17 × 10 −13 cm 3 s −1 is the Hα effective recombination coefficient, V is the volume of the emitting region, φ V its volume filling factor, h the Planck constant, and ν Hα the frequency of the Hα transition. The volume V is the product of the emitting surface and the extent along the line-of-sight l ∼ 4 kpc. Following Boselli et al. (2016) we assume φ V = 0.1 and [Nii]/Hα=0.5. With a mean Hα+ [Nii] surface brightness of 4×10 −18 erg cm −2 s −1 arcsec −2 ) we obtain a mean electron density of n e = 0.05 cm −3 . With an emitting surface of 25 kpc 2 the mass of ionized gas in the filaments is about 10 7 M .
Since the derived mean electron density of the filaments exceeds that derived from the recombination timescale, the hydrogen atoms in the filaments have to be continuously ionized. Boselli et al. (2016) reached the same conclusion for the long, low surface brightness Hα tails of NGC 4569, another ram-pressure stripped galaxy in the Virgo cluster. The source of ionization might be (i) collisional ionization through shocks, (ii) collisional ionization through the thermal electrons of the ambient ICM which confines the filaments, or (iii) radiative ionization by embedded massive stars. We discard the latter hypothesis because the gas density in the filaments is not high enough to permit continuous star formation. Hypothesis (i) is improbable because the direction of filaments is aligned with the ram pressure wind and it is difficult to perceive how ram pressure induced shocks can pervade the entire filaments.
We are thus left with hypothesis (ii): if we assume that the timescale for evaporation by the ICM equals the collisional ionization timescale by ICM-ISM collisions t ICM−ISM = t evap = 10 (N H /10 20 cm −2 ) Myr (Vollmer et al. 2001) and if the recombination timescale is t rec = 0.1 (0.3 cm −3 /1 cm −3 ) −1 Myr=0.3 Myr, we obtain a hydrogen column density of N H ∼ 3 × 10 18 cm −2 . This is well below the VIVA detection limit (∼ 2 × 10 19 cm −2 ; Abramson et al. 2011). We thus suggest that the most probable mechanism for the ionization of the Hα filaments are collisions between the ICM thermal electrons and the hydrogen atoms of the filaments (as in ionized gas filaments in cooling flow clusters; see e.

Summary and conclusions
The edge-on Virgo spiral galaxy NGC 4330 shows signs of ongoing ram pressure stripping: a truncated Hi disk together with a one-sided tail (Chung et al. 2007), a UV tail associated with the Hi tail (Abramson et al. 2011), a truncated Hα disk (Abramson et al. 2011), and an Hα tail with vertical low surface density filaments ; upper panels of Fig. 1). There is no asymmetric ridge of polarized radio continuum emission as observed in other ram-pressure stripped Virgo galaxies because the compression is occurring in the sky plane . A previous dynamical model could reproduce all observed characteristics except the UV color of the stripped gas-free regions of the outer galactic disk . Since the UV color depends on the star formation history of the observed region, the stripping timescale of the dynamical model was not consistent with observations. Fossati et al. (2018) obtained VLT FORS2 deep optical spectra along the major axis of NGC 4330. In addition, they collected literature photometry in 15 bands from the far-UV to the far-IR. Star formation histories in apertures along the major axis of the galaxy were constructed by means of stellar population fitting to the joint observational data, SED and optical spectrum. Based on an analytical model of quenched star formation histories, they found that the outermost radii of the galactic disk were stripped ∼ 500 Myr ago, while the stripping reached the inner 5 kpc in the last 100 Myr.            (1,2,3,5,7,10,15,20,30,50, 100) × 50 µJy/beam. The resolution is (22 × 22 ). Lower panel: emission at 20cm from VIVA ). The contour levels are (1,2,3,5,7,10,15,20,30,50, 100) × 75 µJy/beam. The resolution is (19.4 × 15.4 ).