text:
# Quantum affine reflection algebras of type 𝑑_𝑛⁽¹⁾ and reflection matrices

## 1. Introduction

Solutions to the reflection equation (also known as the boundary Yang-Baxter equation) , called reflection matrices, are required for the construction of quantum integrable models with boundaries . The reflection equation at the boundary is similar to the Yang-Baxter equation in the bulk, the study of which led to the axiomatization of quantum groups . Intertwiners of these quantum groups are solutions of the Yang-Baxter equation (called R-matrices) .

An analogous method to find reflection matrices using intertwiners of reflection equation algebras  embedded as coideal subalgebras in the corresponding Yangian or quantum affine algebra has recently been developed . While there is a general construction of these subalgebras in terms of the universal R-matrix , it is of practical importance to have a simple set of generators for them. Such simple generators have been found for coideal subalgebras of Yangians  and for coideal subalgebras of quantum affine algebras , which we refer to as “quantum affine reflection algebras”. With these expressions it is easy to solve the intertwining condition, which has already led to several new reflection matrices .

The purpose of this letter is to demonstrate the practical utility of these quantum affine reflection algebras by using them to derive hitherto unknown reflection matrices corresponding to the vector representation of $`so(2n)`$. These reflection matrices are needed, for example, to describe the reflection of solitons in $`d_n^{(1)}`$ affine Toda theory off an integrable boundary .

## 2. Quantum affine algebras

In this section we summarize some essential facts about quantum affine algebras. More details can be found in, for example, .

Let $`\widehat{𝔤}`$ be an affine Lie algebra with generalized Cartan matrix $`(a_{ij})_{i,j=0,\mathrm{},n}`$. The quantum affine algebra $`U_h(\widehat{𝔤})`$ is the unital associative algebra over $`[[h]]`$ with generators $`x_i^+,x_i^{},h_i`$, $`i=0,\mathrm{},n`$ and relations

$$[h_i,h_j]=0,[h_i,x_j^\pm ]=\pm a_{ij}x_j^\pm ,$$

(1)
$$[x_i^+,x_j^{}]=\delta _{ij}\frac{q_i^{h_i}q_i^{h_i}}{q_iq_i^1},$$

(4)
$$\underset{k=0}{\overset{1a_{ij}}{}}(1)^k\left[\begin{array}{c}1a_{ij}\\ k\end{array}\right]_{q_i}(x_i^\pm )^kx_j^\pm (x_i^\pm )^{1a_{ij}k}=0ij.$$

Here $`\left[\begin{array}{c}a\\ b\end{array}\right]_q`$ are the $`q`$-binomial coefficients. We have defined $`q_i=e^{d_ih}`$ where the $`d_i`$ are such that $`d_ia_{ij}`$ is a symmetric matrix.

The Hopf algebra structure of $`U_h(\widehat{𝔤})`$ is given by the comultiplication $`\mathrm{\Delta }:U_h(\widehat{𝔤})U_h(\widehat{𝔤})U_h(\widehat{𝔤})`$ defined by

$`\mathrm{\Delta }(h_i)`$ $`=h_i1+1h_i,`$
(5) $`\mathrm{\Delta }(x_i^\pm )`$ $`=x_i^\pm q_i^{h_i/2}+q_i^{h_i/2}x_i^\pm ,`$
and the antipode $`S`$ and counit $`ϵ`$ defined by

(6)
$$S(h_i)=h_i,S(x_i^\pm )=q_i^1x_i^\pm ,ϵ(h_i)=ϵ(x_i^\pm )=0.$$

For any $`(n+1)`$-tuple $`\{\sigma _0,\mathrm{},\sigma _n\}`$ of invertible constants in $`[[h]]`$ there is a Hopf-algebra automorphism $`\sigma `$ of $`U_h(\widehat{𝔤})`$ defined by

(7)
$$\sigma :x_i^\pm \sigma _i^{\pm 1}x_i^\pm ,h_ih_i.$$

Given some finite-dimensional $`U_h(\widehat{𝔤})`$-module $`V`$ with representation map $`\pi :U_h(\widehat{𝔤})\text{End}(V)`$ one can define a spectral parameter dependent representation $`\pi _x:U_h(\widehat{𝔤})\text{End}(V_x)`$ by

(8)
$$\pi _x=\pi \sigma _x$$

where $`\sigma _x`$ is an automorphism of the form (7) with the choice $`\sigma _0=x`$, $`\sigma _i=1,i=1,\mathrm{},n`$. The tensor product modules $`V_{x_1}V_{x_2}\mathrm{}V_{x_n}`$ are irreducible for generic values for the spectral parameters $`x_i`$.

Let $`\stackrel{ˇ}{R}(x)`$ be the unique module homomorphism $`\stackrel{ˇ}{R}(x):V_xV_1V_1V_x`$, i.e., let $`\stackrel{ˇ}{R}(x)`$ be a solution of the intertwining equation

(9)
$$\stackrel{ˇ}{R}(x)(\pi _x\pi _1)(\mathrm{\Delta }(Q))=(\pi _1\pi _x)(\mathrm{\Delta }(Q))\stackrel{ˇ}{R}(x),\text{for all }QU_h(\widehat{𝔤}).$$

Then $`\stackrel{ˇ}{R}(x)`$ automatically satisfies  the Yang-Baxter equation

(10)
$$(\stackrel{ˇ}{R}(y)1)(1\stackrel{ˇ}{R}(xy))(\stackrel{ˇ}{R}(x)1)=(1\stackrel{ˇ}{R}(x))(\stackrel{ˇ}{R}(xy)1)(1\stackrel{ˇ}{R}(y)).$$

## 3. Quantum affine reflection algebras

In this paper we are concerned with subalgebras $`_h(\widehat{𝔤})U_h(\widehat{𝔤})`$ generated by

(11)
$$\widehat{Q}_j=q_j^{h_j/2}\left(x_j^++x_j^{}\right)+\widehat{ϵ}_j\left(q_j^{h_j}1\right),j=0,\mathrm{},n,$$

for certain choices of the parameters $`\widehat{ϵ}_j[[h]]`$. We refer to these algebras as quantum affine reflection algebras. It was shown in  how they are related to Sklyanin’s reflection equation algebras . The algebras $`_h(\widehat{𝔤})`$ are left coideal subalgebras of $`U_h(\widehat{𝔤})`$ because

(12)
$$\mathrm{\Delta }(\widehat{Q}_j)=\widehat{Q}_j1+q_j^{h_j}\widehat{Q}_jU_h(\widehat{𝔤})_h(\widehat{𝔤}).$$

Generically, the irreducible evaluation modules $`V_x`$ of $`U_h(\widehat{𝔤})`$ and their tensor products are also irreducible as modules of $`_h(\widehat{𝔤})`$.

A rescaling automorphism $`\sigma `$ of the form (7) maps the subalgebra $`_h(\widehat{𝔤})`$ to a different subalgebra of $`U_h(\widehat{𝔤})`$ which is isomorphic to $`_h(\widehat{𝔤})`$ as an algebra and as a $`U_h(\widehat{𝔤})`$-coideal. In particular the sign of $`\widehat{ϵ}_j`$ in (11) can be changed by such an automorphism with $`\sigma _j=1`$.

Let us assume that for some constant $`\eta `$ there exists an intertwiner $`K(x):V_{\eta x}V_{\eta /x}`$, i.e., a matrix $`K(x)`$ satisfying

(13)
$$K(x)\pi _{\eta x}(\widehat{Q})=\pi _{\eta /x}(\widehat{Q})K(x),\widehat{Q}_h(\widehat{𝔤}).$$

Then $`K(x)`$ is automatically a solution of the reflection equation

(14)
$$(1K(y))\stackrel{ˇ}{R}(xy)(1K(x))\stackrel{ˇ}{R}(x/y)=\stackrel{ˇ}{R}(x/y)(1K(x))\stackrel{ˇ}{R}(xy)(1K(y))$$

This is so because both sides of the equation are intertwiners between the irreducible modules $`V_{\eta x}V_{\eta y}`$ and $`V_{\eta /x}V_{\eta /y}`$. By Schur’s lemma the two sides are proportional. That they are equal then follows from the fact that they have the same determinant and are equal at $`h=0`$.

## 4. Vector representation of $`d_n^{(1)}`$

We now specialize to the case of the $`2n`$-dimensional representation of $`\widehat{𝔤}=d_n^{(1)}=\widehat{so(2n)}`$. We choose $`d_i=1`$ in (2). The $`U_h(\widehat{𝔤})`$ generators are represented by

(15) $`\pi (x_j^+)`$ $`=E_{j,j+1}E_{\overline{j+1},\overline{j}},`$ $`\pi (x_j^{})`$ $`=E_{j+1,j}E_{\overline{j},\overline{j+1}},`$
$`\pi (h_j)`$ $`=E_{j,j}E_{j+1,j+1}E_{\overline{j},\overline{j}}+E_{\overline{j+1},\overline{j+1}},`$
(16) $`\pi (x_n^+)`$ $`=E_{n1,\overline{n}}E_{n,\overline{n1}},`$ $`\pi (x_n^{})`$ $`=E_{\overline{n},n1}E_{\overline{n1},n},`$
$`\pi (h_n)`$ $`=E_{n,n}E_{\overline{n},\overline{n}}E_{\overline{n1},\overline{n1}}+E_{n1,n1},`$
(17) $`\pi (x_0^+)`$ $`=E_{\overline{2},1}E_{\overline{1},2},`$ $`\pi (x_0^{})`$ $`=E_{1,\overline{2}}E_{2,\overline{1}},`$
$`\pi (h_0)`$ $`=E_{\overline{2},\overline{2}}E_{1,1}E_{2,2}+E_{\overline{1},\overline{1}},`$
where $`E_{j,k}`$ is the matrix with a $`1`$ in the $`j`$-th row and the $`k`$-th column and all other entries $`0`$, and

(18)
$$\overline{j}=2n+1j.$$

The intertwining equation (9) was solved by Jimbo  to give the trigonometric R-matrix

$`\stackrel{ˇ}{R}(x)=`$ $`(xq^2)(x\xi ){\displaystyle \underset{a\overline{a}}{}}E_{a,a}E_{a,a}+(x1)(x\xi )/q{\displaystyle \underset{ab,\overline{b}}{}}E_{a,b}E_{b,a}`$
(19) $`+(1q^2)(x\xi )\left({\displaystyle \underset{a<b,a\overline{b}}{}}+x{\displaystyle \underset{a>b,a\overline{b}}{}}\right)E_{a,a}E_{b,b}`$
$`+{\displaystyle \underset{ab}{}}b_{ab}(x)E_{\overline{b},a}E_{b,\overline{a}},`$
where $`\xi =q^{22n}`$ and

(20)
$$b_{ab}(x)=\{\begin{array}{cc}(q^2x\xi )(x1)\hfill & \text{for }a=b\hfill \\ (q^21)(\xi q^{b^{}a^{}}(x1)\delta _{a,\overline{b}}(x\xi ))\hfill & \text{for }a<b\hfill \\ (q^21)x(q^{b^{}a^{}}(x1)\delta _{a,\overline{b}}(x\xi ))\hfill & \text{for }a>b\hfill \end{array},$$

with

(21)
$$a^{}=\{\begin{array}{cc}a+1/2\hfill & \text{if }an\hfill \\ a1/2\hfill & \text{if }a>n\hfill \end{array}$$

A Hopf-algebra automorphism $`\sigma `$ of the form (7) corresponds to a change of basis and a change of spectral parameter, i.e.,

(22)
$$\pi _x(\sigma (Q))=\mathrm{\Sigma }^1\pi _{\zeta x}(Q)\mathrm{\Sigma }QU_h(\widehat{𝔤}),$$

where

$`\zeta =\sigma _0\sigma _1(\sigma _2\mathrm{}\sigma _{n2})^2\sigma _{n1}\sigma _n,`$
(23) $`\mathrm{\Sigma }=\text{diag}(\mathrm{\Sigma }_1,\mathrm{},\mathrm{\Sigma }_{2n})`$
$`\mathrm{\Sigma }_j=(\sigma _j\mathrm{}\sigma _{n2})^1,\mathrm{\Sigma }_{\overline{j}}=\sigma _j\mathrm{}\sigma _n,j=1,\mathrm{},n2`$
$`\mathrm{\Sigma }_{n1}=1,\mathrm{\Sigma }_n=\sigma _{n1},\mathrm{\Sigma }_{\overline{n}}=\sigma _n,\mathrm{\Sigma }_{\overline{n1}}=\sigma _{n1}\sigma _n.`$
The R-matrix satisfies

(24)
$$(\mathrm{\Sigma }^1\mathrm{\Sigma }^1)\stackrel{ˇ}{R}(x)(\mathrm{\Sigma }\mathrm{\Sigma })=\stackrel{ˇ}{R}(x).$$

## 5. Reflection matrix for $`d_n^{(1)}`$

It is now straightforward to solve the intertwining equation $`(\text{13})`$ for the vector representation of $`d_n^{(1)}`$ to obtain the solutions of the reflection equation (14). One finds that a solution exists only if

(25) $`\eta `$ $`=\pm (1)^nq^{1n},`$ $`\widehat{ϵ}_j`$ $`=\pm i(q^{1/2}q^{1/2})^1,j=0,\mathrm{},n.`$
The signs of the $`\widehat{ϵ}_j`$’s can be changed by a rescaling automorphism, as discussed above, and the sign of $`\eta `$ can be changed by $`xx`$. We thus now choose the upper signs in (25). We write the solution as

(26)
$$K(x)=K_{ab}(x)E_{a,b}.$$

We only need to give the entries on and above the antidiagonal because the others are determined by the symmetry

(27)
$$K_{\overline{a}\overline{b}}(x)=(1)^n\stackrel{~}{K}_{ab}(1/x),$$

where the tilde indicates the change $`q1/q`$. The entries are

(28)
$$K_{ab}(x)=\{\begin{array}{cc}i^{\underset{¯}{a}\underset{¯}{b}}q^{b^{}n}k(x)\hfill & \text{for }a+b2n\text{ and }ab,\overline{b}\hfill \\ i^{\underset{¯}{a}\underset{¯}{b}}q^{b^{}n}(k(x)+(1)^n\stackrel{~}{k}(1/x))(q+1)^1\hfill & \text{for }b=\overline{a}\text{ and }an\hfill \\ i^{\underset{¯}{a}\underset{¯}{b}}q^{b^{}n}k(x)(1)^nq^{1/2}k(x)\hfill & \\ +((1)^nq^{(n1)/2}+q^{(1n)/2})\hfill & \text{for }b=a\hfill \end{array}$$

where

(29)
$$\underset{¯}{a}=\{\begin{array}{cc}a\hfill & \text{for }an\hfill \\ \overline{a}\hfill & \text{for }a>n,\hfill \end{array}$$

the notation $`a^{}`$ was defined in (21), and

(30)
$$k(x)=q^{n/21}(1x).$$

Using the invariance (24) of the R-matrix with respect to the automorphisms (7) we can write an $`n`$-parameter family of solutions $`K_\sigma ^+`$ of the reflection equation (14)

(31)
$$K_\sigma ^+(x)=\mathrm{\Sigma }^1K(x)\mathrm{\Sigma }.$$

A second family $`K_\sigma ^{}(x)`$ of solutions is obtained by choosing the opposite sign for $`\eta `$ in (25), leading to $`K_\sigma ^{}(x)=K_\sigma ^+(x)`$. Furthermore, multiplying any of these solutions by an arbitrary function of $`x`$ gives another solution because the reflection equation is linear.

## 6. Discussion

The main new result of this letter is the hitherto unknown trigonometric reflection matrices (28) for the vector representation of $`so(2n)`$. It proves the practical utility of the quantum affine reflection algebras with their simple generators (11) and makes clear the possibility of generalization to other representations and other algebras, which is the topic of continuing research .

In order to use the reflection matrices (28) to describe soliton reflection in $`d_n^{(1)}`$ affine Toda field theory they need to be transformed to the principal gradation. This is achieved by an automorphism $`\sigma `$ of the form (7) with $`\sigma _0=1,\sigma _i=x,i=1,\mathrm{},n`$. Furthermore the normalisation of the reflection matrix has to be fixed by the requirements of unitarity and crossing symmetry, as was done for $`a_2^{(1)}`$ in .

An equation very similar to the classical limit of our intertwining condition (13) appeared in  in the context of constructing classically integrable boundary conditions in affine Toda field theory.

Coideal subalgebras of non-affine quantized enveloping algebras $`U_h(𝔤)`$ have been studied  in the context of quantum symmetric spaces and $`q`$-orthogonal polynomials. Expressions similar to (11) appear there.

The properties of quantum affine reflection algebras and their representation theory deserve further study. For example a knowledge of how irreducible representations of $`U_h(\widehat{𝔤})`$ branch into irreducible representations of $`_h(\widehat{𝔤})`$ will have physical applications to the spectrum of boundary states in integrable models.


format: text

text:
# 1 Introduction

## 1 Introduction

The decays of charm mesons are currently a new source of information for the study of light meson spectroscopy, with the advantages of having well defined initial state (the $`D`$ meson, a $`0^{}`$ state with defined mass). This new information is complementary to that from scattering experiments and can be particularly relevant to the understanding of the scalar sector, due to the well-known difficulties of scattering experiments to observe light and broad scalar resonances.

Here we present an overview of the results we obtained analysing the decays $`D_s^+\pi ^{}\pi ^+\pi ^+`$ , $`D^+\pi ^{}\pi ^+\pi ^+`$  and $`D^+K^{}\pi ^+\pi ^+`$ , using data collected in 1991/92 by Fermilab experiment E791 from 500 GeV/c $`\pi ^{}`$–nucleon interactions. For details see .

For the $`D^+\pi ^{}\pi ^+\pi ^+`$ decays, we find that a model with only known $`\pi \pi `$ resonances plus a non-resonant (NR) channel is not able to describe the data adequately. Thus we include a new amplitude in our fit function and find strong evidence for the presence of a light and broad scalar resonance, the $`\sigma (500)`$. The channel involving this scalar meson is responsible for half of the decay rate. We measure the mass and the width of this scalar meson to be $`478_{23}^{+24}\pm 17`$ MeV/$`c^2`$ and $`324_{40}^{+42}\pm 21`$ MeV/$`c^2`$, respectively. Other experiments have presented controversial evidence for this low-mass $`\pi \pi `$ resonance in partial wave analyses , with ambiguous interpretations of the characteristics of such particles.

We found a similar situation for the $`D^+K^{}\pi ^+\pi ^+`$ analysis. When we include all known $`K\pi `$ resonant channels plus a NR contribution, we find that this model is not able to describe the data. By including an extra scalar resonant state, with unconstrained mass and width, we obtain a fit which is substantially superior to that without this state. The values for its mass and width are found to be $`797\pm 19\pm 42`$ MeV/c2 and $`410\pm 43\pm 85`$ MeV/c2 respectively. We refer to this state as the $`\kappa `$. We also obtain new measurements for the mass and the width of the $`K_0^{}(1430)`$ resonance.

From our analysis of $`D_s^+\pi ^{}\pi ^+\pi ^+`$ decays we found that the isoscalar intermediate states are dominant in this decay as well. The largest contribution to this final state comes from the decay involving the scalar meson $`f_0(980)`$ whose nature is a long-standing puzzle. It has been described as a $`q\overline{q}`$ state , a $`K\overline{K}`$ molecule, a glueball, and a multiquark state . We obtain new measurements for the $`f_0(980)`$ and $`f_0(1370)`$ masses and widths.

To obtain these results, we had to introduce a new approach in Dalitz plot analysis in order to extract the mass and width of the scalar resonances by using them as floating parameters in the fit. We begin this paper presenting the general method, applied in the $`D^+K^{}\pi ^+\pi ^+`$ Dalitz-plot analysis, then we discuss the $`D_s^+\pi ^{}\pi ^+\pi ^+`$ and $`D^+\pi ^{}\pi ^+\pi ^+`$ studies using the same procedure.

## 2 The $`D^+K^{}\pi ^+\pi ^+`$ Dalitz-plot Analysis

From the original $`2\times 10^10`$ events collected by E791, and after reconstruction and selection criteria, we obtained the $`D^+K^{}\pi ^+\pi ^+`$ sample shown in Figure 1(a). The filled area represents the level of background; besides the combinatorial, the other main source of background comes from the reflection of the decay $`D_s^+K^{}K^+\pi ^+`$ (through $`\overline{K}^{}K^+`$ and $`\varphi \pi ^+`$). The cross-hatched region contains the events selected for the Dalitz-plot analysis. There are 15090 events in this sample, of which 6% are background.

Figure 1(b) shows the Dalitz-plot for these events. The two axes are the squared invariant-mass combinations for $`K\pi `$, and the plot is symmetrized with respect to the two identical pions. The plot presents a rich structure, where we can observe the clear bands from $`\overline{K}^{}(890)\pi ^+`$, and an accumulation of events at the upper edge of the diagonal, due to heavier resonances. To study the resonant substructure, we perform an unbinned maximun-likelihood fit to the data, with probability distribution functions (PDF’s) for both signal and background sources. In particular, for each candidate event, the signal PDF is written as the square of the total physical amplitude $`𝒜`$ and it is weighted for the acceptance across the Dalitz plot (obtained by Monte Carlo (MC)) and by the level of signal to background for each event, as given by the line shape of Figure 1(a). The background PDF’s (levels and shapes) are fixed for the Dalitz-plot fit, according to MC and data studies.

We begin describing our first approach to fit the data, which represents the conventional Dalitz-plot analysis including the known $`K\pi `$ resonant amplitudes ($`𝒜_n,n1`$), plus a constant non-resonant contribution. The signal amplitude, a Breit-Wigner parametrization times a Blatt-Weisskopf damping factor, is described in reference .

Using this model with well-known resonances (Model A), we find contributions from the following channels: the non-resonant, responsible for more than 90% of the decay rate, followed by $`\overline{K}_0^{}(1430)\pi ^+`$, $`\overline{K}^{}(892)\pi ^+`$, $`\overline{K}^{}(1680)\pi ^+`$ and $`\overline{K}_2^{}(1430)\pi ^+`$. The decay fractions and relative phases are shown in Table 1. These values are in accordance with previous results from E691  and E687 . We thus confirm a high non-resonant contribution according to this model, which is totally unusual in $`D`$ decays. Besides, there is an important destructive interference pattern, since all fractions add up to 140 %.

To evaluate the fit quality, we compute a $`\chi ^2`$ from binned, two-dimensional distributions of data and decay model events. The $`\chi ^2`$ comes from the differences in the binned numbers of events between the data and the model (from a fast MC simulation). We obtain $`\chi ^2/\nu =2.7`$ ($`\nu `$ being the number of degrees of freedom), with a corresponding confidence level (CL) of $`10^{11}`$. In Figure 2(a) we show the $`K\pi `$ low and high squared-mass projections for data (error bars) and model (solid line). The discrepancies are evident in the very low-mass region for $`m^2(K\pi )_{\mathrm{𝑙𝑜𝑤}}`$ and near 2.5 (GeV/c2)2 for $`m^2(K\pi )_{\mathrm{ℎ𝑖𝑔ℎ}}`$. These regions of disagreement are the same observed previously by E687 . We thus conclude that a model with the known $`K\pi `$ resonances, plus a non-resonant amplitude, is not able to describe the $`D^+K^{}\pi ^+\pi ^+`$ Dalitz plot satisfactorily.

A similar pattern – bad fit quality with large NR fraction – is found in the analysis of the decay $`D^+\pi ^{}\pi ^+\pi ^+`$ when allowing only the established $`\pi \pi `$ resonances  . There we find that the inclusion of an extra scalar resonance improves the fit substantially, giving strong evidence for the $`\sigma (500)`$. See the section on $`D^+\pi ^{}\pi ^+\pi ^+`$ below. Thus, we are led to try an extra scalar resonance in our fit model here.

This second fit model, Model B, is constructed by the inclusion of an extra scalar state, with unconstrained mass and width. For consistency, the mass and width of the other scalar state, the $`K_0^{}(1430)`$, are also free parameters of the fit. We adopt a better description for these scalar states by introducing gaussian-type form-factors  to take into account the finite size of the decaying mesons. Two extra floating parameters are the meson radii $`r_D`$ and $`r_R`$ introduced for these meson sizes.

Using this model, we obtain the values of $`797\pm 19\pm 42`$ MeV/c2 for the mass and $`410\pm 43\pm 85`$ MeV/c2 for the width of the new scalar state (first error statistical, second error systematic), referred to here as the $`\kappa `$. The values of mass and width obtained for the $`K_0^{}(1430)`$ are respectively $`1459\pm 7\pm 6`$ MeV/c2 and $`175\pm 12\pm 12`$ MeV/c2, appearing heavier and narrower than presented by the PDG . The decay fractions and relative phases for Model B, with statitical errors, are given in Table 1. Compared to the results of Model A (without $`\kappa `$), the non-resonant mode drops from over 90% to 13 $`\pm `$ 6%. The $`\kappa \pi ^+`$ state is now the dominant channel with decay fraction about 50%. The meson radii $`r_D`$ and $`r_R`$ are found to be respectively $`5.0\pm 0.5`$ GeV-1 and $`1.6\pm 1.3`$ GeV-1, in agreement with values used by other groups .

Moreover, the fit quality of Model B is substantially superior to that of Model A. The $`\chi ^2/\nu `$ is now 0.73 with a CL of 95%. The very good agreement between the model and the data can be seen in the projections of Figure 2(b). A number of studies were done to check these results. They are described in detail in reference .

## 3 The $`D_s^+\pi ^{}\pi ^+\pi ^+`$ Results

In Figure 3 we show the $`\pi ^{}\pi ^+\pi ^+`$ invariant mass distribution after reconstruction and selection criteria  for the sample collected by E791. Besides combinatorial background, reflections from the decays $`D^+K^{}\pi ^+\pi ^+`$ , $`D^0K^{}\pi ^+`$ (plus one extra $`\pi ^+`$ track) and $`D_s^+\eta ^{}\pi ^+,\eta ^{}\rho ^0(770)\gamma `$ are all taken into account. The hatched regions in Figure 3 show the samples used for the Dalitz-plot analyses. There are 1686 and 937 candidate events for $`D^+`$ and $`D_s^+`$ respectively, with a signal to background ratio of about 2:1. The Dalitz plots for these events are shown in Figure 4, the axes corresponding to the two $`\pi ^{}\pi ^+`$ invariant-masses squared.

For the $`D_s^+\pi ^{}\pi ^+\pi ^+`$ events in Figure 4(a), the signal amplitude includes all channels with well-known dipion resonances : $`\rho ^0(770)\pi ^+`$, $`f_0(980)\pi ^+`$, $`f_2(1270)\pi ^+`$, $`f_0(1370)\pi ^+`$, $`\rho ^0(1450)\pi ^+`$ and the non-resonant, assumed constant across the Dalitz plot.

For the $`f_0(980)\pi ^+`$ amplitude, instead of a simple Breit-Wigner form, we use a coupled-channel Breit-Wigner function ,

$$BW_{f_0(980)}=\frac{1}{m_{\pi \pi }^2m_0^2+im_0(\mathrm{\Gamma }_\pi +\mathrm{\Gamma }_K)},$$
(1)
$$\mathrm{\Gamma }_\pi =g_\pi \sqrt{m_{\pi \pi }^2/4m_\pi ^2},\mathrm{\Gamma }_K=\frac{g_K}{2}\left(\sqrt{m_{\pi \pi }^2/4m_{K^+}^2}+\sqrt{m_{\pi \pi }^2/4m_{K^0}^2}\right).$$
(2)
The $`D_s^+\pi ^{}\pi ^+\pi ^+`$ Dalitz plot is fit to obtain not only the decay fractions and phases of the possible sub-channels, but also the parameters of the $`f_0(980)`$ state, $`g_\pi `$, $`g_K`$, and $`m_0`$, as well as the mass and width of the $`f_0(1370)`$. The other resonance masses and widths are taken from the PDG. The resulting fractions and phases are given in Table 2. The measured $`f_0(980)`$ parameters are $`m_0=977\pm 3\pm 2`$ MeV/c2, $`g_\pi =`$ 0.09 $`\pm `$ 0.01 $`\pm `$ 0.01 and $`g_K=`$ 0.02 $`\pm `$ 0.04 $`\pm `$ 0.03. By fitting the Dalitz plot using a simple Breit-Wigner function, for the $`f_0(980)`$ we find $`m_0=975\pm 3`$ MeV/c2 and $`\mathrm{\Gamma }_0=44\pm 2\pm 2`$ MeV/c2, and the results for fractions and phases are indistinguishable.

The confidence level of the fit for $`D_s^+\pi ^{}\pi ^+\pi ^+`$ is 35%. In Figure 5 we show the $`\pi ^{}\pi ^+`$ mass-squared projections for data (points) and model (solid lines, from fast-MC).

## 4 The $`D^+\pi ^{}\pi ^+\pi ^+`$ Results

In a first approach, we try to fit the $`D^+\pi ^{}\pi ^+\pi ^+`$ Dalitz plot of Figure 4(b) with the same amplitudes used for the $`D_s^+\pi ^{}\pi ^+\pi ^+`$ analysis. Using this model, the non-resonant, the $`\rho ^0(1450)\pi ^+`$, and the $`\rho ^0(770)\pi ^+`$ amplitudes are found to dominate, as shown in Table 3, and in agreement with previous reported analyses . However, this model does not describe the data satisfactorily, especially at low $`\pi ^{}\pi ^+`$ mass squared, as can be seen from Fig. 6(a). The $`\chi ^2/\nu `$ obtained from the binned Dalitz plot for this model is 1.6, with a CL less than $`10^5`$.

To investigate the possibility that another $`\pi ^{}\pi ^+`$ resonance contributes to the $`D^+\pi ^{}\pi ^+\pi ^+`$ decay, we add an extra scalar resonance amplitude to the signal PDF, with mass and width as floating parameters in the fit.

We find that this model improves our fit substantially. The mass and the width of the extra scalar state are found to be $`478_{23}^{+24}\pm 17`$ MeV/$`c^2`$ and $`324_{40}^{+42}\pm 21`$ MeV/$`c^2`$, respectively. Refering to this state as the $`\sigma (500)`$, we obtain that $`\sigma (500)\pi ^+`$ channel produces the largest decay fraction, as shown in Table 3. The non-resonant amplitude, which is dominant in the model without $`\sigma (500)\pi ^+`$, drops substantially. This model describes the data much better, as can be seen by the $`\pi \pi `$ mass squared projection in Fig. 6(b). The $`\chi ^2/\nu `$ is now 0.9, with a corresponding confidence level of 91%.

## 5 Conclusion

From the data of the Fermilab E791 experiment, we studied the Dalitz plots of the decays $`D^+K^{}\pi ^+\pi ^+`$, $`D_s^+\pi ^{}\pi ^+\pi ^+`$ and $`D^+\pi ^{}\pi ^+\pi ^+`$. In these three final states, the scalar intermediate resonances were found to give the main contribution to the decay rates. We obtained strong evidence for the existence of $`\sigma (500)`$ and $`\kappa `$ scalar mesons, measuring their masses and widths. We also obtained new measurements for masses and widths of the other scalars studied, $`f_0(980)`$, $`f_0(1430)`$ and $`K_0^{}(1430)`$.

The results presented here show the potential of $`D`$ meson decays for the study of light meson spectroscopy, in particular in the scalar sector.


text:
# Solving the Maxwell equations by the Chebyshev method: A one-step finite-difference time-domain algorithm

## I Introduction

Most finite-difference time-domain (FDTD) calculations solve the time-dependent Maxwell equations using algorithms based on a proposal by Yee Taflove ; Kunz ; Yee66 . The Yee algorithm is flexible, fast and easy to implement. A limitation of Yee-based FDTD techniques is that their stability is conditional, meaning that their numerical stability depends on the mesh size used for the spatial discretization and on the time step of the time integration Taflove ; Kunz . In practice, the amount of computational work required to solve the time-dependent Maxwell equations by present FDTD techniques Taflove ; Kunz ; Website ; Zheng1 ; Namiki ; Zheng2 ; Harsh00 ; Kole01 ; Kole02  prohibits applications to a class of important fields such as bioelectromagnetics and VLSI design Taflove ; Gandi ; Houshmand . The basic reason for this is that the time step in the FDTD calculation has to be relatively small in order to maintain a reasonable degree of accuracy in the time integration.

In this paper we describe a one-step algorithm, based on Chebyshev polynomial expansions TAL-EZER0 ; TAL-EZER ; LEFOR ; Iitaka01 ; SILVER ; LOH , to solve the time-dependent Maxwell equations for arbitrarily long times. We demonstrate that the computational efficiency of this one-step algorithm can be orders of magnitude larger than of other FDTD techniques.

## II Algorithm

We consider EM fields in linear, isotropic, nondispersive and lossless materials. The time evolution of EM fields in these systems is governed by the time-dependent Maxwell equations BornWolf . Some important physical symmetries of the Maxwell equations can be made explicit by introducing the fields

$$𝐗(t)\sqrt{\mu }𝐇(t)\mathrm{and}𝐘(t)\sqrt{\epsilon }𝐄(t).$$
(1)
Here, $`𝐇(t)=(H_x(𝐫,t),H_y(𝐫,t),H_z(𝐫,t))^T`$ denotes the magnetic and $`𝐄(t)=(E_x(𝐫,t),E_y(𝐫,t),E_z(𝐫,t))^T`$ the electric field vector, while $`\mu =\mu (𝐫)`$ and $`\epsilon =\epsilon (𝐫)`$ denote, respectively, the permeability and the permittivity. In the absence of electric charges, Maxwell’s curl equations Taflove  read

$$\frac{}{t}\left(\begin{array}{c}𝐗(t)\\ 𝐘(t)\end{array}\right)=\left(\begin{array}{c}𝐗(t)\\ 𝐘(t)\end{array}\right)\frac{1}{\sqrt{\epsilon }}\left(\begin{array}{c}0\\ 𝐉(t)\end{array}\right),$$
(2)
where $`𝐉=(J_x(𝐫,t),J_y(𝐫,t),J_z(𝐫,t))^T`$ represents the source of the electric field and $``$ denotes the operator

$$\left(\begin{array}{cc}0& \frac{1}{\sqrt{\mu }}\times \frac{1}{\sqrt{\epsilon }}\\ \frac{1}{\sqrt{\epsilon }}\times \frac{1}{\sqrt{\mu }}& 0\end{array}\right).$$
(3)
Writing $`𝐙(t)=(𝐗(t),𝐘(t))^T`$ it is easy to show that $``$ is skew symmetric, i.e. $`^T=`$, with respect to the inner product $`𝐙|𝐙^{}_V𝐙^T𝐙^{}𝑑𝐫`$, where $`V`$ denotes the system’s volume. In addition to Eq.(2), the EM fields also satisfy $`(\sqrt{\mu }𝐗(t))=0`$ and $`(\sqrt{\epsilon }𝐘(t))=0`$.

A numerical algorithm that solves the time-dependent Maxwell equations necessarily involves some discretization procedure of the spatial derivatives in Eq. (2). Ideally, this procedure should not change the basic symmetries of the Maxwell equations. We will not discuss the (important) technicalities of the spatial discretization (we refer the reader to Refs. Kole01 ; Kole02 ) as this is not essential to the construction of the one-step algorithm.

On a spatial grid Maxwell’s curl equations (2) can be written in the compact form Kole01 ; Kole02

$$\frac{}{t}𝚿(t)=H𝚿(t)𝚽(t).$$
(4)
The vector $`𝚿(t)`$ is a representation of $`𝐙(t)`$ on the grid. The matrix $`H`$ is the discrete analogue of the operator (3), and the vector $`𝚽(t)`$ contains all the information on the current source $`𝐉`$. The formal solution of Eq. (4) is given by

$$𝚿(t)=U(t)𝚿(0)_0^tU(tu)𝚽(u)𝑑u,$$
(5)
where $`U(t)=e^{tH}`$ denotes the time-evolution matrix. The underlying physical symmetries of the time-dependent Maxwell equations are reflected by the fact that the matrix $`H`$ is real and skew symmetric Kole01 , implying that $`U(t)`$ is orthogonal WILKINSON .

Numerically, the time integration is carried out by using a time-evolution operator $`\overline{U}(t)`$ that is an approximation to $`U(t)=e^{tH}`$. We denote the approximate solution by $`\overline{𝚿}(t)`$. First we use the Chebyshev polynomial expansion to approximate $`U(t)`$ and then show how to treat the source term in Eq. (5). We begin by “normalizing” the matrix $`H`$. The eigenvalues of the skew-symmetric matrix $`H`$ are pure imaginary numbers. In practice $`H`$ is sparse so it is easy to compute $`H_1\mathrm{max}_j_i|H_{i,j}|`$. Then, by construction, the eigenvalues of $`BiH/H_1`$ all lie in the interval $`[1,1]`$ WILKINSON . Expanding the initial value $`𝚿(0)`$ in the (unknown) eigenvectors $`𝐛_j`$ of $`B`$, we find from Eq. (5) with $`𝚽(t)0`$:

$$𝚿(t)=e^{izB}𝚿(0)=\underset{j}{}e^{izb_j}𝐛_j𝐛_j|𝚿(0),$$
(6)
where the $`b_j`$ denote the (unknown) eigenvalues of $`B`$. Although there is no need to know the eigenvalues and eigenvectors of $`B`$ explicitly, the current mathematical justification of the Chebyshev approach requires that $`B`$ is diagonalizable and that its eigenvalues are real. The effect of relaxing these conditions on the applicability of the Chebyshev approach is left for future research. We find the Chebyshev polynomial expansion of $`U(t)`$ by computing the expansion coefficients of each of the functions $`e^{izb_j}`$ that appear in Eq. (6). In particular, as $`1b_j1`$, we can use the expansion ABRAMOWITZ  $`e^{izb_j}=J_0(z)+2_{k=1}^{\mathrm{}}i^kJ_k(z)T_k(b_j),`$ where $`J_k(z)`$ is the Bessel function of integer order $`k`$ to write Eq. (6) as

$$𝚿(t)=\left[J_0(z)I+2\underset{k=1}{\overset{\mathrm{}}{}}J_k(z)\stackrel{~}{T}_k(B)\right]𝚿(0).$$
(7)
Here, $`I`$ is the identity matrix and $`\stackrel{~}{T}_k(B)=i^kT_k(B)`$ is a matrix-valued modified Chebyshev polynomial that is defined by the recursion relations

$$\stackrel{~}{T}_0(B)𝚿(0)=𝚿(0),\stackrel{~}{T}_1(B)𝚿(0)=iB𝚿(0),$$
(8)
and

$$\stackrel{~}{T}_{k+1}(B)𝚿(0)=2iB\stackrel{~}{T}_k(B)𝚿(0)+\stackrel{~}{T}_{k1}(B)𝚿(0),$$
(9)
for $`k1`$. In practice we truncate the sum in Eq. (7), i.e. to obtain the approximation $`\overline{𝚿}(t)`$ we will sum only the contributions with $`kK`$. As $`\stackrel{~}{T}_k(B)_11`$ by construction and $`|J_k(z)||z|^k/2^kk!`$ for $`z`$ real ABRAMOWITZ , the resulting error vanishes exponentially fast for sufficiently large $`K`$. In Fig.1 we show a plot of $`J_n(z=200)`$ as a function of $`n`$ to illustrate this point. From Fig.1 it is clear that the Chebyshev polynomial expansion will only be useful if $`K`$ lies to the right of the right-most extremum of $`J_n(z=200)`$. From numerical analysis it is known that for fixed $`K`$, the Chebyshev polynomial is very nearly the same polynomial as the minimax polynomial NumercalRecipes , i.e. the polynomial of degree $`K`$ that has the smallest maximum deviation from the true function, and is much more accurate than for instance a Taylor expansion of the same degree $`K`$. The coefficients $`J_k(z)`$ should be calculated to high precision and the number $`K`$ is fixed by requiring that $`|J_k(z)|<\kappa `$ for all $`k>K`$. Here, $`\kappa `$ is a control parameter that determines the accuracy of the approximation. For fixed $`\kappa `$, $`K`$ increases linearly with $`z=tH_1`$ (there is no requirement on $`t`$ being small), a result that is essential for the efficiency of the algorithm. Using the recursion relation of the Bessel functions, all $`K`$ coefficients can be obtained with $`𝒪(𝒦)`$ arithmetic operations NumercalRecipes . Clearly this is a neglible contribution to the total computational cost for solving the Maxwell equations.

Performing one time step amounts to repeatedly using recursion (9) to obtain $`\stackrel{~}{T}_k(B)𝚿(0)`$ for $`k=2,\mathrm{},K`$, multiply the elements of this vector by $`J_k(z)`$ and add all contributions. This procedure requires storage for two vectors of the same length as $`𝚿(0)`$ and some code to multiply such a vector by the sparse matrix $`H`$. The result of performing one time step yields the solution at time $`t`$, hence the name one-step algorithm. In contrast to what Eqs. (8) and (9) might suggest, the algorithm does not require the use of complex arithmetic.

We now turn to the treatment of the current source $`𝐉(t)`$. The contribution of the source term to the EM field at time $`t`$ is given by the last term in Eq. (5). One approach might be to use the Chebyshev expansion (7) for $`U(tu)=e^{(tu)H}`$ and to perform the integral in Eq. (5) numerically. However that is not efficient as for each value of $`tu`$ we would have to perform a recursion of the kind Eq. (9). Thus, it is better to adopt another strategy. For simplicity we only consider the case of a sinusoidal source

$$𝐉(𝐫,t)=\mathrm{\Theta }(Tt)𝐬(𝐫)\mathrm{sin}(\mathrm{\Omega }t),$$
(10)
where $`𝐬(𝐫)`$ specifies the spatial distribution and $`\mathrm{\Omega }`$ the angular frequency of the source. The step function $`\mathrm{\Theta }(Tt)`$ indicates that the source is turned on at $`t=0`$ and is switched off at $`t=T`$. Note that Eq. (10) may be used to compose sources with a more complicated time dependence by a Fourier sine transformation.

The formal solution for the contribution of the sinusoidal source (10) reads

$`{\displaystyle _0^t}e^{(tu)H}𝚽(u)𝑑u`$ $`=`$ $`(\mathrm{\Omega }^2+H^2)^1e^{(tT^{})H}\times (\mathrm{\Omega }e^{T^{}H}\mathrm{\Omega }\mathrm{cos}\mathrm{\Omega }T^{}H\mathrm{sin}\mathrm{\Omega }T^{})𝚵`$ (11)
$``$ $`f(H,t,T^{},\mathrm{\Omega })𝚵,`$
where $`T^{}=\mathrm{min}(t,T)`$ and $`𝚽(u)\mathrm{\Theta }(Tt)\mathrm{sin}(\mathrm{\Omega }t)𝚵`$ with $`𝚵`$ a vector of the same length as $`𝚿(0)`$ that represents the time-independent, spatial distribution $`𝐬(𝐫)`$. The coefficients of the Chebyshev polynomial expansion of the formal solution (11) are calculated as follows. First we repeat the scaling procedure described above and substitute in Eq. (11) $`H=ixH_1`$, $`t=z/H_1`$, $`T^{}=Z^{}/H_1`$, and $`\mathrm{\Omega }=\omega H_1`$. Then, we compute the (Fast) Fourier Transform with respect to $`x`$ of the function $`f(x,z,Z^{},\omega )`$ (which is non-singular on the interval $`1x1`$). By construction, the Fourier coefficients $`S_k(tH_1)`$ are the coefficients of the Chebyshev polynomial expansion ABRAMOWITZ .

Taking into account all contributions of the source term with $`k`$ smaller than $`K^{}`$ (determined by a procedure similar to the one for $`K`$), the one-step algorithm to compute the EM fields at time $`t`$ reads

$`\overline{𝚿}(t)`$ $`=`$ $`\left[J_0(tH_1)I+2{\displaystyle \underset{k=1}{\overset{K}{}}}J_k(tH_1)\stackrel{~}{T}_k(B)\right]𝚿(0)`$ (12)
$`+`$ $`\left[S_0(tH_1)I+2{\displaystyle \underset{k=1}{\overset{K^{}}{}}}S_k(tH_1)\stackrel{~}{T}_k(B)\right]𝚵.`$
We emphasize that in our one-step approach the time dependence of the source is taken into account exactly, without actually sampling it.

## III Results

The following two examples illustrate the efficiency of the one-step algorithm. First we consider a system in vacuum ($`\epsilon =\epsilon _0`$ and $`\mu =\mu _0`$) which is infinitely large in the $`y`$\- and $`z`$-direction, hence effectively one dimensional. The current source (10) is placed at the center of a system of length 250.1 and oscillates with angular frequency $`\mathrm{\Omega }=2\pi `$ during the time interval $`0tT=4`$ units . In Table 1 we present results of numerical experiments with two different time-integration algorithms. In general, the error of a solution $`\stackrel{~}{𝚿}(t)`$ as obtained by the FDTD algorithm of Yee Taflove ; Yee66  or the unconditionally stable FDTD algorithm $`T4S2`$ Kole01 ; Kole02  is defined by $`\mathrm{\Delta }(t)\stackrel{~}{𝚿}(t)\overline{𝚿}(t)/\overline{𝚿}(t)`$, where $`\overline{𝚿}(t)`$ denotes the vector of the EM fields as obtained by the one-step algorithm. The error on the Yee-algorithm result vanishes as $`\tau ^2`$ for sufficiently small $`\tau `$ Taflove ; Yee66 . However, as Table I shows, unless $`\tau `$ is made sufficiently small ($`\tau 0.0125`$ in this example), the presence of the source term changes the quadratic behavior to almost linear. The rigorous bound on the error between the exact and $`T4S2`$ results tells us that this error should vanish as $`\tau ^4`$ Kole01 ; DeRaedt87 . This knowledge can be exploited to test if the one-step algorithm yields the exact numerical answer. Using the triangle inequality we can write

$`𝚿(t)\overline{𝚿}(t)𝚿(t)\stackrel{~}{𝚿}(t)`$ $`+`$ $`\stackrel{~}{𝚿}(t)\overline{𝚿}(t)`$
$`\tau ^4tC\left(1+{\displaystyle _0^t}𝐉(u)𝑑u\right)`$ $`+`$ $`\mathrm{\Delta }(t)\overline{𝚿}(t)`$ (13)
where $`C`$ is a positive constant DeRaedt87 . The numerical data in Table 1 (third column) show that $`\mathrm{\Delta }(t)0`$ as $`\tau ^4`$ and, therefore, we can be confident that the one-step algorithm yields the correct answer within rounding errors. Furthermore, since the results of the one-step algorithm are exact within almost machine precision, in general the solution also satisfies $`(\sqrt{\mu }𝐗(t))=0`$ and $`(\sqrt{\epsilon }𝐘(t))=0`$ within the same precision. This high precision also allows us to use the one-step algorithm for genuine time stepping with arbitrarily large time steps, this in spite of the fact that strictly speaking, the one-step algorithm is not unconditionally stable.

From Table 1 it follows that if one finds an error of more than 2.5% acceptable, one could use the Yee algorithm, though we recommend to use the one-step algorithm because then the time-integration error is neglegible. The Yee algorithm is no competition for the $`T4S2`$ algorithm if one requires an error of less than 1%, but the $`T4S2`$ algorithm is not nearly as efficient as the one-step algorithm with respect to the number of required matrix-vector operations.

A more general quantitative analysis of the efficiency can be made using the fact that for an $`n`$th-order algorithm ($`n=2`$ for the Yee algorithm and $`n=4`$ for the $`T4S2`$ algorithm), the error $`\mathrm{\Delta }(t)`$ vanishes no faster with $`\tau `$ than $`\tau ^nt`$. Each time step takes a number $`W(n)`$ of matrix-vector operations (of the type $`𝚿^{}M𝚿`$), e.g. for a three-dimensional system we have $`W(2)=1`$ and $`W(4)=10`$ for the Yee algorithm and the $`T4S2`$ algorithm, respectively. In practice the actual number of floating point operations carried out by our algorithms agrees with these estimates. The total number of matrix-vector operations it takes to obtain the solution at a reference time $`t_r`$ with error $`\mathrm{\Delta }_r(t_r)`$ is then given by $`N_r=W(n)t_r/\tau _r`$ and thus $`\mathrm{\Delta }_r(t_r)W(n)^nt_r^{n+1}/N_r^n`$. The number of operations $`N`$ that it will take to compute the EM fields at time $`t`$ with accuracy $`\mathrm{\Delta }(t)`$ is then calculated from

$$N=N_r\left(\frac{\mathrm{\Delta }_r(t_r)}{\mathrm{\Delta }(t)}\right)^{1/n}\left(\frac{t}{t_r}\right)^{(n+1)/n}.$$
(14)
We note that one numerical reference experiment per $`n`$th-order algorithm is sufficient to determine the parameters $`N_r`$, $`\mathrm{\Delta }_r(t_r)`$, and $`t_r`$. While these parameters may be different for different systems, the scaling of $`N`$ with $`t^{3/2}`$ and with $`t^{5/4}`$, respectively, for second- and fourth-order algorithms, will not be affected. Most importantly, since the number of matrix-vector operations required by the one-step algorithm scales linearly with $`t`$, it is clear that for long enough times $`t`$, the one-step algorithm will be orders of magnitude more efficient than the current FDTD methods. In Fig.2 we show the required number of operations as a function of time $`t`$ taking, as an example, simulation data of 3D systems (discussed below) to fix the parameters $`N_r`$, $`\mathrm{\Delta }_r(t_r)`$, and $`t_r`$. We conclude that for longer times none of the FDTD algorithms can compete with the one-step algorithm in terms of efficiency. For $`t=20`$, the one-step algorithm is a factor of ten faster than the Yee algorithm. Thereby we have disregarded the fact that the Yee algorithm yields results within an error of 0.1% while the one-step algorithm gives the numerically exact solution.

As the second example we use the one-step algorithm to compute the frequency spectrum of a three-dimensional photonic woodpile Lin98 . This structure, shown in the inset of Fig. 3, possesses a large infrared bandgap and is under current experimental and theoretical investigation Lin98 ; Fleming02 . To determine all eigenvalues of the corresponding matrix $`H`$ we follow the procedure described in Refs. Kole01 ; Alben75 ; Hams00 . We use random numbers to initialize the elements of the vector $`𝚿(0)`$. Then we calculate the inner product $`F(t)=𝚿(0)|𝚿(t)`$ as a function of $`t`$ and average $`f(t)=F(t)/F(0)`$ over several realizations of the initial vector $`𝚿(0)`$. The full eigenmode distribution, $`𝒟(\omega )`$, is obtained by Fourier transformation of $`f(t)`$. In Fig. 3 we show $`𝒟(\omega )`$, as obtained by $`T4S2`$ and the one-step algorithm, with a time step $`\tau =0.075`$ (set by the largest eigenvalue of $`H`$), a mesh size $`\delta =0.1`$, and 8192 time steps. For this choice of parameters, the Yee algorithm would be unstable Taflove ; Kunz  and would yield meaningless results. The $`T4S2`$ calculation shows a peak at $`\omega =0`$. This reflects the fact that, in a strict sense, the $`T4S2`$ algorithm does not conserve $`(\sqrt{\mu }𝐗(t))`$ and $`(\sqrt{\epsilon }𝐘(t))`$ Kole01 ; Kole02 . However, the peak at $`\omega =0`$ vanishes as $`\tau ^4`$. Repeating the $`T4S2`$ calculation with $`\tau =0.01`$ yields a $`𝒟(\omega )`$ (not shown) that is on top of the result of the one-step algorithm (see Fig. 3) and is in good agreement with band-structure calculations Lin98 . For $`\tau =0.01`$ the one-step algorithm is 3.5 times more efficient than $`T4S2`$. Note that in this example, the one-step algorithm is used for a purpose for which it is least efficient (time-stepping with relatively small time steps). Nevertheless the gain in efficiency is still substantial. In simulations of the scattering of the EM fields from the same woodpile (results not shown), the one-step algorithm is one to two orders of magnitude more efficient than current FDTD algorithms, in full agreement with the error scaling analysis given above.

## IV Conclusion

We have described a one-step algorithm, based on the Chebyshev polynomial expansions, to solve the time-dependent Maxwell equations with spatially varying permittivity and permeability and current sources. In practice this algorithm is as easy to implement as FDTD algorithms. Our error scaling analysis shows and our numerical experiments confirm that for long times the one-step algorithm can be orders of magnitude more efficient than current FDTD algorithms. This opens possibilities to solve problems in computational electrodynamics that are currently intractable.

###### Acknowledgements.

H.D.R. and K.M. are grateful to T. Iitaka for drawing our attention to the potential of the Chebyshev method and for illuminating discussions.


text:
# On the Stability of Particles in Extrasolar Planetary Systems

## 1. Introduction

The discovery of multi-body extrasolar planetary systems during the past few years has once again confronted astrodynamicists with the old question of the stability of such systems. The discovery of GJ 876 where two planets are locked in a 2:1 mean-motion resonance (Marcy et al. 2001), the confirmation of three planets in orbit about Upsilon Andromedae (Butler et al. 1999), and the discovery of planetary systems around 47 Uma and 55 Cancri with planets in orbits more closely resembling those in the Solar System (Fischer et al. 2002; Marcy et al. 2002), have set the grounds for a deeper look at the problem of the stability of multi-body systems. Here we present a preliminary investigation of the dynamical stability of the 47 Uma and 55 Cancri systems. A more detailed analysis, including similar work on GJ 876 and $`\upsilon `$ And, will be addressed in future work.

## 2. Methodology

We performed dynamical fits to the radial velocity data for the stars 47 Uma and 55 Cancri using a Levenberg-Marquardt minimization algorithm (Press et al. 1992), as in Laughlin & Chambers (2001) and Rivera & Lissauer (2001). We assumed that the planets in each system were coplanar and that the plane of the planets contained the line of sight. We then used the resulting fitted parameters as initial conditions for N-body simulations. In order to give the resulting systems the ability to explore extra degrees of freedom, we artificially added a mutual inclination of one degree. We also performed simulations of each system in which hundreds of (massless) test particles were added. All test particles were started on circular orbits with respect to the central star. No test particles were placed in regions such that a test particle would initially cross the orbit of a planet. The results may be used as an indication of the presence and locations of potential terrestrial planets in these systems.

The simulations were performed with the second-order mixed variable symplectic (MVS) integrator in the MERCURY integration package (Chambers 1999), which is based on the technique pioneered by Wisdom & Holman (1991). This code was modified to include the principal effects of general relativity, as in Lissauer & Rivera (2001).

## 3. 47 Uma

In agreement with Fischer et al. (2002), we find that in dynamical fits to the radial velocity data for 47 Uma, the value of $`\chi _\nu ^2`$ does not change significantly for a fit with the eccentricity of the outer companion in the range 0 to 0.2. When we fit for the eccentricity of the outer planet, the fitting routine converged on a value of 0.19. Note that this value is near the boundary of stability determined by Fischer et al. (2002). A long-term simulation of the system based on this fit is stable over billions of years. In this simulation, the two companions are in a 5:2 mean motion resonance. Also, throughout the simulation the orbits are nearly anti-aligned. Table 1 gives the parameters from the fit at epoch JD 2446959.737, and Figure 1a shows the nearest and farthest points in the companions’ orbits vs. time for the first 1 Gyr of a long-term simulation based on the parameters given in Table 1. Figure 2 shows the eccentricities, inclinations, and the period ratio of the companions over short times.

Figure 1b shows the stability of test particles in a 10 Myr simulation of the same system. It shows the time that a test particle was lost vs. its initial semimajor axis. The dots and lines toward the bottom of the figure indicate the initial semimajor axes and radial excursions of the planets. Note that islands of instability occur at the locations of several mean motion resonances with the planets. Also, exterior to the outer companion, there are islands of stability between mean motion resonances. Since 47 Uma is similar to the sun, the region around 1 AU is of particular interest. Unfortunately, the 3:1 mean motion resonance with the inner companion lies in this region, and there are signs that the region may not be stable. Note that while only the one test particle at 1 AU was lost in less than 10 Myr, proximity to the 3:1 mean motion resonance could endanger the stability of other nearby orbits over longer timescales.

## 4. 55 Cancri

This system proved to be difficult to fit with the Levenberg-Marquardt routine. We obtained two fits with almost the same $`\chi _\nu ^2`$ value with very different parameters for the outer planet. This is an indication of the large uncertainties in this planet’s orbital parameters. Table 1 gives the parameters from one of our fits at epoch JD 2450250.0, and Figure 3a shows the periastron and apastron of each planet vs. time from a simulation based on this fit. The simulation shows that the inner and middle companions are near a 3:1 mean motion resonance. These results are in agreement with Marcy at al. (2002).

Figure 3b shows the stability of test particles in a 5 Myr simulation of the same system. It clearly shows a large stable region between the middle and outer companions. This is in rough agreement with Marcy et al. (2002), in which they showed that a terrestrial planet on a circular orbit at 1 AU would be stable.

## 5. Summary

We have presented results on the stability of the 47 Uma and 55 Cancri systems with and without massless particles. Dynamical fits were used to determine the initial conditions for the simulations. A large area of parameter space is consistent with the observations. The area around 1 AU in the 47 Uma system may not harbor a terrestrial planet, while the opposite is true in the 55 Cancri system. These results are in rough agreement with previous studies (Fischer et al. 2002; Marcy et al. 2002).

### Acknowledgments.

We thank Debra Fischer, Geoff Marcy, and Paul Butler for providing their radial velocity observations.

## References

Butler, R. P., Marcy, G. W., Fischer, D. A., Brown, T. M., Contos, A. R., Korzennik, S. G., Nisenson, P., & Noyes, R. W. 1999, ApJ, 526, 916

Chambers, J. E. 1999, MNRAS, 304, 793

Fischer, D. A., Marcy, G. W., Butler, R. P., Laughlin, G., & Vogt, S. S. 2002, ApJ, 564, 1028

Laughlin, G., & Chambers, J. E. 2001, ApJ, 551, L109

Lissauer, J. J. & Rivera, E. J. 2001, ApJ, 554, 1141

Marcy, G. W., Butler, R. P., Fischer, D. A., Vogt, S. S., Lissauer, J. J., & Rivera, E. J. 2001, ApJ, 556, 296

Marcy, G. W., Butler, R. P., Fischer, D. A., Laughlin, G., Vogt, S. S., Henry, G. W., & Pourbaix, D. 2002, ApJ, submitted

Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. 1992, Numerical Recipes in Fortran (2nd ed.; Cambridge: Cambridge Univ. Press)

Rivera, E. J. & Lissauer, J. J. 2001, ApJ, 558, 392

Wisdom, J. & Holman, M. 1991, AJ, 102, 1528


text:
# Top-charm associated production at hadron colliders in the standard model with large extra dimensions 11footnote 1Supported by National Natural Science Foundation of China.

The precise calculations are carried out on the flavor changing neutral current couplings in the process $`ppggt\overline{c}\left(\overline{t}c\right)`$ at the large hadron collider(LHC) and very large hadron collider(VLHC) in both frameworks of the minimal standard model(MSM) and its extension with extra dimensions. We find that the effects from the large extra dimensions can enhance the total cross section up to about several hundred times as that in the MSM, quantitatively.

PACS: 04.50.+h, 12.60.-i, 13.85.-t

Keywords: Extra dimensions, Top-charm associated production

It is known that the minimal standard model(MSM) has achieved a great success under the electroweak energy scale, $`m_{EW}`$, which is about $`10^2GeV`$ energy scale . But there still exist several problems, such as hierarchy problem, which can not been solved in this theory. There are many extensions of MSM being proposed to solve these problems, such as supersymmetric models. Recently, a new framework for solving the hierarchy problem, which does not rely on either supersymmetry or technicolor, was established by Arkani-Hamed, Dimopoulos and Dvali(ADD) , They predicted that if we confine the matter fields to our 4-dimensional world, the size of extra dimensions could be large enough to be detectable. In this model, the gravitation and gauge interactions are united at the weak scale, $`M_D`$, which is the only fundamental scale in this theory, although it is now a model-dependent parameter, and no new low-energy predictions can be derived. This has the exciting implication that future high-energy collider experiments can directly probe the physics of quantum gravity. Above $`TeV`$ energy scale, completely new phenomena could emerge as resonant production of the Regge recurrences of string theory or excitations of Kaluza-Klein modes of ordinary particles .

There are stringent experimental constraints against the existence of tree-level flavor changing scalar interactions (FCSI’s) involving the light quarks. This leads to the suppression of the flavor changing neutral current(FCNC) couplings, an important feature of MSM, which is explained in terms of the Glashow-Illiopoulos-Maiani(GIM) . The advantage of examining top quark physics than other quark physics is that one can directly determine the properties of top quark itself and does not need to worry about no-perturbative QCD effects which are difficult to attack because there exist no top-flavored hadron states at all. The properties of top quark makes the study on the FCNC couplings to be one of important fields in top physics. The FCNC couplings correlated with top quarks can be probed either in rare decays of $`t`$-quark or via top-charm associated production which appear at loop-level and offer a good place to test quantum effects of the fundamental quantum field theory. In models beyond MSM, new particles may appear in the loop and have significant contributions to flavor changing transitions. Any positive experimental observation of FCNC existence deviated from that in the MSM would unambiguously pronounce the presence of new physics. On the contrary, in case no deviation from the SM is observed, they define in a quantitative way the strategy to obtain lower bounds on the new physics energy scale.

In this letter we shall investigate the process $`ppt\overline{c}(\overline{t}c)`$ at hadron colliders in the MSM extension with large extra dimensions. The cross section of the subprocess $`ppd\overline{d}t\overline{c}(\overline{t}c)`$ is proportional to $`|V_{dt}^{}V_{dc}|^2`$, but that of the subprocess $`ppggt\overline{c}(\overline{t}c)`$ is approximately proportional to $`|V_{bt}^{}V_{bc}|^2`$ for $`|V_{bt}^{}V_{bc}|^2|V_{st}^{}V_{sc}|^2`$. According to the latest data from Ref. , $`\frac{|V_{bt}^{}V_{bc}|^2}{|V_{dt}^{}V_{dc}|^2}2.5\times 10^3`$. The luminosity of gluon from proton at hadron colliders is much larger than that of $`d`$-quark. Although the process $`ggt\overline{c}`$ appears at one-loop level, while the process $`d\overline{d}t\overline{c}`$ does at tree level, we can still conclude that the cross section of the first one will be much larger than that of the later process. We found exactly that the cross section of $`ppggt\overline{c}`$ is almost 100 times large as that of $`ppd\overline{d}t\overline{c}`$ after a precise numerical computing. So we will focus only on the process $`ppggt\overline{c}`$ in the following investigation.

There were some discussions about $`t\overline{c}(\overline{t}c)`$ production at the LHC based on the THDM and the R-violating MSSM. The results of Ref.  show that the signals of the $`t\overline{c}(\overline{t}c)`$ production induced by those models maybe observable, but strongly depends on the chosen parameters. The research on $`t\overline{c}(\overline{t}c)`$ production at $`e^+e^{}`$ in the SM was presented in Ref., in the THDM-III in Ref. and in the MSSM in Ref.. In this work we consider the process $`ppt\overline{c}(\overline{t}c)`$ at hadron colliders in the MSM extension with large extra dimensions. Within the framework of the MSM with extra dimensions, the Feynman diagrams with graviton appear only in the $`s`$-channel and the graviton presents as a propagator in the $`t\overline{c}(\overline{t}c)`$ association production, but not appears in loop. Then in this process the FCNC couplings originates still in CKM mixing matrix rather than other new mechanics. Since the scale of reaction energy at the LHC and VLHC has gone up to $`14TeV`$ and $`100TeV`$, respectively, the effects of the new dimensions will be expected emerge observably. There are eighteen diagrams including the three diagrams involving graviton for the subprocess $`ggt\overline{c}(\overline{t}c)`$ as shown in Fig.1.

Here we present the relevant Feynman rules for our calculation which can be extracted from Ref.. The relevant Feynman rules are listed below,

$$forFig.2\left(a\right),\frac{1}{2}\frac{i}{k^2m^2}\left(\eta _{\mu \mu ^{}}\eta _{\nu \nu ^{}}+\eta _{\mu \nu ^{}}\eta _{\nu \mu ^{}}\frac{2}{D2}\eta _{\mu \nu }\eta _{\mu ^{}\nu ^{}}\right),$$

$$forFig.2\left(b\right),i\frac{\kappa }{\sqrt{V}_\delta }\left[B_{\mu \nu \alpha \beta }m_A^2+\left(C_{\mu \nu \alpha \beta \rho \sigma }C_{\mu \nu \alpha \sigma \beta \rho }\right)k_1^\rho k_2^\sigma \right]\delta ^{ab},$$

$$forFig.2\left(c\right),i\frac{\kappa }{8\sqrt{V}_\delta }\left[\gamma _\mu \left(p_1+p_2\right)_\nu +\gamma _\nu \left(p_1+p_2\right)_\mu 2\eta _{\mu \nu }\left(\text{/}p_1+\text{/}p_22m_\psi \right)\right],$$

respectively, where the flat matric $`\eta _{\mu \nu }=diag\{+1,1,1,1\}`$. $`\delta =D4`$ is the number of the extra dimensions. $`V_\delta `$ is the volume of the compactified space. $`B_{\mu \nu \alpha \beta }`$ and $`C_{\mu \nu \alpha \beta \rho \sigma }`$ are defined as below

$`B_{\mu \nu \alpha \beta }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\eta _{\mu \nu }\eta _{\alpha \beta }\eta _{\mu \alpha }\eta _{\nu \beta }\eta _{\mu \beta }\eta _{\nu \alpha }\right),`$
$`C_{\mu \nu \alpha \beta \rho \sigma }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\eta _{\mu \nu }\eta _{\alpha \beta }\eta _{\rho \sigma }\left(\eta _{\mu \rho }\eta _{\nu \sigma }\eta _{\alpha \beta }+\eta _{\mu \sigma }\eta _{\nu \rho }\eta _{\alpha \beta }+\eta _{\mu \alpha }\eta _{\nu \beta }\eta _{\rho \sigma }+\eta _{\mu \beta }\eta _{\nu \alpha }\eta _{\rho \sigma }\right)\right).`$
For the s-channel virtual graviton exchange diagram, the sum over all Kaluza-Klein modes has to be performed at the amplitude level (The detailed deduction can be found in Refs. ). In order to absorbing the factor $`\frac{\kappa }{\sqrt{V}_\delta }`$ of the interaction vertex for which we have the equation $`\frac{\kappa }{\sqrt{V}_\delta }=\frac{2}{\overline{M}_P}`$, where $`\overline{M}_P=M_P/\sqrt{8\pi }=2.4\times 10^{18}GeV`$ is the reduced Plank mass. The we can define

$$D\left(s\right)=\frac{1}{\overline{M_P}^2}\underset{n}{}\frac{1}{sm_n^2}.$$

Converting the sum into an integral which can be evaluated using dimensional regularization, we have

$$D\left(s\right)=\frac{1}{M_D^{2+\delta }}\frac{S_{\delta 1}}{2}C\mathrm{\Lambda }^{\delta 2}$$

for $`\delta >2`$. Where $`M_D10^2GeV`$ is the only fundamental mass scale in the model, $`S_{\delta 1}`$ is the hypersurface of a unit-radius sphere in $`\delta `$ dimensions, $`\mathrm{\Lambda }`$ is the ultraviolet cutoff, $`C`$ is an unknown coefficient. So $`D(s)`$ is fixed by several uncertain parameter. Here we can postulate $`C1`$, $`\mathrm{\Lambda }M_D`$ and have

$$D\left(s\right)=\frac{1}{M_D^4}\frac{S_{\delta 1}}{2}.$$

In our calculation we take $`M_D`$ in the range of $`100500GeV`$.

In the numerical calculations the following input values of the parameters have been employed

$`m_t=174.3GeV,m_c=1.25GeV,m_b=4.2GeV,m_s=0.12GeV`$
$`m_W=80.42GeV,\alpha ={\displaystyle \frac{1}{128}},\alpha _s=0.117.`$
We take the ”standard” parameterization of CKM matrix  and chose

$$s_{12}=0.221814,s_{23}=0.0459998,s_{13}=0.00500002,\delta _{13}=0.$$

The total cross section for parent process at $`pp`$ collider can be obtained by folding the cross section of subprocess $`\widehat{\sigma }(ggt\overline{c})`$ with the gluon luminosity.

$$\sigma \left(s,ppggt\overline{c}+X\right)=_{\left(m_t+m_c\right)^2/s}^1𝑑\tau \frac{dL_{gg}}{d\tau }\widehat{\sigma }\left(ggt\overline{c}at\widehat{s}=\tau s\right),$$

where $`\sqrt{s}`$ and $`\sqrt{\widehat{s}}`$ are the $`pp`$ and $`gg`$ c.m.s. energies respectively, at the LHC $`\sqrt{s}=14TeV`$ and the VLHC $`\sqrt{s}=100TeV`$, and $`dL_{gg}/d\tau `$ is the distribution function of gluon luminosity, which is defined as

$$\frac{dL_{gg}}{d\tau }=_\tau ^1\frac{dx_1}{x_1}\left[f_g(x_1,Q^2)f_g(\frac{\tau }{x_1},Q^2)\right],$$

where $`\tau =x_1x_2`$, the definitions of $`x_1`$ and $`x_2`$ are from Ref., and in our calculation we adopt the CTEQ5L parton distribution functions . The factorization scale Q is chosen as the average of the final particles masses $`\frac{1}{2}(m_t+m_c)`$. The total cross section of $`ppt\overline{c}(\overline{t}c)+X`$ is obtained by multiplying the cross section of $`ppt\overline{c}+X`$ by factor 2. We have the total cross section in the MSM $`4.23\times 10^3fb`$ at the LHC and $`0.105fb`$ at the VLHC.

The cross sections of parent process $`ppggt\overline{c}+\overline{t}c+X`$ at the LHC as the functions of $`M_D`$ are depicted in Fig.3 with $`\delta =3,5,7`$, respectively. The figure shows that when $`M_D=100GeV`$, the cross sections are about $`1fb`$ and sharply decrease with the increament of $`M_D`$ at first, then approach asymptotically to the cross section in the MSM after the point of $`M_D=400GeV`$. In Fig.5 the curves describe the relations between the cross sections of parent process $`ppggt\overline{c}+\overline{t}c+X`$ and $`M_D`$ at the VLHC, which is almost the same feature as that at the LHC but the cross sections at the VLHC are generally as several hundred times as those at the LHC. So if we take $`M_D=100GeV`$ and assume that the LHC and the VLHC run with integrated luminosity of $`L300[fb]^1`$ in one year, there will be the order of $`10^210^3`$ events per year accumulated at the LHC and $`10^410^5`$ at the VLHC, respectively. In Fig.4 and Fig.6, we present the cross section relative enhancement $`(\sigma \sigma _{SM})/\sigma _{SM}`$ by the virtual graviton exchange as the functions of $`M_D`$ at the LHC and the VLHC, respectively. Since $`M_D`$ appears in denominator of graviton propagator, the cross section decreases with the increment of $`M_D`$ as shown in the figures. We can figure out that the correction from the large extra dimensions to the SM cross section is approximately proportional to $`M_D^4`$. The contribution of $`\delta `$ origins mainly in the hypersurface of a unit-radius sphere in $`\delta `$ dimensions. As we know, $`S_{\delta 1}=\pi ^{\frac{\delta }{2}}/\mathrm{\Gamma }(\frac{\delta }{2})`$. We have $`S_{\delta 1}|_{\delta =3}<S_{\delta 1}|_{\delta =5}<S_{\delta 1}|_{\delta =7}`$, so the cross sections of $`\delta =7`$ are larger than those of $`\delta =5(\delta =3)`$ in the same condition as shown in all figures.

To summarize, we find that the the cross section of top-charm associated production at hadron colliders is largely enhanced due to large extra dimensions within the framework of the MSM with large extra dimensions, and can be detectable both at the LHC and the VLHC with the favorable parameters. This cross section enhancement for the $`t\overline{c}(\overline{t}c)`$ association production process is strongly related to the energy scale $`M_D`$. When $`M_D`$ is above $`400GeV`$, the difference between the cross sections in the MSM with and without large extra dimensions is indistinctly demonstrated. The measurement of the cross section of $`ppggt\overline{c}+\overline{t}c+X`$ process can be used to give the constraint on the fundamental energy scale $`M_D`$.

Acknowledgements:

This work was supported in part by the National Natural Science Foundation of China and a grant from University of Science and Technology of China.

Figure Captions

Fig.1 The Feynman diagrams of the subproecesses $`ggt\overline{c}+\overline{t}c`$.

Fig.2 (a) The propagator of graviton. (b) The coupling between the graviton and gluons. (c) The coupling between the graviton and fermions.

Fig.3 The cross sections of $`ppggt\overline{c}+\overline{t}c+X`$ at the LHC as the functions of the energy scale $`M_D`$.

Fig.4 The enhancement to the cross sections of $`ppggt\overline{c}+\overline{t}c+X`$ at the VLHC as the functions of the energy scale $`M_D`$.

Fig.5 The cross sections of $`ppggt\overline{c}+\overline{t}c+X`$ at the LHC as the functions of the energy scale $`M_D`$.

Fig.6 The enhancement to the cross sections of $`ppggt\overline{c}+\overline{t}c+X`$ at the VLHC as the functions of the energy scale $`M_D`$.


text:
# Local excitations of a spin glass in a magnetic field

## Abstract

We study the minimum energy clusters (mec) above the ground state for the $`3d`$ Edwards-Anderson Ising spin glass in a magnetic field. For fields $`B`$ below $`0.4`$, we find that the field has almost no effect on the excitations that we can probe, of volume $`V64`$. As found previously for $`B=0`$, their energies decrease with $`V`$, and their magnetization remains very small (even slightly negative). For larger fields, both the mec energy and magnetization grow with $`V`$, as expected in a paramagnetic phase. However, all results appear to scale as $`BV`$ (instead of $`B\sqrt{V}`$ as expected from droplet arguments), suggesting that the spin glass phase is destroyed by any small field. Finally, the geometry of the mec is completely insensitive to the field, giving further credence that they are lattice animals, in the presence or the absence of a field.

*Introduction —* Strong disorder or frustration in a system can lead to a low temperature phase with “frozen order”. This is believed to be the case for spin glasses, the archetypes of such systems Mézard et al. (1987); Young (1998). Does such a frozen phase persist under generic perturbations? Here we consider that question in the context of the Edwards-Anderson Edwards and Anderson (1975) Ising spin glass when the perturbation is an external magnetic field. The effect of the field is to align the spins in one direction and this bias may break up the spin glass ordering. Do arbitrarily small fields destroy the frozen order, or can the ordering co-exist with a field as long as it is not too large? This is a long-standing question that has been surprisingly difficult to answer.

From an experimental point of view, this question has been addressed only twice for systems of Ising spins. In the first study Ito and Katori (1994), the (non-equilibrium) properties were interpreted to suggest the presence of spin glass ordering at low fields. Several years later the same sample was re-analyzed and the lines of constant relaxation time determined Mattsson et al. (1995). From that, a scaling analysis was performed, leading to the conclusion that the relaxation time is finite for non-zero fields, and thus the system is paramagnetic. On the theoretical size, the situation remains very controversial. In the mean field picture, where one is guided by the Sherrington-Kirkpatrick Sherrington and Kirkpatrick (1975) or other mean field models, the de Almeida-Thouless line de Almeida and Thouless (1978) $`B_{AT}(T)`$ separates the paramagnetic and spin glass phases, and one has $`B_{AT}>0`$ if $`T<T_c`$ ($`T_c`$ is the critical temperature in zero field). Thus spin glass ordering continues in the presence of a sufficiently small magnetic field in the mean field picture. At present, there is no consensus whether this picture correctly describes what happens in the $`d=3`$ Ising spin glass, nor even when $`d`$ is large but finite. Indeed it has been argued by Fisher and Huse Fisher and Huse (1986) that the magnetic field is always relevant, driving the system at large scales towards paramagnetism Pimentel et al. (2002). Then any non-zero magnetic field will destroy the spin glass ordering and will render the system paramagnetic beyond a length scale $`\mathrm{}_B`$. In the droplet Fisher and Huse (1986) or scaling Bray and Moore (1986) pictures, $`\mathrm{}_B`$ diverges as $`B`$ vanishes like

$$\mathrm{}_BB^{2/(d2\theta )}$$
(1)
where $`\theta `$ is the spin stiffness exponent. Because of this divergence as $`B`$ decreases, the relaxation times also diverge. The difficulty is thus to determine whether these relaxation times diverge only as $`B0`$ or at some positive value $`B_c>0`$. Such a question is nearly impossible to address from simulations because the time scales that can be reached always remain quite small, while equilibrium studies are confined to small systems for the same reason. In view of such hurdles, much numerical effort Houdayer and Martin (1999a); Krzakala et al. (2001) in the last few years has focused instead on the zero temperature limit.

In this paper we also use zero temperature, but we study how minimum energy clusters Lamarcq et al. (2002) above the ground state are affected by a magnetic field. Essentially, we search for a putative critical field $`B_c`$; if there is an AT-line, $`B_c`$ is then simply $`B_{AT}(T=0)`$. Our investigation reduces to finding whether there is a phase transition in $`B`$ at $`T=0`$, and for that we need an order parameter. Since the field breaks the up-down symmetry, the Edwards-Anderson order parameter that gives the spatial average of the square magnetization at equilibrium is not of use. It is thus necessary to characterize the spin glass ordering otherwise, for instance by an infinite (spin glass) susceptibility or by the presence of irreversibility on all time scales. However these quantities are not accessible at $`T=0`$, and so other measures of ordering are necessary. One characteristic of a spin glass order is sensitivity to external perturbations; one can thus see whether the ground state evolves chaotically with $`B`$ up to a critical value $`B_c`$. The study in Houdayer and Martin (1999a) showed that there was no sign of chaos when the lattice size was increased for values of $`B1`$, suggesting that $`0B_c<1.0`$. This should be compared to the mean field value for connectivity $`6`$ random graphs for which $`B_{AT}(T=0)2.1`$. Another study Krzakala et al. (2001) considered how the ground state responded to a domain wall twist; this gave some evidence that the system had frozen order for $`BB_c0.6`$, although $`B_c=0`$ could not be excluded. Our approach here probes the spin glass ordering through its local excitations; the energy to flip a cluster of spins behaves differently in a spin glass phase and in a paramagnetic phase, and we use this to test for spin glass ordering. First we shall describe what these minimum energy clusters are; then we shall consider how their energies and magnetizations depend on the external magnetic field.

*Minimum energy clusters —* We start with the $`3d`$ Edwards-Anderson (EA) model with periodic boundary conditions. The Hamiltonian is defined on a cubic lattice of $`N=L^3`$ spins,

$$H=\underset{<ij>}{}J_{ij}S_iS_j+B\underset{i}{}S_i.$$
(2)
The spins are Ising, i.e., $`S_i=\pm 1`$, and the nearest-neighbor interactions $`\{J_{ij}\}`$ are quenched random variables distributed according to a Gaussian law with zero mean and unit variance.

In this model, when $`B`$ is set to $`0`$, there is numerical evidence Palassini and Caracciolo (1999) of an ordered phase at low temperature where each spin $`i`$ is frozen in a random direction $`S_i0`$. We are interested in local excitations in the putative frozen phase, so following previous work Lamarcq et al. (2002), we define a Minimum Energy Cluster or mec as the connected cluster of spins of lowest energy that contains a specified number of spins and a given, arbitrarily chosen site. The study of such mec in zero field showed Lamarcq et al. (2002) that they are, perhaps surprisingly, fractals whose typical energy decreases with their size, in contrast with the properties expected for the Fisher-Huse droplets. This result is striking and we take it to be a signature of the spin glass phase, allowing us to probe now the case where $`B>0`$.

*Numerical method —* Our measurements are performed on lattices with $`L=6`$ and $`L=10`$ with different values of the field taken in the range $`0B1.5`$; in all cases, we generated $`100`$ disorder samples which was enough to obtain reasonably small statistical errors on our observables.

For each disorder sample, we first compute the ground state of the system using a genetic renormalization algorithm Houdayer and Martin (1999b). Then, we choose an arbitrary “reference” spin and flip it along with a cluster containing $`V1`$ other spins connected to it. The goal is to find the mec of size $`V`$, with the constraint that the reference spin is held flipped and the cluster is always connected. To find that mec, we use non-local Kawasaki dynamics as in Lamarcq et al. (2002) and exchange Monte Carlo Hukushima and Nemoto (1996) with a set of 30 temperatures (between $`T_1=0.5`$ and $`T_{30}=3.0`$). This search for the mec is repeated five independent times to check whether the same cluster is found. (Since our algorithm is heuristic, we obtain an upper bound rather than the true energy of the mec.) In practice, we find the mec quite reliably for $`V32`$ as shown by our tests. Note that our optimization procedure has been improved since the analysis made in Lamarcq et al. (2002) in zero field; for instance, we have been able to find, for $`V=64`$, lower excitations than before. With this improvement, we confirm our earlier conclusions but also we have been able to go to larger mec sizes. In all that follows, we take $`V_n=2^n,1n6`$; even though the true mec are not always found at $`V=64`$, we present the results also for that size as we believe the corresponding bias is small.

*Geometry of mec* Let us look at two geometric properties of mec. Consider first the extension of mec as a function of their volume $`V`$; our data for the mean end-to-end extensions are displayed in the insert of Fig. 1. Given these values, we conclude that mec span the whole system when $`V15`$ if $`L=6`$ and when $`V30`$ if $`L=10`$.

Second, consider the fractal dimension of the mec. When $`B=0`$, we argued Lamarcq et al. (2002) that configurational entropy is sufficient to drive mec towards the lattice animal phase in which they are stringy objects with a fractal dimension $`d_f=2`$. This excitation branch is thus distinct from the Fisher-Huse droplets which are compact objects with a fractal boundary. What happens in the presence of a field? It is useful to think of the large $`B`$ limit which is strongly paramagnetic. The surface energy of a cluster is a random variable $`J_{ij}`$ that is symmetric about zero, while its bulk energy is $`BV`$. Clearly, energy minimization in a random environment will drive the clusters to extend their surface, leading to the lattice animal phase. We can see in Fig. 1 that this geometrical property, just like the extension, is nearly independent of the strength of the magnetic field. This shows that geometrical properties of mec are not relevant to distinguish between a paramagnet and a spin glass, but it also gives further support to the claim that mec are lattice animals even when $`B=0`$.

*Energies of mec* Consider first the typical excitation energy $`\overline{E}`$ of a mec, where $`(\overline{})`$ denotes the average over disorder. In the paramagnetic phase, i.e., for large $`B`$, $`\overline{E}`$ is expected to be self-averaging as $`V\mathrm{}`$. Since there is a non zero magnetization per spin, that energy should grow as $`\overline{E(V)}BV`$, leading therefore to $`\theta =3`$. Note however that because of the self-averaging property, the probability to find $`E=O(1)`$ should go to zero much faster than any power of $`V`$ for large $`V`$; there are thus no low-energy excitations at large scales. Following an Imry-Ma argument, Fisher and Huse Fisher and Huse (1986) argued that this should actually occur for any $`B0`$ for sufficiently large $`V`$’s. The paramagnetic behavior should set in beyond a cross-over volume $`V_B`$, obtained by comparing the Zeeman energy $`B\sqrt{V}`$ of a $`B=0`$ droplet with its excitation energy $`\mathrm{{\rm Y}}V^{\theta /d_f}`$. The result is the cross-over volume $`V_BB^{2d_f/(d_f2\theta )}`$, with the corresponding length being given in Eq. 1. Note that $`d_f`$ is the (possibly fractal) dimension of the excitation. In the usual picture of compact droplets Fisher and Huse (1987), one has $`d_f=3`$ and $`\theta 0.2`$. In the numerical investigation of mec Lamarcq et al. (2002), it was found that the excitation clusters were instead fractal, $`d_f2.0`$, and that the effective $`\theta `$ was small and negative. Since $`\theta `$ is small compared to $`d_f`$ in all cases, we expect roughly $`V_B(\mathrm{{\rm Y}}/B)^2`$. Using the value $`\mathrm{{\rm Y}}6`$ from our numerical results, we obtain a rather large value $`V_B36/B^2`$.

Let us focus on the volume dependence of the mean mec energy at different fields $`B`$. For all our $`B`$ values, this mean initially decreases with volume; furthermore, this decay is compatible with a power law, $`\overline{E(V)}V^{\theta _f/d_f}`$ with $`\theta _f0.13`$ as found in zero field Lamarcq et al. (2002). From the point of view of the mean field picture, the behavior observed in our data is qualitatively as expected: below $`B_c0.5`$, the spin glass ordering leads to a decreasing $`\overline{E(V)}`$ at all $`V`$, while for $`B>B_c`$ we recover the paramagnetic behavior of an increasing energy when $`V`$ is sufficiently large. Note however that the estimate $`B_c0.5`$ is far below the mean field value Houdayer and Martin (1999a) associated with random graphs of connectivity $`6`$, $`B_{MF}2.1`$. The important issue is therefore whether the curves eventually bend upwards at large $`V`$ for all values of $`B>0`$. In order to be more quantitative, we have attempted to scale the data, looking for a collapse when plotting $`\overline{E(V)}/\overline{E(V_B)}`$ as a function of $`V/V_B`$. Surprisingly, the data does not collapse well at all when $`V_BB^2`$ but merges much better if we take $`V_BB^1`$ as shown by the “scaling curve” in the insert of Fig. 2. This result suggests that even the low $`B`$ curves might eventually bend upwards for large enough $`V`$; that could be interpreted as signaling the instability of the spin-glass phase in a field. We will return to this point later.

Note that the $`BV`$ scaling obtained above is rather unexpected since it means that the influence of a small magnetic field on low energy excitations is actually stronger than anticipated by the Fisher-Huse argument. Although we do not have a clear understanding of this result, one could slightly alter the argument as follows. For $`B=0`$ the ground state has a zero magnetization per spin, which means, as recalled above, that the magnetization of a region of volume $`\mathrm{}^3`$ is typically of order $`\mathrm{}^{3/2}`$. For $`B0`$, a non zero magnetization per spin $`\chi B`$ (see Fig. 3) appears; this means that the spatial correlations in the ground state change for distances larger than $`\mathrm{}_B`$ such that $`\mathrm{}_B^3\chi B\mathrm{}_B^{3/2}`$, leading to $`\mathrm{}_B(\chi B)^{2/3}`$. Since the mec are fractal, the excitations built at $`B0`$ are expected to be affected as soon as $`V\mathrm{}_B^{d_f}`$, or, using $`d_f=2`$, at $`V(\chi B)^{4/3}`$, a result that is closer to the scaling reported in Fig. 2.

We have also studied the distribution of energies of mec, for a given volume; the case $`V=32`$ is given in Fig. 4. As expected, for large fields ($`B=1.0`$ and $`1.5`$), the distributions become symmetric around their maximum, and their weights for small excitation energies rapidly tend to zero. However for low fields, $`B0.5`$, we find that these distributions hardly change at all with $`B`$.

*Magnetization of mec* Another observable of interest is the mean magnetization per site $`m(V,B)`$ of mec of volume $`V`$ in the presence of the magnetic field $`B`$. (The mec’s magnetization is defined before it is flipped.) In the droplet picture, the system is driven to the large $`B`$ limit when $`V`$ grows as soon as $`B0`$. Thus, beyond the cross-over scale $`V_B`$ discussed above, $`m(V,B)`$ should converge to a function $`m(B)`$ that is non-zero for all $`B0`$. The situation is different in the mean field picture. Mean field predicts that the equilibrium magnetization is independent of the temperature $`T`$ at low enough $`T`$. Now since the low $`T`$ properties follow from those of the lowest-lying excitations, we see that these excitations cannot have a non-zero mean magnetization. The mean field picture then consists in having the large $`V`$ limit of $`m(V,B)`$ vanish for $`BB_c`$ and grow for $`B>B_c`$. In both pictures of the spin glass phase, $`m(B)`$ acts as an order parameter for the paramagnetic phase. Fig. 5 shows $`m(V,B)`$ as a function of $`V`$. We see that the magnetizations for $`B0.5`$ are indeed very small and in fact are negative rather than positive. Also, we see a clear upturn to positive $`m`$ for $`V=32`$, $`B0.4`$; this upturn seems to be present even for smaller fields, suggesting that $`B_c<0.4`$. This bound is more stringent than what we found from the energies of the mec. Finally, just like the energy, the magnetization appears to approximately rescale as a function of $`BV`$.

*Conclusions —* Our results show that when $`B<0.4`$, the system behaves qualitatively just as in the case $`B=0`$ (at least for the sizes studied here), whereas significant changes arise when $`B0.4`$. If there is a critical field, which appears to us rather unplausible, it must obey $`B_c<0.4`$, a bound that is stronger than that most recently obtained Krzakala et al. (2001). This value is small compared to the value arising in the mean field approximation for which $`B_c2.1`$. Although we cannot completely rule out the mean field picture, the possibility of rescaling the different curves as a function of $`BV`$ points towards the destruction of the spin-glass phase for $`B0`$. On the other hand, our data are also unexpected from the point of view of the droplet model. First, as emphasized before Lamarcq et al. (2002), mec appear as fractal, rather than compact, objects. Second, the Fisher-Huse prediction of a crossover volume $`V_BB^2`$ appears not to be correct. We find that the magnetic field has a much stronger effect, since we obtain $`V_BB^1`$, which might be a consequence of the fractal nature of the mec. Finally, the geometric properties of the mec appear to be $`B`$ independent, and most probably they are lattice animals for all $`B`$. Many of the points raised by the present study would be clarified if the corresponding results in $`d=4`$ were available. Note that recent work Krzakala and Parisi (2003) on Bethe lattices reveals a non zero value of the critical field $`B_c`$ as expected.

*Acknowledgments —* J. L. acknowledges a fellowship from the MENRT and the CEA for computer time on a Compaq SC232. O. M. thanks the SPEC for its hospitality during this work. We thank F. Krzakala for communicating his results prior to publication.