Series solution of the time-dependent Schr\"{o}dinger-Newton equations in the presence of dark energy via the Adomian Decomposition Method

The Schr\"{o}dinger-Newton model is a nonlinear system obtained by coupling the linear Schr\"{o}dinger equation of canonical quantum mechanics with the Poisson equation of Newtonian mechanics. In this paper we investigate the effects of dark energy on the time-dependent Schr\"{o}dinger-Newton equations by including a new source term with energy density $\rho_{\Lambda} = \Lambda c^2/(8\pi G)$, where $\Lambda$ is the cosmological constant, in addition to the particle-mass source term $\rho_m = m|\psi|^2$. The resulting Schr\"{o}dinger-Newton-$\Lambda$ (S-N-$\Lambda$) system cannot be solved exactly, in closed form, and one must resort to either numerical or semianalytical (i.e., series) solution methods. We apply the Adomian Decomposition Method, a very powerful method for solving a large class of nonlinear ordinary and partial differential equations, to obtain accurate series solutions of the S-N-$\Lambda$ system, for the first time. The dark energy dominated regime is also investigated in detail. We then compare our results to existing numerical solutions and analytical estimates, and show that they are consistent with previous findings. Finally, we outline the advantages of using the Adomian Decomposition Method, which allows accurate solutions of the S-N-$\Lambda$ system to be obtained quickly, even with minimal computational resources.

Since the first attempts to create a quantum theory of the gravitational field in the 1930's [1,2], the search for a complete theory of quantum gravity has been one of the major fields of research in theoretical physics. In the pioneering work [1], Bronstein investigated the quantum mechanical measurement of the Γ 0 01 component of the Christoffel symbols. This led to a fundamental limit on the temporal uncertainty intrinsic to any quantum measurement, ∆t ≥ /c 2 Gρ 2 V 1/3 , where V and ρ are the volume and the density of a self-gravitating massive body, respectively. The time uncertainty can then be related to the spatial uncertainty via ∆x ≤ c∆t. By introducing the standard mass-density-volume relation, M = ρV , one also obtains the mass-time-density uncertainty relation, M ≥ /c 2 Gρ (∆t) 3 . This was one of the first generalised uncertainty relations (GURs) and, since then, many others have been proposed in the literature. (See [3,4] for reviews.) Of these, the most widely studied are the generalised uncertainty principle (GUP) and the extended uncertainty principle (EUP). The former aims to incorporate the effects of canonical gravitational attraction between quantum mechanical particles [5,6], while the latter accounts for the repulsive effects of dark energy in the form of a cosmological constant, Λ [7][8][9]. The extended generalised uncertainty principle (EGUP) accounts for both [10][11][12].
Over the years, many theoretical models and diverse approaches to the problem of quantum gravity have been developed. These include postulating the existence of the graviton, the hypothetical spin-2 boson that mediates quantum gravitational interactions, string theory, loop quantum gravity, and noncommutative geometry, to mention just a few of the directions investigated. (For detailed presentations of the different approaches, see [13][14][15][16]. For recent reviews of the present status of quantum gravity research, see [17][18][19][20].) However, recently, important progress in experimental techniques has enabled researchers to cool, control, and measure physical systems in the weak gravity and quantum mechanical regimes, with far greater accuracy than ever before. For the first time, it may be possible to directly observe quantum gravity effects at scales accessible in terrestrial and near-Earth-orbit laboratories, in near-future experiments [21,22].
Nonetheless, the conceptual and technical challenges to the construction of complete theory are manifold. From a quantum field theory perspective, even the Newtonian theory of gravity is problematic. Comparing Newton's law of gravity, Φ(r) = −Gm/r, with Coulomb's law of electrostatics, V (r) = k e q/r, it follows that the gravitational constant has mass dimension −2. Hence, the theory of Newtonian quantum gravity is nonrenormalizable. This result follows from the calculation of the gravitongraviton scattering at energy E, in which the divergent series ∼ 1 + GE 2 + GE 2 2 + ... appears [23]. The nonrenormalizability of such naive quantum gravity theories therefore suggests that new physics should emerge at the Planck scale, M Pl = c/G ≃ 10 19 m proton . Because of these challenges, much research interest has been devoted to the study of semi-classical models. In these models matter fields are quantised while gravity remains classical, or is perturbatively quantised at nextto-leading order in the expansion of the metric. Such an approach to quantum gravity was proposed in [24], which is based on the decomposition of the quantum metric into a classical and a fluctuating part,ĝ µν = g µν + δĝ µν . Further assuming that δĝ µν = K µν = 0, where K µν is a classical tensor, one arrives at an effective gravitational Lagrangian of the form L = L g (ĝ µν ) + √ −gL m (ĝ µν ) ≃ L g + δLg δg µν δĝ µν + √ −gL m + δ( √ −gLm) δg µν δĝ µν , where κ 2 = 8πG/c 4 . The gravitational field equations obtained from this Lagrangian lead to theories that require geometrymatter coupling at the classical level. The coupling is of the kind that also appears in the f (R, T ) type modified gravity theories [25,26] and the cosmological implications of effective field theories with fluctuating metric components were investigated in [27].
However, because gravity is in many ways different from the other fundamental forces, and due to the intrinsic difficulties in its quantization, some researchers have suggested that the gravitational field may be essentially classical, and that it should not and cannot be quantized [28,29]. But, even gravity is not quantum, ordinary matter is. Hence, in order to describe the gravitational dynamics of quantum fields, one must still combine classical gravity with quantized matter. In this scenario quantized matter is coupled to the classical gravitational field by replacing the classical energy momentum tensor, T µν , with the expectation value of the energy-momentum operator, T µν , in Einstein's field equations. The expectation value is constructed by averaging with respect to an appropriately chosen quantum state, Ψ, yielding the semi-classical field equations [30], Equations (1) can be also obtained from the variational principle δS = δ (S g + S ψ ) = 0 [31], where S g = (1/16πG) R √ −gd 4 x is the standard Hilbert-Einstein action of general relativity, while the quantum part of the action, S ψ , is introduced in the form (2) By varying the quantum action (2) with respect to Ψ we obtain the normalization condition Ψ|Ψ = 1 and the following Schrödinger equation for Ψ, in addition to the semi-classical Einstein equations. Note that the Bianchi identities still impose the conservation of the effective energy-momentum tensor, ∇ µ Ψ|T µν |Ψ = 0. The difficulty of building a successful quantum theory of general relativity, as well as the intrinsic problems of treating quantum field theories in curved spacetimes, have also led to the hypothesis that a satisfactory description of quantum gravity could be achieved by unifying quantum mechanics with Newtonian gravity [32]. This corresponds to the weak field limit of Eqs. (1), which reduce to the semi-classical Poisson equation Equation (4) is the basis of the Schrödinger-Newton approach [21] and, in this model, the equation of motion of a self-gravitating massive particle can be formulated as where V is the canonical quantum potential and Φ is the gravitational self-interaction potential, obtained by solving (4). For a system of N non-relativisitc free particles V = 0 and the mass-density operator may be written asρ = This is the standard Schrödinger-Newton equation, whose static and time-dependent solutions have been intensively investigated in . The average value of the self-interaction potential can be obtained as Φ ≃ 2 G/mR 3 and it turns out that the particle behavior is essentially quantum if the condition m 3 R ≪ 2 G ≃ 10 −47 cm g 3 , where R is average radius of the wave function, is satisfied [32].
Due to its extreme nonlinearity, the time-dependent Schrödinger-Newton system has mostly been investigated numerically. An exception is the variational approach, which was considered in [42], where the system of equations (4)-(5) was investigated in the hydrodynamical representation of quantum mechanics. In this formalism the wave function is represented as ψ ( r, t) = √ ρ e iS and the canonical Schrödinger equation reduces to the equations of classical fluid mechanics, in the presence of the quantum potential. The quantum fluid flows with a velocity u = ∇S and the equations of motion can be obtained from the Lagrangian +ρV. (7) By adopting a spherical Gaussian profile for the density, ρ(r, t) = π −3/2 R −3 (t)e −r 2 /R 2 (t) , one can then obtain the gravitational potential and the Lagrangian of the system reduces to L(R,Ṙ) =Ṙ 2 /2 − 1/2R 2 + C/R, where C is an arbitrary constant. The corresponding equation of motion for R isR = 1/R 3 − C/R 2 . Using this formalism, one can obtain the energy eigenvalues, linear frequencies, and nonlinear late-time behavior of the S-N wave packet [42].
More recently, the Schrödinger-Newton system was generalized by considering the effects of dark energy in the form of a cosmological constant Λ [54]. This is consistent with the standard ΛCDM model of cosmology, in which it is assumed that the late-time acceleration of the Universe is driven by a constant vacuum energy density, ρ Λ = Λc 2 /(8πG) ≃ 10 −30 g cm −3 [55]. (For alternative models of dark energy as modified gravity, see [56] and references therein.) The physically interesting regime in which dark energy dominates both gravitational self-attraction and canonical quantum diffusion was investigated numerically and using analytical estimates. It turns out that this takes place for objects with arbitrary mass that are sufficiently delocalized. An estimate of the minimum delocalization width required, of the order of 67 m, was determined, and this prediction was verified by the numerical results.
However, the exact delocalisation radius required for dark energy domination can be much higher for very massive particles. In general, the wave function of a free particle in the S-N-Λ system was found to split into a core region that collapses due to gravitational self-attraction and an outer region that undergoes accelerated diffusion due to presence of dark energy [54]. While the former behaviour is present in the standard S-N model, the latter is unique to the S-N-Λ system. The order of magnitude of the critical radius separating collapse from expansion was found to match analytical estimates of the classical turnaround radius for a massive compact object in the presence of a cosmological constant [57].
The goal of this paper is to further investigate the mathematical and physical properties of the timedependent S-N-Λ system introduced in [54]. In order to obtain a better understanding of the dynamics of the model, we adopt a semi-analytical approach and construct series solutions using the Adomian Decomposition Method (ADM) [58][59][60][61]. This is a powerful method that can be used to obtain accurate series solutions of a large class of nonlinear differential equations, and systems of equations, with applications in diverse fields of science and engineering [62][63][64][65][66][67][68][69][70][71][72][73][74]. Here, we apply it to the S-N-Λ model for the first time.
An essential advantage of this method is that it can be used to obtain analytical approximations to the full numerical solutions, without any need for perturbation theory, closure approximations, linearization, or discretization methods. For many highly nonlinear models, including the S-N and S-N-Λ systems of equations, the use of these methods leads to complicated and time-consuming numerical computations. On the other hand, to obtain even approximate closed-form analytical solutions of a nonlinear problem requires introducing restrictive and simplifying assumptions. The key advantage of the ADM is that it can be used to find the solution of a given equation or system of equations in the form of a rapidly converging power series. Successive terms in the series are obtained via a recursive relation, with the help of a special class of functions known as Adomian polynomials [58][59][60][61]. In most cases the series converges fast, so that the application of this method saves a lot of computational time.
Although the ADM has been used extensively in many areas of engineering and physics, it has been used very little in the study of gravitation and quantum mechan-ics. (For some applications of the method in these fields, see [64,68,70].) In order to apply the method, we must first reformulate the time-dependent S-N-Λ system as a system of two integral equations. We then obtain the series solutions of the system by expanding the nonlinear terms using the Adomian polynomials [58][59][60][61]. To eliminate the unwanted oscillatory behavior of the solution, we represent the Adomian series in terms of their Padé approximants.
After obtaining the recursive relation for the full S-N-Λ system, we test the efficiency of the ADM for a free Gaussian wave packet, in the limit G → 0, Λ → 0. In this case, the canonical Schrödinger equation can be solved exactly, and we show that the ADM recovers the exact solution in just a few simple steps. Next, series solutions are obtained for both the wave function and the gravitational potential, in the presence of gravitational self-interaction and dark energy. The associated probability density is computed with the help of the Padé approximants and we pay special attention to he dark energy dominated regime.
This paper is organized as follows. In Section II we present the basic structure and mathematical formalism of the Adomian Decomposition Method. The S-N-Λ system is reformulated as a system of integral equations in Section III, and the recurrence relations for the series solution of the system are obtained. The method is tested for the case of the canonical Schrödinger equation describing the free propagation of a Gaussian wave packet, and it is shown that the exact solution of this system can be re-obtained in a few simple steps. We obtain the semianalytical solution of the time-dependent S-N-Λ system, for Gaussian initial conditions, in Section IV. The dark energy dominated regime is also considered in detail, and a numerical analysis of the evolution of the probability density is presented. Our results are compared with previous analytical and numerical studies in Sec. V and our discussion and final conclusions are presented in Section VI.

II. THE ADOMIAN DECOMPOSITION METHOD
Let us consider a partial differential equation written in the general form is the linear remainder operator that may contain partial derivatives with respect to x, N [.] is a nonlinear operator, which we assume is analytic, and g is a non-homogeneous term that is independent of u. Equation (8) must be solved with the initial condition u(x, 0) = f (x). We assume thatL t is invertible, so that we can applyL −1 t to both sides, obtaining The ADM posits the existence of a series solution in which u(x, t) is given by while the nonlinear termN [u (x, t)] is decomposed aŝ where {A n } ∞ n=0 are the Adomian polynomials. These are generated according to the rule Substituting the series expansions (10) and (11) into Hence, we can obtain the following recurrence relation, giving the series solution of Eq. (8) as Therefore, we obtain the approximate solution of Eq. (8) as where For a given nonlinearityN [u], the Adomian polynomials are obtained as and so on. The greater the number of terms, the higher the accuracy of the truncated series solution.

III. THE ADOMIAN DECOMPOSITION METHOD FOR THE TIME-DEPENDENT SCHRÖDINGER-NEWTON-Λ SYSTEM
For a single particle of mass m, the time-dependent S-N-Λ system is given by the following two equations, where the last term in the above Poisson equation has been chosen so that the dark energy density is given by its standard form ρ Λ = Λc 2 /(8πG). For spherically symmetric systems, Eqs. (21) and (22) take the form which must be solved with the boundary conditions ψ (r, 0) = Ψ(r) and Φ (r, t) = φ (r), respectively. Introducing the operatorsL t = ∂/∂t andL rr = ∂ 2 /∂r 2 , Eqs. (23) and (24) can be rewritten aŝ These equations can be solved formally to give By taking into account the fact thatL −1 where φ 1 (t) and φ 2 (t) are arbitrary integration functions. Using the Cauchy formula for repeated integration, Eq. (30) can be rewritten as In order to make the gravitational potential finite in the origin r = 0, we chose φ 2 (t) = 0, giving (32) We now determine the series solution of the reformulated system of equations, (29) and (32), by assuming expansions of the form and In addition, we decompose the terms Φ (r, t) ψ (r, t) and |ψ (ξ, t)| 2 in terms of the Adomian polynomials as Substituting the above decompositions into Eqs. (29) and Hence, we obtain the following recursive series solution for the time-dependent S-N-Λ system, The first three Adomian polynomials in each series are obtained as and respectively. For simplicity, we assume that the initial state of the wave function is a spherically symmetric Gaussian, with initial width σ 0 = 1/ √ 2α. The gravitational potential φ(r), corresponding to the time-independent Gaussian wave packet, then satisfies the Poisson equation whose general solution is where erf(z) = (2/ √ π) z 0 e −t 2 dt is the error function, and c 1 and c 2 are arbitrary integration constants. In order to avoid a singularity at the origin, we take c 1 = 0. The initial distribution of the gravitational potential is then finite everywhere, and satisfies the condition lim r→0 φ(r) = c 2 − Gm α/π.

A. Testing the Adomian Decomposition Method
To test the efficiency of the ADM, we consider the evolution of a Gaussian wave packet in the absence of the gravitational interaction, and dark energy, by setting Φ(r, t) = 0. In this case, the evolution of the wave function is given by the canonical Schrödinger equation and must be solved subject to the initial condition ψ(r, 0) = (α/π) 3/4 e −αr 2 (48). The general solution of Eq. (51), satisfying the required initial condition, is given by In order to simplify the mathematical formalism, we introduce a new set of dimensionless variables (τ, θ), defined as where m p denotes the proton mass, and (55) respectively. Moreover, we rescale the wave function so that The rescaled wave function satisfies the dimensionless Schrödinger equation whose general solution, satisfying the initial condition (48), is given bỹ This may be expanded as a power series, with respect to the dimensionless time τ , as To solve Eq. (57) using the ADM, we apply the oper-atorL −1 τ to both sides, giving We then decomposeψ (θ, τ ) into an infinite sum of components, so thatψ (θ, τ ) = ∞ n=0ψ n (θ, τ ), where the componentsψ n (θ, τ ) will be determined recurrently. By substituting the series expansion into Eq. (60), we obtain (61) Thus, we obtain the recursive relations The first few iterations are given bỹ and so on. Clearly, the Adomian series solutionψ (θ, τ ) ≃ ψ 0 (θ, 0) +ψ 1 (θ, τ ) +ψ 2 (θ, τ ) +ψ 3 (θ, τ ) + . . . exactly reproduces the series expansion of the exact solution (59) and, in the limit of an infinite number of iterations, fully recovers it. Hence, we have shown that the ADM gives the exact series representation of the solution of the spherically symmetric, three-dimensional Schrödinger equation describing the time-evolution of a free Gaussian wave packet. The probability distribution P (θ, τ ) is obtained as Here, we approximate the series representation ofP (θ, τ ) by its Padé approximantP [3/4](θ, τ ), which is given bỹ Generally, for a power series of the form f (z) = ∞ z=0 f k z k the Padé approximant of the order (m, n) in the vicinity of the point z = 0 is the rational function Π m.n ∈ R m,n having the property that it takes the closest values to the given series near z = 0. Here, by R m,n , we have denoted the set of rational functions of the form P/Q, where P and Q are polynomials in z of degree p ≤ m and q ≤ n, respectively [75].
The comparison between the exact probability density of the Gaussian quantum wave packet, and its approximation, given by Eq. (67), is represented in Fig. 1. As one can see from this figure, Eq. (67) gives an excellent description of the time-evolution of the wave packet for −1 ≤ τ ≤ 1, and a good approximation for τ outside this range.

IV. SERIES SOLUTION OF THE TIME-DEPENDENT SCHRÖDINGER-NEWTON-Λ SYSTEM
Using the mathematical formalism developed in the previous Section, we now construct explicit series solutions of the S-N-Λ system. In addition, we perform a detailed numerical study. Our main goal is to highlight the effects of self-gravity and dark energy on the timeevolution of a the quantum wave packet.
In the first order of approximation, we obtain (96) and at the second order of approximation, To third order, the probability density can be approximated as Higher order approximations of the probability density of a Gaussian wave packet, evolving under self-gravity in the presence of dark energy, can also be calculated easily with the aid of computer algebra systems. The gravitational self-potential can be approximated as or, in terms of the Padé approximants of the power series, (100)

B. The dark energy dominated regime
We now consider the limiting case in which the dark energy density dominates the matter density, λ ≫ σ ψ (θ, τ ) . The Poisson equation then takes the simple form and can be immediately integrated to givẽ where we have assumed that the background gravitational potential is independent of time. Hence, for the dark energy dominated phase, the Schrödinger equation takes the form and can be formally solved to givẽ By decomposing the wave function asψ(θ, τ ) = ∞ n=0 ψ n (θ, τ ), we obtain the following recurrence relations for the determination of the componentsψ n (θ, τ ): and the first few approximations of the quantum wave packet in the dark energy dominated regime are obtained asψ The probability density of the Gaussian wave packet can be obtained, with the help of the Padé approximants, to different orders of approximation, as and and so on. The analytical expressions for the probability density can also be obtained easily to any desired order of approximation.

C. Numerical analysis
In the final part of this section, we consider the numerical results obtained from the Adomian series solutions of the S-N-Λ system. Our main goal is to highlight the effects of the self-gravitational potential and the dark energy density on the evolution of the probability density associated with the Gaussian quantum wave packet. In Fig. 2, we present the three-dimensional evolution of the rescaled gravitational potential in the absence of dark energy, i.e., with λ = 0, and with σ = 1. For convenience, we takeã = 1 as its initial value. In this case, the gravitational potential can be approximated bỹ and the full solution satisfies the condition lim τ →∞Φ (θ, τ ) = 0. Mathematically, a singularity develops inΦ(θ, τ ) for values of θ satisfying 4θ(4a + σ) − √ 2πσerf √ 2θ = 0. However, to at least third order in the approximation, this equation does not have any real roots, except at θ = 0.
The variation of the self-gravity potential in the presence of dark energy is represented in Fig. 3, for two different values of λ; λ = 0.20 and λ = 0.35. In the large θ limit its behavior can be approximated as Thus, we see that the presence of a positive cosmological constant does have a significant effect on the distribution of the gravitational potential. To at least the considered order of approximation, the condition lim τ →∞Φ (θ, τ ) = 0 still holds. On the other hand, as expected, lim λ→∞Φ (θ, τ ) = −∞. In both cases, thẽ Φ(θ, τ ) has a sharp maximum at the origin of the coordinate system, θ = 0.
The time variation of the probability density of the Gaussian wave packet is represented, for fixed values of the radial coordinate θ, in Figs. 4. There are two significant effects induced by the presence of the dark energy. As one can see from the left-hand panel, for (relatively) small values of the dimensionless radial coordinate θ, the probability density in the presence of Λ > 0 almost coincides with the function describing the evolution with Λ = 0, for −1 ≤ τ ≤ 1, and has the same maximum value. In the absence of Λ, the probability density tends to zero at a finite value of τ . However, dark energy significantly modifies the tail of the Gaussian distribution, which extends in time and induces much higher values of the probability density, as compared to the Λ = 0 case. From the right-hand panel we see that, for larger values of θ, the dark energy has two different effects on the probability density. The first is a significant increase in the amplitude of the probability density, with the maximum increased by a factor of at least two. This indicates the increased probability of finding the wave packet at larger distances from the center, the effect being a direct consequence of the presence of repulsive dark energy. Secondly, at large distances, the probability density tends to zero. However, the decrease is much slower for Λ > 0, and is directly correlated with the increase of the amplitude of the wave. Another interesting effect is related to the change in the shape of the wave function function, which evolves from a single-peaked into a double-peaked symmetric function.
The three-dimensional evolution of the wave packet in the presence of the self-gravitational field and the dark energy density is depicted in Figs. 5. The same effects, as previously mentioned, are also apparent when considering the three-dimensional evolution of the wave packet.
For large values of τ and θ, P (θ, τ ) → 0, but the dynamics of the transition to the asymptotic limit are strongly influenced by the presence of dark energy, whose effect becomes significant at late times and for large values of the radial coordinate.
The behavior of the probability density in the dark energy dominated regime is presented in Fig. 6, for two distinct physical situations, corresponding to a fixed value of θ (left panel), and to a fixed value of τ (right panel). Even though, in this regime, there is a qualitative similarity with the Λ = 0 case, significant differences also appear. The double-peaked shape of the Gaussian distribution is extended in time for fixed θ, and the shape of the Gaussian tail is strongly modified, indicating an increase in the probability of finding the particle at higher values of τ . Moreover, the maximum value of the probability as a function of θ, at a given time, increases dramatically with increasing λ. However, at least at the considered order of approximation, and for the adopted values of λ in the large θ limit, the probability distribution of the initially Gaussian wave packet still tends to zero. Nonetheless, much larger values of λ, in the range λ ∈ 10 2 , 10 3 , would greatly modify the dynamics of the wave packet at infinity.

V. COMPARISON OF THE ADOMIAN METHOD WITH PREVIOUS ANALYTICAL AND NUMERICAL RESULTS
A particle obeying the S-N-Λ equation of motion experiences three tendencies in its dynamics. Both canoni-cal quantum diffusion and dark energy induced acceleration cause its wave function to spread whereas Newtonian self-gravity, represented by the non-linear term, acts to localize the wave packet. In [54], it was argued that the relative strengths of these three tendencies can be esti- mated, at least approximately, by considering the motion of the peak radial probability density, r p (t). This is the position of the spherical shell at which the radial probability density dP/dr = 4πr 2 |ψ| 2 reaches its maximum, that is, the radius at which the particle is most likely to be found at a given time t. It is determined by solving the equation or, alternatively, which is equivalent to setting d 2 P (r, t)/dr 2 = 0 or d 2P (θ, τ )/dθ 2 = 0, respectively. The contributions to the total acceleration experienced by r p (t) due to canonical quantum diffusion, self-gravity, and dark energy are then estimated as and The subscript SE refers to the canonical Schrödinger equation, SN refers to the standard Schrödinger-Newton contribution, and Λ denotes the additional term induced by the dark energy density. In order to determine the regimes in which the different tendencies dominate the dynamics, we consider equality between the absolute magnitudes of the accelerations (117)-(119) in a pair-wise manner, i.e., or, equivalently, plus giving and where λ C (m) = /(mc) is the reduced Compton wavelength of the particle, r S (m) = 2Gm/c 2 is its Schwarzschild radius, l Pl = G /c 3 ≃ 10 −33 cm is the Planck length, and l dS = 3/Λ ≃ 10 28 cm is de Sitter radius. Note that the latter is comparable to the present day radius of the Universe [76] and that we have neglected numerical factors of order unity in all three equations.
The critical value of r p (t) in Eq. (124) is the classical turn-around radius for a spherical compact object in the Schwarzschild-de Sitter spacetime [57], Mass Behavior Below 1 × 10 −18 g Evolution indistinguishable from that of a free particle in canonical quantum mechanics 2 × 10 −18 g to 3 × 10 −17 g The whole wave packet spreads faster than that of a canonical free particle 4 × 10 −17 g to 5 × 10 −17 g The inner core of the wave function spreads slower than the wave function of the canonical free particle while the outer shells spread faster 6 × 10 −17 g to 1 × 10 −16 g The inner core of the wave function collapses under self-gravity while the outer shells spread faster than in canonical quantum mechanics ∼ 2 × 10 −16 g Chaotic Above 3 × 10 −16 g Stationary   TABLE I. Dynamical evolution of a Gaussian wave packet, with initial width σ0 = 7.5 × 10 2 cm, under the S-N-Λ equation, according to the numerical solution obtained in [54]. The comparison is made to a free particle in canonical quantum mechanics, evolving under the canonical Schrödinger equation.
or, equivalently, where m Pl = c/G ≃ 10 −5 g and m dS = ( /c) Λ/3 ≃ 10 −66 g are the Planck mass and the de Sitter mass, respectively. The approximate value of the peak radial probability is and a more careful estimate, accounting accurately for numerical factors, gives r p ≃ 67 m, as shown in [54]. For a Gaussian distribution, the initial peak radial probability is comparable to the initial width of the wave function, r p (0) ≃ σ 0 , and the two are equivalent up to a multiplicative constant of order one for a large class of physically reasonable wave functions [54]. Therefore, Eq. (129) also gives the order of magnitude value of the minimum initial width required, in order for the acceleration due to dark energy to dominate both canonical quantum diffusion and self-gravitation. This is a very clear and somewhat surprising prediction: in the S-N-Λ system, the spreading of any spherically symmetric wave packet with an initial width σ 0 67 m will be dominated by the accelerated expansion of the Universe, due to dark energy, regardless of its initial mass. For particles with masses m 10 −17 g, Eq. (122) implies that the dark energy term always dominates over canonical quantum diffusion, whenever the initial width of the wave packet exceeds this critical value. For heavier particles, we expect the outer shells of the wave packet to undergo accelerated expansion due to dark energy while the inner core region contracts under self-gravity. By Eq. (124), the critical radius marking the division between collapsing and expanding shells should be of the order of the classical turn-around radius (126).
However, these very strong predictions were derived using rather crude analytical techniques and approximations. It is therefore reasonable to ask: can they be trusted? To answer this question, the numerical solution of the S-N-Λ system was presented in [54], for an initially Gaussian wave packet with a range of initial widths and particle masses. Remarkably, the existence of both a critical mass of order 10 −17 g, and of critical initial width of order σ 0 ≃ 6.7 × 10 2 cm, was verified by the numerical results. A summary of the numerical results obtained in [54], for a particle wave function of initial width σ 0 = 7.5 × 10 2 cm, and particle masses in the range 10 −18 kg ≤ m ≤ 10 −16 kg, is given in Table 1.
We note that the chaotic and stationary regimes obtained for larger values of m are artifacts of the numerics, which were unable to probe masses above ∼ 2 × 10 −16 g due to limited computational resources. The critical radius marking the boundary between the collapsing inner core and the expanding outer shells of the wave packet was also verified to be within one order of magnitude of the classical turn-around radius (126), which isn't bad for such a crude analysis [54].
(131) Fig. 7 shows the evolution of θ p (τ ), obtained from the Adomian series solution, for a Gaussian wave packet with a = 1 (describing the effect of the background gravitational field), and for different values of the dimensionless dark energy parameter, λ. The presence of the gravitational field and of the dark energy significantly modifies the behavior of θ p . Although, in the absence of self-gravity, the peak probability density of the Gaussian wave packet satisfies lim τ →∞ θ (0) p = ∞ in the presence of an extremely high dark energy density, corresponding to very large values of λ, the presence of self-gravitational interaction significantly alters the behaviour of θ p , at least at the first order of approximation, which may now tend to zero for finite values of τ . This represents the regime in which the total collapse of the wave function occurs under the action of self-gravitational attraction, which successfully counteracts both dark energy repulsion and canonical quantum diffusion. In the range −1 ≤ τ ≤ 1, the time evolution of θ (1) p closely follows, on a qualitative level, the dynamics of θ (0) p , even though some quantitative differences do appear.
In Figs. 8, the dimensionless radial probability density dP/dθ = 4πθ 2 |ψ| 2 is plotted for fixed θ, and for various values of τ . This clearly shows the formation of a collapsing inner core and an outer shell undergoing accelerated expansion. The critical value of θ that demarcates between the two regions corresponds, to within an order of magnitude, to the classical turn-around radius of the particle mass, and is therefore consistent with the numerical results summarised in Table 1.
Finally, before concluding this section, we note that, since the Compton wavelength of the proton is of order 10 −15 m, lab-based experiments for which the dark energy dominated regime (128)-(129) is accessible require macromolecules with approximately 108 amu. This is two orders of magnitude below the estimated mass required for tests of the standard Schrödinger-Newton equation using opto-mechanical traps [77], which corresponds to the generic estimate for the onset of the semi-classical gravity regime with Λ = 0 [30]. In other words, in terms of the mass parameter, current experiments are sufficiently precise to allow the effects of Λ on the quantum dynamics of a macromolecule to be observed and measured. The associated length scale is σ 0 ≃ 1 − 10 m, though, unfortunately, the associated time-scales may astronomical [54]. However, for macromolecules with ∼ 10 10 amu, the canonical quantum contribution to the peak acceleration is of the same order as the dark energy contribution for σ 0 ≃ 1 m. This raises the intriguing possibility that dark energy effects may be observable in near-future experiments on local quantum systems, though, to date, the preceding order-of-magnitude estimates seem to have been overlooked in the quantum gravity literature. Crucially, the present work shows that we may go beyond such crude estimates, to obtain detailed analytical predictions of the S-N-Λ model under realistic experimental conditions. As a proof-of-concept, our work also shows that we may fruitfully apply the ADM to any number of competing semi-classical gravity models [78][79][80]. This may be useful for a range of experimental tests, including tests of gravitationally-induced wave function collapse [81][82][83].

VI. DISCUSSIONS AND FINAL REMARKS
In the present paper, we have investigated the semi-analytical series solutions of the time-dependent Schrödinger-Newton-Λ (S-N-Λ) system, which describes quantum matter in the presence of a nonlinear self- gravitational interaction and a background dark energy density. For the latter, we adopted for the simple form of a positive cosmological constant, which enters into the mathematical formalism through the modified Poisson equation. In order to solve the coupled system of S-N-Λ equations, we used a powerful mathematical method called the Adomian Decomposition Method (ADM), which provides in a fast and efficient way of obtaining series solutions of strongly nonlinear differential equations. The starting point of this method is the transformation of the given system of differential equations into an equivalent system of integral equations. Then, by positing the existence of series solutions of the integral system, one can obtain sets of recurrence relations for each unknown term in the power series expansion.
Usually, the ADM series converges fast, allowing detailed studies of the solutions of highly nonlinear differential equations using purely analytical methods. The main advantage of the method outlined in this paper is that it is based on a rigorous mathematical procedure, namely, the series expansions of the wave function, and of the nonlinear self-gravity term, while at the same time providing results that are mathematically simple and physically intuitive. This allows the in-depth investigation of the role dark energy may play in the microscopic dynamics of a quantum particle.
In the cosmological context, the dark energy density can be inferred from the critical density of the Universe, given by ρ cr = 3H 2 0 /8πG = 1.88h 2 × 10 −29 g/cm 3 , where H 0 is the present day value of the Hubble constant, and h = H 0 /100 km s −1 Mpc −1 . Since the cosmological data indicates a dark energy density of the order of ρ vac ≃ 0.75ρ cr , it follows that ρ vac ≃ 10 −29 g/cm 3 . On the other hand, the cosmological dark energy can be obtained from physical considerations, once it is interpreted as a vacuum energy, as ρ vac = √ k Pl k dS k dS k 2 + (mc/ ) 2 dk, where k Pl = 2π/l Pl , and k dS = 2π/l dS , where l Pl = G/c 3 is the Planck length and l dS = 3/Λ is the de Sitter length [84][85][86]. This is consistent with the existence of the GUP and EUP [87] and with the recent tentative observational evidence for the granular nature of dark energy on scales of order (k Pl k dS ) −1/2 ≃ 0.1 mm [88][89][90][91][92].
In quantum physics, a quantum fluctuation (also called vacuum fluctuation), is the random variation of the energy at a point in space, due to the creation of virtual particle-antiparticle pairs. These pairs are continuously created in the space, according to the energy-time uncer-tainty principle, ∆E∆t ≥ /2. In our present approach, we describe the effects of these processes on the quantum dynamics of the particle via a constant term. Even though, on a cosmological scale, the vacuum energy may have a very low (but extremely important) numerical value, quantum fluctuations may still have a significant impact on the local particle dynamics, at a microscopic level, over sufficiently long time-scales [54]. However, as a future extension of our current work, it would be interesting to reanalyze the problem using an alternative dark energy ansatz, which captures the oscillating, or 'granular' nature of the dark energy density proposed in recent models [84][85][86][87][88][89][90][91][92].
With or without a dark energy term, the importance of the self-gravitational interaction essentially depends on the mass of the particle. For a particle with a mass of the order of m = 10 10 m p ≃ 10 −14 g, where m p is the proton mass, the dimensionless coefficient σ given by Eq. (71) is of order unity, σ ≃ 1. In this regime, the selfgravitational interaction has a significant effect on the evolution of the quantum wave packet. In the absence of dark energy, λ = 0, it follows from Eq. (72) that this phase corresponds to the standard gravity-dominated regime of the Schrödinger-Newton system.
In summary, the consistency of the Adomian series solutions with the exact numerical solutions obtained in previous studies represents a huge step forward in the study of the S-N-Λ system. Up to now, only very crude and approximate analytical methods could be used to investigate its dynamics. Although useful for developing our physical intuition and providing order-of-magnitude estimates, these are no substitute for accurate quantitative solutions. Conversely, obtaining accurate numerical solutions is resource intensive, requiring long periods of time to develop and run the relevant codes, which are also computationally demanding [54].
By contrast, the same results can be obtained using Adomian decomposition in a fraction of the time, with the help of a relatively simple Mathematica or Maple worksheet. Indeed, in [54], it was stated that "we must deal with a complicated integro-differential equation, with little hope for analytical exploration". We have now shown that this is not the case and that the S-N-Λ system may be investigated analytically, to any degree of desired accuracy, using the right series solution techniques. By applying the Adomian decomposition method to PDEs, it should even be possible to obtain non-spherically symmetric solutions of the S-N-Λ equations. To the best of the author's knowledge, this has not yet been attempted in the existing literature, even numerically.
The preliminary results presented here indicate that the Adomian Decomposition Method can be used to obtain accurate solutions of a wide variety of semi-classical gravity models, subject to a wide range of initial conditions. Ultimately, this should help us to test the predictions of these models in greater detail, under realistic experimental conditions [77].