Unravelling cosmic velocity flows: a Helmholtz-Hodge decomposition algorithm for cosmological simulations
UUnravelling cosmic velocity flows: a Helmholtz-Hodge decomposition algorithm forcosmological simulations
David Vall´es-P´erez a, ∗ , Susana Planelles a , Vicent Quilis a,b a Departament d’Astronomia i Astrof´ısica, Universitat de Val`encia, E-46100 Burjassot (Val`encia), Spain b Observatori Astron`omic, Universitat de Val`encia, E-46980 Paterna (Val`encia), Spain
Abstract
In the context of intra-cluster medium turbulence, it is essential to be able to split the turbulent velocity field in acompressive and a solenoidal component. We describe and implement a new method for this aim, i.e., performing aHelmholtz-Hodge decomposition, in multi-grid, multi-resolution descriptions, focusing on (but not being restricted to)the outputs of AMR cosmological simulations. The method is based on solving elliptic equations for a scalar and a vectorpotential, from which the compressive and the solenoidal velocity fields, respectively, are derived through differentiation.These equations are addressed using a combination of Fourier (for the base grid) and iterative (for the refinement grids)methods. We present several idealised tests for our implementation, reporting typical median errors in the order of1 ‰ -1%, and with 95-percentile errors below a few percents. Additionally, we also apply the code to the outcomes ofa cosmological simulation, achieving similar accuracy at all resolutions, even in the case of highly non-linear velocityfields. We finally take a closer look to the decomposition of the velocity field around a massive galaxy cluster. Keywords: turbulence, large-scale structure of Universe, galaxies: clusters: intracluster medium, galaxies: clusters:general, methods: numerical, Adaptive Mesh Refinement
1. Introduction
Cosmological structures and, in particular, galaxy clus-ters, which constitute the most massive structures whichhave had time to collapse under their own gravity, aredynamically interesting objects from several perspectives.The non-linearity of the evolution of their baryonic com-ponent (i.e., ordinary matter, most of which is in theform of a hot, tenuous plasma known as the intra-clustermedium, ICM) couples the different scales in the evo-lution of cosmic inhomogeneities, producing a plethoraof complex hydrodynamical phenomena, such as shockwaves and turbulence. Turbulence is an intrinsically multi-scale phenomenon, since bulk motions trigger (magneto-)hydrodynamical instabilities (e.g., [1]) which cascadedown to smaller scales until they get dissipated by viscouseffects.Turbulent motions can be [2] and have been recently[3, 4] measured on a number of nearby galaxy clusters fromX-ray surface brightness fluctuations, which have in turnbeen connected to signatures of particle acceleration anddiffuse radio emission in the ICM [5, 6]. Future X-ray fa-cilities, like
Athena , will offer unprecedented insight into ∗ Corresponding author
Email addresses: [email protected] (DavidVall´es-P´erez), [email protected] (Susana Planelles), [email protected] (Vicent Quilis) the dynamics of turbulent motions in the ICM [7]. How-ever, the precise theoretical and numerical description ofturbulence in these vast structures is still matter of on-going research (e.g., [8, 9, 10, 11], just to highlight a fewrecent numerical studies).Splitting the (turbulent) velocity field in its compressiveand solenoidal components, i.e., performing a Helmholtz-Hodge decomposition (HHD), is a crucial step towards ex-ploring the role of turbulence in the ICM, since these com-ponents play fundamental and distinct roles in the evolu-tion of cosmic structures. Thus, while the solenoidal com-ponent is likely the major responsible for the amplificationof cosmic magnetic fields [12, 13], compressive turbulencehas an important role in generating weak shocks whichhave consequential effects on the magnetic and thermalevolution of the cluster (e.g., [13, 14]). These two com-ponents also differ in their spatial distribution, the formertending to be more volume-filling [10, 15, 16], and evenin their spectrum, steeper for the compressive component[8, 17, 18, 19]. The distinction between the solenoidal andcompressive component of the velocity field is also of ut-most importance to model the acceleration of cosmic rayparticles in the ICM, which will become a vibrant field ofobservational research with the advent of the new genera-tion of radio telescopes (e.g., SKA [20]). In this line, a lotof effort has been recently put in modelling the accelerationby compressive and solenoidal modes in magnetohydrody-namics (MHD; [21, 22, 23, 24, 25]).Cosmological hydrodynamical simulations, as many Preprint submitted to Computer Physics Communications February 15, 2021 a r X i v : . [ a s t r o - ph . I M ] F e b ther applications in computer fluid dynamics, need ahuge dynamical range to be resolved in order to, for ex-ample, form realistic galaxies in a cosmological environ-ment [26, 27] or adequately describe turbulence in the ICM(see, e.g., [28] for some graphical examples of the effects ofresolution on the ability to capture instabilities; see also[29, 30], who present detailed studies of stratified, ICM-like turbulence in numerical grids of varying resolution).While Lagrangian codes are inherently adaptive, Euleriancodes based on high-resolution shock-capturing (HRSC)techniques are especially capable of handling shocks andother types of discontinuities, which are pervasive in theformation of cosmological structures ([31], for a review).That is why, among several other options (see, for exam-ple, [32] for a broad review), Adaptive Mesh Refinement(AMR) codes are especially suited for cosmological struc-ture formation.Previous studies of ICM turbulence have already im-plemented HHD algorithms. For example, several worksusing uniform grids (or fixed refinement strategies which,ultimately, allow to resolve the object of interest within aconstant resolution) perform the decomposition in Fourierspace, where it simply reduces to linear algebra projections(e.g., [10, 33]). Additionally, [10] also confront this methodwith solving a Poisson equation (in Fourier space) to finda scalar potential for the compressive velocity component,reporting more accurate results for the first method. How-ever, any of these two procedures, because of their usageof FFT algorithms, can only be applied to regular, uni-form grids, which are not the common use in cosmologicalsimulations.In this paper, we propose, implement and test a newalgorithm for performing the Helmholtz-Hodge decompo-sition in a multi-scale AMR grid, which can therefore beapplied to the outcomes of a full-cosmological simulationwithout the need of performing resimulations of specificobjects or constrained simulations. Our method decom-poses the velocity fields by solving elliptic partial deriva-tive equations (PDEs), which can be addressed iterativelyusing a wide variety of well-known algorithms (e.g., [34]).Nevertheless, although our primary focus is the applica-tion to cosmological structure formation, we emphasizethat the approach presented in this work can be directlyapplied to any application of block-structured AMR, andreadily extended to particle-based simulations through asuitable interpolation scheme.The rest of the manuscript is organised as follows. InSec. 2, we describe our method for performing the decom-position and discuss its numerical implementation. In Sec.3 we present and describe a set of tests to validate the ac-curacy of the algorithm, while in Sec. 4 we apply it to thecomplex velocity field of a cosmological simulation. Last,in Sec. 5 we summarise and present our conclusions.
2. Description of the method
The algorithm is based on the Helmholtz-Hodge decom-position (see, e.g., [35, 36, 37]), which allows to univocallysplit any velocity field in three terms, v = v comp + v rot + v harm (1)where v comp is the compressive (or irrotational, ∇ × v comp = ) velocity field, v rot is the purely rotational (orsolenoidal, ∇ · v rot = 0) velocity field and v harm is the har-monic velocity field (both irrotational and solenoidal, thussatisfying ∇ v harm = and v harm = ∇ χ , with ∇ χ = 0).The harmonic component can be shown to be identi-cally null in any domain with periodic boundary condi-tions, and therefore we will no longer consider it. Be-cause of their defining properties, the compressive compo-nent can be written as the gradient of a scalar potential, v comp = −∇ φ , while the rotational component can be de-rived from a vector potential, v rot = ∇ × A . From this, itis easy to derive that the scalar and vector potentials canbe computed, respectively, as the solutions of the followingelliptic PDEs, ∇ φ = −∇ · v (2) ∇ A = −∇ × v (3)which are formally equivalent to a set of four decoupledPoisson equations (one for φ and one for each cartesiancomponent of A ) whose sources are the divergence andthe components of the curl of the overall velocity field.Once the potentials have been obtained, the compres-sive and rotational components of the velocity field can beobtained through differentiation. Note that, in principle,it would only be necessary to find one of the potentials(either φ or A ), since the other velocity component couldbe then derived by subtracting to the total velocity (Eq. 1with v harm = ). Nevertheless, we have chosen to computeall the potentials in order to keep track of the associatednumerical errors. We have designed a code to perform such decompositionin a multi-resolution, block-structured AMR velocity field.As mentioned before, our code can be easily applied to theoutcomes of any AMR simulation (cosmological or not), oreven to a particle-based code, from which the continuousvelocity field can be defined on an ad-hoc AMR grid struc-ture through a smoothing method (e.g., a particle-mesh,[38]).In our particular implementation, for the base (coarsest)level, (cid:96) = 0, taking advantage of the periodic boundaryconditions, Poisson equations are solved in Fourier space , Note, however, that this is not a requirement of the method. Wesolve Poisson’s equations in Fourier space as cosmological simulationsof sufficiently large volumes typically implement periodic bound- F lmn .2. Poisson’s equation is then solved in Fourier space bymultiplying by the Green’s function, G lmn : ˜ φ lmn = G lmn F lmn . The Green’s function is given by: G lmn = (∆ x/ sin (cid:16) πlN x (cid:17) + sin (cid:16) πmN y (cid:17) + sin (cid:16) πnN z (cid:17) (4)where ∆ x is the cell side length and the domain isdiscretised in N x × N y × N z cells [38].3. The inverse FFT of ˜ φ lmn yields the sought potentialat the base level.In subsequent refinement levels, (cid:96) >
0, Poisson equa-tions have to be solved taking into account the boundaryconditions imposed by the coarser grids the refined patchesare embedded into. In order to do so, we use a successiveover-relaxation procedure (SOR; see, for example, [34]) onthe discretised Poisson equation.Each AMR patch is first initialised (both in the bound-ary and in the interior cells) by linear interpolation fromthe values of the potential at the best-resolved lower levelpatch available. Then, the interior cells are iteratively up-dated in a chessboard pattern as φ new i,j,k = ωφ ∗ i,j,k + (1 − ω ) φ old i,j,k , (5)where φ ∗ i,j,k = (cid:104) φ old i +1 ,j,k + φ old i − ,j,k + φ old i,j +1 ,k + φ old i,j − ,k + φ old i,j,k +1 + φ old i,j,k − − (∆ x (cid:96) ) f i,j,k (cid:105) , (6)being ∆ x (cid:96) the cell size at the given refinement level, f i,j,k the source term and 1 < ω < ω accordingto the Chebyshev acceleration procedure [34]. Aiming toavoid undesirable boundary effects due to the interpola-tion of the potential boundary conditions, we extend the ary conditions and, in these situations, solving Poisson’s equationin Fourier space is much more computationally efficient than usingiterative methods. In any case, for non-periodic domains, the baselevel can be addressed through iterative schemes, just as describedbelow for the refinement levels, if suitable boundary conditions areprovided. In the case of non-periodic domains, however, the har-monic term cannot be dropped in Eq. 1. While we cannot compute v harm by solving elliptic equations, this term can be obtained justby using Eq. 1 to solve for it, once v comp and v rot have been found.This is a good approach, as long as our algorithm precisely recon-structs the compressive and rotational velocity fields from an inputfield, which is tested in Sec. 3 and Sec. 4. patches with 3 fictitious cells in all directions, so that theseboundary conditions are enforced slightly far away fromthe region of interest.Once φ and A are known, the velocity components arefound by finite differencing the potentials as defined be-fore. We compute the derivatives (both of the velocityand of the potentials) using an eighth-order scheme with acentered stencil of up to .
3. Tests
Aiming to assess the robustness of our HHD method andits implementation, we have designed a battery of testsfocused on quantifying to which extent the code is able toaccurately identify and disentangle the compressive androtational velocities. We describe such tests in Sec. 3.1,and examine their results in Sec. 3.2.
For the four tests described below, we have first estab-lished a simple AMR grid structure, i.e., a set of patchesat different refinement levels that could reasonably mockthe ones generated in an actual simulation.We have considered a cubic domain of unit length, withorigin at the center of the box, and we have discretisedit with 128 cells. For each octant, we establish a firstrefinement patch with twice the spatial resolution coveringthe central 1 / < x, y, z < / / < x, y, z < / cells (and likewise forthe remaining 7 octants). Then we add a (cid:96) = 2 patch inthe central 1 / (cid:96) = 1 patch, and continuerecursively.In the first three tests we consider up to (cid:96) max = 10refinement levels, providing a peak resolution of 7 . × − (relative to the box size, which is normalised to 1). Figure1 presents graphically the grid structure employed in thesetests. Below, we describe the velocity fields that we haveseeded on these grids. For the first test, we consider a purely compressive ve-locity field with constant divergence, given in cartesianand spherical coordinates by the analytical expression v = ω ( x ˆu x + y ˆu y + z ˆu z ) ≡ ω r ˆu r . (7)This field has ∇· v = 3 ω and ∇× v = , and hence v = v comp and v rot = 0. It is easy to find an analytic expression The stencil length is shortened as cells get closer to the boundary. .50 0.25 0.00 0.25 0.50 x (arbitrary units) y ( a r b i t r a r y un i t s ) Figure 1: AMR grid structure for the tests in Sec. 3. The figurerepresents a slice at z = 0 .
25, with the colors encoding the highestrefinement level at each point in the x – y plane. The white, dashedlines indicate the cartesian x and y axes. for the scalar potential, φ = − ω (cid:0) x + y + z (cid:1) + C , with C any arbitrary real constant. Likewise, the vector potentialought to be A = ∇ χ , with χ an arbitrary scalar functionor, in particular, A = . We have set ω = 0 . Analogous to the previous one, we have also consideredthe case of a purely uniformly rotating velocity field, whichcan be analytically given in cartesian and cylindrical co-ordinates as: v = ω ( − y ˆu x + x ˆu y ) ≡ ω ρ ˆu φ (8)It is straightforward to show that ∇ · v = 0 and ∇ × v =2 ω ˆu z . This velocity field is generated by the potentials φ = C , being C ∈ R a free constant, and A = − ω ρ ˆu z + ∇ χ , with χ any arbitrary scalar function. As in Test 1, weset ω = 0 . We have designed a third test, aimed to assess the ef-fects of the cross-talk between the compressive and therotational components. We have considered the velocityfield: v = [sin (2 πx ) + sin (4 πy ) + sin (6 πz )] ˆu x ++ [sin (6 πx ) + sin (2 πy ) + sin (4 πz )] ˆu y ++ [sin (4 πx ) + sin (6 πy ) + sin (2 πz )] ˆu z . (9)The compressive part corresponds to the terms of angu-lar frequency 2 π , while the higher frequency ones consti-tute the rotational component. Also in this case, it is easyto find analytical solutions to φ and A : φ = 12 π [cos (2 πx ) + cos (2 πy ) + cos (2 πz )] + C (10) A = (cid:20) − π cos (4 πz ) + 16 π cos (6 πy ) (cid:21) ˆu x ++ (cid:20) − π cos (4 πx ) + 16 π cos (6 πz ) (cid:21) ˆu y ++ (cid:20) − π cos (4 πy ) + 16 π cos (6 πx ) (cid:21) ˆu z + ∇ χ. (11) While the previous tests have checked the ability ofthe code to reconstruct idealised solenoidal, compressiveand mixed velocity fields, we have implemented a last testaimed to assess the ability of the code to capture and re-construct variations on a broad range of spatial frequen-cies. The test is, in part, inspired by the one presentedby [10, App. A.1.2], but with some differences aimed tomix both velocity fields (instead of generating a purelysolenoidal or compressive field), while still having an ana-lytic solution of the HHD to compare with the numericalresults.We have generated our mock, ICM-like velocity field ac-cording to the procedure described below:1. We consider a uniform grid of ( N x × (cid:96) max ) cells. Inthat grid, we compute the velocity field v = v comp + v rot , with: v comp = (cid:88) i = x,y,z N max (cid:88) n = N min A comp n sin (cid:18) πnL x i + ψ comp , in (cid:19) ˆu i (12) v rot = (cid:88) i = x,y,z (cid:88) j (cid:54) = i N max (cid:88) n = N min A rot n sin (cid:18) πnL x j + ψ rot , ijn (cid:19) ˆu i (13)being L the box side length ( L = 1 in our case), A comp n ( A rot n ) the amplitude of the mode of frequency n ofthe compressive (solenoidal) component, and ψ comp , in ( ψ rot , ijn ) the initial phases. For the intial phases, we4ave generated 9 sets of N max − N min +1 random num-bers, uniformly sampled from the interval 0 < ψ < π .The amplitudes are generated so that the compressive(solenoidal) component follows a Burgers [39] (Kol-mogorov [40]) spectrum, E ( k ) ∝ k − ( E ( k ) ∝ k − / ).We set N min = 2 and N max = 1024, so that our mockvelocity field presents solenoidal and compressive fluc-tuations over scales differing almost 3 orders of mag-nitudes. In order for both components to be relevant,we fix the amplitudes so that A rot n = A comp n for n = 64.We note that, while the velocity field generated ac-cording to this procedures is not the most general one(e.g., one could add oblique plane waves), it is chal-lenging enough in order to show the capability of ourcode to handle a broad range of spatial frequencies inclose-to-realistic conditions.2. Then, we compute this total velocity field onto theAMR grid structure defined at the beginning of Sec.3.1. In this case, we maintain N x = 128 as the resolu-tion of the base grid, and limit the number of refine-ment patches to (cid:96) max = 4, which is still enough toshow the multi-scale capabilities of our code. Whencomputing the value of the velocities on the base levelor on AMR levels (cid:96) < (cid:96) max , we average over the uni-form grid cells enclosed in the coarser volume element.This process consistently generates a mixed, solenoidaland compressive, velocity field presenting similar scalingfeatures as the ICM over almost 3 decades in spatial fre-quency. Therefore, it can robustly show the capability ofthe code to handle multi-scale (and multi-frequency) ve-locity signals. For each test, we have validated the performance of thecode by computing a series of error statistics, which wepresent below. Let v be the input velocity field, for whichour algorithm returns its compressive and rotational com-ponents, ˜v comp and ˜v rot . Thus, the algorithm recovers atotal velocity field ˜v = ˜v comp + ˜v rot which might differ fromthe original, v , due to the numerical error in the processesinvolved, namely finite-differencing the velocity field, inte-grating the elliptical equations and finite-differencing thepotentials.In Tests 1 and 2, where only one velocity component(compressive and rotational, respectively) was present, wehave quantified the error in reconstructing these velocityfields by computing the cell-wise relative error as : The limitation is due to the fact that we need to compute thevelocity field in a uniform grid, in the first place. Thus, with N x =128 and (cid:96) max = 4, this uniform grid consists of 2048 cells. The equation below corresponds to the propagation of the vari-ance in { v cell i } to the function v cell = (cid:113)(cid:80) i =1 (cid:0) v cell i (cid:1) , assuming thevelocity components are uncorrelated. r ( v ) Median error5 95% error CI25 75% error CI
Refinement level, | v r o t |/| v | Figure 2: Results from Test 1.
Upper panel : median relative error(solid line; defined as in Eq. 14) and confidence intervals (CIs) inthe reconstruction of the velocity field.
Lower panel : fraction ofrotational velocity misreconstructed by the algorithm. ε r ( v cell ) = (cid:118)(cid:117)(cid:117)(cid:116) (cid:88) i =1 (cid:34)(cid:18) v i cell | v cell | (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) ˜ v i cell − v i cell v i cell + (cid:15) (cid:12)(cid:12)(cid:12)(cid:12)(cid:35) (14)where the subindex i runs over the three cartesian com-ponents. Note that we add a small constant, (cid:15) , in thedenominator of the relative error in v cell i to prevent theoverestimation of the error due to the cells with velocitiesclose to 0. We set (cid:15) ≡ − max ( | v i | ). For each refinementlevel, (cid:96) , we have computed the median error over all thecells, which we use as an estimate of the velocity recon-struction error. We also consider the error percentiles 5,25, 75 and 95 in order to give a confidence interval (CI)for this error.The velocity field in Test 1 (Test 2) has null rotational(compressive) component. We have checked this by com-puting the median value and the 5, 25, 75 and 95 per-centiles of | v rot | / | v | ( | v comp | / | v | ).The results for Test 1 and Test 2 are summarised in Fig-ures 2 and 3, respectively. The upper panels show ε r ( v )as defined above. Both in the constant divergence and theconstant rotational test, the median errors are small (typi-cally lower than 1 ‰ ), and even the 95-percentile errors donot exceed 1% at any refinement level. The base level, forwhich FFT is used, has much more precise results, withmedian relative errors below 10 − . For the AMR levels,we only find a very slight trend to increase the error inmore refined levels.The lower panels in Figs. 2 and 3 show | v rot | / | v | , for5 r ( v ) Median error5 95% error CI25 75% error CI
Refinement level, | v c o m p |/| v | Figure 3: Results from Test 2.
Upper panel : median relative error(solid line; defined as in Eq. 14) and confidence intervals (CIs) inthe reconstruction of the velocity field.
Lower panel : fraction ofcompressive velocity misreconstructed by the algorithm.
Test 1, and | v comp | / | v | , for Test 2. For these highly-idealised scenarios, there is virtually no cross-talk betweenthe different components: the rotational velocity in Test 1accounts for less than 10 − in more than 95% of the cellsat any level. Likewise, more than 95% of the cells in Test2 have compressive velocities less than 10 − relative to thetotal velocity magnitude.In Test 3, as both the compressive and the rotationalvelocity components are present, we have quantified therelative error in disentangling and reconstructing each ofthese by applying Eq. (14) separately to each component.We present its results in Figure 4.For the compressive velocity, the typical relative errorsare below 3 × − , while the rotational velocity presentsmedian relative errors in the range 10 − –10 − . For theAMR levels, (cid:96) ≥
1, the magnitude and behaviour of theerrors resembles what has been seen before for Tests 1 and2. Remarkably, in this example where both componentsare present, the base level does no longer exhibit a muchmore precise result, but it shows errors in the same orderas the ones in the AMR levels.Last, for the exceedingly complex Test 4, since the ac-tual decomposition is still known, we follow the same pro-cedure as in Test 3. The detailed results, presented in thesame way as for the previous test, can be found in Fig.5. Compared to Tests 1 −
3, in this case we find errorsup to an order of magnitude higher, which is not surpris-ing since in this case we have included a truly multi-scalevelocity field, with signal spanning over almost 3 decadesin spatial frequency. In any case, our code performs the r ( v c o m p ) Median error5 95% error CI25 75% error CI
Refinement level, r ( v r o t ) Figure 4: Results from Test 3. Upper ( lower ) panel presents therelative error in the reconstruction of the compressive ( rotational )velocity fields. r ( v c o m p ) Median error5 95% error CI25 75% error CI
Refinement level, r ( v r o t ) Figure 5: Results from Test 4. Upper ( lower ) panel presents therelative error in the reconstruction of the compressive ( rotational )velocity fields.
HHD decomposition with median errors around 1%, andnot exceeding 7% even at 95-percentile level, despite thefact that we have introduced oscillations close to the gridnominal frequency.The results for this last test, together with the ones wepresent below in Sec. 4.2, show the ability of our algo-rithm to perform the HHD in challenging, very non-linearvelocity fields.
4. Application to a cosmological simulation
Last, this section focuses on the results of our HHD codewhen applied to the outcomes of a cosmological simulation,which is described in the paragraphs below. As opposedto the previous highly idealised tests 1 −
3, the velocityfield in a full-cosmological simulation exhibits a plethoraof complex features due to the non-linear nature of theequations governing its evolution (see, e.g., [26, 32, 41] forclassical and recent reviews).
The simulation analysed in this paper has been car-ried out with the cosmological code
MASCLET [42], and has been already employed in a series of previous works[43, 44, 45]. Here we shall introduce the main details ofthe simulation, while some topics which are not intimatelyconnected to the analyses in this paper can be found inmore detail in the aforementioned references.
MASCLET is an Eulerian cosmological code, imple-menting high-resolution shock-capturing techniques for thedescription of the gaseous component and an N -Bodyparticle-mesh for dark matter (DM). Both components arebuilt into an AMR scheme to gain spatial and temporalresolution in the regions of interest.We have simulated a cubic domain of comoving sidelength 40 Mpc, assuming a flat ΛCDM cosmology withthe following values of the cosmological parameters: h ≡ H / (100 km s − Mpc − ) = 0 . m = 0 .
31, Ω b = 0 . Λ = 0 . n s = 0 .
96 and σ = 0 .
82, which are consistentwith the latest values reported by the Planck mission [46].The domain is discretised in a base grid of 128 cells, grant-ing a harsh resolution of ∼
310 kpc at the coarsest level.Regions with large gaseous and/or DM densities can getrecursively refined following the AMR scheme. We allow n (cid:96) = 9 refinement levels, each one halving the cell sidelength with respect to the previous level, providing a peakresolution of ∼
610 pc. The peak DM mass resolution is ∼ × M (cid:12) , equivalent to filling the domain with 1024 of such particles.The simulation started at redshift z = 100, with theinitial conditions set up by a CDM transfer function [47]and generated by a constrained realisation of the gaussianrandom field aimed to produce a massive cluster in thecenter of the computational domain by z ∼ z ∼
0, several massive clusters and groups havebeen formed. Besides gravity, the simulation accounts for abroad variety of feedback mechanisms, which are explainedin greater detail in the cited previous works.
We have run our HHD algorithm over 80 snapshots ofthe simulation described above, ranging from z = 100 to z = 0, and computed the cell-wise error as in Eq. (14)and its corresponding percentiles (5, 25, 75 and 95), asdone in the tests in Sec. 3. Figure 6 presents the overallerror statistics, defined as the median of the error statisticsover all the code outputs. In order to keep track of thedispersion of the error statistics in different snapshots, wehave also plotted the 80 individual error profiles in dottedlines, with the line colours encoding the redshift. In this case, as we do not know beforehand the ‘true’ decomposedvelocity fields to compare with the reconstructed ones, we quantifythe error by comparing the reconstructed total velocity field, ˜v = ˜v comp + ˜v rot to the input one, as defined in Eq. 14. Since, bydefinition, ˜v comp ( ˜v rot ) is the gradient of a scalar field (the rotationalof a vector field), it is irrotational (solenoidal). We have, indeed,checked that our high-order derivatives verify this, typically muchbetter than 1 ‰ . Therefore, as the decomposition is unique, checkingthat ˜v = v proves the validity of the method. Refinement level, r ( v ) Median error5 95% error CI25 75% error CI 1004.032.742.091.510.9580.5050.2340.006
Redshift, z Figure 6: Results from the tests of the HHD algorithm over
MAS-CLET outputs. The blue solid line and the contours present themedian over all the snapshots of the median relative error and theCIs defined as in the figures in Sec. 3. The dotted lines representthe individual error statistics for each snapshot, with the colourscaleencoding the redshift.
The median relative errors in describing the velocityfield as the sum of a compressive and a rotational com-ponent are typically in the order of or slightly less than1% for all refinement levels, with only a small trend toincrease the error with the refinement level. Even at the95% error percentile, the relative errors fall below 5%. Athigh redshift, the errors at low refinement levels ( (cid:96) ≤ In order to exemplify the ability of the code to split thecomponents of a highly complex velocity field, we focus ona massive galaxy cluster , with mass M vir (cid:39) . × M (cid:12) and radius R vir (cid:39) .
99 Mpc, at z (cid:39)
0. We present inFigure 7 slices of gas density (top left panel), total ve-locity (top right), compressive velocity (bottom left) andsolenoidal velocity (bottom right). The velocity mapsshow both magnitude (encoded in color) and direction inthe slice plane (arrows). This same object has been analysed in great detail in [44] (focus-ing on its observational properties) and [45] (exploring its accretionhistory).
The density map shows that, by z ∼
0, the clusteris, indeed, relatively relaxed, sitting in the gravitationalpotential dominated by the dark matter component, inan approximately spherical shape. Several density dis-continuities can be easily discerned, mainly correspond-ing to (internal) merger shocks and (external) accretionshocks. A filament penetrating quite inner radii, of around r ∼ (cid:46) r (cid:46) (cid:39) R vir . In the inner regions ( r (cid:46) . r (cid:38) .
5. Conclusions
In this paper, we have proposed a novel method toperform a Hemholtz-Hodge decomposition in AMR ve-locity fields, or virtually in any description which can besmoothed over an ad-hoc hierarchy of grids. Although ourprimary focus has been cosmological simulations of struc-8 igure 7: Maps around a massive cluster. Each map is a mass-weighted projection of the same region around the cluster, 10 Mpc on each sideand ∼
150 kpc thick, limiting the image resolution to ∼
10 kpc.
Top left:
Gas density (in units of the background density of the Universe).
Top right:
Total velocity magnitude.
Bottom left:
Compressive velocity magnitude.
Bottom right:
Solenoidal velocity magnitude. Thearrows in the velocity maps represent the projection of the corresponding velocity fields in the slice plane. .1 1 r (Mpc) v r ( k m s ) TotalCompressiveSolenoidal
Figure 8: Radial profiles of radial velocity for the total ( red ), com-pressive ( blue ) and solenoidal ( green ) velocity components. ture formation, the method is general and could be easilyextended to any type of hydrodynamical simulation.Previous works in the field of numerical cosmology typ-ically use uniform grids and work straightforwardly inFourier space. However, this procedure requires to per-form constrained simulations (or resimulations) of specificobjects of interest (e.g., a galaxy cluster) in order to beable to describe it with a uniform computational grid ata reasonable computational cost. Our algorithm, instead,can be applied to full-cosmological simulations, withoutthe need of performing resimulations and keeping the fulldescription at the maximum resolution at each position.The performance of the code has been validated in aseries of idealised tests, for which the analytical decompo-sition is known. Our algorithm has shown to succeed indisentangling the compressive and solenoidal velocity com-ponents and reconstructing the input velocity field, withtypical errors in the order of 1 ‰ or below (1% in the morecomplex, ICM-like test). Our errors seem comparable toor even better than the ones displayed by the tests in [10,Appendix A1.2].For exceedingly complex velocity fields, like the onesgenerated by actual cosmological simulations at low red-shifts, where turbulence is fully developed (e.g., [19]) andvelocity fluctuates on many different scales, our tests showthat the decomposition can be brought about with medianerrors below 1%, even resolving scales smaller than the kpcin a domain of several tens of Mpc along each direction.This procedure, whose implementation has been madepublicly available (see Sec. 2.1), will allow us to furtherexplore the dynamics of the turbulent velocity field in theICM of large samples of clusters in future works. Acknowledgements
We gratefully acknowledge the anonymous referees fortheir valuable feedback, which has helped us to improvethe quality of this manuscript. This work has been sup-ported by the Spanish Ministerio de Ciencia e Innovaci´on(MICINN, grants AYA2016-77237-C3-3-P and PID2019-107427GB-C33) and by the Generalitat Valenciana (grantPROMETEO/2019/071). Simulations have been carriedout using the supercomputer Llu´ıs Vives at the Serveid’Inform`atica of the Universitat de Val`encia. This re-search has made use of the following open-source packages:
NumPy [49],
SciPy [50], matplotlib [51] and yt [52]. References [1] J. A. ZuHone, A Parameter Space Exploration of Galaxy Clus-ter Mergers. I. Gas Mixing and the Generation of Cluster En-tropy, ApJ 728 (1) (2011) 54. arXiv:1004.3820 , doi:10.1088/0004-637X/728/1/54 .[2] M. Gaspari, E. Churazov, Constraining turbulence and con-duction in the hot ICM through density perturbations, A&A559 (2013) A78. arXiv:1307.4397 , doi:10.1051/0004-6361/201322295 .[3] I. Zhuravleva, E. Churazov, A. A. Schekochihin, S. W. Allen,P. Ar´evalo, A. C. Fabian, W. R. Forman, J. S. Sanders,A. Simionescu, R. Sunyaev, A. Vikhlinin, N. Werner, Turbulentheating in galaxy clusters brightest in X-rays, Nature 515 (7525)(2014) 85–87. arXiv:1410.6485 , doi:10.1038/nature13830 .[4] F. Hofmann, J. S. Sanders, K. Nandra, N. Clerc, M. Gaspari,Thermodynamic perturbations in the X-ray halo of 33 clustersof galaxies observed with Chandra ACIS, A&A 585 (2016) A130. arXiv:1510.08445 , doi:10.1051/0004-6361/201526925 .[5] D. Eckert, M. Gaspari, F. Vazza, F. Gastaldello, A. Tra-macere, S. Zimmer, S. Ettori, S. Paltani, On the Connec-tion between Turbulent Motions and Particle Acceleration inGalaxy Clusters, ApJL 843 (2) (2017) L29. arXiv:1705.02341 , doi:10.3847/2041-8213/aa7c1a .[6] A. Bonafede, M. Br¨uggen, D. Rafferty, I. Zhuravleva, C. J.Riseley, R. J. van Weeren, J. S. Farnes, F. Vazza, F. Savini,A. Wilber, A. Botteon, G. Brunetti, R. Cassano, C. Ferrari,F. de Gasperin, E. Orr´u, R. F. Pizzo, H. J. A. R¨ottgering,T. W. Shimwell, LOFAR discoveryof radio emission in MACSJ0717.5+3745, MNRAS 478 (3) (2018) 2927–2938. arXiv:1805.00473 , doi:10.1093/mnras/sty1121 .[7] M. Roncarelli, M. Gaspari, S. Ettori, V. Biffi, F. Brighenti,E. Bulbul, N. Clerc, E. Cucchetti, E. Pointecouteau, E. Rasia,Measuring turbulence and gas motions in galaxy clusters viasynthetic Athena X-IFU observations, A&A 618 (2018) A39. arXiv:1805.02577 , doi:10.1051/0004-6361/201833371 .[8] F. Miniati, The Matryoshka Run. II. Time-dependent Turbu-lence Statistics, Stochastic Particle Acceleration, and Micro-physics Impact in a Massive Galaxy Cluster, ApJ 800 (1) (2015)60. arXiv:1409.3576 , doi:10.1088/0004-637X/800/1/60 .[9] W. Schmidt, J. F. Engels, J. C. Niemeyer, A. S. Almgren, Hotand turbulent gas in clusters, MNRAS 459 (1) (2016) 701–719. arXiv:1603.04711 , doi:10.1093/mnras/stw632 .[10] F. Vazza, T. W. Jones, M. Br¨uggen, G. Brunetti, C. Gheller,D. Porter, D. Ryu, Turbulence and vorticity in Galaxy clustersgenerated by structure formation, MNRAS 464 (1) (2017) 210–230. arXiv:1609.03558 , doi:10.1093/mnras/stw2351 .[11] L. Iapichino, C. Federrath, R. S. Klessen, Adaptive mesh re-finement simulations of a galaxy cluster merger - I. Resolv-ing and modelling the turbulent flow in the cluster outskirts,MNRAS 469 (3) (2017) 3641–3655. arXiv:1704.02922 , doi:10.1093/mnras/stx882 .
12] F. Vazza, M. Br¨uggen, C. Gheller, P. Wang, On the amplifi-cation of magnetic fields in cosmic filaments and galaxy clus-ters, MNRAS 445 (4) (2014) 3706–3722. arXiv:1409.2640 , doi:10.1093/mnras/stu1896 .[13] D. H. Porter, T. W. Jones, D. Ryu, Vorticity, Shocks, andMagnetic Fields in Subsonic, ICM-like Turbulence, ApJ 810 (2)(2015) 93. arXiv:1507.08737 , doi:10.1088/0004-637X/810/2/93 .[14] C. Federrath, G. Chabrier, J. Schober, R. Banerjee, R. S.Klessen, D. R. G. Schleicher, Mach Number Dependence of Tur-bulent Magnetic Field Amplification: Solenoidal versus Com-pressive Flows, PRL 107 (11) (2011) 114504. arXiv:1109.1760 , doi:10.1103/PhysRevLett.107.114504 .[15] C. Federrath, R. S. Klessen, W. Schmidt, The Fractal DensityStructure in Supersonic Isothermal Turbulence: Solenoidal Ver-sus Compressive Energy Injection, ApJ 692 (1) (2009) 364–374. arXiv:0710.1359 , doi:10.1088/0004-637X/692/1/364 .[16] L. Iapichino, W. Schmidt, J. C. Niemeyer, J. Merklein, Turbu-lence production and turbulent pressure support in the inter-galactic medium, MNRAS 414 (3) (2011) 2297–2308. arXiv:1102.3352 , doi:10.1111/j.1365-2966.2011.18550.x .[17] C. Federrath, J. Roman-Duval, R. S. Klessen, W. Schmidt,M. M. Mac Low, Comparing the statistics of interstellar tur-bulence in simulations and observations. Solenoidal versus com-pressive turbulence forcing, A&A 512 (2010) A81. arXiv:0905.1060 , doi:10.1051/0004-6361/200912437 .[18] C. Federrath, On the universality of supersonic turbulence,MNRAS 436 (2) (2013) 1245–1257. arXiv:1306.3989 , doi:10.1093/mnras/stt1644 .[19] F. Miniati, The Matryoshka Run: A Eulerian Refinement Strat-egy to Study the Statistics of Turbulence in Virialized Cos-mic Structures, ApJ 782 (1) (2014) 21. arXiv:1310.2951 , doi:10.1088/0004-637X/782/1/21 .[20] J. A. Acosta-Pulido, I. Agudo, A. Alberdi, J. Alcolea, E. J. Al-faro, A. Alonso-Herrero, G. Anglada, P. Arnalte-Mur, Y. As-casibar, B. Ascaso, R. Azulay, R. Bachiller, A. Baez-Rubio,E. Battaner, J. Blasco, C. B. Brook, V. Bujarrabal, G. Busquet,M. D. Caballero-Garcia, C. Carrasco-Gonzalez, J. Casares, A. J.Castro-Tirado, L. Colina, F. Colomer, I. de Gregorio-Monsalvo,A. del Olmo, J. F. Desmurs, J. M. Diego, R. Dominguez-Tenreiro, R. Estalella, A. Fernandez-Soto, E. Florido, J. Font,J. A. Font, A. Fuente, R. Garcia-Benito, S. Garcia-Burillo,B. Garcia-Lorenzo, A. Gil de Paz, J. M. Girart, J. R.Goicoechea, J. F. Gomez, M. Gonzalez-Garcia, O. Gonzalez-Martin, J. I. Gonzalez-Serrano, J. Gorgas, J. Gorosabel, A. Gui-jarro, J. C. Guirado, L. Hernandez-Garcia, C. Hernandez-Monteagudo, D. Herranz, R. Herrero-Illana, Y. D. Hu, N. Hue-lamo, M. Huertas-Company, J. Iglesias-Paramo, S. Jeong,I. Jimenez-Serra, J. H. Knapen, R. A. Lineros, U. Lisenfeld,J. M. Marcaide, I. Marquez, J. Marti, J. M. Marti, I. Marti-Vidal, E. Martinez-Gonzalez, J. Martin-Pintado, J. Masegosa,J. M. Mayen-Gijon, M. Mezcua, S. Migliari, P. Mimica,J. Moldon, O. Morata, I. Negueruela, S. R. Oates, M. Oso-rio, A. Palau, J. M. Paredes, J. Perea, P. G. Perez-Gonzalez,E. Perez-Montero, M. A. Perez-Torres, M. Perucho, S. Planelles,J. A. Pons, A. Prieto, V. Quilis, P. Ramirez-Moreta, C. RamosAlmeida, N. Rea, M. Ribo, M. J. Rioja, J. M. RodriguezEspinosa, E. Ros, J. A. Rubi˜no-Martin, B. Ruiz-Granados,J. Sabater, S. Sanchez, C. Sanchez-Contreras, A. Sanchez-Monge, R. Sanchez-Ramirez, A. M. Sintes, J. M. Solanes,C. F. Sopuerta, M. Tafalla, J. C. Tello, B. Tercero, M. C.Toribio, J. M. Torrelles, M. A. P. Torres, A. Usero, L. Verdes-Montenegro, A. Vidal-Garcia, P. Vielva, J. Vilchez, B. B.Zhang, The Spanish Square Kilometre Array White Book, arXive-prints (Jun. 2015). arXiv:1506.03474 .[21] Y. Fujita, M. Takizawa, C. L. Sarazin, Nonthermal Emissionsfrom Particles Accelerated by Turbulence in Clusters of Galax-ies, ApJ 584 (1) (2003) 190–202. arXiv:astro-ph/0210320 , doi:10.1086/345599 .[22] R. Cassano, G. Brunetti, Cluster mergers and non-thermalphenomena: a statistical magneto-turbulent model, MNRAS 357 (4) (2005) 1313–1329. arXiv:astro-ph/0412475 , doi:10.1111/j.1365-2966.2005.08747.x .[23] G. Brunetti, A. Lazarian, Particle reacceleration by compress-ible turbulence in galaxy clusters: effects of a reduced meanfree path, MNRAS 412 (2) (2011) 817–824. arXiv:1011.1198 , doi:10.1111/j.1365-2966.2010.17937.x .[24] A. Pinzke, S. P. Oh, C. Pfrommer, Turbulence and particleacceleration in giant radio haloes: the origin of seed electrons,MNRAS 465 (4) (2017) 4800–4816. arXiv:1611.07533 , doi:10.1093/mnras/stw3024 .[25] G. Brunetti, F. Vazza, Second-order Fermi ReaccelerationMechanisms and Large-Scale Synchrotron Radio Emission inIntracluster Bridges, PRL 124 (5) (2020) 051101. arXiv:2001.07718 , doi:10.1103/PhysRevLett.124.051101 .[26] E. Bertschinger, Simulations of Structure Formation in theUniverse, ARA&A 36 (1998) 599–654. doi:10.1146/annurev.astro.36.1.599 .[27] T. Naab, J. P. Ostriker, Theoretical Challenges in Galaxy For-mation, ARA&A 55 (1) (2017) 59–109. arXiv:1612.06891 , doi:10.1146/annurev-astro-081913-040019 .[28] F. Vazza, E. Roediger, M. Br¨uggen, Turbulence in the ICMfrom mergers, cool-core sloshing, and jets: results from a newmulti-scale filtering approach, A&A 544 (2012) A103. arXiv:1202.5882 , doi:10.1051/0004-6361/201118688 .[29] R. Mohapatra, C. Federrath, P. Sharma, Turbulence in strat-ified atmospheres: implications for the intracluster medium,MNRAS 493 (4) (2020) 5838–5853. arXiv:2001.06494 , doi:10.1093/mnras/staa711 .[30] R. Mohapatra, C. Federrath, P. Sharma, Turbulent densityand pressure fluctuations in the stratified intracluster medium,MNRAS 500 (4) (2021) 5072–5087. arXiv:2010.12602 , doi:10.1093/mnras/staa3564 .[31] A. M. Bykov, K. Dolag, F. Durret, Cosmological Shock Waves,Space Sci. Rev. 134 (1-4) (2008) 119–140. arXiv:0801.0995 , doi:10.1007/s11214-008-9312-9 .[32] K. Dolag, S. Borgani, S. Schindler, A. Diaferio, A. M. Bykov,Simulation Techniques for Cosmological Simulations, Space Sci.Rev. 134 (1-4) (2008) 229–268. arXiv:0801.1023 , doi:10.1007/s11214-008-9316-5 .[33] A. G. Kritsuk, ˚A. Nordlund, D. Collins, P. Padoan, M. L. Nor-man, T. Abel, R. Banerjee, C. Federrath, M. Flock, D. Lee, P. S.Li, W.-C. M¨uller, R. Teyssier, S. D. Ustyugov, C. Vogel, H. Xu,Comparing Numerical Methods for Isothermal Magnetized Su-personic Turbulence, ApJ 737 (1) (2011) 13. arXiv:1103.5525 , doi:10.1088/0004-637X/737/1/13 .[34] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery,Numerical recipes in FORTRAN. The art of scientific comput-ing, 1992.[35] G. Arfken, H. Weber, F. Harris, Mathematical Methods forPhysicists: A Comprehensive Guide, Elsevier Science, 2013.[36] G. Kowal, A. Lazarian, Velocity Field of Compressible Mag-netohydrodynamic Turbulence: Wavelet Decomposition andMode Scalings, ApJ 720 (1) (2010) 742–756. arXiv:1003.3697 , doi:10.1088/0004-637X/720/1/742 .[37] S. K. Harouna, V. Perrier, Helmholtz-Hodge Decomposition on[0,1]d by Divergence-Free and Curl-Free Wavelets, in: J.-D.Boissonnat, P. Chenin, A. Cohen, C. Gout, T. Lyche, M.-L.Mazure, L. Schumaker (Eds.), Curves and Surfaces, SpringerBerlin Heidelberg, Berlin, Heidelberg, 2012, pp. 311–329.[38] R. W. Hockney, J. W. Eastwood, Computer Simulation UsingParticles, 1981.[39] J. Burgers, Mathematical examples illustratingrelations occur-ring in the theory of turbulentfluid motion. verh. k, Akad. Wet.Amst. Addeel. Nat 17 (1) (1939).[40] A. Kolmogorov, The Local Structure of Turbulence in Incom-pressible Viscous Fluid for Very Large Reynolds’ Numbers,Akademiia Nauk SSSR Doklady 30 (1941) 301–305.[41] S. Planelles, D. R. G. Schleicher, A. M. Bykov, Large-ScaleStructure Formation: From the First Non-linear Objects toMassive Galaxy Clusters, Space Sci. Rev. 188 (1-4) (2015) 93–139. arXiv:1404.3956 , doi:10.1007/s11214-014-0045-7 .
42] V. Quilis, A new multidimensional adaptive mesh refinementhydro + gravity cosmological code, MNRAS 352 (2004) 1426–1438. arXiv:astro-ph/0405389 , doi:10.1111/j.1365-2966.2004.08040.x .[43] V. Quilis, S. Planelles, E. Ricciardelli, Is ram-pressure strip-ping an efficient mechanism to remove gas in galaxies?, MNRAS469 (1) (2017) 80–94. arXiv:1703.09446 , doi:10.1093/mnras/stx770 .[44] S. Planelles, P. Mimica, V. Quilis, C. Cuesta-Mart´ınez, Mul-tiwavelength mock observations of the WHIM in a simulatedgalaxy cluster, MNRAS 476 (4) (2018) 4629–4648. arXiv:1802.09458 , doi:10.1093/mnras/sty527 .[45] D. Vall´es-P´erez, S. Planelles, V. Quilis, On the accretion historyof galaxy clusters: temporal and spatial distribution, MNRAS499 (2) (2020) 2303–2318. arXiv:2009.13882 , doi:10.1093/mnras/staa3035 .[46] Planck Collaboration, Planck 2018 results. VI. Cosmologicalparameters, A&A 641 (2020) A6. arXiv:1807.06209 , doi:10.1051/0004-6361/201833910 .[47] D. J. Eisenstein, W. Hu, Baryonic Features in the Matter Trans-fer Function, ApJ 496 (2) (1998) 605–614. arXiv:astro-ph/9709112 , doi:10.1086/305424 .[48] Y. Hoffman, E. Ribak, Constrained realizations of Gaussianfields - A simple algorithm, ApJ 380 (1991) L5–L8. doi:10.1086/186160 .[49] T. E. Oliphant, A guide to NumPy, Vol. 1, Trelgol PublishingUSA, 2006.[50] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haber-land, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson,W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wil-son, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones,R. Kern, E. Larson, C. Carey, ˙I. Polat, Y. Feng, E. W. Moore,J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Hen-riksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H.Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1. 0 Contribu-tors, SciPy 1.0: Fundamental Algorithms for Scientific Comput-ing in Python, Nature Methods 17 (2020) 261–272. doi:https://doi.org/10.1038/s41592-019-0686-2 .[51] J. D. Hunter, Matplotlib: A 2d graphics environment, Com-puting in Science & Engineering 9 (2007) 90–95. doi:10.1109/MCSE.2007.55 .[52] M. J. Turk, B. D. Smith, J. S. Oishi, S. Skory, S. W. Skillman,T. Abel, M. L. Norman, yt: A Multi-code Analysis Toolkit forAstrophysical Simulation Data, ApJS 192 (1) (2011) 9. arXiv:1011.3514 , doi:10.1088/0067-0049/192/1/9 ..