Separation of scales: dynamical approximations for composite quantum systems

We consider composite quantum-dynamical systems that can be partitioned into weakly interacting subsystems, similar to system–bath type situations. Using a factorized wave function ansatz, we mathematically characterize dynamical scale separation. Specifically, we investigate a coupling régime that is partially flat, i.e. slowly varying with respect to one set of variables, for example, those of the bath. Further, we study the situation where one of the sets of variables is semiclassically scaled and derive a quantum–classical formulation. In both situations, we propose two schemes of dimension reduction: one based on Taylor expansion (collocation) and the other one based on partial averaging (mean-field). We analyze the error for the wave function and for the action of observables, obtaining comparable estimates for both approaches. The present study is the first step towards a general analysis of scale separation in the context of tensorized wavefunction representations.


Introduction
We consider composite quantum-dynamical systems that can be partitioned into weakly interacting subsystems, with the goal of developing effective dynamical descriptions that simplify the original, fully quantum-mechanical formulation. Typical examples are small reactive molecular fragments embedded in a large molecular bath, namely, a protein, or a solvent, all being governed dynamically by quite distinct energy and time scales. To this end, various régimes of intersystem couplings are considered, and a quantum-classical approximation is explored. A key aspect is dimension reduction at the wave function level, without referring to the conventional "reduced dynamics" approaches that are employed in system-bath theories.
1.1. The mathematical setting. The quantum system is described by a timedependent Schrödinger equation (1.1) i∂ t ψ = Hψ ; ψ |t=0 = ψ 0 , governed by a Hamiltonian of the form where the coupling potential W (x, y) is a smooth function, that satisfies growth estimates guaranteeing existence and uniqueness of the solution to the Schrödinger equation (1.1) for a rather general set of initial data, as we shall see later in Section 2. The overall set of space variables is denoted as (x, y) ∈ R n × R d such that the total dimension of the configuration space is n + d. The wave function depends on time t 0 and both space variables, that is, ψ = ψ(t, x, y). We suppose that initially scales are separable, that is, we work with initial data of product form (1.2) ψ(0, x, y) = ψ 0 (x, y) = ϕ x 0 (x)ϕ y 0 (y). In the simple case without coupling, that is, W ≡ 0, the solution stays separated, ψ(t, x, y) = ϕ x (t, x)ϕ y (t, y) for all time, where i∂ t ϕ x = H x ϕ x ; ϕ x |t=0 = ϕ x 0 , i∂ t ϕ y = H y ϕ y ; ϕ y |t=0 = ϕ y 0 , and this is an exact formula. Here, we aim at investigating the case of an actual coupling with ∂ x ∂ y W ≡ 0 and look for approximate solutions of the form ψ app (t, x, y) = ϕ x (t, x)ϕ y (t, y), where the individual components satisfy evolution equations that account for the coupling between the variables. The main motivation for such approximations is dimension reduction, since ϕ x (t, x) and ϕ y (t, y) depend on variables of lower dimension than the initial (x, y). Of crucial importance is the choice of the approximate Hamiltonian H app = H x + H y + W app (t, x, y), that governs the approximate dynamics. We consider two different approximate coupling potentials: One is the time-dependent Hartree mean-field potential, the other one, computationally less demanding, is based on a brute force single-point collocation. Time-dependent Hartree methods have been known for a long time and have earned the reputation of oversimplifying the dynamics of real molecular systems. We emphasize, that our present study does not aim at rejuvenating but at deriving rigorous mathematical error estimates, which seem to be missing in the literature. Surprisingly, our error analysis provides similar estimates for both methods, the collocation and the mean-field approach. We investigate the size of the difference between the true and the approximate solution in the L 2 -norm, ψ(t) − ψ app (t) L 2 = R n+d |ψ(t, x, y) − ψ app (t, x, y)| 2 dxdy and in Sobolev norms. We present error estimates that explicitly depend on derivatives of the coupling potential W (x, y) and on moments of the approximate solution.
As an additional error measure we also consider the deviation of true and approximate expectation values ψ(t), Aψ(t) − ψ app (t), Aψ app (t) , for self-adjoint linear operators A. Roughly speaking, the estimates we obtain for observables depend on one more derivative of the coupling potential W (x, y) than the norm estimates. This means that in many situations expectation values are more accurately described than the wave function itself. Even though rigorous error estimates that quantify the decoupling of quantum subsystems in terms of flatness properties of the coupling potential W (x, y) are naturally important, our results here seem to be the first ones of their kind.
1.2. Relation with previous work. Interacting quantum systems have traditionally been formulated from the point of view of reduced dynamics theories, based on quantum master equations in a Markovian or non-Markovian setting [3]. More recently, also tensorized representations of the full quantum system have been considered, as for example by matrix product states (MPS) [33,34] or within a multiconfiguration time-dependent Hartree (MCTDH) approach [5,15,29,39]. Both wavefunctions (pure states) and density operators (mixed states) can be described in this framework, and wavefunction-based computations can be used to obtain density matrices [36]. In the chemical physics literature, dimension reduction for quantum systems has been proposed in the context of mean-field methods [14,16], and the quantum-classical mean-field Ehrenfest approach [13,2]. Also, quantumclassical formulations have been derived in Wigner phase space [27,21] and in a quantum hydrodynamic setting [17,6,35]. Our present mathematical formulation circumvents formal difficulties of these approaches [10,37,32], by preserving a quantum wavefunction description for the entire system. Previous mathematical work we are aware of is concerned with rather specific coupling models, as for example the coupling of Hartree-Fock and classical equations in [7], or the time-dependent self-consistent field equations in [20], or with adiabatic approximations which rely on eigenfunctions for one part of the system, see for example [38,28]. To the best of our knowledge, the rather general mathematical analysis of scale separation in quantum systems we are developing here is new.
1.3. Partially flat coupling. For a first approximate potential, we consider a brute-force approach, where we collocate partially at a single point, for definiteness we choose the origin, and set In comparison, following the more conventional time-dependent Hartree approach, we set where we perform partial and full averages of the coupling potential, For both approximations, the brute-force and the mean-field approximation, we derive various types of estimates for the error in L 2 -norm. Our key finding is that both methods come with error bounds that are qualitatively the same, since they draw from either evaluations or averages of the function Depending on whether one chooses to control the auxiliary function δW in terms of ∇ x W , ∇ y W or ∇ x ∇ y W , the estimate requires a balancing with corresponding moments of the approximate solution. For example, Proposition 3.3 provides L 2norm estimates of the form x yϕ y (s) L 2 y ds, where const ∈ {1, 2, 4}, depending on whether ψ app (t) results from the brute-force or the mean-field approximation. Example 3.5 discusses important variants of this estimate using different ways of quantifying the flatness of the coupling potential. Proposition 3.9 gives analogous estimates in Sobolev norms. In addition, we analyze the deviation of the true and the approximate expectation values in a similar vein. For the expectation values, we again obtain qualitatively similar error estimates for W bf and W mf . The upper bounds differ from the norm bounds in so far as they involve one more derivative of the coupling potential W and low order Sobolev norms of the approximate solution, see Proposition 3.11. Hence, from the perspective of approximation accuracy, the brute force and the mean-field approach differ only slightly. Therefore, other assessment criteria are needed for explaining the prevalence of the Hartree method in many applications, as we will discuss in Section 3.5.
1.4. Dimension reduction via semiclassical analysis. In the second part of the paper we turn to a specific case of the previous general class of coupled Hamiltonians H ε = H x + H ε y + W (x, y) and consider for one part of the system a semiclassically scaled Schrödinger operator We will discuss in Section 5.1 system-bath Hamiltonians that can be recast in this semiclassical format. The initial data are still a product of the form 1.2, but the y-factor is chosen as that is, ϕ y 0 is a semiclassical wave packet with a smooth and rapidly decaying amplitude function a ∈ S(R d ), and an arbitrary phase space centre (q 0 , p 0 ) ∈ R 2d . We will choose a semiclassical wave packet approximation for ϕ y (t, y) exploring two different choices for the centre (q(t), p(t)). As a first option we consider the classical trajectoryq = p ,ṗ = −∇V 2 (q), and as a second option the corresponding trajectory resulting from the averaged gradient of the potential V 2 , Correspondingly, the approximative factor ϕ x (t, x) is evolved by the partial Hamil- We obtain error estimates in L 2 -norm, see Proposition 5.5 and expectation values, see Proposition 5.8. These estimates are given in terms of the semiclassical parameter ε and derivatives of the coupling potential. Again, both choices for the effective potential differ only slightly in approximation accuracy. Measuring the coupling strength in terms of η = ∇ y W L ∞ , we obtain two-parameter estimates of order √ ε + η/ √ ε in norm, while the ones for the expectation values are of order ε + η. Hence, again the accuracy of quadratic observables is higher than the one for wavefunctions.

Assumptions
We describe here the mathematical setting that will be ours, and discuss it in the context of system-bath Hamiltonians [41,3]. Our Hamiltonian is of the form , where the potentials V 1 (x), V 2 (y) and the coupling potential W (x, y) are all smooth functions, that satisfy growth conditions as given in Assumption 2.1. We will denote V (x, y) = V 1 (x) + V 2 (y) + W (x, y) and abbreviate the Lebesgue spaces for the different variables x, y, and (x, y) by . The initial data ψ 0 (x, y) in (1.2) are products of functions ϕ x 0 ∈ L 2 x and ϕ y 0 ∈ L 2 y , that are square-integrable and typically, Schwartz class, see below.
2.1. Assumptions on regularity and growth of the potentials. We choose a very classical set of assumptions on the regularity and the growth of the potential, since our focus is more on finding appropriate ways to approximate the solution in a standard framework than on treating specific situations. Assumption 2.1. All the potentials that we consider are smooth, real-valued, and at most quadratic in their variables: and, for α ∈ N n , β ∈ N d , We also assume that ∇ y W ∈ L ∞ , but note that this condition can easily be relaxed, see Example 3.5. All the initial date we consider are smooth and rapidly decaying, that is, Schwartz class functions: ϕ x 0 ∈ S(R n ; C), ϕ y 0 ∈ S(R d ; C) (hence ψ 0 ∈ S(R n+d ; C)). Under the above assumption, it is well-known that H x , H y and H are essentially self-adjoint on L 2 (R N ), with N = n, d and n + d, respectively (as a consequence of Faris-Lavine Theorem, see e.g. [ with ω j 0 (possibly anisotropic harmonic potential), β j ∈ R, A j ∈ R n×n positive definite symmetric matrices (Gaussian potential), and V per a smooth potential, periodic along some lattice in R n .
The assumptions on the growth and smoothness of the potentials and the regularity of the initial data call for comments. Remark 2.3. Concerning the growth of V 1 , V 2 and W , the assumption that they are at most quadratic concerns the behavior at infinity and could be relaxed, up to suitable sign assumptions. Local behavior is rather free, for example a local double well is allowed, as long as it is not too confining at infinity. We choose to stick to the at most quadratic case, since bounded second order derivatives simplify the presentation.
Remark 2.4. Concerning the smoothness of our potentials, most of our results still hold assuming only smoothness of W , as long as the operators H x and H y are essentially self-adjoint on an adequate domain included in L 2 (R N ), with N = n, d. For example, V 1 and V 2 could both present Coulomb singularities, and the results of Proposition 3.3 would still hold. In the semiclassical régime, we can also allow a Coulomb singularity for V 1 and prove Proposition 5.5 and Proposition 5.8.
Remark 2.5. Concerning the smoothness and the decay of the initial data, most of our results still hold, if the initial data are contained in one of the spaces Σ k (R N ) containing functions f whose norms are bounded. Note that S(R N ) = ∩ k∈N Σ k . For example, Proposition 3.3 still holds for initial data in Σ 1 , while Proposition 5.5 requires initial data in a semiclassically scaled Σ 3 space.

System-bath
Hamiltonians. An important class of coupled quantum systems are described by system-bath Hamiltonians [41,3].
These are naturally given in the format required by (6.1). In the present discussion, we specify that the bath is described by a harmonic oscillator, V b (y) = 1 2 k 0 2 |y| 2 (or a set of harmonic oscillators in more than one dimension) and the system-bath coupling V sb (x, y) = W (x, y) is of cubic form, such that we obtain in the notation of (6.1), where k 0 2 > 0 and η ∈ R d . The cubic, anharmonic coupling W (x, y) is a non-trivial case, which is employed, e.g., in the description of vibrational dephasing [24,18] and Fermi resonances [4]. It is natural to assume smoothness and subquadratic growth for V s (x). However, the coupling potential W (x, y) clearly fails to satisfy the growth condition of Assumption 2.1. Moreover, it is not clear that in such a framework the total Hamiltonian H is essentially self-adjoint. On the other hand, adding a quartic confining potential, guarantees that H is essentially self-adjoint. Indeed, Young's inequality for products yields for some constants C 1 , C 2 0. Hence, the Faris-Lavine Theorem implies that H x , H y and H are essentially self-adjoint. In the following, we will therefore also provide slight extensions of our estimates to accommodate this specific, but interesting type of coupling (see Remark 3.8).

Partially flat coupling
In this section, we present error estimates that reflect partial flatness properties of the coupling potential W (x, y) in the sense, that quantities like ∇ y W L ∞ or ∇ x ∇ y W L ∞ are small. Depending on the regularity of the initial data, the smallness of these norms could also be relaxed to the smallness of x −σx y −σy ∇ y W L ∞ for some σ x , σ y 0, see Example 3.5. We investigate two approximation strategies, one that is based on brute-force collocation, the other one on spatial averaging. In each case, we prove that the coupling in (x, y) is negligible at leading order with respect to ∇ y W . Throughout this section, ψ = ψ(t, x, y) denotes the solution to the initial value problem (1.1)-(1.2).

It satisfies the equation
with approximative potential W app (x, y) = W (x, 0) + W (0, y) − W (0, 0). The last term Σ ψ controls the error ψ − ψ app , as we will see more precisely below. Saying that the coupling potential W is flat in y means that ∇ y W is small, and we write . This suggests that partial flatness of W implies smallness of the approximation error.
The Taylor expansion which corresponds to the standard normal mode expansion. Hence, choosing (x 0 , y 0 ) such that the maximal singular value of M (x 0 , y 0 ) is minimized, we minimize the error of the brute-force approach.

3.2.
Mean-field approach. Instead of pointwise evaluations of the coupling potential, we might also use partial averages for an approximation. We consider where we have denoted where we have used the fact that the L 2 -norms of φ x and φ y are independent of time, since W is real-valued. Note that (3.3) is the nonlinear system of equations of the time-dependent Hartree approximation. Contrary to the brute-force approach, L 2 regularity does not suffice to define partial averages in general. In view of Assumption 2.1, a fixed point argument (very similar to the proof of e.g. [8, Lemma 13.10]) shows that this system has a unique solution with the phase given by the full average It solves the equation Remark 3.2. The correcting phase e i t 0 W ds seems to be crucial if we want to compute the wave function. On the other hand, since it is a purely time dependent phase factor, it does not affect the usual quadratic observables. The same applies for the phase e itW (0,0) of the brute-force approximation.
3.3. Error estimates for wave functions. We begin with an approximation result at the level of L 2 -norms only. For its proof, see Section 4. Proposition 3.3. Under Assumption 2.1, we have the following error estimates: x yφ y (s) L 2 y ds. We see that the smallness of ∇ y W L ∞ controls the error between the exact and the approximate solution in both approaches.
In this special product case, the crucial source term takes the form and Proposition 3.3 can be augmented by the gradient-free estimate The L ∞ -norms, that provide the upper bounds in Proposition 3.3, separate as y , and it is ∇ y W 2 L ∞ that controls the estimates. Suppose we have W 2 (y) = η|y| 2 with η small: ∇W 2 is not bounded, but we can adapt the proof of Proposition 3.3 to get that is, the extra power of y is transferred to the φ y term.
Remark 3.6. In the spirit of the last observations of Example 3.5, in terms of η := x −σx y −σy ∇ y W L ∞ , for some σ x , σ y 0, we get in Proposition 3.3: See B for details of the argument.
c > 0 and k a positive integer, a typical case where V 1 may be super-quadratic while H x and H remain self-adjoint), then we can estimate xφ x L 2 x uniformly in time.
x t , and the order of magnitude in t is sharp, corresponding to a dispersive phenomenon.
Remark 3.8. The framework of a cubic system-bath coupling W (x, y) = 1 2 η · x|y| 2 as described in Section 2.2 is recovered by taking σ x = σ y = 1 in Example 3.5. In addition, in the presence of a quartic confinement with k 0 4 > 0, in view of Remark 3.7, we also know that y |y|φ y (t) L 2 y is bounded uniformly in t. Adding control on the gradients of the functions ϕ x (t), ϕ y (t) respectively φ x (t), φ y (t), allows also error estimates at the level of the kinetic energy. For a proof, see B.2. Brute-force approach: for ψ app (t, x, y) = ϕ x (t, x)ϕ y (t, y)e itW (0,0) defined by (3.1), Remark 3.10. The strategy used to prove Proposition 3.9 can be iterated to infer error estimates in Sobolev spaces of higher order, H k (R n+d ) for k 2, provided that we consider momenta of the same order k, which explains the interest in the functional spaces Σ k . Error estimates in such spaces can also be obtained by first proving that ψ and the approximate solution(s) remain in Σ k , and then interpolating with the L 2 error estimate from Proposition 3.3.

3.4.
Error estimates for quadratic observables. For obtaining quadratic estimates, we consider observables such as the energy or the momenta, that is, operators that are differential operators of order at most 2 with bounded smooth coefficients. These differential operators have their domain in H 2 (R n+d ), as the operator H. More generally, we could consider pseudo-differential operators B = op(b) associated with a smooth real-valued function b = b(Z) with Z = (z, ζ) ∈ R 2(n+d) , whose action on functions f ∈ S(R n+d ) is given by We assume that b satisfies the Hörmander condition We shall also consider observables that depend only on the variable x or the variable y. The following estimates are proven in B.3.
Proposition 3.11. Under Assumption 2.1, for b ∈ C ∞ (R n+d ) satisfying 3.5 and B = op(b), there exists a constant C b > 0 such that we have the following error estimates: Remark 3.12. The averaging process involved in the action of an observable on a wave function allows to prove estimates like the one in Proposition 3.11, that are more precise than the standard ones stemming from norm estimates, Remark 3.13. We point out that the error is governed by derivatives of second order in W , involving a derivative in the y variable that is supposed to be small. Besides, note that the direct use of an estimate on the wave function itself would have involved H 2 norms of ψ app (s), while this estimate only requires H 1 norms. This first improvement is due to the averaging process present in Egorov Theorem.
3.5. Energy conservation. The error estimates of Proposition 3.3, Proposition 3.9, and Proposition 3.11, do not allow to distinguish between the brute-force single point collocation and the mean-field Hartree approach. However, in computational practice most of the employed methods are of mean-field type. Why? Our previous analysis, that specifically addresses the coupling of quantum systems, does not allow for an answer, and we resort to a more general point of view. Both approximations, the brute-force and the mean-field one, are norm-conserving. However, the mean-field approach is energy-conserving with the same energy as 1.
where the mean-field Hamiltonian is given by Below in B.5 we give an elementary ad-hoc proof of Lemma 3.14. However this conserved value does not correspond to the exact energy of (1.1), but only to an approximation of it.

An exemplary proof
Here we discuss our basic proof strategy and apply it for the norm estimate of Proposition 3.3. The norm estimates of Example 3.5, Proposition 3.9 and Proposition 5.5 and the observable estimates of Proposition 3.11 and Proposition 5.8 are carried out in B and C. where ψ 0 ∈ L 2 (R N ) and Σ ∈ L 1 loc (R + ; L 2 (R N )). Then for all t 0, This standard lemma is our main tool for proving norm estimates. It will be applied with either h = 1 or h = ε as parameter. Its proof is given in A. Now we present the proof of Proposition 3.3.
We note that both approximations and their components are norm-conserving for all times t 0, that is, • In the case of the brute-force approach, we consider the Taylor expansions (3.2) and derive the estimates • In the mean-field approach, we note that for (t, x, y) ∈ R × R n+d , Like we did in the brute-force approach, we may use either of the estimates In the first case, we come up with where we have used Cauchy-Schwarz inequality for the last term. We infer where we have used Young inequality (α + β) 2 and finally x . In the case of the second type approximation for δW , we similarly find x yφ y (t) L 2 y . Proposition 3.3 then follows from Lemma 4.1 with h = 1.

Dimension reduction via semiclassical analysis
In this section, we consider coupled systems, where one part is governed by a semiclassically scaled Hamiltonian, that is, First we motivate such a partial semiclassical scaling in the context of system-bath Hamiltonians and introduce wave packets as natural initial data for the semiclassical part of the system. We explore partial semiclassical wave packet dynamics guided by classical trajectories and by trajectories with averaged potentials. Thus, the partially highly-oscillatory evolution of a PDE in dimension n + d is reduced to a less-oscillatory PDEs in dimensions n, and ODEs in dimension d. The corresponding error estimates in Section 5.5 compare the true and the approximate product solution in norm and with respect to expectation values.
5.1. Semiclassical scaling. We reconsider the system-bath Hamiltonian with cubic coupling of Section 2.2, now formulated in physical coordinates (X, Y ), that is, where the coordinates X and Y of the system and the bath part are prescaled, resulting in the single mass parameters µ 1 , µ 2 for each subsystem and one single harmonic frequency ω 2 for the bath (noting that, alternatively, several harmonic bath frequencies ω 2,j could be introduced, without modifying the conclusions detailed below). The corresponding time-dependent Schrödinger equation reads We perform a local harmonic expansion of the potential V s (X) around the origin X = 0 and assume that it is possible to determine a dominant frequency ω 1 . We then define the natural length scale of the system as Rescaling coordinates as (x, y) = 1 a (X, Y ), we obtain where we have introduced the dimensionless parameters η. The rescaling of the system potential V s and the coupling vector η do not alter their role in the Hamiltonian, whereas the two dimensionless parameters ε and deserve further attention. We now consider the régime where both the mass ratio ε between system and bath and the frequency ratio between bath and system are small, that is, where the system is viewed as "light" and "fast" when compared to the "heavy" and "slow" bath.  [19,40]). The resulting dimensionless mass ratios are given as 12 = µ 1 /µ 2 = 0.109 and 13 = µ 1 /µ 3 = 0.51, and the corresponding frequency ratios are 12 = ω 2 /ω 1 = 0.005 and 13 = ω 3 /ω 1 = 0.006. In the case of the H 2 -Kr relative motion, note that the frequency ratio 13 is indeed small whereas the mass ratio is 13 ∼ 0.5; this shows that the quantum-classical boundary is less clearly defined than in the first example of coupled electronic-nuclear motions. In such cases, different choices can be made in defining the quantum-classical partitioning.
In an idealized setting, where ε is considered as a small positive parameter whose size can be arbitrarily small, we would say that = O(ε) as ε → 0, and view the system-bath Hamiltonian H sb as an instance of a partially semiclassical operator whose potentials V 1 (x) and V 2 (y) are independent of the semiclassical parameter ε and satisfy the growth conditions of Assumption 2.1. As emphasized in Section 2.2, the cubic coupling potential W (x, y) does not satisfy the subquadratic estimate, but can be controlled by additional moments of the approximate solution. A corresponding rescaling of time, t = εω 1 τ , translates the time-dependent Schrödinger equation to its semiclassical counterpart where the physical and the rescaled wave functions are related via ψ ε (t, x, y) = a (n+d)/2 Ψ(τ /(εω 1 ), aX, aY ).
Remark 5.3. Criteria for justifying a semiclassical description are somewhat versatile in the literature. Our scaling analysis shows, that for system-bath Hamiltonians with cubic coupling the obviously small parameter ε, that describes a ratio of reduced masses, has to be complemented by an equally small ratio of frequencies .
Otherwise, the standard form of an ε-scaled Hamiltonian, as it is typically assumed in the mathematical literature, does not seem appropriate.

5.2.
Semiclassical initial data and ansatz. As before, the initial data separate scales, where we now assume that g ε is a semiclassically scaled wave packet, with (q 0 , p 0 ) ∈ R 2d , a smooth and rapidly decreasing, i.e. a ∈ S(R d ; C) = ∩ k∈N Σ k .
In the typical case, where the bath is almost structureless (say, near harmonic), the amplitude a could be chosen as a complex Gaussian, but not necessarily. We now seek an approximate solution of the form ψ ε app (t, x, y) = ψ ε 1 (t, x)ψ ε 2 (t, y), where ψ ε 2 is a semiclassically scaled wave packet for all time, Here, (q(t), p(t)) ∈ R 2d , the phase S(t) ∈ R, and the amplitude u 2 (t) ∈ S(R d , C) must be determined.
Remark 5.4. We note that our approximation ansatz differs from the adiabatic one, that would write the full Hamiltonian as is an operator, that parametrically depends on the "slow" variable y and acts on the "fast" degrees of freedom x. From the adiabatic point of view, one would then construct an approximate solution as ψ ε bo (t, x, y) = Φ(x, y)ψ ε 2 (t, y), where Φ(x, y) is an eigenfunction of the operator H f (y); here, the subscript "bo" stands for Born-Oppenheimer. The result of Corollary C.3 emphasizes the difference between these two points of view.
We denote by the part of the approximate solution that just contains the amplitude. With this notation, The analysis developed in the next two sections allows to derive two different approximations, based on ordinary differential equations governing the semiclassical wave packet part, which are justified in Section 5.5 (see Proposition 5.5).

Approximation by partial Taylor expansion.
Plugging the expression of ψ ε app (t, x, y) into (5.1) and writing y = q(t) + z √ ε in combination with the Taylor expansions we find: where the argument of u ε app and its derivatives are taken in z = y−q(t) √ ε . To cancel the first four terms in the √ ε line, it is natural to require Now cancelling the first four terms in the first line of the right hand side yields In other words, (q(t), p(t)) is the classical trajectory in y, and S(t) is the associated classical action. At this stage, we note that the term z · ∇ y W (x, q)u ε app is not compatible with decoupling the variables x and z (or equivalently, x and y). Using that ∇ y W L ∞ is assumed to be small, the above computation becomes In view of (5.5), we set Equation (5.9) is a Schrödinger equation with a time-dependent harmonic potential: it has a unique solution in L 2 as soon as a ∈ L 2 (R d ). In addition, since a ∈ Σ k for all k ∈ N, u 2 ∈ C(R; Σ k z ) for all k ∈ N. The validity of this approximation is stated in Proposition 5.5 below. Note that if a is a Gaussian state, then u 2 too and its (time-dependent) parameters -width matrix and centre point -can be computed by solving ODEs (see e.g. [26,8,12,23] and references therein).

5.4.
Approximation by partial averaging. Following e.g. [11,31] or [23, Section 2], we write where the averages are with respect to |ψ ε 2 (t, y)| 2 dy. For example, where we anticipate the fact that the L 2 y -norm of ψ ε 2 (t) is independent of time. We almost literally repeat the previous argument and find that , and z is taken as z = (y − q(t))/ √ ε. The parameters satisfy the equations of motioṅ We see that we can now define the approximate solution by: Since the matrix ∇ 2 V 2 y (t) is real-valued, we infer that the L 2 z -norm of u 2 (t) is independent of time, hence ψ 2 (t) L 2 y = u 2 (t) L 2 z = a L 2 . The equation in u 2 is now nonlinear, and can be solved in Σ 1 , since ∇V 2 is at most linear in its argument: u 2 ∈ C(R; Σ 1 z ), and higher Σ k regularity is propagated. Here again, if a is a Gaussian, then so is u 2 and its width and centre can be computed by solving ODEs (see [26,8]). Note also that, differently from the previous setting, u 2 is now ε-dependent via the quantity ∇ 2 V 2 y (t) (see (5.10)). However, this dependence is very weak since a Taylor expansion in (5.10) shows that u 2 is close in any Σ k norm from the solution of the equation (5.9). For this reason, we do not keep memory of this ε-dependence and write u 2 . By contrast, the ε-dependence of ψ ε 1 is strong since it involves oscillatory features in time.

5.5.
The approximation results. The main outcome of the approximations can be stated as follows, and is proved in C.1: Proposition 5.5. Let ψ ε be the solution to (5.1)-(5.2), with g ε given by (5.3). Then with ψ ε app given either like in Section 5.3 or like in Section 5.4, there exist constants K 0 , K 1 independent of ε such that for all t 0, Corollary 5.6. Assume η := ∇ y W L ∞ √ ε, then for all T > 0, Remark 5.7. Using the same techniques as in B.2, one can prove estimates on higher regularity norms, using ε-derivatives in y and standard ones in x. For example, if a ∈ Σ 4 , then there exists K 0 , K 1 independent of ε such that We refer to [9] (see also [8,Chapter 12]) for more detailed computations.
Note that, in both approximations, the evolution of u 2 corresponds to the standard quadratic approximation. In particular, if a is Gaussian, then u 2 is Gaussian at all time, and solving the equation in u 2 amounts to solving ordinary differential equations. However, the equation (5.8) solved by ψ 1 (t) is still quantum, such that a reduction of the total space dimension of the quantum system has been made from n + d to n.
Let us now discuss the approximation of observables that we choose as acting only in the variable y. Due to the presence of the small parameter ε, we choose semiclassical observables and associate with b ∈ C ∞ c (R 2d ) (b smooth and compactly supported) the operator op ε (b) whose action on functions f ∈ S(R d ) is given by y + y 2 , ξ e iξ·(y−y )/ε f (y ) dξdy .
As before, the error estimate is better for quadratic observables than for the wave functions. More specifically, the following result, that is proved in C.2, improves the error estimate from Proposition 5.5 by a factor √ ε.
Then with b ∈ C ∞ c (R 2d ) and ψ ε app given either like in Section 5.3 or in Section 5.4, there exist a constant K independent of ε such that for all t 0, Remark 5.9. Of course, we could have considered a mixed setting consisting of pseudodifferential operators as in Section 3.4 in the variable x, and semiclassical as above in the variable y. One would then obtain estimates mixing those of Proposition 3.11 and Proposition 5.8.

A numerical example
For an illustrative numerical application, we consider a system-bath type Hamiltonian with cubic coupling W (x, y) = 1 2 ηxy 2 , as developed in Section 2.2, i∂ t ψ(t, x, y) = H sb ψ(t, x, y) ; ψ(0, x, y) = ϕ x 0 (x)ϕ y 0 (y), in dimension d = n = 1,  Table 1. Parameters defining the four numerically simulated variations of the system-bath Hamiltonian (6.1). The blue model uses the coupling parameter η ref = −k 0 2 /(3a 1 ) and the frequency ratio ref = 1/100. The red model varies the frequency ratio, the grey and yellow models the coupling strength.
The mass ratio between the system and the bath is moderately small, µ 1 /µ 2 = 0.25. The system potential is a quartic double well, while the bath potential is harmonic, The length scale = 4a 1 of the double well is a multiple of the system's natural harmonic unit a 1 . The initial data ψ 0 (x, y) = ϕ x 0 (x)ϕ y 0 (y) are the ground state of the bivariate harmonic oscillator, that results from the harmonic approximation of V s (x) + V b (y) around the left well (x, y) = (0, 0). The coupling constant η < 0 is negative to ensure that the total Hamiltonian's ground state is localised in the right well (x, y) = (2 , 0), providing a setup with pronounced non-equilibrium dynamics. For such a system-bath model, the gradient-free error estimate (3.4) of Example 3.5 is given by x (s) y 4 y (s) − y 2 2 y (s) ds. Figure 1 presents the results from the following numerical experiment: We identify a frequency ratio ref = 1/100 between bath and system and a coupling parameter η ref = −k 0 2 /(3a 1 ) as generating "reference" system-bath dynamics. For this parameter choice, the mean-field approximation, when compared with a numerically converged MCTDH approximation, results in roughly a 0.1% error after 10 units of the natural harmonic time scale t 1 = 1/ω 1 of the system (blue curve), see Figure 1a. Hence, the Hartree approximation is excellent on the time scale under study. Decreasing the frequency ratio by a factor of four, roughly halves the error (red curve). And, as expected, an increase in the coupling strength also increases the error (grey curve), while decreasing the coupling also decreases the error (yellow curve). The corresponding plot of Figure 1b illustrates that the upper bound of the theoretical error estimate correctly captures the initial slope of all four error curves, while slightly over-estimating the actual error as time evolves. A more detailed assessment of the error estimate (6.2), in particular of its long-time behaviour (up to 150 ps), when the mean-field approximation goes up to errors of the order of 1%, and a more complete screening of physically relevant parameter régimes are work in progress for a numerical companion paper to the present theoretical study.   Table 1. The right plot shows the corresponding upper bound of the error estimate (6.2).

Conclusion and outlook
We have presented quantitative error bounds for the approximation of quantum dynamical wave functions in product form. For both considered approaches, a brute-force single point collocation and the conventional mean-field Hartree approximation, we have obtained similar error estimates in L 2 -norm (Proposition 3.3, Example 3.5), in H 1 -norm (Proposition 3.9), and for quadratic observables (Proposition 3.11). To our knowledge, such general estimates, that quantify decoupling in terms of flatness properties of the coupling potential, are new. The corresponding analysis for semiclassical subsystems (Proposition 5.5, Proposition 5.8) confirms the more general finding, that error estimates for quadratic observables provide smaller bounds than related norm estimates. The single product analysis, as presented here, provides a well-posed starting point for the investigation of more elaborate approximation methods. If the initial data satisfy then we may invoke the linearity of (1.1) to write ψ(t, x, y) = J j=1 ψ j (t, x, y), where each ψ j solves i∂ t ψ j = Hψ j , with ψ j|t=0 = ϕ x 0j ⊗ ϕ y 0j . We approximate each ψ j ≈ ψ j,app individually in one of the ways discussed in the present paper and use the triangle inequality for where the norm can be an L 2 or an energy norm, for instance. However, working on each ψ j instead of ψ directly, seems to prevent control of the limit J → ∞.
Multi-configuration methods therefore use ansatz functions of the form where the families (ϕ (x) j (t)) j 1 and (ϕ (y) (t)) 1 satisfy orthonormality or rank conditions, while gauge constraints lift redundancies for the coefficients a k (t) ∈ C. We view our contributions here as an important first step for a systematic assessment of such multi-configuration approximations in the context of coupled quantum systems. A numerical companion paper, that explores the dynamics of system-bath Hamiltonians with cubic coupling on multiple time-scales with respect to various parameter régimes, is currently in preparation.
and therefore, by the Cauchy-Schwarz inequality, d dt ψ(t) 1 h Σ(t) . Integrating in time, we obtain In the context of observables, refined error estimates will follow from the application of the following lemma.
, and ψ (1) , ψ (2) , solutions to the homogeneous Cauchy problems where ψ 0 ∈ L 2 (R N ). Then, for all t 0 Proof. We denote the unitary evolution operators by U j (t) = e −iAj t/h and calculate Appendix B. Proof of error estimates: partially flat coupling In this section, we prove error estimates in L 2 -norm for general potentials W (Remark 3.6), in H 1 -norm (Proposition 3.9), and for quadratic observables (Proposition 3.11).
B.1. Proof of Remark 3.6. To prove the estimates of Example 3.5, recall that we have denoted The Fundamental Theorem of Calculus yields and we replace the pointwise estimate of δW with |δW (x, x , y, y )| |y − y | x σx + x σx max ( y , y ) σy η.
The estimate on Σ φ becomes , and we conclude by resuming the same estimates as above: and, in view of the inequality |y−y | max y σy , y σy 2 max y σy |y|, y σy |y | 2 y σy |y| + y σy |y | , we find B.2. Proof of Proposition 3.9. To prove error estimates in H 1 (R n+d ), we differentiate (4.1) in space, and two aspects must be considered: (i) In our framework, the operator ∇ x,y does not commute with H. (ii) We must estimate ∇ x,y Σ ψ and ∇ x,y Σ φ . Indeed, we compute In the typical case where V 1 is harmonic, ∇ x V 1 is linear in x, and so xr ψ appears as a source term. Note that in the general setting of Assumption 2.1, |∇ x V 1 (x)| x .
Remark B.1. If ∇ x V 1 and ∇ x W are bounded, then Lemma 4.1 yields The term r ψ (s) L 2 is estimated in Proposition 3.3, and ∇ x Σ ψ (s) L 2 is estimated below.
Multiplying (4.1) by x, we find similarly Energy estimates provided by Lemma 4.1 applied to the equation for ∇ x r ψ and xr ψ then yield a closed system of estimates: where we have used the estimate |∇ x V 1 + ∇ x W | C(1 + |x|), and the uncertainty principle (uncertainty in x, Cauchy-Schwarz in y), The Gronwall Lemma then yields for some C > 0. We compute The first term in controlled by |y| ∇ x ∇ y W L ∞ |ψ app |. The second term is controlled like in Section 3.3, by replacing ψ app with ∇ x ψ app . We can of course resume the same approach when considering ∇ y r ψ , and the analogue of the above first term is now controlled by |x| ∇ x ∇ y W L ∞ |ψ app |.
By Egorov Theorem, see [42,Theorem 11.1], the operator B(σ) is also a pseudodifferential operator, that is, B(σ) = op(b(σ)) for some function b(σ) that satisfies the growth condition 3.5. We have with the notations of B.1. Then, by the direct estimate of Lemma B.2 below, where C b(σ) > 0 depends on derivative bounds for the function b(σ) and We therefore obtain Using the rectangular n × d matrix M (x, y) introduced in Remark 3.1, the gradient of δW (x, y) can be written as We estimate the Sobolev norm by so that integration in time provides where the constant C b = max σ∈[0,t] C b(σ) depends on derivatives of b. In the meanfield case, the approximate Hamiltonian is time-dependent, The difference of the Hamiltonians is also a function, which is now time-dependent,  1 and b = b(z, ζ) be a smooth function on R 2N satisfying the Hörmander growth condition 3.5. Let δW be a smooth function on R N with bounded derivatives. Then, there exist constants C b > 0 and M N > 0 such that Proof. We explicitly write the commutator as We Taylor expand the function δW (z) around the point z , so that Corresponding to the above decomposition, we write [op(b), δW ]ψ(z) = f 1 (z)+f 2 (z) and estimate the two summands separately. We observe that (z − z )e iζ·(z−z ) = −i∇ ζ e iζ·(z−z ) and perform an integration by parts to obtain where the constant C b > 0 depends on derivative bounds of the function b. For the remainder term of the above Taylor approximation we write and obtain that where C b > 0 depends on derivative bounds of b, and M N > 0 depends on the dimension N .
B.5. Proof of energy conservation. Here we provide an elementary ad-hoc proof for energy conservation of the time-dependent Hartree approximation, Lemma 3.14.
Proof. A first observation is that where we have used the self-adjointness of H mf (t) and norm-conservation in the multiplicative components. Secondly, since the approximate energy coincides with the actual energy, and we obtain the aimed for result.
Appendix C. Proof of error estimates in the semiclassical régime Here we present the proofs of the semiclassical estimates given in Proposition 5.5 and Proposition 5.8. C.1. Error estimates for the wave function. In this section, we prove Proposition 5.5, and comment on the constants K 0 , K 1 , which may be analyzed more explicitly in some cases.
Proof. Section 5.3 defines an approximate solution of the from C.1.2. Approximation by partial averaging.
Proof. The semiclassical approximation obtained by partial averaging reads: To estimate the size ofr 1 andr 2 introduced in Section 5.4, we might argue again via Taylor expansion. Indeed, Hence, we have for all averages f = V 2 , ∇V 2 , ∇ 2 V 2 , W (x, ·) that f y (t) = f (q(t)) + O( √ ε), where the error constant depends on moments of |u 2 | 2 . In particular, if u 2 (0) is Gaussian, the odd moments of |u 2 (t, z)| 2 vanish, and the above estimate improves to O(ε). Hence, the L 2 -norm ofr 1 is O( √ ε) close to the L 2 -norm of r 1 , and the L 2 -norm ofr 2 is O(η √ ε), η = ∇ y W L ∞ , close to the L 2 -norm of r 2 (with each time an extra √ ε gain in the above mentioned Gaussian case). In particular, the order of magnitude for the difference between exact and approximate solution is the same as in the previous subsection, only multiplicative constants are affected. We emphasize that the constants C 0 and C 1 from the previous subsection are in general delicate to assess. On the other hand, in specific cases (typically when u 2 is Gaussian and ∇ 2 V 2 is known), they can be computed rather explicitly. C.2. Error estimates for quadratic observables. The proof of Proposition 5.8 is discussed in the next two sections. C.2.1. Approximation by Taylor expansion.
We also assume that some eigenvalue Λ j (t) is separated from the remainder of the spectrum for all t ∈ R and that the initial datum ϕ x 0 is in the eigenspace of h(0) for the eigenvalue Λ j (0): Then adiabatic theory as developed by Kato [22] states that ψ 1 (t) stays in the eigenspace of Λ j (t) on finite time, up to a phase.
Proposition C.2 (Kato [22]). Assume we have (C.1) and that Λ j (0) is a simple eigenvalue of h(0) such that there exists δ 0 > 0 for which Denote by Φ x j (t) a family of normalized eigenvectors of h(t) such that Then, for all T > 0, there exists a constant C T > 0 such that In contrast to the Born-Oppenheimer point of view recalled in Remark 5.4, we obtain the following time-adiabatic extension for our wave-packet approximation: Corollary C.3. In the setting of Proposition 5.5 and Proposition C.2, we obtain the following approximate solution D.1. TDH and MCTDH calculations. For each set of parameters we compared mean-field (TDH) calculations (one single-particle function per mode) to converged Multiconfiguration Time-Dependent Hartree (MCTDH) calculations (three singleparticle functions per mode, i.e., nine configurations). The latter were compared to numerically "exact" standard propagations on a two-dimensional grid, showing virtually no difference. In all cases, we chose for the primitive basis set a direct product of Gauss-Hermite functions (128 × 1024) associated to the local harmonic approximation of the potential energy centered at (x, y) = (0, 0). We used default integrator schemes for both TDH and MCTDH: norm and energy preservation were dutifully checked. The initial datum of each situation was chosen as a quasi-coherent state corresponding to the local harmonic ground state (centered at (x, y) = (0, 0), with no extra momentum, and with natural widths 1/ √ 2 in natural units).