On the numerical experiments of the Cauchy problem for semi-linear Klein-Gordon equations in the de Sitter spacetime

TL;DR

This paper presents numerical experiments analyzing the long-term behavior of solutions to semi-linear Klein-Gordon equations in de Sitter spacetime, highlighting the effects of the Hubble constant and the importance of structure-preserving schemes.

Contribution

It introduces a structure-preserving numerical scheme for simulating the equations and demonstrates its effectiveness in capturing long-term dynamics and diffusion effects.

Findings

Large Hubble constant induces strong diffusion effects.

The scheme preserves the numerically modified Hamiltonian.

Simulations show stable, long-term behavior of solutions.

Abstract

The computational analysis of the Cauchy problem for semi-linear Klein-Gordon equations in the de Sitter spacetime is considered. Several simulations are performed to show the time-global behaviors of the solutions of the equations in the spacetime based on the structure-preserving scheme. It is remarked that the sufficiently large Hubble constant yields the strong diffusion-effect which gives the long and stable simulations for the defocusing semi-linear terms. The reliability of the simulations is confirmed by the preservation of the numerically modified Hamiltonian of the equations.

Figures

Equations

$${G^{\alpha }}_{\beta }+\Lambda {g^{\alpha }}_{\beta }=\mu {T^{\alpha }}_{\beta },$$

$${T^{\alpha }}_{\beta }:=\text{mbox}diag(-\rho c^2,p,\dots ,p).$$

$$-\Lambda {g^{\alpha }}_{\beta }=\mu \text{mbox}diag(-\tilde{\rho }{c}^2,\tilde{p},\dots ,\tilde{p}),\tilde{\rho }:=\frac{\Lambda }{\mu c^2},\tilde{p}:=-\frac{\Lambda }{\mu },$$

$$d{s}^2=g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=-c^{2}d{t}^2+e^{2Ht}{\delta }_{ij}d{x}^{i}d{x}^j,$$

$$\mathcal{L}:=-\frac{\sqrt{-g}}{2}\left\{g^{\alpha \beta }({\nabla }_{\alpha }\varphi )({\nabla }_{\beta }\varphi )+\frac{c^2{m}^2}{{\hslash }^2}{\varphi }^2+2V(\varphi )\right\},$$

$$-\frac{1}{\sqrt{-g}}{\partial }_{\alpha }(\sqrt{-g}{g}^{\alpha \beta }{\partial }_{\beta }\varphi )+\frac{m^2{c}^2}{{\hslash }^2}\varphi +V^{\prime }(\varphi )=0.$$

$${\partial }_t^2\varphi +nH{\partial }_t\varphi -c^2{e}^{-2Ht}{\delta }^{ij}{\partial }_i{\partial }_j\varphi +\frac{m^2{c}^4}{{\hslash }^2}\varphi +c^2{V}^{\prime }(\varphi )=0,$$

$${\partial }_t{e}^0+{\partial }_j{e}^j+e_{\ast }=0,$$

$$e^0$$

$$e^j$$

$$e_{\ast }$$

$$E(t):={\int }_{{\mathbb{R}}^n}{e}^0(t,x)dx,$$

$$E(t)+{\int }_0^t{\int }_{{\mathbb{R}}^n}{e}_{\ast }(s,x)dxds=E(0)$$

$$V(\varphi ):=\frac{\lambda }{p+1}|\varphi {|}^{p+1},V^{\prime }(\varphi )=\lambda |\varphi {|}^{p-1}\varphi$$

$$\left\{\begin{array}[]{l}\left(\partial_{t}^{2}\phi+nH\partial_{t}\phi-c^{2}e^{-2Ht}\delta^{ij}\partial_{i}\partial_{j}\phi+m^{2}c^{4}\phi/\hbar^{2}+c^{2}\lambda|\phi|^{p-1}\phi\right)(t,x)=0\\ \hfill\ \ \ \ \mbox{for}\ \ (t,x)\in[0,T)\times{\mathbb{R}}^{n}\\ \phi(0,\cdot)=\phi_{0}(\cdot)\in H^{1}({\mathbb{R}}^{n}),\ \ \partial_{t}\phi(0,\cdot)=\phi_{1}(\cdot)\in L^{2}({\mathbb{R}}^{n}),\end{array}\right.$$

$$1\leq p\left\{\begin{array}[]{ll}<\infty&\mbox{if}\ \ n=1,2\\ \leq 1+\frac{2}{n-2}&\mbox{if}\ \ n\geq 3.\end{array}\right.$$

$$\begin{gathered} \left\{{\partial }_t^2+(nH-2\kappa ){\partial }_t-c^2{e}^{-2Ht}{\delta }^{ij}{\partial }_i{\partial }_j+\frac{m^2{c}^4}{{\hslash }^2}+\kappa (\kappa -nH)\right\}u \\ +c^2{e}^{\kappa t}{V}^{\prime }(e^{-\kappa t}u)=0, \end{gathered}$$

$$({\partial }_t^2-c^2{e}^{-2Ht}{\delta }^{ij}{\partial }_i{\partial }_j+M^2)u(t,x)+c^2{e}^{nHt/2}{V}^{\prime }(e^{-nHt/2}u(t,x))=0$$

$$\mathcal{L}$$

$$\psi :=\frac{\partial \mathcal{L}}{\partial ({\partial }_t\varphi )}=\frac{e^{nHt}}{c}{\partial }_t\varphi ,$$

$$\mathcal{H}$$

$$H_C(t):={\int }_{{\mathbb{R}}^n}\mathcal{H}(t,x)dx$$

$$\displaystyle\begin{array}[]{lll}\mbox{the kinetic term}&:&K(t):=\displaystyle\int_{{\mathbb{R}}^{n}}\frac{1}{2}(\partial_{t}\phi(t,x))^{2}dx,\\ \mbox{the diffusion term}&:&D(t):=\displaystyle\int_{{\mathbb{R}}^{n}}\frac{1}{2}c^{2}e^{-2Ht}\delta^{ij}(\partial_{i}\phi(t,x))(\partial_{j}\phi(t,x))dx,\\ \mbox{the mass term}&:&M(t):=\displaystyle\int_{{\mathbb{R}}^{n}}\frac{m^{2}c^{4}}{2\hbar^{2}}\phi(t,x)^{2}dx,\\ \mbox{the nonlinear term}&:&N(t):=\displaystyle\int_{{\mathbb{R}}^{n}}c^{2}V(\phi(t,x))dx=\displaystyle\int_{{\mathbb{R}}^{n}}\frac{\lambda c^{2}}{p+1}|\phi|^{p+1}dx\end{array}$$

$$E(t)=K(t)+D(t)+M(t)+N(t).$$

$$\mathcal{H}=\frac{e^{nHt}}{c}{e}_0,$$

$$H_C(t)=\frac{e^{nHt}}{c}(K(t)+D(t)+M(t)+N(t)).$$

$$\displaystyle\begin{array}[]{ll}\partial_{t}\phi&:=\displaystyle{\frac{\delta\mathcal{H}}{\delta\psi}=ce^{-nHt}\psi},\\ \partial_{t}\psi&:=\displaystyle{-\frac{\delta\mathcal{H}}{\delta\phi}=ce^{(n-2)Ht}\delta^{ij}(\partial_{i}\partial_{j}\phi)-\frac{c^{3}m^{2}}{\hbar^{2}}e^{nHt}\phi-ce^{nHt}V^{\prime}(\phi)}.\end{array}$$

$${\partial }_t\mathcal{H}=nH\mathcal{H}+c^2{e}^{-2Ht}{\partial }_i\{\psi {\delta }^{ij}({\partial }_j\varphi )\}-{\mathcal{H}}_{\ast },$$

$${\mathcal{H}}_{\ast }=H\{cn{e}^{-nHt}{\psi }^2+c{e}^{(n-2)Ht}{\delta }^{ij}({\partial }_i\varphi )({\partial }_j\varphi )\},$$

$${\partial }_t{H}_C(t)$$

Full text

††journal: Journal of Computational and Applied Mathematics

1 Introduction

The mathematical structure of partial differential equations in non-flat spacetimes is changed by the variations of curvatures since the differential operators are influenced by the metrics of the spacetimes. The de Sitter spacetime is one of the solutions of the Einstein gravitational equations with the cosmological constant in the vacuum. The space is expanding or contracting along the time in this spacetime. To study the effects of the spatial variation, we consider the semilinear field equation of the Klein-Gordon type, and we carry out the numerical experiments of the solutions and their energy. To investigate the partial differential equations with computational analysis, we need to perform high precise and accurate simulations. There are many methods to make discretized equations.

The Adomian decomposition method [1] can provide analytical approximation solutions of linear and nonlinear differential equations. Since this method dose not require the linearization, it would be appropriate to make dicretization of the nonlinear differential equations. It is reported that the simulations of the nonlinear Klein-Gordon equation in the flat spacetime with the method [5]. In [31], the differential transform method [32] and the variational iteration method [24] are used to perform the nonlinear Klein-Gordon equation in the de Sitter spacetime. The former is one of the power series expansions. The virtue is easy to estimate the numerical errors because the truncation higher order terms are corresponding to the dominant numerical errors. Since the latter is classified as the method of Lagrange multipliers, it would be close to the exact solutions by repeating the iteration if the appropriate conditions are given. In addition, the results with the explicit fourth-order Runge-Kutta method are reported in [4]. The Runge-Kutta method is well-known discretized scheme for the nonlinear differential equations.

The above methods are generic discretization schemes for the differential equations to calculate simulations. Therefore, the properties of the differential equations in continuous level often lost in the discretization, and it makes the numerical errors in the evolutions. On the other hand, we use the structure-preserving scheme (SPS) called as “discrete variational derivative method” in this paper. The concept of this method is to conserve the structures and the properties of the equations in continuous level. This method is widely used to derive the discrete equations of partial differential equations (see e.g., [8, 9, 22]). To make discretized equations using SPS, we need the Lagrangian or the Hamiltonian[12, 25]. This is because the derivations of the discretized equations with SPS are similar with the process of making the continuous equations using the variational principle. In comparison with the construction of the Lagrangian formulation, the one of the Hamiltonian formulation is clear. Especially, since the Hamiltonian means the total energy of the system, this value is regarded as one of the constraints of the system. By monitoring the value in the evolution, we can judge whether the numerical simulations are success or not. Conservation of the constraints is one of the necessary conditions to perform successive simulations. If the constraints dose not conserve in the evolution, the simulations fail (e.g. [21]). Therefore, we use the Hamiltonian formulation in this paper.

We start from the introduction of the de Sitter spacetime. Let $n\geq 1$ be the spatial dimension, $c>0$ be the speed of light. In the following, Greek letters $\alpha,\beta,\gamma,\dots$ run from [math] to $n$, Latin letters $j,k,\ell,\dots$ run from $1$ to $n$. Let $g_{\alpha\beta}dx^{\alpha}dx^{\beta}$ be the metric in $\mathbb{R}^{1+n}$. We denote by $(g_{\alpha\beta})_{0\leq\alpha,\beta\leq n}$ the $(1+n)\times(1+n)$-matrix whose $(\alpha,\beta)$-component is given by $g_{\alpha\beta}$. Put $g:={\rm det}(g_{\alpha\beta})_{0\leq\alpha,\beta\leq n}$, and let $(g^{\alpha\beta})_{0\leq\alpha,\beta\leq n}$ be the inverse matrix of $(g_{\alpha\beta})_{0\leq\alpha,\beta\leq n}$. We use the Einstein rule for the sum of indices of the tensors, for example, ${T^{\alpha}}_{\alpha}:=\sum_{\alpha=0}^{n}{T^{\alpha}}_{\alpha}$ and ${T^{i}}_{i}:=\sum_{i=1}^{n}{T^{i}}_{i}$. The change of upper and lower indices of the tensors is done by $g_{\alpha\beta}$ and $g^{\alpha\beta}$, for example, ${T^{\alpha}}_{\beta}:=g^{\alpha\gamma}T_{\gamma\beta}$.

For a stress-energy tensor ${T^{\alpha}}_{\beta}$, the $(1+n)$-dimensional Einstein equation is defined by

[TABLE]

where ${G^{\alpha}}_{\beta}$ is the Einstein tensor, $\Lambda$ is the cosmological constant, and $\mu$ is the Einstein gravitational constant. When we consider the universe filled with the perfect fluid of the mass density $\rho$ and the pressure $p$, we are able to use the stress-energy tensor ${T^{\alpha}}_{\beta}$ of the perfect fluid given by

[TABLE]

When the cosmological term $\Lambda{g^{\alpha}}_{\beta}$ is transposed to right-hand-side in (1.1), the term is regarded as a part of the stress-energy tensor. Then the term $-\Lambda{g^{\alpha}}_{\beta}$ is rewritten as

[TABLE]

and the cosmological constant $\Lambda>0$ is regarded as the energy which has positive density and negative pressure. Thus, we regard the cosmological constant $\Lambda$ as “the dark energy.” The study of roles of the cosmological constant and the spatial variance is important to describe the history of the universe, especially, the inflation and the accelerating expansion of the universe (see e.g., [13], [14], [17], [18], [19], [20]). One of the solutions of the equation (1.1) is the de Sitter spacetime, and its line-element with the spatial zero-curvature is given by

[TABLE]

where $t:=x^{0}$, $H$ is the Hubble constant, and $\delta_{ij}$ denotes the Kronecker delta. The relation between the cosmological constant $\Lambda$ and the Hubble constant $H$ is given by $H=\sqrt{2\Lambda/(n^{2}-n)}$. We can write $(g_{\alpha\beta})_{0\leq\alpha,\beta\leq n}={\rm diag}(-c^{2},e^{2Ht},\cdots,e^{2Ht})$.

Next, we introduce the semilinear Klein-Gordon equation (KGE) in curved spacetimes. KGE presents the equation of motion for the massive scalar field. Then, to derive the equation of motion described by a real-valued function $\phi$ with the mass $m$ and the potential $V$, we consider the Lagrangian density $\mathcal{L}$ defined by

[TABLE]

where $\nabla_{\mu}$ is the covariant derivative associated with $g_{\mu\nu}$ and $\hbar$ is the Planck constant. If the spacetime is flat, the Lagrangian density $\mathcal{L}$ is consistent with the well-known Lagrangian density of KGE (see e.g., [25]). The Euler-Lagrange equation of the action $\int_{{\mathbb{R}}^{1+n}}\mathcal{L}\,dtdx$ gives KGE in curved spacetime

[TABLE]

In the de Sitter spacetime, KGE is rewritten as

[TABLE]

where $V^{\prime}$ denotes the derivative $V^{\prime}(r)=dV(r)/dr$ for $r\in{\mathbb{R}}$ (see, e.g., [7, 15, 26]). Multiplying $\partial_{t}\phi$ to the both sides of (1.4), we have the divergence formula

[TABLE]

where $e^{\alpha}$, $0\leq\alpha\leq n$, and $e_{\ast}$ mean the energy density and the external energy, respectively. They are defined by

[TABLE]

We define the total energy by

[TABLE]

then we have the energy estimate

[TABLE]

for $t\in[0,\infty)$. Let us consider the defocusing potential ($V\geq 0$), by which $E$ is a positive-valued function for any non-trivial solution $\phi$. When the space is expanding ($H>0$), the estimate (1.6) shows the dissipative property $E(t)<E(0)$ for $t>0$. On the other hand, when the space is contracting ($H<0$), (1.6) shows the anti-dissipative property $E(t)>E(0)$ for $t>0$.

We introduce the known theoretical results on KGE (1.4). Let us consider the scalar potential $V$ of power type and its derivative $V^{\prime}$ given by

[TABLE]

for $\lambda\in{\mathbb{R}}$ and $1<p<\infty$. We consider the Cauchy problem

[TABLE]

where $\phi_{0}$, $\phi_{1}$ are given initial data, $H^{1}(\mathbb{R}^{n})$ denotes the Sobolev space and $L^{2}(\mathbb{R}^{n})$ denotes the Lebesgue space. On the Cauchy problem (1.8) for the positive Hubble constant $H>0$, the following theoretical results have been shown. The first result shows the existence of the time-local solutions for any initial data. The second result shows that the local solutions are time-global solutions if the initial data are sufficiently small. The third result shows the local solutions are global for any data if the semilinear term is defocusing.

Theorem 1.1

([15, Theorems 1.2 and 1.7]) Let $H>0$. Assume $m>nH/2$. Let $p$ satisfy

[TABLE]

Then we have the following results.

(1) For any $\phi_{0}$ and $\phi_{1}$, there exists $T=T(\|\phi_{0}\|_{H^{1}({\mathbb{R}^{n}})}+\|\phi_{1}\|_{L^{2}({\mathbb{R}^{n}})})>0$ such that (1.8) has a unique solution $\phi$ in $C([0,T),H^{1}({\mathbb{R}^{n}}))\cap C^{1}([0,T),L^{2}({\mathbb{R}^{n}}))$.

(2) Let $n\leq 4$. If $\|\phi_{0}\|_{H^{1}({\mathbb{R}^{n}})}+\|\phi_{1}\|_{L^{2}({\mathbb{R}^{n}})}$ is sufficiently small, and $1+4/n\leq p$, then (1.8) has a unique solution $\phi$ in $C([0,\infty),H^{1}({\mathbb{R}^{n}}))\cap C^{1}([0,\infty),L^{2}({\mathbb{R}^{n}}))$.

(3) If $\lambda\geq 0$, then (1.8) has a unique global solution $\phi$ in $C([0,\infty),H^{1}({\mathbb{R}}^{n}))\cap C^{1}([0,\infty),L^{2}({\mathbb{R}}^{n}))$ for any data $\phi_{0}\in H^{1}(\mathbb{R}^{n})$ and $\phi_{1}\in L^{2}(\mathbb{R}^{n})$.

We refer to the corresponding known results on (1.8). D’Ancona and Giuseppe have shown in [2] and [3] global classical solutions for $(\partial_{t}^{2}-a(t)\delta^{ij}\partial_{i}\partial_{j})u+|u|^{p-1}u=0$ with some additional conditions on $a(t)\geq 0$ and $p$ when $n=1,2,3$. Yagdjian has shown in [28] small global solutions for (1.11) with arbitrary $n$ when the nonlinear term $f$ is of power type $p>1$, and the norm of initial data $\|u_{0}\|_{H^{s}({\mathbb{R}^{n}})}+\|u_{1}\|_{H^{s}({\mathbb{R}^{n}})}$ is sufficiently small for some $s>n/2\geq 1$ (see also [29] for the system of the equations). Galstian and Yagdjian has extended this result to the case of the Riemann metric space for each time slices in [11]. In [15], the energy solutions $(s=1)$ in Theorem 1.1 have been shown. This result 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^{\mu\nu}\nabla_{\mu}\nabla_{\nu}+\lambda)u+f(u)=0$ when $f(u)$ is a type of $|u|^{p-1}u$, $p=1+4/(n-1)$, $\lambda>n^{2}/4$, $(u_{0},u_{1})\in H^{1}\oplus L^{2}$, where $g^{\mu\nu}$ gives the asymptotic de Sitter spacetime (see also [6] for the cases $p=5$ with $n=3$, $p=3$ with $n=4$). Blow-up phenomena 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/2$ by the following reason. Since the first equation in (1.8) has the dissipative (or anti-dissipative) term $nH\partial_{t}\phi$, we use the transformation $u=e^{\kappa t}\phi$ for $\kappa\in\mathbb{R}$ to transform the equation to the Klein-Gordon equation. Then the equation is rewritten as

[TABLE]

and we obtain the equation

[TABLE]

when $\kappa=nH/2$, $M^{2}:=c^{4}m^{2}/\hbar^{2}-n^{2}H^{2}/4$. We assume $m$ is large such that $M^{2}>0$ (see [26]). Here, we expect some dissipative effects when $\kappa<nH/2$ in (1.10). However, $\kappa\leq nH/2$ is needed for the energy estimates for (1.11), and $\kappa\geq nH/2$ is needed for the contraction argument to show the existence of the solutions of (1.11). So that, we assume $\kappa=nH/2$, by which we lose the dissipative term $(nH-2\kappa)\partial_{t}u$ in (1.11). Since the first equation in (1.8) has the term $e^{-2Ht}\delta^{ij}\partial_{i}\partial_{j}\phi$ with the non-constant coefficient, it is not easy to obtain the critical theoretical results so far. Instead, in this paper, we carry out some numerical simulations, and show the detailed behaviors to clarify the dissipative effect of the spatial expansion.

2 Hamiltonian formulation of KGE

Let us consider the Hamiltonian density for the Lagrangian density $\mathcal{L}$. The Lagrangian density of KGE (1.3) in the de Sitter spacetime becomes

[TABLE]

Then, we define the momentum $\psi$ by

[TABLE]

the Hamiltonian density $\mathcal{H}$ by the Legendre transformation

[TABLE]

and the Hamiltonian $H_{C}$ by

[TABLE]

for $t\in{\mathbb{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 $(\partial_{t}\phi)^{2}/2$, $c^{2}e^{-2Ht}\delta^{ij}(\partial_{i}\phi)(\partial_{j}\phi)/2$, $m^{2}c^{4}\phi^{2}/(2\hbar^{2})$ and $c^{2}V(\phi)$. Namely,

[TABLE]

for $t\in{\mathbb{R}}$. With (2.7), the total energy (1.5) can be expressed as

[TABLE]

Since the relation between $\mathcal{H}$ and the energy density $e_{0}$ becomes

[TABLE]

$H_{C}$ can be expressed as

[TABLE]

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 the canonical equations

[TABLE]

Since $\mathcal{H}$ satisfies the differential equation

[TABLE]

where we have put

[TABLE]

$H_{C}$ satisfies the equation

[TABLE]

and by using (2.7), (2.13) can be expressed as

[TABLE]

The equation (2.13) (or (2.14)) indicates $H_{C}$ is not conserved in time evolution if the Hubble constant $H$ is not zero. Therefore, we define $\tilde{H}_{C}$ as

[TABLE]

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 spatial expansion and contraction by the Hubble constant $H$. Namely, we are able to understand the properties of the dark energy, which is equivalent to the cosmological constant mathematically, through the asymptotic behaviors of $\phi$, $E$ and $H_{C}$. However, it is not easy to obtain the detailed behaviors of them by theoretical methods since the waves of $\phi$ would be oscillating and interfered by the semilinear term $V(\phi)$, and (1.4) is the equation with the variable coefficient. Thus, we perform the numerical simulations on $\phi$, $E$ and $H_{C}$ to study the detailed dissipative effects of the spatial expansion.

3 Discrete KGE

To investigate KGE by computational analysis, we need the high precise and accurate numerical data. Then, we use SPS to get the high precise and accurate numerical results. With this scheme, the structures of the equations in the continuous case are preserved in the discrete case. In this paper, we use this scheme to make the discrete Hamiltonian formulation. The discrete values of the variable $u$ are defined as $u{}^{(\ell)}_{(k)}:=u(\ell\Delta t,k\Delta x)$, where $\ell$ is the time index, $\Delta t$ is the time mesh width, $k$ is the spatial grid index, and $\Delta x$ is the spatial mesh width. The details of this are in Ref.[9].

For comparison with SPS, we use the Crank-Nicolson scheme (CNS) and the fourth-order Runge-Kutta scheme (RKS) which are widely used for calculating partial differential equations numerically.

3.1 SPS

By using SPS, we can get a set of discrete evolution equations of KGE. We give a discrete Hamiltonian density $\mathcal{H}^{(\ell)}_{(k)}$ as

[TABLE]

where $t_{\ell}$ presents the time of the $\ell$th step and $\widehat{\delta}^{\langle 1\rangle}_{i}u^{(\ell)}_{(k)}$ is the centered space operator defined by

[TABLE]

For $\mathcal{H}^{(\ell)}_{(k)}$, $\phi^{(\ell)}_{(k)}$ and $\psi^{(\ell)}_{(k)}$, we define the values

[TABLE]

by the equation

[TABLE]

Then, a set of discrete KGE with SPS becomes

[TABLE]

The discrete evolution equations (3.9) are corresponding to (2.12) for the continuous case. The nonlinear term of the last term in (3.9) is expressed as

[TABLE]

when $p$ is a natural number. We refer to the expression in the nonlinear Schrödinger equations in Ref.[9] to make above relation. With (3.9), the time difference of $\mathcal{H}^{(\ell)}_{(k)}$ is calculated as

[TABLE]

[TABLE]

[TABLE]

where we use the relation of $t_{\ell+1}=t_{\ell}+\Delta t$. We define the discrete Hamiltonian $H_{C}{}^{(\ell)}$ by

[TABLE]

then (3.11) is calculated as

[TABLE]

where $K^{(\ell)}$, $D^{(\ell)}$, $M^{(\ell)}$ and $N^{(\ell)}$ are defined by

[TABLE]

$K^{(\ell)}$, $D^{(\ell)}$, $M^{(\ell)}$ and $N^{(\ell)}$ are discrete values corresponding to $K(t)$, $D(t)$, $M(t)$ and $N(t)$ defined by (2.7), respectively. We set

[TABLE]

then $\tilde{H}_{C}{}^{(n)}$ is a discretized value of $\tilde{H}_{C}$ defined by (2.15). This value is a guideline to perform precise numerical simulations in the cases of $H\neq 0$.

3.2 CNS and RKS

CNS and RKS are frequently used to perform partial differential equations. The discrete equations of (2.12) with CNS are given by

[TABLE]

where the right hand sides of (3.22) are defined as

[TABLE]

and the second-order central difference operator $\widehat{\delta}^{\langle 2\rangle}_{ij}$ is defined as

[TABLE]

The discrete equations with RKS are give by

[TABLE]

where

[TABLE]

4 Numerical tests

In this section, we carry out the numerical simulations on $\phi$, $E$ and $H_{C}$ (or $\tilde{H}_{C}$). Since the physical background of the Einstein equation is for the $1+3$ dimensional spacetime, we consider the case $n=3$ in the following.

The numerical settings are below. Set $(x^{1},x^{2},x^{3})=(x,y,z)\in\mathbb{R}^{3}$.

Initial condition:

[TABLE]

• 2.

Numerical domains: $0\leq x\leq 1$, $0\leq t\leq 1000$.

• 3.

Boundary condition: periodic.

• 4.

Grids: $\Delta x=1/200$, and $\Delta t=1/500$.

• 5.

Physical settings: mass $m=1$, speed of light $c=1$, Planck constant $\hbar=1$.

• 6.

The number of exponent in the nonlinear term: $p=2,3,4,5,6$.

We consider the solution $\phi(t,x,y,z)$ for $(t,x,y,z)\in[0,\infty)\times\mathbb{R}^{3}$ which is a constant along the variables $y$ and $z$ to give simple and reliable simulations. More general solutions and data will be considered in the future. We perform the simulations by changing the values of the amplitude $A$ in the initial condition, the coefficient value $\lambda$ of the nonlinear term, and the Hubble constant $H$. First, we perform some simulations to study the accuracies with SPS. For comparison with SPS, we use the Crank-Nicolson scheme (CNS) and the Runge-Kutta scheme (RKS) which are widely used for calculating partial differential equations numerically. Because of the nonlinearity of KGE (2.12), the discretized equations are calculated implicitly. Thus, the numerical simulations are performed iteratively. For SPS and CNS, the iterations are nine times and three times, respectively. If we calculate more times iteratively for each scheme, the numerical results become worse in our codes. For the wave equation

[TABLE]

with (1.7), namely, $H=0$, $m=0$ in (1.4) with (1.7), the scaling argument for the solution $\phi_{\rho}(t,x):=\rho^{2/(p-1)}\phi(\rho t,\rho x)$ for $\rho\in\mathbb{R}$ yields the critical exponent of $p(s):=1+4/({3}-2s)$ for the well-posedness of the Cauchy problem of the equation in the Sobolev space $H^{s}(\mathbb{R}^{{3}})$ for $s\in\mathbb{R}$. Especially, $p(-{3/2})$ is called the Fujita exponent, $p(0)$ is called the conformal exponent, $p(1)$ is called the energy critical exponent. This exponent $p(s)$ also plays crucial roles for the Klein-Gordon equation, and it satisfies

[TABLE]

By the simulations for $p=2,3,4,5,6$, we are able to know the behaviors of the solutions for $p$ which is equivalent or close to these critical exponents.

4.1 Case 1: $H=0$

To compare with the numerical results of SPS, CNS and RKS, we set the Hubble constant $H$ as zero in this subsection. In this case, $H_{C}$ is consistent with $E$, and $H_{C}$ is constant with time evolutions theoretically. To show the accuracy of the simulations, we monitor $H_{C}$.

We draw $H_{C}$, and the differences between $H_{C}$ and the initial value $H_{C}(0)$ in the evolutions with SPS, CNS and RKS for $A=4$ and $\lambda=1$ in Fig.1. From the bottom panels (1-b1)–(1-b3), we see the values of $|(H_{C}-H_{C}(0))/H_{C}|$ with CNS and RKS are larger than the ones with SPS.

Next, we show the results for $A=0.9$ and $\lambda=-1$ in Fig.2.

We see the simulations stop before $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. In comparison with Fig.2, the simulations in Fig.1 are robust. Since it is necessary condition 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 numerical scheme as SPS hereafter.

The results of $\phi$ with SPS are shown in Fig.3 and Fig.4. Fig.3 shows the case of $A=4$ and $\lambda=1$. We see that the vibrations occur for $p=5,6$. Fig.4 shows the case of $A=0.9$ and $\lambda=-1$. We see all of the simulations stop before $t=140$.

4.2 Case 2: $H=10^{-3}$

In this subsection and next subsection, we perform the numerical simulations for $H>0$. First, we calculate for $\lambda=1$ and $H=10^{-3}$. When $H\neq 0$, $H_{C}$ is not constant in time evolutions. Since we cannot judge the reliability of simulations by monitoring the values of $H_{C}$, we show $\tilde{H}_{C}$ defined by (2.15) instead of $H_{C}$.

Fig.5 shows $\tilde{H}_{C}$, and the differences between $\tilde{H}_{C}$ and the initial value $\tilde{H}_{C}(0)$ for $A=1,3,4$. We see the value of $p=2$ is the largest and that of $p=6$ is the smallest in the panel (5-t1). Contrarily, in the others of the top panels we 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 other hand, if $A>1$, the value becomes large as $p$ becomes large. From the bottom panels (5-b1)–(5-b3), we see $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ of the case $p=6$ is the largest of them in each panel but all of the values in the bottom panels are enough small.

Fig.6 shows the total energy $E$ defined by (1.5). We see that the energy is exponential decay. Furthermore, the slopes of the lines in the panels of the figure are almost the same. This result means the diffusion effect with time evolutions 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 between the values of $E$ for $A=3,4$. To investigate the reason, we show the components of $E$ defined by (2.7).

The kinetic term $K$, the diffusion term $D$, the mass term $M$ and the nonlinear term $N$ for $A=1$ are shown in Fig.7 and the ones for $A=4$ are in Fig.8. We see that $K$ and $D$ are the dominant terms, and $M$ and $N$ are very small compared with $K$ and $D$ in Fig.7. The decline of the energy-component $N$ is becoming 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 is because there are the differences between $E$s in Fig.6.

Fig.9 shows $\phi$ for $A=4$ and $\lambda=1$, and the numerical settings are the same as Fig.3 except for the value of $H$. The vibrations occur in the only case of $p=6$. Since the vibrations occur in the cases of $p=5,6$ in Fig.3, the Hubble constant seems to affect generations of the vibrations of $\phi$.

Next we show the results of the simulations with the case of $\lambda=-1$.

Fig.10 shows $\tilde{H}_{C}$ and $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ for $A=0.3,0.6,0.9$, $\lambda=-1$ and $H=10^{-3}$. While we see the simulation for $A=0.9$ and $p=2$ stops before $t=200$ in Panel (10-t3) and Panel (10-b3), the values $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ for the other cases are enough small.

Fig.11 shows $E$ for $A=0.3,0.6,0.9$, $\lambda=-1$ and $H=10^{-3}$. It seems there is no difference between the slopes of all the lines in Fig.6 and Fig.11 regardless of the values of $A$, $p$ and $\lambda$. Thus, the dissipation effect of $E$ is not caused by changing in $A$, $p$ and the signature of $\lambda$.

Fig.12 shows $\phi$ for $A=0.9$, $\lambda=-1$ and $H=10^{-3}$. These results are the same as Fig.4 except for the value of $H$. By comparing Fig.4 with Fig.12, the simulation times for $H=10^{-3}$ are longer than that for $H=0$. It means the Hubble constant would influence the stability of the simulation. By comparing Fig.9 and Fig.12, the simulations for $\lambda=1$ are robust against the case of $\lambda=-1$ since the calculation time for $p=2$ and $\lambda=-1$ is short. The tendency is the same as the case of $H=0$.

4.3 Case 3: $H=10^{-2}$

In this subsection, we perform some simulations with $H=10^{-2}$ and $\lambda=1$.

Fig.13 shows $\tilde{H}_{C}$ and $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ for $H=10^{-2}$ and this figure is the same as Fig.5 except for the value of $H$. We see all the values of $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ for $H=10^{-2}$ are enough small.

Fig.14 shows $E$ with the Hubble constant as $10^{-2}$. There are few differences between the changes of $p$ compared with the cases of $H=10^{-3}$ in Fig.6. By comparing Fig.6 and Fig.14, the diffusion effect of $E$ for $H=10^{-2}$ is stronger than that for $H=10^{-3}$. Thus, the diffusion effect is mainly caused by the Hubble constant, and the effect would become strong as $H$ becomes large.

In Fig.15 which shows $\phi$ for $H=10^{-2}$, we see the waveform is damped and there are few vibrations of the waveform. By comparing Fig.3, Fig.9 and Fig.15, the stability of the simulations becomes good as $H$ becomes large since the vibrations of $\phi$ decreases as $H$ becomes large.

Next we show the results for $\lambda=-1$.

In Fig.16, we see the simulations are reliable since all of the values of $|(\tilde{H}_{C}-\tilde{H}_{C}(0))/\tilde{H}_{C}|$ are enough small.

Fig.17 shows the total energy $E$. By comparing Fig.11 and Fig.17, we see the diffusion effect for $\lambda=-1$ and $H=10^{-2}$ is stronger than that for $\lambda=-1$ and $H=10^{-3}$. Thus, considering the results for $\lambda=1$, the diffusion effect is mainly caused by not the signature of the nonlinear term but the Hubble constant.

Fig.18 shows $\phi$ for $\lambda=-1$ and $H=10^{-2}$. The simulation times become long as $H$ becomes large from the results in Fig.4, Fig.12 and Fig.18. These results indicate the simulations become stable as $H$ becomes large.

5 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 proposed the constant value $H_{C}$ when the Hubble constant is zero case and $\tilde{H}_{C}$ when the Hubble constant is nonzero case. With SPS, the Crank-Nicolson scheme (CNS) and the Runge-Kutta scheme (RKS), we performed some simulations in flat spacetime. Then, we showed the superiority of SPS to CNS and RKS in some simulations. We performed some simulations with small $\tilde{H}_{C}$ and showed the influence of the Hubble constant on the numerical stability. Especially, if the signature of the nonlinear term is negative, the simulations stop in some cases. However, with the negative nonlinear term, we showed the enough large value of the Hubble constant gives the long and stable simulations. It is remarkable that the diffusion effect caused by the positive Hubble constant is much stronger than the nonlinear term. Thus, we conclude that we are able to perform stable simulations when the Hubble constant is a sufficiently large. While the diffusion effects are expected from the theoretical point of view since the equation (1.4) has the positive-dissipative term $nH\partial_{t}\phi$ for the positive Hubble constant ($H>0$), the case of the negative Hubble constant ($H<0$) seems to be unstable and requires more delicate consideration since the numerical errors must be rigorously estimated for the blow-up solutions in the unstable case, which will be reported in the subsequent paper.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Number 16H03940 (M.N.).