A toy model for gas sloshing in galaxy clusters

We apply a toy model based on 'pendulum waves' to gas sloshing in galaxy clusters. Starting with a galaxy cluster potential filled with a hydrostatic intra-cluster medium (ICM), we perturb all ICM by an initial small, unidirectional velocity, i.e., an instantaneous kick. Consequently, each parcel of ICM will oscillate due to buoyancy with its local Brunt-V\"ais\"al\"a (BV) period, which we show to be approximately proportional to the cluster radius. The oscillation of gas parcels at different radii with different periods leads to a characteristic, outwards-moving coherent pattern of local compressions and rarefactions; the former form the sloshing cold fronts (SCFs). Our model predicts that SCFs (i) appear in the cluster centre first, (ii) move outwards on several Gyr timescales, (iii) form a staggered pattern on opposite sides of a given cluster, (iv) each move outwards with approximately constant speed; and that (v) inner SCFs form discontinuities more easily than outer ones. These features are well known from idealised (magneto)-hydrodynamic simulations of cluster sloshing. We perform comparison hydrodynamic+N-body simulations where sloshing is triggered either by an instantaneous kick or a minor merger. Sloshing in these simulations qualitatively behaves as predicted by the toy model. However, the toy model somewhat over-predicts the speed of sloshing fronts, and does not predict that inner SCFs emerge with a delay compared to outer ones. In light of this, we identify the outermost cold front, which may be a 'failed' SCF, as the best tracer of the age of the merger that set a cluster sloshing.


INTRODUCTION
Mergers between galaxy clusters leave observable features in the Xray emitting intracluster medium (ICM) of a galaxy cluster.These features include shocks and cold fronts (Markevitch & Vikhlinin 2007).Cold fronts differ from shocks in that the pressure is continuous across a cold front, so that the denser side of the discontinuity is colder than the more diffuse side.
Here we focus on sloshing cold fronts (SCFs), which arise when the ICM of a cluster is perturbed by, e.g., a minor merger as first proposed by Tittley & Henriksen (2005); Ascasibar & Markevitch (2006).They showed that the gravitational disturbance caused by a subcluster passing through the primary cluster is sufficient to cause the ICM of the primary cluster to 'slosh' about the gravitational potential minimum, leading to the familiar arc-shaped 'edges' in Xray surface brightness, wrapped around the cluster core.SCFs have been observed in many galaxy clusters (for a review see ZuHone et al. 2016) and are thought to be ubiquitous in cool-core (CC) clusters (Markevitch et al. 2003;Ghizzardi et al. 2010).
The positive entropy gradient in the ICM leads to the ICM being stable against convection, i.e. after a perturbation a parcel of ICM oscillates rather than keeping rising or sinking.The frequency of such a radial oscillation is known as the Brunt-Väisälä frequency (Cox 1980), and can be written as: where Ω K = √︃    3 is the Keplerian frequency,  = 5/3 is the ratio of specific heats and  =   −2/3 the entropy index.Churazov et al. (2003) and Su et al. (2017) have used the BV period, as an estimate of the sloshing timescale.
In this paper we present a toy model that links local oscillations of ICM parcels with their BV period to the global motion of SCFs, following in broad terms the scenario outlined in Churazov et al. (2003).In essence, our toy model draws an analogy between sloshing fronts and pendulum waves (e.g.Flaten & Parendo 2001).To this end, in Section 2 we describe the basic toy model, summarising its predictions in Section 2.6.Section 2.7 discusses extensions to the basic toy model, in particular variations of the initial perturbation.Section 3 introduces the hydrodynamic simulations of sloshing resulting from an instantaneous kick and from a minor binary cluster merger, and compares the motion of the sloshing fronts in the simulations with the toy model predictions.Sections 4 and 5 discuss limitations and implications, and summarise the results, respectively.A more detailed account of the phenomenology of sloshing in terms of the linear perturbations of a cluster atmosphere is presented in Nulsen et al. (in prep.).

THE BASIC TOY MODEL
As a first step whose result is needed later, we show that the BV period,  BV , in a galaxy cluster is approximately a linear function of radius.To this end, we write the BV period (Equation 2) using the Kepler speed  K = Ω K : In galaxy clusters, it is established empirically that both the Kepler speed  K and the logarithmic derivative of the entropy index  ln  ( )  ln  are approximately constant with radius.If  is a power law  ∝   , its logarithmic derivative is its power index .Theoretically derived and observed entropy power law indices are 1.1 to 1.2 (Tozzi & Norman 2001;Voit 2005;Cavagnolo et al. 2009).Thus, we can write the BV period as or  ≈ 0.15  .
In the last step we made use of the fact that in a hydrostatic cluster, the Kepler speed  K is comparable to the sound speed   .We write the proportionality constant as 1/, as the quantity  will turn out to be the characteristic sloshing front speed.Later in this paper we present hydrodynamic sloshing simulations for a model cluster.In Figure 1, we compare the BV period of our model cluster, calculated by the full equation (Equation 1), with the approximation from Equation 4.

A row of simple harmonic oscillators whose period depends on the position of their equilibrium points
Our model starts with a galaxy cluster potential filled with a hydrostatic ICM.In the simplest version, we imagine the ICM to be perturbed by an instantaneous, unidirectional 'kick', i.e. all ICM is given an initial small velocity in the same direction.For the sake of the toy model, we first focus on the resulting motion of ICM parcels along the cluster radius, , parallel to the kick direction, such that the kick is directed outwards, in positive -direction.As a result of the kick, each ICM parcel along this radius will oscillate radially due to buoyancy with its local Brunt-Väisälä (BV) period (Equation 4).Due to the radial variation of the oscillation period, patterns of compression and rarefaction regions will arise.In this toy model, we assume that the ICM parcels simply oscillate locally without influencing each other like in a pendulum wave experiment (Flaten & Parendo 2001), which is a simplification of the fluid nature of the ICM.  4. Solid black line: Sloshing timescale calculated from the BV period (Equation 1) for our M 200 = 5 × 10 14  ⊙ cluster as a function of its radius out to  200 (1.67 Mpc).Dash-dotted black line: A linear regression fitted to the sloshing timescale which yields a characteristic sloshing front speed of  = 0.133 Mpc/Gyr (see Equations 12 and 4).Dotted green line: the linear approximation from Equation 4 assuming a single sound speed within  200 , which yields a sloshing speed of  = 0.124 Mpc/Gyr.Dotted blue and red lines and shaded region: the linear approximation from Equation 4 using the maximum and minimum values of the sound speed within  200 .
To visualise the resulting patterns, we imagine a row of simple harmonic oscillators along this radius.To start with, we assume that each oscillator is kicked such that all oscillators have the same positive amplitude, .Variations to this perturbation are discussed below.
The displacement of each oscillator away from its equilibrium position at time, , shall be  (); its period is the BV period, assumed to depend linearly on radius as shown in Equation 4. Thus, we can write the displacement, , away from equilibrium of the oscillator with equilibrium position, , at time, , as If indeed all oscillators receive their first kick at the same time,  = 0, the phase is zero ( = 0).In a more general case, the different oscillators could start at different moments in time, and thus have different phases; we discuss this case below.With Equation 4, the displacement, , of the oscillator with equilibrium position, , at time, , becomes (6) its dependence on radius for a fixed time, , is shown in Figure 2. The 1/ term in the argument of the sine function distorts the sine function such that its 'wavelength' decreases with decreasing .The displacement  () approaches zero for  > 4, and its outermost -axis intercept is at  = 2.

Relating sloshing cold fronts to the row of oscillators
Once the oscillators start their oscillation, the variation in period along  will lead to a pattern of enhanced, and reduced, densities   6).
of the oscillating parcels.The locations of enhanced densities of oscillating parcels mark the locations of actual SCFs.These locations occur close to those -intercepts (zeros) of  () where  () has a negative slope.The zeros can be found easily by setting the argument of the sine function to multiples of , i.e. the zeros occur at but only the even-numbered ones are those with a negative slope.The oscillators with equilibrium points just left of these zeros have moved towards the right (outwards), and the oscillators with equilibrium points just right of these zeros have moved towards the left (inwards).Consequently, the density of oscillators is enhanced around those zeros.This pattern of alternating inwards and outwards motion is well-known in hydrodynamic simulations of sloshing; the SCFs are located where outward-moving ICM meets inward-moving ICM, as described above.
A second method to visualise the enhanced density of oscillators is to consider the actual positions of the oscillators at a given time .The position of the oscillator with equilibrium position, , at time, , with respect to the cluster centre is i.e. its equilibrium position plus its local displacement.We show the function () along with  (), both zoomed into a relevant range of , in Figure 3 in the left-hand-side column.SCFs will occur where the density of oscillating parcels is enhanced.This is at the radii where the slope of () has local minima, which is at, or near, the zeros of  () with negative slopes.These radii are marked by cyan lines in Figure 3.However, at smaller radii, neighboring enhanced regions start overlapping, which washes out the enhancements again.This would be apparent if we took the black markers in the bottom left panel of Figure 3 (which are equidistant in their equilibrium position) and project them onto the -axis.The result of such a projection is shown by the band of blue dots on the right of the panel.Instead of projecting the black markers directly onto the -axis, i.e. setting all their -coordinates to zero, we projected them into a small -range in the margin of the plot, giving each a random -coordinate in that range to avoid crowding.There are bands of clearly enhanced densities: these mark the locations of SCFs.
As stated above, due to overlap at inner radii, there are only a limited number of regions with actual enhanced oscillator density.This washing-out of inner SCFs will be reduced with a more realistic initial perturbation, as discussed below.As real ICM parcels cannot cross through each other, this washing-out effect might not happen in reality.The outermost SCFs do not suffer from crowding and should always exist, though they may not be discontinuities, as discussed below.

A more realistic perturbation: constant initial velocity instead of amplitude
The velocity amplitude of the oscillation stated in Equation 6 would be i.e. it would depend on  because the oscillation period depends on .This would mean oscillators at lower  would have a much higher energy if they had equal masses.A more even energy distribution with radius would be more realistic for the ICM of a cluster.Under the impulse approximation, i.e. assuming the perturber passes a region faster than the matter can respond, the perturbation caused is a constant velocity kick if the perturber is an isothermal sphere.
To achieve a constant velocity amplitude, the oscillation amplitude, , needs to be a matching function of radius, : where  is a positive dimensionless parameter (signifying an outwards kick), and consequently the velocity amplitude is constant: The above relationship links the proportionality constant  for () to the initial kick velocity if we want to think of a scenario where the initial condition is a constant speed for all oscillators instead of a constant amplitude.Thus, the position of each oscillator now is This new displacement,  (), and new oscillator position, (), are shown in the right-hand-side column of Figure 3. Again, the oscillators will pile up where the slope of () has its local minima.
Differentiating () with respect to  yields its slope, differentiating again to find the minima of the slope yields , which has zeros where the argument of the sine function equals integer multiples of .Thus, minima and maxima of the slope of () occur again at but only the even-numbered instances are minima, i.e. locations of sloshing fronts.We mark them with cyan lines in Figure 3.
In the blue-dot-band next to the bottom-right panel of Figure 3, we repeat the exercise of visualising the enhanced oscillator densities at SCFs.Making the amplitude a linear (or, more generally, monotonically growing) function of  strongly reduces the washing out of inner SCFs.In the shown example, 7 to 8 SCFs can be identified instead of only 3 fronts in the left-hand-side panels.The left-hand-side panels are for the case of constant oscillator amplitude, the right-hand-side panels illustrate the case where the amplitude grows linearly with radius.The top panels show the displacement  ( ) of the oscillators around their equilibrium position (compare to Figure 2).The bottom panels show the distance of the oscillator positions to the cluster centre at a given time.The bar of blue dots on the right in the bottom panels visualises the enhanced densities of oscillators at certain radii, see text for full description.The cyan lines in all panels mark locations of potential sloshing cold fronts, i.e. locations of enhanced densities of oscillators, i.e. locations of zeros of  ( ) with negative slopes.

The location and motion of sloshing cold fronts
We established that along the radius along which the initial kick was directed outwards, SCFs appear at the even-numbered instances of To consider the case of a kick in the negative -direction, i.e. the other side of the cluster in the case of a cluster-wide unidirectional kick, we need to invert the sign of the amplitude, , in Equations 6 and 7, and of the parameter  in Equation 10.Finding again the locations of minimum slope in () now identifies the odd-numbered instances of Equation 11.Thus, Equation 11 lists the locations of potential SCFs on both sides of the cluster, counting SCFs from the outermost one inwards, alternating sides of the cluster.The sloshing fronts form a staggered pattern.The outermost front is expected on the side of the cluster that experienced the inwards kick.
Each SCF moves outwards with a constant speed, the inner fronts move slower than outer ones, and the pattern of  () remains self-similar.As the characteristic speed, , is much smaller than the sound speed, all SCFs move subsonically.We note that in the toy model framework, sloshing cold fronts are simply a pattern of local enhancements that is travelling through space, similarly to the patterns seen in a pendulum wave experiment (Flaten & Parendo 2001).

True and 'failed' sloshing fronts
At the locations of the potential cold fronts identified by Equation 11, the slope of the function () decreases with decreasing .For the outermost potential fronts, the slope of () can still be positive, but for more inner potential fronts it is negative.In the framework of the toy model, a region of negative / means oscillators have crossed through each other, whereas at the outermost fronts where the slope of () is positive, the oscillators have not changed their order but simply moved closer together.A profile of oscillator density as a function of radius shows a continuous enhancement at an outer potential front when / > 0, but the oscillator density shows a discontinuous enhancement if / < 0 (see Figure 3).It is known from hydrodynamical sloshing simulations that for mild mergers the outermost sloshing 'fronts' can fail to become discontinuities, whereas the inner fronts are discontinuous.Thus, by analogy, in the toy model we identify potential sloshing fronts with / > 0 as 'failed' fronts.True, discontinuous fronts require / < 0.
We note that in a 1-D scenario, ICM parcels would not pass through each other, and with adiabatic processes alone the gaseous ICM would not form discontinuities.However, we know from observations and hydrodynamic simulations that in 3-D discontinuous fronts form.The toy model alone cannot explain the exact process, though.Within the toy model framework, the formation of true, discontinuous fronts occurs where the slope of () not only has a local minimum, but the slope is zero or negative at that minimum.Equivalent considerations for both sides of the cluster lead to the following condition for true, discontinuous fronts: This is progressively easier to fulfil for potential fronts with higher , i.e. closer to the cluster centre.If this condition is not fulfilled, we expect the SCF to not be a discontinuity but only a gradient in ICM density and temperature in the correct direction, and we call this a 'failed' SCF.

Predicted behaviour of sloshing cold fronts from basic toy model
In summary, our toy model made the following assumptions: • We considered the ICM along one diameter in an initially hydrostatic cluster.
• All ICM parcels along this diameter simultaneously receive a kick, i.e. a small, unidirectional, initial velocity.
• As a result, the ICM parcels will oscillate locally as simple harmonic oscillators along the diameter around their equilibrium radius with their local BV period.Their amplitude shall be small compared to the radial range of interest (this defines 'small kick' in the assumption above).
• The BV period depends linearly on cluster radius, see Equation 4.
We showed that the dependence of the oscillation period on radius leads to density enhancements appearing in the ICM.We identified those as (potential) SCFs.
Based on these assumptions, this toy model predicts the following behaviour of SCFs: (i) The radii of SCFs on opposite sides of the cluster make a staggered pattern.(ii) At any given time, SCFs are located at the radii given in Equation 11(the index counts inwards), i.e. their radii keep a self-similar pattern over time.(iii) Each SCF moves outwards with constant, clearly subsonic speed (Equations 12 and 4).(iv) Inner SCFs move slower than outer ones (Equation 12).(v) For every time, , there is an outermost SCF, i.e. despite the assumed cluster-wide perturbation, the SCF pattern will grow from the cluster centre outwards.(vi) Not all sloshing fronts identified by Equation 11 are true fronts.
In particular the outer 'fronts' could fail to become true discontinuities.Forming true discontinuities requires a sufficiently strong initial kick velocity as specified in Equation 13.This condition is easier to fulfil for inner sloshing fronts.However, even a very mild kick velocity will lead to a sloshing-front-like pattern, except that the classic discontinuities are replaced by corresponding slopes in density and temperature.
The qualitative aspects of these predictions are well known SCF features in hydrodynamic simulations.

Variations to the toy model
An obvious question is whether the chosen initial perturbation impacts the prediction.To this end, we discuss some variations to the initial perturbation, namely: • an initial offset instead of an initial velocity, • a constant oscillation amplitude throughout the cluster instead of a constant oscillator velocity or energy, • a non-simultaneous perturbation where the perturber's velocity through the cluster is much faster than the characteristic speed  identified above.
This section reveals that all qualitative conclusions are unaffected, and even quantitative results for sloshing front locations, and speeds, are very similar.The strongest impact could arise from the second point.The constant oscillator velocity throughout the cluster favours the appearance of numerous true fronts, whereas in the constant amplitude case inner cold fronts could be washed out by overlapping each other, although the toy model cannot predict how gas parcels would behave in this scenario.

Initial offset instead of initial kick
If we consider the case of an initial offset instead of an initial kick, the equation describing the position of each oscillator as a function of its equilibrium position (equivalent of Equation 10) becomes Here we have kept a radius-dependent initial offset or amplitude, i.e. the maximum oscillation velocity (and energy) of each oscillator is the same.Considerations equivalent to the ones above reveal that now potential sloshing cold fronts are expected at locations (equivalents to Equations 11 and 12) and their speeds are Again, sloshing fronts are numbered from the outermost one inwards.
Odd numbered fronts appear on the side where the offset perturbation was directed towards the cluster centre, even-numbered fronts on the other side.The difference of front speeds between the two perturbation modes becomes less with increasing , i.e. for inner fronts.For front number  to be a true front, i.e. an ICM discontinuity, the kick velocity needs to obey (equivalent of Equation 13) Thus, all qualitative conclusions remain, and quantitative conclusions change only mildly.

Constant amplitude throughout cluster instead of constant velocity/energy
If the initial perturbation would lead to oscillators at different radii having the same amplitude rather than the same energy, we would expect fewer true cold fronts.Inner cold fronts would wash each other out easily (see the bottom left panel of Figure 3).In the case of constant oscillator energy throughout the cluster, the resulting radial growth of amplitude reduces fronts being washed out at smaller radii, and supports fronts being true fronts at larger radii.

Non-simultaneous perturbation
We return to the case of an initial kick.So far we considered the case that the perturbation occurs at the same time throughout the cluster.
We relax this condition now by expressing the displacement of an oscillator with equilibrium position  at time  as i.e. we use a location-dependent delay time, ().For  < ,  shall be zero.As a simple case, we write the location-dependent delay time () as This describes a perturber arriving at cluster radius  max =    0 at  = 0, which takes the time,  0 , to travel to the cluster centre, travelling with constant speed   inwards.
Inserting the delay function,  (), into the displacement function yields The zeros, i.e. the locations of potential sloshing fronts on alternating sides of the cluster, can be derived as above, and are and their speeds are For a typical cluster, the velocity, , characterising the dependence of BV period of radius, is only about 15% of the sound speed (Equation 4), whereas the infall velocity of a subcluster, i.e. a perturber, is easily 1.5 times the sound speed.Thus, the characteristic sloshing front speed, , is at least 10 times smaller than the typical speed of a perturber crossing the cluster, and /  is small.Thus, if we shift into the time frame t =  −  0 where the perturber arrives in the cluster centre at t = 0, the positions and speeds of the potential cold fronts are only slightly larger compared to the fully instantaneous perturbation approach.The effect is largest (of the order of 10%) for the outermost front.All qualitative conclusions remain the same.However, we note that this scenario still assumes a locally instantaneous perturbation, and not a perturbation over an extended amount of time or region.

BV period not a linear function of radius
The outward motion of SCFs will occur, even if not at constant speed, as long as the BV period is a monotonically increasing function of radius.The positions and speeds of sloshing fronts can be calculated by the same formalism as above but may require a numerical solution.

Simulation Method
In order to test the efficacy and predictive power of the toy model, we perform a set of three highly idealised simulations (dubbed Kick1, Kick2 and Kick3) in addition to an idealised binary merger simulation for comparison.We initialise a spherically symmetric cluster (M 200 = 5 × 10 14  ⊙ , r 200 = 1.67 Mpc) in hydrostatic equilibrium.Details of the method used to generate the cluster used in these simulations can be found in Vaezzadeh et al. (2022) which follows the methods of ZuHone ( 2011).The particles are set up to form a dark matter (DM) halo, as explained in Vaezzadeh et al. (2022), at rest in the grid, i.e the particles are not given any bulk velocity.The gas in our simulation domain is initialised with a uniform initial velocity to the right which we vary between our three simulations.This method is similar to the one used by Churazov et al. (2003) who used a planar shock front running over a cluster to initiated sloshing.
The simulations are run using the hydrodynamic + N-body code, FLASH v4.6 (Fryxell et al. 2000).FLASH is an Eulerian adaptive mesh refinement (AMR) hydrodynamics code which allows us to save computational effort in areas of the simulation domain that are of little interest.We use FLASH's N-body solver with 5 × 10 6 particles to realistically capture the response of the cluster potential to the induced gas sloshing.We use particle density to refine our domain: when the number of particles in a block (16 3 cells) exceeds 1750, the block is refined, and conversely when the number of particles in a block falls below 1500, the block will be de-refined.This allows us to achieve a resolution ranging from ∼ 9.76 kpc within a radius of ∼ 1 Mpc of the cluster core to ∼ 2.44 kpc within a radius of 0.22 Mpc.We run the simulations in a domain of 10 Mpc 3 in size with diode (isolated) boundary conditions.We allow the simulations to run for ∼ 10 Gyr with snapshots produced every 50 Myr.For simplicity we do not take account of cosmological expansion in the simulations, nor do we include radiative cooling or viscosity.
To automatically detect, and thus track, the SCFs in our simulations we use the SCF detection algorithm detailed in Vaezzadeh et al. (2022), interpreting changes in the temperature profile of > 2% over a radial range of 1.3 kpc as SCFs, and discarding those fronts that have a temperature ratio between the start and end points of < 10%.Once SCFs have been detected in this way, we relax the criterion that the temperature ratio across the front be > 10% in order to carefully trace their evolution as far back in time as possible.For our analysis of the simulations we use the Python based library, yt (Turk et al. 2011).

Qualitative evolution of the Kick simulations
Figure 4 shows a series of snapshots for each of the Kick simulations (top 3 rows).The sloshing process evolves very similarly despite the different kick strengths.The motion of the ICM, after the kick along the -axis, results in a sloshing pattern that is clearly orientated along the -axis.The arc-shaped sloshing fronts have an angular extent of approximately 180 • .Because the initial perturbation is perfectly axisymmetric (i.e.contains no angular momentum), the sloshing pattern does not feature the characteristic 'one-arm spiral' pattern so often seen in binary merger simulations.The overall staggered pattern of SCFs is present with approximately the same size and staggering pattern in all three Kick simulations.Indeed, at t max , corresponding SCFs across each simulation have positions within ∼ 20 kpc of one another.
In the case of Kick1, the first discontinuity emerges at ∼ 0.2 Gyr, travelling to the left (in the opposite direction to the kick).However, this front is no longer detected by the detection algorithm beyond ∼ 100 kpc from the core, and is no longer visible by ∼ 0.5 Gyr, as seen in Figure 4.The second front emerges to the right of the core at ∼ 0.4 Gyr (as seen in the left hand column (0.5 Gyr) of Figure 4), and the third emerges at ∼ 0.7 Gyr to the left of the core.The low initial gas velocity (in line with the small perturbation assumed by the toy model) has proven insufficient for these cold fronts to develop into true contact discontinuities, i.e. they are 'failed' fronts, and thus they cease to be detected by the detection algorithm at ∼ 1.85 Gyr, and ∼ 4.8 Gyr, respectively, though the structure is visible beyond these times.Despite these fronts not having sufficient temperature jumps to be detected by the algorithm, we include them in subsequent analysis comparing toy model speed predictions to our simulations in order to maintain consistency between simulations.At t max , Kick1 features eight SCFs (eleven if one includes the three 'failed' outer SCFs).
Kick2, which has an initial gas velocity twice that of Kick1, proceeds in much the same way as Kick1 with the first front emerging at ∼ 0.15 Gyr, and ceasing its evolution at ∼ 0.4 Gyr.The second front emerges to the right of the core at ∼ 0.45 Gyr, with subsequent fronts emerging on alternative sides of the core until  max .A clear, staggered pattern of SCFs about the core is visible by t max , with little to no instability seen affecting the edges of the fronts.There are nine SCFs visible at t max (ten if one includes the initial failed front), with one less small radius front as compared with Kick1.
In the case of Kick3, the first cold front emerges at ∼ 0.1 Gyr, and travels left (opposite direction to the kick) but, as in the other simulations, it is no longer visible by ∼ 0.5 Gyr.The second front then emerges at ∼ 0.6 Gyr to the right of the core.From this point on, fronts continue to emerge on alternating sides of the core which continue to grow until t max , at which point the system features eight SCFs (nine including the initial failed front).The large initial perturbation velocity (approximately half of the ambient sound speed of the cluster) causes sufficient disruption to the cluster that it can no longer be considered a cool-core cluster.One would expect that this would stop SCFs emerging from the core of the cluster; however, fronts do continue to emerge, though there are fewer at small radii than Kick1 and Kick2.The large velocity also leads to clear instabilities along the edges of the cold fronts, with prominent Kelvin-Helmholtz instabilities (KHIs) visible.It is interesting to note that these KHIs do not disrupt the SCFs sufficiently to hinder the fronts' growth and visibility.

Tracking cold front position over time
Because the sloshing occurs along the -axis (due to the perturbation being along the -axis), we limit our analysis of SCF positions to their position along the -axis.The toy model predicts that the first front should emerge in the opposite direction to the perturbation, i.e. the first CF should appear in the negative -direction.Figure 5 shows the positions of SCF detections throughout the duration of each simulation, with each SCF shown by a different colour.As the first cold front in all three simulations fails at an early stage, we do not attempt to track its evolution beyond ∼ 0.5 Gyr.As the next two outermost 'failed' SCFs in Kick1 have ceased being detected automatically beyond ∼ 1.85 Gyr, and ∼ 4.8 Gyr, respectively, but are still visible by eye in temperature slices, we simply add their final positions manually by inspection of the slice images, such that Figure 5 captures their full motion.
SCFs move outwards with almost constant speeds, as predicted by the toy model.The front speeds also decrease with each subsequent front that emerges, also in agreement with the toy model.We then perform a linear regression on each individual SCF in Figure 5 within each simulation in order to extract speeds for each SCF. Figure 6 summarises SCF speeds vs. front number for the different simulation runs and the toy model.
According to the toy model, each SCF should emerge from the core of the cluster (i.e. a radius of zero), and therefore a linear regression should be performed with the stipulation of a null -intercept; however, leaving the -intercepts as free parameters clearly gives the better fits to the overall motion.The implied delayed emergence of inner SCFs is discussed below.

Global features
The toy model correctly predicts the outwards motion of sloshing fronts along the axis of initial perturbation, including the staggered pattern in the direction parallel and anti-parallel to the initial per- turbation.As predicted by the toy model, the sloshing fronts move outwards with approximately constant speed.
A substantial difference to the toy model is that the fronts 'emerge' from the cluster centre one by one, each with a clear delay to the previous one of the order of 1 Gyr, with the delay increasing with each subsequent front.Figure 7 summarises these delay times for all simulations.There is a clear trend in which the emergence time of each SCF is increasingly delayed as the kick strength is increased.In contrast, the toy model predicts that the whole front system arises at once, and moves outwards in a self-similar fashion.If this was the case, all graphs of front radii as a function of time should start from  = 0 at  = 0, but this is not the case in the simulations.The toy model correctly predicts that lower perturbation speeds lead to 'failed' outer cold fronts, i.e. features that are temperature gradients in the correct direction, but are not discontinuities.For example, in the case of Kick1 (the weakest perturbation case) the three outermost fronts are 'failed' fronts, whereas the fronts further inwards are discontinuities.This behaviour is predicted by the toy model qualitatively, and even quantitatively: according to Equation 13, for the characteristic speed  = 134 km/s, and kick speed 100 km/s, only fronts with  ≥ 2 × 134/100 = 2.68 should be discontinuities, which therefore predicts one more successful front than is observed.Performing this calculation for Kick2 results in a prediction of one failed front, which agrees with the simulation.In the case of Kick3, the toy model predicts no failed fronts, but in the simulation the first front fails.

Cold front speeds and positions
Figure 1 shows the sloshing timescale as calculated from the Brunt-Väisälä frequency as a function of radius for our model cluster.Clearly the sloshing timescale is monotonically increasing with radius in an approximately linear fashion within  200 .Departures from linearity in the function will cause deviations from the simple model outlined here.We do not see any SCFs travel beyond  200 , and so we limit our analysis to within this radius, where the linear approximation is appropriate.
From Equation 4, we obtain the characteristic sloshing front speed in this cluster of  ≈ 134 km/s.We can then scale this speed via Equation 12 to predict a unique speed for each of the SCFs that emerge during the course of the simulations.As explained above, the counter, , counts the sloshing fronts from the outermost one inwards.Odd fronts arise in the direction opposite to the kick directions, i.e. along the −-direction, and even-numbered fronts arise in +-direction.The solid black line in Figure 6 shows the CF speed as a function of front number as predicted by the toy model.The coloured lines show the CF speeds derived for the Kick simulations from Figure 5 and for the binary merger from Figure 8.
Figure 6 reveals that, overall, the CF speeds in the Kick simulations are about a factor of 2-3 below the predicted value, and that the front speed depends on front number in a similar power law fashion as predicted (power −1.2 instead of −1).
Differences between the Kick simulations occur for the outermost front, and for fronts beyond number 7. The outermost front in Kick1 was a 'failed' front, and moves slower than the one in the other Kick simulations.The speed of the outermost CF between Kick2 and Kick3 agrees well.The deviation from the predicted pattern at higher CF numbers could arise because the sloshing process changes the inner entropy profile of the cluster, and thus the conditions of the initial state, assumed by the toy model throughout, are not true anymore.This behaviour is increased by the fact that inner cold fronts indeed arise with a delay.More deviations from the initial entropy profile are expected to arise with increasing time, and with increasing perturbation strength -both are seen in Figure 6.
Figure 7 summarises the emergence delay times of SCFs in the different simulations.This behaviour is not predicted by the toy model.These delay times were derived by taking the earliest time at which each SCF could be seen, and aligning these times relative to the perturbation time (0 Gyr in the Kick simulations, and 1.6 Gyr in the binary merger).Given these unpredicted delay times, the toy model alone will not be able to correctly predict the positions of sloshing fronts.While the front speeds could be calibrated, the toy model does not predict the delay in the inner cold fronts emerging.Given their slow speed, this delay has a big impact on the actual CF position at a given time.

Comparison to idealised binary cluster merger
In order to test our model in a more realistic ICM sloshing scenario, we compare our toy model's predictions to a binary cluster merger with mass ratio,  = 1 : 10, in which the primary cluster is the same cluster as in the other simulations presented in this paper, with the exception of it having a 'WHIM'-like atmosphere beyond 2.17 Mpc.This 'WHIM' is a uniform gas background with density, temperature, and pressure values of 1.03×10 −29 g/cm 3 , 1.70 keV, and 2.85 × 10 −14 erg/cm 3 , respectively.The 'WHIM' has the effect and purpose of pre-truncating the atmosphere of the infaller to avoid it carrying too much gas into the primary.The infalling subcluster has a mass of 5×10 13  ⊙ , with an  200 of 777 kpc, and a particle resolution of 5 × 10 5 particles.The subcluster is initialised at 2.45 Mpc north of the primary (the sum of the respective  200 radii), with a radial velocity of −950 km/s, and a tangential velocity of 450 km/s such that the subcluster will pass to the right of the primary with a large pericentre distance, and thus deliver a 'kick' along the same axis as in the simulations presented in the previous section.

Considerations regarding the nature of the perturbation
It is important to note that the key difference between the Kick simulations presented in the previous section, and the binary merger, with regard to the toy model, is the initial perturbation that the cluster receives.The perturbation is continuous, and non-constant in space and time in the binary merger, i.e. the infaller is already perturbing the primary cluster when the simulation begins, and continues to do so as it moves through the primary.Furthermore, the binary merger introduces angular momentum into the host cluster.The perturbation may be a mix of extended kicks and offsets in the primary's ICM.It is also not obvious which direction of perturbation matters moston first approach, the perturber attracts the primary's ICM towards it, but then pushes and pulls the primary's ICM somewhat along its orbit after it passed a given location.This has important implications for the application of the toy model, as the toy model assumes a single 'kick' or offset to a row of oscillating gas parcels which then oscillate independently.In the binary merger case, the perturber first moves approximately parallel to the -axis, but its second passage through the cluster occurs more in a diagonal direction from the −, − quadrant to the +, + quadrant.
We pointed out that the perturbation in a minor merger is continuous.However, for the sake of comparison, we simplify this scenario to the often-invoked model of sloshing in which the first pericentric passage is the key moment of perturbation.As the infaller passes to the right of the primary in our binary merger simulation (therefore pulling the primary to the right at pericentre time), the 'kick' is to the right, and therefore sloshing will occur along the -axis.We note that it is significantly more difficult to trace each SCF for its full evolution than in the case of the highly idealising sloshing simulations.Due to the highly 'messy' nature of the SCFs' evolution in the binary merger, automatic tracking is more difficult, and is therefore augmented by manual tracking of the SCFs.Due to the angular momentum imparted by the infaller, the primary's ICM is swirling at the same time that it is sloshing, and as such fronts that emerge along the -axis to a given side rotate around to the other side of the cluster in some cases.Once the coherent SCF points have been identified, the same procedure of linear regression is applied as described in Section 3.3.

Comparison of binary merger with Kick simulations and toy model
Figure 8 shows the CF positions as a function of time along the direction.Again, to each CF we fit a linear position-time function to determine the speed of each CF. Figure 6 compares the CF speeds from the binary merger to the toy model prediction and the Kick simulations.Figure 7 compares the delay times of the emergence of the SCFs in the binary merger to the delay times for the Kick simulations.
Despite the significantly more complex perturbation in the binary merger, sloshing fronts also arise in a staggered pattern, and move outwards with about constant speed.Similar to the Kick simulations, the inner sloshing fronts emerge with a delay, which is not predicted by the toy model.The dependence of delay time on SCF number is approximately linear, and of a similar order to that seen in Kick3 for early SCFs, and Kick1 for later SCFs.
The SCF speeds (Figure 6) are again within a factor of 2 of the toy model prediction.For SCF numbers of 5 and below, the SCF speeds show a similar dependence on SCF number as the ones in the Kick simulations.For outer SCFs, there is an approximate agreement with the 1/ dependence, although there is significantly more scatter around that trend.From CF 6 onwards, the SCF speeds increase again, even becoming faster than predicted by the toy model.
It is worthwhile noting that, despite the slightly lower overall evolution time in the binary merger (the 'kick' at pericentre occurs at 1.6 Gyr as opposed to 0 Gyr), there is an equal number of SCFs in the binary at  max to in Kick1.The ICM bulk motions in the primary cluster triggered by the merger are about 600 km/s, which is comparable with the kick speed of simulation Kick3.Thus, the perturbation in this binary merger is not a particularly weak perturbation compared to the Kick simulations, especially given that it is a continuous perturbation.

DISCUSSION
We presented a simple toy model for the sloshing process in galaxy clusters.The toy model assumes that sloshing arises as the ICM parcels in an initially hydrostatic cluster start oscillating around their equilibrium radius with their local Brunt-Väisälä (BV) period after an initial perturbation.We showed in Equation 4 that the BV period in galaxy clusters can be approximated as a linear function of radius.The variation of the BV period with radius leads to a characteristic pattern of density enhancements in the ICM that can be linked to sloshing cold fronts.These enhancements travel outwards with about constant speed for a range of initial perturbations (instantaneous kick, instantaneous offset, instantaneous local perturbation that travels through the cluster supersonically).We compared the toy model's prediction in detail to hydrodynamic+N-body simulations of both sloshing initiated by an initial kick to the ICM, and sloshing caused by a binary minor cluster merger.

Successes of the toy model
The toy model correctly predicts several key qualitative characteristics of the sloshing front system in a galaxy cluster: • Sloshing fronts arise from the centre of the cluster, and move outwards, i.e. at any time there is an outermost SCF.The physical way of numbering SCFs should start with this outermost front.
• Sloshing fronts form a staggered pattern, i.e. appear in an alternating fashion on opposite sides of the cluster core.
• Sloshing fronts move outwards with about constant speeds.
• For the outer few SCFs, the front speed decreases with front number, , approximately in proportion to 1/.• Weak perturbations can lead the outermost fronts to be 'failed' sloshing fronts, i.e. they are not discontinuities but have sloshing cold front characteristics in every other aspect.
The toy model quantitatively predicts SCF speeds within a factor of 2 to 3. The characteristic front speed is 14% of the cluster's sound speed (Equation 4).Thus, the toy model predicts subsonic motion of sloshing fronts.These features have been seen in hydrodynamic simulations (e.g.Ascasibar & Markevitch 2006;ZuHone 2011), including the feature of 'failed' SCFs (Roediger et al. 2011), and are reproduced in our idealised Kick simulations presented here as well as our binary merger simulation.
Given that sloshing fronts are a wave phenomenon, we cannot expect sloshing to transport matter over large distances.For example, we do not expect sloshing to transport low entropy gas from cluster centres to the outskirts.

Limits of the toy model
The first major deviation between the toy model and the hydrodynamic simulations is the delay time of inner, i.e. later, cold fronts.In simulations with either an initial instantaneous, kick-like perturbation as well as perturbation by a classic binary minor merger, our analysis showed that the sloshing fronts emerge from the cluster centre with a substantial delay time that grows with cold front number.This means that the position of a given cold front depends not only on its speed, but also on its emergence time.The toy model does not predict this emergence delay, but predicts that the whole sloshing front system emerges together, and grows in a self-similar pattern.Further investigations are required to uncover the origin of the emergence delay of later CFs, and why the toy model still predicts reasonable cold front speeds despite this mismatch.
Secondly, the toy model over-predicts the speed of the sloshing fronts by a factor of 2-3.This is related to the fact that that the toy model considers only radial oscillation modes.In the full treatment (Nulsen et al., in prep.), the oscillation frequency depends on the angle  between the wave vector and the radial direction as  BV sin , thus reducing the characteristic sloshing speed.
There are some further, expected differences between the toy model and the hydrodynamic simulations.Over time, the sloshing process alters the entropy profile of the ICM in the cluster core, thus the motion of the sloshing fronts must change.This aspect is not included in the toy model.The effect is expected to be stronger for stronger perturbations, and later (i.e. more inner) SCFs, due to a stronger resulting modification of the central entropy profile.This is indeed the case.
The version of the toy model presented here considers only instantaneous perturbations, either simultaneously throughout the whole cluster, or a locally instantaneous perturbation moving through the cluster at a speed considerably larger than the characteristic sloshing speed.We have also separated kick and offset perturbations in this version of the toy model.A binary merger causes a more complex perturbation, extended in time and space, so a perfect match cannot be expected.
The current toy model only predicts speeds and locations of cold fronts, but not their strength in terms of density or temperature contrast across them, apart from predicting potentially 'failed' outer cold fronts.

Implications for determining merger ages
There are two difficulties in deriving the age of a merger from an observation of a set of sloshing cold fronts in a given cluster.Most clusters are observed first, and best, in their central regions, i.e. we observe most easily the inner cold fronts of a sloshing front system.If we misidentify them for outer cold fronts, we will generally overestimate their speed, and thus underestimate the age of the merger.If only part of the sloshing front system is known, there is no easy way to know which cold fronts are observed.Thus, merger ages from studies interpreting only sloshing in the cluster centre can only give a lower limit on the cluster's merger age (e.g.Roediger et al. 2011Roediger et al. , 2012)).
Even using a reasonable estimate of the particular sloshing front's speed, simply tracing back the current SCF radius to the cluster centre can strongly underestimate the age of the merger because the SCF emergence delay time of approximately  × 2/3 Gyr is not included.We note that if the age of a particular front is of interest instead of the age of the merger, the simple trace-back method gives good estimates each sloshing front moves with approximately constant speed.We note that the speed of a given cold front is approximately independent of the perturbation only for mild perturbations.For, e.g., stronger mergers, the cold front speed increases (Roediger et al. 2011;Bellomi et al. 2023).
To estimate the age of the merger that caused the sloshing, we need a view of the cluster as a whole, and must identify the outermost sloshing front.This front is affected least by the delay in emergence, and is a direct tracer of the merger age.However, the search for the largest SCF must include looking for 'failed' SCFs, and not simply use the outermost front with a discontinuity.We show in a forthcoming paper that the outermost CF is indeed a good tracer of the merger's age.
The toy model implies a relationship between the orientation of the cold front system, and the merger direction.However, even a single non-head-on merger introduces a rotational component into the ICM that rotates the sloshing direction, introducing a bias in direction.This effect could depend on the impact parameter of the merger.

Sloshing fronts as a wave phenomenon and their resilience against destruction by Kelvin-Helmholtz instability
True SCFs are contact discontinuities with gas of different entropy on either side.However, the identity of the gas on either side of a given front changes with time.ICM that still is on the outside of a given front will be on its inside a while later when the front has moved further outwards.SCFs are like waves moving through the ICM, they do not transport ICM from inner to outer radii over large distances.Thus, SCFs differ in their nature from, e.g., the CF at the upstream edge of a subcluster falling into a host cluster.In this scenario, the gas on the hotter side of the front is always host cluster ICM and the gas on the colder side always subcluster gas.Gas parcels at this kind of cold front can be replaced by flows inside the subcluster or the flow of host cluster ICM around the subcluster atmosphere, but gas from each reservoir does not change to the other side of the front.In contrast, at SCFs, material changes from the hotter side to the colder side as the sloshing front moves over it.
Thus, while we expect shear flows along sloshing fronts to cause Kelvin-Helmholtz instabilities (KHIs), we should not expect these KHIs to be able to fully erase a given SCF.For example, KHIs of perturbation length 10 kpc, arising on an interface of density contrast 2 and shear velocity 300 km/s, have a growth time of 30 Myr, but take about 5 growth times to form the classic KHI rolls (e.g.Roediger et al. 2013a), and would take even longer to erase the interface.A naive interpretation could conclude that a ∼ Gyr old discontinuity should be erased by KHIs as their growth time is much shorter than the age of the front.However, assuming a sloshing front speed of 0.05 kpc/Myr (Figure 6), in 150 Myr, the sloshing front would travel 7.5 kpc, i.e. almost a perturbation length, which is typically larger than the KHI roll height.Thus, KHI growth and sloshing front propagation, and re-formation, take place on similar timescales, which explains why KHIs generally cannot erase sloshing fronts.Sloshing CFs are indeed known to survive KHIs from hydrodynamic simulations.Rather than being erased or washed out, they are only distorted (e.g.Roediger et al. 2013b;ZuHone et al. 2013).Further ICM properties stabilising fronts against KHI, e.g.viscosity (e.g.ZuHone et al. 2010;Roediger et al. 2013b) or magnetic fields (e.g.Vikhlinin & Markevitch 2002;Brzycki & ZuHone 2019;Chadayammuri et al. 2022), are not necessary to ensure front survival, but would help to slow down or even prevent the onset of KHIs.

CONCLUSION AND SUMMARY
We presented a simple toy model for sloshing of the ICM in galaxy clusters that describes sloshing fronts as a coherent pattern arising from ICM parcels oscillating locally with their Brunt-Väisälä period.This period can be approximated by a linear function of radius.The proportionality constant, 1/, is the inverse of the characteristic speed of the resulting sloshing fronts, and is about 14% of the ICM sound speed.The simple model successfully predicts the staggered pattern of sloshing fronts on opposite sides of the cluster, the outwards motion of sloshing fronts with approximately constant speed, and the finite size of the sloshing front pattern.Sloshing fronts should be numbered from the outside inwards.A careful analysis of hydrodynamic simulations reveals that in the hydrodynamic treatment, sloshing fronts emerge near the cluster centre one after the other, with a delay of roughly 0.5 Gyr between them.This effect is not captured by the toy model.
We explained that the best option to derive the age of the merger that triggered sloshing is to trace back the outermost sloshing front.However, such an analysis needs to take into account 'failed' sloshing cold fronts, i.e. those outer sloshing cold fronts that did not form a true discontinuity.The existence of such 'failed' sloshing cold fronts is predicted by the toy model.

Figure 1 .
Figure1.Agreement of the correct sloshing timescale for our model cluster used in Section 3 with the simple approximation from Equation 4. Solid black line: Sloshing timescale calculated from the BV period (Equation1) for our M 200 = 5 × 10 14  ⊙ cluster as a function of its radius out to  200 (1.67 Mpc).Dash-dotted black line: A linear regression fitted to the sloshing timescale which yields a characteristic sloshing front speed of  = 0.133 Mpc/Gyr (see Equations 12 and 4).Dotted green line: the linear approximation from Equation 4 assuming a single sound speed within  200 , which yields a sloshing speed of  = 0.124 Mpc/Gyr.Dotted blue and red lines and shaded region: the linear approximation from Equation 4 using the maximum and minimum values of the sound speed within  200 .

Figure 2 .
Figure 2. Displacement, , away from the equilibrium position of oscillators as a function of their equilibrium position,  (see Equation6).

Figure 3 .
Figure3.Oscillator positions and enhanced densities.The left-hand-side panels are for the case of constant oscillator amplitude, the right-hand-side panels illustrate the case where the amplitude grows linearly with radius.The top panels show the displacement  ( ) of the oscillators around their equilibrium position (compare to Figure2).The bottom panels show the distance of the oscillator positions to the cluster centre at a given time.The bar of blue dots on the right in the bottom panels visualises the enhanced densities of oscillators at certain radii, see text for full description.The cyan lines in all panels mark locations of potential sloshing cold fronts, i.e. locations of enhanced densities of oscillators, i.e. locations of zeros of  ( ) with negative slopes.

Figure 4 .
Figure 4. Time series of temperature slices of the simulations centred on the minimum of gravitational potential.Each row is for a different simulation, from top to bottom: Kick1, Kick2, Kick3, binary merger.The columns show the simulations at 0.5 Gyr, 1.0 Gyr, 5 Gyr after the initial perturbation, and at maximum simulation time.As perturbation time, we take 0 Gyr in the Kick simulations, and pericentre passage (1.6 Gyr) in the binary merger.Each panel is 2.5 Mpc on a side.A movie of the simulations can be found at https://www.youtube.com/shorts/n9GBAzCs5T8.

Figure 5 .
Figure 5. Positions of detected cold fronts throughout each of the simulations Kick1, Kick2, and Kick3.The coloured dots represent distinct cold fronts.Solid lines of the same colour show linear fits to the cold front radii as a function of time.

Figure 6 .Figure 7 .
Figure6.Lineplots of the linear best-fit SCF speeds from each simulation, and the speeds predicted by the toy model as a function of SCF number.The solid black line shows the toy model predictions for each SCF's speed; the solid coloured lines show the linear fit SCF speeds from each of the simulations; and the dash-dotted black line shows a power-law of 1/ 1.2 .

Figure 8 .
Figure 8. Positions of detected cold fronts, along the -axis, throughout the binary merger simulation.The data points begin at 1.6 Gyr, which is the time of first pericentre; the vertical blue line at 5.2 Gyr shows the time of second pericentre.The coloured dots represent different tracked SCFs, with the correspondingly coloured solid lines representing linear fits to those tracked SCFs.