text:
# Chaos in PDEs and Lax Pairs of Euler Equations

## 1. Introduction

Unlike chaos in finite dimensional systems, the area of chaos in partial differential equations has been a long standing open field. During the past ten years, there have been some remarkable developments in the area of chaos in partial differential equations by the author and collaborators. Our main results can be briefly described as follows: A systematic program is developed for proving the existence of homoclinic orbits for perturbed soliton equations. The program involves machineries from integrable theory, dynamical systems, and partial differential equations. The types of homoclinic orbits thus obtained can be either transversal or of Silnikov. With regard to transversal homoclinic orbits, shadowing lemmas are utilized or developed to prove the existence of chaos. With regard to Silnikov homoclinic orbits, Smale horseshoes are constructed to prove the existence of chaos. The main machineries from integrable theory are the isospectral theory and Darboux transformations. The main machineries from dynamical systems are persistence of invariant manifolds, Fenichel fibers, Melnikov measurement, and other measurements. The main machineries from partial differential equations are local well-posedness and regularity.

Another exciting development happened to Euler equations of incompressible inviscid fluids. Lax pairs for both 2D and 3D Euler equations was found by the author and Steve Childress. A Darboux transformation for the Lax pair of 2D Euler equation was also found by Artyom Yurov and the author. The philosophical significance of the existence of Lax pairs for Euler equations is fundamentally important. If one defines integrability of an equation by the existence of a Lax pair, then both 2D and 3D Euler equations are integrable. More importantly, both 2D and 3D Navier-Stokes equations at high Reynolds numbers are near integrable systems. Such a point of view changes our old ideology on Euler and Navier-Stokes equations.

The most important application of the theory on chaos in partial differential equations in theoretical physics will be on the study of turbulence. Our approach of studying turbulence is different from many other studies in which one starts with Stokes equation to prove results on Navier-Stokes equations for small Reynolds number. In our studies, we start with Euler equations and view Navier-Stokes equations at large Reynolds number as singular perturbations of Euler equations. For this goal, we chose the 2D Navier-Stokes equations under periodic boundary conditions to begin a dynamical system study on 2D turbulence. We started with a simple fixed point of 2D Euler equation, and studied the linearized 2D Euler equation at the fixed point. A complete spectral theorem is obtained. In particular, we found unstable eigenvalues. Then, naturally we are interested in the corresponding unstable manifold. Such an unstable manifold forms the prototype of the attractor of Navier-Stokes equations at large Reynolds number. And long term turbulence lives around the attractor. Finding the unstable manifold is not successful yet. Partial success is accomplished in finding the corresponding unstable manifold for some Galerkin truncation. Here the unstable manifold is the 2D surface of an ellipsoid. Stable and unstable manifolds together form a lip-shape configuration. Whether or not the stable and unstable manifolds for 2D Euler equation have the same topology as above is not known. In fact, the existence of stable and unstable manifolds for 2D Euler equation has not been proved.

## 2. Homoclinic Orbits

In terms of proving the existence of a homoclinic orbit, the most common tool is the so-called Melnikov integral method  . This method was subsequently developed by Holmes and Marsden , and most recently by Wiggins . For partial differential equations, this method was mainly developed by Li et al.    .

There are two derivations for the Melnikov integrals. One is the so-called geometric argument     . The other is the so-called Liapunov-Schmitt argument  . The Liapunov-Schmitt argument is a fixed-point type argument which directly leads to the existence of a homoclinic orbit. The condition for the existence of a fixed point is the Melnikov integral. The geometric argument is a signed distance argument which applies to more general situations than the Liapunov-Schmitt argument. It turns out that the geometric argument is a much more powerful machinary than the Liapunov-Schmitt argument. In particular, the geometric argument can handle geometric singular perturbation problems. I shall also mention an interesting derivation in .

In establishing the existence of homoclinic orbits in high dimensions, one often needs other tools besides the Melnikov analysis. For example, when studying orbits homoclinic to fixed points created through resonances in ($`n4`$)-dimensional near-integrable systems, one often needs tools like Fenichel fibers, as presented in previous chapter, to set up geometric measurements for locating such homoclinic orbits. Such homoclinic orbits often have a geometric singular perturbation nature. In such cases, the Liapunov-Schmitt argument can not be applied. For such works on finite dimensional systems, see for example   . For such works on infinite dimensional systems, see for example  .

### 2.1. Silnikov Homoclinic Orbits in Nonlinear Schrödinger (NLS) Equation Under Regular Perturbations

Consider the regularly perturbed nonlinear Schrödinger (NLS) equation ,

(2.1)
$$iq_t=q_{xx}+2[|q|^2\omega ^2]q+i\epsilon [\widehat{}_x^2q\alpha q+\beta ],$$

where $`q=q(t,x)`$ is a complex-valued function of the two real variables $`t`$ and $`x`$, $`t`$ represents time, and $`x`$ represents space. $`q(t,x)`$ is subject to periodic boundary condition of period $`2\pi `$, and even constraint, i.e.,

$$q(t,x+2\pi )=q(t,x),q(t,x)=q(t,x).$$

$`\omega `$ is a positive constant, $`\alpha >0`$ and $`\beta >0`$ are constants, $`\widehat{}_x^2`$ is a bounded Fourier multiplier,

$$\widehat{}_x^2q=\underset{k=1}{\overset{N}{}}k^2\xi _k\stackrel{~}{q}_k\mathrm{cos}kx,$$

$`\xi _k=1`$ when $`kN`$, $`\xi _k=8k^2`$ when $`k>N`$, for some fixed large $`N`$, and $`\epsilon >0`$ is the perturbation parameter. The following theorem was proved in .

###### Theorem 2.1.

There exists a $`\epsilon _0>0`$, such that for any $`\epsilon (0,\epsilon _0)`$, there exists a codimension 1 surface in the external parameter space $`(\alpha ,\beta ,\omega )^+\times ^+\times ^+`$ where $`\omega (\frac{1}{2},1)`$, and $`\alpha \omega <\beta `$. For any $`(\alpha ,\beta ,\omega )`$ on the codimension 1 surface, the regularly perturbed nonlinear Schrödinger equation (2.1) possesses a symmetric pair of Silnikov homoclinic orbits asymptotic to a saddle $`Q_ϵ`$. The codimension 1 surface has the approximate representation given by $`\alpha =1/\kappa (\omega )`$, where $`\kappa (\omega )`$ is plotted in Figure 1.

To prove the theorem, one starts from the invariant plane

$$\mathrm{\Pi }=\{q_xq=0\}.$$

On $`\mathrm{\Pi }`$, there is a saddle $`Q_\epsilon =\sqrt{I}e^{i\theta }`$ to which the symmetric pair of Silnikov homoclinic orbits will be asymptotic to, where

(2.2)
$$I=\omega ^2ϵ\frac{1}{2\omega }\sqrt{\beta ^2\alpha ^2\omega ^2}+\mathrm{},\mathrm{cos}\theta =\frac{\alpha \sqrt{I}}{\beta },\theta (0,\frac{\pi }{2}).$$

Its eigenvalues are

(2.3)
$$\lambda _n^\pm =\epsilon [\alpha +\xi _nn^2]\pm 2\sqrt{(\frac{n^2}{2}+\omega ^2I)(3I\omega ^2\frac{n^2}{2})},$$

where $`n=0,1,2,\mathrm{}`$, $`\omega (\frac{1}{2},1)`$, $`\xi _n=1`$ when $`nN`$, $`\xi _n=8n^2`$ when $`n>N`$, for some fixed large $`N`$, and $`I`$ is given in (2.2). The crucial points to notice are: (1). only $`\lambda _0^+`$ and $`\lambda _1^+`$ have positive real parts, $`\text{Re}\{\lambda _0^+\}<\text{Re}\{\lambda _1^+\}`$; (2). all the other eigenvalues have negative real parts among which the absolute value of $`\text{Re}\{\lambda _2^+\}=\text{Re}\{\lambda _2^{}\}`$ is the smallest; (3). $`|\text{Re}\{\lambda _2^+\}|<\text{Re}\{\lambda _0^+\}`$. Actually, items (2) and (3) are the main characteristics of Silnikov homoclinic orbits.

The unstable manifold $`W^u(Q_\epsilon )`$ of $`Q_\epsilon `$ has a fiber representation . The Melnikov measurement measures the signed distance between $`W^u(Q_\epsilon )`$ and a center-stable manifold $`W_\epsilon ^{cs}`$ . By virtue of a Fenichel Fiber Theorem , one can show that, to the leading order in $`\epsilon `$, the signed distance is given by the Melnikov integral

$`M`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle _0^{2\pi }}[_qF_1(q_0(t))(\widehat{}_x^2q_0(t)\alpha q_0(t)+\beta )`$
$`+_{\overline{q}}F_1(q_0(t))(\widehat{}_x^2\overline{q_0(t)}\alpha \overline{q_0(t)}+\beta )]dxdt,`$
where $`q_0(t)`$ is a homoclinic orbit for the integrable NLS, which can be generated through Darboux transformations ; and $`_qF_1`$ and $`_{\overline{q}}F_1`$ are Melnikov vectors which can be generated through Floquet discriminant in the isospectral theory of the integrable NLS . The zero of the signed distance implies the existence of an orbit in $`W^u(Q_\epsilon )W_\epsilon ^{cs}`$. The stable manifold $`W^s(Q_\epsilon )`$ of $`Q_\epsilon `$ is a codimension 1 submanifold in $`W_\epsilon ^{cs}`$. To locate a homoclinic orbit, one needs to set up a second measurement measuring the signed distance between the orbit in $`W^u(Q_\epsilon )W_\epsilon ^{cs}`$ and $`W^s(Q_\epsilon )`$ inside $`W_\epsilon ^{cs}`$. To set up this signed distance, first one can rather easily track the (perturbed) orbit by an unperturbed orbit to an $`O(\epsilon )`$ neighborhood of $`\mathrm{\Pi }`$, then one needs to prove the size of $`W^s(Q_\epsilon )`$ to be $`O(\epsilon ^\nu )`$ ($`\nu <1`$) with normal form transform. To the leading order in $`\epsilon `$, the zero of the second signed distance is given by

$$\beta \mathrm{cos}\gamma =\frac{\alpha \omega (\mathrm{\Delta }\gamma )}{2\mathrm{sin}\frac{\mathrm{\Delta }\gamma }{2}},$$

where $`\mathrm{\Delta }\gamma =4\vartheta _0`$ and $`\vartheta _0`$ is a phase shift of $`q_0(t)`$. To the leading order in $`\epsilon `$, the common zero of the two second signed distances satisfies $`\alpha =1/\kappa (\omega )`$, where $`\kappa (\omega )`$ is plotted in Figure 1. Then the claim of the theorem is proved by virtue of the implicit function theorem. For rigorous details, see .

### 2.2. Silnikov Homoclinic Orbits in NLS Under Singular Perturbations

Consider the singularly perturbed nonlinear Schrödinger equation ,

(2.4)
$$iq_t=q_{xx}+2[|q|^2\omega ^2]q+i\epsilon [q_{xx}\alpha q+\beta ],$$

where $`q=q(t,x)`$ is a complex-valued function of the two real variables $`t`$ and $`x`$, $`t`$ represents time, and $`x`$ represents space. $`q(t,x)`$ is subject to periodic boundary condition of period $`2\pi `$, and even constraint, i.e.,

$$q(t,x+2\pi )=q(t,x),q(t,x)=q(t,x).$$

$`\omega (1/2,1)`$ is a positive constant, $`\alpha >0`$ and $`\beta >0`$ are constants, and $`\epsilon >0`$ is the perturbation parameter. The following theorem was proved in .

###### Theorem 2.2.

There exists a $`\epsilon _0>0`$, such that for any $`\epsilon (0,\epsilon _0)`$, there exists a codimension 1 surface in the external parameter space $`(\alpha ,\beta ,\omega )^+\times ^+\times ^+`$ where $`\omega (\frac{1}{2},1)/S`$, $`S`$ is a finite subset, and $`\alpha \omega <\beta `$. For any $`(\alpha ,\beta ,\omega )`$ on the codimension 1 surface, the singularly perturbed nonlinear Schrödinger equation (2.4) possesses a symmetric pair of Silnikov homoclinic orbits asymptotic to a saddle $`Q_ϵ`$. The codimension 1 surface has the approximate representation given by $`\alpha =1/\kappa (\omega )`$, where $`\kappa (\omega )`$ is plotted in Figure 1.

The proof of the theorem is also completed through two measurements: the Melnikov measurement and the second measurement. But more powerful machineries are needed .

### 2.3. Silnikov Homoclinic Orbits in Vector NLS Under Perturbations

In recent years, novel results have been obtained on the solutions of the vector nonlinear Schrödinger equations   . Abundant ordinary integrable results have been carried through  , including linear stability calculations . Specifically, the vector nonlinear Schrödinger equations can be written as

$`ip_t+p_{xx}+{\displaystyle \frac{1}{2}}(|p|^2+\chi |q|^2)p=0,`$
$`iq_t+q_{xx}+{\displaystyle \frac{1}{2}}(\chi |p|^2+|q|^2)q=0,`$
where $`p`$ and $`q`$ are complex valued functions of the two real variables $`t`$ and $`x`$, and $`\chi `$ is a positive constant. These equations describe the evolution of two orthogonal pulse envelopes in birefringent optical fibers  , with industrial applications in fiber communication systems  and all-optical switching devices . For linearly birefringent fibers , $`\chi =2/3`$. For elliptically birefringent fibers, $`\chi `$ can take other positive values . When $`\chi =1`$, these equations are first shown to be integrable by S. Manakov , and thus called Manakov equations. When $`\chi `$ is not 1 or 0, these equations are non-integrable. Propelled by the industrial applications, extensive mathematical studies on the vector nonlinear Schrödinger equations have been conducted. Like the scalar nonlinear Schrödinger equation, the vector nonlinear Schrödinger equations also possess figure eight structures in their phase space. Consider the singularly perturbed vector nonlinear Schrödinger equations,

(2.5) $`ip_t+p_{xx}+{\displaystyle \frac{1}{2}}[(|p|^2+|q|^2)\omega ^2]p=i\epsilon [p_{xx}\alpha p\beta ],`$
(2.6) $`iq_t+q_{xx}+{\displaystyle \frac{1}{2}}[(|p|^2+|q|^2)\omega ^2]q=i\epsilon [q_{xx}\alpha q\beta ],`$
where $`p(t,x)`$ and $`q(t,x)`$ are subject to periodic boundary condition of period $`2\pi `$, and are even in $`x`$, i.e.

$$p(t,x+2\pi )=p(t,x),p(t,x)=p(t,x),$$

$$q(t,x+2\pi )=q(t,x),q(t,x)=q(t,x),$$

$`\omega (1,2)`$, $`\alpha >0`$ and $`\beta `$ are real constants, and $`\epsilon >0`$ is the perturbation parameter. We have

###### Theorem 2.3 ().

There exists a $`\epsilon _0>0`$, such that for any $`\epsilon (0,\epsilon _0)`$, there exists a codimension 1 surface in the space of $`(\alpha ,\beta ,\omega )^+\times ^+\times ^+`$ where $`\omega (1,2)/S`$, $`S`$ is a finite subset, and $`\alpha \omega <\sqrt{2}\beta `$. For any $`(\alpha ,\beta ,\omega )`$ on the codimension 1 surface, the singularly perturbed vector nonlinear Schrödinger equations (2.5)-(2.6) possesses a homoclinic orbit asymptotic to a saddle $`Q_ϵ`$. This orbit is also the homoclinic orbit for the singularly perturbed scalar nonlinear Schrödinger equation studied in last section, and is the only one asymptotic to the saddle $`Q_ϵ`$ for the singularly perturbed vector nonlinear Schrödinger equations (2.5)-(2.6). The codimension 1 surface has the approximate representation given by $`\alpha =1/\kappa (\omega )`$, where $`\kappa (\omega )`$ is plotted in Figure 1.

### 2.4. Silnikov Homoclinic Orbits in Discrete NLS Under Perturbations

Consider the following perturbed discrete cubic nonlinear Schrödinger equations ,

(2.7) $`i\dot{q_n}`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}\left[q_{n+1}2q_n+q_{n1}\right]+|q_n|^2(q_{n+1}+q_{n1})2\omega ^2q_n`$
$`+i\epsilon \left[\alpha q_n+{\displaystyle \frac{1}{h^2}}(q_{n+1}2q_n+q_{n1})+\beta \right],`$
where $`i=\sqrt{1}`$, $`q_n`$’s are complex variables,

$$q_{n+N}=q_n,(\text{periodic}\text{condition});\text{and}q_n=q_n,(\text{even}\text{condition});$$

$`h=\frac{1}{N}`$, and

$`N\mathrm{tan}{\displaystyle \frac{\pi }{N}}<\omega <N\mathrm{tan}{\displaystyle \frac{2\pi }{N}},\text{for}N>3,`$
$`3\mathrm{tan}{\displaystyle \frac{\pi }{3}}<\omega <\mathrm{},\text{for}N=3.`$
$`\epsilon [0,\epsilon _1),\alpha (>0),\beta (>0)\text{are}\text{constants.}`$
This is a $`2(M+1)`$ dimensional system, where

$$M=N/2,(N\text{even});\text{and}M=(N1)/2,(N\text{odd}).$$

This system is a finite-difference discretization of the perturbed NLS (2.4). The following theorem was proved in .

Denote by $`\mathrm{\Sigma }_N(N7)`$ the external parameter space,

$`\mathrm{\Sigma }_N`$ $`=`$ $`\left\{(\omega ,\alpha ,\beta )\right|\omega (N\mathrm{tan}{\displaystyle \frac{\pi }{N}},N\mathrm{tan}{\displaystyle \frac{2\pi }{N}}),`$
$`\alpha (0,\alpha _0),\beta (0,\beta _0);`$
$`\text{where}\alpha _0\text{and}\beta _0\text{are}\text{any}\text{fixed}\text{positive}\text{numbers}.\}`$
###### Theorem 2.4.

For any $`N`$ ($`7N<\mathrm{}`$), there exists a positive number $`\epsilon _0`$, such that for any $`\epsilon (0,\epsilon _0)`$, there exists a codimension $`1`$ surface $`E_\epsilon `$ in $`\mathrm{\Sigma }_N`$; for any external parameters ($`\omega ,\alpha ,\beta `$) on $`E_\epsilon `$, there exists a homoclinic orbit asymptotic to a saddle $`Q_\epsilon `$. The codimension $`1`$ surface $`E_\epsilon `$ has the approximate expression $`\alpha =1/\kappa `$, where $`\kappa =\kappa (\omega ;N)`$ is shown in Fig.2.

In the cases ($`3N6`$), $`\kappa `$ is always negative. For $`N7`$, $`\kappa `$ can be positive as shown in Fig.2. When $`N`$ is even and $`7`$, there is in fact a symmetric pair of homoclinic orbits asymptotic to a fixed point $`Q_\epsilon `$ at the same values of the external parameters; since for even $`N`$, we have the symmetry: If $`q_n=f(n,t)`$ solves (2.7), then $`q_n=f(n+N/2,t)`$ also solves (2.7). When $`N`$ is odd and $`7`$, the study can not guarantee that two homoclinic orbits exist at the same value of the external parameters.

### 2.5. Comments on DSII Under Perturbations

Consider the perturbed Davey-Stewartson II equations ,

$`iq_t`$ $`=`$ $`\mathrm{{\rm Y}}q+\left[2(|q|^2\omega ^2)+u_y\right]q+iϵf,`$
$`\mathrm{\Delta }u=4_y|q|^2,`$
where $`q`$ is a complex-valued function of the three variables ($`t,x,y`$), $`u`$ is a real-valued function of the three variables ($`t,x,y`$), $`\mathrm{{\rm Y}}=_{xx}_{yy}`$, $`\mathrm{\Delta }=_{xx}+_{yy}`$, $`\omega >0`$ is a constant. We also consider periodic boundary conditions. The perturbation $`f`$ can be for example a singular perturbation

$$f=\mathrm{\Delta }q\alpha q+\beta ,$$

where $`\alpha >0`$, $`\beta >0`$ are constants, or a regular perturbation by mollifying $`\mathrm{\Delta }`$ into a bounded Fourier multiplier

$$\widehat{\mathrm{\Delta }}q=\underset{kZ^2}{}\beta _k|k|^2\stackrel{~}{q}_k\mathrm{cos}k_1x\mathrm{cos}k_2y,$$

in the case of periods ($`2\pi ,2\pi `$),

$$\beta _k=1,|k|N,\beta _k=|k|^2,|k|>N,$$

for some large $`N`$, $`|k|^2=k_1^2+k_2^2`$.

Under both regular and singular perturbations, the rigorous Melnikov measurement can be established . It turns out that only local well-posedness is necessary for rigorously setting up the Melnikov measurement, thanks to the fact that the unperturbed homoclinic orbits are classical solutions. The obstacle toward proving the existence of homoclinic orbits comes from a technical difficulty in solving a linear system to get the normal form for proving the size estimate of the stable manifold of a saddle. For details, see .

### 2.6. Transversal Homoclinic Orbits in a Periodically Perturbed Sine-Gordon Equation

Transversal homoclinic orbits in continuous systems often appear in two types of systems: (1). periodic systems where the Poincaré period map has a transversal homoclinic orbit, (2). autonomous systems where the homoclinic orbit is asymptotic to a hyperbolic limit cycle.

Consider the periodically perturbed sine-Gordon (SG) equation,

(2.8)
$$u_{tt}=c^2u_{xx}+\mathrm{sin}u+ϵ[au_t+u^3\chi (u)\mathrm{cos}t],$$

where

$$\chi (u)=\{\begin{array}{cc}1,\hfill & uM,\hfill \\ 0,\hfill & u2M,\hfill \end{array}$$

for $`M<u<2M`$, $`\chi (u)`$ is a smooth bump function, under odd periodic boundary condition,

$$u(x+2\pi ,t)=u(x,t),u(x,t)=u(x,t),$$

$`\frac{1}{4}<c^2<1`$, $`a>0`$, $`ϵ`$ is a small perturbation parameter.

###### Theorem 2.5 (, ).

There exists an interval $`I^+`$ such that for any $`aI`$, there exists a transversal homoclinic orbit $`u=\xi (x,t)`$ asymptotic to $`0`$ in $`H^1`$.

### 2.7. Transversal Homoclinic Orbits in a Derivative NLS

Consider the derivative nonlinear Schrödinger equation,

(2.9)
$$iq_t=q_{xx}+2|q|^2q+i\epsilon \left[(\frac{9}{16}|q|^2)q+\mu |\widehat{}_xq|^2\overline{q}\right],$$

where $`q`$ is a complex-valued function of two real variables $`t`$ and $`x`$, $`\epsilon >0`$ is the perturbation parameter, $`\mu `$ is a real constant, and $`\widehat{}_x`$ is a bounded Fourier multiplier,

$$\widehat{}_xq=\underset{k=1}{\overset{K}{}}k\stackrel{~}{q}_k\mathrm{sin}kx,\text{for}q=\underset{k=0}{\overset{\mathrm{}}{}}\stackrel{~}{q}_k\mathrm{cos}kx,$$

and some fixed $`K`$. Periodic boundary condition and even constraint are imposed,

$$q(t,x+2\pi )=q(t,x),q(t,x)=q(t,x).$$

###### Theorem 2.6 ().

There exists a $`\epsilon _0>0`$, such that for any $`\epsilon (0,\epsilon _0)`$, and $`|\mu |>5.8`$, there exist two transversal homoclinic orbits asymptotic to the limit cycle $`q_c=\frac{3}{4}\mathrm{exp}\{i[\frac{9}{8}t+\gamma ]\}`$.

## 3. Existence of Chaos

The importance of homoclinic orbits with respect to chaotic dynamics was first realized by Poincaré . In 1961, Smale constructed the well-known horseshoe in the neighborhood of a transversal homoclinic orbit   . In particular, Smale’s theorem implies Birkhoff’s theorem on the existence of a sequence of structurely stable periodic orbits in the neighborhood of a transversal homoclinic orbit . In 1984 and 1988  , Palmer gave a beautiful proof of Smale’s theorem using a shadowing lemma. Later, this proof was generalized to infinite dimensions by Steinlein and Walther   and Henry . In 1967, Silnikov proved Smale’s theorem for autonomous systems in finite dimensions using a fixed point argument . In 1996, Palmer proved Smale’s theorem for autonomous systems in finite dimensions using shadowing lemma  . In 2002, Li proved Smale’s theorem for autonomous systems in infinite dimensions using shadowing lemma .

For nontransversal homoclinic orbits, the most well-known type which leads to the existence of Smale horseshoes is the so-called Silnikov homoclinic orbit     . Existence of Silnikov homoclinic orbits and new constructions of Smale horseshoes for concrete nonlinear wave systems have been established in finite dimensions   and in infinite dimensions   .

### 3.1. Shift Automorphism

Let $`𝒲`$ be a set which consists of elements of the doubly infinite sequence form:

$$a=(\mathrm{}a_2a_1a_0,a_1a_2\mathrm{}),$$

where $`a_k\{l_1,\mathrm{},l_m\}`$, m labels, and $`k`$. We introduce a topology in $`𝒲`$ by taking as neighborhood basis of

$$a^{}=(\mathrm{}a_2^{}a_1^{}a_0^{},a_1^{}a_2^{}\mathrm{}),$$

the set

$$W_j=\left\{a𝒲\right|a_k=a_k^{}(|k|<j)\}$$

for $`j=1,2,\mathrm{}`$. This makes $`𝒲`$ a topological space. The shift automorphism $`\chi `$ is defined on $`𝒲`$ by

$`\chi `$ $`:`$ $`𝒲𝒲,`$
$`a𝒲,\chi (a)=b,\text{where}b_k=a_{k+1}.`$
The shift automorphism $`\chi `$ exhibits sensitive dependence on initial conditions, which is a hallmark of chaos.

### 3.2. NLS Under Regular Perturbations

We continue from section 2.1 and consider the regularly perturbed nonlinear Schrödinger (NLS) equation (2.1). Starting from the Homoclinic Orbit Theorem 2.1, We have the following theorem on the existence of chaos.

###### Theorem 3.1 (Chaos Theorem, ).

Under certain generic assumptions for the perturbed nonlinear Schrödinger system (2.1), there exists a compact Cantor subset $`\mathrm{\Lambda }`$ of $`H^1`$ (the Sobolev space), $`\mathrm{\Lambda }`$ consists of points, and is invariant under a Poincaré map $`P`$. $`P`$ restricted to $`\mathrm{\Lambda }`$, is topologically conjugate to the shift automorphism $`\chi `$ on four symbols $`1,2,1,2`$. That is, there exists a homeomorphism

$$\varphi :𝒲\mathrm{\Lambda },$$

such that the following diagram commutes:

$$\begin{array}{ccc}𝒲& \stackrel{\varphi }{}& \mathrm{\Lambda }\\ \chi & & P\\ 𝒲& \underset{\varphi }{}& \mathrm{\Lambda }\end{array}$$

Although the symmetric pair of Silnikov homoclinic orbits is not structurally stable, the Smale horseshoes are structurally stable. Thus, the Cantor sets and the conjugacy to the shift automorphism are structurally stable.

### 3.3. NLS Under Singular Perturbations

We continue from section 2.2 and consider the singularly perturbed nonlinear Schrödinger (NLS) equation (2.4). Starting from the Homoclinic Orbit Theorem 2.2, we have

###### Theorem 3.2 (Chaos Theorem, ).

Under certain generic assumptions for the perturbed nonlinear Schrödinger system (2.4), there exists a compact Cantor subset $`\mathrm{\Lambda }`$ of $`H^1`$, $`\mathrm{\Lambda }`$ consists of points, and is invariant under a Poincaré map $`P`$. $`P`$ restricted to $`\mathrm{\Lambda }`$, is topologically conjugate to the shift automorphism $`\chi `$ on four symbols $`1,2,1,2`$. That is, there exists a homeomorphism

$$\varphi :𝒲\mathrm{\Lambda },$$

such that the following diagram commutes:

$$\begin{array}{ccc}𝒲& \stackrel{\varphi }{}& \mathrm{\Lambda }\\ \chi & & P\\ 𝒲& \underset{\varphi }{}& \mathrm{\Lambda }\end{array}$$

### 3.4. Discrete NLS Under Perturbations

We continue from section 2.4. Starting from the Homoclinic Orbit Theorem 2.4, we have the theorem on the existence of chaos for the perturbed discrete nonlinear Schrödinger equation (2.7).

###### Theorem 3.3 (Chaos Theorem, ).

Under certain generic assumptions for the perturbed discrete nonlinear Schrödinger system (2.7), there exists a compact Cantor subset $`\mathrm{\Lambda }`$ of $`^{2(M+1)}`$, $`\mathrm{\Lambda }`$ consists of points, and is invariant under a Poincaré map $`P`$. $`P`$ restricted to $`\mathrm{\Lambda }`$, is topologically conjugate to the shift automorphism $`\chi `$ on two symbols $`0,1`$. That is, there exists a homeomorphism

$$\varphi :𝒲\mathrm{\Lambda },$$

such that the following diagram commutes:

$$\begin{array}{ccc}𝒲& \stackrel{\varphi }{}& \mathrm{\Lambda }\\ \chi & & P\\ 𝒲& \underset{\varphi }{}& \mathrm{\Lambda }\end{array}$$

### 3.5. Numerical Simulation of Chaos

The finite-difference discretization of both the regularly and the singularly perturbed nonlinear Schrödinger equations (2.1) and (2.4) leads to the same discrete perturbed nonlinear Schrödinger equation (2.7).

In the chaotic regime, typical numerical output is shown in Figure 3. Notice that there are two typical profiles at a fixed time: (1). a breather type profile with its hump located at the center of the spatial period interval, (2). a breather type profile with its hump located at the boundary (wing) of the spatial period interval. These two types of profiles are half spatial period translate of each other. If we label the profiles, with their humps at the center of the spatial period interval, by “C”; those profiles, with their humps at the wing of the spatial period interval, by “W”, then

(3.1)
$$\mathrm{`}\mathrm{`}W\text{}=\sigma \mathrm{`}\mathrm{`}C\text{},$$

where $`\sigma `$ is the symmetry group element representing half spatial period translate. The time series of the output in Figure 3 is a chaotic jumping between “C” and “W”, which we call “chaotic center-wing jumping”. We interpret the chaotic center-wing jumping as the numerical realization of the shift automorphism $`\chi `$ on symbols. We can make this more precise in terms of the phase space geometry. The figure-eight structure of the integrable NLS projected onto the plane of the Fourier component $`\mathrm{cos}x`$ is illustrated in Figure 4, and labeled by $`L_C`$ and $`L_W`$. From the symmetry, we know that

$$L_W=\sigma L_C.$$

$`L_C`$ has the spatial-temporal profile realization as in Figure 5 with the hump located at the center. $`L_W`$ corresponds to the half spatial period translate of the spatial-temporal profile realization as in Figure 5, with the hump located at the boundary (wing). An orbit inside $`L_C`$, $`L_{Cin}`$ has a spatial-temporal profile realization as in Figure 6. The half period translate of $`L_{Cin}`$, $`L_{Win}`$ is inside $`L_W`$. An orbit outside $`L_C`$ and $`L_W`$, $`L_{out}`$ has the spatial-temporal profile realization as in Figure 7. $`S_l`$ and $`S_{l,\sigma }`$ are two phase blocks. From these figures, we can see clearly that the chaotic center-wing jumping (Figure 3) is the realization of the shift automorphism on symbols in $`S_lS_{l,\sigma }`$.

### 3.6. Shadowing Lemma and Chaos in Finite-D Periodic Systems

Since its invention , shadowing lemma has been a useful tool for solving many dynamical system problems. Here we only focus upon its use in proving the existence of chaos in a neighborhood of a transversal homoclinic orbit. According to the type of the system, the level of difficulty in proving the existence of chaos with a shadowing lemma is different.

A finite-dimensional periodic system can be written in the general form

$$\dot{x}=F(x,t),$$

where $`x^n`$, and $`F(x,t)`$ is periodic in $`t`$. Let $`f`$ be the Poincaré period map.

###### Definition 3.4.

A doubly infinite sequence $`\{y_j\}`$ in $`^n`$ is a $`\delta `$ pseudo-orbit of a $`C^1`$ diffeomorphism $`f:^n^n`$ if for all integers $`j`$

$$|y_{j+1}f(y_j)|\delta .$$

An orbit $`\{f^j(x)\}`$ is said to $`\epsilon `$-shadow the $`\delta `$ pseudo-orbit $`\{y_j\}`$ if for all integers $`j`$

$$|f^j(x)y_j|\epsilon .$$

###### Definition 3.5.

A compact invariant set $`S`$ is hyperbolic if there are positive constants $`K`$, $`\alpha `$ and a projection matrix valued function $`P(x)`$, $`xS`$, of constant rank such that for all $`x`$ in $`S`$

$$P(f(x))Df(x)=Df(x)P(x),$$

$$|Df^j(x)P(x)|Ke^{\alpha j},(j0),$$

$$|Df^j(x)(IP(x))|Ke^{\alpha j},(j0).$$

###### Theorem 3.6 (Shadowing Lemma ).

Let $`S`$ be a compact hyperbolic set for the $`C^1`$ diffeomorphism $`f:^n^n`$. Then given $`\epsilon >0`$ sufficiently small there exists $`\delta >0`$ such that every $`\delta `$ pseudo-orbit in $`S`$ has a unique $`\epsilon `$-shadowing orbit.

The proof of this theorem by Palmer  is overall a fixed point argument with the help of Green functions for linear maps.

Let $`y_0`$ be a transversal homoclinic point asymptotic to a saddle $`x_0`$ of a $`C^1`$ diffeomorphism $`f:^n^n`$. Then the set

$$S=\{x_0\}\{f^j(y_0):jZ\}$$

is hyperbolic. Denote by $`A_0`$ and $`A_1`$ the two orbit segments of length $`2m+1`$

$$A_0=\{x_0,x_0,\mathrm{},x_0\},A_1=\{f^m(y_0),f^{m+1}(y_0),\mathrm{},f^{m1}(y_0),f^m(y_0)\}.$$

Let

$$a=(\mathrm{},a_1,a_0,a_1,\mathrm{}),$$

where $`a_j\{0,1\}`$, be any doubly infinite binary sequence. Let $`A`$ be the doubly infinite sequence of points in $`S`$, associated with $`a`$

$$A=\{\mathrm{},A_{a_1},A_{a_0},A_{a_1},\mathrm{}\}.$$

When $`m`$ is sufficiently large, $`A`$ is a $`\delta `$ pseudo-orbit in $`S`$. By the shadowing lemma (Theorem 3.6), there is a unique $`\epsilon `$-shadowing orbit that shadows $`A`$. In this manner, Palmer  gave a beautiful proof of Smale’s horseshoe theorem.

###### Definition 3.7.

Denote by $`\mathrm{\Sigma }`$ the set of doubly infinite binary sequences

$$a=(\mathrm{},a_1,a_0,a_1,\mathrm{}),$$

where $`a_j\{0,1\}`$. We give the set $`\{0,1\}`$ the discrete topology and $`\mathrm{\Sigma }`$ the product topology. The Bernoulli shift $`\chi `$ is defined by

$$[\chi (a)]_j=a_{j+1}.$$

###### Theorem 3.8.

Let $`y_0`$ be a transversal homoclinic point asymptotic to a saddle $`x_0`$ of a $`C^1`$ diffeomorphism $`f:^n^n`$. Then there is a homeomorphism $`\varphi `$ of $`\mathrm{\Sigma }`$ onto a compact subset of $`^n`$ which is invariant under $`f`$ and such that when $`m`$ is sufficiently large

$$f^{2m+1}\varphi =\varphi \chi ,$$

that is, the action of $`f^{2m+1}`$ on $`\varphi (\mathrm{\Sigma })`$ is topologically conjugate to the action of $`\chi `$ on $`\mathrm{\Sigma }`$.

Here one can define $`\varphi (a)`$ to be the point on the shadowing orbit that shadows the midpoint of the orbit segment $`A_{a_0}`$, which is either $`x_0`$ or $`y_0`$. The topological conjugacy can be easily verified. For details, see . Other references can be found in  .

### 3.7. Shadowing Lemma and Chaos in Infinite-D Periodic Systems

An infinite-dimensional periodic system defined in a Banach space $`X`$ can be written in the general form

$$\dot{x}=F(x,t),$$

where $`xX`$, and $`F(x,t)`$ is periodic in $`t`$. Let $`f`$ be the Poincaré period map. When $`f`$ is a $`C^1`$ map which needs not to be invertible, shadowing lemma and symbolic dynamics around a transversal homoclinic orbit can both be established   . Other references can be found in    . There exists also a work on horseshoe construction without shadowing lemma for sinusoidally forced vibrations of buckled beam .

### 3.8. Periodically Perturbed Sine-Gordon (SG) Equation

We continue from section 2.6, and use the notations in section 3.6. For the periodically perturbed sine-Gordon equation (2.8), the Poincaré period map is a $`C^1`$ diffeomorphism in $`H^1`$. As a corollary of the result in last section, we have the theorem on the existence of chaos.

###### Theorem 3.9.

There is an integer $`m`$ and a homeomorphism $`\varphi `$ of $`\mathrm{\Sigma }`$ onto a compact Cantor subset $`\mathrm{\Lambda }`$ of $`H^1`$. $`\mathrm{\Lambda }`$ is invariant under the Poincaré period-$`2\pi `$ map $`P`$ of the periodically perturbed sine-Gordon equation (2.8). The action of $`P^{2m+1}`$ on $`\mathrm{\Lambda }`$ is topologically conjugate to the action of $`\chi `$ on $`\mathrm{\Sigma }:P^{2m+1}\varphi =\varphi \chi `$. That is, the following diagram commutes:

$$\begin{array}{ccc}\mathrm{\Sigma }& \stackrel{\varphi }{}& \mathrm{\Lambda }\\ \chi & & P^{2m+1}\\ \mathrm{\Sigma }& \underset{\varphi }{}& \mathrm{\Lambda }\end{array}$$

### 3.9. Shadowing Lemma and Chaos in Finite-D Autonomous Systems

A finite-dimensional autonomous system can be written in the general form

$$\dot{x}=F(x),$$

where $`x^n`$. In this case, a transversal homoclinic orbit can be an orbit asymptotic to a normally hyperbolic limit cycle. That is, it is an orbit in the intersection of the stable and unstable manifolds of a normally hyperbolic limit cycle. Instead of the Poincaré period map as for periodic system, one may want to introduce the so-called Poincaré return map which is a map induced by the flow on a codimension 1 section which is transversal to the limit cycle. Unfortunately, such a map is not even well-defined in the neighborhood of the homoclinic orbit. This poses a challenging difficulty in extending the arguments as in the case of a Poincaré period map. In 1996, Palmer  completed a proof of a shadowing lemma and existence of chaos using Newton’s method. It will be difficult to extend this method to infinite dimensions, since it used heavily differentiations in time. Other references can be found in    .

### 3.10. Shadowing Lemma and Chaos in Infinite-D Autonomous Systems

An infinite-dimensional autonomous system defined in a Banach space $`X`$ can be written in the general form

$$\dot{x}=F(x),$$

where $`xX`$. In 2002, the author  completed a proof of a shadowing lemma and existence of chaos using Fenichel’s persistence of normally hyperbolic invariant manifold idea. The setup is as follows,

* Assumption (A1): There exist a hyperbolic limit cycle $`S`$ and a transversal homoclinic orbit $`\xi `$ asymptotic to $`S`$. As curves, $`S`$ and $`\xi `$ are $`C^3`$.
* Assumption (A2): The Fenichel fiber theorem is valid at $`S`$. That is, there exist a family of unstable Fenichel fibers $`\{^u(q):qS\}`$ and a family of stable Fenichel fibers $`\{^s(q):qS\}`$. For each fixed $`qS`$, $`^u(q)`$ and $`^s(q)`$ are $`C^3`$ submanifolds. $`^u(q)`$ and $`^s(q)`$ are $`C^2`$ in $`q,qS`$. The unions $`_{qS}^u(q)`$ and $`_{qS}^s(q)`$ are the unstable and stable manifolds of $`S`$. Both families are invariant, i.e.

$$F^t(^u(q))^u(F^t(q)),t0,qS,$$

$$F^t(^s(q))^s(F^t(q)),t0,qS,$$

where $`F^t`$ is the evolution operator. There are positive constants $`\kappa `$ and $`\widehat{C}`$ such that $`qS`$, $`q^{}^u(q)`$ and $`q^+^s(q)`$,

$$F^t(q^{})F^t(q)\widehat{C}e^{\kappa t}q^{}q,t0,$$

$$F^t(q^+)F^t(q)\widehat{C}e^{\kappa t}q^+q,t0.$$
* Assumption (A3): $`F^t(q)`$ is $`C^0`$ in $`t`$, for $`t(\mathrm{},\mathrm{})`$, $`qX`$. For any fixed $`t(\mathrm{},\mathrm{})`$, $`F^t(q)`$ is a $`C^2`$ diffeomorphism on $`X`$.

###### Remark 3.10.

Notice that we do not assume that as functions of time, $`S`$ and $`\xi `$ are $`C^3`$ , and we only assume that as curves, $`S`$ and $`\xi `$ are $`C^3`$.

Under the above setup, a shadowing lemma and existence of chaos can be proved . Another crucial element in the argument is the establishment of a $`\lambda `$-lemma (also called inclination lemma) .

### 3.11. A Derivative Nonlinear Schrödinger Equation

We continue from section 2.7, and consider the derivative nonlinear Schrödinger equation (2.9). The transversal homoclinic orbit given in Theorem 2.6 is a classical solution. Thus, Assumption (A1) is valid. Assumption (A2) follows from the standard arguments in   . Since the perturbation in (2.9) is bounded, Assumption (A3) follows from standard arguments. Thus there exists chaos in the derivative nonlinear Schrödinger equation (2.9) .

## 4. Lax Pairs of Euler Equations of Inviscid Fluids

The governing equations for the incompressible viscous fluid flow are the Navier-Stokes equations. Turbulence occurs in the regime of high Reynolds number. By formally setting the Reynolds number equal to infinity, the Navier-Stokes equations reduce to the Euler equations of incompressible inviscid fluid flow. One may view the Navier-Stokes equations with large Reynolds number as a singular perturbation of the Euler equations.

Results of T. Kato show that 2D Navier-Stokes equations are globally well-posed in $`C^0([0,\mathrm{});H^s(R^2)),s>2`$, and for any $`0<T<\mathrm{}`$, the mild solutions of the 2D Navier-Stokes equations approach those of the 2D Euler equations in $`C^0([0,T];H^s(R^2))`$ . 3D Navier-Stokes equations are locally well-posed in $`C^0([0,\tau ];H^s(R^3)),s>5/2`$, and the mild solutions of the 3D Navier-Stokes equations approach those of the 3D Euler equations in $`C^0([0,\tau ];H^s(R^3))`$, where $`\tau `$ depends on the norms of the initial data and the external force  . Extensive studies on the inviscid limit have been carried by J. Wu et al.    . There is no doubt that mathematical study on Navier-Stokes (Euler) equations is one of the most important mathematical problems. In fact, Clay Mathematics Institute has posted the global well-posedness of 3D Navier-Stokes equations as one of the one million dollars problems.

V. Arnold  realized that 2D Euler equations are a Hamiltonian system. Recently, the author found Lax pair structures for Euler equations     . Understanding the structures of solutions to Euler equations is of fundamental interest. Of particular interest is the question on the global well-posedness of 3D Navier-Stokes and Euler equations. Our number one hope is that the Lax pair structures can be useful in investigating the global well-posedness. Our secondary hope is that the Darboux transformation  associated with the Lax pair can generate explicit representation of homoclinic structures .

The philosophical significance of the existence of Lax pairs for Euler equations is even more important. If one defines integrability of an equation by the existence of a Lax pair, then both 2D and 3D Euler equations are integrable. More importantly, both 2D and 3D Navier-Stokes equations at high Reynolds numbers are singularly perturbed integrable systems. Such a point of view changes our old ideology on Euler and Navier-Stokes equations.

### 4.1. A Lax Pair for 2D Euler Equation

The 2D Euler equation can be written in the vorticity form,

(4.1)
$$_t\mathrm{\Omega }+\{\mathrm{\Psi },\mathrm{\Omega }\}=0,$$

where the bracket $`\{,\}`$ is defined as

$$\{f,g\}=(_xf)(_yg)(_yf)(_xg),$$

$`\mathrm{\Omega }`$ is the vorticity, and $`\mathrm{\Psi }`$ is the stream function given by,

$$u=_y\mathrm{\Psi },v=_x\mathrm{\Psi },$$

and the relation between vorticity $`\mathrm{\Omega }`$ and stream function $`\mathrm{\Psi }`$ is,

$$\mathrm{\Omega }=_xv_yu=\mathrm{\Delta }\mathrm{\Psi }.$$

###### Theorem 4.1 (Li, ).

The Lax pair of the 2D Euler equation (4.1) is given as

(4.2)
$$\{\begin{array}{c}L\phi =\lambda \phi ,\hfill \\ _t\phi +A\phi =0,\hfill \end{array}$$

where

$$L\phi =\{\mathrm{\Omega },\phi \},A\phi =\{\mathrm{\Psi },\phi \},$$

and $`\lambda `$ is a complex constant, and $`\phi `$ is a complex-valued function.

### 4.2. A Darboux Transformation for 2D Euler Equation

Consider the Lax pair (4.2) at $`\lambda =0`$, i.e.

(4.3) $`\{\mathrm{\Omega },p\}=0,`$
(4.4) $`_tp+\{\mathrm{\Psi },p\}=0,`$
where we replaced the notation $`\phi `$ by $`p`$.

###### Theorem 4.2 ().

Let $`f=f(t,x,y)`$ be any fixed solution to the system (4.3, 4.4), we define the Gauge transform $`G_f`$:

(4.5)
$$\stackrel{~}{p}=G_fp=\frac{1}{\mathrm{\Omega }_x}[p_x(_x\mathrm{ln}f)p],$$

and the transforms of the potentials $`\mathrm{\Omega }`$ and $`\mathrm{\Psi }`$:

(4.6)
$$\stackrel{~}{\mathrm{\Psi }}=\mathrm{\Psi }+F,\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }+\mathrm{\Delta }F,$$

where $`F`$ is subject to the constraints

(4.7)
$$\{\mathrm{\Omega },\mathrm{\Delta }F\}=0,\{\mathrm{\Delta }F,F\}=0.$$

Then $`\stackrel{~}{p}`$ solves the system (4.3, 4.4) at $`(\stackrel{~}{\mathrm{\Omega }},\stackrel{~}{\mathrm{\Psi }})`$. Thus (4.5) and (4.6) form the Darboux transformation for the 2D Euler equation (4.1) and its Lax pair (4.3, 4.4).

###### Remark 4.3.

For KdV equation and many other soliton equations, the Gauge transform is of the form ,

$$\stackrel{~}{p}=p_x(_x\mathrm{ln}f)p.$$

In general, Gauge transform does not involve potentials. For 2D Euler equation, a potential factor $`\frac{1}{\mathrm{\Omega }_x}`$ is needed. From (4.3), one has

$$\frac{p_x}{\mathrm{\Omega }_x}=\frac{p_y}{\mathrm{\Omega }_y}.$$

The Gauge transform (4.5) can be rewritten as

$$\stackrel{~}{p}=\frac{p_x}{\mathrm{\Omega }_x}\frac{f_x}{\mathrm{\Omega }_x}\frac{p}{f}=\frac{p_y}{\mathrm{\Omega }_y}\frac{f_y}{\mathrm{\Omega }_y}\frac{p}{f}.$$

The Lax pair (4.3, 4.4) has a symmetry, i.e. it is invariant under the transform $`(t,x,y)(t,y,x)`$. The form of the Gauge transform (4.5) resulted from the inclusion of the potential factor $`\frac{1}{\mathrm{\Omega }_x}`$, is consistent with this symmetry.

Our hope is to use the Darboux transformation to generate homoclinic structures for 2D Euler equation .

### 4.3. A Lax Pair for Rossby Wave Equation

The Rossby wave equation is

$$_t\mathrm{\Omega }+\{\mathrm{\Psi },\mathrm{\Omega }\}+\beta _x\mathrm{\Psi }=0,$$

where $`\mathrm{\Omega }=\mathrm{\Omega }(t,x,y)`$ is the vorticity, $`\{\mathrm{\Psi },\mathrm{\Omega }\}=\mathrm{\Psi }_x\mathrm{\Omega }_y\mathrm{\Psi }_y\mathrm{\Omega }_x`$, and $`\mathrm{\Psi }=\mathrm{\Delta }^1\mathrm{\Omega }`$ is the stream function. Its Lax pair can be obtained by formally conducting the transformation, $`\mathrm{\Omega }=\stackrel{~}{\mathrm{\Omega }}+\beta y`$, to the 2D Euler equation ,

$$\{\mathrm{\Omega },\phi \}\beta _x\phi =\lambda \phi ,_t\phi +\{\mathrm{\Psi },\phi \}=0,$$

where $`\phi `$ is a complex-valued function, and $`\lambda `$ is a complex parameter.

### 4.4. Lax Pairs for 3D Euler Equation

The 3D Euler equation can be written in vorticity form,

(4.8)
$$_t\mathrm{\Omega }+(u\mathrm{NewA})\mathrm{\Omega }(\mathrm{\Omega }\mathrm{NewA})u=0,$$

where $`u=(u_1,u_2,u_3)`$ is the velocity, $`\mathrm{\Omega }=(\mathrm{\Omega }_1,\mathrm{\Omega }_2,\mathrm{\Omega }_3)`$ is the vorticity, $`\mathrm{NewA}=(_x,_y,_z)`$, $`\mathrm{\Omega }=\mathrm{NewA}\times u`$, and $`\mathrm{NewA}u=0`$. $`u`$ can be represented by $`\mathrm{\Omega }`$ for example through Biot-Savart law.

###### Theorem 4.4.

The Lax pair of the 3D Euler equation (4.8) is given as

(4.9)
$$\{\begin{array}{c}L\varphi =\lambda \varphi ,\hfill \\ _t\varphi +A\varphi =0,\hfill \end{array}$$

where

$$L\varphi =(\mathrm{\Omega }\mathrm{NewA})\varphi ,A\phi =(u\mathrm{NewA})\varphi ,$$

$`\lambda `$ is a complex constant, and $`\varphi `$ is a complex scalar-valued function.

###### Theorem 4.5 ().

Another Lax pair of the 3D Euler equation (4.8) is given as

(4.10)
$$\{\begin{array}{c}L\phi =\lambda \phi ,\hfill \\ _t\phi +A\phi =0,\hfill \end{array}$$

where

$$L\phi =(\mathrm{\Omega }\mathrm{NewA})\phi (\phi \mathrm{NewA})\mathrm{\Omega },A\phi =(u\mathrm{NewA})\phi (\phi \mathrm{NewA})u,$$

$`\lambda `$ is a complex constant, and $`\phi =(\phi _1,\phi _2,\phi _3)`$ is a complex 3-vector valued function.

Our hope is that the infinitely many conservation laws generated by $`\lambda `$ can provide a priori estimates for the global well-posedness of 3D Navier-Stokes equations, or better understanding on the global well-posedness . For more informations on the topics, see  .

## 5. Dynamical System Studies on 2D Euler Equation

### 5.1. Linearized 2D Euler Equation at a Fixed Point

To begin an infinite dimensional dynamical system study, we investigate the linearized 2D Euler equation at a fixed point   . We consider the 2D Euler equation (4.1) under periodic boundary condition in both $`x`$ and $`y`$ directions with period $`2\pi `$. Expanding $`\mathrm{\Omega }`$ into Fourier series,

$$\mathrm{\Omega }=\underset{k^2/\{0\}}{}\omega _ke^{ikX},$$

where $`\omega _k=\overline{\omega _k}`$, $`k=(k_1,k_2)^T`$, and $`X=(x,y)^T`$. The 2D Euler equation can be rewritten as

(5.1)
$$\dot{\omega }_k=\underset{k=p+q}{}A(p,q)\omega _p\omega _q,$$

where $`A(p,q)`$ is given by,

(5.5) $`A(p,q)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[|q|^2|p|^2](p_1q_2p_2q_1)`$
$`=`$ $`{\displaystyle \frac{1}{2}}[|q|^2|p|^2]\left|\begin{array}{cc}p_1\hfill & \hfill q_1\\ p_2\hfill & \hfill q_2\end{array}\right|,`$
where $`|q|^2=q_1^2+q_2^2`$ for $`q=(q_1,q_2)^T`$, similarly for $`p`$.

Denote $`\{\omega _k\}_{k^2/\{0\}}`$ by $`\omega `$. We consider the simple fixed point $`\omega ^{}`$ :

(5.6)
$$\omega _p^{}=\mathrm{\Gamma },\omega _k^{}=0,\text{if}kp\text{or}p,$$

of the 2D Euler equation (5.1), where $`\mathrm{\Gamma }`$ is an arbitrary complex constant. The linearized two-dimensional Euler equation at $`\omega ^{}`$ is given by,

(5.7)
$$\dot{\omega }_k=A(p,kp)\mathrm{\Gamma }\omega _{kp}+A(p,k+p)\overline{\mathrm{\Gamma }}\omega _{k+p}.$$

###### Definition 5.1 (Classes).

For any $`\widehat{k}^2/\{0\}`$, we define the class $`\mathrm{\Sigma }_{\widehat{k}}`$ to be the subset of $`^2/\{0\}`$:

$$\mathrm{\Sigma }_{\widehat{k}}=\left\{\widehat{k}+np^2/\{0\}\right|n,p\text{is specified in (}\text{5.6}\text{)}\}.$$

See Fig.8 for an illustration of the classes. According to the classification defined in Definition 5.1, the linearized two-dimensional Euler equation (5.7) decouples into infinitely many invariant subsystems:

(5.8) $`\dot{\omega }_{\widehat{k}+np}`$ $`=`$ $`A(p,\widehat{k}+(n1)p)\mathrm{\Gamma }\omega _{\widehat{k}+(n1)p}`$
$`+A(p,\widehat{k}+(n+1)p)\overline{\mathrm{\Gamma }}\omega _{\widehat{k}+(n+1)p}.`$
Let $`_{\widehat{k}}`$ be the linear operator defined by the right hand side of (5.8), and $`H^s`$ be the Sobolev space where $`s0`$ is an integer and $`H^0=\mathrm{}_2`$.

###### Theorem 5.2 ( ).

The eigenvalues of the linear operator $`_{\widehat{k}}`$ in $`H^s`$ are of four types: real pairs ($`c,c`$), purely imaginary pairs ($`id,id`$), quadruples ($`\pm c\pm id`$), and zero eigenvalues.

###### Definition 5.3 (The Disk).

The open disk of radius $`|p|`$ in $`^2/\{0\}`$, denoted by $`D_{|p|}`$, is defined as

$$D_{|p|}=\left\{k^2/\{0\}\right||k|<|p|\},$$

and the closure of $`D_{|p|}`$, denoted by $`\overline{D}_{|p|}`$, is defined as

$$\overline{D}_{|p|}=\left\{k^2/\{0\}\right||k||p|\}.$$

###### Theorem 5.4 (The Spectral Theorem,  ).

We have the following claims on the spectrum of the linear operator $`_{\widehat{k}}`$:

1. If $`\mathrm{\Sigma }_{\widehat{k}}\overline{D}_{|p|}=\mathrm{}`$, then the entire $`H^s`$ spectrum of the linear operator $`_{\widehat{k}}`$ is its continuous spectrum. See Figure 9, where $`b=\frac{1}{2}|\mathrm{\Gamma }||p|^2\left|\begin{array}{cc}p_1& \widehat{k}_1\\ p_2& \widehat{k}_2\end{array}\right|.`$ That is, both the residual and the point spectra of $`_{\widehat{k}}`$ are empty.
2. If $`\mathrm{\Sigma }_{\widehat{k}}\overline{D}_{|p|}\mathrm{}`$, then the entire essential $`H^s`$ spectrum of the linear operator $`_{\widehat{k}}`$ is its continuous spectrum. That is, the residual spectrum of $`_{\widehat{k}}`$ is empty. The point spectrum of $`_{\widehat{k}}`$ is symmetric with respect to both real and imaginary axes. See Figure 10.

Denote by $`L`$ the right hand side of (5.7), i.e. the whole linearized 2D Euler operator, the spectral mapping theorem holds.

###### Theorem 5.5 ().

$$\sigma (e^{tL})=e^{t\sigma (L)},t0.$$

Let $`\zeta `$ denote the number of points $`q^2/\{0\}`$ that belong to the open disk of radius $`|p|`$, $`D_{|p|}`$, and such that $`q`$ is not parallel to $`p`$.

###### Theorem 5.6 ().

The number of nonimaginary eigenvalues of $`L`$ (counting the multiplicities) does not exceed $`2\zeta `$.

Another interesting discussion upon the discrete spectrum can be found in .

Since the introduction of continued fractions for calculating the eigenvalues of steady fluid flow, by Meshalkin and Sinai , this topics had been extensively explored        . Rigorous justification on the continued fraction calculation was given in  . As an example, we take $`p=(1,1)^T`$. When $`\mathrm{\Gamma }0`$, the fixed point has $`4`$ eigenvalues which form a quadruple. These four eigenvalues appear in the invariant linear subsystem labeled by $`\widehat{k}=(3,2)^T`$. One of them is :

(5.9)
$$\stackrel{~}{\lambda }=2\lambda /|\mathrm{\Gamma }|=0.24822302478255+i0.35172076526520.$$

See Figure 11 for an illustration. The essential spectrum (= continuous spectrum) of $`_{\widehat{k}}`$ with $`\widehat{k}=(3,2)^T`$ is the segment on the imaginary axis shown in Figure 11, where $`b=\frac{1}{4}\mathrm{\Gamma }`$. The essential spectrum (= continuous spectrum) of the linear 2D Euler operator at this fixed point is the entire imaginary axis.

A rather well-known open problem is proving the existence of unstable, stable, and center manifolds. The main difficulty comes from the fact that the nonlinear term is non-Lipschitzian.

### 5.2. Models

To simplify our study, we study only the case when $`\omega _k`$ is real, $`k^2/\{0\}`$, i.e. we only study the cosine transform of the vorticity,

$$\mathrm{\Omega }=\underset{k^2/\{0\}}{}\omega _k\mathrm{cos}(kX).$$

To further simplify our study, we will study a concrete dashed-line model based upon the line of fixed points (5.6) with the mode $`p=(1,1)^T`$ parametrized by $`\mathrm{\Gamma }`$, and the only unstable invariant linear subsystem labeled by $`\widehat{k}=(3,2)^T`$. The so-called dashed-line model is given by ,

$`\dot{\omega }_n`$ $`=`$ $`ϵ_{n1}A_{n1}\omega _p\omega _{n1}ϵ_{n+1}A_{n+1}\omega _p\omega _{n+1},`$
$`\dot{\omega }_p`$ $`=`$ $`{\displaystyle \underset{nZ}{}}ϵ_nϵ_{n1}A_{n1,n}\omega _{n1}\omega _n,`$
where

(5.11)
$$\omega _n=\omega _{\widehat{k}+np},A_n=A(p,\widehat{k}+np),A_{m,n}=A(\widehat{k}+mp,\widehat{k}+np),$$

$$ϵ_n=\{\begin{array}{cc}1,\hfill & \text{if }n5j,jZ,\hfill \\ ϵ,\hfill & \text{if }n=5j,\text{ for some }jZ.\hfill \end{array}$$

The model is designed to model the hyperbolic structure of the 2D Euler equation, connected to the line of fixed points (5.6) with $`p=(1,1)^T`$. Figure 12 illustrates the collocation of the modes in this model, which has the “dashed-line” nature leading to the name of the model.

When $`\epsilon =0`$, one can determine the stable and unstable manifolds of the fixed point $`\omega ^{}`$:

(5.12)
$$\omega _p=\mathrm{\Gamma },\omega _n=0(n).$$

Using the polar coordinates:

$$\omega _1=r\mathrm{cos}\theta ,\omega _4=r\mathrm{sin}\theta ;\omega _2=\rho \mathrm{cos}\vartheta ,\omega _3=\rho \mathrm{sin}\vartheta ;$$

we have the following explicit expressions for the stable and unstable manifolds of the fixed point $`\omega ^{}`$ and its negative $`\omega ^{}`$ represented through perverted heteroclinic orbits:

$`\omega _p`$ $`=`$ $`\mathrm{\Gamma }\mathrm{tanh}\tau ,`$
$`r`$ $`=`$ $`\sqrt{{\displaystyle \frac{A_2}{A_2A_1}}}\mathrm{\Gamma }\text{sech}\tau ,`$
(5.13) $`\theta `$ $`=`$ $`{\displaystyle \frac{A_2}{2\kappa }}\text{ln}\mathrm{cosh}\tau +\theta _0,`$
$`\rho `$ $`=`$ $`\sqrt{{\displaystyle \frac{A_1}{A_2}}}r,`$
(5.16) $`\theta +\vartheta `$ $`=`$ $`\{\begin{array}{cc}\mathrm{arcsin}\left[\frac{1}{2}\sqrt{\frac{A_2}{A_1}}\right],\hfill & (\kappa >0),\hfill \\ \pi +\mathrm{arcsin}\left[\frac{1}{2}\sqrt{\frac{A_2}{A_1}}\right],\hfill & (\kappa <0),\hfill \end{array}`$
where $`\tau =\kappa \mathrm{\Gamma }t+\tau _0`$, $`(\tau _0,\theta _0)`$ are the two parameters parametrizing the two-dimensional stable (unstable) manifold, and

$$\kappa =\sqrt{A_1A_2}\mathrm{cos}(\theta +\vartheta )=\pm \sqrt{A_1A_2}\sqrt{1+\frac{A_2}{4A_1}}.$$

The two auxilliary variables $`\omega _0`$ and $`\omega _5`$ have the expressions:

$`\omega _0`$ $`=`$ $`{\displaystyle \frac{\alpha \beta }{1+\beta ^2}}\text{sech}\tau \left\{\mathrm{sin}[\beta \text{ln}\mathrm{cosh}\tau +\theta _0]{\displaystyle \frac{1}{\beta }}\mathrm{cos}[\beta \text{ln}\mathrm{cosh}\tau +\theta _0]\right\},`$
$`\omega _5`$ $`=`$ $`{\displaystyle \frac{\alpha \beta }{1+\beta ^2}}\text{sech}\tau \left\{\mathrm{cos}[\beta \text{ln}\mathrm{cosh}\tau +\theta _0]+{\displaystyle \frac{1}{\beta }}\mathrm{sin}[\beta \text{ln}\mathrm{cosh}\tau +\theta _0]\right\},`$
where

$$\alpha =A_1\mathrm{\Gamma }\kappa ^1\sqrt{\frac{A_2}{A_2A_1}},\beta =\frac{A_2}{2\kappa }.$$

The explicit expression (5.13) of the perverted heteroclinic orbits shows that the unstable manifold $`W^u(\omega ^{})`$ of the fixed point $`\omega ^{}`$ is the same as the stable manifold $`W^s(\omega ^{})`$ of the negative $`\omega ^{}`$ of $`\omega ^{}`$, and the stable manifold of the fixed point $`\omega ^{}`$ is the same as the unstable manifold of the negative $`\omega ^{}`$ of $`\omega ^{}`$. Both $`W^u(\omega ^{})`$ and $`W^s(\omega ^{})`$ have the shapes of “painted eggs” (Figure 13). $`W^u(\omega ^{})`$ and $`W^s(\omega ^{})`$ together form the “lip” (Figure 13) which is a higher dimensional generalization of the heteroclinic connection on plane.

As the first step toward understanding the degeneracy v.s. nondegeneracy of the hyperbolic structure of 2D Euler equation, we are interested in the $`\epsilon `$-homotopy deformation of such hyperbolic structure for the dashed-line model. We have calculated the Melnikov functions to study the “breaking” or “persistence” of such hyperbolic structure when $`\epsilon `$ is small. It turns out that both the first and the second order Melnikov functions are identically zero .

## 6. Conclusion

Results on chaos in partial differential equations are reported. Results on Lax pairs of Euler equations of incompressible inviscid fluids are also reported. Preliminary results on dynamical system studies of 2D Euler equation are summarized.


format: text

text:
# CMB power spectrum estimation and map reconstruction with the Expectation - Maximization algorithm

## 1 Introduction

Undergoing CMB experiments like BOOMERanG (Netterfield et al. 2002, Rulh et al. 2002), MAXIMA (Hanany et al. 2000), DASI (Halverson et al. 2002), VSA (Grainge et al. 2002), CBI (Mason et al. 2002), ACBAR (Kuo et al. 2002), Archeops (Benoit et al. 2002) and MAP as well as future ones (Planck) will reveal with unprecedented quality the primordial matter fluctuations responsible for the CMB. Very recently the first detection of the E-mode polarization in the CMB has been claimed (DASI, Kovac et al. 2002), which independently supports the structure formation models via gravitational instability. Thanks to these experimental results the fundamental cosmological parameters can be determined with good accuracy.

These experiments will measure also the emission coming from our own Galaxy at mm frequencies as well as the emission due to other galaxies in the same wavebands. On the other hand, galaxy clusters will distort the CMB radiation with an intensity which is proportional to the total mass of the cluster times its temperature. All these components will generate a confusion limit which, together with the intrinsic detector noise, will make very difficult to disentangle which percentage of the observed emission is due to the CMB or to the other components. Several component separation methods have been proposed so far in the literature: multifrequency Wiener filter (MWF, Tegmark and Efstathiou 1996, Bouchet and Gispert 1999), maximum entropy methods (MEM, Hobson et al. 1998, Vielva et al. 2001b; extended to the sphere in Stolyarov et al. 2002, Barreiro et al. 2003), independent component analysis (ICA, Baccigalupi et al. 2000, Maino et al. 2002), blind Bayes (Snoussi et al. 2001).
Although there are some similarities and differences between the previous methods, all of them share a common aspect: they try to separate all the components simultaneously. For this purpose, these methods usually need to assume some a priori information about the components. Thus, if the different components have different frequency dependencies, then, by examining the data at different frequencies it is possible to distinguish (at a certain level) the different components. Or, if we know the correlation function of the components (or their power spectrum) and each one is significantly different from each other, then, it is also possible to distinguish (again at a certain level) the components since their modes will behave differently in the Fourier space. If one knows both, the frequency dependence and the power spectrum, then the component separation improves dramatically since now it is possible to combine the information in the different channels by correlating the Fourier modes with a correlation given essentially by the frequency dependence of each component and its power spectrum. This approach has proved to be very useful if the frequency dependence and power spectrum of the components are known.
It is interesting to explore what kind of information can be obtained when no prior information is assumed about none of the components. In some cases we simply do not know the priors (as it happens in the case of the free-free emission or the spinning dust). On the other hand, if we assume something about the components which is not accurate but an approach, then we are intrinsically introducing some bias in our final result. If, after the component separation, we end up with a biased estimate on, at least, one of the components, this will have an effect on some (if not all) of the other components which must absorb that bias in order to obey the constraint that the sum of all the components must be equal to the original data. This very last point is one of the risks of the simultaneous component separation methods.
Another assumption usually made is that the components are, one from each other, statistically independent. This assumption is not true for several of the components (like the Galactic ones). This assumption can be a source of systematic errors. Other typical assumption is that the probability distribution function (pdf) of the individual components is a Gaussian. This is false for components like the point sources, the Sunyaev-Zeldovich (SZ) and the Galactic emissions where the true pdf has a bell-shape with a long tail towards positive values (negative for the SZ for frequencies below 217 GHz). This assumtion will be somewhat relaxed in this work by assuming that the combined foreground emission plus instrumental noise is Gaussian distributed in Fourier space, instead of assuming it for each individual component. We will study the effect of this assumption at different scales. At those scales where the assumption is less justified solutions will be proposed to reduce the effect.

In a recent paper (Delabrouille et al. 2002) the authors have used the Expectation-Maximization (EM) algorithm to minimize a spectral matching criterion. They consider the problem of component separation with four components, CMB, SZ effect, dust and instrumental noise. In that work, the authors have shown that this is an interesting alternative to estimate jointly the frequency dependence and spatial power spectra of the components. This method allows to introduce certain degrees of freedom or free parameters in the problem. These free parameters can be determined after iterating an expectation-maximization process.

In the work presented here, we will explore this direction but in a much more simplified manner focusing only on the estimation of the power spectrum and the map reconstruction of the CMB but including more components than in the previous work and doing no assumptions at all about the power spectrum or frequency dependence of any of the components. In addition to the diffuse foregrounds, we also will study the effect of point sources which were not considered in the previous work.

The approach presented in this paper can be extended to determine jointly the CMB and the SZ components. The physics of these two emissions is well known, in particular their frequency dependence. Therefore, no prior assumptions about their frequency behaviour is required. In addition, these two components –together with the point sources one– are the most important signals from the cosmological point of view. Several works have already been presented to determine the SZ (Herranz et al. 2001b,c; Diego et al. 2002) and the point sources emissions (Tegmark & Oliveira-Costa 1998; Cayón et al. 2000; Sanz et al. 2001; Vielva et al. 2001, 2003; Chiang et al. 2002; Hobson & McLachlan 2002) from microwave images. Our aim for the future will be to develop a method based on the EM algorithm combined with a point source detection technique, to recover simultaneously the three cosmological emissions.

The paper is organized as follows. In Section 2, we summarize the key ideas behind the EM algorithm. In Section 3 we apply the EM algorithm to the problem of determining an estimator of the CMB power spectrum. The main results are shown in Section 4. These results are compared with the ones which would be obtained in the ideal case when only CMB and white noise are considered. Finally, the conclusions are given in Section 5.

## 2 The Expectation-Maximization algorithm

In this section we will present an alternative method which will be useful to get a robust estimate of the power spectrum of the CMB. We will also present a single iterative expression for the new estimator of the CMB power spectrum. This estimator can be used directly with multifrequency microwave data to give a fast, accurate and robust estimate of the power spectrum of the CMB.
EM is an algorithm useful when there is a many-to-one mapping. That is, several variables are combined together to give one observed quantity. This is exactly our problem where the complete data are just the different components (Galactics, extragalactics, CMB and noise) and the observed quantity is just the projection of all of them in one of the channels (data on the receiver). EM allows to introduce degrees of freedom (or free parameters) in the pdf of the complete data and then look for the best parameters by maximising the expected value of that pdf given the observations.
The advantage of this algorithm in our problem is obvious. Since part of the information about the components is unknown, we can parametrise this unknown information as free parameters in the pdf of the complete data. Then, the maximization process will give us the best set of free parameters given the data.

The EM algorithm (Dempster et al. 1977) provides an iterative procedure for computing maximum likelihood estimates (MLEs) in situations where the observed data vector, $`𝐝`$, is viewed as being incomplete. $`𝐝`$ is an observable function of the so-called complete data $`𝐱`$ (see e.g. McLachlan and Krishnan 1997). Let’s denote by $`𝐟_𝐢(𝐝|𝐩)`$ the pdf’s of the incomplete data $`𝐝`$ and by $`𝐟(𝐱|𝐩)`$ the pdf of $`𝐱`$. The vector $`𝐩=(𝐩_\mathrm{𝟏},\mathrm{}.,𝐩_𝐍)`$ will contain the $`N`$ unknown parameters. The complete-data log likelihood function that could be formed for $`𝐩`$ if $`𝐱`$ were fully observable is given by

$$\mathrm{log}L(𝐩)=\mathrm{log}𝐟(𝐱|𝐩).$$
(1)
Instead of observing $`𝐱`$, we observe $`𝐝`$, with many-to-one mapping from $`𝐱`$ to $`𝐝`$ (i.e $`𝐝=𝐝(𝐱)`$). It follows that

$$𝐟_𝐢(𝐝|𝐩)=_{𝐗_𝐝}𝐟(𝐱|𝐩)\mathrm{𝐝𝐱}$$
(2)
where $`𝐗_𝐝`$ is the sub-space of $`𝐱`$ defined by $`𝐝=𝐝(𝐱)`$.
The EM algorithm approaches the problem of solving the incomplete-data likelihood, $`\mathrm{log}(𝐩)=\mathrm{log}𝐟_𝐢(𝐝|𝐩)`$, indirectly by proceeding iteratively in terms of the complete-data log likelihood function, $`logL(𝐩)`$. As it is unobservable, it is replaced by its conditional expectation given $`𝐝`$, using the current fit for $`𝐩`$. More specifically, the EM algorithm consists of two steps: Expectation (E-step) and Maximization (M-step). On the j+1 iteration, these steps are as follows:

$``$ E-step: Calculate the quantity $`Q(𝐩|𝐩^𝐣)`$, defined as $`Q(𝐩|𝐩^𝐣)=E\{\mathrm{log}L(𝐩)|𝐝,𝐩^𝐣\}`$.

$``$ M-step: Choose $`𝐩^{𝐣+\mathrm{𝟏}}`$ to be the value that maximises $`Q(𝐩|𝐩^𝐣)`$; that is, $`Q(𝐩^{𝐣+\mathrm{𝟏}}|𝐩^𝐣)=\mathrm{max}[Q(𝐩|𝐩^𝐣)]`$.

That is, in the first step we compute the expected value of the log of the complete pdf given the data and an estimate of the parameters ($`𝐩^𝐣`$). In the second step we look for the values of the parameters which maximise the previous expected value. At this step, we have to maximise with respect to each one of the parameters ($`p^j(i)`$). This can be done by just setting the first derivative of $`Q(𝐩;𝐩^𝐣)`$ with respect to that parameter to $`0`$ and solving for $`p^j(i)`$.
This is an iterative process which can be started with an arbitrary choice for the parameters in the first step, $`𝐩^𝐨`$. The EM algorithm assure us that at every step the likelihood $`(𝐩)`$ is increased. The method will converge to the optimal set of parameters, $`𝐩`$, in a number of steps which will depend on the nature of the problem and on the initial election for the first iteration, $`𝐩^𝐨`$.

In this paper, we will consider a simple case where we allow the covariance in Fourier space of the CMB signal (or power spectrum) to be a free parameter.
We will apply the EM algorithm in order to determine that free parameter (the power spectrum of the signal). This case can be extended to include more free parameters in the analysis as it is shown in Delabrouille et al. (2002).

## 3 EM applied to CMB experiments

If we know the frequency dependence of one of the components, then we can express the signal at a given frequency as:

$$𝐒(𝐧,\nu )=𝐀(𝐧,\nu )\left(𝐟(\nu )s(𝐧)\right)$$
(3)
where $`s(𝐧)`$ is the spatial pattern of the signal we want to estimate, $`𝐟(\nu )`$ is the frequency dependence of the signal (including the band-with), and $`𝐀(𝐧,\nu )`$ accounts for the convolution with the antenna beam of the experiment and its frequency response.
The data at different frequencies, $`𝐝(\nu )`$, can be expressed as a sum of the signal plus some other contributions.

$$𝐝(𝐧,\nu )=𝐒(𝐧,\nu )+\xi (𝐧,\nu ).$$
(4)
The residual $`\xi (𝐧,\nu )`$ includes all the other components, i.e. Galactic and extragalactic components convolved with the corresponding beam plus the corresponding instrumental noise. Due to the antenna and frequency convolution in $`𝐀(𝐧,\nu )`$, it is easier to work in Fourier space where this convolution is just a product of the Fourier modes which are uncorrelated provided the field is homogeneous and isotropic. Therefore, the previous equation should be rather expressed in Fourier space as:

$$d_\nu (𝐤)=A_\nu (k)f_\nu s(𝐤)+\xi _\nu (𝐤).$$
(5)
If the experiment has $`m`$ different channels, then we have an observation for each channel. That is, we have the vector of observations $`𝐝=(d_1,d_2,\mathrm{},d_m)`$ for each Fourier mode $`𝐤`$. And we can write:

$$𝐝=𝐑s+\xi $$
(6)
where, $`𝐑`$ is the response vector containing $`m`$ elements. Each element of $`𝐑`$ is just the product of the antenna and the frequency response of the instrument times the frequency dependence of the signal,

$$R_\nu =A_\nu f_\nu .$$
(7)
In the case in which the signal is the CMB then $`f_\nu =1\nu `$. From eqn. 6 we can express the residual as

$$\xi =𝐝𝐑s.$$
(8)
The basic starting point when applying EM is to define the pdf of the complete data. In our case the complete data are all the components plus the multifrequency observations $`𝐝`$. We will focus on the problem in which the complete data can be divided in just two elements. The two elements are the CMB signal we want to estimate $`s`$ and the observed data $`𝐝`$ (or equivatently the residual $`\xi `$, noise plus rest of components; see below). By doing this, the nature of our problem reduces drastically since we only have to make assumptions about two elements.
Furthermore, the residual and the signal can be considered as independent. In terms of the probability, the pdf of the complete data is:

$$𝐟(𝐱|𝐩)=𝐏(s,𝐝|𝐩)=𝐏(s,\xi |𝐩)=𝐏(s|𝐩_𝐀)𝐏(\xi |𝐩_𝐁)$$
(9)
where $`𝐩=𝐩_𝐀𝐩_𝐁`$ with $`𝐩_𝐀`$ and $`𝐩_𝐁`$ unknown parameter sets for the pdf’s of $`s`$ and $`\xi `$. The second equality follows from the transformation $`(s,𝐝)(s,\xi )`$ whose Jacobian is equal to one (see equation 5). It is important to note that, in oppossition to other methods, our assumption about independence of the elements ($`s`$ and $`\xi `$) is well justified. No relation is expected between the CMB and the other foregrounds (except maybe the SZ effect which could be weakly related to the CMB through the ISW effect).

Once we have established that the two elements are independent, we need a clue about the specific form of the individual probabilities, $`𝐏(s|𝐩_𝐀)`$ and $`𝐏(\xi |𝐩_𝐁)`$. If the CMB is close to a Gaussian variable (and we know it is), then the pdf of the CMB in Fourier space can be approached by:

$$𝐏(s|𝐩_𝐀)𝐏(s)P(k)^1exp(\frac{s^2}{P(k)})$$
(10)
where $`P(k)`$ is the power spectrum of the CMB: the parameter to be estimated.
The pdf of the residual ($`\xi `$) can be modelled as

$$𝐏(\xi |𝐩_𝐛)𝐏(\xi )exp(\xi C^1\xi ^{})$$
(11)
where $`C^1`$ is the inverse of the correlation matrix of the residual. $`C`$ is an $`n\times n`$ matrix ($`n=`$ number of channels). This correlation matrix is not diagonal if there are some correlations in the residuals between the different frequency channels. The elements of $`C`$ are the power (diagonal) and cross-power (off-diagonal) spectra of the residuals at each channel which can be given as a function of the $`P(k)`$ (see below).
We have assumed in equation (11) that the residuals are Gaussian distributed. This assumption is not as strong as assuming that each component is Gaussian distributed. The sum of several non-Gaussian distributions tends to a more Gaussian one. In the limit of a sum of infinite independent components the central limit theorem assures that the resulting distribution is Gaussian. In any case Gaussianity can be a good approach at scales where the instrumental noise dominates the residuals. The most critical situations for the validity of this assumption appear at large and small scales since, at least for some channels, they are dominated by the Galactic components and the compact sources, respectively. Solutions for these problems will be proposed in the next sections.

In the previous relations for the individual pdf’s appear the power spectra of the signal and the residual. In a real situation these quantities are unknown. This is the reason why we apply EM. We can introduce the unknown CMB power spectrum, $`P(k)`$, as the uncertainty in our problem. In order to apply EM and estimate this free parameter we have to calculate the expected value $`Q(𝐩|𝐩^𝐣)`$ which is just an integral:

$$Q(𝐩|𝐩^𝐣)=log𝐏(s,𝐝|𝐩)𝐏(s|𝐝,𝐩^𝐣)𝑑s.$$
(12)
The only thing we need to know to calculate $`Q(𝐩|𝐩^𝐣)`$ are the probabilities $`𝐏(s|𝐝)`$ and $`𝐏(s,𝐝)`$. By Bayes we know that:

$$𝐏(s|𝐝)=\frac{𝐏(s,𝐝)}{𝐏(𝐝)}$$
(13)
where $`𝐏(𝐝)`$ can be obtained just by marginalising $`𝐏(s,𝐝)`$ over $`s`$:

$$𝐏(𝐝)=𝐏(s,𝐝)𝑑s$$
(14)
and $`𝐏(s,𝐝)`$ can be obtained through the chain rule:

$$𝐏(s,𝐝)=𝐏(s,\xi )J\left(\frac{s,\xi }{s,𝐝}\right).$$
(15)
By taking equation (5), one can prove that the Jacobian in the previous expression is equal to one, as it was commented above, so we have:

$$𝐏(s,𝐝)=𝐏(s,\xi )$$
(16)
where (see equation 9),

$$𝐏(s,\xi )P(k)^1exp(\frac{s^2}{P(k)})exp(\xi C^1\xi ^{})$$
(17)
where the term, $`P(k)`$, accounts for our free parameter in $`𝐏(s|𝐩_𝐀)`$ (see equation 10) and the second term was taken from equation (11). We remark that the covariance matrix $`C`$ also depends on the free parameter $`P(k)`$, as it will be discussed below.

Now it is easy to compute the terms, $`𝐏(𝐝)`$, $`𝐏(s|𝐝)`$, and finally $`Q(P(k)|P(k)^j)`$:

$$𝐏(𝐝)\frac{exp\left(𝐝C^1𝐝^{}+(𝐝C^1𝐑^{})^2/F_1\right)}{P(k)|C|^{1/2}F_1}$$
(18)
where we have expressed $`\xi `$ in terms of $`𝐝`$ and $`s`$ (see eqn. 8).

$$𝐏(s|𝐝)F_1exp\left(F_1\left|s\frac{𝐝C^1𝐑^{}}{F_1}\right|^2\right)$$
(19)
where,

$$F_1=𝐑C^1𝐑^{}+P(k)^1.$$
(20)
Finally the expected value of the log of the complete pdf (equation 12) is:

$`Q(P(k)|P(k)^j)`$ $`=`$ $`log(P(k)){\displaystyle \frac{1}{2}}log|C|𝐝C^1𝐝^{}`$
$``$ $`F_1<|s|^2>+2𝐝C^1𝐑^{}<s>+const.`$
where :

$$<s>=\frac{𝐝C^1𝐑^{}}{F_1}$$
(22)
$$<|s|^2>=\frac{1}{F_1}\left(1+\frac{|𝐝C^1𝐑^{}|^2}{F_1}\right).$$
(23)
Equation (22) is nothing more than multifrequency Wiener filter computed with the CMB power spectrum at iteration $`j`$. This is an expected result since Gaussianity for the signal and residuals have been assumed. Similarly, the signal variance at each $`k`$ $`<|s|^2>`$ is also computed with $`P(k)^j`$ at iteration $`j`$. At this point it is important to note that in eq. (3) the dependence on $`P(k)`$ appears explicitly in the first term and also in the fourth one through the function $`F_1`$. Moreover, since we do not know the residuals the $`P(k)`$ dependence is also present in the other terms through the residual covariance matrix,

$$C^{\nu \mu }=P_\xi ^{\nu \mu }(k)=P_𝐝^{\nu \mu }(k)R^\nu R^\mu P(k)$$
(24)
where $`P_\xi ^{\nu \mu }(k)`$ and $`P_𝐝^{\nu \mu }(k)`$ are the residual and data cross-power spectra after convolution with the corresponding antenna response whereas $`P(k)`$ is the unconvolved CMB power spectrum to be determined.

In the maximisation step, the value of $`P(k)^{j+1}`$ that maximises $`Q(P(k)|P(k)^j)`$ should be found. Finding the maximum is not a trivial task given the complicated dependence on $`P(k)`$. In order to facilitate the finding we consider the residual covariance matrix, $`C^{\nu \mu }`$, constant in the M-step (but changing with $`P(k)`$ in the E-step). With this simplification, we can easily find the analitical solution for the maximum by differentiating $`Q(P(k)|P(k)^j)`$ with respect to $`P(k)`$ and equating to zero. We find the not so surprising result:

$$P(k)^{j+1}=<|s|^2>.$$
(25)
The previous equation is the main result of EM. $`<|s|^2>`$ is given by equation (23) computed with the value of $`P(k)`$ at iteration $`j`$. By iterating equation (25), which only depends on the data, the response vector and the CMB power spectrum obtained in the previous iteration, we can get an estimate of the power spectrum of the underlying signal. The final step after the EM estimation of the signal power spectrum at each $`k`$, is to recover the CMB map by applying eqn. 22, the multifrequency Wiener filter derived within the EM framework.

Finally, this method can be easily extended to determine jointly the CMB and the SZ emissions, since they are almost uncorrelated and the frequency dependence of the SZ is also very well known. In addition, the method can be implemented with more realistic (non-Gaussian) pdfs to model the signal distributions as well as the residuals. We will study these extensions of the method in a future work.

## 4 Results

In this section, we present the results that have been obtained for the CMB power spectrum determination. Firstly, we briefly describe the simulated data that have been used in this work, accounting for the most important contributions to the microwave sky. Then, we present our CMB power spectrum estimation. An improvement at small scales is achived by using a very useful tool for subtracting the brightest point sources: the Mexican Hat Wavelet.

Finally, we applied our estimator to an ideal data set, where the Gaussian assumption for the residuals is satisfied. This exercise is highly interesting, since it shows the robustness of our estimator: the results achieved in this idealistic case are only slightly better than the ones obtained in the reslistic situation.

### 4.1 Simulated data

To show the power of our approach we will apply our estimator (eqn. 25) to realistic simulations. They consider the expected levels of the instrumental white noise and resolution characteristics of the 10 Planck mission channels at 30, 44, 70, 100 (Low and High Frequency Instrument, LFI and HFI), 143, 217, 353, 545 ad 857 GHz (see Table 1). In addition, we use state of the art simulations of the different galactic and extragalactic components, together with the pure CMB signal.

Within the first ones, we have taken into account the synchrotron, free-free, spinning and thermal dust contributions. The first one has been simulated using the all-sky pattern provided by Giardino et al. (2002) including both, temperature and spectral index templates. The free-free emission is very poorly known. Recent experiments focusing on the H-$`\alpha `$ detection, like Southern H-$`\alpha `$ Sky Survey (SHASSA, Reynolds & Haffner, 2000) and the Wisconsin H-$`\alpha `$ Mapper project (WHAM, Gaustad et al. 2001), will produce all-sky surveys that could be used as free-free templates. However, since these data are not available at the present time, we have used the correlation between dust and free-free emissions proposed by Bouchet et al. (1996) as spatial template, and a power law, $`I_\nu \nu ^{0.16}`$, to describe the frequency dependence. We have also included the spinning dust emission –proposed by Draine & Lazarian (1998) as a possible explanation for the anomalous Galactic emission found by CMB experiments like COBE (Kogut 1999) or Saskatoon (Oliveira-Costa et al. 1997). Very recently, a tentative confirmation of that emission has been claimed (Finkbeiner et al. 2002). The simulation of that emission has been carried out using the frequency dependence proposed by Draine & Lazarian (1998)<sup>1</sup><sup>1</sup>1Data provided by the authors. and using the thermal dust component as an spatial template, since both dust emissions are strongly correlated through the neutral hydrogen column density ($`N_H`$). Finally, thermal dust emission has been also simulated, using the best model proposed by Finkbeiner et al. (1999) to fit the FIRAS, IRAS and DIRBE data, consisting of two grey-bodies with mean temperatures of 16.2 K and 9.4 K and emissivities 2.70 and 1.67, respectively.

The simulated extragalactic foregrounds are the thermal Sunyaev-Zel’dovich effect (SZ) and radio and infrared point sources. The SZ was performed following the Diego et al. (2001) model for a flat $`\mathrm{\Lambda }`$CDM Universe with $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = $`0.7`$. The point sources correspond to the Toffolatti et al. (1998) model for the same Universe, including radio flat-spectrum and infrared sources. We refer to that paper for details.

Finally, the CMB signal was simulated using the $`C_l`$’s provided by the CMBFAST code (Seljak & Zaldarriaga, 1996) for the same cosmological model and assuming Gaussian fluctuations. We would like to remark that synchrotron, thermal dust and the point source emissions have been simulated using a frequency dependence which varies with the sky position. The simulated maps covering a patch of the sky of $`12.8^o\times 12.8^o`$, as seen by the 10 Planck channels, are shown in figure 1.

### 4.2 A simple choice for the correlation matrix

In order to estimate the power spectrum of the CMB through eqn. (25) we still have to iterate a relatively complex expression containing a cocient of polynomials of order twice the number of channels and involving many terms. To further simplify the convergence of that equation we will only consider the diagonal terms of the covariance matrix, making zero the other non-diagonal elements. This assumption does not take into account the correlations among the residuals at different channels due to the Galactic components and compact sources. The former are reflected at the low $`k`$-modes and their effect can be reduced by performing a simple Galactic subtraction as the one proposed by Diego et al. (2002). The later have their impact at intermediate (e.g. radio sources observed in the lowest resolution channels) and large $`k`$-modes (e.g. infrared sources). However, as it will be shown below, once the brightest point sources are subtracted from the maps, the contribution of the remanent point sources is low enough to make their effect significantly smaller.

In figure 2 we present our results for the 10 Planck channels using the data obtained from the combination of the simulated components as described above (figure 1). The comparison of the recovered spectrum with the true one shows a good agreement until $`k`$ values larger than the one corresponding to the best antenna channels and the results start to be dominated by the deconvolved noise.

In order to better appreciate the differences, we show in figure 3 the relative error, $`(P_r(k)P_t(k))/P_t(k)`$ %, where $`P_r(k)`$ is the estimator of equation 25.

In both large and relatively small scale range (small and relatively large $`k`$) our estimator finds a reasonably good estimation of the true power spectrum. The error is small (less than 10 %) up to scales $`k55`$ (which corresponds to $`l1500`$). The effect of only considering the diagonal terms of the residual covariance matrix is small (see also the comparison with the ideal case below). This fact proves the applicability of our assumption for the microwave data since our method is not very sensitive to the cross-correlations of the different channels.

At very small scales ($`k>45`$ or $`l>1300`$) the signal is below the noise level and the error bars grow quickly. Our estimator shows a clear bias toward higher values (we recover more power than the true one). This implies that some of our assumptions are wrong at these small scales. The Gaussianity assumption is a good approximation in Fourier space at large $`k`$-modes (because of the dilution in many modes of the clear non-Gaussian features appearing in the real space). However, the uncorrelation assumption among frequency channels is a source of bias at these scales, since the bright point sources are strongly correlated. We have checked this point by removing the brightest point sources (but leaving the galaxy clusters and the weak point sources in the residual). We have used the Mexican Hat Wavelet (MHW) technique proposed by Vielva et al. (2001a) for removing point sources above the fluxes determined by the so-called $`50\%`$ error criterion used in that paper. These fluxes correspond to those detection limits above which the maximum percentage of spurious detection is $`5\%`$, being a detection spurious if the error in the amplitude estimation is larger than $`50\%`$. After the subtraction, we applied our estimator to these new point-source-free maps. The final result is shown in figure 4. The effect of removing the brightest point sources is evident now. This result clearly shows how the bias can be corrected by removing the compact sources. The relative error does not change significantly after removing the point sources but the estimator is able to go up to scales $`k=60`$ ($`l1700`$) with relative error $`<10\%`$. The error is smaller than $`50\%`$ below $`k=75`$ ($`l2100`$).

### 4.3 Comparing with an ideal situation

It is interesting to answer whether or not the above results are close to the most optimistic case where all our assumptions are right (Gaussianity and $`C`$ diagonal). We will check this point in the simple (but unrealistic) case where our data consist only of CMB signal plus the corresponding instrumental noise for each channel. In this simple case, our assumption of Gaussianity is right since both CMB and noise are Gaussian. The correlation matrix of the residual, $`C`$, will be exactly diagonal with the elements being the power spectrum of the data minus the CMB power spectrum convolved with the antenna corresponding to each channel.

In figures 3 and 4 we show the performance of the EM estimator in this ideal case (CMB plus noise) after combining the 10 Planck channels. The result is surprisingly very similar to the one obtained in the realistic case where the brightest point sources have been subtracted (figure 4). This test shows how the EM estimator is extremely robust and can give results very close to the most optimal one.

### 4.4 Map reconstruction

Using the expression given in equation 22 and the estimated power spectrum $`P_r(k)`$, we are able to recover the CMB map. As it was pointed out in Section 3, due to the Gaussianity assumptions, equation 22 is nothing more than multifrequency Wiener filter (MWF). Hence, our map reconstruction is a MWF focused on the CMB recovery. Let us remark that this approach differs from the traditional MWF already used in the microwave sky recovery (e.g., Tegmark & Efstathiou 1996, Bouchet & Gispert 1999 and Hobson et al 1998). These works were focused on the all components separation, whereas our approach is devoted to reconstruct just the CMB signal. Even more, whereas the traditional approaches require previous knowledge of the power spectra (not only the CMB one, but also the power spectra of the foregrounds), the CMB power spectrum is a result of the blind EM method. This and the non requirement of any prior knowledge about the frequency dependence of the components, are the strongest points of the blind methods. In figure 5 we present the CMB map reconstruction (left) together with the input one (right). Both maps have been convolved with the best Planck resolution (FWHM = 5’). The correlation between the input and recovered CMB maps is very good $`r=0.98`$ (figure 6), whereas the slope of the best straightline fit is 0.97. The residual CMB is dominated by the Galactic dust emission, since it is the most important foreground at low $`k`$-modes. Hence, by performing a simple dust subtraction as the one proposed by Diego et al. (2002) –which basically consists in using the 857 GHz channels as a dust template to be subtracted from the others channels– a significant improvement in the CMB map reconstruction is achieved, being the RMS of the residual map lower than $`15\%`$.

## 5 Conclusion

We have presented an efficient estimator (eqn. 25) of the CMB power spectrum. Our estimator is based on the EM algorithm and does not make use of any prior information whatsoever about any of the components. This renders satisfactory results, since we are able to recover the power spectrum up to $`l1500`$ with less than 10% error. This limit can be improved if we remove the brightest compact sources, for instance applying the MHW technique poposed by Vielva et al. (2001a). After point source removal the estimator can recover the power spectrum up to $`l1700`$ with less than 10% error and up to $`l2100`$ with less than 50% error.

Our results are very close to the optimal case when all our assumptions are satisfied (only CMB and instrumental noise are included in the simulations), as can be seen in figure 4. There are several CMB power spectrum estimations in the literature, given for different component separation techniques, like Bayesian methods (MEM: Hobson et al. 1998 and Stolyarov et al. 2002; and MWF: Tegmark & Efstathiou 1996 and Bouchet & Gispert 1999) and blind source separation algorithms (ICA and FastICA: Baccigalupi et al. 2000 and Maino et al. 2002; and Multi-Detector Multi-Component analysis: Delabrouille et al. 2002). This work provides a CMB power spectrum that is comparable to or even better than the previous ones, being more robust than those, since previous knowledge about the underlying components is not requiered (contrary to the Bayesian methods) and simulations used in previous works are significantly more idealised than the ones used in this work (e.g., lack of spatial variation of the frequency dependence for several components; absence of some of the components, like point sources or some Galactic foregrounds). On the other hand, our simulations account for the previous limitations including all the major microwave emissions and allowing for spatial variation of the frequency dependence. Finally, our method is able to reconstruct the CMB map using the estimated power spectrum and the MWF given in equation 22. Let us remark that the MWF is obtained when the Gaussianity hypothesis is assumed, whereas the EM algorithm provides a more general procedure for component separation. That hypothesis will be explored in a future work by including more realistic distributions for the different components. In addition, we are working on a jointly determination of the cosmological signals (CMB, SZ and point sources) using the EM algorithm and a point source detection tool. Finally, an extension of this algorithm to the sphere is straightforward and will be studied soon.

## Acknowledgements

We thank J. Delabrouille for a critical reading of a previous version of the manuscript which resulted in a substantial improvement. We thank L. Toffolatti for providing us with the point source simulation, and B. T. Draine and A. Lazarian for providing us with the emissivity predicted by their spinning dust model. We appreciate interesting discussions with J. L. Sanz, R. B. Barreiro, D. Herranz and V. Stolyarov. We thank the RTN of the EU project HPRN-CT-2000-00124 for partial finantial support. EMG and PV acknowledge partial financial support from the Spanish MCYT projects ESP2001-4542-PE and ESP2002-04141-C03-01. This research has been supported by a Marie Curie Fellowship of the European Community programme Improving the Human Research Potential and Socio-Economic knowledge under contract number HPMF-CT-2000-00967. PV acknowledges support from a fellowship of Universidad de Cantabria.