{"id":"53615da1-a694-4496-b1b2-5859ab37ed28","arxiv_id":"2507.14076","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A multi-Gaussian Wigner-function variational method simulates open bosonic dynamics and finds Ising-like critical slowing down in a driven-dissipative Bose-Hubbard lattice.","lead":"This paper introduces a variational method that approximates quantum dynamics of many interacting bosons using a sum of Gaussian components in Wigner phase space, with equations of motion computed by automatic differentiation. It demonstrates the method on driven-dissipative Bose-Hubbard lattices up to 12x12 sites and reports critical slowing down consistent with the 2D quantum Ising universality class.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The many-body gap result lacks convergence checks: finite-time VMG trajectories with NG=16 are used to extract a Liouvillian gap without verifying convergence in Gaussian number or integration time at the critical point.","rationale":"The reader's conditional verdict identifies the same weak point, and I agree. The variational method is elegant, and the single-mode benchmarks (Figs. 2 and 3) provide genuine evidence that the ansatz can represent strong Wigner negativity with controlled convergence. The steady-state parity comparison with CSR (Fig. 5a) independently validates the static variational manifold for 6×6. However, the dynamical critical claim is a much stronger statement: it requires that the variational flow, not just the steady state, is accurate on the 12×12 lattice in the regime where the true gap is smallest. This is precisely where any variational approximation is most fragile: near criticality the slow mode is a collective fluctuation whose Wigner function is highly non-Gaussian, and the VMG trajectory may miss it entirely. The paper provides no many-body convergence test in NG, no error bars on λ, and no description of the gap-extraction fit. The fact that β, ν, and z are imported from the 2D quantum Ising model and Gc is read off the same finite-size data means the collapse in Fig. 7(c) cannot by itself demonstrate that the method extracts those exponents; it shows consistency under the assumption that they apply. The proposed check—doubling NG and T at the critical point and verifying that λ is stable—directly targets the load-bearing assumption. If it passes, the conditional should be lifted; if it fails, the critical-slowing-down claim and the statement that the paper provides the first calculation of the dynamical critical exponent would need to be substantially weakened. The verdict remains CONDITIONAL, matching the reader's assessment.","tokens_in":27906,"tokens_out":3709,"duration_ms":499727,"concrete_test":"Repeat the 12×12 lattice simulation at G=Gc≈1.0γ with NG=32 (and, if feasible, NG=64) and extend the integration time to T=240γ⁻¹; extract λ from two independent time windows (e.g., [40,80]γ⁻¹ and [120,240]γ⁻¹) and check that the inferred λ stabilizes with respect to both NG and T. If λ shifts by more than the estimated uncertainty, or if the late-time decay is not a clean single exponential, the reported gap and the Fig. 7(c) collapse are not converged and the critical-slowing-down claim is unsupported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central dynamical claim rests on the unverified assumption that the VMG trajectory (Eq. 8, T·dθ/dt = V) on the multi-Gaussian manifold (Eq. 5) reproduces the true late-time decay of the 12×12 lattice with NG=16 Gaussians, and that an exponential extracted from trajectories of duration up to 120γ⁻¹ (Sec. VI) equals the asymptotic Liouvillian gap. Section V demonstrates exponential convergence with NG for the single-mode Kerr oscillator (Fig. 3c), but no analogous NG-convergence is shown for the lattice, where the variational manifold is far more constrained (block-diagonal covariance, only single-mode squeezing per Appendix D). The gap-extraction procedure is not described: no fit window, no fitting function, no uncertainty estimates. At G=Gc, Fig. 7(a) shows parity still far from steady state at the end of the trajectory, so the 'asymptotic' rate is inferred from a transient. Additionally, the collapse in Fig. 7(c) imports β, ν, and z from the literature and uses Gc ≈ 1.0γ estimated from the same finite-size data, making it a consistency check under assumed exponents rather than an independent extraction. If the variational projection or the finite-time fit misses the slow mode that closes the gap, the reported λ and the resulting z-collapse would reflect the approximation, not the Liouvillian.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper introduces a variational method for open quantum bosonic systems based on a Wigner phase-space representation. The ansatz, Eq. (5), is a sum of complex Gaussian components, and the time evolution of its parameters is obtained from the Dirac-Frenkel principle, Eq. (8), with the Wigner quantum geometric tensor and Liouvillian gradient evaluated analytically through Gaussian moments and automatic differentiation. The method is benchmarked on a single-mode driven-dissipative Kerr parametric oscillator against exact diagonalization, including Wigner-function snapshots and exponential convergence with the number of Gaussians (Section V). It is then applied to two-dimensional Bose-Hubbard lattices with two-boson driving and losses (Section VI): the steady-state parity is compared with corner-space renormalization for a 6x6 lattice and computed for lattices up to 12x12, and the Liouvillian gap is extracted from relaxation dynamics (Section VII). The central physics claim is that the finite-size scaling of the gap reveals critical slowing down with dynamical exponents of the 2D quantum Ising universality class, using beta=0.32641871, nu=0.62997097, and z=2.0235.","tokens_in":28177,"tokens_out":5305,"duration_ms":57835,"significance":"The methodological core is valuable and largely well executed. The single-mode benchmarks in Section V, especially the exponential control of the error with the number of Gaussians shown in Fig. 3(c), are convincing and demonstrate that the ansatz can capture strong Wigner negativity. The steady-state comparison against the corner-space renormalization method in Fig. 5(a) provides an independent check for a genuinely many-body lattice. The analytical Gaussian-moment framework combined with Taylor-mode automatic differentiation is a genuine technical contribution, and the reported scalability to a 12x12 lattice with 6920 variational parameters is impressive. If the lattice dynamics and critical-scaling claims are substantiated, this would be a significant advance for open bosonic many-body simulation. The main weakness is that the central physics claim rests on a gap-extraction and scaling-collapse analysis that currently lacks convergence checks, error bars, and a clear fitting protocol; the collapse is a consistency check under externally fixed exponents rather than an independent extraction.","major_comments":[{"comment":"The extraction of the Liouvillian gap lambda is not described: no fitting function, fit window, or uncertainty estimates are given, and Fig. 7(a) shows that at G=Gc the parity has not reached its steady state within the simulation time for the larger lattices. The reported lambda is therefore inferred from a transient, and the identification of this rate with the asymptotic Liouvillian gap is unsupported. Please specify the fitting procedure, show that a single-exponential decay describes the data over a stable window, and demonstrate convergence of lambda with respect to the final integration time.","section":"Section VII, Fig. 7"},{"comment":"No convergence in the number of Gaussians NG is shown for the lattice dynamics. The exponential convergence of Fig. 3(c) is demonstrated only for the single-mode Kerr oscillator, whereas the 12x12 results use NG=16 with the restricted parametrization of Appendix D, whose covariance matrices contain only single-mode squeezing and no inter-mode blocks. Because this ansatz manifold is far more constrained in the lattice case, the authors should show, at least for a 6x6 or 8x8 lattice, that the steady-state parity and the extracted relaxation rate are stable under increasing NG, for example NG=4, 8, 16, 32.","section":"Section VI and Appendix D"},{"comment":"The finite-size collapse is a consistency check under assumed exponents rather than an independent extraction: beta, nu, and z are imported from the 2D quantum Ising literature, and Gc about 1.0 gamma is estimated from the same finite-size data that is then rescaled. The manuscript should quantify the sensitivity of the collapse to the choice of Gc and to the assumed exponents, and should report error bars on lambda and Gc; without this, the claim that the method extracts dynamical exponents of the 2D quantum Ising universality class is overstated.","section":"Section VII, Fig. 7(c)"},{"comment":"The main text states that the reader is referred to Appendix D for a detailed description of the initial conditions used for the calculated dynamics, but Appendix D contains only the variational parametrization and does not specify the initial conditions. This information is necessary to reproduce the relaxation dynamics and to assess whether the extracted long-time rate depends on the initial preparation.","section":"Section VI and Appendix D"}],"minor_comments":[{"comment":"There are several typos that should be corrected: 'systems systems' in the introduction, 'powerfull' and 'architechtures' in Section I, and 'Lindbald' in Appendix B.","section":"Section I and Appendix B"},{"comment":"The onsite sum in Eq. (B8) runs to NG rather than to M, which is inconsistent with the main-text Hamiltonian in Eq. (17); this appears to be a typo.","section":"Appendix B1, Eq. (B8)"},{"comment":"The Dirac-Frenkel equations require solving a linear system with the quantum geometric tensor T, but the manuscript does not discuss the conditioning or possible rank-deficiency of T, nor any regularization used in the numerical solution; a brief statement on how Eq. (8) is inverted in practice would aid reproducibility.","section":"Section III, Eq. (8)"},{"comment":"No data or code availability statement is provided; given the complexity of the implementation, releasing the code or at least the data underlying Figs. 5 and 7 would strengthen the paper.","section":"Section V, Fig. 3 and Section VII"},{"comment":"The phrase 'first calculation of the dynamical critical exponent' should be tempered: since the collapse uses literature exponents and an estimated Gc, the analysis is better described as a consistency check, not a standalone determination of z.","section":"Section VII"}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"This is a solid method paper with one overreach. What is genuinely new is the combination of a multi-Gaussian Wigner ansatz, analytic generalized Gaussian moments, and Taylor-mode automatic differentiation to evaluate the Dirac-Frenkel equations. The single-mode Kerr oscillator benchmark is strong: the VMG results track exact diagonalization through deep Wigner negativity, and the relative error decays exponentially with the number of Gaussians. The steady-state parity comparison against corner-space renormalization on a 6x6 lattice is an independent check, not a circular test, and it works. That part deserves credit.\n\nThe soft spot is the critical-dynamics section. The convergence that is demonstrated for the single-mode case is not shown for the lattice. The 12x12 calculation uses 16 Gaussians with block-diagonal covariance, and nothing in the paper certifies that this manifold captures the slow mode that closes the Liouvillian gap. The gap extraction procedure is also underspecified: no fit window, no fitting function, no uncertainty. Figure 7(a) shows the parity still far from steady state at the end of the trajectory at Gc, so the 'asymptotic' rate is inferred from a transient. And the collapse in Fig. 7(c) imports beta, nu, and z from the literature while estimating Gc from the same finite-size data. This is a consistency check under assumed exponents, not an independent extraction. The conclusions say 'extracted critical exponents,' which overstates what the data support.\n\nNone of this undermines the variational machinery. The method is benchmarked externally, the technical development is real, and the single-mode results are convincing. But the many-body dynamical claim is not yet load-bearing. The paper needs NG-convergence checks at the critical point, time-convergence checks, documented fits with error bars, and language that separates 'consistent with' from 'extracted.' If released with those changes, the critical claim could be taken seriously. As written, it is a promising but unproven application of an otherwise sound method.\n\nI would send this to peer review. The variational framework deserves serious referee time, and the critical-dynamics section needs exactly the scrutiny a good referee would give it. My recommendation: engage with it, but ask for the convergence and uncertainty analysis before accepting the dynamical exponent claim.","headline":"A genuinely useful variational method for open bosonic systems whose critical-dynamics section overclaims: the 'extracted' dynamical exponents are imported from the literature and the gap extraction lacks documented convergence and fitting details.","tokens_in":28685,"tokens_out":1987,"would_cite":true,"duration_ms":25643,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper claims that a variational ansatz of complex-centered Gaussians for the Wigner function, evolved by the Dirac-Frenkel principle and evaluated by automatic differentiation, matches exact diagonalization and scales to…","keywords":["open quantum systems","Wigner function","variational method","Dirac-Frenkel principle","automatic differentiation","Bose-Hubbard model","critical slowing down","Liouvillian spectral gap"],"falsifier":"Compute the Liouvillian gap for a small 2D lattice, such as 4x4, both with the VMG method and with exact diagonalization in a truncated Fock basis; if the VMG gap does not converge to the exact value as the number of Gaussians increases, the finite-time relaxation estimate is not trustworthy. Additionally, refit the 12x12 gap data with twice as many Gaussian components; if the extracted $z$ shifts outside the reported uncertainty, the critical-exponent claim is not converged.","tokens_in":27683,"feed_emoji":"⚛️","tokens_out":4948,"duration_ms":57223,"temperature":0.7,"pith_summary":"This paper introduces a variational method for open quantum bosonic systems that approximates the Wigner function as a sum of complex-centered Gaussian components. It derives equations of motion from the Dirac-Frenkel principle and evaluates all phase-space integrals analytically via automatic differentiation, so accuracy improves systematically as the number of Gaussians is increased. The central claim is that this method matches exact diagonalization on a single-mode benchmark, including negative Wigner fringes, and scales to a 12x12 driven-dissipative Bose-Hubbard lattice. From the relaxation dynamics the authors extract the Liouvillian spectral gap and report a finite-size collapse with the critical exponents of the 2D quantum Ising universality class, implying critical slowing down in the thermodynamic limit. If correct, this provides a scalable route to strongly correlated open bosonic dynamics beyond existing tensor-network and phase-space methods.","feed_headline":"Variational Wigner method tracks 144-mode quantum critical decay","feed_subtitle":"Multi-Gaussian ansatz plus automatic differentiation matches exact diagonalization and finds 2D Ising critical exponents.","key_machinery":"The central object is the Variational Multi-Gaussian ansatz $W_\\theta(\\xi)=\\sum_i \\mathrm{Re}[G(\\xi;\\theta_i)]$ for the Wigner function, where each $G$ is a normalized Gaussian with complex center $\\mu = \\alpha + i\\beta$. Complex centers produce oscillatory modulations of the Gaussian envelope and capture negative interference fringes; the even-parity ansatz adds mirrored Gaussians. The dynamics is governed by the Dirac-Frenkel equations $T\\,d\\theta/dt = V$, where $T$ is the Wigner-space quantum geometric tensor (overlap of parameter derivatives) and $V$ the Liouvillian gradient; both reduce to generalized Gaussian moments. The key technical move is to evaluate these moments as derivatives of a closed-form generating function $Z[J,\\tilde J]$, using Taylor-mode automatic differentiation to handle derivatives up to order six, which keeps the method linear in the number of modes and Gaussians.","core_discovery":"On the paper's own terms, the discovery is that a Variational Multi-Gaussian ansatz for the Wigner function, evolved with the Dirac-Frenkel principle and evaluated with Taylor-mode automatic differentiation, gives controlled and scalable open quantum bosonic dynamics. The ansatz represents negative Wigner regions through complex-centered Gaussians, and the paper shows exponential reduction of observable error with the number of Gaussian components in the single-mode Kerr parametric oscillator. Applied to a two-dimensional Bose-Hubbard lattice with two-boson driving and losses, the method reproduces steady-state parity from corner-space renormalization and extends to lattice sizes up to 12x12, i.e. 144 modes. The paper's strongest claim is the finite-size scaling collapse of the Liouvillian gap using exponents $\\beta = 0.32641871$, $\\nu = 0.62997097$, and $z = 2.0235$, identifying the transition as 2D quantum Ising and demonstrating critical slowing down.","pith_inferences":["Beyond the paper's benchmarks, the same machinery should extend to spin or fermionic phase-space representations, since the Liouvillian differential structure is generic; the authors themselves point toward such extensions.","The reported $z = 2.0235$ is taken from a five-loop epsilon expansion of the 2D Ising model, not derived by the variational method; a direct extraction of $z$ by fitting the VMG gap data without fixing the Ising value would be a stronger independent test.","One could test the method on quench dynamics that create cat states or on regimes with multiple steady states; the paper does not explore these, and they would probe the ansatz's expressivity beyond the studied parameter window."],"forward_implications":["The method provides a systematically convergent variational route for open bosonic systems deep in the quantum regime, including transient Wigner negativities.","For the studied driven-dissipative Bose-Hubbard model, the Liouvillian gap vanishes in the thermodynamic limit with 2D quantum Ising exponents, so the phase transition shows critical slowing down.","The approach extends steady-state benchmarks to lattices of 144 modes, beyond previously accessible corner-space or tensor-network limits.","In the single-mode benchmark, increasing the number of Gaussians by a factor of four reduces the relative observable error by roughly four orders of magnitude, indicating controlled convergence with $N_G$.","Because the Liouvillian acts as a polynomial differential operator in phase space, the same automatic-differentiation machinery can be applied to other analytical ansatze and phase-space representations."],"supporting_citations":[{"why":"Supplies the Dirac-Frenkel variational principle used to derive the equations of motion $T d\\theta/dt = V$.","marker":"[54]"},{"why":"Introduces the corner-space renormalization method used as a benchmark for steady-state observables.","marker":"[23]"},{"why":"Provides the steady-state finite-size scaling of the quadratically driven photonic lattice that the paper reproduces and extends.","marker":"[24]"},{"why":"Establishes that complex-centered Gaussian functions reproduce the negative interference fringes of Schrodinger-cat states, justifying the ansatz.","marker":"[45]"},{"why":"Supplies Taylor-mode automatic differentiation used to compute the generalized Gaussian moments.","marker":"[53]"},{"why":"Provides the exact-diagonalization numerical library used for the single-mode benchmark.","marker":"[59]"},{"why":"Justifies identifying the asymptotic decay rate with the Liouvillian spectral gap and its closing at a dissipative phase transition.","marker":"[68]"},{"why":"Source of the 2D quantum Ising critical exponents $\\beta$ and $\\nu$ used in the finite-size collapse.","marker":"[61–63]"},{"why":"Source of the dynamical critical exponent $z = 2.0235$ used in the rescaling of the Liouvillian gap.","marker":"[71]"}],"fun_headline_variants":["Variational multi-Gaussian method hits 144-mode quantum criticality","Auto-diff accelerates variational Wigner dynamics for open Bose-Hubbard","Multi-Gaussian ansatz captures 2D Ising critical exponents in quantum lattice","Variational Wigner method scales to 144 modes with auto-diff","Critical decay in 144-mode open quantum lattice via variational multi-Gaussian"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The method's central assumption is that the variational equations of motion stay accurate over very long times on the large lattice, so the decay rate read off from finite trajectories equals the true asymptotic Liouvillian gap; convergence is verified explicitly only on the single-mode test case, not on the 12x12 lattice used for the critical exponents.","fun_headline_variants_meta":{"raw":{"variants":["Variational multi-Gaussian method hits 144-mode quantum criticality","Auto-diff accelerates variational Wigner dynamics for open Bose-Hubbard","Multi-Gaussian ansatz captures 2D Ising critical exponents in quantum lattice","Variational Wigner method scales to 144 modes with auto-diff","Critical decay in 144-mode open quantum lattice via variational multi-Gaussian"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000969,"raw_usage":{"total_tokens":4104,"prompt_tokens":908,"completion_tokens":3196,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":524,"completion_tokens_details":{"reasoning_tokens":3099}},"tokens_in":524,"tokens_out":3196,"duration_ms":25768,"temperature":1.0,"reasoning_tokens":3099,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T16:01:20.613391+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the Liouvillian gap for a small 2D lattice, such as 4x4, both with the VMG method and with exact diagonalization in a truncated Fock basis; if the VMG gap does not converge to the exact value as the number of Gaussians increases, the finite-time relaxation estimate is not trustworthy. Additionally, refit the 12x12 gap data with twice as many Gaussian components; if the extracted $z$ shifts outside the reported uncertainty, the critical-exponent claim is not converged.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the Dirac-Frenkel variational principle used to derive the equations of motion $T d\\theta/dt = V$."},{"cited_title":"Kilda, A","cited_arxiv_id":null,"evidence_quote":"Provides the steady-state finite-size scaling of the quadratically driven photonic lattice that the paper reproduces and extends."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Establishes that complex-centered Gaussian functions reproduce the negative interference fringes of Schrodinger-cat states, justifying the ansatz."},{"cited_title":"Vogel and H","cited_arxiv_id":null,"evidence_quote":"Supplies Taylor-mode automatic differentiation used to compute the generalized Gaussian moments."},{"cited_title":"Griewank, J","cited_arxiv_id":null,"evidence_quote":"Provides the exact-diagonalization numerical library used for the single-mode benchmark."},{"cited_title":"Tosca, M","cited_arxiv_id":null,"evidence_quote":"Justifies identifying the asymptotic decay rate with the Liouvillian spectral gap and its closing at a dissipative phase transition."}],"review_version":1}