text:
# 1 Introduction

## 1 Introduction

We begin by recalling the variational principle for Lagrange mechanics. Suppose $`Q`$ denotes the configuration space with coordinates $`q^i`$ and $`TQ`$ the tangent bundle with coordinates $`(q^i,\dot{q}^i)`$, $`i=1,2,\mathrm{},n`$. In the following, we will work with the extended configuration space $`R\times Q`$ with coordinates $`(t,q^i)`$ and the extended tangent bundle $`R\times TQ`$ with coordinates $`(t,q^i,\dot{q}^i)`$. Here $`t`$ denotes the time \[References\].

Consider a Lagrangian $`L:R\times TQR`$. The corresponding action functional is defined by

$`S((t,q^i(t)))={\displaystyle _a^b}L(t,q^i(t),\dot{q}^i(t))𝑑t,`$ (1.1)
where $`q^i(t)𝒞^2([a,b])`$, which is the set of all $`C^2`$ curves in $`Q`$ defined on $`[a,b]`$.

The variational principle seeks the curves $`q^i(t)`$ in $`𝒞^2([a,b])`$, for which the action functional $`S`$ is stationary under variations of $`q^i(t)`$ with fixed endpoints.

Consider a vector field on $`R\times Q`$

$`V=\xi (t,q^j){\displaystyle \frac{}{t}}+{\displaystyle \underset{i=1}{\overset{n}{}}}\varphi ^i(t,q^j){\displaystyle \frac{}{q^i}}.`$ (1.2)
Let $`F^ϵ`$ be the flow of $`V`$. For $`(t,q^j)R\times Q`$, we have $`F^ϵ(t,q^j)=(\stackrel{~}{t},\stackrel{~}{q}^i)`$.

$`\stackrel{~}{t}=f(ϵ,t,q^j),`$ (1.3)
$`\stackrel{~}{q}^i=g^i(ϵ,t,q^j),`$ (1.4)
where

$`\begin{array}{cc}& {\displaystyle \frac{d}{dϵ}}|_{ϵ=0}f(ϵ,t,q^j)=\xi (t,q^j),\hfill \\ & {\displaystyle \frac{d}{dϵ}}|_{ϵ=0}g^i(ϵ,t,q^j)=\varphi ^i(t,q^j).\hfill \end{array}`$ (1.5)
The transformation (1.3-1.4) transforms a curve $`q^j(t)`$ in $`Q`$ into another curve $`\stackrel{~}{q}^i(ϵ,\stackrel{~}{t})`$ in $`Q`$ determined by

$`\stackrel{~}{t}=f(ϵ,t,q^j(t)),`$
$`\stackrel{~}{q}^i=g^i(ϵ,t,q^j(t)).`$
Before calculating the variation of $`S`$, we consider the first order prolongation of $`V`$

$`\text{pr}^1V=\xi (t,q^j){\displaystyle \frac{}{t}}+{\displaystyle \underset{i=1}{\overset{n}{}}}\varphi ^i(t,q^j){\displaystyle \frac{}{q^i}}+{\displaystyle \underset{i=1}{\overset{n}{}}}\alpha ^i(t,q^j,\dot{q}^j){\displaystyle \frac{}{\dot{q}^i}},`$ (1.6)
where $`\text{pr}^1V`$ denote the first order prolongation of $`V`$ and

$`\alpha ^i(t,q^j,\dot{q}^j)=D_t\varphi ^i(t,q^j)\dot{q}^iD_t\xi (t,q^j),`$ (1.7)
where $`D_t`$ denotes the total derivative \[References\].

Now we calculate the variation of $`S`$ directly

$`\delta S=`$ $`{\displaystyle \frac{d}{dϵ}}|_{ϵ=0}S((\stackrel{~}{t},\stackrel{~}{q}^i(ϵ,\stackrel{~}{t})))`$
$`=`$ $`{\displaystyle \frac{d}{dϵ}}|_{ϵ=0}{\displaystyle _{\stackrel{~}{a}}^{\stackrel{~}{b}}}L(\stackrel{~}{t},\stackrel{~}{q}^i(ϵ,\stackrel{~}{t}),{\displaystyle \frac{d}{d\stackrel{~}{t}}}\stackrel{~}{q}^i(ϵ,\stackrel{~}{t}))𝑑\stackrel{~}{t}`$
$`=`$ $`{\displaystyle _a^b}\left[\left({\displaystyle \frac{L}{t}}+{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{L}{\dot{q}^i}}\dot{q}^iL\right)\right)\xi +\left({\displaystyle \frac{L}{q^i}}{\displaystyle \frac{d}{dt}}{\displaystyle \frac{L}{\dot{q}^i}}\right)\varphi ^i\right]𝑑t`$
$`+\left[\left(L{\displaystyle \frac{L}{\dot{q}^i}}\dot{q}^i\right)\xi +{\displaystyle \frac{L}{\dot{q}^i}}\varphi ^i\right]|_a^b.`$ (1.8)
If $`\xi (a,q^j(a))=\xi (b,q^j(b))=0`$ and $`\varphi ^i(a,q^j(a))=\varphi ^i(b,q^j(b))=0`$, the requirement of $`\delta S=0`$ yields the energy evolution equation

$`{\displaystyle \frac{L}{t}}+{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{L}{\dot{q}^i}}\dot{q}^iL\right)=0`$ (1.9)
from $`\xi `$ and the Euler-Lagrange equation

$`{\displaystyle \frac{L}{q^i}}{\displaystyle \frac{d}{dt}}{\displaystyle \frac{L}{\dot{q}^i}}=0`$ (1.10)
from $`\varphi ^i`$.111Note that Hamilton’s variational principle holds in fact for the vertical variation, while the horizontal variation gives rise to the conservation relation between the Euler-Lagrange equation and energy conservation, see, for example, \[References\]. Here, in order to transfer to discrete case more directly, this variational requirement is employed.

If $`L`$ does not depend on $`t`$ explicitly, i.e., $`L`$ is conservative and $`\frac{L}{t}=0`$, then (1.9) becomes the energy conservation law

$`{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{L}{\dot{q}^i}}\dot{q}^iL\right)=0.`$ (1.11)
If we drop the requirement

$`\xi (a,q^j(a))=\xi (b,q^j(b))=0,`$
$`\varphi ^i(a,q^j(a))=\varphi ^i(b,q^j(b))=0,`$
we can define from the second term in (1.8) the extended Lagrangian one form

$`\theta _L=\left(L{\displaystyle \frac{L}{\dot{q}^i}}\dot{q}^i\right)dt+{\displaystyle \frac{L}{\dot{q}^i}}dq^i.`$ (1.12)
We define the extended Lagrangian two form to be

$$\omega _L=d\theta _L.$$

Suppose $`g^i(t,v_𝒒)`$ is a solution of (1.10) depending on the initial condition $`v_𝒒TQ`$. Restricting $`\stackrel{~}{q}^i(ϵ,\stackrel{~}{t})`$ to the solution space of (1.10) and using the same method in \[References\], we can prove

$`(\text{pr}^1g)^{}\omega _L=\omega _L,`$ (1.13)
where $`\text{pr}^1g^i(s,v_𝒒))=(s,g^i(s,v_𝒒),\frac{d}{ds}g^i(s,v_𝒒))`$ denotes the first order prolongation of $`g^i(s,v_𝒒)`$.

By defining a discrete Lagrangian with fixed steps and constructing a discrete action, Veselov develops a kind of discrete mechanics \[References, References\]. See also \[References, References\]. The resulting discrete Euler-Lagrange equation (variational integrator) preserves a discrete Lagrange two form. However the discrete versions of (1.9) and (1.11) are not reflected in Veselov discrete mechanics. Using variable steps and a defined discrete energy, Kane, Marsden and Ortiz obtain symplectic-energy integrators \[References\]. Chen, Guo and Wu develop a discrete total variation and prove the symplectic-energy integrators preserve a discrete extended Lagrange two form \[References\]. Guo and Wu also present the difference variational approach with variable step-lengths \[References\].

The connection between Veselov discrete mechanics and finite element methods was first suggested in \[References\]. Symplectic and multisymplectic structures in simple finite element methods are explored in \[References\]. As is known, it is difficult to obtain high order discretizations in Veselov discrete mechanics. We find that the finite element methods provide a natural way to obtain a kind of discrete Lagrangian, including those of high order. The purpose of this paper is to develop such kind of discrete mechanics based on finite element methods.

This paper is organized as follows. In the next section, we will present the discrete mechanics based on finite element methods. Section 3 is devoted to discussing the role of time steps. We finish this paper by drawing some conclusions in Section 4.

## 2 Discrete mechanics based on finite element methods

### 2.1 Generalized Veselov discrete mechanics

We first present Veselov discrete mechanics in a more general formulation that allows for variable time steps. In this generalized Veselov discrete mechanics, $`M\times M`$ is used for the discrete version of $`R\times TQ`$. Here $`M=R\times Q`$. A point $`(t_k,q_k;t_{k+1},q_{k+1})M\times M`$ corresponds to a tangent vector $`\frac{q_{k+1}q_k}{t_{k+1}t_k}`$. Here for brevity of notations, we restrict ourselves to one-dimensional configuration space $`Q`$. A discrete Lagrangian is defined to be $`𝕃:M\times MR`$ and the corresponding action to be

$`𝕊={\displaystyle \underset{k=0}{\overset{N1}{}}}𝕃(t_k,q_k,t_{k+1},q_{k+1})(t_{k+1}t_k).`$ (2.1)
The discrete variational principle is to extremize $`𝕊`$ for variations of $`q_k`$ holding the endpoints $`q_0`$ and $`q_N`$ fixed. In this discrete variational principle, $`t_k`$ play the role of parameters (They can also play the role of variables, see next section). This discrete variational principle determines a discrete flow $`\mathrm{\Phi }:M\times MM\times M`$ by

$`\mathrm{\Phi }(t_{k1},q_{k1},t_k,q_k)=(t_k,q_k,t_{k+1},q_{k+1}).`$ (2.2)
Here $`q_{k+1}`$ is found from the following discrete Euler-Lagrange equation (variational integrator)

$`(t_{k+1}t_k)D_2𝕃(t_k,q_k,t_{k+1},q_{k+1})+(t_kt_{k1})D_4𝕃(t_{k1},q_{k1},t_k,q_k)=0`$ (2.3)
for all $`k\{1,2,\mathrm{},N1\}`$. Here $`D_i`$ denotes the partial derivative of $`𝕃`$ with respect to the $`i`$th argument. Notice that $`t_k`$ are parameters here. The reason we still call them arguments is that we want to keep the notations consistent since in next section they will play the role of variables.

Now we prove that the discrete flow $`\mathrm{\Phi }`$ preserves a discrete version of the Lagrange two form

$`\omega _L^v=d\theta _L^v,\theta _L^v={\displaystyle \frac{L}{\dot{q}}}dq.`$ (2.4)
We do this directly from the variational point of view, consistent with the continuous case \[References\].

As in continuous case, we calculate $`d𝕊`$ for variations with varied endpoints.

$`d𝕊(q_0,q_1,\mathrm{},q_N)(\delta q_0,\delta q_1,\mathrm{},\delta q_N)`$
$`={\displaystyle \frac{d}{dϵ}}|_{ϵ=0}𝕊(q_0+ϵ\delta q_0,q_1+ϵ\delta q_1,\mathrm{},q_N+ϵ\delta q_N)`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}\left(D_2𝕃(t_k,q_k,t_{k+1},q_{k+1})\delta q_k+D_4𝕃(t_k,q_k,t_{k+1},q_{k+1})\delta q_{k+1}\right)\left(t_{k+1}t_k\right)`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}D_2𝕃(t_k,q_k,t_{k+1},q_{k+1})\left(t_{k+1}t_k\right)\delta q_k+{\displaystyle \underset{k=1}{\overset{N}{}}}D_4𝕃(t_{k1},q_{k1},t_k,q_k)\left(t_kt_{k1}\right)\delta q_k`$
$`={\displaystyle \underset{k=1}{\overset{N1}{}}}\left(D_2𝕃(t_k,q_k,t_{k+1},q_{k+1})\left(t_{k+1}t_k\right)+D_4𝕃(t_{k1},q_{k1},t_k,q_k)\left(t_kt_{k1}\right)\right)\delta q_k`$
$`+D_2𝕃(t_0,q_0,t_1,q_1)\left(t_1t_0\right)\delta q_0+D_4𝕃(t_{N1},q_{N1},t_N,q_N)\left(t_Nt_{N1}\right)\delta q_N.`$ (2.5)
We can see that the last two terms in (2.5) come from the boundary variations. Based on the boundary variations, we can define two one forms on $`M\times M`$

$`\theta _𝕃^v(t_k,q_k,t_{k+1},q_{k+1})=D_2𝕃(t_k,q_k,t_{k+1},q_{k+1})(t_{k+1}t_k)dq_k,`$ (2.6)
and

$`\theta _𝕃^{v+}(t_k,q_k,t_{k+1},q_{k+1})=D_4𝕃(t_k,q_k,t_{k+1},q_{k+1})(t_{k+1}t_k)dq_{k+1}.`$ (2.7)
Here we have used the notation in \[References\]. We regard the pair $`(\theta _𝕃^v,\theta _𝕃^{v+})`$ as being the discrete version of the Lagrange one form $`\theta _L^v`$ in (2.4).

Now we parameterize the solutions of the discrete variational principle by the initial condition $`(q_0,q_1)`$, and restrict $`𝕊`$ to that solution space. Then Eq. (2.5) becomes

$`d`$ $`𝕊(q_0,q_1,\mathrm{},q_N)(\delta q_0,\delta q_1,\mathrm{},\delta q_N)`$
$`=\theta _𝕃^{}(t_0,q_0,t_1,q_1)(\delta q_0,\delta q_1)+\theta _𝕃^+(t_{N1},q_{N1},t_N,q_N)(\delta q_{N1},\delta q_N)`$
$`=\theta _𝕃^{}(t_0,q_0,t_1,q_1)(\delta q_0,\delta q_1)+(\mathrm{\Phi }^{N1})^{}\theta _𝕃^+(t_0,q_0,t_1,q_1)(\delta q_0,\delta q_1).`$ (2.8)
From (2.8), we can obtain

$`d𝕊=\theta _𝕃^v+(\mathrm{\Phi }^{N1})^{}\theta _𝕃^{v+}.`$ (2.9)
The Eq. (2.9) holds for arbitrary $`N>1`$. Taking N=2 leads to

$`d𝕊=\theta _𝕃^v+\mathrm{\Phi }^{}\theta _𝕃^{v+}.`$ (2.10)
By taking exterior differentiation of (2.10), we obtain

$`\mathrm{\Phi }^{}(d\theta _𝕃^{v+})=d\theta _𝕃^v.`$ (2.11)
From the definition of $`\theta _𝕃^{}`$ and $`\theta _𝕃^+`$, we know that

$`\theta _𝕃^v+\theta _𝕃^{v+}=d((t_{k+1}t_k)𝕃),`$
which means $`d\theta _𝕃^{v+}=d\theta _𝕃^v`$.

Defining

$`\omega _𝕃^v=:d\theta _𝕃^{v+}=d\theta _𝕃^v,`$ (2.12)
we show from (2.11) that the discrete flow $`\mathrm{\Phi }`$ preserves the discrete Lagrange two form $`\omega _𝕃^v`$:

$`\mathrm{\Phi }^{}(\omega _𝕃^v)=\omega _𝕃^v.`$ (2.13)
Let us consider an example. For the classical Lagrangian

$`L(t,q,\dot{q})={\displaystyle \frac{1}{2}}\dot{q}^2V(q),`$ (2.14)
we choose the discrete Lagrangian $`𝕃(t_k,q_k,t_{k+1},q_{k+1})`$ as

$`𝕃(t_k,q_k,t_{k+1},q_{k+1})={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}\right)^2V\left({\displaystyle \frac{q_{k+1}+q_k}{2}}\right).`$ (2.15)
The discrete Euler-Lagrange equation (2.3) becomes

$`\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}{\displaystyle \frac{q_kq_{k1}}{t_kt_{k1}}}\right)+{\displaystyle \frac{V^{}(\overline{q}_k)(t_{k+1}t_k)+V^{}(\overline{q}_{k1})(t_kt_{k1})}{2}}=0,`$ (2.16)
which preserves the Lagrange two form

$`\left({\displaystyle \frac{1}{t_{k+1}t_k}}+{\displaystyle \frac{t_{k+1}t_k}{4}}V^{\prime \prime }(\overline{q}_k)\right)dq_{k+1}dq_k,`$ (2.17)
where $`\overline{q}_k=\frac{q_k+q_{k+1}}{2},\overline{q}_{k1}=\frac{q_{k1}+q_k}{2}.`$

If we take fixed variables $`t_{k+1}t_k=t_kt_{k1}=h`$, then (2.16) becomes

$`{\displaystyle \frac{q_{k+1}2q_k+q_{k1}}{h^2}}+{\displaystyle \frac{V^{}(\overline{q}_k)+V^{}(\overline{q}_{k1})}{2}}=0,`$
which preserves the Lagrange two form

$`\left({\displaystyle \frac{1}{h}}+{\displaystyle \frac{h}{4}}V^{\prime \prime }(\overline{q}_k)\right)dq_{k+1}dq_k.`$
### 2.2 Discrete mechanics based on finite element methods

Now we consider discrete mechanics based on finite element methods. Let us go back to the variation problem of the action factional (1.1). The finite element method is an approximate method for solving the variation problem. Instead of solving the variation problem in the space $`𝒞^2([a,b])`$, the finite element method solves the problem in a subspace $`V_{h^m}([a,b])`$ of $`𝒞^2([a,b])`$. $`V_{h^m}([a,b])`$ consists of piecewise $`m`$-degree polynomials interpolating the curves $`q(t)𝒞^2([a,b])`$.

First, we consider the piecewise linear interpolation. Given a partition of \[a, b\]

$`a=t_0<t_1<\mathrm{}<t_k<\mathrm{}<t_{N1}<t_N=b,`$
the intervals $`I_k=[t_k,t_{k+1}]`$ are called elements. $`h_k=t_{k+1}t_k`$. $`V_h([a,b])`$ consists of piecewise linear function interpolating $`q(t)`$ at $`(t_k,q_k),k=0,1,\mathrm{},N`$. Now we derive the expressions of $`q_h(t)V_h([a,b])`$. First we construct the basis functions $`\phi _k(t)`$, which are piecewise linear function on $`[a,b]`$ satisfying $`\phi _k(t_i)=\delta _k^i,i,k=0,1,\mathrm{},N`$.

$`\phi _0(t)=\{\begin{array}{cc}1\frac{tt_0}{h_0},\hfill & t_0tt_1;\hfill \\ 0,\hfill & \text{otherwise};\hfill \end{array}\phi _N(t)=\{\begin{array}{cc}1+\frac{tt_N}{h_{N1}},\hfill & t_{N1}tt_N;\hfill \\ 0,\hfill & \text{otherwise};\hfill \end{array}`$ (2.18)
and for $`k=1,2,\mathrm{},N1`$,

$`\phi _k(t)=\{\begin{array}{cc}1+\frac{tt_k}{h_{k1}},\hfill & t_{k1}tt_k;\hfill \\ 1\frac{tt_k}{h_k},\hfill & t_ktt_{k+1};\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}`$ (2.19)
Using these basis functions, we obtain the expression $`q_hV_h([a,b])`$:

$`q_h(t)={\displaystyle \underset{k=0}{\overset{N}{}}}q_k\phi _k(t).`$ (2.20)
In the space $`V_h([a,b])`$, the action functional (1.1) becomes

$`S((t,q_h(t)))`$ $`={\displaystyle _a^b}L(t,q_h(t),\dot{q}_h(t))𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}{\displaystyle _{t_k}^{t_{k+1}}}L(t,{\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t)),{\displaystyle \frac{d}{dt}}{\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t)))𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}𝕃(t_k,q_k,t_{k+1},q_{k+1})(t_{k+1}t_k),`$ (2.21)
where

$`𝕃(t_k,q_k,t_{k+1},q_{k+1})=`$ $`{\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}L(t,{\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t)),{\displaystyle \frac{d}{dt}}{\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t)))𝑑t`$
$`=`$ $`{\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}L(t,{\displaystyle \underset{i=k}{\overset{k+1}{}}}(q_i\phi _i(t)),{\displaystyle \frac{d}{dt}}{\displaystyle \underset{i=k}{\overset{k+1}{}}}(q_i\phi _i(t)))𝑑t.`$ (2.22)
Therefore, restricting to the subspace $`V_h([a,b])`$ of $`𝒞^2([a,b])`$, the original variational problem reduces to the extremum problem of the function (2.21) in $`q_k,k=0,1,\mathrm{},N`$. Notice that (2.21) is just one of the discrete actions (2.1). Thus, what remains to be done is just to perform the same calculations on (2.21) as on (2.1). We can then obtain the discrete Euler-Lagrange equation (2.3) that preserves the discrete Lagrange two form (2.12). Therefore, discrete mechanics based on finite element methods consists of two steps: first, use finite element methods to obtain a kind of discrete Lagrangian; second, use the method of Veselov mechanics to obtain the variational integrators.

Let us consider the previous example again. For the classical Lagrangian (2.14), we choose the discrete Lagrangian $`𝕃(t_k,q_k,t_{k+1},q_{k+1})`$ as

$`𝕃(t_k,q_k,t_{k+1},q_{k+1})`$
$`={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{d}{dt}}{\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t))\right)^2V\left({\displaystyle \underset{i=0}{\overset{N}{}}}(q_i\phi _i(t))\right)\right)𝑑t`$
$`={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}\right)^2V\left({\displaystyle \frac{t_{k+1}t}{t_{k+1}t_k}}q_k+{\displaystyle \frac{tt_k}{t_{k+1}t_k}}q_{k+1}\right)\right)𝑑t`$
$`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}\right)^2F(q_k,q_{k+1}),`$ (2.23)
where

$`F(q_k,q_{k+1})={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}V\left({\displaystyle \frac{t_{k+1}t}{t_{k+1}t_k}}q_k+{\displaystyle \frac{tt_k}{t_{k+1}t_k}}q_{k+1}\right)𝑑t.`$ (2.24)
The discrete Euler-Lagrange equation (2.3) becomes

$`\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}{\displaystyle \frac{q_kq_{k1}}{t_kt_{k1}}}\right)+{\displaystyle \frac{F(q_k,q_{k+1})}{q_k}}\left(t_{k+1}t_k\right)+{\displaystyle \frac{F(q_{k1},q_k)}{q_k}}\left(t_kt_{k1}\right)=0,`$ (2.25)
which preserves the Lagrange two form

$`\left({\displaystyle \frac{1}{t_{k+1}t_k}}+(t_{k+1}t_k){\displaystyle \frac{^2F(q_k,q_{k+1})}{q_kq_{k+1}}}\right)dq_{k+1}dq_k.`$ (2.26)
Again, if we take fixed time steps $`t_{k+1}t_k=t_kt_{k1}=h`$, (2.25) becomes

$`{\displaystyle \frac{q_{k+1}2q_k+q_{k1}}{h^2}}+{\displaystyle \frac{F(q_k,q_{k+1})}{q_k}}+{\displaystyle \frac{F(q_{k1},q_k)}{q_k}}=0,`$
which preserves the Lagrange two form

$`\left({\displaystyle \frac{1}{h}}+h{\displaystyle \frac{^2F(q_k,q_{k+1})}{q_kq_{k+1}}}\right)dq_{k+1}dq_k.`$
Suppose $`q_k`$ is the solution of (2.25) and $`q(t)`$ is the solution of

$`{\displaystyle \frac{d^2q}{dt^2}}+{\displaystyle \frac{V(q)}{q}}=0,`$ (2.27)
then from the convergence theory of finite element methods \[References, References\], we have

$`q(t)q_h(t)Ch^2,`$ (2.28)
where $``$ is the $`L^2`$ norm. $`q_h(t)=_{k=0}^Nq_k`$, $`h=\underset{𝑘}{\text{max}}\{h_k\}`$ and $`C`$ is a constant independent of $`h`$.

If we use mid-point numerical integration formula in (2.24), we obtain

$`F(q_k,q_{k+1})=`$ $`{\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}V\left({\displaystyle \frac{t_{k+1}t}{t_{k+1}t_k}}q_k+{\displaystyle \frac{tt_k}{t_{k+1}t_k}}q_{k+1}\right)𝑑t`$
$``$ $`V\left({\displaystyle \frac{q_k+q_{k+1}}{2}}\right).`$
In this case, (2.25) is just (2.16). We can also use trapezoid formula or Simpson formula and so on to integrate (2.24) numerically and obtain another kind of discrete Lagrangian.

### 2.3 Discrete mechanics with Lagrangian of high order

Now we consider piecewise quadratic polynomial interpolation, which will result in a kind of discrete Lagrangian of high order. To this aim, we add an auxiliary node $`t_{k+\frac{1}{2}}`$ to each element $`I_k=[t_k,t_{k+1}]`$. There are two kinds of quadratic basis functions: $`\varphi _k(t)`$ for nodes $`t_k`$ and $`\varphi _{k+\frac{1}{2}}(t)`$ for $`t_{k+\frac{1}{2}}`$ that satisfy

$`\varphi _k(t_i)=\delta _i^k,\varphi _k(t_{i+\frac{1}{2}})=0,`$
$`\varphi _{k+\frac{1}{2}}(t_{i+\frac{1}{2}})=\delta _i^k,\varphi _{k+\frac{1}{2}}(t_i)=0,i,k=0,1,\mathrm{},N.`$
We list the basis functions as follows:

$`\varphi _0(t)=\{\begin{array}{cc}(\frac{2(tt_0)}{h_0}1)\left(\frac{tt_0}{h_0}1\right),\hfill & t_0tt_1;\hfill \\ 0,\hfill & \text{otherwise};\hfill \end{array}`$ (2.29)
$`\varphi _N(t)=\{\begin{array}{cc}(\frac{2(t_Nt)}{h_{N1}}1)\left(\frac{t_Nt}{h_{N1}}1\right),\hfill & t_{N1}tt_N;\hfill \\ 0,\hfill & \text{otherwise};\hfill \end{array}`$ (2.30)
and for $`k=1,2,\mathrm{},N1`$,

$`\varphi _k(t)=\{\begin{array}{cc}(\frac{2(t_kt)}{h_{k1}}1)\left(\frac{t_kt}{h_{k1}}1\right),\hfill & t_{k1}tt_k;\hfill \\ (\frac{2(tt_k)}{h_k}1)\left(\frac{tt_k}{h_k}1\right),\hfill & t_ktt_{k+1};\hfill \\ 0,\hfill & \text{otherwise};\hfill \end{array}`$ (2.31)
and for $`k=0,1,\mathrm{},N1`$,

$`\varphi _{k+\frac{1}{2}}(t)=\{\begin{array}{cc}4\frac{tt_k}{h_k}\left(1\frac{tt_k}{h_k}\right),\hfill & t_ktt_{k+1};\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}`$ (2.32)
Using these basis functions, we construct subspace $`V_{h^2}([a,b])`$ of $`𝒞^2([a,b])`$:

$`q_{h^2}(t)={\displaystyle \underset{k=0}{\overset{N}{}}}q_k\varphi _k(t)+{\displaystyle \underset{k=0}{\overset{N1}{}}}q_{k+\frac{1}{2}}\varphi _{k+\frac{1}{2}}(t),q_{h^2}(t)V_{h^2}([a,b]).`$ (2.33)
In the space $`V_{h^2}([a,b])`$, the action functional (1.1) becomes

$`S((t,q_{h^2}(t)))`$ $`={\displaystyle _a^b}L(t,q_{h^2}(t),\dot{q}_{h^2}(t))𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}{\displaystyle _{t_k}^{t_{k+1}}}L(t,q_{h^2}(t),\dot{q}_{h^2}(t))𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})(t_{k+1}t_k),`$ (2.34)
where

$`𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}L(t,q_{h^2}(t),\dot{q}_{h^2}(t))𝑑t.`$ (2.35)
For the discrete action (2.34), we have

$`d𝕊(q_0,q_{\frac{1}{2}},q_1,\mathrm{},q_{N1+\frac{1}{2}},q_N)(\delta q_0,\delta q_{\frac{1}{2}},\delta q_1,\mathrm{},\delta q_{N1+\frac{1}{2}},\delta q_N)`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}\left(D_2𝕃\left(𝒘_k\right)\delta q_k+D_3𝕃\left(𝒘_k\right)\delta q_{k+\frac{1}{2}}+D_5𝕃\left(𝒘_k\right)\delta q_{k+1}\right)\left(t_{k+1}t_k\right)`$
$`={\displaystyle \underset{k=0}{\overset{N1}{}}}\left(D_2𝕃\left(𝒘_k\right)\delta q_k+D_3𝕃\left(𝒘_k\right)\delta q_{k+\frac{1}{2}}\right)\left(t_{k+1}t_k\right)+{\displaystyle \underset{k=1}{\overset{N}{}}}D_5𝕃\left(𝒘_{k1}\right)\left(t_kt_{k1}\right)\delta q_k`$
$`={\displaystyle \underset{k=1}{\overset{N1}{}}}\left(D_2𝕃\left(𝒘_k\right)\left(t_{k+1}t_k\right)+D_5𝕃\left(𝒘_{k1}\right)\left(t_kt_{k1}\right)\right)\delta q_k+{\displaystyle \underset{k=0}{\overset{N1}{}}}D_3𝕃\left(𝒘_k\right)\delta q_{k+\frac{1}{2}}\left(t_{k+1}t_k\right)`$
$`+D_2𝕃\left(𝒘_0\right)\left(t_1t_0\right)\delta q_0+D_5𝕃\left(𝒘_{N1}\right)\left(t_Nt_{N1}\right)\delta q_N,`$ (2.36)
where $`𝒘_k=(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1}),k=0,1,\mathrm{},N1`$. From (2.36), we obtain the discrete Euler-Lagrange equation

$`D_2𝕃(𝒘_k)(t_{k+1}t_k)+D_5𝕃(𝒘_{k1})(t_kt_{k1})=0,`$ (2.37)
$`D_3𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})=0,`$ (2.38)
$`D_3𝕃(t_{k1},q_{k1},q_{k1+\frac{1}{2}},t_k,q_k)=0.`$ (2.39)
We can from (2.38) and (2.39) solve for $`q_{k+\frac{1}{2}}`$ and $`q_{k1+\frac{1}{2}}`$ respectively, then substitute them into (2.37) and finally solve for $`q_{k+1}`$. Therefore, the discrete Euler-Lagrange equation (2.37-2.39) still determines a discrete flow

$`\mathrm{\Psi }:M\times MM\times M`$
$`\mathrm{\Psi }(t_{k1},q_{k1},t_k,q_k)=(t_k,q_k,t_{k+1},q_{k+1}).`$
From (2.36), we can define two one forms

$`\mathrm{\Theta }_𝕃^v(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})=D_2𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})(t_{k+1}t_k)dq_k,`$
and

$`\mathrm{\Theta }_𝕃^{v+}(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})=D_5𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})(t_{k+1}t_k)dq_{k+1}.`$
Using the same method as before, we can prove that

$`\mathrm{\Psi }^{}(d\mathrm{\Theta }_𝕃^{v+})=d\mathrm{\Theta }_𝕃^v.`$ (2.40)
From the definition of $`\mathrm{\Theta }_𝕃^v`$ and $`\mathrm{\Theta }_𝕃^{v+}`$, we have

$`\mathrm{\Theta }_𝕃^v+\mathrm{\Theta }_𝕃^{v+}=d((t_{k+1}t_k)L)D_3L(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})dq_{k+\frac{1}{2}}.`$ (2.41)
From (2.38), we obtain $`D_3L(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})=0`$. Thus,

$`\mathrm{\Theta }_𝕃^v+\mathrm{\Theta }_𝕃^{v+}=d((t_{k+1}t_k)L),`$
which means

$`d\mathrm{\Theta }_𝕃^{v+}=d\mathrm{\Theta }_𝕃^v.`$ (2.42)
From (2.40) and (2.42), we arrive at

$`\mathrm{\Psi }^{}(\mathrm{\Omega }_𝕃^v)=\mathrm{\Omega }_𝕃^v,`$ (2.43)
where $`\mathrm{\Omega }_𝕃^v=d\mathrm{\Theta }_𝕃^{v+}`$.

For the classical Lagrangian (2.14), From (2.33) and (2.35), we obtain

$`𝕃(t_k,q_k,q_{k+\frac{1}{2}},t_{k+1},q_{k+1})`$
$`={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}\left({\displaystyle \frac{1}{2}}(\dot{q}_{h^2}(t))^2V(q_{h^2}(t))\right)𝑑t`$
$`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{3}}a^2(t_{k+1}^2+t_kt_{k+1}+t_k^2)+ab(t_k+t_{k+1})+b^2\right)`$
$`G(q_k,q_{k+\frac{1}{2}},q_{k+1}),`$ (2.44)
where

$`a={\displaystyle \frac{4}{h_k^2}}\left(q_k+q_{k+1}2q_{k+\frac{1}{2}}\right),`$
$`b={\displaystyle \frac{1}{h_k^2}}\left(4(t_k+t_{k+1})q_{k+\frac{1}{2}}(3t_k+t_{k+1})q_{k+1}(t_k+3t_{k+1})q_k\right),`$
and

$`G(q_k,q_{k+\frac{1}{2}},q_{k+1})={\displaystyle \frac{1}{t_{k+1}t_k}}{\displaystyle _{t_k}^{t_{k+1}}}V\left(q_kf_k(t)+q_{k+1}f_{k+1}(t)+q_{k+\frac{1}{2}}f_{k+\frac{1}{2}}(t)\right)𝑑t,`$
where

$`f_k(t)=({\displaystyle \frac{2(tt_k)}{h_k}}1)\left({\displaystyle \frac{tt_k}{h_k}}1\right),`$
$`f_{k+1}(t)=({\displaystyle \frac{2(t_{k+1}t)}{h_k}}1)\left({\displaystyle \frac{t_{k+1}t}{h_k}}1\right),`$
$`f_{k+\frac{1}{2}}(t)=4{\displaystyle \frac{tt_k}{h_k}}\left(1{\displaystyle \frac{tt_k}{h_k}}\right).`$
For the discrete Lagrangian (2.44), the discrete Euler-Lagrangian equation (2.37-2.39) becomes

$`a_1q_{k1}+a_2q_k+a_3q_{k+1}+a_4q_{k\frac{1}{2}}+a_4q_{k+\frac{1}{2}}d_1h_kd_2h_{k1}=0,`$ (2.45)
$`{\displaystyle \frac{8}{3h_k^2}}\left(q_k+q_{k+1}2q_{k+\frac{1}{2}}\right){\displaystyle \frac{G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}}}=0,`$ (2.46)
$`{\displaystyle \frac{8}{3h_{k1}^2}}\left(q_{k1}+q_k2q_{k1+\frac{1}{2}}\right){\displaystyle \frac{G(q_{k1},q_{k1+\frac{1}{2}},q_k)}{q_{k1+\frac{1}{2}}}}=0,`$ (2.47)
where

$`a_1={\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{h_{k1}}},a_2={\displaystyle \frac{7}{3}}\left({\displaystyle \frac{1}{h_{k1}}}+{\displaystyle \frac{1}{h_k}}\right),a_3={\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{h_k}},`$
$`a_4={\displaystyle \frac{8}{3}}{\displaystyle \frac{1}{h_{k1}}},a_5={\displaystyle \frac{8}{3}}{\displaystyle \frac{1}{h_k}},`$
$`d_1={\displaystyle \frac{G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_k}},d_2={\displaystyle \frac{G(q_{k1},q_{k1+\frac{1}{2}},q_k)}{q_k}}.`$
The solution of (2.45-2.47) preserves the Lagrange two form

$`\left({\displaystyle \frac{1}{3h_k}}h_k{\displaystyle \frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_kq_{k+1}}}M\right)dq_kdq_{k+1},`$ (2.48)
where

$$M=\frac{\left(\frac{16}{3h_k}+h_k\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}q_k}\right)\left(\frac{16}{3h_k}+h_k\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}q_k}\right)}{\left(\frac{32}{3h_k}h_k\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}^2}\right)}.$$

If we take the fixed time steps $`h_{k1}=h_k=h`$, then (2.45-2.47)) become

$`{\displaystyle \frac{q_{k1}8q_{k\frac{1}{2}}+14q_k8q_{k+\frac{1}{2}}+q_{k+1}}{3h^2}}d_1h_kd_2h_{k1}=0,`$ (2.49)
$`{\displaystyle \frac{8}{3h^2}}\left(q_k+q_{k+1}2q_{k+\frac{1}{2}}\right){\displaystyle \frac{G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}}}=0,`$ (2.50)
$`{\displaystyle \frac{8}{3h^2}}\left(q_{k1}+q_k2q_{k1+\frac{1}{2}}\right){\displaystyle \frac{G(q_{k1},q_{k1+\frac{1}{2}},q_k)}{q_{k1+\frac{1}{2}}}}=0,`$ (2.51)
which preserve

$`\left({\displaystyle \frac{1}{3h}}h{\displaystyle \frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_kq_{k+1}}}M\right)dq_kdq_{k+1},`$ (2.52)
where

$$M=\frac{\left(\frac{16}{3h_k}+h\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}q_k}\right)\left(\frac{16}{3h}+h\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}q_k}\right)}{\left(\frac{32}{3h}h\frac{^2G(q_k,q_{k+\frac{1}{2}},q_{k+1})}{q_{k+\frac{1}{2}}^2}\right)}.$$

Suppose $`q_k`$ is the solution of (2.45-2.47)) and $`q(t)`$ is the solution of (2.27), then from the convergence theory of finite element methods \[References, References\], we have

$`q(t)q_{h^2}(t)Ch^3,`$ (2.53)
where

$$q_{h^2}(t)=\underset{k=0}{\overset{N}{}}q_k\varphi _k(t)+\underset{k=0}{\overset{N1}{}}q_{k+\frac{1}{2}}\varphi _{k+\frac{1}{2}}(t),$$

$`h=\underset{𝑘}{\text{max}}\{h_k\}`$ and $`C`$ is a constant independent of $`h`$.

## 3 Time steps as variables

In §2, the time steps $`t_k`$ play the role of parameters. They are determined beforehand according to some requirements. In fact, we can also regard $`t_k`$ as variables and the variation of the discrete action with respect to $`t_k`$ yields the discrete energy conservation law. This fact was first observed by Lee \[ References, References, References\]. The symplecticity of the resulting integrators was investigated in \[References, References\]. These results also apply to the discrete mechanics based on finite element methods.

We regard $`t_k`$ as variables and calculate the variation of the discrete action (2.1) as follows

$`d𝕊(t_0,q_0,\mathrm{},t_N,q_N)(\delta t_0,\delta q_0,\mathrm{},\delta t_N,\delta q_N)`$
$`={\displaystyle \frac{d}{dϵ}}|_{ϵ=0}𝕊(t_0+ϵ\delta t_0,q_0+ϵ\delta q_0,\mathrm{},t_N+ϵ\delta t_N,q_N+ϵ\delta q_N)`$
$`={\displaystyle \underset{k=1}{\overset{N1}{}}}\left[D_2𝕃(𝒘_k)(t_{k+1}t_k)+D_4𝕃(𝒘_{k1})(t_kt_{k1})\right]\delta q_k`$
$`+{\displaystyle \underset{k=1}{\overset{N1}{}}}\left[D_1𝕃(𝒘_k)(t_{k+1}t_k)+D_3L(𝒘_{k1})(t_kt_{k1})+L(𝒘_{k1})L(𝒘_k)\right]\delta t_k`$
$`+D_2L(𝒘_0)(t_1t_0)\delta q_0+D_4L(𝒘_{N1})(t_Nt_{N1})\delta q_N`$
$`+\left[D_1L(𝒘_0)(t_1t_0)L(𝒘_0)\right]\delta t_0`$
$`+\left[D_3L(𝒘_{N1})(t_Nt_{N1})+L(𝒘_{N1})\right]\delta t_N,`$ (3.1)
where $`𝒘_k=(t_k,q_k,t_{k+1},q_{k+1}),k=0,1,\mathrm{},N1`$. From (3.1), we see that the variation $`\delta q_k`$ yields the discrete Euler-Lagrange equation

$`D_2𝕃(𝒘_k)(t_{k+1}t_k)+D_4𝕃(𝒘_{k1})(t_kt_{k1})=0`$ (3.2)
and the variation $`\delta t_k`$ yields the discrete energy evolution equation

$`D_1𝕃(𝒘_k)(t_{k+1}t_k)+D_3L(𝒘_{k1})(t_kt_{k1})+L(𝒘_{k1})L(𝒘_k)=0,`$ (3.3)
which is a discrete version of (1.9). For a conservative $`L`$, (3.3) becomes the discrete energy conservation law.

From the boundary terms in (3.1), we can define two one-forms

$`\theta _𝕃^{}(𝒘_k)=(D_1𝕃(𝒘_k)(t_{k+1}t_k)𝕃(𝒘_k))dt_k+D_2𝕃(𝒘_k)(t_{k+1}t_k)dq_k,`$ (3.4)
and

$`\theta _𝕃^+(𝒘_k)=(D_3𝕃(𝒘_k)(t_{k+1}t_k)+𝕃(𝒘_k))dt_{k+1}+D_4𝕃(𝒘_k)(t_{k+1}t_k)dq_{k+1}.`$ (3.5)
These two one forms are the discrete version of the extended Lagrange one form (1.12).

Unlike the continuous case, the solution of (3.2) does not satisfy (3.3) in general. Therefore, we must solve (3.2) and (3.3) simultaneously. Using the same method in §2, we can show that the coupled integrator

$`\begin{array}{cc}& D_2𝕃(𝒘_k)(t_{k+1}t_k)+D_4𝕃(𝒘_{k1})(t_kt_{k1})=0,\hfill \\ & D_1𝕃(𝒘_k)(t_{k+1}t_k)+D_3L(𝒘_{k1})(t_kt_{k1})+L(𝒘_{k1})L(𝒘_k)=0\hfill \end{array}`$ (3.6)
preserves the extended Lagrange two form $`\omega _𝕃=d\theta _𝕃^+`$.

For the discrete Lagrangian (2.23), (3.6) becomes

$`\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}{\displaystyle \frac{q_kq_{k1}}{t_kt_{k1}}}\right)+{\displaystyle \frac{F\left(𝒘_k\right)}{q_k}}h_k+{\displaystyle \frac{F\left(𝒘_{k1}\right)}{q_k}}h_{k1}=0,`$
$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{q_{k+1}q_k}{t_{k+1}t_k}}\right)^2+F\left(𝒘_k\right){\displaystyle \frac{F\left(𝒘_k\right)}{t_k}}h_k={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{q_kq_{k1}}{t_kt_{k1}}}\right)^2+F\left(𝒘_{k1}\right)+{\displaystyle \frac{F\left(𝒘_{k1}\right)}{t_k}}h_{k1}.`$
For the kind of high order discrete Lagrangian, we can obtain similar formulas.

## 4 Conclusions

In this paper, we have presented the discrete mechanics based on finite element methods. The finite element methods provide a natural way to obtain a kind of discrete Lagrangian. Using piecewise quadratic functions, we have constructed a high order discrete Lagrangian and derived the corresponding variational integrator (the discrete Euler-Lagrange equation). The symplecticity of the high order variational integrator has also been proved.

The time steps can play the role of parameters or the role of variables. We can obtain discrete energy conservation law when the time steps are regarded as variables.

Finite element methods are well-established in numerical mathematics. Many results of finite element methods can be applied to the discrete mechanics presented in this paper such as convergence theory, error estimates as well as the solving method of the variational integrators.

Recently, it has been proved \[References\] that the simplectic structure preserves not only on the phase flow but also on the flows with respect to symplectic vector fields as long as certain cohomological condition is satisfied in both continuous and discrete cases. This should be able to extend to the cases in this paper. We will study this issue elsewhere.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (grant Nos. 90103004, 10171096) and the National Key Project for Basic Research of China (G1998030601).


format: text

text:
# Spatial distribution of emission in Unidentified Infrared Bands from Midcourse Space Experiment Survey

## 1 Introduction

The near to mid infrared spectrum originating from the interstellar medium of the Galactic star forming regions consists of various features in addition to a continuum. The continuum emission is attributed to the thermal emission from interstellar dust component and almost all narrow emission lines have been identified with atomic nebular, molecular, etc transitions originating from the interstellar gas component. In addition, several broader features have been detected which have been identified with features due to the solid state material of the dust (e.g. silicate absorption features at $``$ 10 and $``$ 18 $`\mu `$m etc). There exists a class of broad emission features, sometimes called “Unidentified Infrared emission Bands” (UIBs; at 3.3, 6.2, 7.7, 8.6, 11.3, 12.7 $`\mu `$m), identity of whose carriers and the emission mechanisms are still a subject of active research. Some of these bands are widely believed to be characteristic of the bending and stretching modes of C=C and C–H bonds in aromatic molecules (e.g. fluorescent emission of Polycyclic Aromatic Hydrocarbons or PAHs; Leger & Puget 1984, Allamandolla et al. 1985). However, other contenders also exist in literature (e.g. amorphous materials with aromatic hydrocarbon; Sakata et al. 1984, Borghesi et al. 1987). The study of large scale distribution of emission in the UIBs in the Galactic plane in general and selected star forming regions in particular would be important in understanding details of their emission mechanism.

The recent mission Infrared Space Observatory (ISO), in particular the imaging camera ISOCAM, has made it possible to study selected Galactic star forming regions in the UIBs. The ISOCAM instrument had several filters with the passbands selected to cover these UIBs (and also neighbouring continuum) so that emission in these individual features can be measured very precisely. However, it is unreasonable to expect Galactic plane surveys in UIB emission using the ISOCAM, since its primary objective was to achieve best possible (nearly diffraction limited) angular resolution in studies of individual astrophysical sources. As a result, largest single image from ISOCAM covers $``$ 3′$`\times `$ 3′.

With the advent of the Midcourse Space Experiment (MSX), new possibilities have emerged. The SPIRIT III instrument onboard MSX spacecraft, has surveyed the entire Galactic Plane ($`|b|<5`$°) in four mid infrared bands centered around 8.3, 12.13, 14.65 and 21.34 $`\mu `$m with an angular resolution $``$ 18″ (Price et al. 2001). These four bands are referred to as $`A`$, $`C`$, $`D`$ and $`E`$ spectral bands of MSX respectively (the SPIRIT III instrument also had two additional narrower bands at 4.29 and 4.35 $`\mu `$m, called $`B_1`$ and $`B_2`$ bands). The usefulness of the MSX survey to the study of diffuse interstellar medium and global characteristics has already been demonstrated (Cohen & Green 2001; Cohen 1999). The MSX band $`A`$ includes the dominant UIB features at 6.2, 7.7 and 8.7 $`\mu `$m. Similarly the MSX band $`C`$ includes the UIB features at 11.3 and 12.7 $`\mu `$m.

Whereas, the ISOCAM provides imaging capability in narrower spectral bands at higher angular resolution (3″ or 6″ ) of selected regions, the MSX survey covers the entire Galactic plane in four broader bands (two of these covering several UIBs in addition to the continuum). Making use of this complementarity of ISOCAM vis a vis MSX, the following scheme has been explored to study large scale emission in the UIBs in the Galaxy :

* Model each picture element of MSX images with an integrated emission in UIB features superposed on a gray body continuum spectrum under some reasonable assumptions. The best fit solution (for each pixel) provides a measure of the UIB emission locally.
* Test the reliability of the above scheme by comparing the results of some selected Galactic star forming regions which have been studied using the ISOCAM and whose emission in individual UIBs have been quantified (e.g. Zavagno & Ducci 2001).
* The comparison should cover not only qualitative (e.g. structural details / morphology) aspects, but also quantitative correlation between the integrated estimate of UIB emission from our scheme and the ISOCAM results.

The section 2 describes the modelling scheme in detail, and the section 3 presents the results for the sample of six Galactic star forming regions using the MSX survey data. In section 4, a comparison between our results have been made with those from the ISOCAM by Zavagno and Ducci (2001). The conclusions are summarised in section 5.

## 2 The Scheme

The publicly available MSX Galactic plane survey radiance images (Infrared Processing and Analysis Center at http://irsa.ipac.caltech.edu/applications/MSX) in the four bands at 8.3 ($`A`$), 12.1 ($`C`$), 14.7 ($`D`$) and 21.3 ($`E`$) $`\mu `$m are gridded in 6″$`\times `$ 6″ pixel, though the true angular resolution is $``$ 18.3″ (Price et al. 2001), in the unit of $`Wm^2sr^1`$. The zodiacal background has already been subtracted out from these MSX survey maps. The spectrum emitted from each pixel is assumed to be a combination of a thermal continuum (modified Planck function or gray body) and the total radiance due to the relevant UIB features within the MSX band.

$$R_i=r_i^{UIBs}+(1e^{\tau _\nu })\times B_\nu (T)\times RSR_i(\nu )𝑑\nu $$
$`(1)`$
$$i=A,C,D,E$$

where $`R_i`$ are the measured radiances in the MSX bands. $`B_\nu (T)`$ is the Planck function, and the term in parenthesis emulates the gray body spectrum emitted by the dust grains. $`\tau _\nu `$ is the optical depth due to the interstellar dust component at the frequency $`\nu `$. The functions $`RSR_i(\nu )`$ represent the normalized relative spectral responses of the four MSX bands (Egan et al. 1999).

Since the range of frequencies covered by the MSX bands is limited, we assume a power law dependence of the dust emissivity on frequency, viz.,

$$\tau _\nu =\tau _{10}\times \left(\frac{\nu }{\nu _{10}}\right)^\beta .$$
$`(2)`$
here $`\tau _{10}`$ is the optical depth at 10 $`\mu `$m and $`\nu _{10}`$ the frequency corresponding to wavelength 10 $`\mu `$m. The value of $`\beta `$ is a constant determined from the type of dust assumed. The effect of varying $`\beta `$ is discussed later. The $`r_i^{UIBs}`$ are the modelled total radiances in UIBs within the $`i`$-th MSX band. Since there are no known UIBs within the bands $`D`$ and $`E`$ of MSX, $`r_D^{UIBs}`$ = 0 and $`r_E^{UIBs}`$ = 0. In addition, it has been assumed that the total radiance due to the UIBs in band $`C`$ is proportional to that in band $`A`$, viz.,

$$r_C^{UIBs}=\alpha \times r_A^{UIBs}.$$
$`(3)`$
Here $`\alpha `$ is held fixed to a reasonable value (based on available observational data and understanding) for all pixels of all star forming regions studied, though the effect of changing the value of $`\alpha `$ has also been discussed later.

For each pixel on the sky with sufficient signal to noise ratio in each of the four MSX bands (implemented by map dynamic range cuts), we solve the set of four equations (eq. 1), for the three unknown variables viz., $`T`$, $`\tau _{10}`$ and $`r_A^{UIBs}`$. A non-linear chi-square minimization scheme based on the finite difference Levenberg-Marquardt algorithm has been used for this purpose. The integrals in equation 1 are evaluated numerically using a cautious adaptive Romberg extrapolation method. In order to ensure that the best solution obtained indeed corresponds to a global minimum of chi-square, the computations are repeated for 125 different sets of initial guesses comprising of 5 values each of $`T`$, $`\tau _{10}`$ and $`r_A^{UIBs}`$. The grid of initial guess values for these three variables have been selected to cover a wide range of physical situations (e.g. $`T`$ ranging between 50 and 800 K; $`\tau _{10}`$ between $`10^6`$ and $`10^2`$; $`r_A^{UIBs}`$ between $`10^7`$ and $`10^3`$ $`Wm^2sr^1`$).

Invariably, the same solution is obtained starting from almost all different sets of initial guesses. The procedure is repeated for all pixels of the MSX map resulting in spatial distributions of these three physical variables. Here we extensively use the map of $`r_A^{UIBs}`$, which can be compared with the measurements from ISOCAM.

## 3 Results

### 3.1 Results from ISOCAM

#### 3.1.1 Available results from the literature & sample selection

From the literature we have selected the work of Zavagno and Ducci (2001; hereafter ZD) which is based on ISOCAM measurements, for detailed comparison with results from our scheme of extracting emission in UIBs from the MSX data. The reason for selecting ZD is that they very uniformly studied (using same set of filters) a reasonable sample size comprising of five Galactic compact H II regions, viz., Sharpless(Sh)-61, Sh-138, Sh-152, Sh-156 and Sh-186. They studied the entire 3–12 $`\mu `$m wavelength range accessible through ISOCAM. The ZD sample of compact H II regions are bright in IRAS 12 $`\mu `$m band and are known to be strong emitters of UIBs at 7.7, 8.6 and 11.3 $`\mu `$m but show no silicate absorption feature at 10 or 18 $`\mu `$m. In addition, they represent a sequence in equivalent stellar type of the main exciting star.

Results of ZD are based on imaging with 3″$`\times `$ 3″ pixel mode of ISOCAM covering 87″$`\times `$ 87″ regions around Sh-61, Sh-152 & Sh-186 and 174″$`\times `$ 87″ around Sh-138 & Sh-156. Based on images in SW1 (centre $`\lambda `$ : 3.57 $`\mu `$m; $`\lambda `$ range : 3.05–4.10 $`\mu `$m; Césarsky et al. 1996), SW2 (3.30; 3.20–3.40), LW4 (6.00; 5.50–6.50), LW6 (7.75; 7.00–8.50), LW8 (11.4; 10.7–12.0) filters and five selected wavelengths using the CVF, they have quantified the UIB fluxes (actually radiances) in the 3.3, 6.2, 7.7 and 11.2 $`\mu `$m features integrated over the mapped regions.

#### 3.1.2 Our estimation of 7.7 $`\mu `$m UIB feature emission

Using the publicly available ISOCAM data of the ZD sample sources (ISO Postcards from ISO Data Archive for General Users; http://www.iso.vilspa.esa.es), we have extracted the spatial distribution of emission in the 7.7 $`\mu `$m UIB feature, using a method similar (but not identical) to ZD. First of all the underlying continuum at 7.7 $`\mu `$m has been estimated from power law interpolation using the CVF images at 6.91 and 8.22 $`\mu `$m. Next, this continuum has been subtracted from the LW6 image and the resulting emission has been attributed to the 7.7 $`\mu `$m UIB (hereafter $`UIB_{7.7}`$ map). The resulting maps are presented and discussed in the next section.

### 3.2 Results obtained from the MSX survey

The MSX images of the five sources from ZD sample were processed following the scheme described above. The dust emissivity power law index $`\beta `$ has been taken to be 1.0, which is commonly used in literature for general interstellar grains in the mid infrared wavelength region relevant to MSX bands (Scoville and Kwan 1976, Savage and Mathis 1979: Mathis 1990). However, other values of $`\beta `$ have also been explored to ensure that the results obtained here are not sensitive to the choice of $`\beta `$, as discussed later.

#### 3.2.1 Choice of $`\alpha `$

The value of $`\alpha `$ represents the ratio of radiances in the UIBs within the passbands of MSX bands $`C`$ and $`A`$ (see eq. 3). These correspond to the features at 11.2 and 12.7 $`\mu `$m (in band $`C`$) and at 6.2, 7.7 and 8.6 $`\mu `$m (in band $`A`$). The value of $`\alpha `$ used here is based on the recent work by Verstraete et al. (2001) using the ISO-SWS. They have studied the 2.4 – 25 $`\mu `$m spectra of three selected bright Galactic interstellar regions where dense molecular gas is illuminated by stellar radiation. Their wavelength range covers all the four MSX bands adequately. In addition, their choice of the three regions, viz., NGC 2023, Orion bar and M17-SW spans a wide range of excitation parameter (flux as well as hardness of the radiation field). The spectral resolution of the ISO-SWS measurements used by them was either 200 or 500, significantly large to resolve individual UIBs. The average value of $`\alpha `$ for the above three regions representing different physical conditions, has been estimated by us from Fig.1 of Verstraete et al. (2001) to be 0.32 (with a very small dispersion). The effect of choosing different values of $`\alpha `$ have been studied by us, which is discussed later.

From the study of Verstraete et al. (2001) it is clear that the contribution of forbidden ionic lines and molecular rotational lines, to the radiance within the MSX bands is negligible compared to the UIBs and the underlying continuum. Hence presence of these narrow lines should not affect our scheme of extracting the radiance due to the UIBs from the MSX data.

#### 3.2.2 Maps of extracted UIB radiance

The extent of the angular region considered around each target was determined by the available dynamic ranges in the four MSX bands (i.e. only pixels satisfying the dynamic range condition in each band, were modelled). The usable dynamic range, $`UDR`$, for each band was defined from the frequency distribution of the radiance values, $`f(R)`$, in the corresponding map in the following manner : $`UDR`$ = $`R_{max}`$ / ($`R_{median}`$ \+ $`R_{1/2}^{}`$), where $`R_{max}`$ is the brightest pixel value, $`R_{median}`$ is the median value of $`R`$ as determined from $`f(R)`$, and $`R_{1/2}^{}`$ represents the brightness value satisfying $`f(R_{1/2}^{})`$ = $`f(R_{median})/2`$ and $`R_{1/2}^{}<R_{median}`$. For a purly Gaussian distribution, the above translates to using the brightest $``$ 15 % pixels of the full sample. It may be noted here that any particular choice of dynamic range only changes the outer boundary of the region where our scheme is applied by either including or excluding these pixels, without affecting any numerical results for other pixels. The details of the dynamic range used for the five regions are presented in Table 1.

The resulting integrated UIB radiance maps extracted by us from the MSX data for Sh-61, Sh-138, Sh-152, Sh-156 and Sh-186 are displayed in Figures 1 to 5 respectively as isophot contour plots. The sizes of these maps are selected based on the regions covered by the study by ZD. The peak UIB radiance values are also listed in Table 1.

Morphologically the extracted UIB radiance maps are very similar to the corresponding maps of thermal continuum emission from the dust grains, in general (at the scale of MSX resolution). However, there are differences between the spatial distribution of UIB emission and the other modelled parameters, viz., the dust temperature and the optical depth corresponding to the thermal continuum emission. As one example, the map of dust optical depth ($`\tau _{10}`$, at 10 $`\mu `$m) for the Sharpless 152 region is presented in Figure 6, which can be compared with Figure 3.

## 4 Comparison between MSX and ISOCAM results

### 4.1 Morphological similarities

It is instructive to compare our maps from modelling of MSX data with those based on the ISOCAM data. Most relevant ISOCAM filters for this comparison are LW4, LW6 and LW8 covering the 6.2, 7.7 and 11.2 $`\mu `$m UIBs. The former two features contribute to the UIB radiance in MSX band $`A`$ and the last one in MSX band $`C`$.

The ISOCAM based $`UIB_{7.7}`$ maps have been superposed (grey scale) on our integrated UIB map extracted from MSX maps, viz. $`r_A^{UIBs}`$, (hereafter $`UIB_A`$), in Figures 7 to 11 corresponding to Sh-61, Sh-138, Sh-152, Sh-156 and Sh-186 respectively. For a more general comparison involving other UIBs, we refer to Figure 2 of ZD. In principle, our $`UIB_{7.7}`$ map for each source must resemble LW6 map of ZD, which indeed is the case.

A marked morphological similarity between the contour plots corresponding to all these three ISOCAM filters (LW4, LW6 & LW8) for all the five sources, support our assumption that $`r_C^{UIBs}`$ and $`r_A^{UIBs}`$ are correlated (equation 3). Next, we comment qualitatively about structural similarities between the spatial distribution of emission in individual UIBs as obtained from ISOCAM (our $`UIB_{7.7}`$; LW4, LW6 & LW8 maps of ZD), vis a vis the total emission due to all UIBs within the MSX band $`A`$ ($`UIB_A`$). It may be noted that whereas ZD’s maps represent intensity while ours are radiance. For qualitative comparison of structural details, it may be acceptable, however quantitive comparison is made in identical units later in this subsection.

#### 4.1.1 Sharpless 61

The peak position for Sh-61 in $`UIB_A`$ map matches exactly with that of $`UIB_{7.7}`$ (Fig. 7 here) / LW6 and LW4 (Fig. 2 of ZD). General structure of the $`UIB_A`$ contours is also very similar to the maps LW4, LW6 as well as LW8 of ZD, with extension along E-W. We have detected an additional emission region due S-W from the main peak in the $`UIB_A`$ map, which is perhaps barely below the lowest contour displayed in LW4 or LW6. However, our $`UIB_{7.7}`$ map gives some hint of the same despite being near the edge of the ISOCAM field.

#### 4.1.2 Sharpless 138

For Sh-138 (Fig. 8), the $`UIB_{7.7}`$ map has been able to resolve two components along N-S, and the $`UIB_A`$ map shows the peak position as well as the contour shape consistent with that (note that the former has an angular resolution is $``$ 18″ while ISOCAM’s is 3″). The contours in LW4, LW6 and LW8 maps are remarkably similar with a main peak and extension due north.

#### 4.1.3 Sharpless 152

Our $`UIB_{7.7}`$ (Fig. 9) as well as the LW4, LW6 and LW8 maps for Sh-152 show very rich (but again similar between these filters) structure at an angular scale beyond the scope of MSX’s resolution. Still, our $`UIB_A`$ map for Sh-152 reproduces the major structures consistent with the above, viz., the curved shape extending approximately along S-E to N-W direction. The entire $`UIB_A`$ map seems to be shifted by $``$ 10″ with respect to $`UIB_{7.7}`$ map. Additional interesting structures are detected in our map (see Fig. 3; due north and due south with respect to the main peak) which lie outside the region mapped by ISOCAM.

#### 4.1.4 Sharpless 156

In the case of Sh-156 also the $`UIB_{7.7}`$, LW4, LW6 and LW8 maps show structural details at smaller angular scale which are unresolved in our $`UIB_A`$ map (Fig. 10). But the curvature and density of contours in the latter indicate a morphology consistent with the above ISOCAM maps. The positions of respective peaks match very well in all these maps. UIB Emission is detected from a larger extended region all around the region mapped by ISOCAM, with one isolated secondary peak due west of the main peak (Fig. 4).

#### 4.1.5 Sharpless 186

The $`UIB_{7.7}`$ as well as the $`UIB_A`$ maps have resolved the Sh-186 region into two major peaks separated mainly along S-E to N-W line (Fig. 11). The stronger component being the S-E one in all maps. Similar structure is seen in the LW4, LW6 and LW8 maps too. The position of the peaks and other structural features also match very well (within $`<`$ 10″) between the $`UIB_A`$ and all ISOCAM maps. The $`UIB_A`$ map shows much more extended emission well outside the region mapped by ZD (Fig. 5).

### 4.2 Quantitative correlation

In order to make quantitative comparison, the $`UIB_A`$ radiance maps of the five regions have been integrated over the same corresponding regions as imaged by ZD, to get I(UIB$`{}_{}{}^{MSX}{}_{A}{}^{}`$). ZD have tabulated the solid angle integrated UIB fluxes F(3.3$`\mu `$m), F(6.2$`\mu `$m), F(7.7$`\mu `$m) and F(11.2$`\mu `$m) in their Table 4. We have compared their \[F(6.2$`\mu `$m)+F(7.7$`\mu `$m)\]= I(ISOCAM4+6) for each region with our integrated $`UIB_A`$ radiances (see Table 2). Surprizingly, an extremely tight linear correlation has been found between I(ISOCAM4+6) and I(UIB$`{}_{}{}^{MSX}{}_{A}{}^{}`$)! It has been found that the ratio, $`_A`$ = I(UIB$`{}_{}{}^{MSX}{}_{A}{}^{}`$) / I(ISOCAM4+6) = 2.29 $`\pm `$ 0.07 (mean value and the error on the mean). This is despite several simplifying assumptions made in our analysis of the MSX data. This is indeed a remarkable finding, considering the very complex microscopic as well as macroscopic details that must go into deciding the amount of emission in the UIBs. Whereas the UIBs at 3.3, 8.6, 11.3 and 12.7 $`\mu `$m originate from vibrational modes of the aromatic C–H bond, the bands at 6.2 and 7.7 $`\mu `$m arise from the aromatic C–C bonds. There is strong evidence from the laboratory work on PAHs (carriers of UIBs), that the relative strengths of these features are very sensitive to the ionization state of the PAHs (Allamandola, Hudgins & Sandford 1999). Several observational studies of the spatial distribution of UIBs in well resolved H II regions support the same. Joblin et al. (1996) found that the ratio of 8.6 $`\mu `$m to 11.3 $`\mu `$m emission increases with the local far-UV (FUV) flux, which can ionize the PAH molecules. In a more recent study, Crété et al. (1999) found similar variations in the ratio of UIB features for the M 17 complex. One possible explanation for the above correlation is as follows : (i) the angular sizes of the sample of sources considered here are such that the intrinsic resolution of MSX ($``$ 18″) spatially averages out various excitation effects expected close to the source of FUV radiation field; and (ii) the sample spans a somewhat limited range of FUV luminosity ($``$ 2). It may be interesting to extend the present work to nearby star forming regions for which the MSX resolution is adequate to probe the UIB emitting regions close to the exciting FUV source, and also covering a larger range of $`L_{FUV}`$.

In any case, our strong empirical correlation has many important and useful implications. For example, spatial distribution of emission in the UIBs for the entire Galactic plane can be studied following our method and the MSX survey with an angular resolution $``$ 20″. Of course the ISOCAM data provides the higher angular resolution information of selected regions and also the very important calibration factor above.

Since the band $`A`$ of MSX includes the UIB at 8.6 $`\mu `$m in addition to the ones at 6.2 and 7.7 $`\mu `$m, it is natural to obtain a value for $`_A`$ greater than unity. The mean value of $`_A`$ so obtained can be interpreted in terms of the relative strength of 8.6 $`\mu `$m feature vis a vis 6.2 $`\mu `$m + 7.7 $`\mu `$m features.

A similar correlation between $`UIB_C`$ (which is just a scaled down value of $`UIB_A`$) and the 11.2 $`\mu `$m feature emission measured using ISOCAM (LW8) by ZD has been explored for our sample of five sources. The ratio, $`_C`$ = I(UIB$`{}_{C}{}^{}{}_{}{}^{MSX}`$)/I(ISOCAM8), turns out to be 2.98 $`\pm `$ 0.28, which shows a much larger dispersion. In view of the fact that band $`C`$ of MSX includes the UIB feature at 12.7 $`\mu `$m also, the above is not surprizing. The average value of $`_C`$ should provide information regarding relative strengths of the features at 11.2 and 12.7 $`\mu `$m. The larger dispersion in the values of $`_C`$ among the five regions considered here, is perhaps indicative of the observed variability of the strength of 12.7 $`\mu `$m feature vis-a-vis other UIBs at shorter wavelengths (e.g. among the sources whose ISO-SWS spectra have been presented by Verstraete et al. 2001).

Next, we discuss, how sensitive is the value of $`_A`$ to the various assumptions made in analysing / modelling the MSX data. One of the parameters while modelling out the continuum from the MSX radiance in band $`A`$, is $`\beta `$, the dust emissivity index. Let us consider a few most popular types of dust grains, viz., Draine and Lee type (DL; Draine & Lee 1984), and Mathis, Mezger and Panagia (MMP; Mathis et al. 1983). One of the most popular dust size distribution is due to Mathis, Rumple and Nordsieck (MRN; Mathis et al. 1977). The MRN size distribution averaged values of absorption cross sections for DL and MMP type dust, for silicate and graphite grains have been computed earlier (Mookerjea & Ghosh 1999; Mookerjea et al. 1999). Using these for the MSX bands $`A`$ to $`E`$, we find the effective value of $`\beta `$ to be in the range 0.56 to 0.89. In any case, all calculations were repeated for the values of $`\beta `$ = 0 and 0.5 (in addition to 1.0). This had insignificant effect on the value of $`_A`$, and also the tightness of correlation between I(ISOCAM4+6) and I(UIB$`{}_{}{}^{MSX}{}_{A}{}^{}`$). It may be noted here that in case strong silicate absorption feature at $``$ 9.8 $`\mu `$m was important, the above correlation could have been lost. The sample of ZD had been selected such that silicate feature is not visible in the IRAS LRS spectrum.

Another very important parameter in our scheme is $`\alpha `$ (see equation 3), which has been held at a value of 0.32 on the basis of result of Verstraete et al. (2001). In order to study the sensitivity of our results on the numerical value of $`\alpha `$, we have repeated the calculations for a range of its values between 0 and 0.7. We find the following, for the value of $`\alpha `$ between 0.15 and 0.35, the correlation between I(ISOCAM4+6) and I(UIB$`{}_{}{}^{MSX}{}_{A}{}^{}`$) remains very tight, though the numerical value of $`_A`$ changes slightly between 2.0 and 2.3. For values of $`\alpha `$ outside this range (0.15–0.35), the correlation becomes much poorer and also the value of $`_A`$ decreases on $`either`$ side of this range ! All of the above can be understood, if the total UIB feature strength in MSX bands $`A`$ and $`C`$ are really proportional. Hence we conclude that value of $`_A`$ determined here is physically meaningful and should help quantifying the UIB emission in the Galactic plane in general.

## 5 Summary

A scheme has been developed to extract the contribution of Unidentified Infrared emission Bands (UIBs) from the mid infrared Galactic plane survey carried out by the SPIRIT III instrument onboard Midcourse Space Experiment (MSX) satellite in four bands. The scheme models the observations with a combination of thermal emission (gray body) from interstellar dust and the UIB emission from the gas component, under reasonable assumptions. Thus the spatial distribution of emission in the UIBs with an angular resolution $``$ 20″ (intrinsic to MSX survey) has been extracted.

In order to verify the reliability of this scheme, a detailed comparison has been made with the results obtained by Zavagno & Ducci (2001) using the ISOCAM instrument onboard Infrared Space Observatory (ISO), which has superior spectral and spatial resolutions than that of the MSX survey. Five Galactic star forming regions, viz., Sharpless 61 (Sh-61), Sh-138, Sh-152, Sh-156 and Sh-186, studied by Zavagno & Ducci (2001) have been used in this comparison.

The following results have been found :

* the UIB emission extracted from the MSX data is able to reproduce all major structural / morphological details detected by ISOCAM, consistent with its angular resolution;
* a very tight linear correlation has been found between the integrated UIB emission extracted from MSX data and that obtained from ISOCAM data (hence MSX based UIB radiance estimates can be calibrated well); The above validates our assumption regarding emissivity law for interstellar dust grains and is self consistent with the assumed proportionality of UIB emission between the bands $`A`$ and $`C`$ of MSX. The correlation may be a result of spatial averaging over regions with different excitation details.
* the numerical value of the ratio of UIB radiances obtained from MSX and ISOCAM is an important information for understanding relative strengths of UIBs.

Hence, we conclude that our empirical scheme along with the MSX Galactic plane survey can be a powerful tool to study large scale spatial distribution of UIB carriers.

###### Acknowledgements.

It is a pleasure to thank the referee Dr. M. P. Egan whose suggestions have improved the scientific content of this paper. This research made use of data products from the Midcourse Space Experiment. Processing of the data was funded by the Ballistic Missile Defense Organization with additional support from NASA Office of Space Science. This research has also made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The present work is based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA.


text:
# Towards Quantum Information Theory in Space and Time

## 1 Introduction

Remarkable experimental and theoretical results obtained in quantum computing, teleportation and cryptography (these topics sometimes are considered as belonging to quantum information theory) are based on the investigation of basic properties of quantum mechanics. Especially important are properties of nonfactorized entangled states introduced by Einstein, Podolsky and Rosen which were named by Schrodinger as the most characteristic feature of quantum mechanics.

Ideas of Shannon‘s classical information theory are important for the modern quantum information theory as well as the notions of qubit, quantum relative entropy, quantum channel, and entangled states , see for example  \- .

The spacetime dependence is not explicitly indicated in this approach. As a result, many important achievements in modern quantum information theory have been obtained for an idealized situation when the spacetime dependence of quantum phenomena is neglected.

We emphasize the importance of the investigation of quantum information effects in space and time. 111The importance of the investigation of quantum information effects in space and time and especially the role of relativistic invariance in classical and quantum information theory was stressed in the talk by the author at the First International Conference on Quantum Information which was held at Meijo University, Japan, November 4-8, 1997. Transmission and processing of (quantum) information is a physical process in spacetime. Therefore a formulation of such basic notions in quantum information theory as composite systems, entangled states and the channel should include the spacetime variables .

Ultimately, quantum information theory should become a part of quantum field theory (perhaps, in future, a part of superstring theory) since quantum field theory is our most fundamental physical theory.

Quantum field theory  is not just an abstract mathematical theory of operators in a Hilbert space. Basic equations of quantum field theory such as the Maxwell, Dirac, Yang–Mills equations are differential equations for operator functions defined on the spacetime. The nonrelativistic Schrodinger equation is also a differential equation in spacetime. Therefore a realistic quantum information theory should be based on the study of the solutions of these equations propagated in spacetime.

One could suggest to define a context described in  as a boundary condition for a differential equation. Then we would derive the contextual dependence of probabilities from the study of the dependence of solutions of the equation on the boundary conditions.

In modern quantum information theory the basic notion is the two dimentional Hilbert space, i.e. qubit. We suggest that in a relativistic quantum information theory, when the existense of spacetime is taken into account, the basic notion should be a notion of an elementary quantum system, i.e. according to Wigner (see ) it is an infinite dimensional Hilbert space $`H`$ invariant under an irreducible representation of the Poincare group labeled by $`[m,s]`$ where $`m0`$ is mass and $`s=0,1/2,1,\mathrm{}`$ is spin (helicity).

Entangled states, i.e. the states of two particles with the wave function which is not a product of the wave functions of single particles, have been studied in many theoretical and experimental works starting from the paper of Einstein, Podolsky and Rosen, see e.g. .

J. Bell proved  that there are quantum spin correlation functions in entangled states that can not be represented as classical correlation functions of separated random variables. Bell’s theorem reads, see :

$$\mathrm{cos}(\alpha \beta )E\xi _\alpha \eta _\beta $$

where $`\xi _\alpha `$ and $`\eta _\beta `$ are two random processes such that $`|\xi _\alpha |1`$$`|\eta _\beta |1`$ and $`E`$ is the expectation. Here the function $`\mathrm{cos}(\alpha \beta )`$ describes the quantum mechanical correlation of spins of two entangled particles. Bell‘s theorem has been interpreted as incompatibility of the requirement of locality with the statistical predictions of quantum mechanics . For a recent discussion of Bell’s theorem and Bell‘s inequalities see, for example  \-  and references therein.

However if we want to speak about locality in quantum theory then we have to localize somehow our particles. For example we could measure the density of the energy or the position of the particles simultaneously with the spin. Only then we could come to some conclusions about a relevance of the spin correlation function to the problem of locality.

The function $`\mathrm{cos}(\alpha \beta )`$ describes quantum correlations of two spins in the two qubit Hilbert space when the spacetime dependence of the wave functions of the particles is neglected. Let us note however that the very formulation of the problem of locality in quantum mechanics prescribes a special role to the position in ordinary three-dimensional space. It is rather strange therefore that the problem of local in space observations was neglected in discussions of the problem of locality in relation to Bell’s inequalities .

Let us stress that we discuss here not a problem of interpretation of quantum theory but a problem of how to make correct quantum mechanical computations describing an experiment with two detectors localized in space. Recently it was pointed out  that if we make local observations of spins then the spacetime part of the wave function leads to an extra factor in quantum correlations and as a result the ordinary conclusion from the Bell theorem about the nonlocality of quantum theory fails.

We present a modification of Bell‘s equation which includes space and time variables. The function $`\mathrm{cos}(\alpha \beta )`$ describes the quantum mechanical correlation of spins of two entangled particles if we neglect the spacetime dependence of the wave function. It was shown in  that if one takes into account the space part of the wave function then the quantum correlation describing local observations of spins in the simplest case will take the form $`g\mathrm{cos}(\alpha \beta )`$ instead of just $`\mathrm{cos}(\alpha \beta )`$. Here the parameter $`g`$ describes the location of the system in space and time. In this case one gets a modified equation

$$g\mathrm{cos}(\alpha \beta )=E\xi _\alpha \eta _\beta $$

One can prove that if the distance between detectors is large enough then the factor $`g`$ becomes small and there exists a solution of the modified equation. We will show that in fact at large distances all reasonable quantum states become disentangled. This fact leads also to important consequences for quantum teleportation and quantum cryptography, . Bell’s theorem constitutes an important part in quantum cryptography. In  it is discussed how one can try to improve the security of quantum cryptography schemes in space by using a special preparation of the space part of the wave function.

It is important to study also a more general question: which class of functions $`f(s,t)`$ admits a representation of the form

$$f(s,t)=Ex_sy_t$$

where $`x_s`$ and $`y_t`$ are bounded stochastic processes and also analogous question for the functions of several variables $`f(t_1,\mathrm{},t_n).`$

Such considerations could provide a noncommutative generalization of von Neumann‘s spectral theorem.

In this paper entangled states in space and time are considered. We point out a simple but the fundamental fact that the vacuum state $`\omega _0`$ in a free quantum field theory is a nonfactorized (entangled) state for observables belonging to spacelike separated regions:

$$\omega _0(\phi (x)\phi (y))\omega _0(\phi (x))\omega _0(\phi (y))0$$

Here $`\phi (x)`$ is a free scalar field in the Minkowski spacetime and $`(xy)^2<0`$. Hence there is a statistical dependence between causally disconnected regions.

However one has an asymptotic factorization of the vacuum state for large separations of the spacelike regions. Moreover one proves that in quantum field theory there is an asymptotic factorization for any reasonable state and any local observables. Therefore at large distances any reasonable state becomes disentangled. We have the relation

$$\underset{|l|\mathrm{}}{lim}[\omega (A(l)B)\omega (A(l))\omega (B)]=0$$

Here $`\omega `$ is a state from a rather wide class of the states which includes entangled states, $`A`$ and $`B`$ are two local observables, and $`A(l)`$ is the translation of the observable $`A`$ along the 3 dim vector $`l`$. As a result a violation of Bell‘s inequalities (see below) can be observed without inconsistency with principles of relativistic quantum theory only if the distance between detectors is rather small. We suggest a further experimental study of entangled states in spacetime by studying the dependence of the correlation functions on the distance between detectors.

There is no a factorization of the expectation value $`\omega _0(\phi (x)\phi (y))`$ even for the space-like separation of the variables $`x`$ and $`y`$ if the distance between $`x`$ and $`y`$ is not large enough. However we will prove that there exist a representation of the form

$$\omega _0(\phi (x)\phi (y))=E\xi (x)\xi ^{}(y)$$

which is valid for all $`x`$ and $`y`$. Here $`\xi (x)`$ is a classical (generalized) complex random field and $`E`$ is the expectation value. Therefore the quantum correlation function is represented as a classical correlation function of separated random fields. This representation can be called a local realistic representation by analogy with the Bell approach to the spin correlation functions.

In the next section Bell‘s theorem is discussed and a slight generalization of the known CHSH result is proved. In Sect.3 the locality in space is considered for entangled states and the asymptotic factorization of the states is proved. A hidden variable representation for quantum correlation which is local in the space is also obtained. Noncommutative spectral theory and local realism are considered in Sect.4. Finally the disentanglement at large distances in quantum field theory is considered in Sect.5.

## 2 Bell’s Theorem

### 2.1 Bell‘s Theorem and Stochastic Processes

In the presentation of Bell’s theorem we will follow  where one can find also more references. Bell’s theorem reads:

$$\mathrm{cos}(\alpha \beta )E\xi _\alpha \eta _\beta $$
(1)
where $`\xi _\alpha `$ and $`\eta _\beta `$ are two random processes such that $`|\xi _\alpha |1`$$`|\eta _\beta |1`$ and $`E`$ is the expectation. In more details:

Theorem 1. There exists no probability space $`(\mathrm{\Lambda },,d\rho (\lambda ))`$ and a pair of stochastic processes $`\xi _\alpha =\xi _\alpha (\lambda ),\eta _\beta =\eta _\beta (\lambda ),0\alpha ,\beta 2\pi `$ which obey $`|\xi _\alpha (\lambda )|1`$$`|\eta _\beta (\lambda )|1`$ such that the following equation is valid

$$\mathrm{cos}(\alpha \beta )=E\xi _\alpha \eta _\beta $$
(2)
for all $`\alpha `$ and $`\beta `$.

Here $`\mathrm{\Lambda }`$ is a set, $``$ is a sigma-algebra of subsets and $`d\rho (\lambda )`$ is a probability measure, i.e. $`d\rho (\lambda )0,𝑑\rho (\lambda )=1.`$ The expectation is

$$E\xi _\alpha \eta _\beta =_\mathrm{\Lambda }\xi _\alpha (\lambda )\eta _\beta (\lambda )𝑑\rho (\lambda )$$

One can write Eq. (2) as an integral equation

$$\mathrm{cos}(\alpha \beta )=_\mathrm{\Lambda }\xi _\alpha (\lambda )\eta _\beta (\lambda )𝑑\rho (\lambda )$$
(3)
We say that the integral equation (3) has no solutions $`(\mathrm{\Lambda },,d\rho (\lambda ),\xi _\alpha ,\eta _\beta )`$ with the bound $`|\xi _\alpha |1`$$`|\eta _\beta |1.`$

We will prove the theorem below. Let us discuss now the physical interpretation of this result.

Consider a pair of spin one-half particles formed in the singlet spin state and moving freely towards two detectors. If one neglects the space part of the wave function then one has the Hilbert space $`C^2C^2`$ and the quantum mechanical correlation of two spins in the singlet state $`\psi _{spin}C^2C^2`$ is

$$D_{spin}(a,b)=\psi _{spin}|\sigma a\sigma b|\psi _{spin}=ab$$
(4)
Here $`a=(a_1,a_2,a_3)`$ and $`b=(b_1,b_2,b_3)`$ are two unit vectors in three-dimensional space $`R^3`$, $`\sigma =(\sigma _1,\sigma _2,\sigma _3)`$ are the Pauli matrices,

$$\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\sigma a=\underset{i=1}{\overset{3}{}}\sigma _ia_i$$

and

$$\psi _{spin}=\frac{1}{\sqrt{2}}\left(\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)\right)$$

If the vectors $`a`$ and $`b`$ belong to the same plane then one can write $`ab=\mathrm{cos}(\alpha \beta )`$ and hence Bell’s theorem states that the function $`D_{spin}(a,b)`$ Eq. (4) can not be represented in the form

$$P(a,b)=\xi (a,\lambda )\eta (b,\lambda )𝑑\rho (\lambda )$$
(5)
i.e.

$$D_{spin}(a,b)P(a,b)$$
(6)
Here $`\xi (a,\lambda )`$ and $`\eta (b,\lambda )`$ are random fields on the sphere, $`|\xi (a,\lambda )|1`$$`|\eta (b,\lambda )|1`$ and $`d\rho (\lambda )`$ is a positive probability measure, $`𝑑\rho (\lambda )=1`$. The parameters $`\lambda `$ are interpreted as hidden variables in a realist theory. It is clear that Eq. (6) can be reduced to Eq. (1).

### 2.2 CHSH Inequality

To prove Theorem 1 we will use the following theorem which is a slightly generalized Clauser-Horn-Shimony-Holt (CHSH) result.

Theorem 2. Let $`f_1,f_2,g_1`$ and $`g_2`$ be random variables (i.e. measured functions) on the probability space $`(\mathrm{\Lambda },,d\rho (\lambda ))`$ such that

$$|f_i(\lambda )g_j(\lambda )|1,i,j=1,2.$$

Denote

$$P_{ij}=Ef_ig_j,i,j=1,2.$$

Then

$$|P_{11}P_{12}|+|P_{21}+P_{22}|2.$$

Proof of Theorem 2. One has

$$P_{11}P_{12}=Ef_1g_1Ef_1g_2=E(f_1g_1(1\pm f_2g_2))E(f_1g_2(1\pm f_2g_1))$$

Hence

$$|P_{11}P_{12}|E(1\pm f_2g_2)+E(1\pm f_2g_1)=2\pm (P_{22}+P_{21})$$

Now let us note that if $`x`$ and $`y`$ are two real numbers then

$$|x|2\pm y|x|+|y|2.$$

Therefore taking $`x=P_{11}P_{12}`$ and $`y=P_{22}+P_{21}`$ one gets the bound

$$|P_{11}P_{12}|+|P_{21}+P_{22}|2.$$

The theorem is proved.

The last inequality is called the CHSH inequality. By using notations of Eq. (5) one has

$$|P(a,b)P(a,b^{})|+|P(a^{},b)+P(a^{},b^{})|2$$
(7)
for any four unit vectors $`a,b,a^{},b^{}`$.

Proof of Theorem 1. Let us denote

$$f_i(\lambda )=\xi _{\alpha _i}(\lambda ),g_j(\lambda )=\eta _{\beta _j}(\lambda ),i,j=1,2$$

for some $`\alpha _i,\beta _j.`$ If one would have

$$\mathrm{cos}(\alpha _i\beta _j)=Ef_ig_j$$

then due to Theorem 2 one should have

$$|\mathrm{cos}(\alpha _1\beta _1)\mathrm{cos}(\alpha _1\beta _2)|+|\mathrm{cos}(\alpha _2\beta _1)+\mathrm{cos}(\alpha _2\beta _2)|2.$$

However for $`\alpha _1=\pi /2,\alpha _2=0,\beta _1=\pi /4,\beta _2=\pi /4`$ we obtain

$$|\mathrm{cos}(\alpha _1\beta _1)\mathrm{cos}(\alpha _1\beta _2)|+|\mathrm{cos}(\alpha _2\beta _1)+\mathrm{cos}(\alpha _2\beta _2)|=2\sqrt{2}$$

which is greater than 2. This contradiction proves Theorem 1.

It will be shown below that if one takes into account the space part of the wave function then the quantum correlation in the simplest case will take the form $`g\mathrm{cos}(\alpha \beta )`$ instead of just $`\mathrm{cos}(\alpha \beta )`$ where the parameter $`g`$ describes the location of the system in space and time. In this case one can get a representation

$$g\mathrm{cos}(\alpha \beta )=E\xi _\alpha \eta _\beta $$
(8)
if $`g`$ is small enough. The factor $`g`$ gives a contribution to visibility or efficiency of detectors that are used in the phenomenological description of detectors.

## 3 Local Observations

### 3.1 Modified Bell‘s equation

In the previous section the space part of the wave function of the particles was neglected. However exactly the space part is relevant to the discussion of locality. The Hilbert space assigned to one particle with spin 1/2 is $`C^2L^2(R^3)`$ and the Hilbert space of two particles is $`C^2L^2(R^3)C^2L^2(R^3).`$ The complete wave function is $`\psi =(\psi _{\alpha \beta }(𝐫_1,𝐫_2,t))`$ where $`\alpha `$ and $`\beta `$ are spinor indices, $`t`$ is time and $`𝐫_1`$ and $`𝐫_2`$ are vectors in three-dimensional space.

We suppose that there are two detectors (A and B) which are located in space $`R^3`$ within the two localized regions $`𝒪_A`$ and $`𝒪_B`$ respectively, well separated from one another. If one makes a local observation in the region $`𝒪_A`$ then this means that one measures not only the spin observable $`\sigma _i`$ but also some another observable which describes the localization of the particle like the energy density or the projection operator $`P_𝒪`$ to the region $`𝒪`$. We will consider here corelation functions of the projection operators $`P_𝒪`$.

Quantum correlation describing the localized measurements of spins in the regions $`𝒪_A`$ and $`𝒪_B`$ is

$$\omega (\sigma aP_{𝒪_A}\sigma bP_{𝒪_B})=\psi |\sigma aP_{𝒪_A}\sigma bP_{𝒪_B}|\psi $$
(9)
Let us consider the simplest case when the wave function has the form of the product of the spin function and the space function $`\psi =\psi _{spin}\varphi (𝐫_1,𝐫_2)`$. Then one has

$$\omega (\sigma aP_{𝒪_A}\sigma bP_{𝒪_B})==g(𝒪_A,𝒪_B)D_{spin}(a,b)$$
(10)
where the function

$$g(𝒪_A,𝒪_B)=_{𝒪_A\times 𝒪_B}|\varphi (𝐫_1,𝐫_2)|^2𝑑𝐫_1𝑑𝐫_2$$
(11)
describes correlation of particles in space. It is the probability to find one particle in the region $`𝒪_A`$ and another particle in the region $`𝒪_B`$.

One has

$$0g(𝒪_A,𝒪_B)1.$$
(12)
If $`𝒪_A`$ is a bounded region and $`𝒪_A(l)`$ is a translation of $`𝒪_A`$ to the 3-vector $`l`$ then one has

$$\underset{|l|\mathrm{}}{lim}g(𝒪_A(l),𝒪_B)=0.$$
(13)
Since

$$\psi _{spin}|\sigma aI|\psi _{spin}=0$$

we have

$$\omega (\sigma aP_{𝒪_A}I)=0.$$

Therefore we have proved the following proposition which says that the state $`\psi =\psi _{spin}\varphi (𝐫_1,𝐫_2)`$ becomes disentangled at large distances.

Proposition. One has the following property of the asymptotic factorization (disentanglement) at large distances:

$$\underset{|l|\mathrm{}}{lim}[\omega (\sigma aP_{𝒪_A(l)}\sigma bP_{𝒪_B})\omega (\sigma aP_{𝒪_A(l)}I)\omega (I\sigma bP_{𝒪_B})]=0$$
(14)
or

$$\underset{|l|\mathrm{}}{lim}\omega (\sigma aP_{𝒪_A(l)}\sigma bP_{𝒪_B})=0.$$

Now one inquires whether one can write a representation

$$\omega (\sigma aP_{𝒪_A(l)}\sigma bP_{𝒪_B})=\xi (a,𝒪_A,\lambda )\eta (b,𝒪_B,\lambda )𝑑\rho (\lambda )$$
(15)
where $`|\xi (a,𝒪_A(l),\lambda )|1,|\eta (b,𝒪_B,\lambda )|1`$.

Remark. A local modified equation reads

$$|\varphi (𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐},t)|^2\mathrm{cos}(\alpha \beta )=E\xi (\alpha ,𝐫_\mathrm{𝟏},t)\eta (\beta ,𝐫_\mathrm{𝟐},t).$$

If we are interested in the conditional probability of finding the projection of spin along vector $`a`$ for the particle 1 in the region $`𝒪_A(l)`$ and the projection of spin along the vector $`b`$ for the particle 2 in the region $`𝒪_B`$ then we have to divide both sides of Eq. (15) by $`g(𝒪_A(l),𝒪_B)`$.

Note that here the classical random variable $`\xi =\xi (a,𝒪_A(l),\lambda )`$ is not only separated in the sense of Bell (i.e. it depends only on $`a`$) but it is also local in the 3 dim space since it depends only on the region $`𝒪_A(l)`$. The classical random variable $`\eta `$ is also local in 3 dim space since it depends only on $`𝒪_B`$. Note also that since the eigenvalues of the projector $`P_𝒪`$ are 0 or 1 then one should have $`|\xi (a,𝒪_A)|1.`$

Due to the property of the asymptotic factorization and the vanishing of the quantum correlation for large $`|l|`$ there exists a trivial asymptotic classical representation of the form (15) with $`\xi =\eta =0.`$

We can do even better and find a classical representation which will be valid uniformly for large $`|l|`$.

If $`g`$ would not depend on $`𝒪_A`$ and $`𝒪_B`$ then instead of Eq (2) in Theorem 1 we could have a modified equation

$$g\mathrm{cos}(\alpha \beta )=E\xi _\alpha \eta _\beta $$
(16)
The factor $`g`$ is important. In particular one can write the following representation  for $`0g1/2`$:

$$g\mathrm{cos}(\alpha \beta )=_0^{2\pi }\sqrt{2g}\mathrm{cos}(\alpha \lambda )\sqrt{2g}\mathrm{cos}(\beta \lambda )\frac{d\lambda }{2\pi }$$
(17)
Therefore if $`0g1/2`$ then there exists a solution of Eq. (16) where

$$\xi _\alpha (\lambda )=\sqrt{2g}\mathrm{cos}(\alpha \lambda ),\eta _\beta (\lambda )=\sqrt{2g}\mathrm{cos}(\beta \lambda )$$

and $`|\xi _\alpha |1,|\eta _\beta |1.`$ If $`g>1/\sqrt{2}`$ then it follows from Theorem 2 that there is no solution to Eq. (16). We have obtained

Theorem 3. If $`g>1/\sqrt{2}`$ then there is no solution $`(\mathrm{\Lambda },,d\rho (\lambda ),\xi _\alpha ,\eta _\beta )`$ to Eq. (16) with the bounds $`|\xi _\alpha |1`$$`|\eta _\beta |1.`$ If $`0g1/2`$ then there exists a solution to Eq. (16) with the bounds $`|\xi _\alpha |1`$$`|\eta _\beta |1.`$

Remark. Further results on solutions of the modified equation have been obtained by A.K. Guschchin, S. V. Bochkarev and D. Prokhorenko. Local variable models for inefficient detectors are presented in .

Let us take now the wave function $`\varphi `$ of the form $`\varphi =\psi _1(𝐫_1)\psi _2(𝐫_2)`$ where

$$_{R^3}|\psi _1(𝐫_1)|^2𝑑𝐫_1=1,_{R^3}|\psi _2(𝐫_2)|^2𝑑𝐫_2=1$$

In this case

$$g(𝒪_A(l),𝒪_B)=_{𝒪_A(l)}|\psi _1(𝐫_1)|^2𝑑𝐫_1_{𝒪_B}|\psi _2(𝐫_2)|^2𝑑𝐫_2$$

There exists such $`L>0`$ that

$$_{B_L}|\psi _1(𝐫_1)|^2𝑑𝐫_1=ϵ<1/2,$$

where $`B_L=\{𝐫R^3:|𝐫|L\}.`$ Let us make an additional assumption that the classical random variable has the form of a product of two independent classical random variables $`\xi (a,𝒪_A)=\xi _{space}(𝒪_A)\xi _{spin}(a)`$ and similarly for $`\eta .`$ We have the following

Theorem 4. Under the above assumptions and for large enough $`|l|`$ there exists the following representation of the quantum correlation function

$$g(𝒪_A(l),𝒪_B)\mathrm{cos}(\alpha \beta )=(E\xi _{space}(𝒪_A)(l))(E\eta _{space}(𝒪_B))E\xi _{spin}(\alpha )\xi _{spin}(\beta )$$

where all classical random variables are bounded by 1.

Proof. To prove the theorem we write

$$g(𝒪_A(l),𝒪_B)\mathrm{cos}(\alpha \beta )=_{𝒪_A(l)}\frac{1}{ϵ}|\psi _1(𝐫_1)|^2𝑑𝐫_1_{𝒪_B}|\psi _2(𝐫_2)|^2𝑑𝐫_2ϵ\mathrm{cos}(\alpha \beta )$$

$$=(E\xi _{space}(𝒪_A(l))(E\eta _{space}(𝒪_B))E\xi _{spin}(\alpha )\xi _{spin}(\beta )$$

Here $`\xi _{space}(𝒪_A(l))`$ and $`\eta _{space}(𝒪_B)`$ are random variables on the probability space $`B_L\times R^3`$ with the probability measure

$$dP(𝐫_1,𝐫_2)=\frac{1}{ϵ}|\psi _1(𝐫_1)|^2|\psi _2(𝐫_2)|^2d𝐫_1d𝐫_2$$

of the form

$$\xi _{space}(𝒪_A(l),𝐫_1,𝐫_2)=\chi _{𝒪_A(l)}(𝐫_1),\eta _{space}(𝒪_B,𝐫_1,𝐫_2)=\chi _{𝒪_B}(𝐫_2)$$

where $`\chi _𝒪(𝐫)`$ is the characteristic function of the region $`𝒪.`$ We assume that $`𝒪_A(l)`$ belongs to $`B_L.`$ Further $`\xi _{spin}(\alpha )`$ is a random process on the circle $`0\phi 2\pi `$ with the probability measure $`d\phi /2\pi `$ of the form

$$\xi _{spin}(\alpha ,\phi )=\sqrt{2ϵ}\mathrm{cos}(\alpha \phi )$$

The theorem is proved.

### 3.2 Expansion of Wave Packet

Let us remind that there is a well known effect of expansion of wave packets due to the free time evolution. If $`ϵ`$ is the characteristic length of the Gaussian wave packet describing a particle of mass $`M`$ at time $`t=0`$ then at time $`t`$ the characteristic length $`ϵ_t`$ will be

$$ϵ_t=ϵ\sqrt{1+\frac{\mathrm{}^2t^2}{M^2ϵ^4}}.$$
(18)
It tends to $`(\mathrm{}/Mϵ)t`$ as $`t\mathrm{}`$. Therefore the locality criterion is always satisfied for nonrelativistic particles if regions $`𝒪_A`$ and $`𝒪_B`$ are far enough from each other.

### 3.3 Relativistic Particles

We can not immediately apply the previous considerations to the case of relativistic particles such as photons and the Dirac particles because in these cases the wave function can not be represented as a product of the spin part and the spacetime part. Let us show that the wave function of photon can not be represented in the product form. Let $`A_i(k)`$ be the wave function of photon, where $`i=1,2,3`$ and $`kR^3`$. One has the gauge condition $`k^iA_i(k)=0`$ . If one supposes that the wave function has a product form $`A_i(k)=\varphi _if(k)`$ then from the gauge condition one gets $`A_i(k)=0`$. Therefore the case of relativistic particles requires a separate investigation (see below).

## 4 Noncommutative Spectral Theory and Local Realism

As a generalisation of the previous discussion we would like to suggest here a general relation between quantum theory and theory of classical stochastic processes which expresses the condition of local realism. Let $``$ be a Hilbert space, $`\rho `$ is the density operator, $`\{A_\alpha \}`$ is a family of self-adjoint operators in $``$. One says that the family of observables $`\{A_\alpha \}`$ and the state $`\rho `$ satisfy to the condition of local realism if there exists a probability space $`(\mathrm{\Lambda },,d\rho (\lambda ))`$ and a family of random variables $`\{\xi _\alpha \}`$ such that the range of $`\xi _\alpha `$ belongs to the spectrum of $`A_\alpha `$ and for any subset $`\{A_i\}`$ of mutually commutative operators one has a representation

$$Tr(\rho A_{i_1}\mathrm{}A_{i_n})=E\xi _{i_1}\mathrm{}\xi _{i_n}$$

The physical meaning of the representation is that it describes the quantum-classical correspondence. If the family $`\{A_\alpha \}`$ would be a maximal commutative family of self-adjoint operators then for pure states the previous representation can be reduced to the von Neumann spectral theorem . In our case the family $`\{A_\alpha \}`$ consists from not necessary commuting operators. Hence we will call such a representation a noncommutative spectral representation. Of course one has a question for which families of operators and states a noncommutative spectral theorem is valid, i.e. when we can write the noncommutative spectral representation. We need a noncommutative generalization of von Neumann‘s spectral theorem.

It would be helpful to study the following problem: describe the class of functions $`f(t_1,\mathrm{},t_n)`$ which admits the representation of the form

$$f(t_1,\mathrm{},t_n)=Ex_{t_1}\mathrm{}z_{t_n}$$

where $`x_t,\mathrm{},z_t`$ are random processes which obey the bounds $`|x_t|1,\mathrm{},|z_t|1`$.

From the previous discussion we know that there are such families of operators and such states which do not admit the noncommutative spectral representation and therefore they do not satisfy the condition of local realism. Indeed let us take the Hilbert space $`=C^2C^2`$ and four operators $`A_1,A_2,A_3,A_4`$ of the form (we denote $`A_3=B_1,A_4=B_2`$)

$$A_1=\left(\begin{array}{cc}\mathrm{sin}\alpha _1& \mathrm{cos}\alpha _1\\ \mathrm{cos}\alpha _1& \mathrm{sin}\alpha _1\end{array}\right)I,A_2=\left(\begin{array}{cc}\mathrm{sin}\alpha _2& \mathrm{cos}\alpha _2\\ \mathrm{cos}\alpha _2& \mathrm{sin}\alpha _2\end{array}\right)I$$

and

$$B_1=I\left(\begin{array}{cc}\mathrm{sin}\beta _1& \mathrm{cos}\beta _1\\ \mathrm{cos}\beta _1& \mathrm{sin}\beta _1\end{array}\right),B_2=I\left(\begin{array}{cc}\mathrm{sin}\beta _2& \mathrm{cos}\beta _2\\ \mathrm{cos}\beta _2& \mathrm{sin}\beta _2\end{array}\right)$$

Here operators $`A_i`$ correspond to operators $`\sigma a`$ and operators $`B_i`$ corresponds to operators $`\sigma b`$ where $`a=(\mathrm{cos}\alpha ,0,\mathrm{sin}\alpha ),b=(\mathrm{cos}\beta ,0,\mathrm{sin}\beta ).`$ Operators $`A_i`$ commute with operators $`B_j,[A_i,B_j]=0,i,j=1,2`$ and one has

$$\psi _{spin}|A_iB_j|\psi _{spin}=\mathrm{cos}(\alpha _i\beta _j),i,j=1,2$$

We know from Theorem 2 that this function can not be represented as the expected value $`E\xi _i\eta _j`$ of random variables with the bounds $`|\xi _i|1`$$`|\eta _j|1.`$

However, as it was discussed above, the space part of the wave function was neglected in the previous consideration. We suggest that in physics one could prepare only such states and observables which satisfy the condition of local realism. Perhaps we should restrict ourself in this proposal to the consideration of only such families of observables which satisfy the condition of relativistic local causality. If there are physical phenomena which do not satisfy this proposal then it would be important to describe quantum processes which satisfy the above formulated condition of local realism and also processes which do not satisfy this condition.

## 5 Quantum Probability and Quantum Field Theory

In quantum probability (see ) we are given a * - algebra $`𝒜`$ and a state (i.e. a linear positive normalized functional) $`\omega `$ on $`𝒜`$. Elements from $`𝒜`$ are called random variables. Two random variables $`A`$ and $`B`$ are called (statistically) independent if $`\omega (AB)=\omega (A)\omega (B)`$.

First we will prove the following

Proposition. There is a statistical dependence between two spacelike separated regions in the theory of free scalar quantum field.

Proof. Let us consider a free scalar quantum field $`\phi (x)`$:

$$\phi (x)=\frac{1}{(2\pi )^{3/2}}_{R^3}\frac{d𝐤}{\sqrt{2k^0}}(e^{ikx}a^{}(𝐤)+e^{ikx}a(𝐤))$$

Here $`kx=k^0x^0\mathrm{𝐤𝐱}`$, $`k^0=\sqrt{𝐤^2+m^2},m0`$ and $`a(𝐤)`$ and $`a^{}(𝐤)`$ are annihilation and creation operators,

$$[a(𝐤),a^{}(𝐤^{})]=\delta (𝐤𝐤^{})$$

The field $`\phi (x)`$ is an operator valued distribution acting in the Fock space $``$ with the vacuum $`|0>`$,

$$a(𝐤)|0>=0$$

The vacuum expectation value of two fields is

$$\omega _0(\phi (x)\phi (y))=<0|\phi (x)\phi (y)|0>=W_0(xy,m^2)$$

where

$$W_0(xy,m^2)=\frac{1}{(2\pi )^3}_{R^3}\frac{d𝐤}{2k^0}e^{ik(xy)}$$

The statistical independence of two spacelike separated regions in particular would lead to the relation

$$\omega _0(\phi (x)\phi (y))\omega _0(\phi (x))\omega _0(\phi (y))=0$$

if $`(xy)^2<0`$. But since $`\omega _0(\phi (x))=0`$ in fact we have

$$\omega _0(\phi (x)\phi (y))\omega _0(\phi (x))\omega _0(\phi (y))=W_0(xy,m^2)0$$

However the violation of the statistical independence vanish exponentially with the spacial separation of $`x`$ and $`y`$ since for large $`\lambda =m\sqrt{x^2}`$ the function $`W_0(x,m^2)`$ behaves like

$$\frac{m^2}{4\pi \lambda }\left(\frac{\pi }{2\lambda }\right)^{1/2}e^\lambda $$

Let us prove that any polynomial state is asymptotically disentangled (factorized) for large spacelike distances. Let $`𝒜`$ be the algebra of polinomials in the Fock space $``$ at the field $`\phi (f)`$ with the test functions $`f`$ . Let $`C𝒜`$ and $`|\psi >=C|0>.`$ Denote the state $`\omega (A)=\psi |A|\psi /\psi ^2`$ for $`A𝒜`$.

Theorem 5. One has the following asymptotic disentanglement property

$$\underset{|l|\mathrm{}}{lim}[\omega (A(l)B)\omega (A(l))\omega (B)]=0$$

Here $`A`$ and $`B`$ belong to $`𝒜`$ and $`A(l)`$ is the translation of $`A`$ along the 3 dim vector $`l`$. One has also

$$\underset{|l|\mathrm{}}{lim}[\omega (A(l))0|A(l)|0]=0$$

The proof of the theorem is based on the Wick theorem and the Riemann-Lebesgue lemma.

Similar theorems take place also for the Dirac and the Maxwell fields. In particular for the Dirac field $`\psi (x)`$ one can prove the asymtotic factorization for the local spin operator

$$𝐒(𝒪)=_𝒪\psi ^{}𝚺\psi 𝑑x$$

Here $`\mathrm{\Sigma }`$ is made from the Dirac matrices.

Finally let us show that some correlation functions in the relativistic quantum field theory can be represented as mathematical expectations of the classical (generalized) random fields.

Theorem 6. If $`\phi (x)`$ is a scalar complex quantum field then one has a representation

$$0|\phi (x_1)\mathrm{}\phi (x_n)\phi ^{}(y_1)\mathrm{}\phi ^{}(y_n)|0=E\xi (x_1)\mathrm{}\xi (x_n)\xi ^{}(y_1)\mathrm{}\xi ^{}(y_n).$$

Here $`\xi (x)`$ is a complex random field.

The proof of the theorem follows from the positivity of the quantum correlation functions. It is interesting that we have obtained a functional integral representation for the quantum correlation functions in real time. Similar representation is valid also for the 2-point correlation function of an interacting scalar field. It follows from the Kallen-Lehmann representation.

## 6 Conclusions

We have discussed some problems in quantum information theory which requires the inclusion of spacetime variables. In particular entangled states in space and time were considered. A modification of Bell‘s equation which includes the spacetime variables is suggested and investigated. A general relation between quantum theory and theory of classical stochastic processes was proposed which expresses the condition of local realism in the form of a noncommutative spectral theorem. Applications of this relation to the security of quantum key distribution in quantum cryptography was mentioned.

There are many interesting open problems in the approach to quantum information in space and time discussed in this paper. Some of them related with the noncommutative spectral theory and theory of classical stochastic processes have been discussed above. It would be useful if the local algebraic approach to quantum theory  will be developed in this direction.

Entangled states in space and time are considered. It is shown that any reasonable state in relativistic quantum field theory becomes disentangled (factorizable) at large spacelike distances if one makes local observations. As a result a violation of Bell‘s inequalities can be observed without inconsistency with principles of relativistic quantum theory only if the distance between detectors is rather small. We suggest a further experimental study of entangled states in spacetime by studying the dependence of the correlation functions on the distance between detectors.

## Acknowledgments

I am grateful to R. Gill, A. Holevo, A. Khrennikov, T.Ishiwatari, S. Lloyd, S. Molotkov, M. Ohya, R. Roschin, A. Suarez and K.A. Valiev for helpful discussions and remarks. This work is supported in part by RFFI and by INTAS 99-00545 grants.