{"id":"f12effb6-5f96-4da1-906b-9487ea795e6b","arxiv_id":"2501.00846","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A thread-safe, high-order lattice Boltzmann model with Allen-Cahn interface tracking is shown to reproduce bubble rising and droplet collision dynamics at density and viscosity ratios up to 1000 and 100.","lead":"This paper extends a thread-safe lattice Boltzmann method to simulate two-phase flows with large density and viscosity differences, like water and air, and tests it on rising bubbles, colliding droplets, and raindrops hitting a surface. It matters because the method is designed to run efficiently on GPU supercomputers, which could make realistic multiphase simulations cheaper without machine learning shortcuts.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The printed Laplacian stencil in Eq. 19 does not annihilate constant fields, so the method as written has a spurious source in every interface force.","rationale":"The reader's weakest assumption was the force splitting in Eqs. 16-17. I agree that the balance of pressure and viscous forces is important, but I found a more elementary and easily checked defect upstream of it: the discrete Laplacian used to form the chemical potential and the surface-tension force is not a Laplacian as printed. This affects both the benchmark validations and the claimed generality at density ratio 1000. The concern is not that the method disagrees with the literature; the graphical comparisons in Figs. 4-13 are plausible evidence if the code implements the equations correctly. It is that the equations in the paper cannot be the equations the code runs unless Eq. 19 is a typographical error. Because this is fixable and does not contradict the benchmark agreement, it raises the same need for revision the reader already expressed; I therefore leave the conditional verdict unchanged.","tokens_in":13345,"tokens_out":8966,"duration_ms":89085,"concrete_test":"Apply Eq. 19 to a constant scalar field Ψ=1 on a D3Q27 lattice. If the result is nonzero, the printed stencil is invalid. Then re-run a static droplet test at density ratio 1000 with an isotropic Laplacian that annihilates constants and check that spurious currents and interface drift converge to zero with grid refinement; if they do not, the high-density-contrast force balance remains unverified.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The most load-bearing issue is not the force balance itself but the discrete derivative operators it relies on. Section II.C states that gradients and Laplacians are computed with Eqs. 18-19, and Eq. 19 is used for the chemical potential's second derivative and for the Allen-Cahn diffusion term in Eq. 20. As printed, Eq. 19 is not a second-derivative operator: for a constant field Ψ=C, the sum over i≠0 of wi equals 1-w0, so the right-hand side equals (1-2w0)C/c_s^2, which is (11/27)C/c_s^2 for D3Q27 (w0=8/27), not zero. A valid isotropic Laplacian must annihilate constants; the standard form is (2/c_s^2)(Σ_i wi Ψ(x+ci) - Ψ(x)) or an equivalent operator. Consequently, every simulation run with the printed stencil contains a spurious constant-field source in the chemical potential and in the Allen-Cahn diffusion term, independent of density ratio. This is an internal inconsistency in the method description: the central high-density-contrast claim is unverifiable from the paper as written unless the code silently substitutes a different, correct operator.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript presents a hybrid lattice Boltzmann method for two-component flows with high density and viscosity contrasts, combining a thread-safe, high-order regularized collision-streaming scheme with a conservative Allen-Cahn interface-tracking equation solved by finite differences. The surface tension, pressure-gradient, and viscous forces are added as external forcings, with the pressure and viscous contributions intended to extend the method to density and viscosity ratios up to 1000 and 100. Validation is carried out against reference simulations for rising bubbles and against experiments for head-on and off-axis droplet collisions, and GPU strong/weak scaling is reported on the Leonardo supercomputer. An application to successive raindrop impacts on a solid surface is also presented.","tokens_in":13581,"tokens_out":5293,"duration_ms":53333,"significance":"If the method performs as claimed, it would be a useful contribution to GPU-oriented multiphase LBM, since the thread-safe formulation removes race conditions while retaining third-order Hermite reconstruction. The paper has concrete strengths: it benchmarks against external experiments and independent simulations rather than only self-comparisons, it demonstrates the challenging air-water density ratio, and it reports scaling data on a modern GPU cluster. However, the manuscript as written contains an inconsistency in the discrete Laplacian operator used in several load-bearing equations, and the validation is largely qualitative, with no grid-convergence study and at least one visible discrepancy with the cited experiment. Because the central claim depends on these points, the paper cannot be accepted in its present form.","major_comments":[{"comment":"The discrete second-derivative operator printed in Eq. (19) does not annihilate constant fields. For a constant field Ψ=C, the right-hand side evaluates to (1-2w0)C/c_s^2, which is (11/27)C/c_s^2 for the D3Q27 weights (w0=8/27), rather than zero. Since this operator is used for the chemical potential in Eq. (3) and for the Allen-Cahn diffusion term in Eq. (20), the method as written contains a spurious constant-field source in every time step, independent of density ratio. The standard isotropic Laplacian has the form (2/c_s^2)(Σ_i w_i Ψ(x+c_i)-Ψ(x)) or an equivalent expression that annihilates constants. The authors should either correct Eq. (19) or explicitly define the operator actually implemented, and they should show that this operator passes the constant-annihilation test. This is load-bearing because the high-density-contrast claim relies on the discrete forces being balanced at the interface.","section":"II.C, Eq. (19)"},{"comment":"The validation of the rising-bubble case is qualitative: the center-of-mass velocity comparison with Adelsberger et al. is presented graphically with no error metrics, and the statement that 'good agreement' holds 'provided that the interface width is adequately resolved' is not supported by a grid-convergence or interface-width-convergence study. Given that the main novelty is the ability to handle density ratio 1000 and viscosity ratio 100, the absence of a resolution study is a significant gap: it leaves open whether the observed agreement arises from correct force balance or from sufficiently diffuse interfaces masking errors.","section":"III.A, Fig. 4 and surrounding text"},{"comment":"The paper acknowledges that the simulation produces a satellite droplet after thread breakup at We=25, while the cited experiments of Ashgriz and Poo do not report one. The explanation invoking Huang et al. is plausible, but it is not accompanied by a quantitative comparison (satellite size, pinch-off time, or a resolution study showing that this satellite is converged). In addition, Fig. 7 attributes a thin-film breakup at We=40 to insufficient resolution. Together these points undermine the abstract's claim of results 'in agreement with experiments' and need to be addressed with either a systematic resolution/parameter study or a more nuanced statement of agreement.","section":"III.B, We=25 case"},{"comment":"The force-splitting construction assumes that the discrete pressure-gradient force F_p and viscous force F_ν, together with the surface tension force, reproduce the target momentum equation at density ratio 1000. No test is reported for spurious currents at a stationary droplet or for the discrete balance of these forces at a sharp interface. Such a test is directly relevant because the discrete stencils of Eqs. (18)-(19) determine whether the pressure-gradient and viscous forces are accurately represented. The authors should add a quantitative assessment, e.g., maximum spurious velocity as a function of grid resolution and interface width for a stationary droplet at high density ratio.","section":"II.C, Eqs. (15)-(17)"}],"minor_comments":[{"comment":"The left-hand side of Eq. (19) has two indices (∂α∂βΨ) while the right-hand side is a scalar; please clarify whether this equation is intended as the Laplacian or as the full Hessian, and supply the correct lattice weights and central term.","section":"II.C, Eq. (19)"},{"comment":"The caption of Fig. 5 states 'We=23' whereas the text discussing this case refers to We=25; the inconsistency should be reconciled.","section":"III.B, Fig. 5 caption"},{"comment":"The text spells the name 'Jaqmin' but the reference [22] is by Jacqmin; please correct the spelling.","section":"II.C, Eq. (3)"},{"comment":"The paragraph claiming that machine learning is 'non-essential' for studying complex flows is an opinion not supported by the scaling benchmarks or by any comparison with ML-based methods; it should be removed or supported with quantitative evidence.","section":"II.E, Performance section"},{"comment":"The collision map reports only the numerical simulation outcomes as triangles; overlaying the experimental regime boundaries from [34] would make the claimed agreement quantitatively assessable.","section":"III.B, Fig. 13"},{"comment":"The equilibrium profile in Eq. (5) uses absolute values of coordinate differences, which is not the standard planar interface profile and is introduced without derivation; please clarify this choice and its relation to the coordinate axes.","section":"II.A, Eq. (5)"}],"recommendation":"major_revision","confidential_remarks":"The concern in Eq. (19) may well be a typesetting error rather than an error in the code, since the authors cite Thampi et al. for isotropic discrete Laplacians and the surrounding formulas otherwise follow standard practice. However, as printed, the manuscript is internally inconsistent and the central validation claims cannot be checked. I would ask the authors to supply the corrected operator, a brief verification that it annihilates constants, and a quantitative resolution study for at least one high-density-ratio case before reconsidering the paper. The novelty relative to the authors' earlier TSLB papers is incremental but real, and the droplet-collision benchmarks are appropriate for the journal's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: the paper extends the thread-safe LBM framework to air-water-like density and viscosity contrasts, and the benchmark results are mostly consistent with experiments and prior simulations. It also gives concrete GPU scaling numbers, which is more than many LBM papers do. But there is a load-bearing problem in the printed discrete operators. Eq. 19, used for the chemical potential's second derivative and the Allen-Cahn diffusion term, is not a Laplacian. For a constant field it returns (11/27)C/c_s^2 on D3Q27, not zero. That means every simulation described in the paper, if run with the printed stencil, contains a spurious constant-field source in the interface forces. The code almost certainly substitutes the correct isotropic stencil, but a reader cannot verify the central claim from the manuscript. That is a reproducibility issue, not a minor typo.\n\nWhat is genuinely new: the thread-safe, Hermite-regularized collision scheme has been applied to single-phase and unit-density multicomponent flows before; here it is adapted to density ratio 1000 and viscosity ratio 100, with a force splitting that looks physically reasonable. The rising-bubble regime map, droplet collision outcomes, and raindrop impact sequences are visually in line with the references. The A100 strong and weak scaling tables are honest and useful.\n\nSoft spots beyond Eq. 19: validation is largely qualitative. There is no grid-convergence study and no error metric against the benchmark data. The We=25 collision produces a satellite droplet that the original 1990 experiment did not report; the authors attribute it to resolution and cite a newer experiment seeing the same, which is fair, but it remains a discrepancy in a central test case. Also, the Allen-Cahn mobility D and the coefficients κ and β are defined in two places with different notation, and it takes work to see which definition is used in the code. Minor, but worth cleaning up.\n\nThe paper is for people building GPU multiphase LBM solvers. They will get real value from the force-splitting details and the performance data. It deserves a serious referee, not a desk reject. The core idea is sound and the benchmarks mostly support it, but the stencil must be corrected and a convergence study added before it can be trusted as a reproducible method.","headline":"A genuinely useful extension of thread-safe LBM to high density contrasts, but the printed Laplacian stencil in Eq. 19 does not annihilate constants, so the method as written is not reproducible until that is fixed.","tokens_in":14098,"tokens_out":1861,"would_cite":false,"duration_ms":19545,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["76M28","76T10"],"pacs":["47.11.-j","47.55.D-"],"model":"deepseek-v4-flash","headline":"A thread-safe, third-order lattice Boltzmann method reproduces droplet and bubble dynamics at density ratios up to 1000 and viscosity ratios up to 100.","keywords":["lattice Boltzmann method","thread-safe GPU computing","Allen-Cahn phase field","high density contrast multiphase flow","rising bubble","droplet collision","regularized lattice Boltzmann","air-water simulation"],"falsifier":"A static planar interface with $\\rho_L/\\rho_G=1000$ and no body force should remain perfectly quiescent: if the spurious velocity produced by the unbalanced discrete pressure and viscous force terms does not vanish, or does not decrease as the interface width is refined, the high-density-ratio capability claim is false. The same test with $\\nu_L/\\nu_G=100$ isolates the viscous correction.","tokens_in":13181,"feed_emoji":"💧","tokens_out":13281,"duration_ms":115242,"temperature":0.7,"pith_summary":"This paper claims that a thread-safe lattice Boltzmann implementation, which reconstructs distribution populations from macroscopic fields to avoid race conditions, can be extended to two-fluid flows with density ratios up to 1000 and viscosity ratios up to 100 by adding an Allen-Cahn interface tracker and three extra force terms. The resulting solver reproduces the rising-bubble deformation regimes on the Eötvös–Galilei map for an air-water-like contrast and matches experimental sequences of head-on and off-axis water-droplet collisions, including satellite formation. It also reports weak-scaling efficiency near 0.95 up to 32 GPUs, so high-contrast multiphase simulations become practical on shared-memory and GPU architectures. A sympathetic reader would care because this pushes diffuse-interface lattice Boltzmann methods into the parameter range of realistic environmental, microfluidic, and biological flows.","feed_headline":"Lattice Boltzmann now handles 1000:1 density and 100:1 viscosity","feed_subtitle":"Thread-safe GPU solver matches droplet and bubble experiments at air-water density contrast.","key_machinery":"The central object is the fused streaming-collision update $f_i(\\mathbf{x}+\\mathbf{c}_i\\Delta t,t+\\Delta t)=f_i^{eq}(\\mathbf{x},t)+(1-\\omega)f_i^{neq}(\\mathbf{x},t)+S_i(\\mathbf{x},t)$, in which both the equilibrium and non-equilibrium parts are reconstructed from macroscopic fields at each node rather than read from memory, using the recursive properties of Hermite polynomials up to third order on a D3Q27 (27-velocity cubic) stencil. This reconstruction is what removes race conditions, letting streaming and collision be fused into one kernel. The density and viscosity contrast is then injected through the three force terms, while the conservative Allen-Cahn equation $\\partial_t\\phi+u_\\alpha\\partial_\\alpha\\phi=D\\partial_\\alpha\\partial_\\alpha\\phi-\\kappa\\partial_\\alpha(\\phi(1-\\phi)n_\\alpha)$ is advanced by a forward-time centred-space finite-difference scheme to track the interface.","core_discovery":"The paper's central claim is that the thread-safe lattice Boltzmann update—rebuilding post-collision populations from macroscopic moments instead of reading neighbouring populations—can carry a high-contrast multiphase model if the non-equilibrium part is reconstructed to third order in Hermite polynomials and the density/viscosity contrast is injected through three explicit forces: surface tension $F_s=\\mu_\\phi\\partial_\\alpha\\phi$, pressure correction $F_p=-p^*c_s^2\\partial_\\alpha\\rho$, and viscous correction $F_\\nu=-\\frac{\\nu\\omega}{c_s^2\\Delta t}[\\sum_i(f_i-f_i^{eq})c_{i\\alpha}c_{i\\beta}]\\partial_\\alpha\\rho$. With this construction, the paper demonstrates that $\\rho_L/\\rho_G=1000$ and $\\nu_L/\\nu_G=100$ are reachable on a 27-velocity cubic lattice, recovering the rising-bubble regime map and reproducing experimental water-droplet collision outcomes at $\\mathit{Oh}=0.0044$.","pith_inferences":["The high-density-ratio capability ultimately depends on the balance of the pressure and viscous force corrections at the sharpest interface cells; a direct test the paper does not report is to compute the discrete momentum balance across a static interface at density ratio 1000 before running any dynamic simulation.","The paper's assertion that machine learning is non-essential for complex multiphase flow is a broader opinion than the benchmarks prove; the measurements establish scalability and fidelity, not a comparison with learned surrogate models.","The raindrop-train simulation points toward a testable extension: quantifying how impact spacing, surface wettability, and roughness control the breakup statistics and the ejection of microplastics from the lamella rim.","The paper flags dual-grid resolution of the phase field as future work; if implemented, the interface width could be controlled independently of flow resolution, which would further relax the resolution constraint seen in the thin-film breakup case."],"forward_implications":["At density ratio 1000 and viscosity ratio 100, the method reproduces the air-water rising-bubble deformation regimes on the Eötvös–Galilei map, including axisymmetric, hat-shaped, peripheral-breakup, central-breakup, and oscillatory regimes, matching the benchmark reference when the interface is resolved.","Head-on and off-axis water-droplet collisions with Weber numbers from 10 to 96 and Reynolds numbers from roughly 700 to 2200 reproduce the experimental coalescence and separation boundaries and the formation of satellite droplets at Ohnesorge number 0.0044.","Because equilibrium and non-equilibrium populations are reconstructed from macroscopic fields, the fused stream-and-collision kernel is free of read-after-write race conditions, which is what permits weak-scaling efficiency near 0.95 up to 32 GPUs.","A train of five raindrops impacting a solid substrate reproduces the observed lamella spreading, recoil, and growth of the liquid film with successive impacts, giving a direct numerical route to studying microplastic and pollutant dispersion by rain.","Third-order recursive regularization raises the stability ceiling of the thread-safe scheme, so a droplet diameter of 100 lattice units with a Cahn number of 0.04 is enough to capture reflexive separation and satellite droplets."],"supporting_citations":[{"why":"Supplies the base thread-safe fused streaming-collision scheme that the multiphase extension builds on.","marker":"[18]"},{"why":"Provides the earlier hybrid thread-safe lattice Boltzmann model with interface capturing at unit density ratio and the accLB code base.","marker":"[20]"},{"why":"Supplies the recursive Hermite regularization used to reconstruct non-equilibrium moments up to third order.","marker":"[21]"},{"why":"Provides the forcing scheme that encodes the trapezoidal rule in the collision step.","marker":"[29]"},{"why":"Relates the second-order moment of the distribution to the deviatoric stress, used to build the viscous correction.","marker":"[30]"},{"why":"Supplies the isotropic discrete gradient and Laplacian stencils used for the force terms and interface tracking.","marker":"[31]"},{"why":"Classifies the Eötvös–Galilei deformation regimes of a rising bubble used for comparison.","marker":"[32]"},{"why":"Provides the benchmark rising-bubble center-of-mass velocity data at density ratios 10 and 1000.","marker":"[33]"},{"why":"Supplies the experimental droplet-collision sequences and regime data for head-on and off-axis comparisons.","marker":"[34]"},{"why":"Provides the experimental comparison for a successive train of droplets impacting a solid surface.","marker":"[38]"}],"fun_headline_variants":["Thread-safe LBM handles 1000:1 density and 100:1 viscosity","GPU-ready LBM matches droplet experiments at 1000:1 density","High-contrast multiphase LBM validated on bubble and droplet tests","Thread-safe solver reaches extreme density and viscosity contrasts"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The method's high-density-ratio reach rests on the assumption that the pressure-gradient and viscous force corrections (Eqs. 16–17) exactly recover the target momentum equation at the discrete interface, using the stencils of Eqs. 18–19; if those forces are not balanced at the sharpest interface cells, the claim of density ratio 1000 collapses even though the thread-safe streaming itself is sound.","fun_headline_variants_meta":{"raw":{"variants":["Thread-safe LBM handles 1000:1 density and 100:1 viscosity","GPU-ready LBM matches droplet experiments at 1000:1 density","High-contrast multiphase LBM validated on bubble and droplet tests","Thread-safe solver reaches extreme density and viscosity contrasts"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000605,"raw_usage":{"total_tokens":2787,"prompt_tokens":874,"completion_tokens":1913,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":490,"completion_tokens_details":{"reasoning_tokens":1835}},"tokens_in":490,"tokens_out":1913,"duration_ms":13087,"temperature":1.0,"reasoning_tokens":1835,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T22:41:19.452550+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"A static planar interface with $\\rho_L/\\rho_G=1000$ and no body force should remain perfectly quiescent: if the spurious velocity produced by the unbalanced discrete pressure and viscous force terms does not vanish, or does not decrease as the interface width is refined, the high-density-ratio capability claim is false. The same test with $\\nu_L/\\nu_G=100$ isolates the viscous correction.","supporting_citations":[{"cited_title":"Montessori, G","cited_arxiv_id":null,"evidence_quote":"Supplies the base thread-safe fused streaming-collision scheme that the multiphase extension builds on."},{"cited_title":"Montessori, M","cited_arxiv_id":null,"evidence_quote":"Provides the earlier hybrid thread-safe lattice Boltzmann model with interface capturing at unit density ratio and the accLB code base."},{"cited_title":"A high-performance lattice Boltzmann model for multicomponent turbulent jet simulations","cited_arxiv_id":"2403.15773","evidence_quote":"Supplies the recursive Hermite regularization used to reconstruct non-equilibrium moments up to third order."},{"cited_title":"Grad, Note on n-dimensional hermite polynomials, Communications on Pure and Applied Mathematics 2 (4) (1949) 325–330","cited_arxiv_id":null,"evidence_quote":"Provides the forcing scheme that encodes the trapezoidal rule in the collision step."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Relates the second-order moment of the distribution to the deviatoric stress, used to build the viscous correction."},{"cited_title":"Kr¨ uger, F","cited_arxiv_id":null,"evidence_quote":"Supplies the isotropic discrete gradient and Laplacian stencils used for the force terms and interface tracking."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Classifies the Eötvös–Galilei deformation regimes of a rising bubble used for comparison."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the benchmark rising-bubble center-of-mass velocity data at density ratios 10 and 1000."},{"cited_title":"Adelsberger, P","cited_arxiv_id":null,"evidence_quote":"Supplies the experimental droplet-collision sequences and regime data for head-on and off-axis comparisons."},{"cited_title":"Amani, N","cited_arxiv_id":null,"evidence_quote":"Provides the experimental comparison for a successive train of droplets impacting a solid surface."}],"review_version":1}