We introduce adaptive particle refinement for compressible smoothed particle hydrodynamics SPH. SPH calculations have the natural advantage that resolution follows mass but this is not always optimal. Our implementation allows the user to specify local regions of the simulation that can be more highly resolved. We test our implementation on practical applications including a circumbinary disc a planet embedded in a disc and a flyby. By comparing with equivalent globally high resolution calculations we show that our method is accurate and fast with errors in the mass accreted onto sinks of less than 9 percent and speed ups of 1.076.62 times for the examples shown. Our method is adaptable and easily extendable for example with multiple refinement regions or derefinement.
['Astrophysics - Instrumentation and Methods for Astrophysics']
Adaptive particle refinement for compressible smoothed particle hydrodynamics
{'No Header': '<!-- image -->', 'Adaptive particle refinement for compressible smoothed particle hydrodynamics': 'Rebecca Nealon 1,2 and Daniel J. Price 3  \n- 1 Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Road, CV4 7AL Coventry, UK\n- 2 Department of Physics, University of Warwick, Gibbet Hill Road, CV4 7AL Coventry, UK\n- 3 School of Physics and Astronomy, Monash University, Clayton, Vic 3800, Australia\n- Author for correspondence: R. Nealon, Email: rebecca.nealon@warwick.ac.uk.', 'Abstract': "We introduce adaptive particle refinement for compressible smoothed particle hydrodynamics (SPH). SPH calculations have the natural advantage that resolution follows mass, but this is not always optimal. Our implementation allows the user to specify local regions of the simulation that can be more highly resolved. We test our implementation on practical applications including a circumbinary disc, a planet embedded in a disc and a flyby. By comparing with equivalent globally high resolution calculations we show that our method is accurate and fast, with errors in the mass accreted onto sinks of less than 9% and speed ups of 1.07 - 6.62 × for the examples shown. Our method is adaptable and easily extendable, for example with multiple refinement regions or derefinement.  \nKeywords: hydrodynamics, methods: numerical  \nThis paper considers the best way to locally adapt resolution in the simulation method known as smoothed particle hydrodynamics (SPH, Lucy 1977; Gingold and Monaghan 1977). Because resolution follows mass in SPH, the best resolved region often corresponds to the region of interest. However, this may be inefficient, as high densities correlate with short timesteps which are expensive. A potential solution is Adaptive Particle Refinement (APR): splitting and merging of particles. This method allows multiple resolutions to co-exist in the same simulation (Monaghan and Varnas 1988; Meglicki, Wickramasinghe, and Bicknell 1993; Kitsionas and Whitworth 2002, 2007; Vacondio et al. 2013; López, Roose, and Recarey Morfa 2013; Barcarolo et al. 2014; Chiron et al. 2018; Gao, Qiu, and Fu 2022). APR is conceptually similar to Adaptive Mesh Refinement (AMR, Berger and Colella 1989; Truelove et al. 1997) commonly used in simulations employing Eulerian meshes.  \nAlthough less widely exploited than mesh refinement, the fundamentals of particle refinement and its applications are well developed. Monaghan and Varnas (1988) and Meglicki, Wickramasinghe, and Bicknell (1993) first applied particle splitting in SPH. They merged particles in their simulations to create fewer, more massive particles when the density in a cell exceeded a critical density threshold. Meglicki, Wickramasinghe, and Bicknell (1993) additionally included particle splitting in low density regions to increase the local resolution. In the context of cosmological simulations, based on the method outlined in Porter (1985), Katz and White (1993) (and later Navarro and White 1994) scaled the mass of the particles in spherical, nested layers throughout the computational domain such that the innermost layer had the smallest mass and thus highest resolution. Thacker and Couchman (2000) found that there was a small amount of noise generated at the boundary between the mass layers but that it did not affect the structures identified  \nin their simulations. Bertschinger (2001) produced a method to accurately initialise simulations with multiple masses for a Gaussian random field. Klypin et al. (2001) demonstrated numerical convergence with this simulation style. This 'multiple mass' or 'zoom' method, the approach of nested regions is a core feature of cosmological simulations (e.g. Springel 2005).  \nBørve, Omang, and Trulsen (2001) implemented an alternative multi-resolution scheme they called Regularised-SPH (RSPH). In this method, the smoothing length h was set to a piece-wise constant in steps of 2 with the contribution from neighbouring particles interpolated using a grid. In its most modern invocation (Børve, Omang, and Trulsen 2005), this method used auxiliary particles that exist at the boundary between different h regions and which are passively evolved. Børve, Omang, and Trulsen (2006) tested RSPH on a multidimensional MHD shock, showing excellent shock capturing properties. Børve, Speith, and Trulsen (2009) also successfully applied RSPH to planet-disc interactions in two dimensions.  \nKitsionas and Whitworth (2002) considered the collapse of clumps in self gravitating filaments. They increased the resolution in high density regions to ensure that the Jean's criteria was met, finding their method was robust and comparable to results from AMR calculations. Importantly they used 13 children for each parent ( n child = 13 where we refer to more massive and less massive particles as 'parents' and 'children', respectively) to ensure spherical symmetry in 3D when splitting particles and established an empirical test to determine the radius of the sphere that the children are located on ( r sep). Kitsionas and Whitworth (2007) showed that using a mass weighted kernel (with 50 times the mass of the largest particle) was preferable to volume weighted.  \nFeldman and Bonet (2007) introduced a general splitting procedure, placing child particles symmetrically in hexagons or  \ntriangles around their parent to conserve angular momentum. Importantly, they quantified the error that is generated when a parent particle is split into children; when the children have equal mass this error is due to the new particle arrangement and is a function of the distance between the children and their parent.  \nLastiwka, Quinlan, and Basa (2005) showed that this error could be reduced by slightly shifting the children. López, Roose, and Recarey Morfa (2013) extended this idea by measuring the error between the parent and child distributions as a function of the separation between the children when they are newly placed, r sep. Solutions to this error include shuffling the particles until the error is minimised (Vacondio et al. 2013), blending the region between the parent and children regions (Barcarolo et al. 2014), particle disordering (Chiron et al. 2018) and particle regularisation (Gao, Qiu, and Fu 2022). Vacondio et al. (2013) first introduced particle merging in incompressible applications but found it too expensive for practical use. Chiaki and Yoshida (2015) showed that the noise introduced by a split could be reduced by placing the children particles on the vertices of a Voronoi mesh.  \nMost recently, Franchini, Lupi, and Sesana (2022) used particle splitting in a meshless-finite-volume method (Hopkins 2014) to better resolve circumstellar discs within the inner cavity of a circumbinary disc. Adopting the method from Anglés-Alcázar et al. (2021) they used n child = 2 and r sep = min(0.25 × h parent, 0.35 × d neighbour ) (where d neighbour is the distance to the nearest neighbour) to ensure that split particles were not inadvertently placed too close to existing particles (i.e. that the fluid elements did not overlap). The robustness of their method was demonstrated in Duffell et al. (2024) where their splitting method found agreement with comparable codes including P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc (Price, Wurster, et al. 2018).  \nWhile the above implementations are promising, practical limitations remain. First, while most authors state that merging is possible, in practice it is not used - and never in simulations of compressible flow. On the rare occasion that merging is included it is found to be computationally prohibitive (Vacondio et al. 2013) which restricts what problems APR can be practically applied to. Second, typical implementations require hard-wiring of parameters like r sep in order to reduce noise when particles are split, but these values are derived from simple tests (e.g. Kitsionas and Whitworth 2002). Third, some of the above applications in astrophysics allow particles of different refinement levels to mix which may result in numerical instabilities (e.g. Chiron et al. 2018).  \nIn this paper we introduce an APR implementation into the smoothed particle hydrodynamics code P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc that includes both splitting and merging, separation of refinement levels and which is accurate and fast. In Section 1 we outline the method. Section 2 applies our method to examples and practical calculations compared to equivalent global high resolution simulations. We measure speed, accuracy and disc storage. In Sections 3 and 4 we discuss and draw conclusions. Basic tests for the interested reader are summarised in Appendix 1.", '1. Method': 'Here we describe the core of our implementation including splitting, merging, relaxing and the order in which these are completed. The splitting and merging processes are summarised in Figure 1.', '1.1 Overall procedure': "We assume spherical refinement and de-refinement zones. Their size and location may be either fixed or dynamic and comoving with another particle. In some of our later applications the refinement zone is centred on a moving point mass particle but it may also be set by a particle property like density. Here we show examples with spherical zones but other volumes are possible. The main points of our implementation can be captured with five 'rules':  \n- 1. We assign all particles a refinement level, ℓ , determined exclusively from the spatial position of the particle. The refinement level represents the number of refinements above the base resolution.\n- 2. When a particle enters a new refinement zone with a given refinement level\n- · it is split if the particle's current refinement level is less than that of the zone.\n- · it is merged if the particle's current refinement level is greater than that of the zone.\n- 3. In simulations where the goal is to locally increase the resolution the default refinement level is set to ℓ = 0. In simulations where we want to locally decrease the resolution the default refinement level is instead set to ℓ = ℓ max.\n- 4. We restrict the difference in mass between adjacent refinement levels to be strictly a factor of two.\n- 5. As P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc stores the mass by particle type, the refinement level is also used to relate the refined mass mp ( ℓ ) from the largest particle mass mp (0) with the following relation:  \nmp ( ℓ ) = mp (0) 2 ℓ , (1)  \nwhere ℓ >= 1. When we compare particles with different refinement levels we refer to the more highly resolved particles as 'children' and the lower resolution particles as 'parents'. In our scheme then, two children merge to form one parent and on a split one parent produces two children.", '1.2 Splitting': "A parent particle with refinement level ℓ is split into two children when it crosses the boundary into a refinement zone with a higher refinement level. Once a particle is confirmed to be split;  \n- 1. A new child particle is made with identical properties to the parent.\n- 2. The parent itself is reassigned to be a second child particle.  \nFigure 1. Schematic showing our refining and derefining process. The flow of the fluid is le/f\\_t to right and parent particles are larger and blue, children particles smaller and pink and the particle size is proportional to its mass. Particles that are split or merged in the time step shown are indicated with outlines and ℓ shows the refinement level. Le/f\\_t: As a parent particle crosses the boundary it is split into two children particles, their common centre of mass is at the parent's location and they are split tangentially to the boundary. Right: Children particles are paired according to our modified k -d tree grouping (indicated with the grey boxes). When the centre of mass of a cell crosses the boundary the particles are merged, with the parent adopting the average velocity and position of the children. Further details are in Section 1.  \n<!-- image -->  \n- 3. Both children have their refinement level increased and their smoothing lengths rescaled (assuming constant density)  \nh ℓ +1 = ( 1 2 ) 1/3 h ℓ . (2)  \n- 4. Children are placed on opposite sides of the parent at a distance r sep = 0.2 h ℓ .\n- 5. Children are additionally placed tangential to the boundary of the refinement region (as in Figure 1).  \nThe scaling value in r sep has been determined empirically to reduce noise following the test from Kitsionas and Whitworth (2002), however our inclusion of relaxing (see Section 1.4) means that our implementation is not particularly sensitive to this choice (as demonstrated in Appendix 1.3). Additionally, as P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc calculates the density and smoothing lengths iteratively (2.1.4, Price, Wurster, et al. 2018), the smoothing lengths of surrounding particles are accurately adjusted to accommodate the new particles straight after the split (and later, merge) has occurred.  \nTo ensure particles are split tangentially to the boundary in 3D we identify the vector ⃗ v between the original particle and the centre of the refinement zone. A vector ⃗ w represents the perpendicular bisector of ⃗ v and is always tangential to the boundary of the refinement zone. We rotate ⃗ w around ⃗ v through a random angle and place the children on opposite sides of the parent along ± ⃗ w . As in Franchini, Lupi, and Sesana (2022), if the distance shifted is more than 0.35 × the distance to the closest neighbour we adopt this value instead, however we have found that this only affects a handful of particles even in very well resolved simulations.  \nFigure 2. Flowchart summarising the APR routine implemented in P/H.s/A.s/N.s/T.s/O.s/M.s.  \n<!-- image -->", '1.3 Merging': "Wemake use of a modified k -d tree in P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc to make rapid on-the-fly merging possible. We have altered the existing k -d tree so that it always returns leaf nodes with exactly two particles (this modification is only used in the call from the APR routine and is not used elsewhere in the calculations). Every time the domain is split through the centre of mass we enforce an even number of particles on either side of the new 'subdomains'. When a split occurs that results in an odd number of particles in each new subdomain, we move the location of the split across one particle in the direction that makes the most balanced split by number of particles. The tree continues until finally all leaves have precisely two particles.  \nAt each time-step where merging occurs;  \n- 1. All particles at the same refinement level ℓ are isolated and then paired them using our modified k -d tree.\n- 2. For each leaf in this tree, we calculate what refinement level the leaf has based on it's centre of mass and compare it to the refinement level of the child particles, ℓ .\n- 3. The particles are merged when the refinement level based on the position of the leaf is less than the refinement level of the child particles in the leaf.\n- 4. One of these children particle is removed from the simulation and the remaining child particle is reassigned as a parent.\n- 5. The parent is given the position and velocity of the centre of mass of the leaf (i.e. the two children), the refinement level is decreased and the smoothing length is scaled by  \nh ℓ -1 = 2 1/3 h ℓ . (3)  \n- 6. If an odd number of particles are considered for merging, we discard the last particle that is listed at that timestep. In the next timestep when particles are again considered for merging this particle is then included.  \nThe right panel of Figure 2 schematically shows the merging process. With our implementation it is important to note that although the boundary for where merging occurs is well defined, because we are merging based on the children particle distribution the practical boundary is fuzzy and a transition region typically about 10% of the larger smoothing length naturally forms. In well resolved simulations (e.g. Section 2) this region is negligible in width. A major benefit of our approach is we have at most a maximum of one spare particle per refinement level per timestep. As we do not consider ℓ = 0 particles for merging, the maximum error e in the total mass of these spare particles for ℓ max refinement levels is  \ne = mp (0) ℓ max ∑ ℓ =1 1 2 ℓ , (4)  \nwhich is conveniently always less than mp (0), the largest particle mass.", '1.4 Relaxing': "When splitting or merging a significant number of particles - for example the first time the splitting/merging routine is called - there is a discrepancy between the original and refined density distribution. This effect was identified by Feldman and Bonet (2007) and is due to how well the new particle arrangement can represent the original density distribution.  \nTo combat this we have implemented a relaxing procedure which shuffles the introduced particles until they more accurately represent the original distribution. When splitting or merging an original set of particles into a new set, our algorithm  \n- 1. Calculates the accelerations of the new set of particles at their current locations, a new.\n- 3. Shifts each particle by ∆ x calculated from\n- 2. Calculates the accelerations at the locations of the new set of particles, interpolated from the original reference particles, a ref .  \n∆ x = 0.5 ∆ t 2 ( a new a ref ), (5)  \nwhere ∆ t is calculated from the sound speed and smoothing length by ∆ t = 0.3 h / cs on the new particle.  \nFor any given particle the magnitude of the shuffle is also capped to be less than the particles smoothing length. At each shuffle step, we estimate a 'kinetic energy' defined as the magnitude of the shift divided by the particle's timestep squared and summed across all the shuffled particles. For multiple shuffles we find that the kinetic energy decreases in an exponential fashion with each subsequent shuffle smaller than the previous. Our shuffling process is repeated until either the total kinetic energy for a shuffle has decreased to 0.5% of its original value or 50 shuffles have been completed. These limits are set to strike a balance between accuracy and computational expense, as we have found that more shuffles does not continue to dramatically improve the particle distribution. How quickly the kinetic energy limit is reached does depend on the resolution, with higher resolution simulations achieving it in fewer shuffles.", '1.5 With individual timesteps': "A key part of our implementation is compatibility with individual particle timesteps as these are used widely in P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc applications (see Price, Wurster, et al. 2018). Each particle is assigned to a timestep bin; these are arranged such that the zeroth bin has ∆ t = ∆ t max (which is limited to be the time between outputs) and particles are arranged in bins according to their local timestep constraint, where ∆ t decreases by factor of two in each bin.  \nEach bin is then evolved separately, synchronising with the bin above when the appropriate timesteps are synchronised. Particles are then defined as 'active' (as in, they will be moved by the stepping routine) when their bin is being evolved.  \nWhen implemented in conjunction with APR we restrict the splitting and merging procedure to occur over active particles only. This does mean that a merging particle's closest  \nTable 1. Summary of the APR simulations shown in Section 2. Columns state the Name of the simulation, the radius of the central refinement region r ℓ , the steps into each refinement zone dr ℓ and the average number of particles used N . The error, speed up and storage are all compared to the high resolution reference cases. The method to measure the error is described in the text for each simulation and is measured according to Equation 6. The storage is calculated as a fraction compared to the high resolution reference calculation. All simulations have 3 levels of increased refinement except F6 which has ℓ = 6 . Here 'B' refers to binary, 'P' to planet-disc and 'F' to flyby simulations.  \nneighbour could be on a different time-step, forcing the particle to merge with an active particle further away. In practice this is rare because the particles have been paired using our modified k -d tree. This pairing takes advantage of the fact that spatially associated particles tend to have the same acceleration and so are naturally on the same timestep bin such that the closest neighbour is also on the same timestep. Additionally, when a particle is either split, merged or a new particle introduced through either of those procedures it is conservatively assigned the shortest timestep (corresponding to the maximum timestep bin).", '1.6 Every timestep': 'Figure 2 summarises the process undertaken at each timestep on active particles. First the location, size and number of refinement regions are updated. If relaxing is required at this step the properties of all particles are also saved as a reference.  \nWethen consider whether any particles should be split. Using the particles current position we calculate what refinement level it should have and compare this to its current refinement level. When the particle has a refinement level less than what it should have we split it according to Section 1.2. We repeat this process across for all active particles and for each refinement level, starting from the lowest refinement level.  \nAfter all the particles have been split we then consider whether any should be merged. We isolate all particles that are at a particular refinement level and put these particles into our modified tree routine, merging as per Section 1.3. Mirroring the splitting procedure, this process is repeated from the highest to lowest refinement level. This order means that particles are only ever changed one refinement level each split/merge which we have found leads to a smoother distribution (e.g.  \nAppendix 1.4). Only once all the particles have been split and merged do we consider the relaxing process, using the original particle distribution as the reference set. With the particles refinement levels updated the program then returns to the normal time stepping routine and continues.', '2. Example applications': "We demonstrate the accuracy, speed up and typical use cases of our APR implementation with three example applications: a circumbinary disc, planet disc interactions and a stellar flyby. In each case we include a low and high resolution comparison simulation, where the high resolution simulation has the same global resolution as the highest APR zone. Our examples also demonstrate the adaptability of the implementation and provide typical values for the size of the refinement region and the step width. Our examples include gravity from sink particles (stars and planets) as well as accretion onto those sinks. We also assume a vertically isothermal temperature profile for these applications, but we refer to Appendix 1 for tests with an adiabatic equation of state. Radiative cooling or source terms that are sensitive to temperature are left to future work. Our results are summarised in Table 1. In all cases we find that our implementation is accurate, fast and requires less storage space when compared to the high resolution reference cases.  \nIn the examples shown here we choose to use n child = 2 and nested refinement levels and with 'relaxing' only implemented in the first step (see Appendix 1). We additionally ensure that the velocity of the fluid across the boundary is small (see Appendix 1.5) by choosing the location of our refinement regions carefully. These examples are presented with a Wendland C2 kernel with a kernel radius of 2 h and h fact = 1.3 (Wendland 1995; Price, Wurster, et al. 2018). As we are considering parti-  \nFigure 3. Columndensity of the HD142527 simulations from Price, Cuello, et al. (2018) at t = 1137 years. The le/f\\_t and right panels have no APR and are separated by a factor of two in spatial resolution. The middle panel shows simulation B4 with three levels of refinement, locally matching the evolution of the high resolution reference case. The refinement zone is shown with the light green circles where the highest resolution region is indicated with a solid line and the nested regions with dotted lines. The similarities in the streamers and circumprimary disc confirm our APR implementation is capable of accurately locally increasing the resolution of a simulation. A movie of these simulations is available online.  \n<!-- image -->  \nFigure 4. Mass of each sink in the HD142527 simulations with di/fferent combinations of refinement region sizes. The mass accretion rate reflects the properties of the disc surrounding the sink, confirming the structure of circumprimary disc in the APR simulations is the same as the high resolution reference calculation (as seen in Figure 3).  \n<!-- image -->  \ncles of different masses co-existing in the same simulation it is prudent to consider kernels other than the typical M 4 cubic spline (e.g. Dehnen and Aly 2012). Our comparison of M 4 , Wendland C2, Wendland C4 and Wendland C6 is shown in Appendix 1. Here we find that the M 4 cubic spline is robust to pairing in these tests but that the higher order splines offer a slightly smoother transition between the refinement zones. We thus chose the Wendland C2 as a balance between computational cost and accuracy but note that the cubic spline is likely sufficient for most purposes.", '2.1 Circumbinary disc': 'We repeat the concluding simulation from Price, Cuello, et al. (2018) of a circumbinary disc nominally representing the protoplanetary disc HD 142527 (but see Nowak et al. 2024). The star and disc parameters are the same as used in these works and our low resolution reference simulation has N = 1 × 10 6 particles, as in their work. Our high resolution reference simulation is the same as the low except it has double the spatial resolution with N = 8 × 10 6 particles. In our APR versions (simulations B1-B4) the refinement region has three nested levels of refinement ( ℓ = 3) and is centred on the centre of mass of the two stars so that the accretion streams and circumstellar discs are more highly resolved. As the gas has a negative radial velocity, we find that this application mostly represents a test of splitting - little if any merging occurs in these simulations. This is similar to the published applications of APR in GIZMO (Franchini, Lupi, and Sesana 2022; Duffell et al. 2024).  \nFigure 3 shows a representative APR simulation of HD142527 book-ended by reference calculations. At t = 1137 years in the low resolution reference calculation we find a poorly resolved circumprimary disc consistent with Price, Cuello, et al. (2018), but this disc is recovered in the high resolution case on the far right. Our APR simulation shows the same disc structure around the primary star as in the high resolution case. The spiral arms at the inner edge of the disc are well recovered when they are within the  \nFigure 5. Column density of the region surrounding a planet embedded in a disc shown in the corotating frame. The le/f\\_t and right panels have no APR and are separated by a factor of two in spatial resolution. The middle panel shows simulation P6 with three levels of refinement, matching the circumprimary disc structure and spiral arms in the high resolution reference simulation. The refinement zone is shown as in Figure 3. This test demonstrates that structures like spiral arms are faithfully reproduced even when they cross the refinement boundaries. A movie of these simulations is available online.  \n<!-- image -->  \nFigure 6. Mass and radius of the planet in the planet-disc simulation with di/fferent combinations of refinement region sizes. The APR simulations and high resolution simulations are again distinct from the low resolution reference simulation.  \n<!-- image -->  \nhighest refinement region (e.g. simulations B3 and B4). The consistency between the APR and high resolution case is maintained through the end of the simulation at t = 3283 years.  \nFigure 4 shows how the mass of the two sink particles changes throughout the reference simulations and simulations B1-B4. The evolution of the mass for the APR simulations  \nclosely follows the high resolution reference cases where the refinement region covers the cavity, indicating that the disc structure and thus the accretion rate is the same across these simulations. The low resolution simulation with the poorly resolved circumprimary disc accretes faster, resulting in a larger final mass. Simulations B1 and B2 show distinct behaviour after t ∼ 1000 years with up to ∼ 3 times more mass accreted onto both sinks due to their small refinement regions; the spiral arms and even the secondary sink are not consistently resolved inside the highest refinement region. For both B1 and B2 the lower resolution of the inner edge of the cavity then leads to about 30% more mass falling onto the sinks and thus the higher accretion rate for both.  \nWe calculate the error, E , using the largest percentage difference in mass between the APR versions and high resolution reference case on each sink as  \nE = max [ m APR ( t ) -m ref ( t ) m ref ( t ) × 100 ] . (6)  \nWerepeat this calculation across both sinks and take the largest of these as the error in our simulation when compared to the high resolution reference case, summarised in Table 1. We find the largest error of 9.0 × 10 -3 %lies with the secondary sink of the B1 simulation, possibly because it ventures closer to the edge of the refinement region each orbit. In contrast to our other example applications, our APR simulations here show a speed-up compared to the globally high resolution case of ∼ 1.1 - 2.3 × faster (Table 1). The speed up is inversely proportional to the size of the refinement zone because it is centred on the binary; the tightest and most computationally expensive orbits are located here and our APR simulations add more particles to this region.', '2.2 Planet-disc interaction': 'We simulate the interaction of a planet embedded in a disc where the region around the planet is refined using APR. The disc extends from R in = 1 au to R out = 10 au with a total  \ndisc mass of 0.05M ⊙ . We adopt typical protoplanetary disc values with a surface density profile Σ ∝ (1 -√ R in / R ) R -p with p = 1, a sound speed c s ∝ R -q with q = 0.25, a disc aspect ratio of H / R = 0.05 at the inner edge R = R in where H ≡ c s( R )/ Ω ( R ) is the scale-height, Ω ( R ) = √ GM ∗ / R 3 is the Keplerian angular velocity, M ∗ = 1M ⊙ is the central mass and R is the cylindrical radius. We also assume a Shakura and Sunyaev (1973) viscosity of α = 0.005, modelled with the method outlined in Lodato and Price (2010). Into this disc we add a planet with 0.1M J initially located at R p = 5 au, which is large enough to generate spirals in the disc but not so massive that the wakes of the spirals interact strongly with the planet. We set the accretion radius of the sink to be 0.25 R H = 0.04 au where R H is the Hill radius (Nealon et al. 2018). Our low and high reference simulations use N = 1 × 10 6 and N = 8 × 10 6 respectively. We conduct six APR simulations (P1-P6, see Table 1) varying the size of the refinement region and the size of the step around it, all with a maximum refinement level of ℓ = 3 (corresponding to the high resolution reference).  \nOur simulations with a planet show the formation of spiral arms, indications of a circumplanetary disc and a shallow gap being carved. Figure 5 compares the column density of our low and high resolution reference case to the representative example P6. For all of these simulations we find that the spiral arms generated by the planet are recovered faithfully even as they cross multiple refinement levels. The column density contrast of the circumplanetary disc and the spiral arms are visually almost identical between our APR simulations and the high resolution reference case.  \nFigure 6 shows the mass of the planet and radius measured from the central star. As with the HD 142527 simulations, the high resolution reference and the APR versions are qualitatively distinct from the low resolution case. Using Equation 6 we calculate the error from the mass but also equivalently from the radius, taking the largest percentage difference as the error reported in Table 1. We find the largest difference of 9.1% in the mass and 3.06% in the radius, both for simulation P5. While this error is the largest we recover with our APR implementation in the examples shown here it is still about 1/3 of the largest difference between the low and high resolution reference cases. Additionally, Nealon et al. (2018) showed that the choice of accretion radius has a stronger impact on the planet location and mass (their Figure A1) than the difference we have found here. In general we find that the error increases with both the size of and the width of the refinement region.', '2.3 Flyby encounter': "Smallwood et al. (2024) showed that there is a direct and robust relationship between the orientation of the disc formed from captured material in a flyby encounter and the original disc. With a broad suite of simulations, their results showed that the disc that formed around the perturber was always tilted twice as much as the original inclination between the perturber and the disc around the primary star. However the perturber disc in their simulations was generally underesolved (see their Figure A2). Smallwood et al. (2024) conducted a resolution  \nstudy with 5 × 10 5 and 4 × 10 6 particles and determined that in order to well resolve the disc around the perturber they would need to model the original disc with ∼ 144 × 10 6 particles.  \nWe repeat the calculation from Smallwood et al. (2024) at 5 × 10 5 particles for our low resolution reference and 4 × 10 6 as our high resolution reference. We then conduct four APR variations (F1 - F4) with the refinement region centred on the perturber star. The upper panel of Figure 7 shows the column density of the two reference cases and simulations F1 with a refinement region radius of 20 au and F4 with 40 au, about 300 years after pericentre passage (as in Smallwood et al. 2024). The two APR cases and the high resolution share the same disc structure but the tidal stream is different. In F1 with the smallest refinement region the tidal stream more closely resembles that of the low resolution reference simulation, in line with our expectations.  \nAs before, Figure 8 shows the mass of the perturber sink throughout the simulations. We find that the high resolution reference case and APR versions describe the same mass accretion pathway and that it is distinct from the low resolution reference case. The error measured from the mass of the perturber sink is ≲ 10 -3 %for all combinations of APR used here and with speed ups of between 5.28 - 6.62 × faster than the high resolution reference case.  \nAs in Smallwood et al. (2024) we also measure the resolution of the perturber disc for the different simulations. Figure 9 shows the ⟨ h ⟩ / H resolution at t = 2500 years after the pericentre encounter where ⟨ h ⟩ / H is the shell averaged smoothing length. We find that the high resolution reference case and APR versions produce the same profile, confirming that they are the same local resolution. In particular, the profiles for F1 and F2 deviate around R = 25au and 30 au respectively, near to the edge of the refinement region for these simulations. The factor of 1/2 in ⟨ h ⟩ / H between the APR simulations and the low resolution case demonstrates that the resolution has doubled (as expected from three refinement levels).  \nThese example applications so far have demonstrated the accuracy and speed of the APR implementation by comparing to a globally high resolution reference case, but this is not the most powerful or intended use of APR. Here we perform an additional two simulations using APR on the flyby problem: simulation F5 which initially has N = 4 × 10 6 particles and uses three levels of refinement ( ℓ = 3) and simulation F6 which initially has N = 5 × 10 5 particles and uses six levels of refinement ( ℓ = 6). Around the perturber these simulations achieve the same local resolution in terms of local particle number.  \nThe lower row of Figure 7 shows the discs formed around the perturber star for F5 and F6, finding the same disc structure as in the previous cases. These two simulations are also included in Figures 8 and 9. The mass accreted onto the sink in simulations F5 and F6 is 9.2 × 10 -4 %lower than the high resolution reference case. This difference is about half that between the low and high resolution reference cases (1.89 × 10 -3 %) and consistent with our previous observation that increasing the resolution decreases the mass accretion rate. Figure 9 demonstrates that the linear resolution of these two simulations is double the previous F1 - F4 simulations. We have thus shown  \nFigure 7. Comparison of discs formed from captured material around the perturber in the flyby simulation (see Smallwood et al. 2024). The 'Low resolution' and 'High resolution' panels do not have APR and are separated by a factor of two in resolution. Simulations F1 and F4 initially have N = 5 × 10 5 particles with 3 levels of APR, F5 has N = 4 × 10 6 particles with 3 levels of refinement and F6 has N = 5 × 10 5 particles with 6 levels of refinement. The refinement zone is shown as in Figure 3. The disc structure is similar irrespective of the base resolution of the simulation but the tidal stream onto the disc depends on the size of the refinement region and the number of levels. A movie of these simulations is available online.  \n<!-- image -->  \nFigure 8. Mass of the perturber star in the flyby simulation for di/fferent combinations of refinement regions and maximum refinement levels. As before, increasing the resolution using APR leads to a lower mass accretion rate and the similarities between the APR and high resolution reference evidence the similarity in the disc structure around the perturber star.  \n<!-- image -->  \nFigure 9. Quantifying the resolution of the perturber disc in the flyby calculations using ⟨ h ⟩ / H as a function of R . Three levels of refinement corresponds to a factor of two in linear resolution and this is recovered here by the ⟨ h ⟩ / H decreasing by about half for three refinement levels. Colour scheme is the same as in Figure 8.  \n<!-- image -->  \nthat with just 5.7 × 10 6 (F5) and 3.0 × 10 6 particles we can achieve an equivalent resolution around the perturber of simulation that uses 32 × 10 6 particles globally.", '3. Discussion': "Our results demonstrate that our method is accurate and fast (see Table 1) with the most rapid speed up seen in the flyby simulations and the slowest in the circumbinary disc simulations. In the flyby case the majority of the splitting/merging occurred during the actual flyby encounter and the time-step constraining orbits were around the primary star which we did not change in our APR simulations. Together, these led to a significant speed up in the simulation. By contrast, in the circumbinary disc we increased the resolution in the region of the simulation where particles had the largest acceleration and thus the shortest time-steps. We still found a speed up in that case because of our use of individual time-steps (Price, Wurster, et al. 2018). In the planet-disc interaction simulations we found constant splitting and merging occurred as the variation in the orbital speed meant particles were continually entering and exiting the refinement region. For those simulations the size of the refinement region generally dictated the speed up, with smaller regions corresponding to faster speed ups.  \nThe circumbinary disc and flyby tests demonstrated the importance of where to set the refinement region. In the former, the largest error was measured when the inner edge of the circumbinary disc was not adequately resolved (simulations B1 and B2) such that material was entering the cavity with the same rate as the low resolution reference case and occasionally, the secondary star popped out of the refinement region. But when the inner edge of the circumbinary disc and the complete orbit of the secondary was captured in the refinement region (B3 and B4) the error decreased significantly as compared to the high resolution reference case. This effect was also noted by Franchini, Lupi, and Sesana (2022), who recommend setting the inner edge of the refinement region to 4 × the semi-major axis.  \nOur implementation is also designed to be as adaptable as possible. Our use of the k -d tree in P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc for merging means that our APR method can naturally be used for derefinement, where the APR zone has a lower resolution than the global simulation. This is particularly useful when increases in density result in very small time-steps that effectively kill a simulation - commonly experienced in simulations including self-gravity (e.g. Longarini et al. 2023; Lau et al. 2022; Hall, Forgan, and Rice 2017). Our method is also easily extendable to multiple refinement regions (e.g. refining separately around two planets in a disc).  \nAn important and desirable property of SPH is its conservation properties. Conservation of kinetic energy and both linear and angular momentum are respected with APR when the children particles inherit the velocity of the parent (Feldman and Bonet 2007; López, Roose, and Recarey Morfa 2013). While our inclusion of relaxing means that children may not necessarily be placed symmetrically around the parent particle this only affects the conservation of angular momentum. Our  \nmethod perfectly conserves mass and in simulations with accretion (where we may get an odd number of particles to be merged) our implementation caps the number of unmerged particles to be ℓ max - 1 (Equation 4). In addition to this, APR in P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc is subjected to the regular conservation checks (see Section 2.2, Price, Wurster, et al. 2018).  \nThe current limitation of our method is a ∼ 5% blip in density that occurs every time particles are either split or merged. In our 3D simulations this 'blip' in density is cosmetic (it can just be seen in the middle panel of Figure 5) but is more obvious in our tests on small amplitude linear sound waves outlined in Appendix 1 (although it does not seem to affect the propagation of the wave itself ). We tested different methods to mitigate this feature and found that reducing n child was ultimately the most effective measure, setting our limitation of n child = 2 and thus having nested refinement regions. Importantly, the 'blip' feature is also visibly present in the GIZMO implementation of APR (see density renderings in Duffell et al. 2024) which demonstrates that it is not a function of our method, but is inherent to the process of changing the particle mass.  \nThe most promising method to remove this blip is the blending method outlined by Barcarolo et al. (2014) which allows particles to be introduced in a way that removes the noise before they contribute to the density summation. Although successful for incompressible flows, blending is practically very difficult to implement with compressible flows due to mass conservation in the blending zone. For example, we could successfully incorporate blending as in Barcarolo et al. (2014) in our initial conditions, but once particles started moving through the blending region we found that the solution deteriorated rapidly. We additionally implemented other forms of relaxing (e.g. Lind et al. 2012; Diehl et al. 2015; Sun et al. 2017) and other forms of blending (Chiron et al. 2018; Gao, Qiu, and Fu 2022) but did not find a satisfactory improvement on the blip found in the simulations.", '4. Conclusion': 'We introduce a live adaptive particle refinement implementation into the smoothed particle hydrodynamics code P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc. We have considered example applications of a circumbinary disc, planet-disc interaction and a flyby to demonstrate our method. With these examples we have shown that our implementation is  \n- 1. Accurate: We measure the mass accreted onto sink particles as a proxy of the disc structure in our simulations. For the circumbinary and flyby examples, we find that the APR calculations are accurate to < 0.1%. For the planet-disc interaction example we measure both the mass and the radial location of the planet and find that it is accurate to at least 9% (but is generally more accurate than this).\n- 2. Fast: Every APR simulation showed a speed up, offering between 1.07-6.62 × than when compared to a simulation with the same resolution but globally. The speed up is application-dependent, with the flyby example being the most rapid and the circumbinary disc the least amenable to speedup in this manner.  \n- 3. Uses less storage: Because our APR simulations use fewer particles in total, they require between 15 - 27% of the storage of an equivalent globally resolved calculation.  \nOur example applications suggested optimal sizes of the refinement region as a guide for future calculations. We found accuracy of the implementation depends sensitively on whether or not key features were uniformly resolved. We also demonstrated that the location of the refinement region can be dynamic and note that derefinement is possible. Finally we showed that APR can increase the resolution of simulations at low cost; for the flyby example we achieved a local resolution of 32 million particles with totals of either 3.0 or 5.7 million particles.  \nWe showed examples limited to hydrodynamical simulations, but P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc also includes dust (Guillaume Laibe and Daniel J. Price 2012; G. Laibe and D. J. Price 2014), magnetohydrodynamics (Tricco and Price 2012), self-gravity, general relativity (Liptai and Price 2019) and can be coupled with MCFOST (Pinte et al. 2006; Pinte et al. 2009; Nealon, Price, and Pinte 2020). For brevity we leave dedicated testing of our APR implementation with these features to future works.', 'Acknowledgement': 'The authors thank James Wurster for significant contributions in the early stages of this work and the referee for helpful suggestions. RN acknowledges Matthew Bate, Richard Booth, Terry Tricco and Alessia Franchini for discussions and Sahl Rowther for aesthetics assistance.  \nFunding Statement R.N. acknowledges support from UKRI/EPSRC through a Stephen Hawking Fellowship (EP/T017287/1). This work was performed using Avon, the HPC clusters at the University of Warwick. DP is grateful for Australian Research Council funding via DP220103767 and DP240103290. We thank the Monash-Warwick alliance for sponsorship of the /P.sc/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc users workshop where this work commenced and Nicolás Cuello for his hospitality.  \nData Availability Statement The data to reproduce figures in this study are openly available on Zenodo at http://doi.org/10.5281/zenodo.11581005. All data from the simulations will be made available upon request to the corresponding author. The software used to create and visualise the simulations are publicly available:  \nP/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc:  \nhttps://github.com/danieljprice/phantom (Price, Wurster, et al. 2018) https://github.com/ttricco/sarracen  \nS/A.sc/R.sc/R.sc/A.sc/C.sc/E.sc/N.sc:  \n(Harris and Tricco 2023)', 'References': "Anglés-Alcázar, Daniel, Eliot Quataert, Philip F. Hopkins, Rachel S. Somerville, Christopher C. Hayward, Claude-André Faucher-Giguère, Greg L. Bryan, Dušan Kereš, Lars Hernquist, and James M. Stone. 2021. Cosmological Simulations of Quasar Fueling to Subparsec Scales Using Lagrangian Hyper-refinement. ApJ 917, no. 2 (August): 53. https : / / doi . org / 10 . 3847 / 1538 - 4357 / ac09e8. arXiv: 2008 . 12303 [astro-ph.GA] .  \nBarcarolo, D.A., D. Le Touzé, G. Oger, and F. de Vuyst. 2014. Adaptive particle refinement and derefinement applied to the smoothed particle hydrodynamics method. Journal of Computational Physics 273:640-657. /I.sc/S.sc/S.sc/N.sc: 0021-9991. https://doi.org/https://doi.org/10.1016/j.jcp. 2014.05.040. https://www.sciencedirect.com/science/article/pii/ S0021999114004033.  \nBerger, M. J., and P. Colella. 1989. Local Adaptive Mesh Refinement for Shock Hydrodynamics. Journal of Computational Physics 82, no. 1 (May): 64-84. https://doi.org/10.1016/0021-9991(89)90035-1.  \nBertschinger, Edmund. 2001. Multiscale Gaussian Random Fields and Their Application to Cosmological Simulations. ApJS 137, no. 1 (November): 1-20. https://doi.org/10.1086/322526. arXiv: astro- ph/0103301 [astro-ph] .  \nBørve, S., M. Omang, and J. Trulsen. 2001. Regularized Smoothed Particle Hydrodynamics: A New Approach to Simulating Magnetohydrodynamic Shocks. ApJ 561, no. 1 (November): 82-93. https://doi.org/10. 1086/323228.  \n. 2005. Regularized smoothed particle hydrodynamics with improved multi-resolution handling. Journal of Computational Physics 208, no. 1 (September): 345-367. https://doi.org/10.1016/j.jcp.2005.02.018.  \n. 2006. Multidimensional MHD Shock Tests of Regularized Smoothed Particle Hydrodynamics. ApJ 652, no. 2 (December): 1306-1317. https: //doi.org/10.1086/508454.  \nBørve, S., R. Speith, and J. Trulsen. 2009. Numerical Dissipation in RSPH Simulations of Astrophysical Flows with Application to Protoplanetary Disks. ApJ 701, no. 2 (August): 1269-1282. https://doi.org/10.1088/ 0004-637X/701/2/1269.  \nChiaki, Gen, and Naoki Yoshida. 2015. Particle splitting in smoothed particle hydrodynamics based on Voronoi diagram. MNRAS 451, no. 4 (August): 3955-3963. https://doi.org/10.1093/mnras/stv1227. arXiv: 1504.07347 [astro-ph.IM] .  \nChiron, L., G. Oger, M. de Leffe, and D. Le Touzé. 2018. Analysis and improvements of adaptive particle refinement (apr) through cpu time, accuracy and robustness considerations. Journal of Computational Physics 354:552-575. /I.sc/S.sc/S.sc/N.sc: 0021-9991. https://doi.org/https://doi.org/10.1016/ j.jcp.2017.10.041. https://www.sciencedirect.com/science/article/pii/ S0021999117308082.  \nDehnen, W., and H. Aly. 2012. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. MNRAS 425 (September): 1068-1082. https://doi.org/10.1111/j.1365-2966.2012. 21439.x. arXiv: 1204.2471 [astro-ph.IM] .  \nDiehl, S., G. Rockefeller, C. L. Fryer, D. Riethmiller, and T. S. Statler. 2015. Generating optimal initial conditions for smoothed particle hydrodynamics simulations. Publications of the Astronomical Society of Australia 32:e048. https://doi.org/10.1017/pasa.2015.50.  \nDuffell, Paul C., Alexander J. Dittmann, Daniel J. D'Orazio, Alessia Franchini, Kaitlin M. Kratter, Anna B. T. Penzlin, Enrico Ragusa, et al. 2024. The Santa Barbara Binary-Disk Code Comparison. arXiv e-prints (February): arXiv:2402.13039. https://doi.org/10.48550/arXiv.2402.13039. arXiv: 2402.13039 [astro-ph.SR] .  \nFeldman, J., and J. Bonet. 2007. Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems. International Journal for Numerical Methods in Engineering 72, no. 3 (October): 295324. https://doi.org/10.1002/nme.2010.  \nFranchini, Alessia, Alessandro Lupi, and Alberto Sesana. 2022. Resolving Massive Black Hole Binary Evolution via Adaptive Particle Splitting. ApJL 929, no. 1 (April): L13. https://doi.org/10.3847/2041-8213/ac63a2. arXiv: 2201.05619 [astro-ph.HE] .  \nGao, Tianrun, Huihe Qiu, and Lin Fu. 2022. A block-based adaptive particle refinement sph method for fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering 399:115356. /I.sc/S.sc/S.sc/N.sc: 00457825. https://doi.org/https://doi.org/10.1016/j.cma.2022.115356. https: //www.sciencedirect.com/science/article/pii/S0045782522004376.  \nGingold, R. A., and J. J. Monaghan. 1977. Smoothed particle hydrodynamics Theory and application to non-spherical stars. MNRAS 181 (November): 375-389.  \nHall, Cassandra, Duncan Forgan, and Ken Rice. 2017. Identifying and analysing protostellar disc fragments in smoothed particle hydrodynamics simulations. MNRAS 470, no. 3 (September): 2517-2538. https : / / doi . org / 10 . 1093 / mnras / stx1244. arXiv: 1705 . 06690 [astro-ph.EP] .  \nHarris, Andrew, and Terrence S. Tricco. 2023. Sarracen: a python package for analysis and visualization of smoothed particle hydrodynamics data. Journal of Open Source Software 8 (86): 5263. https://doi.org/10.21105/ joss.05263. https://doi.org/10.21105/joss.05263.  \nHopkins, P. F. 2014. GIZMO: A New Class of Accurate, Mesh-Free Hydrodynamic Simulation Methods. ArXiv e-prints (September). arXiv: 1409.7395.  \nKatz, Neal, and Simon D. M. White. 1993. Hierarchical Galaxy Formation: Overmerging and the Formation of an X-Ray Cluster. ApJ 412 (August): 455. https://doi.org/10.1086/172935.  \nKitsionas, S., and A. P. Whitworth. 2002. Smoothed Particle Hydrodynamics with particle splitting, applied to self-gravitating collapse. MNRAS 330, no. 1 (February): 129-136. https://doi.org/10.1046/j.13658711.2002.05115.x. arXiv: astro-ph/0203057 [astro-ph] .  \n. 2007. High-resolution simulations of clump-clump collisions using SPH with particle splitting. MNRAS 378, no. 2 (June): 507-524. https: //doi.org/10.1111/j.1365-2966.2007.11707.x. arXiv: astro-ph/0703177 [astro-ph] .  \nKlypin, Anatoly, Andrey V. Kravtsov, James S. Bullock, and Joel R. Primack. 2001. Resolving the Structure of Cold Dark Matter Halos. ApJ 554, no. 2 (June): 903-915. https://doi.org/10.1086/321400. arXiv: astroph/0006343 [astro-ph] .  \nLaibe, G., and D. J. Price. 2014. Dusty gas with one fluid in smoothed particle hydrodynamics. MNRAS 440 (May): 2147-2163. https://doi.org/10. 1093/mnras/stu359. arXiv: 1402.5249 [astro-ph.EP] .  \nLaibe, Guillaume, and Daniel J. Price. 2012. Dusty gas with smoothed particle hydrodynamics - I. Algorithm and test suite. MNRAS 420, no. 3 (March): 2345-2364. https://doi.org/10.1111/j.1365-2966.2011.20202.x. arXiv: 1111.3090 [astro-ph.IM] .  \nLastiwka, Martin, Nathan Quinlan, and Mihai Basa. 2005. Adaptive particle distribution for smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids 47, nos. 10-11 (April): 1403-1409. https: //doi.org/10.1002/fld.891.  \nLau, Mike Y. M., Ryosuke Hirai, Miguel González-Bolívar, Daniel J. Price, Orsola De Marco, and Ilya Mandel. 2022. Common envelopes in massive stars: towards the role of radiation pressure and recombination energy in ejecting red supergiant envelopes. MNRAS 512, no. 4 (June): 54625480. https://doi.org/10.1093/mnras/stac049. arXiv: 2111.00923 [astro-ph.SR] .  \nLind, S.J., R. Xu, P.K. Stansby, and B.D. Rogers. 2012. Incompressible smoothed particle hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. Journal of Computational Physics 231 (4): 1499-1523. /I.sc/S.sc/S.sc/N.sc: 0021-9991. https://doi.org/https://doi.org/10.1016/j.jcp.2011.10.027. https: //www.sciencedirect.com/science/article/pii/S0021999111006279.  \nLiptai, David, and Daniel J. Price. 2019. General relativistic smoothed particle hydrodynamics. MNRAS 485, no. 1 (May): 819-842. https://doi.org/ 10.1093/mnras/stz111. arXiv: 1901.08064 [astro-ph.IM] .  \nLodato, G., and D. J. Price. 2010. On the diffusive propagation of warps in thin accretion discs. MNRAS 405 (June): 1212-1226. https://doi.org/10. 1111/j.1365-2966.2010.16526.x. arXiv: 1002.2973 [astro-ph.HE] .  \nLongarini, Cristiano, Philip J. Armitage, Giuseppe Lodato, Daniel J. Price, and Simone Ceppi. 2023. The role of the drag force in the gravitational stability of dusty planet-forming disc - II. Numerical simulations. MNRAS 522, no. 4 (July): 6217-6235. https://doi.org/10.1093/mnras/stad1400. arXiv: 2305.03659 [astro-ph.EP] .  \nLópez, Yaidel, Dirk Roose, and Carlos Recarey Morfa. 2013. Dynamic particle refinement in sph: application to free surface flow and non-cohesive soil simulations. Computational Mechanics 51 (May). https://doi.org/10.1007/ s00466-012-0748-0.  \nLucy, L. B. 1977. A numerical approach to the testing of the fission hypothesis. AJ 82 (December): 1013-1024. https://doi.org/10.1086/112164.  \nMeglicki, Z., D. Wickramasinghe, and G. V. Bicknell. 1993. ThreeDimensional Structure of Truncated Accretion Discs in Close Binaries. MNRAS 264 (October): 691. https://doi.org/10.1093/mnras/264.3.691.  \nMonaghan, J. J., and J. C. Lattanzio. 1985. A refined particle method for astrophysical problems. AAp 149 (August): 135-143.  \nMonaghan, J. J., and S. R. Varnas. 1988. The dynamics of interstellar cloud complexes. MNRAS 231 (March): 515-534. https://doi.org/10.1093/ mnras/231.3.515.  \nNavarro, Julio F., and Simon D. M. White. 1994. Simulations of dissipative galaxy formation in hierarchically clustering universes-2. Dynamics of the baryonic component in galactic haloes. MNRAS 267, no. 2 (March): 401-412. https://doi.org/10.1093/mnras/267.2.401.  \nNealon, Giovanni Dipierro, Richard Alexander, Rebecca G. Martin, and Chris Nixon. 2018. Warping a protoplanetary disc with a planet on an inclined orbit. MNRAS 481 (November): 20-35. https://doi.org/10.1093/mnras/ sty2267. arXiv: 1808.05693 [astro-ph.EP] .  \nNealon, Rebecca, Daniel J. Price, and Christophe Pinte. 2020. Rocking shadows in broken circumbinary discs. MNRAS 493, no. 1 (March): L143L147. https://doi.org/10.1093/mnrasl/slaa026. arXiv: 2002.02983 [astro-ph.EP] .  \nNowak, M., S. Rowther, S. Lacour, F. Meru, R. Nealon, and D. J. Price. 2024. The orbit of HD 142527 B is too compact to explain many of the disc features. arXiv e-prints (February): arXiv:2402.03595. https://doi.org/10. 48550/arXiv.2402.03595. arXiv: 2402.03595 [astro-ph.SR] .  \nPinte, C., T. J. Harries, M. Min, A. M. Watson, C. P. Dullemond, P. Woitke, F. Ménard, and M. C. Durán-Rojas. 2009. Benchmark problems for continuum radiative transfer. High optical depths, anisotropic scattering, and polarisation. AAp 498 (May): 967-980. https://doi.org/10.1051/ 0004-6361/200811555. arXiv: 0903.1231 [astro-ph.SR] .  \nPinte, C., F. Ménard, G. Duchêne, and P. Bastien. 2006. Monte Carlo radiative transfer in protoplanetary disks. AAp 459 (December): 797-804. https: //doi.org/10.1051/0004-6361:20053275. eprint: astro-ph/0606550.  \nPorter, D. H. 1985. A Study of Hierarchical Clustering of Galaxies in an Expanding Universe. PhD diss., University of California, Berkeley, January.  \nPrice, N. Cuello, C. Pinte, D. Mentiplay, S. Casassus, V. Christiaens, G. M. Kennedy, et al. 2018. Circumbinary, not transitional: on the spiral arms, cavity, shadows, fast radial flows, streamers, and horseshoe in the HD 142527 disc. MNRAS 477 (June): 1270-1284. https://doi.org/10.1093/ mnras/sty647. arXiv: 1803.02484 [astro-ph.SR] .  \nPrice, Daniel J., David Liptai, Ilya Mandel, Joanna Shepherd, Giuseppe Lodato, and Yuri Levin. 2024. Eddington envelopes: The fate of stars on parabolic orbits tidally disrupted by supermassive black holes. arXiv e-prints (April): arXiv:2404.09381. https://doi.org/10.48550/arXiv.2404.09381. arXiv: 2404.09381 [astro-ph.HE] .  \nPrice, James Wurster, Terrence S. Tricco, Chris Nixon, Stéven Toupin, Alex Pettitt, Conrad Chan, et al. 2018. Phantom: A Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code for Astrophysics. PASA 35 (September): e031. https://doi.org/10.1017/pasa.2018.25. arXiv: 1702.03930 [astro-ph.IM] .  \nShakura, N.I., and R.A. Sunyaev. 1973. Black Holes in Binary-Systems - Observational Appearance. AAp (17, AVE DU HOGGAR, PA COURTABOEUF, BP 112, F-91944 LES ULIS CEDEX A, FRANCE) 24 (3): 337-355. /I.sc/S.sc/S.sc/N.sc: 0004-6361.  \nSmallwood, Jeremy L., Rebecca Nealon, Nicolás Cuello, Ruobing Dong, and Richard A. Booth. 2024. Formation of misaligned second-generation discs through fly-by encounters. MNRAS 527, no. 2 (January): 20942109. https://doi.org/10.1093/mnras/stad3057. arXiv: 2310.03327 [astro-ph.EP] .  \nSpringel, Volker. 2005. The cosmological simulation code GADGET-2. MNRAS 364, no. 4 (December): 1105-1134. https://doi.org/10.1111/j.13652966.2005.09655.x. arXiv: astro-ph/0505010 [astro-ph] .  \nSun, P.N., A. Colagrossi, S. Marrone, and A.M. Zhang. 2017. The δ plussph model: simple procedures for a further improvement of the sph scheme. Computer Methods in Applied Mechanics and Engineering 315:2549. /I.sc/S.sc/S.sc/N.sc: 0045-7825. https://doi.org/https://doi.org/10.1016/j.cma. 2016.10.028. https://www.sciencedirect.com/science/article/pii/ S0045782516309112.  \nThacker, R. J., and H. M. P. Couchman. 2000. Implementing Feedback in Simulations of Galaxy Formation: A Survey of Methods. ApJ 545, no. 2 (December): 728-752. https : / / doi . org / 10 . 1086 / 317828. arXiv: astro-ph/0001276 [astro-ph] .  \nTricco, Terrence S, and Daniel J Price. 2012. Constrained hyperbolic divergence cleaning for smoothed particle magnetohydrodynamics. Journal of Computational Physics 231 (21): 7214-7236.  \nTruelove, J. Kelly, Richard I. Klein, Christopher F. McKee, II Holliman John H., Louis H. Howell, and Jeffrey A. Greenough. 1997. The Jeans Condition: A New Constraint on Spatial Resolution in Simulations of Isothermal Self-gravitational Hydrodynamics. ApJL 489, no. 2 (November): L179-L183. https://doi.org/10.1086/310975.  \nVacondio, R., B. D. Rogers, P. K. Stansby, P. Mignosa, and J. Feldman. 2013. Variable resolution for SPH: A dynamic particle coalescing and splitting scheme. Computer Methods in Applied Mechanics and Engineering 256 (April): 132-148. https://doi.org/10.1016/j.cma.2012.12.014.  \nWendland, H. 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in computational Mathematics 4 (1): 389-396.  \nYang, Xiufeng, and Song-Charng Kong. 2019. Adaptive resolution for multiphase smoothed particle hydrodynamics. Computer Physics Communications 239:112-125. /I.sc/S.sc/S.sc/N.sc: 0010-4655. https://doi.org/https://doi.org/10. 1016/j.cpc.2019.01.002. https://www.sciencedirect.com/science/article/ pii/S0010465519300037.", 'Appendix 1. Simple tests': 'To demonstrate some of the key choices of our implementation we use the straightforward example of a sound wave in a periodic box simulated with P/H.sc/A.sc/N.sc/T.sc/O.sc/M.sc. The density is set by  \nρ ( x ) = ρ 0 + A sin ( 2 π x L ) , (7)  \nwhere A = 0.02 is the amplitude of the perturbation and ρ 0 = 1. The thermal energy is perturbed in the same manner with the same amplitude A but scaled by P 0 / ρ 0 , the unperturbed pressure and density. The fluid has an unperturbed sound speed of c s,0 = 1.0 and we adopt an adiabatic equation of state. Unless otherwise stated the velocity is set to vx = 0.001, vy = vz = 0. The particles are initially set on a close-packed lattice with the above density distribution. We add noise to the lattice by randomly adjusting the x and y positions of the particles by 0.0001 × the length of the box. We use 64 particles across the box of width dx = 1.0 and dy = 1.0, corresponding to a total of N = 56832 particles initially.  \nThe refinement region is a circle centred at x = 0.5, y = 0.5 with a radius of r = 0.10. In contrast to our example applications, these simulations only make use of one refinement level ( ℓ = 1) which corresponds to a doubling of the particle number inside this region. The first split occurs at t = 0.5 to allow the particles to settle a bit first and the simulation is then performed until t = 2.0, corresponding to five wave-crossings in total. After the split the simulation has around N = 63500 particles.  \nOur results in this section are all shown with the cross section of the density, ρ ( x ), across the box divided by the initial average density, ρ IC , to avoid small variations in the background density caused by kernel bias. For clarity we also restrict the visualisation to only include particles between 0.4 < y < 0.6, essentially making it a cross section centred on y = 0.5. On each plot we also display the average density as a function of position across the box for both the APR simulation and a reference simulation without any refinement. This allows us to easily compare how much noise is introduced as a result of the splitting/merging and how accurate these simulations are compared to the simulations that do not include APR. Movies of some of the simulations in this section are available online.', "Appendix 1.1 Density discontinuity - the 'blip'": "A common feature in all of our simulations is a density discontinuity at the refinement boundary (i.e. 'the blip' between the different resolution regions). It appears as a slight increase in the density just outside the refinement region and a slight decrease immediately inside. This discontinuity does not introduce spurious waves, has a constant amplitude and does not appear to affect the density profile inside the refinement region. Its apparent prominence is because of our use of a comparatively small amplitude sound wave as our test case we chose it specifically so that any spurious features would be easily identifiable. This feature maintains a constant amplitude (here, of < 5% of the original density) and in our examples in Section 2 is barely visible in renderings (e.g. Figure 5),  \nsimilar to the density profiles seen in examples with GIZMO (Anglés-Alcázar et al. 2021; Franchini, Lupi, and Sesana 2022; Duffell et al. 2024). Importantly, the blip does not seem to affect the propagation of the wave itself.  \nThis density blip is a manifestation of a known issue in APR, where the splitting and merging procedure introduces small errors into the density distribution. Feldman and Bonet (2007) formally established that an error in the density distribution is introduced each time a parent particle is replaced by children particles when the children all have the same mass. When a split occurs and the children particles are introduced, they should be placed symmetrically around the parent to ensure conservation of angular momentum (e.g. López, Roose, and Recarey Morfa 2013). However, the introduction of these new particles also introduces this error, the magnitude of which is a function of the distance the children are placed at r sep. In turn, this means the error can be minimised by adjusting the locations of the children.  \nLópez, Roose, and Recarey Morfa (2013) extended this idea to consider the error introduced in the derivative of the density, suggesting that this was the more relevant error to minimise for subsequent calculations. Further works looked at how best to minimise this error that is introduced every time particles are split or merged with various r sep (e.g. López, Roose, and Recarey Morfa 2013), particle shuffling/relaxing (e.g. Lind et al. 2012; Diehl et al. 2015; Sun et al. 2017), the use of ghost particles and blending techniques (e.g. Barcarolo et al. 2014; Chiron et al. 2018; Gao, Qiu, and Fu 2022) to mitigate this at the boundaries of refinement regions. Barcarolo et al. (2014) eliminated this density blip by employing a blending zone between each refinement level, where both children and parents co-exist but their contribution to the density summation is graduated by their distance across the blending zone. This approach was also adopted by Gao, Qiu, and Fu (2022) but only as new particles moved into the higher refinement level, allowing them to regularise first. While promising, blending is difficult to implement in compressible flows due to mass conservation in the blending zone.  \nThe alternative is particle relaxation (shuffling): to minimise the error introduced when particles were split, Diehl et al. (2015) used a WVT shuffling method to rearrange the particles to a lower error state compared to the original parent distribution. Yang and Kong (2019) employed multiple, stepped refinement regions to control the ratio of the smoothing lengths across the boundaries and mitigate instabilities (similar to Børve, Omang, and Trulsen 2001). In our implementation we employ relaxation, but find that it is only important when a large number of particles are split or merged at once.", 'Appendix 1.2 Kernel choice': "Figure 10 shows the final snapshot of our test simulation comparing the use of the Cubic spline (Monaghan and Lattanzio 1985), Wendland C2, Wendland C4 and Wendland C6 kernels (Wendland 1995), summarised in Table 2. The density profile in the refinement region (blue line) is in excellent agreement  \nFigure 10. Wave in a box test showing the e/ffect of di/fferent kernels with APR. The blue points show the particles within the y = 0.5 cross-section across the x dimension of the box at t = 2.0 . The dark blue line shows their average density, while the red line shows the average density of the same simulation without APR. The density discontinuity at the refinement boundary maintains a constant amplitude across all of the simulations.  \n<!-- image -->  \nwith the expected profile (red line) for all kernel choices. As might be expected, the kernels with more neighbours have less noise inside the refinement region and at the boundary. In particular, the cubic has the largest noise at the refinement region boundary, the largest spread in density inside the refinement region and has the least accurate density profile (although the difference here is marginal).  \nFigure 11 shows our planet disc interaction test at t = 10 orbits but with different combinations of r sep and either a Cubic, Wendland C2 or Wendland C4 kernel. Here we have zoomed in to see the boundaries of the refinement zone to examine the effect of these choices. Moving between kernels, the WC4 in the lowest row is smoothest while the cubic has the largest noise at the boundaries. For any of the kernels, increasing r sep beyond 0.2 also corresponds to an increase in the noise at the boundary (seen most clearly in the WC4 case). Both of these tests confirm that a higher order kernel is more effective at mitigating the density 'blip' but that the difference is marginal. For our example applications we adopt the Wendland C2 kernel as a balance between accuracy and computational expense.", 'Appendix 1.3 Adjusting the split': 'Figure 12 demonstrates the effectiveness of both particle relaxing (as per Section 1.4) and directional splitting, where particles are split tangentially to the refinement boundary. Importantly, once the particles have settled after their initial split (Figure 12, top panel) there is no significant difference between the methods in the density profiles. However, their time to compute is quite different with the simulations that include relaxing  \nFigure 11. Planet disc interaction test showing the e/ffect of di/fferent kernel choices and r sep on the noise introduced at the refinement boundary in a full simulation. Refinement boundaries are indicated in the same way as Figure 3. The boundary is smoothest for the Wendland C4 kernel when r sep = 0.2 but we note the di/fference is marginal.  \n<!-- image -->  \nTable 2. Summary of the h fact employed and the average number of neighbours N neigh for the kernels tested in Figure 10 (e.g. Price, Wurster, et al. 2018).  \ntaking much longer. In practice, we thus adopt directional splitting as it is fast, gives an improved initial density estimate and prevents particles from being immediately placed across the boundary into a low resolution region when they are split. To improve computational time by default we only employ the relaxing routine when refinement regions are activated and there is a large number of particles that are split or merged at once.', 'Appendix 1.4 Number of children': 'In our tests so far we have used n child = 2, splitting parents into two children and merging two children to become a parent. In Figure 13 we show the impact of this choice by considering n child = 4 and n child = 8. We compare this to simulations with n child = 2 and ℓ = 2, 3 respectively which have equivalent local resolutions.  \nTosplit one parent into four children we follow the method outlined in Section 1.2 and treat those two children as the diagonal corners of a square. We then add an additional two particles on the opposing diagonal which ensures that the face of the square of four children particles is parallel to the refinement boundary. For eight children we create two squares that are offset by 0.1 × the smoothing length of the original parent, with the centre of mass of the resultant cube centred on the original parent. For our merging routines we simply edit our modified k -d tree to return cells with n child = 4 and n child = 8.  \nFigure 13 shows the planet disc interaction test at t = 4 orbits for these choices, where the centrally refined zone has the same width ( r ℓ = 0.35au and the steps into the refinement  \nFigure 12. Wave in a box test showing the e/ffect of di/fferent particle placement options when a split or merge occurs shown at the initial split (upper, t = 0.5 ) and at the end of the simulation (lower, t = 2.0 ). The colours are the same as in Figure 10. Relaxing is most successful when the split first occurs but makes negligible di/fference during the course of the simulation. Directional splitting also makes little di/fference in the long term but does prevent particles from splitting across the boundary.  \n<!-- image -->  \nzone at increments of 0.10au). The noise introduced by the nested zones is visibly lower in both cases; in the n child = 8 case the noise is the largest and is also prominent across the whole refined region, even after it has had the opportunity to settle after the first split. We note that while this test may be improved by relaxing at each step, in practice this becomes computationally extremely expensive. This implies that using a smaller number of children - even if it means having nested refinement regions - is the preferred approach. In other words, even four children is too many.  \nFigure 13. The planet disc interaction test, examining the e/ffect of having di/fferent n child . The le/f\\_t column has n child = 2 with nested refinement levels, the upper right has n child = 4 and one level of refinement, the lower right n child = 8 and one level. The refinement zones are shown as in Figure 3. Rows have the same local resolution around the planet. While the nested refinement zones do add noise at the boundaries, this is demonstrably less than bigger numbers of children.  \n<!-- image -->', 'Appendix 1.5 Size of the APR zone': 'We test the effect of the size of the refinement region using the jet test presented in Price et al. (2024), Appendix E, Figure 14. In this test a jet of gas shoots in one direction from a fixed location in an empty domain, the gas flares out as it expands slightly. Based on the results of Price et al. (2024) we set the Mach number of the gas to 10 to ensure a fairly narrow jet. The particles are injected as a cylinder with radius R = 1 16 particles in each layer. The rest of the parameters in this simulation are scale free.  \nWe apply four different APR zones to this jet test with ℓ = 1 for all but r ℓ = 3, 5, 10 and 15. The refinement boundary of each different zone has the same position of y = 5 so that particles are always split at the same location but they are merged at different distances from the launching point of the jet. Figure 14 shows this with the refinement zone indicated with a green dashed circle. We chose an injection velocity of vy = 37.5 and calculate how many sound crossing times the particles will experience between refinement and derefinement to be < 1, 1.3, 2.7 and 4.0.  \nFigure 14 shows the column density of the jets at t = 50 with the refinement zones superimposed. For the tests where particles have fewer sound crossing times we find that there is a slight narrowing of the jet within the refinement zone but for the final test where particles have several sound crossing times this narrowing is not apparent. This test demonstrates that particles require several sound crossing times between the refinement boundaries to allow them to relax before they are derefined.  \nFigure 14. Testing the size of the refinement zone, where the four simulations are characterised by the number of sound crossings that can occur between the refinement and derefinement boundaries. The refinement region is indicated with a green circle in the figure.  \n<!-- image -->'}