{"id":"86489064-014b-496a-8515-da1da54b7e9b","arxiv_id":"2610.28727","paper_version":1,"verdict":"ACCEPT","confidence":"HIGH","novelty_score":7.0,"correctness_risk":"low","formal_verification":"none","parameter_count":8,"one_line_summary":"A new hybrid framework combines collisional stellar dynamics with moving-mesh magnetohydrodynamics, enabling simulations that resolve close encounters and ejections during star cluster formation.","lead":"Arepo-CD couples star-by-star collisional N-body dynamics to the moving-mesh MHD code Arepo, using unsoftened star-star gravity and a hierarchy of integrators for compact stellar groups. Benchmark tests against the established N-body code PeTar show close agreement in cluster evolution, and application to a turbulent giant molecular cloud produces a population of runaway stars ejected at speeds above 100 km/s that softened-gravity simulations miss.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek/deepseek-v4-flash","headline":"Parameter Δt_group (≈10 yr) determines which stars receive few-body integration; its sensitivity on cluster dynamics and runaway production is not tested, raising uncertainty about the robustness of the central results.","rationale":"The reader identified the group-identification parameters as the weakest assumption. I fully agree: Δt_group, N_FSI_max, and N_SDAR_max are free parameters that directly control which stars enter the collisional few-body solver. Among these, Δt_group is the most sensitive because it sets the dynamical timescale threshold for pair identification; N_FSI_max and N_SDAR_max primarily affect computational cost and are less physically decisive as long as they are large enough to capture the relevant subsystems. The paper presents compelling benchmarks against PeTar and a plausible GMC application, but these use a single parameter set. Without a sensitivity analysis, the central claim that Arepo-CD produces a substantial runaway population is contingent on the particular choice of Δt_group. Given that the paper otherwise meets the standard for acceptance (method clear, benchmarks convincing, application physically interesting), I recommend conditional acceptance pending the addition of a sensitivity test (as described above) or a clear physical argument for why the default choice is robust. This does not reject the paper but asks for a modest strengthening that would significantly raise confidence in the quantitative results.","tokens_in":31794,"tokens_out":3597,"duration_ms":39401,"concrete_test":"Run the same GMC simulation (Sect. 4) with Δt_group set to 5 yr and 20 yr (factor-of-two variations around the fiducial 10 yr). Compare the following at t = 2.61 Myr (1.91 t_ff): (a) the number of isolated stars with v > 30 km/s, (b) the cumulative velocity distribution above 10 km/s, and (c) the radial density profile between 0.1 and 10 pc. If these quantities agree within the stochastic scatter of the simulation (estimated by comparing two realizations at the fiducial value), the threshold choice is not load-bearing. A cheaper alternative: vary Δt_group in the N=10^4 IMF cluster benchmark (Sect. 3.2) and compare the velocity distribution at 10 Myr.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The group-identification procedure (Sects. 2.4.1, 2.4.2) uses Δt_group (typically 10 yr) to define collisional pairs. Pairs with Δt_ff or Δt_fb < Δt_group are merged into groups that receive specialized few-body integration. This threshold directly controls which stars are treated with unsoftened gravity and high-order integrators, while the rest retain softened gravity and lower-order integration. The paper's central claims—that Arepo-CD reproduces collisional benchmarks and, in a GMC simulation, produces a substantial runaway population—depend on this choice. However, no sensitivity study is presented: the effect of varying Δt_group on core collapse, density profiles, or the runaway tail is not explored. The benchmark tests (Sect. 3) use a single Δt_group value; it is not varied even in the simpler N-body setups. Because the threshold is user-defined and may interact with the numerical timestep hierarchy and the group-size limits (N_FSI_max, N_SDAR_max), the absence of a sensitivity test leaves a load-bearing assumption unvalidated. A different choice could shift the boundary between group-integrated and tree-integrated stars, altering both energy conservation and the production of high-velocity escapers. Without evidence that the results are robust to plausible variations in Δt_group, the quantitative predictions (e.g., runaway fraction) cannot be considered firmly established.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek/deepseek-v4-flash","summary":"This paper presents Arepo-CD, a new computational framework that couples star-by-star collisional N-body dynamics with the moving-mesh magnetohydrodynamics code Arepo. Unsoftened gravity is used for star-star interactions, while softened gravity is retained for interactions involving gas and collisionless components. Compact stellar subsystems are identified on the fly and evolved with a hybrid integrator combining a fourth-order forward symplectic integrator (FSI) and slow-down algorithmic regularization (SDAR). The framework is benchmarked against the collisional N-body code PeTar, showing good agreement for equal-mass and realistic-IMF clusters in terms of core collapse timescales, radial density profiles, and velocity distributions. The method is then applied to a self-consistent simulation of a 10^5 M_sun giant molecular cloud (GMC) coupled with the Rigel star formation and feedback model. The collisional run and a softened-gravity version show similar star formation histories and efficiencies, but the collisional run produces a substantial population of run away and walk away stars (v > 30 km/s) with a tail ex tending beyond 100 km/s, as well as a more extended distribution of heated gas due to the redistribution of massive stars.","tokens_in":32265,"tokens_out":2311,"duration_ms":21446,"significance":"If the central result holds, Arepo-CD provides a practical and efficient tool for studying the coupled evolution of stars, gas, and feedback in star-cluster formation, addressing a key limitation of current hyrodynamic codes that rely on softened gravity. The benchmark tests againt an independent collisional code (PeTar) lend credence to the method. The application to a GMC demonstrate that close enounters during the embedded phase can generate ob servable runaway populations and alter the spatial distribution of feedback, which has implications for understanding star cluster dynamics and the impact of massive stars on their environment.","major_comments":[{"comment":"The group-identification threshold Δt_group (typically set to 10 yr) is a user-defined parameter that determines which stellar pairs receive specialized few-body integration. The paper do not present any sensitivity test for this parameter. The benchmark tests (Sect. 3) use a single value, and the GMC application (Sect. 4) also uses a single value. While the results agree with PeTar for that choice, the robustness of the method to larger or smaller Δt_group is unknown. For example, a too large Δt_group could include many stars in groups and degrade efficiency; a too small value could miss important close enounters and reduce accurancy. The authors should either perfor m a sensitivity study varying Δt_group (e.g., by a factor of 2 or 5) and show that core collapse, density profiles, and the runaway population remain consistent, or provide a strong phyical argument why the results are insensitive to the pre cise value. As this is a load-bearing assumption, a major revision is required.","section":"Sect. 2.4.1, Eq. (24) and (25)","authors_acknowledge":null,"authors_acknowledge_p":null},{"comment":"The benchmark tests are performed only for a single set of group-identification parameters (including Δt_group, r_search, and limiting sizes). Since the central claim is that Arepo-CD reproduces collisional N-body results, the dependence of this agreement on the user-defined parameters should be established. I recommend mu ltiplying the benchmark runs with at least one alternative setting (e.g., Δt_group = 5 yr and 20 yr) to demonstrate that the match to PeTar is not a coincidental result of that particular parameter choice.","section":"Sect. 3.1 and 3.2","authors_acknowledge":null,"authors_acknowledge_p":null}],"minor_comments":[{"comment":"The phrasing 'different in their treatments on collisional gravity' is awkward. Consider revising to 'differ in their treatment of collisional gravity' or 'differ in how they treat collisional gravity'.","section":"Sect. 1 (Introduction), line 'treatments on collisional gravity'","authors_acknowledge":null,"authors_acknowledge_p":null},{"comment":"The parameter r_search is introduced but no typical value is given until Section 4.1 (0.01 pc). It would be helpful to mention the typical scale (e.g., ~0.01 pc) already in the methods section to give the reader immediate context.","section":"Sect. 2.4.1","authors_acknowledge":null,"authors_acknowledge_p":null},{"comment":"The subscript 'L' in the kick operators K^{(G_k)}_L is used both for the long-range kicks of the outer splitting and for internal kicks within the group. While the meaning is clear from context, consider using a different subscript (e.g., K^{(G_k)}_{\rm int}) for clarity to avoid potential confusion.","section":"Sect. 2.2.4, Eq. (19)","authors_acknowledge":null,"authors_acknowledge_p":null},{"comment":"The y-axis label in Figure 5 is given as 'Number of isolated stars N_isol' but the plot appears to show a histogram of velocities; please verify that the label is consistent with the data (count per bin).","section":"Table 3 (not present in this draft) or Fig. 5 caption","authors_acknowledge":null,"authors_acknowledge_p":null},{"comment":"The statement that the fiducial GMC simulation requires ~3000 CPU hours is useful, but the paper also notes that a systematic scaling analysis remains to be performed. This is fine. However, to help the reader assess operational cost, please also indicate the approximate wall-clock time and number of cores used.","section":"Sect. 5.1","authors_acknowledge":null,"authors_acknowledge_p":null}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek/deepseek-v4-flash","letter":"This is a solid paper. The authors couple collisional stellar dynamics to the moving-mesh MHD code Arepo, using a three-level Hamiltonian splitting that combines leapfrog, a fourth-order forward symplectic integrator, and slow-down algorithmic regularization. That combination is new in the Arepo ecosystem, and the description of the method is detailed and readable. The benchmark tests against PeTar are convincing: core collapse times, Lagrangian radii, density profiles, and velocity distributions all agree within realization scatter. The GMC application, while exploratory, shows a genuine difference in the high-velocity tail of escapers that simply does not appear in softened-gravity runs. That is a physically interesting result and demonstrates the code's value.\n\nThe weak point is the group-identification threshold Δt_group (~10 yr). The paper does not test sensitivity to this parameter, either in the N-body benchmarks or in the GMC simulation. The benchmark clusters are simple (isolated, no gas, no binaries), so the group assignment there is trivial. In the GMC run, the choice of Δt_group directly controls which stars get few-body integration and which stay in the tree. A different value could shift the boundary and affect the runaway fraction. The authors list Δt_group as a free parameter but do not vary it. This is a limitation, but not a fatal one for a methods paper—it is common to present a first implementation with default choices. Still, a cautious reader will want to see a sensitivity test before trusting the quantitative predictions of the runaway tail.\n\nOther minor points: the code is not public; the cost estimate (~3000 CPU-hours) is given without a scaling analysis; the GMC simulation is a single realization, and the authors note that multiple realizations are important. These are acknowledged or can be deferred.\n\nWho this is for: anyone doing star cluster formation simulations with hydrodynamics. The method is practical and the paper compares well to existing frameworks (Torch, Asura+Bridge, Enzo-N). The scientific results on runaway stars and feedback redistribution are suggestive, not definitive, but they motivate further work.\n\nThe paper deserves peer review and will, after minor revisions (tone down claims about the runaway fraction, add a paragraph on parameter sensitivity or at least flag it as future work), be a useful contribution.","headline":"A well-executed hybrid collisional N-body + MHD method with clear benchmarks and a plausible GMC application, though the parameter sensitivity of the group-identification scheme is not tested.","tokens_in":32734,"tokens_out":1129,"would_cite":true,"duration_ms":12563,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek/deepseek-v4-flash","headline":"Arepo-CD couples collisional N-body dynamics with moving-mesh MHD to resolve close stellar encounters in star cluster formation simulations.","keywords":["star cluster formation","collisional N-body dynamics","magnetohydrodynamics","hybrid integration scheme","runaway stars","giant molecular cloud","stellar feedback"],"falsifier":"Run a $10^{5}$ M⊙ GMC simulation identical to the fiducial model but vary Δt_group from 1 yr to 100 yr; if the runaway star fraction changes by less than the stochastic scatter between realizations, the parameter choice is robust. Conversely, a strong dependence would show that the results are not converged.","tokens_in":33327,"feed_emoji":"🌟","tokens_out":1923,"duration_ms":20845,"temperature":0.7,"pith_summary":"The paper presents Arepo-CD, a hybrid framework that embeds collisional N-body dynamics into the moving-mesh magnetohydrodynamics code Arepo. It treats star-star gravity without softening while keeping gas and dark matter interactions softened, allowing accurate modeling of close stellar encounters, binary formation, and dynamical ejections within self-consistent star formation simulations. By identifying compact stellar subsystems on-the-fly and evolving them with a fourth-order forward symplectic integrator and slow-down algorithmic regularization, Arepo-CD achieves high accuracy for few-body dynamics without imposing short timesteps on the whole simulation. Benchmark tests against the collisional N-body code PeTar show agreement in core collapse, radial density profiles, and velocity distributions. In a giant molecular cloud simulation, Arepo-CD produces a substantial population of runaway stars with speeds exceeding 100 km/s and redistributes radiative feedback into the surrounding lower-density gas, effects that softened-gravity simulations miss.","feed_headline":"Hybrid integrator catches runaway stars in star cluster formation","feed_subtitle":"Arepo-CD couples collisional stellar dynamics to MHD, producing a high-velocity tail missing in softened-gravity simulations.","key_machinery":"The key machinery is the hybrid Hamiltonian-splitting scheme based on the particle-particle particle-tree (P3T) method, which separates the Hamiltonian into short-range and long-range components. Stellar particles in compact groups (identified by free-fall and fly-by timescales below a user-defined threshold, typically 10 yr) are integrated internally with a fourth-order forward symplectic integrator, with the tightest subgroups handled by slow-down algorithmic regularization. The external gravitational influence from particles outside the group is applied as a long-range kick via the tree gravity solver, keeping the global timestep hierarchy efficient.","core_discovery":"Arepo-CD demonstrates that a three-level hybrid integration scheme (second-order leapfrog for large-scale evolution, fourth-order forward symplectic integrator for compact stellar groups, and slow-down algorithmic regularization for the stiffest interactions) can efficiently and accurately couple collisional stellar dynamics to moving-mesh MHD star formation simulations. The framework reproduces core collapse and post-collapse expansion in isolated cluster tests, and when applied to a self-consistent GMC simulation, close few-body encounters generate runaway stars and alter the spatial distribution of feedback, showing that collisional dynamics matters during the embedded phase of cluster formation.","pith_inferences":["The paper's group identification parameters (Δt_group = 10 yr, N_FSI_max ∼ 100, N_SDAR_max ≤ 4) are not tested for sensitivity, so the accuracy and cost of the method may depend on these choices in regimes with different densities or binary fractions.","Because Arepo-CD treats all stellar particles without a mass threshold, it could be applied to lower-mass star formation simulations where interactions among low-mass stars also matter, not only massive stars.","The three-level hierarchy suggests that similar hybrid integrator designs could be adopted in other astrophysical codes that need to span wide dynamical ranges, such as simulations of galactic centers or dense stellar systems with black holes."],"forward_implications":["If Arepo-CD is correct, simulations of star cluster formation that use softened gravity systematically underestimate the production of runaway and walkaway stars during the embedded phase.","Ejected massive stars carry ionizing radiation, winds, and supernovae into lower-density gas, which may increase the escape fraction of ionizing photons from molecular clouds and affect galaxy-scale feedback.","The hybrid integration scheme makes it computationally practical (∼3000 CPU hours for a 10^5 M⊙ cloud) to run multiple realizations and parameter surveys for studying cluster formation statistics.","Arepo-CD provides a foundation for future extensions such as stellar collisions, primordial binaries, and binary stellar evolution, enabling predictive studies of compact-object populations in young clusters."],"supporting_citations":[{"why":"Supplies the slow-down algorithmic regularization method (SDAR) used for the tightest few-body subsystems.","marker":"L. Wang et al. 2020b"},{"why":"Provide the fourth-order forward symplectic integrator (FSI) formalism adopted for compact stellar groups.","marker":"A. Rantala et al. 2021, 2023"},{"why":"Describe the Arepo moving-mesh MHD code and its tree-gravity solver, on which Arepo-CD is built.","marker":"V. Springel 2010; R. Weinberger et al. 2020"},{"why":"Describe the Rigel star formation and feedback model that Arepo-CD couples to for the GMC simulations.","marker":"Y. Deng et al. 2024b, 2026b"},{"why":"Documents the PeTar collisional N-body code used as the benchmark for validation tests.","marker":"L. Wang et al. 2020a"}],"fun_headline_variants":["Hybrid code reveals runaway stars in star formation","Runaway stars emerge in coupled MHD-star dynamics","New integrator tracks high-velocity stellar ejections","Collisional stellar dynamics reshape cluster feedback","Unsoftened gravity yields stellar speed tail above 100 km/s"],"cache_read_input_tokens":0,"weakest_assumption_plain":"The group identification relies on user-defined parameters (free-fall timescale threshold Δt_group = 10 yr and maximum group sizes N_FSI_max ≈ 100, N_SDAR_max ≤ 4) whose impact on accuracy and computational cost has not been tested through a sensitivity study.","fun_headline_variants_meta":{"raw":{"variants":["Hybrid code reveals runaway stars in star formation","Runaway stars emerge in coupled MHD-star dynamics","New integrator tracks high-velocity stellar ejections","Collisional stellar dynamics reshape cluster feedback","Unsoftened gravity yields stellar speed tail above 100 km/s"]},"model":"deepseek/deepseek-v4-flash","effort":"low","cost_usd":0.000104,"raw_usage":{"cost":0.000104234592,"is_byok":false,"cost_details":{"upstream_inference_cost":0.000104234592,"upstream_inference_prompt_cost":6.9889512e-05,"upstream_inference_completions_cost":3.434508e-05},"total_tokens":1163,"prompt_tokens":974,"completion_tokens":189,"prompt_tokens_details":{"audio_tokens":0,"video_tokens":0,"cached_tokens":256,"cache_write_tokens":0},"completion_tokens_details":{"audio_tokens":0,"image_tokens":0,"reasoning_tokens":113}},"tokens_in":974,"tokens_out":189,"cost_source":"openrouter_usage","duration_ms":2520,"image_pages":[],"temperature":1.0,"billed_attempts":1,"reasoning_tokens":113,"transport_provider":"openrouter","cache_read_input_tokens":0,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-10-22T14:50:38.218589+00:00","model_set":{"reader":"deepseek/deepseek-v4-flash"},"falsifier":"Run a $10^{5}$ M⊙ GMC simulation identical to the fiducial model but vary Δt_group from 1 yr to 100 yr; if the runaway star fraction changes by less than the stochastic scatter between realizations, the parameter choice is robust. Conversely, a strong dependence would show that the results are not converged.","supporting_citations":[{"cited_title":"2021, MNRAS, 502, 5546, doi: 10.1093/mnras/stab057","cited_arxiv_id":null,"evidence_quote":"Provide the fourth-order forward symplectic integrator (FSI) formalism adopted for compact stellar groups."},{"cited_title":"Monthly Notices of the Royal Astronomical Society , author =","cited_arxiv_id":null,"evidence_quote":"Describe the Arepo moving-mesh MHD code and its tree-gravity solver, on which Arepo-CD is built."}],"review_version":1}