Numerical Bogoliubov approximation of bosonic many-body quantum dynamics

1. Introduction.

• The problem is approximating many-boson Schrödinger dynamics, whose exact solution is generally unavailable; the setting assumes bounded, even interactions in the mean-field regime.The presentation sets the external trapping potential to zero.Source: §1 · pp. 1–2

1.1. Existing numerical schemes.

• Direct MCTDH has cost O(r^N), while second-quantized MCTDHB still scales polynomially with N, motivating methods suited to large particle numbers.The cited MCTDHB works are.Source: §1.1 · p. 2

1.2. Elements of Bogoliubov theory and approximation ansatz.

• The method uses a Bogoliubov-inspired ansatz of the condensate plus a few excitation functions, a time-dependent strategy the authors say prior numerical work had not used.The paper describes its approximation as first order around the Hartree product and notes it is not intended to match other methods at small or moderate N.Source: §1.2 · p. 3

1.3. Overview & main ideas.

• The paper organizes the task around evolving the ansatz continuously and updating it efficiently over discrete time steps.Its parameterization replaces the full N-particle wave function with scalar coefficients and the two-particle excitation function χ2.Source: §1.3 · p. 4, §1.3.1 · p. 4

1.3.1. How to evolve 𝑢in continuous time? Our starting point is the approximation ansatz (11), viz.

• The Dirac–Frenkel variational principle derives parameter equations while the condensate follows the Hartree equation; the resulting parameter dynamics vary slowly, uniformly in N.The coefficient evolution is generated by a Hermitian matrix with slowly varying components uniformly bounded in N.Source: §1.3.1 · pp. 4–5
• The discretization aims to keep computational effort uniformly bounded in N while addressing large-variable brakets, many interaction operators, permutation sums, kinetic steps, and the costly two-particle function.The overview identifies these as the main computational hurdles and pairs symmetry-based reductions with operator splitting and low-rank approximation.Source: §1.3.1 · p. 5

1.4. Outline.

• The paper proceeds from notation and motion equations to complexity reduction, then numerical methods and experiments on conservation and Bogoliubov predictions at finite N.Source: §1.4 · p. 5

1.5. Notation.

• The notation section establishes conventions for time derivatives, inner products, marginals, projections, and symmetric function spaces used in the derivations.Source: §1.5 · p. 6

2.1. Derivation.

• The equations of motion for the ansatz coefficients and excitation functions are derived using the Dirac–Frenkel time-dependent variational principle, without imposing tensorization.The ansatz and condensate function together define an approximate solution to the Schrödinger equation.Source: §2.1 · pp. 7–8, §2 · p. 7

2.2. Tensorization and variational principle, and a variational crime.

• Under tensorization constraints, the Dirac–Frenkel principle requires inverting a nontrivial operator whose inverse may be arbitrarily ill-conditioned or nonexistent.The authors therefore replace higher excitation functions in the equation for χ₂ with their tensorized forms, yielding differential equations for the coefficients and χ₂ and algebraic equations for higher excitations.Source: §2.2 · p. 11

2.3. Behaviour with respect to 𝑁.

• For Hartree-driven φ, excitation derivatives are uniformly bounded in N, while the coefficient vector c follows a smooth unitary evolution.The coefficient evolution is governed by a Hermitian matrix whose entries are slowly varying in time and uniformly bounded with respect to N.Source: §2.3 · pp. 11, 16
• Permutation symmetry reduces the effective kinetic and potential operators to m + 2 and m(m − 1)/2 + 2m + 3 operators, respectively, uniformly in N.The unreduced operators involve N one-body terms and N(N − 1)/2 two-body terms.Source: §2.3 · pp. 12–13
• The matrix ˇH has diagonal and off-diagonal elements uniformly bounded in N and is smooth in time because it depends on smooth φ and excitation functions χ₂ⱼ.Source: §2.3 · p. 16
• Choosing μφ(t) as in (7) eliminates the highly oscillatory phase proportional to N that another choice of μ would produce in the coefficient vector c.Source: §2.3 · p. 17
• Theorems 2 and 3 show that the equations of motion (23) have right-hand sides bounded as N increases, and the next section addresses computing them with complexity uniformly bounded in N.Source: §2.3 · p. 17

3. Computational complexity of the terms in the equations of motion

• The cost of evaluating the equations-of-motion terms is uniformly bounded with respect to the number N of particles.Source: §3 · p. 17

3.1. Reduction of terms by a pigeonhole argument.

• A braket vanishes when 2p > 2q + 2 or 2q > 2p + m + 2, though these conditions are sufficient rather than necessary.Source: §3.1 · p. 17

3.2. Reducing the number of operators.

• The kinetic terms reduce to 2 one-body operators for the coefficient equation and 3 for the excitation equation; the potential terms reduce to 3 and 4 two-body operators, respectively.Source: §3.2 · p. 18

3.3. Reducing the number of permutations and variables in the braket.

• Symmetry and orthogonality reduce permutation counts to a number independent of N and limit braket computation to at most 2p + 2 variables instead of N − m.Source: §3.3 · pp. 19–21

4.1. Time and space discretizations.

• The discretization samples a computational box on a uniform trapezoidal-rule grid with M points per dimension, applies periodic boundaries, and advances the solution at times t_k = kΔt.The box is chosen large enough to contain the region where the solution is not negligibly small.Source: §4.1 · p. 21

4.2. Numerical solution with 𝑛> 1.

• For n > 1, higher excitation functions are held fixed while 𝜒2 is integrated, then reconstructed by tensorization at each time-step end.This avoids the difficult or potentially ill-conditioned operator inversion while preserving the tensorized ansatz after every step.Source: §4.2 · pp. 21–22

4.3. Operator splitting.

• Strang splitting separates free kinetic evolution from potential evolution in three substeps, yielding a time-reversible, second-order method with local error O(Δt^3).The free Schrödinger substeps are integrated exactly in Fourier space, avoiding the stepsize restriction associated with the unbounded Laplacian.Source: §4.3 · p. 22

4.4. Low-rank approximation of 𝜒2𝑗.

• The second-order rank-adaptive BUG integrator evolves the low-rank two-particle excitation while remaining robust to arbitrarily small singular values and adapting rank to at most 4r per step.It updates the basis from the range of an evolved factor, then advances the coefficient matrix in the new basis; a modification preserves complex symmetry.Source: §4.4 · pp. 23–24

4.4.1. Using low-rank approximations to further reduce computational complexity.

• Low-rank factors and orthonormality let the method skip bracket evaluations whose paired marginals contribute inner products equal to one, further reducing computational cost.The factorization expresses the discretized symmetric 𝜒2 as UBU⊤, and the symmetry and low-rank structure are preserved during time integration.Source: §4.4 · p. 23, §4.4.1 · pp. 24–25

4.5. Dealing with the potential using convolutions.

• Low-rank excitation functions reduce potential calculations to one-particle functions, with FFT-based convolutions computed in O(d M^d log M) and the BUG right-hand side used only in matrix products.The method does not treat that right-hand side as a dense matrix after space discretization.Source: §4.5 · p. 26
• The treatment using convolutions avoids the challenging projection of the marginals of F when computing χ2's equation of motion.Source: §4.5 · p. 27

4.6. Time integration of the equations of motion.

• The potential step uses Strang splitting, integrates only χ2, and reconstructs higher excitations from χ2 and the angles.The χ2 equation is integrated with a complex-symmetric BUG variant; splitting the differential system introduces error dependent on the scheme and time step Δt.Source: §4.6 · pp. 27–28
• The split advances the coefficients, angles, and φ, restores χ2's orthogonality to the current φ, and then integrates χ2.With other quantities held fixed, the coefficient evolution is unitary, and φ is integrated exactly.Source: §4.6 · p. 28

4.6.1. Preserving the orthogonality.

• After φ evolves, χ2 is replaced by a nearby function satisfying the orthogonality constraint, obtained by projection and normalization.The low-rank expression is made unitary using an Autonne–Takagi decomposition.Source: §4.6.1 · pp. 28–29

4.7. Algorithm.

• Algorithm 1 summarizes the numerical method.Source: §4.7 · p. 29

5. Numerical experiments

• Both contact-interaction and soft-Coulomb experiments use n=6 and a rank-adaptive BUG(2r) integrator with r_max=15 to investigate the numerical method itself, not its relevance for precise physics applications.Source: §5 · pp. 29, 31

5.1. Double well with contact interaction.

• Energy conservation reflects competition between finite-N and timestep errors: at N=10^2 the N-error dominates for Δt∈{10^-2,10^-3}, whereas at large N the timestep error dominates and conservation is second-order in time.At N=10^6, decreasing Δt improves energy conservation but not up to second order. With Δt=10^-2, results for N=10^6 and N=10^10 are quantitatively very similar, unlike those for N=10^2.Source: §5.1 · pp. 31–32
• For the double-well contact interaction, the smallest eigenvalue of B is relatively robust to N and slightly smaller for smaller Δt.The authors relate the N-robustness to the interaction size being parameterized by λ, independent of N.Source: §5.1 · p. 32
• The initial ordering |c_2| > |c_4| > ··· > |c_2n| is approximately preserved; increasing |c_2|^2, …, |c_2n|^2 with a sharp decrease in |c_0|^2 shows the Hartree product is a very poor approximation.The values of |c_2j|^2 have little dependence on N.Source: §5.1 · p. 32

5.2. Harmonic potential with Coulomb interaction. √︃

• For the soft-Coulomb experiment, energy conservation again reflects competing N-error and Δt-error, while the Hartree product remains a reasonable approximation and the smallest eigenvalue of B remains robust to N.The excitation-coefficient ordering is approximately preserved, with faster decay across j than in the contact-interaction experiment.Source: §5.2 · p. 32

6. Conclusion

• The ansatz represents the wave function using the condensate function and a few excitation functions, so the full condensate need not be explicitly formed and computational effort is uniformly bounded in N.The condensate evolves by the Hartree equation; coefficients and excitation functions are the quantities approximated numerically.Source: §6 · p. 32
• The method splits the Hamiltonian and equations of motion, derives uniformly N-smooth equations through the Dirac–Frenkel principle, and evolves rank-r excitation approximations with a dynamically updated BUG integrator.The Laplacian part is solved exactly, while the potential part is approximated numerically.Source: §6 · p. 32
• The resulting algorithm has computational cost independent of particle number, and numerical experiments show second-order accuracy and satisfying conservation properties.Figure 1 describes the double-well contact-interaction experiment with r_final=r_max=15.Source: §6 · p. 32, Figure 1 · p. 33

AI disclosure

• The authors used Vibe for proofreading and copyediting support.The disclosure names Mistral and the GLM model by Z.ai.Source: AI disclosure · p. 34
• Vibe is identified as the Mistral assistant using the GLM model by Z.ai.Source: AI disclosure · p. 34
• The authors state that no AI-based tool was used to implement the algorithm or numerical experiments, apart possibly from external libraries they could not check.Source: AI disclosure · p. 34

Appendix A. Technical lemmas

• The appendix collects technical lemmas used throughout the work.Source: Appendix A · p. 34
• The authors explicitly disclaim novelty for the collected results.Source: Appendix A · p. 34
• The authors report that they did not find references for these results.Source: Appendix A · p. 34

Appendix B. Components of the tridiagonal Hermitian matrix H

• Bogoliubov orthogonality makes the kinetic matrix T diagonal, while H is split into kinetic and potential matrices T and V.The component calculations are attributed to Lemmas 2, 3, and 10.Source: Appendix B · pp. 36–37

Appendix C. BUG integrator for complex symmetric matrices

• At the end of a step, the integrator can truncate the basis by removing the smallest singular values of S under a chosen smallness condition.The preceding implementation discussion focuses on differential equations for complex symmetric matrices.Source: Appendix C · p. 37, §C.1 · p. 38

C.1. Classical BUG integrator.

• The classical BUG integrator evolves a low-rank matrix using unitary bases U and V and an invertible core matrix S.Its update proceeds through K-, L-, and S-steps, followed by orthonormalization of the updated bases.Source: §C.1 · p. 38

C.2. Preserving symmetry constraints during the BUG integration.

• The integrator must retain a complex-symmetric factorization at every substep because the differential-equation right-hand side accepts only symmetric arguments.These restrictions follow from deriving the equations under the assumption that χ2 is symmetric.Source: §C.2 · p. 38
• Setting V=U makes the L-substep equivalent to the K-substep, so only K- and S-substeps are needed.This uses the required symmetric form U S U⊤ throughout the integration.Source: §C.2 · p. 38
• An Autonne–Takagi representation of χ2 lets the K-step evaluate the right-hand side on symmetric arguments while updating the low-rank factors.The method rewrites U S U⊤ as Ũ Σ Ũ⊤ and uses that form for the BUG update.Source: §C.2 · pp. 38–39

C.3. Enforcing a diagonal core matrix.

• Diagonalizing the core at explicit Runge–Kutta substeps makes right-hand-side evaluations efficient when the core is diagonal.The updated basis remains unitary after diagonalizing the intermediate core.Source: §C.3 · p. 39

C.4. Autonne–Takagi factorization.

• Every complex symmetric matrix A can be factored as W Σ W⊤, where W is unitary and Σ is diagonal.The section presents this as the symmetric equivalent of the SVD.Source: §C.4 · p. 39
• The diagonal entries of Σ are the singular values of A.This identifies the diagonal factor in the Autonne–Takagi decomposition.Source: §C.4 · p. 39
• The section points to [15, Corollary 4.4.4] for a proof and [16, Section 4] for an algorithm based on the usual SVD.These references support the factorization result and its computation, respectively.Source: §C.4 · p. 39

D.1. Computation of 𝜒2𝑘.

• Low-rank structure reduces storage for χ2k from M^2k values to a core tensor of size r^2k, though that core may still be too large.The passage motivates computing values without storing the full core tensor.Source: §D.1 · p. 40
• With a diagonal core matrix, only index permutations pairing equal values contribute, requiring every index value to occur an even number of times.The contributing permutations pair indices within subgroups of equal values.Source: §D.1 · p. 40
• The counting groups contributing permutations by equal-index subgroups and accounts for rearrangements of the paired terms.The derivation counts permutations within each subgroup before adjusting for permutations that reorder product terms.Source: §D.1 · p. 40
• When the multi-index is sorted, each right-hand-side value can be computed in O(k) time in one iteration and evaluated on the fly.The norm of the associated basis tensor can be computed once and stored whenever b is updated.Source: §D.1 · p. 40

• Core tensors b are obtained through tensorization, and a summand is nonzero only when both factors in its product are nonzero.Source: §D.2 · p. 41
• A core-tensor entry can be nonzero only when every index in i′ occurs an even number of times; G_l counts occurrences of index l.Source: §D.2 · p. 41
• When the two sets differ, P_C,iQ_C is zero; when they match, the sum in (60) can be restricted to multi-indices i′.Source: §D.2 · p. 41
• Replacing the iteration over P in (60) with a sum over R_s in [1:r]^s substantially improves computation for large |P| or large r.Source: §D.2 · p. 41