On the numerical experiments of the Cauchy problem for semi-linear Klein-Gordon equations in the de Sitter spacetime
OOn the numerical experiments of the Cauchy problemfor semi-linear Klein-Gordon equations in the de Sitterspacetime
Takuya Tsuchiya a, ∗ , Makoto Nakamura b, ∗∗ a Center for Liberal Arts and Sciences, Hachinohe Institute of Technology, 88-1, Obiraki,Myo, Hachinohe, Aomori 031-0814, JAPAN. b Faculty of Science, Yamagata University, Kojirakawa-machi 1-4-12, Yamagata 990-8560,JAPAN.
Abstract
The computational analysis of the Cauchy problem for semi-linear Klein-Gordonequations in the de Sitter spacetime is considered. Several simulations are per-formed to show the time-global behaviors of the solutions of the equations inthe spacetime based on the structure-preserving scheme. It is remarked thatthe sufficiently large Hubble constant yields the strong diffusion-effect whichgives the long and stable simulations for the defocusing semi-linear terms. Thereliability of the simulations is confirmed by the preservation of the numericallymodified Hamiltonian of the equations.
Keywords: semilinear Klein-Gordon equation, Cauchy problem, de Sitterspacetime, structure-preserving scheme
1. Introduction
The mathematical structure of partial differential equations in non-flat space-times is changed by the variations of curvatures since the differential operatorsare influenced by the metrics of the spacetimes. The de Sitter spacetime is oneof the solutions of the Einstein gravitational equations with the cosmologicalconstant in the vacuum. The space is expanding or contracting along the timein this spacetime. To study the effects of the spatial variation, we considerthe semilinear field equation of the Klein-Gordon type, and we carry out thenumerical experiments of the solutions and their energy. To investigate the par-tial differential equations with computational analysis, we need to perform highprecise and accurate simulations. There are many methods to make discretized ∗ [email protected] ∗∗ [email protected] Preprint submitted to Journal of Computational and Applied Mathematics May 23, 2019 a r X i v : . [ g r- q c ] M a y quations. The Adomian decomposition method [1] can provide analytical ap-proximation solutions of linear and nonlinear differential equations. Since thismethod dose not require the linearization, it would be appropriate to makedicretization of the nonlinear differential equations. It is reported that the sim-ulations of the nonlinear Klein-Gordon equation in the flat spacetime with themethod [5]. In [31], the differential transform method [32] and the variational it-eration method [24] are used to perform the nonlinear Klein-Gordon equation inthe de Sitter spacetime. The former is one of the power series expansions. Thevirtue is easy to estimate the numerical errors because the truncation higher or-der terms are corresponding to the dominant numerical errors. Since the latteris classified as the method of Lagrange multipliers, it would be close to the ex-act solutions by repeating the iteration if the appropriate conditions are given.In addition, the results with the explicit fourth-order Runge-Kutta method arereported in [4]. The Runge-Kutta method is well-known discretized scheme forthe nonlinear differential equations.The above methods are generic discretization schemes for the differentialequations to calculate simulations. Therefore, the properties of the differentialequations in continuous level often lost in the discretization, and it makes thenumerical errors in the evolutions. On the other hand, we use the structure-preserving scheme (SPS) called as “discrete variational derivative method” inthis paper. The concept of this method is to conserve the structures andthe properties of the equations in continuous level. This method is widelyused to derive the discrete equations of partial differential equations (see e.g.,[8, 9, 22]). To make discretized equations using SPS, we need the Lagrangianor the Hamiltonian[12, 25]. This is because the derivations of the discretizedequations with SPS are similar with the process of making the continuous equa-tions using the variational principle. In comparison with the construction of theLagrangian formulation, the one of the Hamiltonian formulation is clear. Espe-cially, since the Hamiltonian means the total energy of the system, this value isregarded as one of the constraints of the system. By monitoring the value in theevolution, we can judge whether the numerical simulations are success or not.Conservation of the constraints is one of the necessary conditions to performsuccessive simulations. If the constraints dose not conserve in the evolution, thesimulations fail (e.g. [21]). Therefore, we use the Hamiltonian formulation inthis paper.We start from the introduction of the de Sitter spacetime. Let n ≥ c > α, β, γ, . . . run from 0 to n , Latin letters j, k, (cid:96), . . . run from 1 to n . Let g αβ dx α dx β be the metric in R n . We denote by ( g αβ ) ≤ α,β ≤ n the (1+ n ) × (1+ n )-matrix whose ( α, β )-component is given by g αβ . Put g := det( g αβ ) ≤ α,β ≤ n ,and let ( g αβ ) ≤ α,β ≤ n be the inverse matrix of ( g αβ ) ≤ α,β ≤ n . We use the Einsteinrule for the sum of indices of the tensors, for example, T αα := (cid:80) nα =0 T αα and T ii := (cid:80) ni =1 T ii . The change of upper and lower indices of the tensors is doneby g αβ and g αβ , for example, T αβ := g αγ T γβ .For a stress-energy tensor T αβ , the (1 + n )-dimensional Einstein equation is2efined by G αβ + Λ g αβ = µ T αβ , (1.1)where G αβ is the Einstein tensor, Λ is the cosmological constant, and µ is theEinstein gravitational constant. When we consider the universe filled with theperfect fluid of the mass density ρ and the pressure p , we are able to use thestress-energy tensor T αβ of the perfect fluid given by T αβ := diag( − ρc , p, . . . , p ) . When the cosmological term Λ g αβ is transposed to right-hand-side in (1.1), theterm is regarded as a part of the stress-energy tensor. Then the term − Λ g αβ isrewritten as − Λ g αβ = µ diag( − (cid:101) ρc , (cid:101) p, . . . , (cid:101) p ) , (cid:101) ρ := Λ µc , (cid:101) p := − Λ µ , and the cosmological constant Λ > ds = g αβ dx α dx β = − c dt + e Ht δ ij dx i dx j , (1.2)where t := x , H is the Hubble constant, and δ ij denotes the Kronecker delta.The relation between the cosmological constant Λ and the Hubble constant H isgiven by H = (cid:112) / ( n − n ). We can write ( g αβ ) ≤ α,β ≤ n = diag( − c , e Ht , · · · , e Ht ).Next, we introduce the semilinear Klein-Gordon equation (KGE) in curvedspacetimes. KGE presents the equation of motion for the massive scalar field.Then, to derive the equation of motion described by a real-valued function φ with the mass m and the potential V , we consider the Lagrangian density L defined by L := − √− g (cid:26) g αβ ( ∇ α φ )( ∇ β φ ) + c m (cid:126) φ + 2 V ( φ ) (cid:27) , (1.3)where ∇ µ is the covariant derivative associated with g µν and (cid:126) is the Planckconstant. If the spacetime is flat, the Lagrangian density L is consistent withthe well-known Lagrangian density of KGE (see e.g., [25]). The Euler-Lagrangeequation of the action (cid:82) R n L dtdx gives KGE in curved spacetime − √− g ∂ α ( √− gg αβ ∂ β φ ) + m c (cid:126) φ + V (cid:48) ( φ ) = 0 . In the de Sitter spacetime, KGE is rewritten as ∂ t φ + nH∂ t φ − c e − Ht δ ij ∂ i ∂ j φ + m c (cid:126) φ + c V (cid:48) ( φ ) = 0 , (1.4)3here V (cid:48) denotes the derivative V (cid:48) ( r ) = dV ( r ) /dr for r ∈ R (see, e.g., [7, 15,26]). Multiplying ∂ t φ to the both sides of (1.4), we have the divergence formula ∂ t e + ∂ j e j + e ∗ = 0 , where e α , 0 ≤ α ≤ n , and e ∗ mean the energy density and the external energy,respectively. They are defined by e := 12 (cid:26) ( ∂ t φ ) + c e − Ht δ ij ( ∂ i φ )( ∂ j φ ) + m c (cid:126) φ + 2 c V ( φ ) (cid:27) ,e j := − c e − Ht δ ij ( ∂ i φ )( ∂ t φ ) ,e ∗ := H { n ( ∂ t φ ) + c e − Ht δ ij ( ∂ i φ )( ∂ j φ ) } . We define the total energy by E ( t ) := (cid:90) R n e ( t, x ) dx, (1.5)then we have the energy estimate E ( t ) + (cid:90) t (cid:90) R n e ∗ ( s, x ) dx ds = E (0) (1.6)for t ∈ [0 , ∞ ). Let us consider the defocusing potential ( V ≥ E is a positive-valued function for any non-trivial solution φ . When the space isexpanding ( H > E ( t )
0. On the other hand, when the space is contracting (
H < E ( t ) > E (0) for t > V of power type and its derivative V (cid:48) given by V ( φ ) := λp + 1 | φ | p +1 , V (cid:48) ( φ ) = λ | φ | p − φ (1.7)for λ ∈ R and 1 < p < ∞ . We consider the Cauchy problem (cid:0) ∂ t φ + nH∂ t φ − c e − Ht δ ij ∂ i ∂ j φ + m c φ/ (cid:126) + c λ | φ | p − φ (cid:1) ( t, x ) = 0for ( t, x ) ∈ [0 , T ) × R n φ (0 , · ) = φ ( · ) ∈ H ( R n ) , ∂ t φ (0 , · ) = φ ( · ) ∈ L ( R n ) , (1.8)where φ , φ are given initial data, H ( R n ) denotes the Sobolev space and L ( R n ) denotes the Lebesgue space. On the Cauchy problem (1.8) for thepositive Hubble constant H >
0, the following theoretical results have beenshown. The first result shows the existence of the time-local solutions for anyinitial data. The second result shows that the local solutions are time-globalsolutions if the initial data are sufficiently small. The third result shows thelocal solutions are global for any data if the semilinear term is defocusing.4 heorem 1.1. ([15, Theorems 1.2 and 1.7]) Let
H > . Assume m > nH/ .Let p satisfy ≤ p (cid:26) < ∞ if n = 1 , ≤ n − if n ≥ . (1.9) Then we have the following results.(1) For any φ and φ , there exists T = T ( (cid:107) φ (cid:107) H ( R n ) + (cid:107) φ (cid:107) L ( R n ) ) > suchthat (1.8) has a unique solution φ in C ([0 , T ) , H ( R n )) ∩ C ([0 , T ) , L ( R n )) .(2) Let n ≤ . If (cid:107) φ (cid:107) H ( R n ) + (cid:107) φ (cid:107) L ( R n ) is sufficiently small, and /n ≤ p , then (1.8) has a unique solution φ in C ([0 , ∞ ) , H ( R n )) ∩ C ([0 , ∞ ) , L ( R n )) .(3) If λ ≥ , then (1.8) has a unique global solution φ in C ([0 , ∞ ) , H ( R n )) ∩ C ([0 , ∞ ) , L ( R n )) for any data φ ∈ H ( R n ) and φ ∈ L ( R n ) . We refer to the corresponding known results on (1.8). D’Ancona and Giuseppehave shown in [2] and [3] global classical solutions for ( ∂ t − a ( t ) δ ij ∂ i ∂ j ) u + | u | p − u = 0 with some additional conditions on a ( t ) ≥ p when n = 1 , , n when the nonlinear term f is of power type p >
1, and the norm of initialdata (cid:107) u (cid:107) H s ( R n ) + (cid:107) u (cid:107) H s ( R n ) is sufficiently small for some s > n/ ≥ s = 1) in Theorem 1.1 have been shown. Thisresult has been extended to the case of general Friedmann-Lemaˆıtre-Robertson-Walker spacetime in [10]. Baskin has shown in [7] small global solution for( g µν ∇ µ ∇ ν + λ ) u + f ( u ) = 0 when f ( u ) is a type of | u | p − u , p = 1 + 4 / ( n − λ > n /
4, ( u , u ) ∈ H ⊕ L , where g µν gives the asymptotic de Sitter space-time (see also [6] for the cases p = 5 with n = 3, p = 3 with n = 4). Blow-upphenomena are considered in [27]. See also the references in the summary [30]by Yagdjian.In Theorem 1.1, we have assumed the condition m > nH/ nH∂ t φ , we use the transformation u = e κt φ for κ ∈ R to transform theequation to the Klein-Gordon equation. Then the equation is rewritten as (cid:26) ∂ t + ( nH − κ ) ∂ t − c e − Ht δ ij ∂ i ∂ j + m c (cid:126) + κ ( κ − nH ) (cid:27) u + c e κt V (cid:48) ( e − κt u ) = 0 , (1.10)and we obtain the equation( ∂ t − c e − Ht δ ij ∂ i ∂ j + M ) u ( t, x ) + c e nHt/ V (cid:48) ( e − nHt/ u ( t, x )) = 0 (1.11)when κ = nH/ M := c m / (cid:126) − n H /
4. We assume m is large such that M > κ < nH/ κ ≤ nH/ κ ≥ nH/ κ = nH/
2, by which we lose the5issipative term ( nH − κ ) ∂ t u in (1.11). Since the first equation in (1.8) has theterm e − Ht δ ij ∂ i ∂ j φ with the non-constant coefficient, it is not easy to obtainthe critical theoretical results so far. Instead, in this paper, we carry out somenumerical simulations, and show the detailed behaviors to clarify the dissipativeeffect of the spatial expansion.
2. Hamiltonian formulation of KGE
Let us consider the Hamiltonian density for the Lagrangian density L . TheLagrangian density of KGE (1.3) in the de Sitter spacetime becomes L = − c e nHt (cid:26) − c ( ∂ t φ ) + e − Ht δ ij ( ∂ i φ )( ∂ j φ ) + c m (cid:126) φ + 2 V ( φ ) (cid:27) . (2.1)Then, we define the momentum ψ by ψ := ∂ L ∂ ( ∂ t φ ) = e nHt c ∂ t φ, the Hamiltonian density H by the Legendre transformation H := ψ ( ∂ t φ ) − L = ce nHt (cid:26) e − nHt ψ + e − Ht δ ij ( ∂ i φ )( ∂ j φ ) + c m (cid:126) φ + 2 V ( φ ) (cid:27) , and the Hamiltonian H C by H C ( t ) := (cid:90) R n H ( t, x ) dx (2.2)for t ∈ R . To investigate the properties of H C , we denote the kinetic term,the diffusion term, the mass term and the nonlinear term by the integration of( ∂ t φ ) / c e − Ht δ ij ( ∂ i φ )( ∂ j φ ) / m c φ / (2 (cid:126) ) and c V ( φ ). Namely,the kinetic term : K ( t ) := (cid:90) R n
12 ( ∂ t φ ( t, x )) dx, the diffusion term : D ( t ) := (cid:90) R n c e − Ht δ ij ( ∂ i φ ( t, x ))( ∂ j φ ( t, x )) dx, the mass term : M ( t ) := (cid:90) R n m c (cid:126) φ ( t, x ) dx, the nonlinear term : N ( t ) := (cid:90) R n c V ( φ ( t, x )) dx = (cid:90) R n λc p + 1 | φ | p +1 dx (2.3)for t ∈ R . With (2.3), the total energy (1.5) can be expressed as E ( t ) = K ( t ) + D ( t ) + M ( t ) + N ( t ) . (2.4)6ince the relation between H and the energy density e becomes H = e nHt c e , (2.5) H C can be expressed as H C ( t ) = e nHt c ( K ( t ) + D ( t ) + M ( t ) + N ( t )) . If the spacetime is flat ( H = 0) and the speed of light c is set as 1, then H C is consistent with E . The Hamiltonian formulation of (1.4) is given by thecanonical equations ∂ t φ := δ H δψ = ce − nHt ψ,∂ t ψ := − δ H δφ = ce ( n − Ht δ ij ( ∂ i ∂ j φ ) − c m (cid:126) e nHt φ − ce nHt V (cid:48) ( φ ) . (2.6)Since H satisfies the differential equation ∂ t H = nH H + c e − Ht ∂ i { ψδ ij ( ∂ j φ ) } − H ∗ , where we have put H ∗ = H { cne − nHt ψ + ce ( n − Ht δ ij ( ∂ i φ )( ∂ j φ ) } ,H C satisfies the equation ∂ t H C ( t ) = nHH C ( t ) − (cid:90) R n H ∗ dx, (2.7)and by using (2.3), (2.7) can be expressed as ∂ t H C ( t ) = Hc e nHt {− nK ( t ) + ( n − D ( t ) + nM ( t ) + nN ( t ) } . (2.8)The equation (2.7) (or (2.8)) indicates H C is not conserved in time evolution ifthe Hubble constant H is not zero. Therefore, we define ˜ H C as˜ H C ( t ) := H C ( t ) − (cid:90) t ∂ t H C ( s ) ds, (2.9)which is a constant value in time evolution independent of the value of H .The behaviors of E and H C show the properties of the effects of the spatialexpansion and contraction by the Hubble constant H . Namely, we are ableto understand the properties of the dark energy, which is equivalent to thecosmological constant mathematically, through the asymptotic behaviors of φ , E and H C . However, it is not easy to obtain the detailed behaviors of them bytheoretical methods since the waves of φ would be oscillating and interfered bythe semilinear term V ( φ ), and (1.4) is the equation with the variable coefficient.Thus, we perform the numerical simulations on φ , E and H C to study thedetailed dissipative effects of the spatial expansion.7 . Discrete KGE To investigate KGE by computational analysis, we need the high preciseand accurate numerical data. Then, we use SPS to get the high precise andaccurate numerical results. With this scheme, the structures of the equationsin the continuous case are preserved in the discrete case. In this paper, we usethis scheme to make the discrete Hamiltonian formulation. The discrete valuesof the variable u are defined as u ( (cid:96) )( k ) := u ( (cid:96) ∆ t, k ∆ x ), where (cid:96) is the time index,∆ t is the time mesh width, k is the spatial grid index, and ∆ x is the spatialmesh width. The details of this are in Ref.[9].For comparison with SPS, we use the Crank-Nicolson scheme (CNS) and thefourth-order Runge-Kutta scheme (RKS) which are widely used for calculatingpartial differential equations numerically. By using SPS, we can get a set of discrete evolution equations of KGE. Wegive a discrete Hamiltonian density H ( (cid:96) )( k ) as H ( (cid:96) )( k ) = ce nHt (cid:96) (cid:26) e − nHt (cid:96) ( ψ ( (cid:96) )( k ) ) + e − Ht (cid:96) δ ij ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) )( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) )( k ) ) + c m (cid:126) ( φ ( (cid:96) )( k ) ) + 2 λp + 1 | φ ( (cid:96) )( k ) | p +1 (cid:27) , (3.1)where t (cid:96) presents the time of the (cid:96) th step and (cid:98) δ (cid:104) (cid:105) i u ( (cid:96) )( k ) is the centered spaceoperator defined by (cid:98) δ (cid:104) (cid:105) i u ( (cid:96) )( k ) := u ( (cid:96) )( k +1) − u ( (cid:96) )( k ) + u ( (cid:96) )( k − x i . For H ( (cid:96) )( k ) , φ ( (cid:96) )( k ) and ψ ( (cid:96) )( k ) , we define the values (cid:98) δ H / (cid:0)(cid:98) δ ( φ ( (cid:96) +1)( k ) , φ ( (cid:96) )( k ) ) (cid:1) and (cid:98) δ H / (cid:0)(cid:98) δ ( ψ ( (cid:96) +1)( k ) , ψ ( (cid:96) )( k ) ) (cid:1) by the equation H ( (cid:96) +1)( k ) − H ( (cid:96) )( k ) = (cid:98) δ H (cid:98) δ ( φ ( (cid:96) +1)( k ) , φ ( (cid:96) )( k ) ) ( φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ) + (cid:98) δ H (cid:98) δ ( ψ ( (cid:96) +1)( k ) , ψ ( (cid:96) )( k ) ) ( ψ ( (cid:96) +1)( k ) − ψ ( (cid:96) )( k ) ) . (3.2)8hen, a set of discrete KGE with SPS becomes φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ∆ t = (cid:98) δ H (cid:98) δ ( ψ ( (cid:96) +1)( k ) , ψ ( (cid:96) )( k ) )= 14 c ( e − nHt (cid:96) +1 + e − nHt (cid:96) )( ψ ( (cid:96) +1)( k ) + ψ ( (cid:96) )( k ) ) ,ψ ( (cid:96) +1)( k ) − ψ ( (cid:96) )( k ) ∆ t = − (cid:98) δ H (cid:98) δ ( φ ( (cid:96) +1)( k ) , φ ( (cid:96) )( k ) )= 14 c ( e ( n − Ht (cid:96) +1 + e ( n − Ht (cid:96) ) δ ij (cid:98) δ (cid:104) (cid:105) j (cid:98) δ (cid:104) (cid:105) i ( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) − c m (cid:126) ( e nHt (cid:96) +1 + e nHt (cid:96) )( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) − λc p + 1) ( e nHt (cid:96) +1 + e nHt (cid:96) ) | φ ( (cid:96) +1)( k ) | p +1 − | φ ( (cid:96) )( k ) | p +1 φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) . (3.3)The discrete evolution equations (3.3) are corresponding to (2.6) for the contin-uous case. The nonlinear term of the last term in (3.3) is expressed as | φ ( (cid:96) +1)( k ) | p +1 − | φ ( (cid:96) )( k ) | p +1 φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) = {| φ ( (cid:96) +1)( k ) | p + | φ ( (cid:96) +1)( k ) | p − | φ ( (cid:96) )( k ) | + · · · + | φ ( (cid:96) +1)( k ) || φ ( (cid:96) )( k ) | p − + | φ ( (cid:96) )( k ) | p } | φ ( (cid:96) +1)( k ) | − | φ ( (cid:96) )( k ) | φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) (3.4)when p is a natural number. We refer to the expression in the nonlinearSchr¨odinger equations in Ref.[9] to make above relation. With (3.3), the time9ifference of H ( (cid:96) )( k ) is calculated as H ( (cid:96) +1)( k ) − H ( (cid:96) )( k ) ∆ t = c { ( ψ ( (cid:96) +1)( k ) ) + ( ψ ( (cid:96) )( k ) ) } e − nHt (cid:96) +1 − e − nHt (cid:96) ∆ t + c e − nHt (cid:96) +1 + e − nHt (cid:96) )( ψ ( (cid:96) +1)( k ) + ψ ( (cid:96) )( k ) ) ψ ( (cid:96) +1)( k ) − ψ ( (cid:96) )( k ) ∆ t + c δ ij { ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) +1)( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) +1)( k ) ) + ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) )( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) )( k ) ) } e ( n − Ht (cid:96) +1 − e ( n − Ht (cid:96) ∆ t + c e ( n − Ht (cid:96) +1 + e ( n − Ht (cid:96) ) δ ij (cid:98) δ (cid:104) (cid:105) i ( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ∆ t + c m (cid:126) { ( φ ( (cid:96) +1)( k ) ) + ( φ ( (cid:96) )( k ) ) } e nHt (cid:96) +1 − e nHt (cid:96) ∆ t + c m (cid:126) ( e nHt (cid:96) +1 + e nHt (cid:96) )( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ∆ t + cλ p + 1) ( | φ ( (cid:96) +1)( k ) | p +1 + | φ ( (cid:96) )( k ) | p +1 ) e nHt (cid:96) +1 − e nHt (cid:96) ∆ t + cλ p + 1) ( e nHt (cid:96) +1 + e nHt (cid:96) ) | φ ( (cid:96) +1)( k ) | p +1 − | φ ( (cid:96) )( k ) | p +1 φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ∆ t = c { ( ψ ( (cid:96) +1)( k ) ) + ( ψ ( (cid:96) )( k ) ) } e − nHt (cid:96) +1 − e − nHt (cid:96) ∆ t + c
16 ( e − nHt (cid:96) +1 + e − nHt (cid:96) )( e ( n − Ht (cid:96) +1 + e ( n − Ht (cid:96) ) · δ ij (cid:98) δ (cid:104) (cid:105) i { ( ψ ( (cid:96) +1)( k ) + ψ ( (cid:96) )( k ) ) (cid:98) δ (cid:104) (cid:105) j ( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) } + c δ ij { ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) +1)( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) +1)( k ) ) + ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) )( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) )( k ) ) } e ( n − Ht (cid:96) +1 − e ( n − Ht (cid:96) ∆ t + c m (cid:126) { ( φ ( (cid:96) +1)( k ) ) + ( φ ( (cid:96) )( k ) ) } e nHt (cid:96) +1 − e nHt (cid:96) ∆ t + cλ p + 1) ( | φ ( (cid:96) +1)( k ) | p +1 + | φ ( (cid:96) )( k ) | p +1 ) e nHt (cid:96) +1 − e nHt (cid:96) ∆ t − cnH { ( ψ ( (cid:96) +1)( k ) ) + ( ψ ( (cid:96) )( k ) ) } e − nHt (cid:96) + (cid:98) δ (cid:104) (cid:105) i (cid:26) c
16 ( e − nHt (cid:96) +1 + e − nHt (cid:96) )( e ( n − Ht (cid:96) +1 + e ( n − Ht (cid:96) ) · { ( ψ ( (cid:96) +1)( k ) + ψ ( (cid:96) )( k ) ) δ ij (cid:98) δ (cid:104) (cid:105) j ( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) } (cid:27) + c ( n − H δ ij { ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) +1)( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) +1)( k ) ) + ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) )( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) )( k ) ) } e ( n − Ht (cid:96) + c m nH (cid:126) { ( φ ( (cid:96) +1)( k ) ) + ( φ ( (cid:96) )( k ) ) } e nHt (cid:96) + cλnH p + 1) ( | φ ( (cid:96) +1)( k ) | p +1 + | φ ( (cid:96) )( k ) | p +1 ) e nHt (cid:96) + O (∆ t ) , (3.5)where we use the relation of t (cid:96) +1 = t (cid:96) + ∆ t . We define the discrete Hamiltonian H C ( (cid:96) ) by H C ( (cid:96) ) := (cid:88) k ∈ Z H ( (cid:96) )( k ) (cid:89) ≤ i ≤ n ∆ x i , (3.6)then (3.5) is calculated as H C ( (cid:96) +1) − H C ( (cid:96) ) ∆ t = H c e nHt (cid:96) {− n ( K ( (cid:96) +1) + K ( (cid:96) ) ) + ( n − D ( (cid:96) +1) + D ( (cid:96) ) )+ n ( M ( (cid:96) +1) + M ( (cid:96) ) ) + n ( N ( (cid:96) +1) + N ( (cid:96) ) ) } + O (∆ t ) , (3.7)where K ( (cid:96) ) , D ( (cid:96) ) , M ( (cid:96) ) and N ( (cid:96) ) are defined by K ( (cid:96) ) := (cid:88) k ∈ Z c e − nHt (cid:96) ( ψ ( (cid:96) )( k ) ) (cid:89) ≤ i ≤ n ∆ x i ,D ( (cid:96) ) := (cid:88) k ∈ Z c e − Ht (cid:96) δ ij ( (cid:98) δ (cid:104) (cid:105) i φ ( (cid:96) )( k ) )( (cid:98) δ (cid:104) (cid:105) j φ ( (cid:96) )( k ) ) (cid:89) ≤ i ≤ n ∆ x i ,M ( (cid:96) ) := (cid:88) k ∈ Z m c (cid:126) ( φ ( (cid:96) )( k ) ) (cid:89) ≤ i ≤ n ∆ x i ,N ( (cid:96) ) := (cid:88) k ∈ Z c V ( φ ( (cid:96) )( k ) ) (cid:89) ≤ i ≤ n ∆ x i = (cid:88) k ∈ Z λc p + 1 | φ ( (cid:96) )( k ) | p +1 (cid:89) ≤ i ≤ n ∆ x i . (3.8) K ( (cid:96) ) , D ( (cid:96) ) , M ( (cid:96) ) and N ( (cid:96) ) are discrete values corresponding to K ( t ), D ( t ), M ( t )and N ( t ) defined by (2.3), respectively. We set˜ H C ( (cid:96) ) := H C ( (cid:96) ) − (cid:88) ≤ i ≤ (cid:96) − H c e nHt (cid:96) {− n ( K ( i +1) + K ( i ) ) + ( n − D ( i +1) + D ( i ) )+ n ( M ( i +1) + M ( i ) ) + n ( N ( i +1) + N ( i ) ) } , (3.9)then ˜ H C ( n ) is a discretized value of ˜ H C defined by (2.9). This value is a guidelineto perform precise numerical simulations in the cases of H (cid:54) = 0.11 .2. CNS and RKS CNS and RKS are frequently used to perform partial differential equations.The discrete equations of (2.6) with CNS are given by φ ( (cid:96) +1)( k ) − φ ( (cid:96) )( k ) ∆ t = Φ (cid:16) ( ψ ( (cid:96) +1)( k ) + ψ ( (cid:96) )( k ) ) / , ( t (cid:96) +1 + t (cid:96) ) / (cid:17) ,ψ ( (cid:96) +1)( k ) − ψ ( (cid:96) )( k ) ∆ t = Ψ (cid:16) ( φ ( (cid:96) +1)( k ) + φ ( (cid:96) )( k ) ) / , ( t (cid:96) +1 + t (cid:96) ) / (cid:17) , (3.10)where the right hand sides of (3.10) are defined asΦ( ψ, t ) := ce − nHt ψ, Ψ( φ, t ) := ce ( n − Ht δ ij (cid:98) δ (cid:104) (cid:105) ij φ − c m (cid:126) e nHt φ − λce nHt | φ | p − φ, (3.11)and the second-order central difference operator (cid:98) δ (cid:104) (cid:105) ij is defined as (cid:98) δ (cid:104) (cid:105) ij u ( (cid:96) )( k ) := (cid:40) ( u ( (cid:96) )( k +1) − u ( (cid:96) )( k ) + u ( (cid:96) )( k − ) / (∆ x i ) , ( i = j ) (cid:98) δ (cid:104) (cid:105) i (cid:98) δ (cid:104) (cid:105) j u ( (cid:96) )( k ) . ( i (cid:54) = j ) (3.12)The discrete equations with RKS are give by φ ( (cid:96) +1)( k ) = φ ( (cid:96) )( k ) + ( h + 2 h + 2 h + h )∆ t/ ,ψ ( (cid:96) +1)( k ) = ψ ( (cid:96) )( k ) + ( s + 2 s + 2 s + s )∆ t/ , (3.13)where h = Φ( ψ ( (cid:96) )( k ) , t (cid:96) ) ,s = Ψ( φ ( (cid:96) )( k ) , t (cid:96) ) ,h = Φ( ψ ( (cid:96) )( k ) + s ∆ t/ , t (cid:96) + ∆ t/ ,s = Ψ( φ ( (cid:96) )( k ) + h ∆ t/ , t (cid:96) + ∆ t/ ,h = Φ( ψ ( (cid:96) )( k ) + s ∆ t/ , t (cid:96) + ∆ t/ ,s = Ψ( φ ( (cid:96) )( k ) + h ∆ t/ , t (cid:96) + ∆ t/ ,h = Φ( ψ ( (cid:96) )( k ) + s ∆ t, t (cid:96) + ∆ t ) ,s = Ψ( φ ( (cid:96) )( k ) + h ∆ t, t (cid:96) + ∆ t ) . (3.14)
4. Numerical tests
In this section, we carry out the numerical simulations on φ , E and H C (or˜ H C ). Since the physical background of the Einstein equation is for the 1 + 3dimensional spacetime, we consider the case n = 3 in the following.The numerical settings are below. Set ( x , x , x ) = ( x, y, z ) ∈ R .12 Initial condition: φ = φ ( x, y, z ) = A cos(2 πx ) , φ = φ ( x, y, z ) = 2 πA sin(2 πx ) . • Numerical domains: 0 ≤ x ≤
1, 0 ≤ t ≤ • Boundary condition: periodic. • Grids: ∆ x = 1 / t = 1 / • Physical settings: mass m = 1, speed of light c = 1, Planck constant (cid:126) = 1. • The number of exponent in the nonlinear term: p = 2 , , , , φ ( t, x, y, z ) for ( t, x, y, z ) ∈ [0 , ∞ ) × R which is a con-stant along the variables y and z to give simple and reliable simulations. Moregeneral solutions and data will be considered in the future. We perform thesimulations by changing the values of the amplitude A in the initial condition,the coefficient value λ of the nonlinear term, and the Hubble constant H . First,we perform some simulations to study the accuracies with SPS. For compari-son with SPS, we use the Crank-Nicolson scheme (CNS) and the Runge-Kuttascheme (RKS) which are widely used for calculating partial differential equa-tions numerically. Because of the nonlinearity of KGE (2.6), the discretizedequations are calculated implicitly. Thus, the numerical simulations are per-formed iteratively. For SPS and CNS, the iterations are nine times and threetimes, respectively. If we calculate more times iteratively for each scheme, thenumerical results become worse in our codes. For the wave equation ∂ t φ − c δ ij ∂ i ∂ j φ + c V (cid:48) ( φ ) = 0with (1.7), namely, H = 0, m = 0 in (1.4) with (1.7), the scaling argument forthe solution φ ρ ( t, x ) := ρ / ( p − φ ( ρt, ρx ) for ρ ∈ R yields the critical exponentof p ( s ) := 1 + 4 / (3 − s ) for the well-posedness of the Cauchy problem of theequation in the Sobolev space H s ( R ) for s ∈ R . Especially, p ( − /
2) is calledthe Fujita exponent, p (0) is called the conformal exponent, p (1) is called theenergy critical exponent. This exponent p ( s ) also plays crucial roles for theKlein-Gordon equation, and it satisfies1 < p (cid:18) − (cid:19) = 53 < < p (0) = 73 < p (cid:18) (cid:19) = 3 < < p (1) = 5 < . By the simulations for p = 2 , , , ,
6, we are able to know the behaviors of thesolutions for p which is equivalent or close to these critical exponents. H = 0To compare with the numerical results of SPS, CNS and RKS, we set theHubble constant H as zero in this subsection. In this case, H C is consistent13 H C Time p =2 p =3 p =4 p =5 p =6 40060080010001200 0 200 400 600 800 1000 H C Time p =2 p =3 p =4 p =5 p =6 40060080010001200 0 200 400 600 800 1000 H C Time p =2 p =3 p =4 p =5 p =6 (1-t1): H C with SPS (1-t2): H C with CNS (1-t3): H C with RKS -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 (1-b1): | ( H C − H C (0)) /H C | with SPS (1-b2): | ( H C − H C (0)) /H C | with CNS (1-b3): | ( H C − H C (0)) /H C | with RKS Figure 1: The Hamiltonian H C , and the difference between H C and the initial data H C (0)for A = 4 and λ = 1. The horizontal axis is time. The top panels are the values of H C andthe bottom panels are the values of log | ( H C − H C (0)) /H C | . The left panels are drawn withSPS, the middle panels are with CNS, and the right panels are with RKS. with E , and H C is constant with time evolutions theoretically. To show theaccuracy of the simulations, we monitor H C . We draw H C , and the differencesbetween H C and the initial value H C (0) in the evolutions with SPS, CNS andRKS for A = 4 and λ = 1 in Fig.1. From the bottom panels (1-b1)–(1-b3), wesee the values of | ( H C − H C (0)) /H C | with CNS and RKS are larger than theones with SPS.Next, we show the results for A = 0 . λ = − t = 140. From the bottom panels, we see the values of | ( H C − H C (0)) /H C | with SPS are smaller than the ones with CNS and RKS. Incomparison with Fig.2, the simulations in Fig.1 are robust. Since it is necessarycondition of performing correct simulations that the value of | ( H C − H C (0)) /H C | is small, KGE have to be performed with SPS. Therefore, we use the numericalscheme as SPS hereafter. The results of φ with SPS are shown in Fig.3 andFig.4. Fig.3 shows the case of A = 4 and λ = 1. We see that the vibrationsoccur for p = 5 ,
6. Fig.4 shows the case of A = 0 . λ = −
1. We see all ofthe simulations stop before t = 140. H = 10 − In this subsection and next subsection, we perform the numerical simulationsfor
H >
0. First, we calculate for λ = 1 and H = 10 − . When H (cid:54) = 0, H C is notconstant in time evolutions. Since we cannot judge the reliability of simulationsby monitoring the values of H C , we show ˜ H C defined by (2.9) instead of H C .Fig.5 shows ˜ H C , and the differences between ˜ H C and the initial value ˜ H C (0)for A = 1 , ,
4. We see the value of p = 2 is the largest and that of p = 6is the smallest in the panel (5-t1). Contrarily, in the others of the top panels14 H C Time p =2 p =3 p =4 p =5 p =6 16.0816.1016.1216.1416.1616.1816.20 0 20 40 60 80 100 120 140 H C Time p =2 p =3 p =4 p =5 p =6 16.0816.1016.1216.1416.1616.1816.20 0 20 40 60 80 100 120 140 H C Time p =2 p =3 p =4 p =5 p =6 (2-t1): H C with SPS (2-t2): H C with CNS (2-t3): H C with RKS -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 20 40 60 80 100 120 140 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 20 40 60 80 100 120 140 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -16.00-14.00-12.00-10.00-8.00-6.00-4.00-2.000.00 0 20 40 60 80 100 120 140 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 (2-b1): | ( H C − H C (0)) /H C | with SPS (2-b2): | ( H C − H C (0)) /H C | with CNS (2-b3): | ( H C − H C (0)) /H C | with RKS Figure 2: The same as Fig.1 except for the values of A and λ . These results are set as A = 0 . λ = −
1. These calculations stop before t = 140. π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (3-1): A = 4 and p = 2. (3-2): A = 4 and p = 3. (3-3): A = 4 and p = 4. π π φ p =5 x Time φ π π φ p =6 x Time φ (3-4): A = 4 and p = 5. (3-5): A = 4 and p = 6. Figure 3: φ for A = 4 and λ = 1 with SPS. The vibrations for p = 5 , t = 500 and t = 100, respectively. π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (4-1): A = 0 . p = 2. (4-2): A = 0 . p = 3. (4-3): A = 0 . p = 4. π π φ p =5 x Time φ π π φ p =6 x Time φ (4-4): A = 0 . p = 5. (4-5): A = 0 . p = 6. Figure 4: The same as Fig.3 except for the values of A and λ . These results are set as A = 0 . λ = −
1. All of the simulations stop before t = 140. ˜ H C Time p =2 p =3 p =4 p =5 p =6 180200220240260280300 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 3004005006007008009001000110012001300 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 (5-t1): ˜ H C of A = 1 (5-t2): ˜ H C of A = 3 (5-t3): ˜ H C of A = 4 -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 (5-b1): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 1 (5-b2): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 3 (5-b3): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 4 Figure 5: ˜ H C , and the difference between ˜ H C and the initial data ˜ H C (0) for A = 1 , , λ = 1 and H = 10 − . The horizontal axis is time. The top panels show the values of ˜ H C , andthe bottom panels show the values of log | ( ˜ H C − ˜ H C (0)) / ˜ H C | . The left panels are the caseof A = 1, the middle panels are A = 3, and the right panels are A = 4. The lines of Panel(5-b1) are almost overlapping.
16e see the largest value is the case of p = 6. The reason is A = 1 or not. If A = 1, the nonlinear term becomes small as p becomes large. On the otherhand, if A >
1, the value becomes large as p becomes large. From the bottompanels (5-b1)–(5-b3), we see | ( ˜ H C − ˜ H C (0)) / ˜ H C | of the case p = 6 is the largestof them in each panel but all of the values in the bottom panels are enoughsmall. Fig.6 shows the total energy E defined by (1.5). We see that the energy -0.500.000.501.001.502.002.503.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -0.500.000.501.001.502.002.503.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -0.500.000.501.001.502.002.503.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 (6-1) E of A = 1 (6-2) E of A = 3 (6-3) E of A = 4 Figure 6: The total energy E for λ = 1 and H = 10 − . The horizontal axis is time, and thevertical axis is log ( E ). The left panel is the case of A = 1, the middle panel is A = 3 andthe right panel is A = 4. The lines of Panel (6-1) are almost overlapping. is exponential decay. Furthermore, the slopes of the lines in the panels of thefigure are almost the same. This result means the diffusion effect with timeevolutions is not caused by changing in A or p . While the lines in Panel (6-1)which is the case of A = 1 are almost overlapping, there are differences betweenthe values of E for A = 3 ,
4. To investigate the reason, we show the componentsof E defined by (2.3). -5-4-3-2-101 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN -5-4-3-2-101 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN -5-4-3-2-101 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN (7-1): A = 1 and p = 2 (7-2): A = 1 and p = 3 (7-3): A = 1 and p = 4 -5-4-3-2-101 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN -5-4-3-2-101 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN (7-4): A = 1 and p = 5 (7-5): A = 1 and p = 6 Figure 7: The kinetic term K , the diffusion term D , the mass term M , and the nonlinearterm N are shown for A = 1, λ = 1 and H = 10 − . The horizontal axis is time and thevertical axis is the logarithmic values of K , D , M and N . The kinetic term K , the diffusion term D , the mass term M and the non-linear term N for A = 1 are shown in Fig.7 and the ones for A = 4 are in Fig.8.17 l og ( v a l u e s ) Time
KDMN -2-10123 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN -2-10123 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN (8-1): A = 4 and p = 2 (8-2): A = 4 and p = 3 (8-3): A = 4 and p = 4 -2-10123 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN -2-10123 0 200 400 600 800 1000 l og ( v a l u e s ) Time
KDMN (8-4): A = 4 and p = 5 (8-5): A = 4 and p = 6 Figure 8: The same as Fig.7 except for the value of A . These results are set as A = 4. We see that K and D are the dominant terms, and M and N are very smallcompared with K and D in Fig.7. The decline of the energy-component N isbecoming faster than K , D and M as the power p increases. On the other hand,the proportion of N to E becomes large as p becomes large in Fig.8. This isbecause there are the differences between E s in Fig.6. Fig.9 shows φ for A = 4 π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (9-1): A = 4 and p = 2. (9-2): A = 4 and p = 3. (9-3): A = 4 and p = 4. π π φ p =5 x Time φ π π φ p =6 x Time φ (9-4): A = 4 and p = 5. (9-5): A = 4 and p = 6. Figure 9: The same as Fig.3 except for the value of H . These results are set as H = 10 − .The vibrations occur for p = 6 which is Panel (9-5). and λ = 1, and the numerical settings are the same as Fig.3 except for the valueof H . The vibrations occur in the only case of p = 6. Since the vibrations occurin the cases of p = 5 , φ . 18ext we show the results of the simulations with the case of λ = −
1. Fig.10 ˜ H C Time p =2 p =3 p =4 p =5 p =6 7.1607.1657.1707.1757.1807.1857.1907.1957.2007.205 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 16.0816.1016.1216.1416.1616.1816.20 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 (10-t1) ˜ H C of A = 0 . H C of A = 0 . H C of A = 0 . -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 -9.00-8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( ˜ H C − ˜ H C ( )) / ˜ H C | Time p =2 p =3 p =4 p =5 p =6 (10-b1) ( ˜ H C ( t ) − ˜ H C (0)) / ˜ H C ( t )of A = 0 . H C ( t ) − ˜ H C (0)) / ˜ H C ( t )of A = 0 . H C ( t ) − ˜ H C (0)) / ˜ H C ( t )of A = 0 . Figure 10: The same as Fig.5 except for the values of A and λ . These results are set as A = 0 . , . , . λ = −
1. The simulation for p = 2 and A = 0 . t = 200. shows ˜ H C and | ( ˜ H C − ˜ H C (0)) / ˜ H C | for A = 0 . , . , . λ = − H = 10 − .While we see the simulation for A = 0 . p = 2 stops before t = 200 in Panel(10-t3) and Panel (10-b3), the values | ( ˜ H C − ˜ H C (0)) / ˜ H C | for the other casesare enough small. Fig.11 shows E for A = 0 . , . , . λ = − H = 10 − . -2.00-1.50-1.00-0.500.000.501.001.50 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -2.00-1.50-1.00-0.500.000.501.001.50 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -2.00-1.50-1.00-0.500.000.501.001.50 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 (11-1) E of A = 0 . E of A = 0 . E of A = 0 . Figure 11: The same as Fig.6 except for the values of A and λ . These results are calculatedfor A = 0 . , . , . λ = −
1. The lines are almost overlapping in each panel except forthe case of p = 2 and A = 0 .
9. The simulation of p = 2 and A = 0 . t = 200. It seems there is no difference between the slopes of all the lines in Fig.6 andFig.11 regardless of the values of A , p and λ . Thus, the dissipation effect of E is not caused by changing in A , p and the signature of λ . Fig.12 shows φ for A = 0 . λ = − H = 10 − . These results are the same as Fig.4 exceptfor the value of H . By comparing Fig.4 with Fig.12, the simulation times for H = 10 − are longer than that for H = 0. It means the Hubble constant wouldinfluence the stability of the simulation. By comparing Fig.9 and Fig.12, thesimulations for λ = 1 are robust against the case of λ = − p = 2 and λ = − π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (12-1): A = 0 . p = 2 (12-2): A = 0 . p = 3 (12-3): A = 0 . p = 4 π π φ p =5 x Time φ π π φ p =6 x Time φ (12-4): A = 0 . p = 5 (12-5): A = 0 . p = 6. Figure 12: The same as Fig.4 except for the values of H . These results are set as H = 10 − .The simulation of the case p = 2 stops before t = 180. H = 0. H = 10 − In this subsection, we perform some simulations with H = 10 − and λ = 1.Fig.13 shows ˜ H C and | ( ˜ H C − ˜ H C (0)) / ˜ H C | for H = 10 − and this figure isthe same as Fig.5 except for the value of H . We see all the values of | ( ˜ H C − ˜ H C (0)) / ˜ H C | for H = 10 − are enough small. Fig.14 shows E with the Hubbleconstant as 10 − . There are few differences between the changes of p comparedwith the cases of H = 10 − in Fig.6. By comparing Fig.6 and Fig.14, thediffusion effect of E for H = 10 − is stronger than that for H = 10 − . Thus,the diffusion effect is mainly caused by the Hubble constant, and the effect wouldbecome strong as H becomes large. In Fig.15 which shows φ for H = 10 − , wesee the waveform is damped and there are few vibrations of the waveform. Bycomparing Fig.3, Fig.9 and Fig.15, the stability of the simulations becomes goodas H becomes large since the vibrations of φ decreases as H becomes large.Next we show the results for λ = −
1. In Fig.16, we see the simulations arereliable since all of the values of | ( ˜ H C − ˜ H C (0)) / ˜ H C | are enough small. Fig.17shows the total energy E . By comparing Fig.11 and Fig.17, we see the diffusioneffect for λ = − H = 10 − is stronger than that for λ = − H = 10 − .Thus, considering the results for λ = 1, the diffusion effect is mainly caused bynot the signature of the nonlinear term but the Hubble constant. Fig.18 shows φ for λ = − H = 10 − . The simulation times become long as H becomeslarge from the results in Fig.4, Fig.12 and Fig.18. These results indicate thesimulations become stable as H becomes large.20 ˜ H C Time p =2 p =3 p =4 p =5 p =6 180200220240260280300 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 300400500600700800900100011001200 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 (13-t1): ˜ H C of A = 1 (13-t2): ˜ H C of A = 3 (13-t3): ˜ H C of A = 4 -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 (13-b1): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 1 (13-b2): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 3 (13-b3): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 4 Figure 13: The same as Fig.5 except of the Hubble constant H . These results are set as H = 10 − . The lines in Panel (13-b1) are almost overlapping. -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.004.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.004.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.004.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 (14-1): E of A = 1 (14-2): E of A = 3 (14-3): E of A = 4 Figure 14: The same as Fig.6 except for the Hubble constant H . These results are set as H = 10 − . The lines in Panel (14-1) and Panel (14-2) are almost overlapping in each panel. π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (15-1): A = 4 and p = 2. (15-2): A = 4 and p = 3. (15-3): A = 4 and p = 4. π π φ p =5 x Time φ π π φ p =6 x Time φ (15-4): A = 4 and p = 5. (15-5): A = 4 and p = 6. Figure 15: The same as Fig.3 and Fig.9 except for the Hubble constant H . These results areset as H = 10 − . .7941.7951.7961.7971.7981.7991.800 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 7.1607.1657.1707.1757.1807.1857.1907.1957.2007.2057.210 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 16.0816.1016.1216.1416.1616.1816.20 0 200 400 600 800 1000 ˜ H C Time p =2 p =3 p =4 p =5 p =6 (16-t1): ˜ H C of A = 0 . H C of A = 0 . H C of A = 0 . -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 -8.00-7.00-6.00-5.00-4.00-3.00-2.00-1.000.00 0 200 400 600 800 1000 l og | ( H C − H C ( )) / H C | Time p =2 p =3 p =4 p =5 p =6 (16-b1): | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 0 . | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 0 . | ( ˜ H C − ˜ H C (0)) / ˜ H C | of A = 0 . Figure 16: The same as Fig.13 except for the value of λ . These results are set as λ = − -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 -14.00-12.00-10.00-8.00-6.00-4.00-2.000.002.00 0 200 400 600 800 1000 l og ( E ) Time p =2 p =3 p =4 p =5 p =6 (17-1): E of A = 0 . E of A = 0 . E of A = 0 . Figure 17: The same as Fig.11 except for H . These results are set as H = 10 − . The linesare almost overlapping each other. π π φ p =2 x Time φ π π φ p =3 x Time φ π π φ p =4 x Time φ (18-1): A = 0 . p = 2. (18-2): A = 0 . p = 3. (18-3): A = 0 . p = 4. π π φ p =5 x Time φ π π φ p =6 x Time φ (18-4): A = 0 . p = 5. (18-5): A = 0 . p = 6. Figure 18: The same as Fig.4 and Fig.12 except for H . These results are set as H = 10 − . . Concluding remark In this paper, we made the Hamiltonian formulation of the semilinear Klein-Gordon equation, and we derived the discrete equation with the structure-preserving scheme (SPS). To show the reliability of the simulations, we proposedthe constant value H C when the Hubble constant is zero case and ˜ H C when theHubble constant is nonzero case. With SPS, the Crank-Nicolson scheme (CNS)and the Runge-Kutta scheme (RKS), we performed some simulations in flatspacetime. Then, we showed the superiority of SPS to CNS and RKS in somesimulations. We performed some simulations with small ˜ H C and showed theinfluence of the Hubble constant on the numerical stability. Especially, if thesignature of the nonlinear term is negative, the simulations stop in some cases.However, with the negative nonlinear term, we showed the enough large valueof the Hubble constant gives the long and stable simulations. It is remarkablethat the diffusion effect caused by the positive Hubble constant is much strongerthan the nonlinear term. Thus, we conclude that we are able to perform stablesimulations when the Hubble constant is a sufficiently large. While the diffusioneffects are expected from the theoretical point of view since the equation (1.4)has the positive-dissipative term nH∂ t φ for the positive Hubble constant ( H >
H <
0) seems to be unstableand requires more delicate consideration since the numerical errors must berigorously estimated for the blow-up solutions in the unstable case, which willbe reported in the subsequent paper.
Acknowledgements
This work was supported by JSPS KAKENHI Grant Number 16H03940(M.N.).
References [1] G. Adomian, A review of the decomposition method in applied mathe-matics, J. Math. Anal. Appl. 135 (2) (1988) 501–544.[2] P. D’Ancona, A note on a theorem of J¨orgens, Math. Z. 218 (1995) 239–252.[3] P. D’Ancona, A. Di Giuseppe, Global existence with large data for anonlinear weakly hyperbolic equation, Math. Nachr. 231 (2001) 5–23.[4] A. Balogh, J. Banda, K. Yagdjian, High-performance implementation of aRunge-Kutta finite-difference scheme for the Higgs boson equation in thede Sitter spacetime, Commun. Nonlinear. Sci. Numer. Simul. 68 (2019)15–30.[5] K.C. Basak, P.C. Ray, R.K. Bera, Solution of non-linear Klein-Gordonequation with a quadratic non-linear term by Adomian decompositionmethod Commun. Nonlinear. Sci. Numer. Simul. 14 (2009) 718–723.236] D. Baskin, A Strichartz estimate for de Sitter space, The AMSI-ANUWorkshop on Spectral Theory and Harmonic Analysis, 97–104, Proc. Cen-tre Math. Appl. Austral. Nat. Univ., 44, Austral. Nat. Univ., Canberra,2010.[7] D. Baskin, Strichartz Estimates on Asymptotically de Sitter Spaces, An-nales Henri Poincar´e 14 (2) (2013) 221–252.[8] D. Furihata, Finite difference schemes for ∂u∂t = (cid:0) ∂∂x (cid:1) α δGδuα δGδu