Galactic seismology: joint evolution of impact-triggered stellar and gaseous disc corrugations

The initial conditions (particle positions and velocities) of each of the components making up our Galaxy model are created with the Action-based GAlaxy Modelling Architecture software package (agama; Vasiliev, 2019). We refer the reader to Tepper-Garcia et al. (2021, their section 3) for a detailed description of agama’s self-consistent modelling module for the collisionless components in our model (DM halo, bulge, stellar disc). We have complemented the agama galaxy model framework to include the gas phase (cf. Deg et al., 2019); our methodology is described elsewhere (Tepper-García et al. 2022b, in prep.). Here we provide only a brief description.

Our approach follows Wang et al. (2010) who give a prescription for setting up isothermal gas discs in equilibrium. In brief, the gas disc is isothermal and initially axisymmetric, with a surface density profile described by a radially declining exponential, $\propto\exp[-R/R_{g}]$, with a scalelength $R_{g}=7$ kpc. Its vertical structure follows a $\operatorname{sech}^{2}[z/z_{0}]$ profile, appropriate for a gas distribution in vertical hydrostatic equilibrium, with a scaleheight²²2We estimate the scaleheight by measuring the root-mean-square of the height of the gas with respect to the galactic plane as a function of $R$ $z_{0}$ that varies with cylindrical radius from $z_{0}\approx 20$ pc near the centre to $z_{0}\approx 160$ pc at $R=20$ kpc (a ‘flaring’ disc). The radially increasing thickness of the gas disc is a consequence of the isothermal constraint and the fact that the (vertical) potential weakens with radius from the galaxy’s centre. The azimuthal velocity profile of the gas disc ensures rotational support against radial instabilities (see Wang et al., 2010, their equation 13).

We track the evolution of the following models: 1) the isolated hybrid model; and 2) the interaction hybrid model. The former is intended as a reference simulation to gauge the response of the synthetic galaxy to the interaction, and to assess the impact of any artefacts (e.g. numerical noise). In the interaction hybrid simulation, the point mass representing Sgr is placed at $\vec{r}\approx(-11.1,0,28.6)$ kpc with a velocity $\vec{v}=(-145.4,0,-220.2)$ km s^-1 with respect to the galactic centre.

For both simulations, the synthetic galaxy and the infalling point mass are evolved with the ramses code (Teyssier, 2002) which incorporates adaptive mesh refinement (AMR). The system is placed into a cubic box spanning 600 kpc on a side. In addition to the galaxy and the perturber, the simulation box is initially filled with a cold ($T=10^{3}$ K), tenuous, uniform gas. We do not attempt to set the ambient background gas in hydrostatic equilibrium, which is a non-trivial but feasible task (e.g. Grønnow et al., 2021). To deter the ambient gas from collapsing in a timescale shorter than our simulation time span, we set its density to an extreme threshold value ($10^{-20}$ cm^-3).

We stress that our simulations are strictly isothermal, i.e. the gas is artificially kept at a temperature of $10^{3}$ K throughout. But to ensure this equation of state does not influence our dynamical results, we have run one additional (expensive) simulation to $t\approx 1$ Gyr with a different equation of state for the gas. This test simulation starts from the same initial conditions and using a setup identical to those of our isolated hybrid model, with the exception of the following: 1) The ambient gas is initially hot ($T=10^{6}$ K); 2) we include an idealised, uniform ultraviolet (UV) background (appropriate for a redshift $z=0$) that heats the gas; 3) we allow the gas to cool radiatively, assuming it is composed of hydrogen, helium, and a mixture of metals resulting in a metallicity equivalent to 10 percent of the solar value. To deter the gas from cooling indefinitely, we set a temperature floor at $T=10^{3}$ K. In addition to the effect of UV photo-heating, the gas temperature increases as a result of compression.

We have compared the evolution of the discs with cooling/heating and the strict isothermal condition, and while there are differences between them, none are important for our dynamical study. This may seem contradictory to our assertion that the setup of the gas disc must be strictly isothermal. However, due to the existence of a temperature floor, the high gas density across much of the gas disc as well as its relatively high metallicity, the gas disc to a large extent retains its equilibrium temperature (and thus its structure) throughout, while the ambient medium barely cools during the time span of the simulation. In other words, the simulation starts and ends up (after $\sim 2$ Gyr) with virtually two distinct gas phases: a cold, dense disc, and a hot, tenuous background, connected by an insignificant intermediate temperature phase.

Given the similarity between the simulation with cooling/heating and the isothermal condition, and also because the former is significantly more expensive to run, we choose to investigate the corrugations in stars and gas with our isothermal simulations.

In addition to neglecting any cooling/heating processes, we ignore the presence of a hot halo around the Galaxy (e.g. Miller & Bregman, 2013), or the presence of magnetic fields (e.g. Jansson & Farrar, 2012). We also neglect galactic processes such as star formation or stellar feedback of any kind. Our aim is to construct the simplest possible model to isolate the dynamical effect of the interaction from other, potentially obfuscating processes. More realistic simulations incorporating the physics we have neglected will be explored in future studies.

The total simulation time is about 2 Gyr for both the isolated hybrid model and the interaction hybrid model. At runtime, the AMR grid is maximally refined up to level 14, implying a limiting spatial resolution of 600 kpc / $2^{14}\approx 36$ pc. As a result, the vertical structure of the gas disc is barely resolved at $R\lesssim 4$ kpc.

A requirement of numerical experiments aimed at investigating the response of a system to a prescribed stimulus is for the system to be in near-perfect equilibrium when left in isolation, and for it to maintain that state for a time span that is long compared to the actual experiment. We have shown in our previous work that the isolated N-body model preserves its initial state for at least 4 Gyr (BT21). Since the isolated hybrid model follows virtually an identical setup, we expect this model to retain the same degree of stability. To this end, we analyse the overall stability of isolated hybrid model by comparing the initial structural and kinematic properties of the stellar disc and of the gas disc with its state after 2 billion years of evolution. A detailed discussion of this analysis is deferred to the Appendix. Here it suffices to say that the synthetic galaxy, and in particular both the stellar and gas discs, maintain their initial state throughout their evolution in isolation, implying that our setup yields a stable model from the outset.

To summarise: In what follows, we discuss three different simulations: 1) isolated hybrid simulation; 2) interaction hybrid simulation; and 3) interaction N-body simulation. Simulations 1 and 2 are new and are presented for the first time in the present paper. Simulation 3 was introduced and discussed at length in BT21.

While the stellar and gas discs are fundamentally different in character, we have already seen that they share a common evolution at the outset but evolve separately at later times. To understand the mechanism behind this divergent behaviour, we explore whether there exists a measurable correlation between the vertical response of the stellar and gas discs to the perturbation induced by the satellite. We consider a quantitative analysis of $z$ and $V_{z}$ simultaneously by mapping the spatial and the kinematic data to the vertical phase space $(z,V_{z})$. Consider a pair of $\overline{z}$ and $\overline{V_{z}}$ maps of the synthetic galaxy at a given epoch (Figs.  1 and  2). But instead of looking at $\overline{z}$ or $\overline{V_{z}}$ separately, we combine these two quantities into a vector, an element of $(z,V_{z})$ space. However, because $z$ and $V_{z}$ carry different units, we apply a scale transformation to bring them to a common dimension. The natural choice is to scale $\overline{z}$ by the vertical frequency $\nu$, defined by $\nu^{2}(R)=d^{2}\Phi/dz^{2}\,,$ where $\Phi(R,z)$ is the total potential at $(R,z)$, and the derivative is evaluated at $z=0$.

Thus we define the normalised phase space coordinates

and the corresponding phase-space vectors, $w_{g}\equiv\{\tilde{z},\tilde{V}_{z}\}_{g}$ and $w_{\star}\equiv\{\tilde{z},\tilde{V}_{z}\}_{\star}\,,$ where the subscript is used to distinguish between the vectors corresponding either to the gas or to the stars (Fig. 3). These vectors are elements of the plane spanned by $(\tilde{z},\tilde{V}_{z})$, a transformed version of the vertical phase space. Note that this transformation renders the distribution of stars on this plane roughly circular, regardless of their location in cylindrical radius $R$.

Within the transformed phase space, we measure the angle $\alpha$ between the vectors $w$ geometrically via

where $|w|$ indicates the length of the vector $w$ to the centroid of the distribution.

Evaluation of Eq. (2) over a pair of maps such as the left panels of Figs. 1 and 2 yields values in the range $[-1,1]$; non-negative values correspond to angles (measured in radian) in the range $|\alpha|\leq\pi/2$, and thus indicate that stars and gas at $p(x,y)$ oscillate vertically in phase. Negative values indicate the opposite, i.e. the gas and stars are out of phase. We further categorise the in-phase corrugations in stars and gas to be strongly correlated if $|\alpha|\leq\pi/4$ (corresponding to $\cos\alpha\gtrsim 0.71$) or weakly correlated if $\pi/4\leq|\alpha|\leq\pi/2$ (corresponding to $0\leq\cos\alpha\lesssim 0.71$).

Using the approach just described, we calculate the angle $\alpha$ and the correlation at each time step ($\Delta t\approx 10$ Myr) in our simulation from $t=0$ Gyr to $t\approx 2$ Gyr, i.e. from roughly -100 Myr prior to the impact to roughly 1800 Myr after the impact. Thus we are able to quantify the evolution of the correlation angle $\alpha$ between the vertical motion of stars and gas across the galaxy. In order to quantify the occurrence of by-chance correlations, we repeat the analysis, but this time rotating each of the $\overline{z}_{\star}$ and $\overline{V_{z}}_{\star}$ maps anti-clockwise by $\pi/2$ prior to calculating the vectors $w_{\star}$. We refer to the latter as our ‘control data’, and to the former as our ‘simulation data’.

We further quantify the correlation globally (at each time step) by calculating the arithmetic average of the pixel-wise values delivered by Eq. (2),

where $\left(\cos\alpha\right)_{i}$ is the value of the angle $\alpha$ at pixel (or bin) $i$, and $M$ is the total number of bins in the map (see beginning of Sec. 3). Analogous to the angle $\alpha$, positive values of $\mathcal{C}$ indicate a positive correlation between stars and gas across the map, while negative values indicate the opposite.

The results of our approach are displayed in Fig. 4. There we present a series of panels, each corresponding to a face-on projection of the synthetic galaxy at a different epoch in the system’s evolution from its initial state roughly 100 Myr prior to impact (top left), to about 850 Myr after the impact (bottom right). The colour coding indicates the value of $\cos\alpha$: blue and green correspond to $\cos\alpha\gtrsim 0.71$ and $0\leq\cos\alpha\lesssim 0.71$, i.e. to strong and weak in-phase oscillations, respectively, while a grey colour corresponds to negative values, and thus indicate an out of phase oscillation. In each panel, the simulation time is indicated in the top right corner, and the value of $\mathcal{C}$ in the bottom right corner.

Initially, i.e. prior to the impact (top left panel), there is no correlation between the vertical response of the stars and the gas because an equilibrium distribution has no centroid shift in the transformed phase space.

At impact time ($t=95$ Myr, central panel, top row), we see a signal consistent with a strong in-phase vertical response between stars and gas at $R\gtrsim 5$ kpc. The correlation index increases up to about 150 Myr after impact; its value climbs from roughly 0.6 to 0.8 (top row, right panel). After that, the oscillations in the stars and gas start to move out of phase, as indicated by a weakening of the overall signal (bottom row), a trend that continues until the end of the simulation.

It is interesting that the difference in correlation angle between stars and gas reveals a great deal of substructure in the disc, notably the spiral structure. Presumably, the contrast between stars and gas in this space arises because the gas is partially compressed by the stellar density wave, and is thus displaced downstream, causing it to move out of phase with respect to the stars along the density wave and generating a ‘galactic shock’ (Fujimoto, 1968; Roberts, 1969; Lubow et al., 1986; Baba et al., 2015).

The global behaviour of the phase correlation is captured in the top panel of Fig. 5. The black symbols (+) trace the behaviour of the global correlation angle between the corrugations of the stellar disc and the gas disc (Eq. 3) over the full simulation. The red symbols ($\bullet$) indicate the corresponding result for our control data. The background has been colour-coded following the same convention followed in Fig. 4. The vertical dotted line flags the impact epoch. In this form, it becomes apparent that the corrugations of stars and the gas in the disc are strongly correlated ($\mathcal{C}\gtrsim 0.71$, blue shaded area) for roughly 150 Myr after the impact, and only weakly correlated after that ($0\lesssim\mathcal{C}\lesssim 0.71$, green shaded area), consistent with the results presented in Fig. 4. It is worth noting that the level of correlation reached at $t\approx 1.3$ Gyr is comparable to the level reached at the end of the isolated hybrid simulation, likely seeded by numerical noise (see Fig. 5, bottom panel), and even drops below this level at later times. This is consistent with the dynamical response of the discs triggered by the interaction dominating over numerical noise effects for a significant fraction (if not all) of the simulation time span.

Our findings suggest that, in the absence of any further, strong perturbations, estimates of the degree of correlation between the corrugations observed in stars and gas can be used to age-date the phenomenon, regardless of whether the perturbation is extrinsic (Elmegreen et al., 1995) or intrinsic (Weinberg, 1991) to the galaxy, provided the origin of the perturbation does not lie more than $\sim 1$ Gyr in the past.

Our impulsive simulation provides the ideal framework to study the vertical response of the stellar and gaseous discs to a one-time interaction with an external perturber (cf. Chakrabarti & Blitz, 2009), without the obfuscating effects of subsequent disc crossings (e.g. Chequers et al., 2018; Laporte et al., 2018b).

In order to analyse the vertical response of the stellar and gaseous discs throughout their evolution, we proceed as explained next. We stress that our approach is preferred over using Fourier methods because we find a full modal analysis is not warranted here, at least not in our initial analysis.

At each time step during the evolution of the interaction hybrid model, the absolute value of $z$ and $V_{z}$ is mapped from the Cartesian $(x,y)$ plane onto the polar $(R,\phi)$ plane, where $R$ and $\phi$ are the cylindrical radius and the azimuthal angle, respectively. We do the same for the isolated hybrid model, which we then use to correct³³3In practice, we subtract one map from the other before taking the absolute value; our results are virtually identical if we skip this step. each $(R,\phi)$ map (i.e. at each time step) for any vertical displacement or velocity that results from numerical noise. Then, for each corrected $(R,\phi)$ map, we calculate the 50 percentile (median) and 90 percentile (here considered a proxy for the maximum value) at each radial bin along the full range in azimuth. Thus we obtain an estimate of the median and maximum vertical displacement and vertical velocity at each radius $R$ across the disc, at each time $T$, yielding two $|z|$ maps and two $|V_{z}|$ maps, one for the stellar disc and one for the gas disc, on the $(T,R)$ plane. Furthermore, to assess the potential effect of the gas disc on the stellar disc, we repeat the analysis, this time for the stellar disc in the interaction N-body model. The results are displayed in Fig. 6.

The first two rows display the $|z|$ maps, the next two rows the $|V_{z}|$ maps. The first and third rows display the 50 percentile in $|z|$ and $|V_{z}|$, respectively, while the second and fourth rows the 90 percentile of each quantity. The left column displays the results for the stellar disc in interaction N-body model, while the central column and right columns corresponds to the interaction hybrid model; the central column displays the result for the stellar disc, and the last column, the result for the gas disc.

Prior to the impact (whose epoch is flagged by the white vertical line at $t\approx 100$ Myr), neither the stellar disc (in either simulation) nor the gas disc reveal any significant vertical response, consistent with the fact that up to this point the synthetic galaxy can be regarded as an isolated system. But the situation changes dramatically after the interaction: clearly, both the stellar and the gas disc undergo a strong vertical perturbation, a bending wave.

The perturbation, as revealed either by the change in $|z|$ or $|V_{z}|$, is apparently structured, displaying a striation pattern with ridges running from top to bottom that become more slanted with time. The ridges correspond to small $|z|$ or small $|V_{z}|$ amplitudes, and are thus associated with the ‘nodes’ of the bending wave, i.e. the regions across the disc where the bending wave intersects the midplane. Their inclination (with respect to, say, the $R$-axis) on the $(T,R)$ plane is directly related to the degree of windup of both the kinematic density wave (‘spiral arms’) and the bending wave; the higher the degree of wind-up, the more inclined are these ridges.⁴⁴4 A near-perfectly circular and stationary bending wave node – an example of extreme windup – corresponds to a near-perfectly horizontal ridge in the $(T,R)$ plane.

Apparently, the perturbation induced by the satellite does not cover all radii, but extends only from $R=20$ kpc into $R\approx 7-9$ kpc, leading to the formation of a lower ‘envelope’ on the $(T,R)$ plane. This is a natural consequence of the fact that the restoring force of the disc, or equivalently, its resilience to any external perturbation, increases towards the centre (cf. Laporte et al., 2018b, their figure 17). Moreover, the width of the swath along $R$ on the $(T,R)$ plane decreases with time, and does so in an in-out fashion, implying that the coherent vertical motion of the stars, and of the gas, washes away with time, in agreement with Poggio et al. (2021). We will return to this point later in Sec.3.2.3.

We now look at some important differences across panels, i.e. differences between the spatial and kinematic responses of a stellar disc in the absence or presence of a gas disc, and the difference between in the spatial and kinematic responses of a gas disc compared to a stellar disc.

We note first that the response in the vertical displacement as given by $|z|$ is offset by $\pi/2$ from the vertical velocity, i.e. maxima in the vertical displacement at any radius and any time are coincident with minima in the vertical velocity, and vice-versa. To illustrate this point, we include in the panel, at the centre of the top row, a series of red lines roughly tracing the first four minima in the vertical displacement (50 percentile) of the stellar disc in the interaction hybrid model. We insert an identical set of lines in the central panel of the third row. Clearly, these lines trace the maxima in $|V_{z}|$. The same is true at the 90 percentile level, and it is also true for the gas disc. This result is not trivial and indicates that the corrugations induced by the satellite in both the stellar and the gaseous discs behave predominantly as travelling ripple patterns, initially in phase before drifting out of phase, as discussed earlier (cf. Poggio et al., 2021, for the same effect in a pure N-body simulation). This interpretation of Fig. 6 justifies our earlier method in using a correlation angle as a way of tracking the relative evolution of the joint star-gas corrugations (Sec. 3.2.1).

We now focus on the stellar disc in both the interaction N-body simulation and the interaction hybrid simulation (left and central columns). It is reassuring that the behaviour of the stellar disc in these two simulations is very similar. This is true both for $|z|$ (rows 1 and 2) and $|V_{z}|$ (rows 3 and 4), and both at the 50 percentile level (rows 1 and 3) and at the 90 percentile level (rows 2 and 4). The response in each case are obviously not identical, but the differences appear to be minor. We note that the comparison between these two simulations is justified, as the stellar discs in them are virtually identical to one another by design (see Fig. 11).

The similarity between the vertical response of the stellar discs in these different simulations indicates that a light, cold gas disc as we have included in our hybrid model does not significantly affect the dynamical behaviour of the stellar disc. This is important because it suggests that the onset and development of the phase spiral detected by Gaia is not strongly influenced by the presence of gas. We address this issue in a future paper.

We now compare the vertical response of the stellar and gaseous discs in the interaction hybrid simulation (central and right columns). Within the first 700 Myr after the impact, the two components appear to respond in a similar way to the perturbation, both in terms of $|z|$ (rows 1 and 2) and $|V_{z}|$ (rows 3 and 4). The median vertical displacement is $|z|\approx 1$ kpc and the median vertical velocity $|V_{z}|\approx 10$ km s^-1. At the 90 percentile level, the vertical displacement reaches values of up a few kpc, while the corresponding vertical velocity may reach values of several tens of km s^-1. These ranges are in good agreement with the corresponding values observed in the Galaxy both in the stellar disc (Widrow et al., 2012) or the gas layer (Alves et al., 2020), or in the gaseous layers of external galaxies (Sánchez-Gil et al., 2015; Narayan et al., 2020).

After $t\approx 700$ Myr, however, we observe a significant difference between the ongoing response of the stellar disc and that of the gas disc. Indeed, while both discs appear to settle with time, the gas disc does so more quickly, as indicated by the fading amplitudes in $|z|$ and $|V_{z}|$ and both at the 50 percentile and 90 percentile levels. This fact is consistent with, and may actually explain, the decrease in the correlation between the corrugations in stars and gas discussed in Sec. 3.2.1. Likely, the gas waves dampen at a faster rate compared to the stellar waves because the gas is subject to compression, which transforms kinetic energy into internal energy that is radiated away in order to satisfy the isothermal constraint imposed on the gas.

We note finally the near-vertical line apparent in the maps displaying to the response of the gas disc (right column) at around $t\approx 200$ Myr. This corresponds to the external perturber along its polar orbit and which – as it crossed the disc – accreted gas and thus becomes visible in this plane. The additional structure visible very close to the centre $R\lesssim 2$ kpc is caused by the collapse of ambient gas along and around the spin axis towards the galactic plane (see Appendix A).

To sum up, we do not find any significant differences between the vertical response of the stellar disc in a pure N-body simulation and a simulation that includes a light, cold gas disc. We find that the median and maximum vertical displacement and vertical velocity both in the stellar disc and the gaseous disc are consistent with observations during a broad time span after the impact, but that they do dampen with time. The damping of the corrugations is strongest with decreasing galactocentric radius, and it is stronger for the gas waves compared to the stellar waves.

An alternative – and complementary – view to the evolution of the stellar and gas corrugations in terms of their spatial and kinematic amplitudes is based on the analysis of the vertical energy involved in the process. We now pursue this complementary approach.

First, note that squared-length of the vector $w$ defined earlier (Sec. 3.2.1) corresponds to the total energy of a simple harmonic oscillator with frequency $\nu$ and motion constrained to the $z$-axis,

This expression is approximately equal to the vertical specific energy of an ensemble of stars. Based on the results discussed in Sec. 3.2.2, we assume that it can also be applied to the gas phase. We stress that when applying the above equation we are implicitly considering only the bulk motions of temporarily co-spatial ensembles of stars or patches of gas found in a small volume. We are also ignoring internal random motions or any momentary changes due to compression or expansion processes in the gas.⁵⁵5 Recall that the gas is kept isothermal and thus heating / cooling processes are perfectly balanced, by design. The distinction between the energy involved in the corrugation (i.e. the bulk motions), and the actual orbital energy (of stars) is important, as will become clear later.

Using the same approach used to construct the $\overline{z}$ and $\overline{V_{z}}$ maps earlier (Sec. 3.1), we construct separate $\overline{\epsilon_{z}}$ maps for the stars in the disc and for the gas disc. Figs. 7 and 8 display examples of such maps at selected snapshots throughout the evolution of synthetic galaxy in the interaction hybrid model. At $t=0$ Myr (top left panel in either figure), the vertical energy is very low both across the stellar disc and the gas disc. In the absence of an external perturbation, the energy remains virtually constant in the stellar disc, and reasonably constant in the gas disc, in particular in the outer disc (cf. Fig. 15 and 16). The behaviour of the stellar disc in the interaction N-body simulation is qualitatively similar to that of the stellar disc in the interaction hybrid simulation and its corresponding $\overline{\epsilon_{z}}$ maps are therefore omitted here.

But the interaction with a massive perturber, on the other hand, leads to a dramatic increase in the energy across both discs. The stars and the gas near the impact site at $x\approx-18$ kpc experience a significant energy transfer in the vertical direction as a result of the direct collision with the perturber. The signature in the energy increase is consistent with the onset of an $m=1$ bending mode. The stars and the gas on the opposite side also display a significant increase in their energy, as a result of the early perturbation induced by the satellite as it its on its way towards the galactic plane. Overall, the behaviour in the increase of the vertical energy across the discs mimics the behaviour we observe in the vertical displacement and in the vertical velocity (cf. Figs. 1 and 2), which is not surprising – but consistent, given the relation between $\epsilon_{z}$, $z$, and $V_{z}$ (Eq. 4).

In addition to the average vertical energy at each location on the stellar and gaseous discs, we estimate the evolution of the average vertical energy by measuring the median value of the vertical energy across the entire disc, separately for the stellar disc and for the gas disc, at each time step over the full time span of the simulation. This is useful to quantify the global, long-term behaviour of the stellar and gas corrugations. The result is presented in Fig. 9 (top panel). The blue continuous curve indicates the evolution of the energy of the stellar disc in the interaction hybrid model, while the orange continuous curve indicates the corresponding result for the gas disc. For reference, we include the evolution of vertical energy of each disc in the isolated hybrid simulation (dashed curves). In this simulation, the vertical energy of the stellar disc remains virtually constant. The vertical energy of the gas disc is initially practically zero (consistent with the fact that the disc is initially in vertical hydrostatic equilibrium), and increases slowly with time - as a result of the collapse of gas onto the disc plane (see Appendix A), but remains very low overall.

Comparing the blue continuous curve to the blue dashed curve, and the continuous orange curve to the orange dashed curve, we see that there is a clear excess in vertical energy in both discs as a result of the impact (flagged by the vertical dotted line) with respect to a no impact scenario, in line with Edelsohn & Elmegreen (1997)’s findings.

To asses the potential effect of the gaseous disc on the stellar disc, we include in the top panel of Fig. 9 the evolution of the vertical energy in the stellar disc for the interaction N-body model (black continuous curve), and the corresponding result for the isolated N-body simulation (black dashed curve).

We note that in all cases the vertical energy decreases with time, but at different rates, as suggested by the difference between the continuous curves towards the end of the simulation. By fitting a simple exponential $\propto\exp{[-t/\tau]}$ to the evolution of the vertical specific energy of each of the discs from $t=600$ Myr onward⁶⁶6This choice is motivated by the fact that all three curves intersect roughly at $t=600$ Myr. we learn that the discs each lose vertical energy on a typical time scale $\tau\approx 1350\pm 17$ Myr (stellar disc, interaction N-body model), $\tau\approx 1051\pm 14$ Myr (stellar disc, interaction hybrid model), and $\tau\approx 802\pm 15$ (gas disc). These results are summarised in Tab. 2. Thus the corrugations in gas disc disc dampen at a faster pace, compared to the stellar disc. Specifically, between $t=600$ Myr and $t=1800$ Myr, the energy of the corrugations in the interaction hybrid simulation have declined by roughly 68 percent (stellar disc) and roughly 78 percent (gas disc).