text:
# Dynamic Point-Formation in Dielectric Fluids

## Abstract

We use boundary-integral methods to compute the time-dependent deformation of a drop of dielectric fluid immersed in another dielectric fluid in a uniform electric field $`E`$. Steady state theory predicts, when the permittivity ratio, $`\beta `$, is large enough, a conical interface can exist at two cone angles, with $`\theta _<(\beta )`$ stable and $`\theta _>(\beta )`$ unstable. Our numerical evidence instead shows a dynamical process which produces a cone-formation and a transient finite-time singularity, when $`E`$ and $`\beta `$ are above their critical values. Based on a scaling analysis of the electric stress and the fluid motion, we are able to apply approximate boundary conditions to compute the evolution of the tip region. We find in our non-equilibrium case where the electric stress is substantially larger than the surface tension, the ratio of the electric stress to the surface tension in the newly-grown cone region can converge to a $`\beta `$ dependent value, $`\alpha _c(\beta )>1`$, while the cone angle converges to $`\theta _<(\beta )`$. This new dynamical solution is self-similar.

preprint: APS/123-QED

The formation of conical ends on fluid-fluid interfaces in strong electric/magnetic fields has been seen in various electrospraying and ferrofluid experiments zeleny1917 ; taylor1964 ; garton1964 ; fernandez1992 ; nagel2000 ; bacri1982b . Building on the work of Taylor taylor1964 , Li et al. and Ramos & Castellanos studied the electrostatics of an infinite cone with semi-vertical angle $`\theta _0`$ formed between two dielectric fluids with permittivity ratio $`\beta `$. In spherical coordinates, the electric stress $`\sigma _er^{2(\nu 1)}`$. In an equilibrium cone, this stress must be balanced by the surface tension, so that $`\nu `$ must be $`1/2`$. According to their analysis li1994 ; ramos1994 , there are two such solutions of $`\theta _0`$, $`\theta _<(\beta )`$ and $`\theta _>(\beta )`$, which will occur for $`\beta >\beta _c=17.59`$. The former is said to be stable, the latter unstable li1994 .

In contrast, this letter describes a dynamical fixed point in which a cone is formed transiently. At the fixed point, the cone angle is $`\theta _<(\beta )`$ so that surface stress and electric stress have the same scaling in the cone, with the ratio of the two being constant, but different from unity. The electric stress is the larger of the two, so that the total surface stress always acts to elongate the pointed region.

We compute the time-dependent deformation of a dielectric drop (fluid 1) freely suspended in another dielectric fluid (fluid 2) in a uniform electric field. Both fluids are incompressible and have the same viscosity $`\eta `$. There is surface tension with coefficient $`\gamma `$ between the two fluids. The drop is axially symmetric and has round tips, with its shape represented by the radius function $`h(z)`$ in cylindrical coordinates $`(r,z)`$. $`r^{}`$ denotes the radius of curvature at the tip. An electric field with strength $`E`$ is applied in the z direction (Fig. 1). Suppose the initial radius of the drop is $`a`$. Respectively we use $`a`$, $`\gamma /\eta `$, $`\gamma /a`$, $`(\gamma /aϵ_oϵ_2)^{\frac{1}{2}}`$ and $`(\gamma ϵ_o/aϵ_2)^{\frac{1}{2}}`$ to scale length, velocity, stress, electric field and surface charge density sherwood1988 .

Following Sherwood sherwood1988 , we study a situation in which Reynolds number is small so that the fluid flows via Stokes equation and the charge distributions are determined by electrostatics. The surface charge density $`\rho `$ can be expressed in the form of a boundary integral equation

$`{\displaystyle \frac{(\beta +1)}{2(\beta 1)}}\rho (x)`$ $`=`$ $`{\displaystyle _{L_y}}g(x,y)\rho (y)h(y)(1+h^{^{}}(y)^2)^{\frac{1}{2}}𝑑y`$ (1)
$`{\displaystyle \frac{Eh^{^{}}(x)}{(1+h^{^{}}(x)^2)^{\frac{1}{2}}}},`$
where $`\rho (x)`$ is the surface charge density at $`(h(x),x)`$, $`g`$ is a Green function; $`\beta `$ denotes the permittivity ratio $`ϵ_1/ϵ_2`$; $`L_y`$ is the range of z axis occupied by the drop sherwood1988 . From the surface charge density $`\rho `$, the normal and tangential component of the electric field can be calculated to obtain the jump in the electric stress across the interface. Sherwood also uses a boundary integral to determine the interface velocity

$$u_i(x)=\frac{1}{8\pi }_{L_y}G_{ij}(x,y)f_j(y)h(y)(1+h^{^{}}(y)^2)^{\frac{1}{2}}𝑑y,$$
(2)
where $`i`$ and $`j`$ refer to the $`z`$ or $`r`$ component, $`f_j(y)`$ is the $`j`$ component of the total surface stress, and $`G`$ denotes a Green function rallison1978 . The velocity, u, is then used to update the interface position. In our simulation, we apply a boundary element method with many details similar to that described by Sherwood. We distribute mesh points in proportion to the local curvature, and use a cubic spline to interpolate the interface between mesh points, a quartic polynomial to interpolate the surface charge density. The derived linear algebraic equations are solved by LU decomposition. A fourth-order Runge-Kutta scheme is applied to update the interface position.

There exists a critical electric field $`E_c(\beta )`$ for $`\beta >\beta _c`$. When $`E<E_c`$, the drop can reach equilibrium with round tips. We start our simulation from a sphere and apply a sufficiently small electric field. If the maximum velocity on the interface decreases to a value below $`10^4`$ following an exponential decay, we consider that the drop will reach equilibrium. After equilibrium is reached by a numerical extrapolation, the field is increased by a small amount. Through increasing the electric field step by step we find the critical electric field.

When $`E>E_c`$, a finite time singularity develops. From now on, we use $`\beta =18.5`$ as our example, which has $`E_c^2=0.4085\pm 0.0003`$. For instance, when we choose the initial shape to be the equilibrium shape at $`E^2=0.400<E_c^2`$ and suddenly apply $`E^2=0.410>E_c^2`$, the drop forms conical-like ends. The velocity at the tip dramatically increases as a critical time is approached (Fig. 2). Here we can at most obtain about twelve decades of data in one calculation, due to the increasing number of mesh points required and roundoff error.

Figure 3 shows how the shape and stresses evolve as the finite time singularity develops. The slope plot suggests we can partition the interface into three regions, the tip region, the conical-like region and the macroscopic region. The conical-like region is the intermediate region with a small variation in the slope. Figure 3(A) shows as $`r^{}`$ decreases, the macroscopic region and the established part of the cone region almost remain intact, while part of the interface which used to be in the tip region now grows conical-like. This shows that in the course of $`r^{}0`$, only the tip region changes in time, while the established part of the cone region nearly remains independent of time. A careful examination of Fig. 3(A) shows the conical-like region is not exactly a cone, because the slope of the newly-grown cone changes as $`r^{}`$ decreases. As we shall see in more detail later, this slope approaches the value set by $`\theta _<`$. Figure 3(B) shows the ratio of the electric stress to the surface tension is larger than one in the tip region and the conical-like region. The stress ratio also does not change in the part of interface whose shape remains as the conical-like region expands. Hence the numerical evidence says that the shape of each part of the almost conical region and the stress ratio within that part remain frozen as the tip gets smaller. However, as $`r^{}`$ changes the slope and the stress ratio of the newly-grown part change too. So Figure 2 and 3 may show an approach to a fixed point, but they do not show a fixed behavior themselves. Because there is a slow and not-quite uniform convergence to a fixed point, it is hard to estimate the critical time, $`t_c`$, from our raw data. For this reason, we shall henceforth plot our results against tip radius, $`r^{}`$ instead of trying to use $`t_ct`$.

A scaling study, sometimes called an order of magnitude analysis, enables us to estimate the sizes of different contributions, it indicates how the different regions affect one another.

These estimates show that except for a uniform advection, the stresses in the tip region determines the flow within that region and the subsequent shape of the tip. Specifically, the deformation of the tip region is dominantly caused by the local stress jump. The axial strain rate defined as $`u_z(x)/x`$ measures how fast the interface deforms due to the axial velocity. Using (2), we can express the contribution to $`u_z(x)/x`$ from the three regions. Respectively the tip region, the conical region and the macroscopic region have length scales $`r^{}`$, $`h(z)`$ and $`1`$ ($`r^{}h(z)1`$). The electric stresses in the three regions have an order of magnitude

$$\sigma _e\rho ^2\{\begin{array}{cc}E^2r^{2(\nu 1)}\hfill & \text{tip region}\hfill \\ E^2h(z)^{2(\nu 1)}\hfill & \text{conical-like region,}\hfill \\ E^2\hfill & \text{macroscopic region}\hfill \end{array}$$
(3)
with $`0<\nu <1`$. A similar result applies to the surface tension $`\sigma _s`$, but with $`\nu =1/2`$. An argument like that of Lister and Stone shows that the forces in the intermediate and macroscopic region simply advect the tip region without significant contribution to the strain rate lister1998 .

A followup study yang2002  will describe in more detail how the scaling analysis of the electric stress works. For the present purposes, it suffices to say that the shape in the tip region mostly determines the electric stress within that region, except for a coefficient which only depends on the shape in the other regions and the applied electric field. If we change the shape in the other regions, the electric stress in the entire tip region will be changed by a factor which is independent of the shape in the tip region. Changing the applied electric field will have the same effect. So after we reshape the rest part of the interface, we can restore the electric stress in the tip region by applying a different electric field of certain strength.

The scaling study permits us to construct approximate boundary conditions which then permits the accurate determination of the subsequent behavior of the tip. Basically whenever we are about to run out of mesh points, we cut off the part of interface far away from the tip and replace it by a new shape profile which takes fewer mesh points. Then we restore the electric stress in the entire tip region, which we can accomplish by adjusting the applied electric field to restore the electric stress at the tip to its value prior to the truncation. The tip regions of the prescribed new drop and the original drop will subsequently evolve in the same way, because the deformation of the tip region is primarily driven by the local stress jump. We define the rescaled axial distance $`\xi `$ and radius function $`H(\xi )`$ as

$$\xi =(zz_{tip})/r^{},H(\xi )=h(z)/r^{}.$$
(4)
We at least keep the part of interface with $`\xi 10^4`$ and typically match a spherical band to the center region, requiring the slope to be continuous at the truncation points. The center of the spherical band, which locates on the z axis, coincides with the center of the prescribed new drop. The error will be smaller if the truncation point is farther away from the tip. The “truncate and prescribe” idea was invented by Zhang and Lister zhang1999 .

Using the same initial condition as in Fig. 2, we calculate the evolution of the tip region for eighty decades of $`r^{}`$ with approximate boundary conditions. We truncate the drop for twenty-four times, starting at $`r^{}=10^8`$. Figure 4 shows that the approximate boundary condition produces the same result as the exact boundary condition without truncation at $`10^{12}r^{}10^8`$. The point with $`\xi =100`$ is pretty close to the cone region, so the slope and the stress ratio there can respectively reflect the angle and the stress ratio of the newly-grown cone. Out of the eighty decades of data, the curves at small $`r^{}`$ can be adequately fitted as $`c+b(r^{})^p`$ with $`p>0`$. Thus $`H^{^{}}(\xi =100)`$, $`\alpha (\xi =100)`$ and $`dr^{}/dt`$ each approach limiting values as $`r^{}0`$. Just below we shall show those limits are the same for different initial conditions.

Further simulations show there exists a fixed point behavior: the stress ratio of the newly-grown cone converges to a fixed value larger than unity. At $`\beta =18.5`$, $`\alpha (\xi =100)`$ will converge to $`1.130`$ as $`r^{}0`$, if it is close to $`1`$ when the drop starts to develop conical ends, regardless of the initial shapes. For example in Fig. 5, simulation (ii) in solid lines and simulation (iii) in dashed lines have different initial conditions. The dotted lines and dashed lines can all be excellently fitted by $`c+b(r^{})^p`$. As $`r^{}0`$, the two simulations give the same limits $`H^{^{}}(\xi =100)_c=0.4726`$, $`\alpha (\xi =100)_c=1.130`$ and $`(dr^{}/dt)_c=0.0198`$. As you may notice, we have obtained the same limits in the simulation in Fig. 4. The fitted values of the power are very close to each other in all the simulations, and we get $`p=0.013\pm 0.001`$. In simulation (i), we purposely choose a particular initial condition to let $`\alpha (\xi =100)`$ equal to $`\alpha (\xi =100)_c`$ at an early stage in the cone formation. We see that $`\alpha (\xi =100)`$ stably stays at $`\alpha (\xi =100)_c`$ for many decades of $`r^{}`$ with $`H^{^{}}(\xi =100)`$ and $`dr^{}/dt`$ equal to the limits obtained from the curve fitting. We have concrete numerical evidence that $`\alpha (\xi =100)_c=1.130`$ is a stable fixed point at $`\beta =18.5`$. This fixed point is finally approached in a power law of $`r^{}`$. And simulation (i) gives the solution at this fixed point.

At this fixed point, the tip region is self-similar and the intermediate region is a cone with the cone angle $`\theta _<(\beta )`$. The first twelve decades of data calculated with the exact boundary condition in simulation (i) reveal the solution at this fixed point. We find the following properties: (a) The shape profiles of the tip region are self-similar after we rescale them by $`r^{}`$, for example $`H(\xi )`$ at $`\xi =10`$ is a constant for small $`r^{}`$ \[Fig. 6(A)\]. (b) The ever expanding intermediate region has a constant slope equal to $`Tan\theta _<(\beta )`$ \[Fig. 6(B)\]. (c) The velocity at the tip increases logarithmically \[Fig. 6(C)\]. (d) The stress ratio in the intermediate region is a constant substantially larger than one, which we call it $`\alpha _c`$ \[Fig. 6(D)\]. The values of $`\alpha _c`$ and $`\alpha (\xi =100)_c`$ are very close to each other. (e) $`dr^{}/dt`$ is a constant \[Fig. 5(C)\], which indicates that $`r^{}`$ scales like $`t_ct`$. This self-similar solution has some qualitative similarity with a scaling solution found by Lister and Stone lister1998 .

We find similar fixed points at other values of $`\beta >\beta _c`$ such as $`19.0`$. To summarize, our numerical evidence shows a cone with the smaller cone angle $`\theta _<(\beta )`$ can be formed transiently in a non-equilibrium case where the electric stress is not balanced by the surface tension. The angle is approached as the stress ratio in the newly-grown cone region converges to a fixed value $`\alpha _c(\beta )>1`$. The dynamical solution at this fixed point is self-similar.

This project would be simply impossible without Leo P. Kadanoff’s guidance and support. I am very grateful to H. A. Stone and W. W. Zhang for their codes on viscous pinchoff. I also want to thank W. W. Zhang for helpful discussions. This research was supported in part by the DOE ASCI-FLASH program, the NSF grant DMR-0094569 to L. P. Kadanoff, and the MRSEC program of the National Science Foundation under Award No. DMR-9808595.


format: text

text:
# Isospin relations for the tau decay modes

## 1 Introduction

In 1953, at the Bagnères de Bigorre conference , which Roberto Salmeron attended as member of the Manchester group, Dalitz showed the following inequality :

$$1/4\frac{\tau ^+\pi ^+\pi ^0\pi ^0}{\tau ^+\pi ^+\pi ^+\pi ^{}}1.$$
(1)
The $`\tau ^+`$ in Eq. 1, partner of the $`\theta ^+`$ in the celebrated puzzle, is nowadays known as $`K^+`$. The hypotheses leading to Eq. 1 were the existence of an isospin triplet ($`\tau ^+,\tau ^0,\tau ^{}`$) and the conservation of isospin in the $`\tau `$ decay. None of them was founded. Nevertheless, owing to the $`|\mathrm{\Delta }\stackrel{}{I}|=1/2`$ rule, Eq. 1 survived the introduction of the Gell-Mann Nishijima scheme  and stimulated the discovery of the $`K^+\pi ^+\pi ^0\pi ^0`$ decay mode .

More than thirty years later, the same inequality was written in a paper  devoted to the “calculation of exclusive decay modes of the tau”, but the $`\tau `$ studied in the paper was the $`\tau `$ lepton. This coincidence can serve to illustrate the longevity and the generality of the isospin relations. However, in the following pages, we will discuss only their applications  to the decay modes of the $`\tau `$ lepton.

## 2 Isospin relations in hadronic final states

The proof of an equation like Eq. 1 starts from the identification of the different possible isospin states for the hadronic system. The amplitudes for a given charge-configuration ($`\pi ^+\pi ^+\pi ^{}`$ and $`\pi ^+\pi ^0\pi ^0`$ for three $`\pi `$’s with $`Q=1`$) are linear combinations of the isospin amplitudes; thus the partial widths are linear combinations of the squared amplitudes and their interference terms.

Generally, a large fraction of the interference terms are killed by the integration over the phase-space. The remaining terms are bounded by the Schwartz inequality . The resulting constraints can be geometrically represented by an allowed convex domain in the space of the charge-configuration fractions, $`f_{cc}=\mathrm{\Gamma }_{cc}/\mathrm{\Gamma }`$. If all the interferences vanish, the domain is a polyhedron, convex hull of the points that describe each isospin state.

Hence, the first step in the construction of the allowed domain is the setting up of a basis for the isospin states, adapted to the implementation of the Pauli principle.

## 3 Isospin states of $`n\pi `$ systems

Such a basis was constructed by Pais  with the object of studying the many pion systems produced in $`\overline{p}p`$ and $`\overline{p}n`$ annihilations. The construction is based on two simple remarks : i) the representations of the isospin group SU(2) relevant for $`n\pi `$ systems are also representations of SO(3), ii) the group of $`3\times 3`$ orthogonal unimodular matrices, SO(3), is a subgroup of the group of $`3\times 3`$ unimodular matrices, SL(3). Thus, if $`V`$ is the three-dimensional space of the isospin states of one pion, the space of the states of $`n`$ pions, $`V^n`$, supports representations of both SL(3) and the symmetric group $`𝒮_n`$, which acts on $`V^n`$ by permuting its factors. Standard properties of the representations of linear groups  imply the decomposition

$$V^n=\stackrel{\lambda }{}E_\lambda F_\lambda ,$$
(2)
where the symbol111Also denoted  “symmetry class”. $`\lambda =(\lambda _1,\lambda _2,\lambda _3)`$, with $`\lambda _1\lambda _2\lambda _3`$ and $`\lambda _1+\lambda _2+\lambda _3=n`$, is associated to a three-row Young diagram; $`F_\lambda `$ is the irreducible representation of $`𝒮_n`$ determined by the diagram, and $`E_\lambda `$ an irreducible representation of SL(3). For a given $`n`$, the correspondence between the representation of SL(3) and $`\lambda `$ is one–one. The representation $`E_\lambda `$ is characterized by the two numbers: $`\lambda _1\lambda _2`$ and $`\lambda _2\lambda _3`$.

As a representation of SO(3), $`E_\lambda `$ is no longer irreducible because of the invariance under SO(3) of the contraction operation . Its decomposition into irreducible representations of SO(3) reads

$$E_\lambda =\stackrel{I}{}N_I(\lambda )D_I,$$
(3)
where $`D_I`$ is the $`(2I+1)`$-dimensional irreducible representation (integer isospin I). The multiplicity $`N_I(\lambda )`$ was computed by Racah . It can be written:

$$N_I(\lambda )=\varphi (\lambda _1\lambda _3I+2)\varphi (\lambda _2\lambda _3I+1)\varphi (\lambda _1\lambda _2I+1),$$
(4)
where $`\varphi (x)`$ is the greatest integer contained in $`x/2`$ for $`x>0`$, and $`0`$ for $`x0`$. For systems of two pions, we get:

|  | $`\lambda `$ | $`N_0`$ | $`N_1`$ | $`N_2`$ |
| --- | --- | --- | --- | --- |
|  | (1,1,0) | 0 | 1 | 0 |
|  | (2,0,0) | 1 | 0 | 1 |

and, for three pions,

|  | $`\lambda `$ | $`N_0`$ | $`N_1`$ | $`N_2`$ | $`N_3`$ |
| --- | --- | --- | --- | --- | --- |
|  | (1,1,1) | 1 | 0 | 0 | 0 |
|  | (2,1,0) | 0 | 1 | 1 | 0 |
|  | (3,0,0) | 0 | 1 | 0 | 1 |

For any $`\lambda `$, Eq. 4 implies the relation

$$N_0(\lambda )+N_1(\lambda )=1,$$
(5)
which can also be obtained directly from the symmetry ($`\mathrm{}\mathrm{}`$) of a contraction. If both $`\lambda _1\lambda _2`$ and $`\lambda _2\lambda _3`$ are even, $`N_0=1`$, otherwise $`N_1=1`$.

From the equations (2), (3) and (5), we see that the sole degeneracy of the $`n`$-$`\pi `$ states $`|I_3,\lambda `$ with $`I1`$, ($`I_3`$ and $`\lambda `$ fixed) is due to the permutation symmetry. These states form an irreducible representation of $`𝒮_n`$. Since a permutation preserves the charge-configuration ($`n^+,n^0,n^{}`$) of a $`n`$-$`\pi `$ system, the irreducibility implies that the Clebsch-like coefficients used to write the isospin states as combinations of charge-configuration states are determined by $`I_3`$ and $`\lambda `$ only222The elements of the group algebra used to build an orthogonal basis of the representation give also orthogonal charge configuration states .. In other words , the coefficients of the charge-configurations are a class property.

Furthermore, since the permutation symmetry properties of the momentum and isospin amplitudes are the same, because of the Pauli principle, integrating over the phase-space kills all the interference terms for a $`n`$-$`\pi `$ system with $`I1`$; the allowed domain in the space of the charge-configuration fractions is a polyhedron.

Let’s take the simple example of the $`Q=1`$, $`I=1`$ three-$`\pi `$ system alluded to in the introduction. For SO(3), an isospin one is a vector, and any vector made of three vectors can be written

$$\alpha (\stackrel{}{b}\stackrel{}{c})\stackrel{}{a}+\beta (\stackrel{}{c}\stackrel{}{a})\stackrel{}{b}+\gamma (\stackrel{}{a}\stackrel{}{b})\stackrel{}{c}.$$

As an isospin function, $`(\stackrel{}{a}\stackrel{}{b})\stackrel{}{c}`$ describes an $`I=1`$ state with the two first $`\pi `$ in an $`I=0`$ state: $`(\pi ^+\pi ^{}\pi ^0\pi ^0+\pi ^{}\pi ^+)\pi ^+`$. It is then straightforward to write the other terms by cyclic permutations and get the ratio

$$R=\frac{\pi ^+\pi ^0\pi ^0}{\pi ^+\pi ^+\pi ^{}}=\frac{|\alpha |^2+|\beta |^2+|\gamma |^2}{|\beta +\gamma |^2+|\alpha +\gamma |^2+|\alpha +\beta |^2}.$$
(6)
The symmetry class $`\lambda =(3,0,0)`$, is associated to the one-dimensional space of completely symmetric states ($`\alpha =\beta =\gamma `$) for which the ratio is $`R=1/4`$; the class $`\lambda =(2,1,0)`$ to the two-dimensional orthogonal space ($`\alpha +\beta +\gamma =0`$), with $`R=1`$. An example of state in the $`(2,1,0)`$ class is given by a $`\rho \pi `$ system, which is represented by $`(\stackrel{}{a}\stackrel{}{b})\stackrel{}{c}`$, where the vector product ($``$) is interpreted as the combination of two isospins one into an isospin one.

Thus the charge-configuration fractions are:

$$f_{\pi ^+\pi ^0\pi ^0}=\frac{1}{2}W_{(210)}+\frac{1}{5}W_{(300)},f_{\pi ^+\pi ^+\pi ^{}}=\frac{1}{2}W_{(210)}+\frac{4}{5}W_{(300)},$$
(7)
with $`W_{(210)}+W_{(300)}=1`$.

The weights $`W_\lambda `$ depend on the dynamics. In the Fermi statistical model , they are proportional to the dimension of the representation $`F_\lambda `$ of the permutation group $`𝒮_n`$; here $`W_{(210)}^{stat}=2/3`$ and $`W_{(300)}^{stat}=1/3`$.

If $`I=2`$ states are allowed, they have to share the symmetry properties of $`I=0`$ or $`I=1`$ states, so that interference terms must be taken into account. It can be checked on Eq. 4 that, for $`n<6`$, $`N_2(\lambda )1`$ for all $`\lambda `$. Thus, for $`n<6`$ and $`I2`$, the charge-configuration coefficients are determined by $`I`$, $`I_3`$ and $`\lambda `$ only.

Tables of the coefficients can be found in the literature . They are computed by explicitly constructing tensors with the required symmetry, as in the previous example, or by more sophisticated methods .

## 4 Semileptonic decays of the $`\tau `$ lepton

The possible hadronic systems333The $`\eta `$ channels have to be treated separately because $`\eta `$ decay violates isospin. ($`h`$) in the semileptonic decay $`\tau \nu h`$ are: $`n\pi `$, $`\eta n\pi `$, $`Kn\pi `$, $`K\eta n\pi `$, and $`K\overline{K}n\pi `$.

The properties of the charged weak current imply that the total isospin is 1 for the strangeness zero final states, and 1/2 for the others. Thus the isospin of the $`n\pi `$ system is 1 for $`h=n\pi `$ and $`h=\eta n\pi `$; 0 or 1 for $`h=Kn\pi `$; 0, 1 or 2 for $`h=K\overline{K}n\pi `$.

For all the cases but $`h=K\overline{K}n\pi `$ the isospin amplitudes can be labelled by the symmetry class $`\lambda `$ only and the interference terms are killed by the integration over phase-space. Thus the partial width for a given charge-configuration ($`cc`$) can be written:

$$\mathrm{\Gamma }^{cc}=_\lambda C_\lambda ^{cc}\mathrm{\Gamma }^\lambda ,$$
(8)
where the coefficients $`C_\lambda ^{cc}`$ can be found in tables .

For positive G-parity systems ($`h=2n\pi `$, $`h=\eta (2n+1)\pi `$), the weak current is related by an isospin rotation to the electromagnetic current, hence relations can be established between the $`\tau `$ partial widths and the cross sections of $`e^+e^{}`$ annihilations.

For $`h=K\overline{K}n\pi `$, a more detailed analysis taking into account the interferences is needed.

### 4.1 $`\tau \nu (2n+1)\pi `$

For a $`3\pi `$ system, there are two possible symmetry classes and the isospin constraint is Eq. 1.

Experimentally, the ratio $`\pi ^{}\pi ^{}\pi ^+/\pi ^{}\pi ^0\pi ^0`$ is nearly 1 because of the dominance of the $`\rho \pi `$ intermediate state. A detailed analysis of the final state , taking into account the isospin symmetry breaking caused by the difference of the $`\pi ^0`$ and $`\pi ^\pm `$ masses, predicts a ratio 0.985, in good agreement with the measurements .

For a $`5\pi `$ final state, three charge-configurations and four symmetry classes are present. The production of $`\eta `$ contributes to the $`2\pi ^{}\pi ^+2\pi ^0`$ and $`\pi ^{}4\pi ^0`$ final states.

The comparison of the measured branching ratios  (after $`\eta `$ subtraction) with the allowed domain is made in Fig. 1. It shows the dominance of $`\lambda =(2,2,1)`$, which is due to the $`\omega \pi ^{}\pi ^0`$ intermediate state.

### 4.2 $`\tau \nu 2n\pi `$

The final states of an even number of pions are produced by the vector current, which is related by an isospin rotation to the electromagnetic current. Since the isospin rotation commutes with the permutations, the relation between $`\tau `$ partial widths and $`e^+e^{}`$ cross sections  can be written for each symmetry class $`\lambda `$:

$$\frac{1}{\mathrm{\Gamma }_{\nu \overline{\nu }e}}\frac{d\mathrm{\Gamma }_{\nu 2n\pi }^\lambda }{dm^2}=\frac{3\mathrm{cos}\theta _c^2}{2\pi \alpha ^2m_\tau ^8}m^2(m_\tau ^2m^2)^2(m_\tau ^2+2m^2)\sigma _{e^+e^{}2n\pi }^\lambda (m^2).$$
(9)
There is only one class for $`2\pi `$ systems. For $`4\pi `$ final states the possible classes are $`\lambda =(3,1,0)`$ and $`\lambda =(2,1,1)`$. Since there are two charge-configurations for both $`\tau `$ decays and $`e^+e^{}`$ annihilations, the correspondence between cross sections and partial widths is very simple:

$`\sigma _{2\pi ^+2\pi ^{}}`$ $``$ $`2\mathrm{\Gamma }_{\pi ^{}3\pi ^0}`$
$`\sigma _{\pi ^+\pi ^{}2\pi ^0}`$ $``$ $`\mathrm{\Gamma }_{\pi ^{}\pi ^{}\pi ^+\pi ^0}\mathrm{\Gamma }_{\pi ^{}3\pi ^0}.`$
Thorough comparisons  of $`\tau `$-decay and $`e^+e^{}`$-annihilation data for the two- and four-pion channels, including the consideration of isospin symmetry breaking, have been made recently in order to improve the theoretical determination of the muon anomalous magnetic moment $`a_\mu `$. They show some discrepancies between the $`e^+e^{}`$ and $`\tau `$ data, as well as between different $`e^+e^{}`$ experiments.

Four classes can contribute to the $`6\pi `$ states production, and only three charge-configurations are possible in $`\tau `$ decays as well as in $`e^+e^{}`$ annihilations. Therefore, even with complete measurements it would not be possible to predict the partial widths from the cross sections or conversely. Nevertheless, with two measurements in $`e^+e^{}`$ annihilations ($`\sigma _{3\pi ^+3\pi ^{}}`$, $`\sigma _{2\pi ^+2\pi ^{}2\pi ^0}`$) and in $`\tau `$ decays ($`B_{3\pi ^{}2\pi ^+\pi ^0}`$, $`B_{2\pi ^{}\pi ^+3\pi ^0}`$), it is possible to determine the contributions of the four classes, or, at least, to check the consistency of the different measurements.

Figure 2 displays the allowed region in the plane of the charge-configuration fractions and the ratio $`B_{2\pi ^{}\pi ^+3\pi ^0}/B_{3\pi ^{}2\pi ^+\pi ^0}`$ of the measurements , after $`\eta `$ subtraction (grey area).

The large $`\sigma _{2\pi ^{}2\pi ^+2\pi ^0}/\sigma _{3\pi ^{}3\pi ^+}`$ ratio observed in $`e^+e^{}`$ annihilations implies the dominance of $`\lambda =(3,2,1)`$ and/or $`\lambda =(4,1,1)`$, quantitatively shown  by the diagonal hatched region in Fig. 2. The location of the intersection of the two regions corresponds to $`\lambda =(3,2,1)`$, in agreement with the observed  importance of the

$`\omega `$ production.

Discrepancies appear when the total cross section for $`e^+e^{}`$ annihilations into $`6\pi `$ is taken into account. It allows an estimation  of the total $`B_{6\pi }`$ branching ratio, which together with the measured $`B_{3\pi ^{}2\pi ^+\pi ^0}`$ gives the vertical hatched band in Fig. 2, clearly incompatible with the other estimations.

Two hypotheses can be contemplated: either the $`e^+e^{}`$ cross sections are overestimated by a factor of roughly four, or the $`e^+e^{}`$ annihilations into $`6\pi `$ receive a large contribution from $`I=0`$, $`\eta 3\pi `$ final states. In the second hypothesis, the cross section for $`e^+e^{}\pi ^+\pi ^{}+neutrals`$ would be three times larger than $`\sigma _{2\pi ^{}2\pi ^+2\pi ^0}`$.

The contribution of the $`6\pi `$ channel to the estimation of $`a_\mu `$ is small , however the second hypothesis, if true, could have a not completely negligible impact on the estimation.

### 4.3 $`\tau \nu Kn\pi `$

For $`\tau \nu K\pi `$, the isospin 1/2 implies the ratio $`K^{}\pi ^0/\overline{K}^0\pi ^{}=1/2`$, to be compared with the experimental value

$$B_{K^{}\pi ^0}/B_{\overline{K}^0\pi ^{}}=0.51\pm 0.04.$$

For $`Kn\pi `$ systems, the number of charge-configurations is greater than the number of symmetry classes. The resulting relations between branching ratios are

$`B_{K^{}\pi ^+\pi ^{}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}B_{\overline{K}^0\pi ^{}\pi ^0}+2B_{K^{}\pi ^0\pi ^0}`$ (10)
$`B_{\overline{K}^0\pi ^{}\pi ^{}\pi ^+}`$ $`=`$ $`B_{\overline{K}^0\pi ^{}\pi ^0\pi ^0}+2B_{K^{}\pi ^0\pi ^0\pi ^0}`$ (11)
for three and four hadron final states. The data for the $`K2\pi `$ final states are shown in Fig. 3. The location of the experimental point is consistent with the observed dominance of the intermediate states $`K^{}\pi `$ and $`K\rho `$.

### 4.4 $`\tau \nu K\overline{K}n\pi `$

In a $`K\overline{K}`$ system, the possible values of the $`K\overline{K}`$ isospin are 0 and 1. For $`I_{K\overline{K}}=0`$ the isospin of the $`n\pi `$ system is $`I_{n\pi }=1`$; it can be 0, 1 or 2 for $`I_{K\overline{K}}=1`$. Both axial and vector currents contribute to the decay, leading to $`G=+1`$ (V) and $`G=1`$ (A) for the G-parity of the hadronic system, and $`G_{K\overline{K}}=(1)^nG`$.

Since the $`J^P`$ quantum numbers are different for the axial and vector currents, there is no V-A interference term in the partial widths.

The G-parity of a $`K\overline{K}`$ system is related to its isospin and orbital momentum by $`G_{K\overline{K}}=(1)^{I_{K\overline{K}}+l_{K\overline{K}}}`$, thus for a given (V or A) current the values of $`l_{K\overline{K}}`$ associated to $`I_{K\overline{K}}=0`$ and $`I_{K\overline{K}}=1`$ are different and there is no interference term between $`I_{K\overline{K}}=0`$ and $`I_{K\overline{K}}=1`$ amplitudes.

This imply , for any charge-configuration ($`cc`$) of the $`n\pi `$ system, the equality.

$$\mathrm{\Gamma }_{K^0\overline{K}^0(n\pi ,cc)}=\mathrm{\Gamma }_{K^+K^{}(n\pi ,cc)}.$$
(12)
From CPT invariance, $`\mathrm{\Gamma }_{K_SK_S(n\pi ,cc)}=\mathrm{\Gamma }_{K_LK_L(n\pi ,cc)}`$, but the ratio $`\mathrm{\Gamma }_{K_SK_L(n\pi ,cc)}/\mathrm{\Gamma }_{K_SK_S(n\pi ,cc)}`$ is a free parameter depending on dynamics and the respective contributions of V and A currents .

Figure 4 shows the agreement of the data  with Eq. 12 in the case of the $`K\overline{K}\pi `$ final state. The ratio $`K^+K^{}\pi ^+/K^0K^{}\pi ^0`$ is found equal to 1 in agreement with the observed dominance of the intermediate state $`K^{}\pi `$.

The isospin amplitudes for a $`K\overline{K}n\pi `$ final state can be labelled by $`I_{K\overline{K}}`$, $`I_{n\pi }`$ and the symmetry class $`\lambda `$ of the $`n\pi `$ system.

If no $`I_{n\pi }=2`$ amplitude is associated with $`\lambda `$, there is no possible interference term and the class is described by a point in the space of the charge-configuration fractions.

If $`I_{n\pi }=2`$ is possible, there is one interference term, but only one since $`m_\tau <2m_K+6m_\pi `$ (section 3), thus the class is described by a two-dimensional elliptic domain .

The allowed domain in the space of the charge-configuration fractions is the convex hull of the points and ellipses associated with the symmetry classes of the $`n\pi `$ system. Fig. 5 shows a projection of this multidimensional domain in the case of $`n=2`$. It implies the following inequality :

$$B_{K^0K^{}\pi ^0\pi ^0}\frac{3}{4}(B_{\overline{K}^0K^+\pi ^{}\pi ^{}}+B_{K^0K^{}\pi ^+\pi ^{}}),$$
(13)
which, together with the relation $`B_{\overline{K}^0\pi ^{}\pi ^0\pi ^0}B_{\overline{K}^0\pi ^{}\pi ^{}\pi ^+}`$ (Eq. 11), gives the constraint

$$B_{K_Sh^{}\pi ^0\pi ^0}<B_{K_Sh^{}h^{}h^+}$$
(14)
on topological branching ratios.

## 5 Conclusion

Isospin relations have been used in $`\tau `$ lepton physics for nearly twenty years. Their first applications  were the estimation of $`\tau `$ branching ratios from $`e^+e^{}`$ annihilation data, the bounding of the contributions of unobserved channels, and the elucidation of the “one prong problem” .

Today, a large number of decay modes are tabulated , the order of magnitude of the smallest measured branching ratios is $`10^4`$, and the data from $`\tau `$ decays are used  to complement and correct the information given by the $`e^+e^{}`$ annihilations.

The observed discrepancies can only be solved by experiment. However, quoting from Blackett’s closing remark at the Bagnères conference : “if the history of scientific discovery is any guide, the same increase of accuracy which will serve to settle our present controversies will equally, surely, itself bring to birth new controversies by leading of some discoveries.”


text:
# Astrophysical components separation of 𝐶𝑂𝐵𝐸-DMR 4yr data with FastICA

## 1 Introduction

Current and future CMB space missions such as $`WMAP`$ (launched in June 2001, Bennett et al. 1996) and Planck (scheduled for launch in February 2007, Tauber 2000) will map the microwave sky emission over the entire sky with an unprecedented combination of sensitivity and angular resolution. In the meantime a large variety of ground-based and balloon-borne experiments will provide accurate data over various regions of the sky.

The main goal of such experiments is to produce a map cleaned of any contributions from other sources between us and the Last Scattering Surface and thus containing only genuine CMB anisotropy. From such maps the angular power spectrum is evaluated and used to determine cosmological parameters with high accuracy and thereby probe various scenarios for structure formation. Potential foreground contaminants include emission from our Galaxy (mainly synchrotron, free-free and dust emission), compact galactic and extra-galactic sources, and the thermal and kinematic Sunyaev-Zel’dovich effect from clusters of galaxies. In order to exploit the cosmological information encoded into the CMB angular power spectrum, it is crucial to be able to identify and remove such signal components with high accuracy and reliability.

Much recent work has been performed in this area by several authors (de Oliveira-Costa & Tegmark 1999, Tenorio et al.1999, Hobson et al. 1998, 1999, Stolyarov et al. 2001, Prunet et al. 2001) employing a range of techniques from the classical Wiener Filter to maximum entropy methods. These algorithms are referred to as “non-blind” in the sense that they require some a-priori information on the signal to be separated (e.g. spatial templates and frequency dependences of the underlying components, although see Barreiro et al. 2003 for a variation on this approach using a maximum entropy based technique). Recently a blind approach has been proposed (Baccigalupi et al. 2001, Maino et al. 2002) and applied to simulated sky maps similar to those that Planck will produce. The approach appears to be very promising in that it is fast and does not need priors about the underlying signals (at least for high sensitivity missions such as $`WMAP`$ or Planck). It is therefore possible to use such blind algorithms to obtain priors that can be fed into classical Bayesian separation algorithms.

However, to-date these blind algorithms have not been applied to real CMB data to assess and validate their reliability in a real world situation. We present here the first application of FastICA to real CMB measurements from the $`COBE`$-DMR 4yr data.

This paper is organised as follows: in Section 2 we briefly present the component separation problem and the main simplifications required for the application of FastICA. In Section 2.1 and 2.2 we summarise the basics of the FastICA algorithm and its application to $`COBE`$-DMR data. Section 3 shows the impact of the noise amplitude and distribution on FastICA outputs which drives the choice of optimal separation technique. The main CMB-related results are presented in Section 4 while Section 6 considers the implications for foreground emission. A critical discussion follows in Section 7.

## 2 Component Separation Problem

Let us suppose that the observed sky radiation is the superposition of $`N`$ different physical processes and its frequency and spatial dependences can be factorized into two separated terms:

$$\stackrel{~}{x}(\text{r},\nu )=\underset{j=1}{\overset{N}{}}\overline{s}_j(\text{r})f_j(\nu ).$$
(1)
This signal is, in general, observed by an experiment with $`M`$-frequencies through an optical system, whose beam pattern is in general modelled at each frequency as a shift-invariant point spread function $`B(\text{r},\nu )`$. Let us further suppose that $`B(\text{r},\nu )`$ is frequency-independent at least within each frequency bandwidth $`t(\nu )`$. In addition any real experiment adds some instrumental noise to the output $`ϵ_\nu (\text{r})`$. Following our assumptions the observed signal at a frequency $`\nu `$ is given by:

$`x_\nu (\text{r})`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}B_\nu (\text{r})\overline{s}_j(\text{r}){\displaystyle t_\nu (\nu ^{})f_j(\nu ^{})𝑑\nu ^{}}+ϵ_\nu (\text{r})`$ (2)
$`=`$ $`B_\nu (\text{r}){\displaystyle \underset{j=1}{\overset{N}{}}}a_{\nu j}\overline{s}_j(\text{r})+ϵ_\nu (\text{r})`$
where $``$ denotes convolution, and

$$a_{\nu j}=t_\nu (\nu ^{})f_j(\nu ^{})𝑑\nu ^{}.$$
(3)
Our data model can be further simplified assuming that the radiation pattern of the telescope is frequency independent i.e. $`B_\nu (\text{r})=B(\text{r})`$. In this case, Eq. (2) can be written in vector form as:

$$𝐱(\text{r})=𝐀\overline{𝐬}(\text{r})B(\text{r})+\mathit{ϵ}(\text{r})=\mathrm{𝐀𝐬}(\text{r})+\mathit{ϵ}(\text{r}),$$
(4)
where each component, $`s_j`$, of the vector $`𝐬`$ is the corresponding source function convolved with the $`B`$ beam pattern. The matrix $`𝐀`$ is the mixing matrix with elements given by the $`a_{\nu j}`$ coefficients.

We present this derivation in order to stress data model assumptions. It is worth noting that one or more of these might be not justified in real cases. For instance both the $`WMAP`$ and Planck experiments will observe the sky radiation through a telescope with multi-frequency receivers and the resulting beam pattern is strongly frequency dependent. This is up to now one of the main limitation of the FastICA approach to astrophysical components separation and forced Maino et al. (2002) to further convolve simulated Planck sky maps in order to obtain similar beam functions. Furthermore the noise term $`ϵ(\text{r})`$ is usually assumed to be additive, signal-independent, white, Gaussian, stationary and uniformly distributed on the sky. These are additional strong assumptions since the noise spectrum of a real experiment may contain a low-frequency tail due to the so-called $`1/f`$ noise and different scanning (observing) strategies which distribute the integration time on the sky in a non-uniform way.

### 2.1 FastICA algorithm

We briefly summarise here the FastICA algorithm. The problem of obtaining both the mixing matrix $`𝐀`$ and the signals $`𝐬`$ from observed data $`𝐱`$ is unsolvable if additional information is not provided. The ICA approach assumes that

* the signals $`𝐬`$ are independent random processes on the map domain;
* all the signal, but at most one, have non-Gaussian distribution.

The strategy we exploit here is described in detail in Hyvärinen & Oja (1997) and Hyvärinen (1999) and its application in an astrophysical context can be found in Maino et al. (2002): independent components are extracted maximising a suitable measure of non-Gaussianity. Indeed the central limit theorem states that a variable which is a mixture of independent variables is “more Gaussian” than the original ones. Therefore we have to find a transformation such that the Gaussianity of the variables is reduced: this is equivalent to finding a set of transformed variables that are “more independent” than the original ones. Furthermore since our data model is a noisy, we have to define a measure of non-Gaussianity that is robust against noise, this is the the so-called neg-entropy.

Approximations to neg-entropy have been given by Hyvärinen & Oja (2000) and Hyvärinen (1999) and, if the noise has the properties assumed previously and its covariance matrix is known, the Gaussian moments of the transformed variables $`𝐲`$ = $`\mathrm{𝐖𝐱}`$ are shown to be robust estimates of the desired functions. Here the matrix $`𝐖`$ is the separation matrix such that the $`𝐲`$ components are in fact independent.

The algorithm needs a preprocessing step (Hyvärinen 1999) in which the input maps are “quasi whitened”. This reduces the number of unknowns in the problem. Let us assume that we know the covariance matrix $`𝚺`$ of the instrumental noise; at each frequency, the mean value is removed from the data (the offsets of each independent component can be recovered at the end of the separation process) and their covariance matrix $`𝐂`$ is evaluated by computing the following expectation value:

$$𝐂=\mathrm{𝐱𝐱}^T.$$
(5)
A modified covariance matrix and quasi-whitened data sets are respectively given by:

$$\widehat{𝚺}=(𝐂𝚺)^{1/2}𝚺(𝐂𝚺)^{1/2},$$
(6)
$$\widehat{𝐱}=(𝐂𝚺)^{1/2}𝐱.$$
(7)
The separation matrix $`𝐖`$ is estimated row by row i.e. one component at a time. Let $`𝐰`$ be a $`M`$-vector such that $`𝐰^T\widehat{𝐱}`$ gives one component of the transformed vector $`𝐲`$ ($`𝐰^T`$ is a row of the separation matrix $`𝐖`$). In order to find an estimation of the transformed vector $`𝐲`$ that is robust against noise, the following iterative algorithm is applied, together with a convergence criterion:

1. choose an initial vector $`𝐰`$
2. update it by means of:

$$𝐰_{\mathrm{new}}=\widehat{𝐱}g(𝐰^T\widehat{𝐱})(I+\widehat{𝚺})g^{}(𝐰^T\widehat{𝐱})$$

where $`g`$ is a regular non-quadratic function i.e. $`g(u)=u^3`$, $`g(u)=\mathrm{tanh}(u)`$ and $`g(u)=u\mathrm{exp}(u^2)`$.
3. normalise $`𝐰_{\mathrm{new}}`$ to be a unit vector
4. compare $`𝐰_{\mathrm{new}}`$ with the old computed value; if not converged come back to (ii); if converged, begin another row.

This procedure maximises the non-Gaussianity of the component $`𝐰^T\widehat{𝐱}`$. If $`k`$ rows of the matrix $`𝐖`$ are found at a given time, the $`k+1`$ row is searched for in a sub-space orthogonal to the first $`k`$ rows. To this purpose a orthogonalisation procedure (e.g. by mean of Gram-Schmidt rule) is inserted between point (ii) and (iii). Once the separation matrix $`𝐖`$ is obtained the underlying components are derived by using:

$$𝐱=𝐖^\mathrm{𝟏}𝐲$$
(8)
This equation allows us to derive the frequency scalings for each independent component: the scaling between $`\nu `$ and $`\nu ^{}`$ of the $`j^{\mathrm{th}}`$ component is given by the ratio of $`W_{\nu j}^1/W_{\nu ^{}j}^1`$.

It is also possible to recover signal-to-noise ratio for the reconstructed components. Since the noise covariance $`𝚺`$ is supposed to be known, noise constrained realization $`𝐧_𝐱`$ for each frequency channels can be built. Once we have the separation matrix $`𝐖`$, the noise realization in FastICA, outputs are given by $`\mathrm{𝐖𝐧}_𝐱`$. Noise is transformed like signals and an estimation of the noise in the reconstructed component is achieved. This is quite useful: suppose that we have $`M`$ frequency channels but only $`N<M`$ components have to be separated. Using all $`M`$ channels we can reconstruct $`M`$ components and the “fictitious” ones will be characterised by a signal-to-noise ratio lower than unity.

### 2.2 Application to $`COBE`$-DMR 4yr data

The DMR experiment has observed full-sky microwave emission at three frequencies (31.5, 53 and 90 GHz). At each frequency, two separate channels (denoted A and B) measured the temperature difference between two horn antennas of angular resolution $`7^{}`$ pointing in directions separated by $`60^{}`$ on the sky. This differencing information was then used to reconstruct full-sky maps of the microwave sky.

The instrumental noise is almost white and Gaussian but, due to the scanning strategy, it is not uniformly distributed on the sky: variations by a factor of $`4`$ is typically observed at 90 and 53 GHz while factor of $`6`$ is found at 31.5 GHz. This larger factor reflects the fact that part of the original 31 GHz data have been discarded due to a well known systematic effect (Kogut et al. 1996). In this way the 31.5 GHz channels shows a noise distribution on the sky which is substantially different from that at 53 and 90 GHz. Furthermore the $`S/N`$ ratio is quite poor: for 10 effective resolution smoothed maps this is $`0.5,1.5`$ and 1 when combining A and B channels at 31.5, 53 and 90 GHz respectively, and becomes $`2`$ when combining together all the frequencies. However, no extra-smoothing is required before applying FastICA since the angular resolution is the same for each channel.

We have applied the FastICA algorithm to the $`COBE`$-DMR 4yr data in HEALPix format (Górski et al. 1999) with a resolution parameter $`N_{\mathrm{side}}`$ = 32 corresponding to 12288 pixels in the sky with size $`1.83^{}`$. The A and B radiometers at each frequency are combined to obtain a sum map (A+B)/2.

Maino et al.(2002) have shown that the optimal CMB reconstruction, both in terms of frequency scaling and offset normalisation, can be achieved by exploiting data at high galactic latitudes. We use the extended custom galactic cut (Banday et al. 1997) to select regions useful for this CMB reconstruction. Before applying the code the best fit monopole, dipole and quadrupole – computed using only those pixels surviving the Galactic cut – are subtracted from the maps. FastICA “naturally” removes a monopole from the input sky maps but does not do the same for dipole and quadrupole distributions. These modes are removed since we expect that structures on the largest angular scales are contaminated most significantly by a galactic emission which is characterised by a steeply falling power spectrum. Since one figure of merit for the separation performed by FastICA is the frequency scaling, this procedure minimises the contamination of the derived CMB scaling by foreground emission.

We work with all the non-quadratic functions described before: $`g(u)=u^3`$, $`g(u)=\mathrm{tanh}(u)`$ and $`g(u)=u\mathrm{exp}(u^2)`$ identified by p, t and g, respectively.

Finally, we make use of the $`S/N`$ ratio and frequency scaling of each reconstructed component to assess its physical nature.

## 3 Impact of noise distribution and amplitude

One of the assumptions of FastICA is that the noise is Gaussian, stationary and uniformly distributed on the sky. The latter is manifestly not true for DMR. In order to validate the impact of noise amplitude and distribution on the reconstruction procedure, we performed $`3000`$ simulations of fake CMB skies (characterised by a Harrison-Zeldovich initial power law spectrum with an rms quadrupole normalisation of $`18\mu `$K) and instrumental noise. Specifically we consider three cases:

1. DMR nominal noise amplitude and distribution
2. uniform noise distribution (i.e. each pixel is observed the same number of times) with the noise amplitude matched to the mean noise amplitude in the actual DMR data (0.322, 0.100 and 0.139 mK at 31.5, 53 and 90 GHz, respectively)
3. uniform noise distribution with amplitude reduced by a factor of 4 with respect the nominal values

For each realisation we performed a component separation on the cut sky with the p, g and t functions. We have considered the cases where all three frequencies are included, or just the 90 and 53 GHz data (since we know that the 31.5 GHz channel shows a different noise pattern on the sky). In each case we evaluate the frequency scaling (between 90 and 53 GHz) of the reconstructed CMB component as well as its angular power spectrum on the cut sky exploiting the technique described in Górski et al. (1994) and in Górski et al. (1996, hereafter G96).

Table 1 reports the results for the derived frequency scaling for the p function (similar results are obtained with g and t). The results are almost identical for noise cases ii) and iii) with either 2 and 3 DMR frequencies, though naturally, with higher sensitivity a higher accuracy is achieved for the derived CMB frequency scaling. However, the situation is quite different for the nominal noise distribution and amplitude (case i). Although the mean value is almost the same and consistent with the expected theoretical value, the rms is somewhat degraded being 0.477 and 0.205 for 3 and 2 frequency channels separation, respectively. Therefore a first indication of the relevance of the noise distribution is derived from an inspection of the frequency scalings.

A similar situation is found in the analysis of the reconstructed CMB component power spectrum as shown in Figure 1. The CMB component recovered when all three DMR sum maps are inputs to the FastICA algorithm shows an angular power spectrum (dot-dashed line) that is suppressed with respect to that derived from only 90 and 53 GHz inputs (solid line), although they do agree at the 1-sigma level (dark grey and light grey shaded regions for all DMR frequencies and only 90 and 53 GHz respectively). The reason for this behaviour is probably associated with the distinct noise distribution of the 31.5 GHz channel relative to the other two frequencies, which have considerably more consistent noise patterns. This drives the FastICA algorithm to identify two components: one which has the expected CMB frequency scaling and a second with an unphysical (negative) frequency scaling (and interestingly a higher $`S/N`$ ratio) into which part of the CMB power has been aliased.

We conclude that differences in the noise patterns between the DMR channels are a potential source of unphysical components reconstructed with FastICA, and can have a notable effect on the frequency scaling of the reconstructed CMB component and its angular power spectrum. Therefore in what follows we restrict analysis to the 90 and 53 GHz channels only.

## 4 Results – CMB

We have applied FastICA to the 90 and 53 GHz frequency channels. The results are almost identical for the three non-quadratic functions considered. Table 2 summarises the reconstructed CMB frequency scalings and the percentage deviation from the theoretical CMB frequency scaling. Results from the p- and g-functions are consistent in terms of the correct frequency scaling and $`S/N`$ ratio. Results determined using the t-non-quadratic function reveal a non-optimal frequency scaling for the CMB component. Nevertheless the correct normalisation could be derived.

Figure 2 shows the angular power spectrum of the reconstructed CMB component for the p-, g- and t-functions compared with that computed from the minimum-variance noise combination of the 53 and 90 GHz channels. All results are consistent although the t-function shows a small excess of power on all scales. At large $`\mathrm{}`$ this corresponds to a larger noise contribution.

Using the same linear combination coefficients that produce the reconstructed CMB components from the input 90 and 53 GHz maps, it is possible to derive the rms noise properties per pixel for the reconstructions. This noise prescription is used as an input to the likelihood analysis (see Górski et al. 1994 and G96) solving for $`Q_{rmsPS}`$ and $`n`$. Table 3 reports the results of these fits, whilst Figure 3 plots the corresponding likelihood contours.

The results are consistent for p and g while larger $`Q_{rmsPS}`$ and $`n`$ values are recovered for t. Nevertheless, all are consistent with the fit derived from the optimal minimum-variance combination of the 90 and 53 GHz sky maps without foreground corrections applied.

## 5 Adding “priors” on foregrounds

In the previous section we detailed the results of our analysis of FastICA as applied to the 53 and 90 GHz sky maps alone, and stressed the consistency of the cosmological results with those determined in G96 without the application of any galactic foreground correction. The relevance of this comparison can be understood as follows: with only two input sky maps, the FastICA algorithm can only reconstruct two outputs, and we have determined that these correspond to CMB together with an unphysical, noise-related component, i.e. the method is not able to reconstruct the galactic emission at high galactic latitudes. Indeed, we have verified by means of simulations – which include a galactic contribution modelled on a DIRBE 140$`\mu `$m template (scaled to the DMR frequencies according to the correlation coefficients from G96) – that the $`S/N`$ ratio of the galactic emission in the DMR data is too low to allow FastICA to recover it. Only when the noise rms is reduced by a factor of 10 relative to the actual DMR values can some galactic emission be separated by FastICA, although the separation remains less than optimal.

In order to attempt to extract and better separate the CMB signal from that foreground emission which is certainly present in the data, we have proceeded to add additional channels (priors) using foreground templates such as the Haslam map of the diffuse galactic radio emission at 408 MHz (Haslam et al. 1982) and the DIRBE map of galactic diffuse infra-red emission at 140$`\mu `$m. When applying FastICA to the 53 and 90 GHz DMR data in combination with the Haslam and DIRBE templates, we have assumed that the noise contribution to the template maps is small compared to that in the DMR data, which is certainly the case. Since we now work with 4 inputs, FastICA returns 4 derived components. One is again almost consistent with noise ($`S/N`$ lower that 1), one has the CMB frequency scaling and the expected $`S/N`$ ratio while the other two show features respectively of the Haslam and DIRBE maps. It is important to note that FastICA does not completely mix the Haslam and DIRBE signals into a single foreground component, but is able to keep them largely separated (see 6). As before, we fit the CMB component for $`Q_{rmsPS}`$ and $`n`$. Table 4 reports the results: there is a general agreement between our results, especially those determined with the p- and g-functions, and those reported in G96 after galactic emission has been taken into account using the template Haslam and DIRBE maps.

It is clear that p and g almost give the same results while large differences are observed with the t form. This may be due to the low $`S/N`$ ratio of the DMR data: Maino et al. (2002) did not observe such behaviour when working with higher fidelity Planck simulations. Inspection of the extracted components shows that the frequency scalings of the CMB component are almost the same as when only the DMR data were used as inputs (variations are $`\genfrac{}{}{0pt}{}{<}{}\mathrm{\hspace{0.17em}1}\%`$) and that the extracted signals have also the same $`S/N`$ ratio. We also computed the final rms in the extracted maps and compare them with the previous DMR-only case. The rms figures are almost the same for p and g (deviations within 8%) while for the t form we observe a variation in the map rms of about 30%. This implies that the t is unable to adequately separate the underlying components and identify them as CMB and galactic signals.

## 6 Results – Foreground emissions

We have analysed the two additional signal components extracted by FastICA. These appear to closely resemble the Haslam and DIRBE maps, but this is not optimal (the Spearman correlation coefficient $`r_\mathrm{S}90\%`$). Furthermore, visual inspection reveals that some structure which appears in the input Haslam map template appears as “ghosts” (with negative signal) in the reconstructed dust map, and viceversa as shown in Figure 4. Such behaviour almost certainly reflects a more complex spectral signature of the foreground components than the method (with only 4 input data sets) is able to disentangle. Undoubtedly, the significant large angular scale spatial correlation between foregrounds in various regions of the sky and lower fidelity of the foregrounds due to both the $`S/N`$ and large angular smoothing of the DMR data also play a role. In fact, the extent of the signal aliasing depends on the non-quadratic form considered in the analysis, and is particularly evident for the t form. The spectral indices for these galactic independent components derived are not robust to the choice of non-quadratic form and furthermore are not consistent with what expected from previous studies (e.g. G96). However, the CMB component is reconstructed well and remains clean from galactic emission as results in Table 4 show.

Interestingly, this means that FastICA can provide information on foregrounds, not directly from the (spectral) properties of the reconstructed galactic components, but rather from information gleaned in the reconstruction of the CMB component. Each ICA component is obtained by a weighted linear combination of the input data, even though the combination coefficients (actually the rows of the matrix $`𝐖`$) are computed by a non-linear process. This means that we can interpret the derived weights effectively as correlation coefficients between the foreground templates and the DMR data. Indeed, for the previous analysis with the Haslam and DIRBE maps, the recovered values are similar to that found in G96 using an explicit (minimum $`\chi ^2`$) correlation analysis. Furthermore, the consistency of the results between p, g and, to some extent, t, suggest that the method allows a robust estimation of the correlation of the DMR data with the considered templates.

Recently Banday et al. (2003) have derived correlation coefficients between the DMR data and three foreground templates: the Haslam map at 408 MHz, the DIRBE map at 140$`\mu `$m and a new map of H$`\alpha `$ emission (Dickinson, Davies & Davis 2003) that could be used to determine a reliable free-free template. Furthermore when thermal dust and free-free emission are properly accounted for, they found a strong anomalous correlation between data at 31.5 GHz and the DIRBE map. This component is found to be well fitted by a spectral shape of the form $`\nu ^\beta `$ with $`\beta 2.5`$.

We attempt to verify here whether FastICA is able to arrive at similar results, and we therefore extend our analysis to include the H$`\alpha `$ template. An initial consideration is that the non-quadratic forms of the FastICA algorithm are in general sensitive to different features of the astrophysical components, therefore it is possible that a given form may or may not be particularly well matched to the detection of certain signals on the sky. We therefore calibrate which form is best suited to our current studies, by using results from Banday et al. (2003) to construct a fake skies with CMB, foreground emissions (whose correlation coefficients are known a-priori) and $`10000`$ realisations of the DMR instrumental noise. In each run a single DMR channel is considered together with foreground emission templates and each of p, g and t is employed to evaluate the correlation coefficients. These simulations also provide an indication of the uncertainty in the reconstruction process. Simulations show that p and g provide consistent results.

Table 5 reports our estimates of the correlation coefficients at each of the DMR frequencies.

Figure 5 plots our results (filled circles) together with those from Banday et al. (2003) (open circles). There is good general agreement between the two sets of results. FastICA also clearly detects foreground emission correlated with the DIRBE dust map at 140$`\mu `$m and identifies this as the major foreground emission at 31.5 GHz.

In order to further probe the nature of this anomalous dust correlated component, we subtract from the DMR data model 7 of thermal dust emission from Finkbeiner, Davis and Schlegel (1999, hereafter FDS). This comprises two dust components, one at a temperature of 9.6 K with dust emissivity $`\alpha =1.5`$ and the other at 16.4 K with $`\alpha =2.6`$. The results are shown in Table 6.

From Table 6 and Fig. 5 it is clear that, if FDS model 7 provides an accurate prediction for the thermal dust emission at DMR frequencies, then it must account for essentially the entire dust correlation at 90 GHz. At 31.5 GHz a strong anomalous correlation with the DIRBE 140$`\mu `$m dust template persists after removal of the thermal dust model, whilst the results at 53 GHz seem to indicate the presence of both a thermal and anomalous dust correlated component which would have roughly equal contributions at a frequency of $``$ 60 GHz.

We investigate the spectral behaviour of this anomalous dust correlated component by fitting our data points (the correlations coefficients with DIRBE 140$`\mu `$m after FDS model 7 subtraction) by a power law of the form $`𝐀_{\mathrm{norm}}(\nu /\nu _0)^\beta `$. $`𝐀_{\mathrm{norm}}`$ has units of $`\mu `$K/(MJy sr<sup>-1</sup>) and is the amplitude of the emission at the reference frequency $`\nu _0`$ of 31.5 GHz. Table 7 shows the results from this analysis: it is clear that the frequency behaviour is consistent with a spectrum steeper that the typical value of $``$ 2.15 expected for pure free-free emission (although they are still compatible with this kind of emission).

Such conclusions are consistent with those of by Banday et al. (2003), and are particularly satisfactory due to the independence of the two techniques employed, and the fact that the ICA method makes no assumptions about the cosmological signal present in the data.

## 7 Critical Discussion and Conclusion

In this paper we have applied for the first time the FastICA algorithm to real CMB data and specifically to the $`COBE`$-DMR 4-yr data.

One of the requirements for applicability of FastICA is that noise has to be Gaussian and uniform distributed on the sky. Unfortunately, this not the case for the DMR observing strategy. Furthermore part of the original data at 31.5 GHz were removed to eliminate potential contamination from a particular systematic effect (Kogut et al. 1996), resulting in a noise distribution that is notably different than that at 53 and 90 GHz. Therefore, before applying the technique, we have tested the effect of such a non-uniform noise distribution on both the frequency scaling and power spectrum reconstruction of the recovered CMB component. We run $`3000`$ MC simulations and considered three scenarios of different noise distribution in the data: one nominal, one with uniform noise distribution and nominal mean amplitude and the third with uniform distribution and reduced amplitude. Results indicate that the different noise distribution of the 31.5 GHz channel results in the identification of a spurious non-physical component, and the degradation of the frequency scaling of the reconstructed CMB component. Furthermore, the angular power spectrum of the latter has suppressed amplitude with respect to the case in which only the 90 and 53 GHz channels are considered, resulting from the aliasing of CMB power into the non-physical component associated with the noise asymmetry.

Further analysis was therefore restricted to include the 90 and 53 GHz channels only. In this case we obtained a CMB reconstruction with a frequency scaling consistent with the expected theoretical one. We fit the cosmological component with a power-law model parametererised by the rms-quadrupole normalisation $`Q_{rmsPS}`$ and spectral index $`n`$, and find that the results are robust against the choice of non-quadratic function, completely consistent with those derived in G96 from the optimal minimum-variance combination of the 90 and 53 GHz channels, without foreground corrections. This is a very important outcome since equivalent results have been obtained with independent approaches.

Simulations with different noise rms properties have shown that foreground emission in the DMR data has too low a $`S/N`$ ratio to be detected and separated with FastICA. Only with a $`S/N`$ a factor of 10 or better does it become possible to identify foreground emission. However, after adding two additional channels – the Haslam map or radio emission at 408 MHz and the DIRBE map of 140$`\mu `$mgalactic infra-red emission – FastICA is now able to recognise foreground emission in the DMR data and separate it from CMB signal. The CMB frequency scaling for different FastICA forms are consistent with that found without priors on foregrounds. More importantly, we derived values for $`Q_{\mathrm{rms}}`$ and $`n_\mathrm{s}`$ for this new cleaned CMB reconstruction in very good agreement with those reported in G96 after foreground correction.

The situation with the foreground component reconstruction is less satisfying, however. We found that the foreground spectral indices derived from the ICA foreground components are not robust against the choice of the FastICA non-quadratic form and this could be due to several reasons. First of all the different forms tend to preferentially discriminate between different morphological aspects of the foreground emission, with the residuals tending to be mixed together. Secondly, the $`S/N`$ ratio for foregrounds at high galactic latitudes is too low for FastICA to obtain faithful information on these emissions. Finally, and quite realistically, it is likely that the true foreground emission at microwave frequencies is different from our priors which are based on observations at considerably different wavelengths. It is worth noting, however, that Maino et al. (2002) have demonstrated that, for higher $`S/N`$ observations and with an exact correspondence between “priors” and foreground signals, the algorithm was able to properly extract the underlying foreground contributions.

Despite the fact that the spectral behaviour of the derived foreground components is not properly recovered, our results do still allow some information on the non-CMB emission to be determined. Since the ICA components are linear combinations of the input sky maps, the weightings between the various input data used to reconstruct the CMB sky can be interpreted as correlation coefficients between the DMR data at specific frequencies and the foreground templates, in the spirit of the analysis in G96. Indeed, the results in Table 4 are in substantial agreement with this. Following the recent work of Banday et al. (2003) we have further considered an extra channel using the H$`\alpha `$ template provided by Dickinson, Davies & Davis (2002). Using these three templates, FastICA is able to derive the correlation coefficients between the DMR data and foreground templates in very good agreement with their results.We found that after subtraction of the dust thermal emission as accounted for by model 7 of FDS, a strong dust correlated anomalous component remains, and consitutes the strongest foreground emission at 31.5 GHz. This component is found to be well fitted by a power law spectrum with $`\beta 2.56`$.

Blind component separation algorithms such as FastICA are promising and can be applied to real data even if their assumptions are not completely satisfied. The results obtained so far are encouraging and future application to CMB datasets (e.g. $`WMAP`$ and Planck) are foreseen.

## Acknowledgements

It is a pleasure to thank E. Salerno for useful discussion. DM thanks A.J. Banday for warm hospitality at MPA and AJB thanks D. Maino for kind hospitality at OAT in Trieste. We thank Clive Dickinson for providing us with the H$`\alpha `$ template prior to publication.


text:
# Asset trees and asset graphs in financial markets

## I Introduction

Mantegna suggested studying the clustering of companies using the correlation matrix of asset returns Man1 . The correlations are transformed into distances, a subset of which is selected using the minimum spanning tree (MST) criterion. In the resulting tree, the distances are the edges used to connect the nodes, or companies, and thus a taxonomy of the financial market is formed. This method was later studied by Bonanno et al. Gio2 , while other studies on clustering in the financial market are Kull ; Gio ; Kas ; Lal ; Mar , and those specifically on market crashes Lillo ; Lillo2 .

Recently, we have studied the properties of a set of asset trees created using the methodology introduced by Mantegna in jpo ; short ; long , and applied it in the crash context in bali . In these studies, the obtained multitude of trees was interpreted as a sequence of evolutionary steps of a single ‘dynamic asset tree’, and different measures were used to characterize the system, which were found to reflect the state of the market. The economic meaningfulness of the emerging clustering was also discussed and the dynamic asset trees were found to have a strong connection to portfolio optimization.

In this paper, we introduce a modified methodology which, in general, does not not lead to a tree but a graph, or possibly even several graphs that do not need to be inter-connected. Here we limit ourselves to studying only one type of ’dynamic asset graph’, which in terms of its size is compatible and thus comparable with the dynamic asset tree. Although in graph theory a tree is defined as a type of graph, the terms asset graph and asset tree are used here to refer to the two different approaches and as concepts are mutually exclusive and noninterchangeable. The aims of this paper are to introduce this modified approach and to demonstrate some of its similarities and differences with the previous approach. Further considerations of topology and taxonomy of the financial market obtained using dynamic asset graphs are to be presented later.

## II Constructing graphs and trees

In this paper, the term financial market refers to a set of asset price data commercially available from the Center for Research in Security Prices (CRSP) of the University of Chicago Graduate School of Business. Here we will study the split-adjusted daily closure prices for a total of $`N=477`$ stocks traded at the New York Stock Exchange (NYSE) over the period of 20 years, from 02-Jan-1980 to 31-Dec-1999. This amounts a total of 5056 price quotes per stock, indexed by time variable $`\tau =1,2,\mathrm{},5056`$. For analysis and smoothing purposes, the data is divided time-wise into $`M`$ *windows* $`t=1,\mathrm{\hspace{0.17em}2},\mathrm{},M`$ of width $`T`$, where $`T`$ corresponds to the number of daily returns included in the window. Several consecutive windows overlap with each other, the extent of which is dictated by the window step length parameter $`\delta T`$, which describes the displacement of the window and is also measured in trading days. The choice of window width is a trade-off between too noisy and too smoothed data for small and large window widths, respectively. The results presented in this paper were calculated from monthly stepped four-year windows, i.e. $`\delta T=250/1221`$ days and $`T=1000`$ days. We have explored a large scale of different values for both parameters, and the cited values were found optimal jpo . With these choices, the overall number of windows is $`M=195`$.

In order to investigate correlations between stocks we first denote the closure price of stock $`i`$ at time $`\tau `$ by $`P_i(\tau )`$ (Note that $`\tau `$ refers to a date, not a time window). We focus our attention to the logarithmic return of stock $`i`$, given by $`r_i(\tau )=\mathrm{ln}P_i(\tau )\mathrm{ln}P_i(\tau 1)`$ which for a sequence of consecutive trading days, i.e. those encompassing the given window $`t`$, form the return vector $`𝒓_i^t`$. In order to characterize the synchronous time evolution of assets, we use the equal time correlation coefficients between assets $`i`$ and $`j`$ defined as

$$\rho _{ij}^t=\frac{𝒓_i^t𝒓_j^t𝒓_i^t𝒓_j^t}{\sqrt{[𝒓_{i}^{t}{}_{}{}^{2}𝒓_i^t^2][𝒓_{j}^{t}{}_{}{}^{2}𝒓_j^t^2]}},$$
(1)
where $`\mathrm{}`$ indicates a time average over the consecutive trading days included in the return vectors. Due to Cauchy-Schwarz inequality, these correlation coefficients fulfill the condition $`1\rho _{ij}1`$ and form an $`N\times N`$ correlation matrix $`𝐂^t`$, which serves as the basis for graphs and trees discussed in this paper.

For the purpose of constructing asset graphs and asset trees, we define a distance between a pair of stocks. This distance is associated with the edge connecting the stocks and it is expected to reflect the level at which the stocks are correlated. We use a simple non-linear transformation $`d_{ij}^t=\sqrt{2(1\rho _{ij}^t)}`$ to obtain distances with the property $`2d_{ij}0`$, forming an $`N\times N`$ symmetric distance matrix $`𝐃^t`$. Now two alternative approaches may be adopted. The first one leads to asset trees, and the second one to asset graphs. In both approaches the trees (or graphs) for different time windows are not independent of each other, but form a series through time. Consequently, this multitude of trees or graphs is interpreted as a sequence of evolutionary steps of a single *dynamic asset tree* or *dynamic asset graph*.

In the first approach we construct an asset tree according to the methodology by Mantegna Man1 . This approach requires an additional hypothesis about the topology of the metric space, namely, the so-called ultrametricity hypothesis. In practice, it leads to determining the minimum spanning tree (MST) of the distances, denoted $`𝐓^t`$. The spanning tree is a simply connected acyclic (no cycles) graph that connects all $`N`$ nodes (stocks) with $`N1`$ edges such that the sum of all edge weights, $`_{d_{ij}^t𝐓^t}d_{ij}^t`$, is minimum. We refer to the minimum spanning tree at time $`t`$ by the notation $`𝐓^t=(V,E^t)`$, where $`V`$ is a set of vertices and $`E^t`$ is a corresponding set of unordered pairs of vertices, or edges. Since the spanning tree criterion requires all $`N`$ nodes to be always present, the set of vertices $`V`$ is time independent, which is why the time superscript has been dropped from notation. The set of edges $`E^t`$, however, does depend on time, as it is expected that edge lengths in the matrix $`𝐃^t`$ evolve over time, and thus different edges get selected in the tree at different times.

In the second approach we construct asset graphs. This consists of extracting the $`N(N1)/2`$ distinct distance elements from the upper (or lower) triangular part of the distance matrix $`𝐃^t`$, and obtaining a sequence of edges $`d_1^t,d_2^t,\mathrm{},d_{N(N1)/2}^t`$, where we have used a single index notation. The edges are then sorted in an ascending order and form an ordered sequence $`d_{(1)}^t,d_{(2)}^t,\mathrm{},d_{(N(N1)/2)}^t`$. Since we require the graph to be representative of the market, it is natural to build the graph by including only the strongest connections in it. The number of edges to include is, of course, arbitrary. Here we include only $`N1`$ shortest edges in the graph, thus giving $`V^t=\{d_{(1)}^t,d_{(2)}^t,\mathrm{},d_{(N1)}^t\}`$. This is motivated by the fact that the asset tree also consists of $`N1`$ edges, and this choice renders the two methodologies comparable, and possibly even similar to one another. The presented mechanism for constructing graphs defines them uniquely and, consequently, no additional hypotheses about graph topology are required. It is important to note that both the set of vertices $`V^t`$ and the set of edges $`E^t`$ are time dependent, and thus we denote the graph with $`𝐆^t=(V^t,E^t)`$. The choice to include only the $`N1`$ shortest edges in the graph means that the *size* of the graph, defined as the number of its edges, is fixed at $`N1`$. However, the *order* of the graph, defined as the number of its vertices, is not fixed for the graph but varies as a function of time. This is due to the fact that even a small set of vertices may be strongly inter-connected, and thus may use up many of the available edges. This may also lead to the formation of loops in the graph. These aspects are clearly different from the tree approach, where the order is always fixed at $`N`$ and no loops are allowed by definition.

## III Market characterization

We start comparisons between the two approaches by visually examining edge length distribution plots. The three plots are the distribution of (i) distance elements $`d_{ij}^t`$ contained in the distance matrix $`𝐃^t`$ (Figure 1), (ii) distance elements $`d_{ij}`$ contained in the asset (minimum spanning) tree $`𝐓^t`$ (Figure 2), and (iii) distance elements $`d_{ij}`$ contained in the asset graph $`𝐆^t`$ (Figure 3). In all three plots, but most prominently in Figure 1, there appears to be a discontinuity in the distribution between roughly 1986 and 1990. It seems that part has been cut out, pushed to the left and made flatter. This anomaly is a manifestation of Black Monday (October 19, 1987), and its length along the time axis is related to the choice of window width $`T`$ bali ; long .

We can now reconsider the tree and graph construction mechanism described earlier. Starting from the distribution of the $`N(N1)/2`$ distance matrix elements in Figure 1, for the asset tree we pick the shortest $`N1`$ of them, subject to the constraint that all vertices are spanned by the chosen edges. For the purpose of building the graph, however, this constraint is dropped and we pick the shortest elements from the distribution in Figure 1. Therefore, the distribution of graph edges in Figure 3 is simply the left tail of the distribution of distance elements in Figure 1, and it seems that the asset graph rarely contains edges longer than about 1.1, and the largest distance element is $`d_{max}=1.1375`$. In contrast, in the distribution of tree edges in Figure 2 most edges included in the tree seem to come from the area to the right of the value 1.1, and the largest distance element is now $`d_{max}=1.3549`$.

Instead of using the edge length distributions as such, we can characterize the market by studying the location (mean) of the edge length distribution by defining a simple measure, the *mean distance*, as

$$\overline{d}(t)=\frac{1}{N(N1)/2}\underset{d_{ij}^t𝐃^t}{}d_{ij}^t,$$
(2)
where $`t`$ denotes the time at which the tree is constructed, and the denominator is the number of distinct elements in the matrix. We mention in passing that one could instead use the *mean correlation coefficient*, defined as

$$\overline{\rho }(t)=\frac{1}{N(N1)/2}\underset{\rho _{ij}^t𝐂^t}{}\rho _{ij}^t,$$
(3)
which would lead to the same conclusions, as the mean distance and mean correlation coefficient are mirror images of one another and, consequently, it suffices to examine either one of them. In a similar manner, we can characterize the asset tree and the asset graph, which are both simplified networks representing the market, but use only $`N1`$ distance elements $`d_{ij}^t`$ out of the available $`N(N1)/2`$ in the distance matrix $`𝐃^t`$. We define the *normalized tree length* for the asset tree as

$$L_{\text{mst}}(t)=\frac{1}{N1}\underset{d_{ij}^t𝐓^t}{}d_{ij}^t,$$
(4)
and the *normalized graph length* for the asset graph as

$$L_{\text{graph}}(t)=\frac{1}{N1}\underset{d_{ij}^t𝐆^t}{}d_{ij}^t,$$
(5)
where $`t`$ again denotes the time at which the tree or graph is constructed, and $`N1`$ is the number of edges present. The three measures are depicted in Figure 4. Three observations are made. First, all three behave very similarly, which is also reflected by their level of mutual correlation. Pearson’s linear and Spearman’s rank-order correlation coefficients between the mean distance and normalized tree length are 0.98 and 0.92, respectively. The same measures between the mean distance and the normalized graph length are 0.96 and 0.87. Thus, the normalized tree length seems to track the market slightly better. Second, the measures decrease in moving from the mean distance, via the normalized tree length, to the normalized graph length, and their averages are 1.29, 1.12 and 1.00, respectively. Also, the normalized tree length is always higher than the the normalized graph length, and the difference is 0.13 on average. It seems that the asset tee, due to the minimum spanning tree criterion, is forced to “accept” edge lengths that are far less optimal (longer) than the asset graph, resulting in a higher average value for the normalized tree length than for the normalized graph length. Third, the normalized graph length tends to exaggerate the depression caused by the crash, which can be traced back to the graph construction mechanism. Although all of the above measures are very simple, we have earlier studied the normalized tree length and found it to be descriptive of the overall market state. It is closely related to market diversification potential, i.e. the scope of the market to eliminate specific risk of the minimum risk Markowitz portfolio jpo ; short . The fact that the normalized distance and normalized tree length behave so similarly suggests that they are, at least to some extent, interchangeable measures.

Based on how the asset graphs and trees are constructed, the edge length distributions in Figures 2 and 3, and the normalized graph and tree length plots in Figure 4, it is evident that they consist of different edges and, consequently, have different topologies. We can learn about the overall topological differences between the asset graph $`𝐆^t=(V_G^t,E_G^t)`$ and the asset tree $`𝐓^t=(V_T,E_T^t)`$ by studying the overlap of edges present in both as a function of time. The relative overlap is given by the quantity $`\frac{1}{N1}|E_G^tE_T^t|`$ where $``$ is the intersection operator and $`|\mathrm{}|`$ gives the number of elements in the set. As can be see from the plot in Figure 5, the quantity is quite stable over time. On average, the asset graph and asset tree seem to share some 25% of edges.

It is also of interest to study how this overlap of edges changes in the process of determining the asset graph and tree. Assume we use Kruskal’s algorithm to obtain the minimum spanning tree when constructing the asset tree. We would first arrange the edges in nondecreasing order by their length as a list, and select the shortest unexamined edge for inclusion in the tree, with the condition that it does not form a cycle. If it does, we discard it, and move on to the next unexamined edge on the list. Except for the constraint on cycles, the algorithm is identical to the way asset graphs are generated. If we denote the size of graph *in construction* by $`n`$, where $`n=1,2,\mathrm{},N1`$, then at least for small values of $`n`$ asset graphs and asset trees should contain the same set of edges, i.e., $`E(n)_G^t=E(n)_T^t`$ and, therefore, be identical in topology. It is expected that, starting from some value of $`n=n_c`$, the above equality no longer holds. This indicates the formation of the first cycle and, consequently, for all $`nn_c`$ the asset graph and tree are topologically different. This is demonstrated in Figure 6, where the relative overlap of edges, $`\frac{1}{N1}|E_G^tE_T^t|`$, has been plotted as a function of normalized number of edges, $`\frac{n}{N1}`$, and the quantity has been averaged over time. The function decreases rapidly for small values of $`\frac{n}{N1}`$, indicating that for the current set of data with $`N=477`$, only a few edges can be added before the first cycle is formed.

## IV Evolution of asset graphs and trees

The robustness of asset graph and asset tree topology can be studied through the concept of *single-step survival ratio*, defined as the fraction of edges found common in two consecutive graphs or trees at times $`t`$ and $`t1`$ as

$$\sigma (t)=\frac{1}{N1}|E^tE^{t1}|.$$
(6)
As before, $`E^t`$ refers to the set of edges of the graph or the tree at time $`t`$, $``$ is the intersection operator and $`|\mathrm{}|`$ gives the number of elements in the set. Although it has not been explicitly indicated in the definition, $`\sigma (t)`$ is dependent on the two parameters, namely, the window width $`T`$ and the step length $`\delta T`$. Figure 7 shows the plots of single-step survival for both the graph (upper curve) and the tree (lower curve) for $`\delta T=250/1221`$ days and $`T=1000`$ days.

The most evident observation is that the graph seems to have a higher survival ratio than the tree. For the graph, on average, some 94.8% of connections survive, whereas the corresponding number for the tree is about 82.6%. In addition, the fluctuation of the single-step survival ratio, as measured by its standard deviation, is smaller for the graph at 5.3% than for the tree at 6.2%. In general, both curves fluctuate together, meaning that the market events causing rewiring in the graph also cause re-wirings in the tree. This is very clear in the two sudden dips in both curves, which result from the re-wirings related to Black Monday bali . Although both curves fall drastically, the one for the asset graph falls less, indicating that the graph is more stable than the tree also under extreme circumstances, such as market crashes. The higher survival ratio for the asset graph is not caused by any particular choice of parameters, but is reproduced for all examined parametric values. Some indication of the sensitivity of the single-step survival ratio on the window width parameter is given in Table 1. The fact that asset graphs are more stable than asset trees is related to their construction mechanism. The spanning tree constraint practically never allows choosing the shortest available edges for the tree and, consequently, the ensuing structure is more fragile. A short edge between a pair of stocks corresponds to a very high correlation between their returns. This may result from the companies having developed a cooperative relationship, such as a joint venture, or then may be simply incidental. In the first case the created bond between the two stocks is likely to be longer lasting than if the MST criterion forced us to include a weaker bond between two companies.

We can expand the concept of single-step survival ratio to cover survival over several consecutive time steps $`\delta T`$. Whereas the single-step survival ratio was used to study short-term survival, or robustness, of graphs and trees, the *multi-step survival ratio* is used to study their long-term survival. It is defined as

$$\sigma (t,k)=\frac{1}{N1}|E^tE^{t1}\mathrm{}E^{tk+1}E^{tk}|,$$
(7)
where only those connections that have persisted for the whole time period of length $`k\delta T`$ without any interruptions are taken into account. According to this formula, when a bond between two vertices breaks even once within $`k`$ steps and then reappears, it is not counted as a survived connection. A closely related concept is that of graph or tree *half-life* $`t_{1/2}`$, defined as the time in which half the number of initial connections have decayed, i.e., $`\sigma (t,t_{1/2})=0.5`$. The multi-step survival ratio is plotted in Figure 8, where the half-life threshold is indicated by the dashed horizontal line.

The time axis can be divided into two regions based on the nature of the decay process, and these regions are located somewhat differently for the graph and the tree. The precise locations of the regions are, of course, subject to speculation, but for the purpose carrying out fits and analysis they need to be fixed. For the asset graph the first and second regions, discretized according to $`\delta T=\frac{1}{12}`$ year as mentioned before, are defined on the intervals $`(\frac{1}{2},4)`$ and $`(4\frac{1}{12},16\frac{1}{6})`$, respectively, both given in years. Within the first region, the graph exhibits clean exponential decay, as witnessed by the fitted straight line on lin-log scale in the inset of Figure 8. Somewhere in between the two regions there is a cross-over to power-law behavior, which is evident within the second region, resulting in a straight line on the log-log scale plot of the same figure. For the asset tree the regions are defined on $`(\frac{1}{12},1\frac{1}{2})`$ and $`(1\frac{7}{12},11\frac{1}{4})`$. Within the first region the asset tree decays faster than exponentially, as can be verified by comparison with the straight line decay of the graph in the inset. Similarly to the tree, there is a cross-over to a power-law, although the slope is faster than for the graph. If we write the power law decay as $`\sigma (t,k)_tk^\gamma `$, the fits yield for the asset graph $`\gamma 1.39`$, whereas for the asset tree we have $`\gamma 1.19`$.

The finding concerning the slower decay of the asset graph within the first region is fully compatible with the results obtained with single-step survival ratio. Since the graph shows higher survival ratio over a single-step, it is to be expected that, at least in the early time horizon (within the first region), graphs should decay more slowly than trees. The half-lives for both the graph and the tree occur within the first region, and thus it is not surprising that graph half-life is much longer than tree half-life. For the graph, we obtained $`t_{1/2}1.71`$ years, and for the tree $`t_{1/2}0.47`$. Although the half-lives depend on the value of window width $`T`$, the differences between them persists for different parametric values. When measured in years, for window widths of $`T=2`$, $`T=4`$ and $`T=6`$, the corresponding half-lives for the asset tree are 0.22, 0.47 and 0.75 years, whereas for the asset graph they are 0.88, 1.71 and 2.51 years, respectively. Interestingly, the situation seems to be reversed within the second region, where both decay as power-law. Here the higher exponent of $`\gamma `$ for the asset graph indicates that it actually decays faster than the asset tree. This finding could, of course, be influenced by our choice of the window width $`T`$. Explorations with that parameter revealed another interesting phenomenon. Whereas the slope for the asset tree seems to be independent of window width, as discussed in long , that of asset graphs is not. For $`T=500`$, $`T=1000`$ and $`T=1500`$ we obtained for the asset tree the exponents $`\gamma =1.15`$, $`\gamma =1.19`$ and $`\gamma =1.17`$, respectively, which, within the error bars, are to be considered equal long . For the asset graph, we obtain the values of $`\gamma =2.07`$, $`\gamma =1.39`$ and $`\gamma =1.55`$. Although no clear trend can be detected in these values, a matter that calls for further exploration, it is clear that the value of $`\gamma `$ is higher for the asset graph than for the asset tree. Therefore, the asset graph decays more slowly than the asset tree within the first region, while within the second region the situation is just the opposite.

## V Distribution of vertex degrees in asset graphs and asset trees

As the asset graph and asset tree are representative of the financial market, studying their structure can enhance our understanding of the market itself. Recently Vandewalle et al. Van  found scale free behavior for the asset tree in a limited time window. They proposed the distribution of the vertex degrees $`f(k)`$ to follow a power law of the following form

$$f(k)k^\alpha ,$$
(8)
with the exponent $`\alpha 2.2`$. Later, we studied this phenomenon further with a focus on asset tree dynamics long . We found that the asset tree exhibits, most of the time, scale free properties with a rather robust exponent $`\alpha 2.1\pm 0.1`$ during times of normal stock market operation. In addition, within the error limits, the exponent was found to be constant over time. However, during crash periods when the asset tree topology is drastically affected, the exponent changes to $`\alpha 1.8\pm 0.1`$, but nevertheless the asset tree maintains its scale free character. The interesting question is whether asset graphs also display similar scale free behavior and if so, are there are differences in the value of the exponent. As Figures 9 and 10 make clear, for the observed data does not fit as well with scale free behavior for the asset graph as it does for the asset tree. The obtained average value for the exponent of the asset graph is significantly lower, i.e. $`\alpha 0.9\pm 0.1`$. In addition, the exponent for the aseet graph varies less as a function of time and does not show distinctively different behavior between normal and crash markets.

In the case of the asset tree, there were sometimes clear outliers, as one node typically had a considerably higher vertex degree than the power-law scaling would predict. This outlier was used as a *central node*, a reference against which some tree properties were measured. However, the fact that one node often had “too high” vertex degree provided further support for using one of the nodes as the center of the tree, as discussed in detail in long . In the dynamic asset graph, these types of outliers are not present. This observation merely reflects upon the differences between the topologies produced by the two different methodologies but does not, as such, rule out the possibility of using one of the nodes as a central node seli .

In long , we estimated the overall goodness of power-law fits by calculating the $`R^2`$ coefficient of determination, a measure which indicates the fraction of the total variation explained by the least-squares regression line. Averaged over all the time windows, we obtained the values $`R^20.93`$ and $`R^20.86`$, with and without outliers excluded, respectively. Since there were no outliers in the data for asset graphs, it was used as such to give an average of $`R^20.75`$. This indicates that scale free behavior is not evident in this case.

## VI Summary and conclusion

In summary, we have introduced the concept of dynamic asset graph and compared some of its properties to the dynamic asset tree, which we have studied recently. Comparisons between edge length distributions reveal that the asset tree, due to the minimum spanning tree criterion, is forced to “accept” edge lengths that are far less optimal (longer) than the asset graph. This results in a higher average value for the normalized tree length than for the normalized graph length although, in general, they behave very similarly. However, the latter tends to exaggerate market anomalies and, consequently, the normalized tree length seems to track the market better. The asset graph was also found to exhibit clearly higher single-step survival ratio than the asset tree. This is understandable, as the spanning tree criterion does not allow the shortest connections, which were conjectured to have the longest lives, to be included in the tree and, thus, their omission leads to a more fragile structure. This is also witnessed by studying the multi-step survival ratio, where it was found that in the early time horizon the asset graph shows exponential decay, but the asset tree decays faster than exponential. Later on, however, both decay as a power-law, but here the situation is reversed and the asset tree decays more slowly than the asset graph. We also studied the vertex degree distributions produced by the two alternative approaches. Earlier we have found asset trees to exhibit clear scale free behavior, but for the asset graph scale free behavior is not so evident. Further, the values obtained for the scaling exponent are very different from our earlier studies with asset trees.

###### Acknowledgements.

J.-P. O. is grateful to European Science Foundation for REACTOR grant to visit Hungary, the Budapest University of Technology and Economics for the warm hospitality. Further, the role of Harri Toivonen at the Department of Accounting, Helsinki School of Economics, is acknowledged for carrying out CRSP database extractions. J.-P. O. is also grateful to the Graduate School in Computational Methods of Information Technology (ComMIT), Finland. The authors are also grateful to R. N. Mantegna for very useful discussions and suggestions. This research was partially supported by the Academy of Finland, Research Center for Computational Science and Engineering, project no. 44897 (Finnish Center of Excellence Program 2000-2005) and OTKA (T029985).