text:
# Toward Making the Constraint Hypersurface an Attractor in Free Evolution

## I Introduction

The data analysis needs of the LIGO and LISA experiments have driven an interest in finding numerical solutions to the Einstein equations in astrophysically relevant situations. Wave forms from such numerical simulations can, once generated, be used as templates for finding gravitational wave signals in the experimental noise.

In order to find numerical solutions, however, one must make a 3+1 space-time split of the Einstein equations following some variation of the ADM formalism Arnowitt et al. (1962). When using free evolution techniques, which is the most common practice in numerical relativity today, a plethora of authors have found that the exact form of the evolution equations, in particular how the constraints are substituted into the evolution equations, plays a critical role in how stable the evolution will be (e.g. Shibata and Nakamura (1995); Baumgarte and Shapiro (1999); Alcubierre et al. (2000); Kidder et al. (2001); Laguna and Shoemaker (2002)). Although making such modifications to the evolution equations remains more of an art than a science, some authors have proposed well defined methods for making stability-friendly modifications.

Historically, Detweiler Detweiler (1987) seems to be one of the first authors to consider explicitly a method for adding terms proportional to the constraints to the evolution equations with the goal of improving simulation stability. He looked at the evolution of the constraints, and tried to add terms to the evolution equations of the fundamental variables that would drive the constraints to zero as the evolution progressed. This approach is very similar to the approach presented here.

Where the approach differs, however, is significant. Whereas Detweiler needed to add terms such that a negative-definite operator appeared on the right hand side of the evolution equation for the constraints, the method proposed here generates a general term of the correct sign which is added to terms of unknown sign that come from the underlying formalism. Under special circumstances, Detweiler showed that his method does in fact generate a negative-definite operator by making use of a free parameter introduced in his correction terms. The method proposed here also uses a free parameter in arguing that the constraints can be evolved toward zero, with the difference that the argument here is more general but without a guarantee that the overall sign of the the constraint evolution has the correct sign.

The $`\lambda `$-system approach, proposed by Brodbeck, Frittelli, Hübner, and Reula, embeds the Einstein equations into a larger system of symmetric hyperbolic equations. The Einstein equations are coupled to the larger system by adding terms to the Einstein evolution equations which are zero on the constraint hypersurface, with the extra degrees of freedom defined such that they are forever zero if the Einstein constraints and constraints on the new degrees of freedom are satisfied in the initial data. A linear order analysis of this system suggests that for “small enough” deviations from the constraint hypersurface, the system should asymptote to a solution of the Einstein equations Brodbeck et al. (1999).

Finally, Yoneda and Shinkai Yoneda and Shinkai (2001) performed detailed eigenvalue calculations to determine the theoretical signs of the constraint evolution when linearizing around a Minkowski background spacetime. While their program was intensive and thorough, their results depend on the validity of a perturbation expansion around a fixed spacetime. In that limit, they evaluate the sign of the eigenvalues of their expansion matrix, and use those signs to predict the goodness of correction terms which are strictly linear in the constraints and spatial derivatives of the constraints. Although one of their more promising terms was used to good effect in a numerical experiment Yo et al. (2002), their analysis remains strictly valid only in the perturbation regime around Minkowski space, and, at present, they have only examined terms linear in the constraints and derivatives of the constraints.

The method here, in contrast, attempts to make the constraint surface an attractor for the evolution equations without appealing to perturbation theory in any way and without extending the system of equations. The results will show it successful in reaching this goal for the simple case of the Maxwell equations, but cannot confirm with certainty that the results will generalize to non-linear general relativity. The primary problem is that the correction terms proposed here have the form

$$\frac{d}{dt}\mathrm{constraint}^2=(\mathrm{orig}\mathrm{rhs})(\mathrm{correction})$$

with a correction term of definite sign subtracted from a term of indefinite sign coming from the underlying theory. The success of the method in the case of the linear Maxwell equations makes a similar result in the linearized Einstein equations likely, but does not guarantee success for full, non-linear general relativity. This will be discussed more completely in Section IV.

In this paper, following earlier work by Knapp, Walker, and Baumgarte (KWB) Knapp et al. (2002), I will apply the method to Maxwell’s equations as a proof-of-concept test. The work by KWB showed that some qualitative features of two popular ADM formulations of general relativity are manifest in the simpler framework of electromagnetism, making this a natural place to demonstrate the feasibility of the new methods presented here.

## II The Theory

I will introduce the idea here by first considering the ordinary differential equation example of the simple harmonic oscillator, and then considering the more interesting case of partial differential equations.

### II.1 ODE: Harmonic Oscillator

The equations of motion for a harmonic oscillator (with $`\omega =1`$) are

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`v`$ (1)
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`x`$ (2)
and the constraint 111This is actually a conserved quantity. For this simple ODE case, there is no difference, but when we pass to PDEs, there will be a difference between conserved quantities and constraints. is the shifted energy function

$$E(t)=x^2+v^2E_0$$
(3)
with $`E_0=x^2(0)+v^2(0)`$, which should be zero-valued at all times. In a numerical evolution, one would like to evolve the fundamental quantities $`x`$ and $`v`$, and, after doing so, calculate the constraint $`E(t)`$ from $`x`$ and $`v`$. Since $`E`$ is derived, on the other hand, it makes sense to look at its time evolution in terms of the time evolution of $`x`$ and $`v`$ using the chain rule. Furthermore, $`E(t)`$ is not bounded on either side, so consider instead $`E^2(t)`$, which should also always be zero. Applying the chain rule gives the elementary result that

$$\frac{dE^2}{dt}F(x,v)=\frac{E^2}{x}\frac{dx}{dt}+\frac{E^2}{v}\frac{dv}{dt}.$$
(4)
Now we cannot change the two partial derivative factors since they are fixed by the form of the energy expression, but, as noted before, we are free to change the equations of motion without changing the physics by adding multiples of the constraints. Let $`K`$ be a positive constant. If we make the changes

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`vK{\displaystyle \frac{E^2}{x}}`$ (5)
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`xK{\displaystyle \frac{E^2}{v}}`$ (6)
to the equations of motion, then

$$\frac{dE^2}{dt}=F(x,v)K\left[\left(\frac{E^2}{x}\right)^2+\left(\frac{E^2}{v}\right)^2\right]$$
(7)
shows how the constraint evolves under the new equations of motion. For large $`K`$, the term in brackets will dominate the evolution of the constraint, constantly pushing it toward its minimum value of zero.

Numerical results, integrating forward in time with the Euler method, are in agreement with this analytic result.

### II.2 PDEs

Having seen the method in the simple case of ODEs, consider now the the more interesting case of PDEs with one constraint equation. Let $`s`$ be the state vector for a system. Define

$$S_m(t,𝐱)=\frac{s_m}{t}$$
(8)
to be the right hand side of the evolution equation in the unmodified formalism. A general constraint $`C`$ will depend on $`s`$ and its spatial derivatives. Furthermore, the constraint should be satisfied at every point in space.

These considerations motivate looking at the integrated, squared constraint $`\overline{C^2}=C^2d^Nx`$, and, instead of taking partial derivatives with respect to the fields, we need to take variational derivatives of the integrated constraint when considering the analogies of Equations (56) so that the dependence of $`C`$ on the spatial derivatives of the fields is treated properly.

Following this prescription, the appropriate modification to the equation of motion for the state vector is

$$\frac{s_m}{t}=S_m(t,𝐱)K_{mn}(t,𝐱)\frac{\delta \overline{C^2}}{\delta s_n(t,𝐱)}$$
(9)
for some positive-definite matrix-valued function $`K_{mn}`$. Under this change,

$$\frac{d\overline{C^2}}{dt}=D[s]\left(\frac{\delta \overline{C^2}}{\delta s_m}\right)K_{mn}\left(\frac{\delta \overline{C^2}}{\delta s_n}\right)d^Nx$$
(10)
gives the evolution of the constraint in the modified theory. Here $`D[s]`$ gives the functional form of the right hand side of the constraint’s evolution equation in the unmodified theory. For the cases considered here, I chose the $`K_{mn}`$ diagonal and constant.

For systems with $`M`$ constraint equations, the method is easily modified by taking the grand constraint functional to be

$$\overline{C_G^2}=w_{IJ}(t,𝐱)C_IC_Jd^Nx$$
(11)
with any positive definite matrix $`w_{IJ}`$. Like $`K_{mn}`$, the matrix $`w_{IJ}`$, can, in principle, be a function of both space and time, and can have off-diagonal entries. In practice, however, I have only used diagonal and constant matrices. Even in the diagonal case, the matrix is necessary because there is no a priori reason to believe that the constraints have the same dimensions, and it also allows the different constraints to be treated with different relative importance. In addition, since there is no natural scale for the grand constraint, one may always set one of the coefficients $`w_{IJ}`$ in (11) to unity.

## III Application to Maxwell

This section takes Maxwell’s equations as a concrete example of a system of partial differential equations subject to a pointwise constraint. Following KWB Knapp et al. (2002), I consider two ways of writing Maxwell’s equations for the vacuum in terms of the vector potential $`A_i`$. The first system, called System I, uses the evolution equations

$`_tA_i`$ $`=`$ $`E_i_i\psi `$ (12)
$`_tE_i`$ $`=`$ $`_j_jA_i+_i_jA_j`$ (13)
and the constraint

$$C_E=_iE_i=0.$$
(14)
The second system introduces the additional field $`\mathrm{\Gamma }`$ defined by

$$\mathrm{\Gamma }=_iA_i$$
(15)
to eliminate mixed derivatives in (13). The evolution equations for System II are

$`_tE_i`$ $`=`$ $`_j_jA_i+_i\mathrm{\Gamma }`$ (16)
$`_t\mathrm{\Gamma }`$ $`=`$ $`_i_i\psi `$ (17)
and (12). Both systems use a gauge consistent with

$$_t\psi =_iA_i=\mathrm{\Gamma }$$
(18)
using the first equality for System I and the second for System II.

### III.1 System I Evolution Equations

Having defined the systems, I would now like to calculate the terms required for applying the constraint finding method to System I. Here there is only one constraint, $`C=C_E=_iE_i`$, which is zero-valued. It depends only on the first derivatives of the electric field, therefore I need only to calculate

$$\frac{\delta \overline{C^2}}{\delta E_i(𝐱)}=2_iC_E$$
(19)
which modifies (13). The new evolution equation for the electric field is

$$_tE_i=_k_kA_i+_i_kA_k+2K_E_iC_E$$
(20)
for an arbitrary positive $`K_E`$, while the other System I evolution equations (12) and (18) remain unchanged.

### III.2 System II Evolution Equations

System II, unlike System I, has two constraints that should be enforced, the original constraint given by (14) plus the definition of $`\mathrm{\Gamma }`$ in (15), rewritten as

$$C_\mathrm{\Gamma }=_iA_i\mathrm{\Gamma }=0$$
(21)
to make it zero-valued. This provides more freedom in constructing the grand constraint

$$C^2=C_E^2+wC_\mathrm{\Gamma }^2$$
(22)
where one does not necessarily have to treat the constraints on equal footing. In this case, the total constraint depends additionally on $`\mathrm{\Gamma }`$ and first derivatives of the vector potential. In addition to (19), which is still valid, I need

$`{\displaystyle \frac{\delta \overline{C^2}}{\delta A_i(𝐱)}}`$ $`=`$ $`2w_iC_\mathrm{\Gamma }`$ (23)
$`{\displaystyle \frac{\delta \overline{C^2}}{\delta \mathrm{\Gamma }(𝐱)}}`$ $`=`$ $`2wC_\mathrm{\Gamma }`$ (24)
to enforce the definition of $`\mathrm{\Gamma }`$.

Applying these correction terms to the evolution equations gives the new equations of motion

$`_tA_i`$ $`=`$ $`E_i_i\psi +2wK_A_iC_\mathrm{\Gamma }`$ (25)
$`_tE_i`$ $`=`$ $`_k_kA_i+_i\mathrm{\Gamma }+2K_E_iC_E`$ (26)
$`_t\mathrm{\Gamma }`$ $`=`$ $`_k_k\psi +2wK_\mathrm{\Gamma }C_\mathrm{\Gamma }`$ (27)
which, combined with (18), form a complete system. The constants $`K_A`$ and $`K_\mathrm{\Gamma }`$ are arbitrary but positive.

### III.3 Propagation of Constraints

In Knapp et al. (2002), KWB examined the evolution equation for $`C_E`$ in both systems. They showed that for System I, the constraint does not evolve in time, and that in System II, the constraint obeys a wave equation. They did not have reason, however, to consider the time evolution of the secondary constraint $`C_\mathrm{\Gamma }`$. I extend their results here by showing that the secondary constraint also satisfies a wave equation in the unmodified case. Furthermore, I demonstrate the improved behavior of the constraints under the modifications proposed here.

Calculating the first time derivatives of the constraints is easily accomplished by taking the time derivatives of (14) and (21) and replacing the time derivatives that appear on the right hand sides of the equations by the evolution equations (20), (25), and (27). This gives the results

$`_tC_E`$ $`=`$ $`_i^2\left[C_\mathrm{\Gamma }2K_EC_E\right]`$ (28)
$`_tC_\mathrm{\Gamma }`$ $`=`$ $`C_E+2w\left[K_A_i^2C_\mathrm{\Gamma }K_\mathrm{\Gamma }C_\mathrm{\Gamma }\right]`$ (29)
which can be viewed as valid for both systems if $`C_\mathrm{\Gamma }`$ is taken to be identically zero for System I.

To see that KWB have a wave equation for the secondary constraint, set all of the $`K`$s to zero in (2829) to eliminate the modifications, and take the time derivative of (29). The result

$$_t^2C_\mathrm{\Gamma }=_i^2C_\mathrm{\Gamma }$$
(30)
follows immediately.

Of greater interest here, however, is an analysis on the modified equations (2829) in their first derivative form. The equations are linear, so they admit a Fourier analysis by substituting a plane wave solution $`e^{ikx}`$ into the right hand sides. After substituting, the resulting equations

$`_tC_E`$ $`=`$ $`k^2C_\mathrm{\Gamma }2K_Ek^2C_E`$ (31)
$`_tC_\mathrm{\Gamma }`$ $`=`$ $`C_E2w\left[K_Ak^2+K_\mathrm{\Gamma }\right]C_\mathrm{\Gamma }`$ (32)
retain the terms that gave the KWB wave equations for the constraints, but have additional terms that look like they provide exponential decay. This system of equations is, in fact, simple enough for Maple to solve analytically for general values of the $`K`$s in one dimension. The solution, which is too long to display in detail here, consists of a sum of terms with the form

$$\mathrm{exp}\left[\left(f_1^+\pm \sqrt{\sigma \left[\left(f_1^{}\right)^2k^2\right]}\right)t\right]f_2$$
(33)
where

$$f_1^\pm (k^2)=(K_E\pm wK_A)k^2\pm wK_\mathrm{\Gamma }$$
(34)
$`\sigma =\pm 1`$, and $`f_2=f_2(k,C_i(0,x))`$ is some simple function of $`k`$ and the initial values of the constraints. Since $`f_1^+(k^2)>0`$ is manifestly positive and the radical is either positive or pure imaginary, the only term of this form that could cause anything other than exponential decay is $`f_1^++\sqrt{\sigma [(f_1^{})^2k^2]}`$ for parameters where $`\sigma [(f_1^{})^2k^2]>0`$. Simple algebraic analysis, however, shows that even this term has an overall minus sign, giving exponential decay.

In order to make the argument more concrete, I present here the $`k=1`$ solution

$`C_E(t,x)`$ $`=`$ $`e^{3t}\left[C_E(0,x)+S(x)t\right]`$ (35)
$`C_\mathrm{\Gamma }(t,x)`$ $`=`$ $`e^{3t}\left[C_\mathrm{\Gamma }(0,x)S(x)t\right]`$ (36)
for which all of the $`K`$s are set equal to one. Here $`S(x)=C_E(0,x)+C_\mathrm{\Gamma }(0,x)`$ is a short hand. The general solution can be requested from the author in the form of a Maple worksheet.

### III.4 Numerical Results

My numerical experiments on these modified systems of equations used an ICN integration scheme Teukolsky (2000), and a Courant factor of 1/2. The spatial domain ran from -6 to +6, with data stored on 99 points in each coordinate direction. On the evolved fields ($`E_i`$, $`A_i`$, $`\psi `$, and $`\mathrm{\Gamma }`$), I imposed outgoing wave boundary conditions (see Alcubierre et al. (2003) for implementation details), and on the constraints, I imposed $`C_I=0`$ for applicable $`I`$. All runs were performed on a 500 MHz Digital Personal Workstation with 1.5 GB of RAM.

For each system, I ran three parameter sets, one of which reproduced the equations used by KWB. The full definitions of all of the parameter sets are found in Table 1.

Sensible values of the parameters were easily determined by trial and error. Choosing the values too small, as expected, makes little difference in the evolution, while choosing the values too large leads to numerical instabilities during the transient period of the evolution 222One possible explanation for the instability associated with large parameters is that the added terms modify the dispersion relationships for the various Fourier modes, as seen in (33). Large values of the parameters may require adjustments to the Courant condition, which has not been analyzed here..

I followed KWB in using the analytic solution

$`A^\varphi `$ $`=`$ $`0`$ (37)
$`E^\varphi `$ $`=`$ $`8𝒜\lambda ^2re^{\lambda r^2}\mathrm{sin}\theta `$ (38)
of a toroidal dipole to generate the initial data. The other components of the fields are zero. I chose $`\lambda =𝒜=1`$, and the conversion from spherical to Cartesian coordinates was made in the code.

The results of the System I runs are summarized in Figure 1, which shows a plot of $`C_E_2`$ vs. $`t`$.

The I-0 (control) curve reproduces the findings of KWB that, after an initial transient, the $`C_E`$ constraint does not evolve in time. The I-1 and I-2 curves, representing different values (see Table 1) of the parameter $`K_E`$, on the other hand, show a modulated exponential decay. Eventually, around $`t200`$, the I-1 case also stops decaying as rapidly, while the rapid exponential decay continues through the end of the run for I-2.

That the constraint in the modified case I-1 ceases to evolve at some point is consistent with (10), which implies that the constraints will cease to evolve when the first term balances with the second term. From (10), one expects that this balance will be achieved for smaller constraint violation when the driving parameter is larger, which is consistent with the results shown in Figure 1.

Looking at two dimensional slices of the constraint data at various times, also suggests that the source of the modulation in the decay demonstrated by I-1 and I-2 is fluctuations at the boundary, possibly caused by the simple boundary condition applied there on the constraints. Significantly, these fluctuations are unable to penetrate the interior of the computational domain, unlike many scenarios seen in numerical relativity where noise from the boundary noticeably propagates inward from the outer boundary, eventually killing the simulation. Because I am only interested in the Maxwell equations as a test-bed for the method, and because the modified System I equations already perform orders of magnitude better than their unmodified counter part, I have not pursued this point further.

Figure 2 shows a plot of $`C_E_2`$ vs. $`t`$ for the System II case.

Here again, the control run (II-0) reproduces the results of KWB, this time showing exponential decay in the primary constraint. Even with such an ideal result in the unmodified case, the modified runs II-1 and II-2 show improvement. They represent runs with non-zero values (see Table 1) of the various forcing parameters, and in these cases the constraint decays exponentially, but with a smaller characteristic time.

Figure 3 is a plot of $`C_\mathrm{\Gamma }_2`$ versus time $`t`$, and shows similar behavior to the primary constraint in all three cases.

The secondary constraint shows exponential decay in the unmodified II-0 case, while showing a faster decay in the two modified runs, II-1 and II-2. The jump in the graph at $`t=0`$ occurs because the secondary constraint is exactly satisfied in the initial data by construction.

It should be noted that one likely explanation for the extremely favorable performance of the unmodified System II is that, since the constraints satisfy a wave equation, constraint violations propagate off of the grid. This is supported by the data presented by KWB, showing that decay rate of the constraint decreases when the outer boundary is moved farther out Knapp et al. (2002). The method presented here benefits from constraint violations propagating off of the grid as well, but, in addition, attempts to damp the constraint violation throughout the grid at all times.

In both systems, there is a practical limit to how large the forcing parameters may be chosen without creating an instability. I believe that its origin can be seen in (33), which shows how the dispersion relation for the Fourier components of the solutions is modified by the new terms. This modification is effectively a modification of the propagation speed of those components, which, for a fixed Courant factor, can lead to numerical instabilities in the numerical scheme. When the forcing parameters are too large, the code crashes. The constraint blows up rapidly immediately before such a crash, and such crashes usually occur during the transient period of the numerical scheme.

## IV Conclusions and Discussion

The method proposed here does not provide a formalism for any physical system. It attempts to take any formalism and makes that formalism better. The analysis does not depend on the form of the equations and the arguments leading to the modifications do not depend on the validity of any perturbative expansion. The results in the case of the Maxwell equations are very encouraging and consistent with the analytic predictions.

Should the results derived from Maxwell’s equations generalize to the Einstein equations, furthermore, this could prove extremely beneficial to the long term stability of simulations in numerical relativity. Fourier analysis of the ADM system suggests that instabilities seen in various 3+1 formulations are triggered by exponentially growing modes not seen in the initial data Calabrese et al. (2002). On the other hand, at least for the case of the Maxwell equations, the terms added via the method proposed here lead to exponential damping terms in the evolution equations for all Fourier modes.

It should be noted, however, that in some sense the Maxwell equations were an optimal case. Notice, by counting the number of integrations by parts required to calculate the correction terms to the evolution equations, that the highest order of spatial derivatives in the new evolution equations will be the larger of (1) order of the original equations and (2) twice the order of the constraints. This means that, while the order of the Maxwell evolution equations was not changed, the Einstein equations will change from second order to fourth order in space. How this change of order effects the numeric solutions to the equations must be studied as this method is applied in the context of numerical relativity.

###### Acknowledgements.

I would like to thank C. Misner, C. Schiff, T. Jacobson, D. Mattingly, and V. Prieto-Gortcheva for comments on the ideas presented here and on the manuscript. The code used in the numerical simulations made use of the Cactus computational toolkit and a modified piece of code originally written by E. Schnetter for the Maya project. This work was supported in part by NSF grant PHY-0071020.


format: text

text:
# Diffusion-limited aggregation in channel geometry

## Abstract

We performed extensive numerical simulation of diffusion-limited aggregation in two dimensional channel geometry. Contrary to earlier claims, the measured fractal dimension $`D=1.712\pm 0.002`$ and its leading correction to scaling are the same as in the radial case. The average cluster, defined as the average conformal map, is similar but not identical to Saffman-Taylor fingers.

Diffusion-limited aggregation (DLA) has attracted considerable attention since its introduction by Witten and Sander in 1981 Witten and Sander (1981). In this model an aggregate or cluster grows by capturing diffusing particles which irreversibly attach to it on first contact. This is the discrete model of a wide variety of physical systems in the Laplacian growth class. This class can be modeled with a field, satisfying the Laplace equation outside of a growing cluster, where the cluster grows in proportion to the gradient of the field at the boundary. A generalization of DLA, known as the dielectric breakdown model (DBM) Niemeyer et al. (1984), allows growth proportional to the field gradient to the exponent $`\eta `$. DLA is regained for $`\eta =1`$. The majority of our current knowledge about DLA is numerical (mostly in two dimensions and radial geometry), although progress has been made in the theoretical front as well (see eg. Ball and Somfai (2002, 2003)).

One of the controversies surrounding DLA in channel geometry has been that the fractal dimension might differ from that in the radial case. This claim has been based on small size simulations Meakin and Family (1986); Argoul et al. (1988); Evertsz (1990) or small size calculations Kol and Aharony (2000, 2001). In this paper, based on extensive numerical simulations of off-lattice DLA in channel, we show that the fractal dimension is asymptotically the *same* in the two geometries.

One of the most important differences between the two geometries is that for the channel the continuum version of the problem (Laplacian growth) has a stable solution without tip splitting instability and finger competition. These stationary translating solutions —called Saffman-Taylor fingers Saffman and Taylor (1958)— have been studied in viscous fingering experiments in Hele–Shaw cells Thome et al. (1989). In this paper we compare them with average profiles of DLA clusters.

Our first method for generating DBM in a channel is to use iterated conformal maps. The conformal mapping method for the radial case is described in Hastings and Levitov (1998) and Hastings (2001). For radial DLA, a map is created from the unit circle in the $`w`$ (“mathematical”) plane to the unit circle with a bump at a randomly chosen angle in the physical plane; the composition of such maps is a map from the unit circle the DLA cluster. To produce DBM clusters, instead of choosing the angles at random, we use a monte-carlo method to select bump sites with the correct distribution Hastings (2001).

To adapt this method to a channel, we modify the map of Stepanov and Levitov Stepanov and Levitov (2001) for both periodic and reflective boundary conditions by requiring the map to be symmetric about the real axis. The map is given by

$$f_{\mathrm{\Lambda },\theta }(w)=ln\left[g^1(\stackrel{~}{f}_\mathrm{\Lambda }(g(w)))\right],$$
(1)
where $`g(w)=\frac{w1}{w+1}`$ is a map from the unit circle to the imaginary axis, and from the exterior of the unit circle to the positive-real half-plane. The function

$$\stackrel{~}{f}_\mathrm{\Lambda }(w)=\frac{w+\gamma \sqrt{(wxi)^2+\mathrm{\Lambda }^2}+\gamma \sqrt{(w+xi)^2+\mathrm{\Lambda }^2}}{1+\gamma \sqrt{(1xi)^2+\mathrm{\Lambda }^2}+\gamma \sqrt{(1+xi)^2+\mathrm{\Lambda }^2}}$$
(2)
adds two bumps at symmetric points $`x^{}i`$ and $`x^{}i`$. The denominator in Eq. (2) forces $`\stackrel{~}{f}_\mathrm{\Lambda }(w)`$ to map $`w=1`$ to $`1`$, so that $`f_{\mathrm{\Lambda },\theta }(w)`$ maps infinity to infinity. The parameter $`\mathrm{\Lambda }`$ controls the size of the bump, and $`\gamma `$ controls the aspect ratio of the bump. Since $`g(w)`$ maps $`1`$ to infinity, we choose $`\theta `$ at random between $`0`$ and $`\pi `$; and for $`\theta >\frac{\pi }{2}`$, $`f_{\mathrm{\Lambda },\theta }(w)`$ becomes

$$f_{\mathrm{\Lambda },\theta }(w)=\overline{g^1(\stackrel{~}{f}_\mathrm{\Lambda }(g(\overline{w})))}$$
(3)
where the bar denotes complex conjugation.

The actual bump positions are not at $`\pm xi`$, but are off by a small factor determined by $`x`$ and $`\mathrm{\Lambda }`$. The bump size is also dependent on $`x`$ and $`\mathrm{\Lambda }`$. In order to get bumps at angle $`\theta `$, $`\stackrel{~}{f}_\mathrm{\Lambda }`$ must place bumps at $`g(e^{\pm i\theta })=\frac{sin\theta }{1+cos\theta }i`$. We do this by an approximation method. To keep all the particles in the cluster the same size, the bump size on the unit circle is varied according to the first derivative of the composite conformal map in the original version of the conformal map technique Hastings and Levitov (1998) This assumes that the higher order derivatives are negligible, which is not true deep inside a fjord. Thus particles added in a fjord can end up being very large, and sometimes can partially fill the channel. To combat this effect, we measure the bump area at each step, and iteratively correct the size parameter $`\mathrm{\Lambda }`$ if the area is outside of a preset tolerance (10% for the results in this paper); compare Stepanov and Levitov (2001).

While conformal mapping allows one to grow DBM for any $`\eta `$, and directly produces a conformal map for the cluster boundary, it is computationally intensive. A more efficient numerical algorithm for generating off-lattice DLA (that is, $`\eta =1`$) is a simple adaptation of hierarchical maps Ball and Brady (1985) to channel geometry. This method enables close to linear dependence of computing resources on cluster size. We used in total $`1.7\times 10^{11}`$ particles for the dimension calculations, and (including probes) $`4\times 10^{11}`$ particles for the average profile.

Both periodic and reflective boundary conditions have been implemented on the sides of the channel (the periodic boundary condition is sometimes referred to as “cylindrical”). The reflective boundary conditions is achieved as above: the cluster is grown in a channel of double width and periodic boundary conditions, and for each deposited particle we deposit also its mirror image. At the end one of the images was discarded. By the conventions used in this paper the channel is given by the range $`w/2<y<w/2`$ and the clusters grow (macroscopically) in the positive $`x`$ direction.

*Fractal dimension of channel DLA.* The fractal dimension is measured through the density. The average density scales with the width $`w`$ of the channel, with exponent given by the co-dimension:

$$\rho (w)w^{D2}$$
(4)
To avoid transients, we discarded the first and last $`3w`$ long section of the clusters, and measured the density (number of particle centers per area) on the remaining middle section.

We generated clusters of $`8\times 10^6`$ to $`32\times 10^6`$ particles in channels of width $`w=50,100,200,500,1000,2000,5000`$ particle diameters. For each width the number of clusters grown ranged from a few hundred to a few thousand, with more and larger clusters necessary for large widths, to achieve comparable statistical confidence in the average density.

Figure 1 shows the width ($`w`$) dependence of the effective fractal dimension $`D_{\text{eff}}=2+d\mathrm{ln}\rho /d\mathrm{ln}w`$. The fractal dimension tends to $`D=1.712\pm 0.002`$, independent of the choice of boundary conditions. Off-lattice noise reduction Ball et al. (2002) does not change the dimension, but accelerates the convergence to its asymptotic value.

*Average profile.* Analytical solutions for unbranched Laplacian growth in a channel with reflective boundary conditions have been known for a long time. One can find solutions which translate a fixed profile along the channel in time. In the absence of surface tension, these solutions form a one-parameter family, called Saffman-Taylor (ST) fingers Saffman and Taylor (1958). They are parametrized by the asymptotic ratio $`\lambda `$ of the widths of the finger ($`w_{\text{finger}}`$) and the channel ($`w`$), and have the profile:

$$x(y)=\frac{w(1\lambda )}{2\pi }\mathrm{log}\left[\frac{1}{2}\left(1+\mathrm{cos}\frac{2\pi y}{\lambda w}\right)\right]$$
(5)
Of these solutions $`\lambda =1/2`$ is the most important, because in related experiments Saffman and Taylor (1958) this profile has been observed in the limit of vanishing surface tension. Analytical calculations Combescot et al. (1986); Shraiman (1986); Hong and Langer (1986) show that surface tension—a singular perturbation—selects a discrete set of finger solutions (only one of which is linearly stable), which all converge to the $`\lambda =1/2`$ ST-finger in the limit of zero surface tension.

It has been suggested Arneodo et al. (1989) that the $`\lambda =1/2`$ ST-finger solution also models the average profile both of the unstable (highly branched) Hele–Shaw fingering and of the DLA growth in a channel with the corresponding reflective boundary conditions. The profile was defined as a level set of the ensemble averaged mass density, and for the experimental Hele–Shaw profiles half the maximum level was used. For DLA growth it was later shown Arneodo et al. (1996) that the level set at 0.5 maximum matches the width of $`\lambda =0.56`$, whilst the best match to the $`\lambda =1/2`$ profile came from the level set at 0.6 maximum. Outside that range the authors of Ref. Arneodo et al. (1996) concluded they could match level sets only to finger widths but not to the full shape of any ST-finger.

Here we use a different kind of finger averaging, which does not have any fitting parameter (eg. height of level set), as follows. In the conformal map method, we directly average the map. That is, we choose a set of points on the unit circle in the mathematical plane, and repeatedly map to the physical plane. The position of the image points averaged over the different maps, that is, over the different clusters that we have generated, is a reasonable alternative to the ensemble average of Arneodo et al. (1989, 1996).

In Fig. 2(a) we show the average conformal map generated this way. We see that the average map for DLA, $`\eta =1`$, does not correspond to the ST result for $`\lambda =1/2`$, but we get a good match to it for the DBM growth at $`\eta =1.2`$. This is an interesting result, especially in the context of recently proposed equivalences between DBM models with generalized local spatial cutoff. In that framework Ball and Somfai (2002, 2003) a higly ramified viscous finger with simple surface tension cutoff corresponds to standard (fixed size cutoff) DBM with $`\eta 1.2`$. Here we observe that the *non-branching* ST-finger solution is very similar to the *conformal average* of DBM clusters of the same $`\eta `$.

To use DLA grown with random walking particles in a channel, we need only construct the conformal map from the complex unit circle to the perimeter of each cluster, and take the average of these maps, as above. The conformal map is obtained numerically by the following method Somfai et al. (1999). We send $`M`$ probe particles to the frozen cluster, record their impact position, and discard them. These points correspond to $`M`$ uniformly distributed points on the unit circle. The landing positions of the probe particles are labeled topologically as one encounters them when tracking the perimeter of the cluster. Finally the $`m`$-th point is assigned to the angle $`2\pi m/M`$ of the unit circle. This angle has an error of the order of $`M^{1/2}`$, which vanishes for large $`M`$.

We measured the average conformal map on channels with reflective boundary conditions, widths ranging from 10 to 2000 particle diameters. For each width, we grew $`10^5`$ short clusters (only about $`10w`$ long), and probed each with $`10^5`$ test particles. In addition, for a few selected widths we probed $`10^4`$ clusters with $`10^6`$ probes each.

The average map for a wide channel generated this way is shown in Fig. 2(b). The curvature of the tip is consistent with that of the $`\lambda =1/2`$ ST-finger (the measured curvature of the DLA profile for $`w=1000`$ or $`2000`$ corresponds to $`\lambda =0.51\pm 0.03`$). The asymptotic width, however, is larger. The average map significantly differs also from the ST-finger of matching asymptotic width.

The average conformal map, rescaled onto a unit wide channel, shows strong dependence on the channel width. This is shown on Fig. 3 as a function of reduced $`x`$ position relative to the tip, $`\xi =(xx_{\text{tip}})/w`$. Details of the tip and tail regions show that these are clearly not consistent with a common asymptotic ST-finger shape.

We performed a finite size scaling on the $`w`$ dependence of the finger width. The filling ratio of the finger $`w_{\text{finger}}/w=[y_+(x)y_{}(x)]/w`$ was measured at selected $`\xi `$ values: $`\xi =0.5`$, $`1`$, and $`1.5`$, and is plotted on Fig. 4. The finite size scaling exponent was found to be $`\nu =0.33`$, same as for other quantities Somfai et al. (1999); Ball et al. (2002). The most interesting is the second extrapolation: $`\xi \mathrm{}`$. We only have 3 points for this, so it is only reasonable to give a lower bound: $`w_{\text{finger}}/w0.62`$.

In summary, using large scale random walker based simulations we have shown that the fractal dimension of DLA in a channel—with either periodic or reflective boundary conditions—is the same as in radial geometry. This is a great simplification compared to earlier claims of boundary condition (geometry) dependent fractal dimension. Second, using both iterated conformal maps and random walker based simulations, we measured the average profile of the clusters, defined by the average conformal map, and compared them to ST-finger solutions of the corresponding continuum problem. The averaged DLA profile is reminiscent but distinct from the ST-fingers, while the average profile of DBM clusters with $`\eta =1.2`$ are rather similar to the ST-finger with $`\lambda =1/2`$.

###### Acknowledgements.

We are indebted to Dave Kessler for the suggestion of measuring the average conformal map. This research has been supported by the EC under Contract No. HPMF-CT-2000-00800. The computing facilities were provided by the Centre for Scientific Computing of the University of Warwick, with support from the JREI.


text:
# Riemann zeros, prime numbers and fractal potentials

## I Introduction

In the early work of Wu et. al wu1  a general numerical algorithm was developed which could reconstruct a one-dimensional (1D) local potential from an essentially arbitrary set of levels (i.e., bound-states). The basic idea behind the method (see Sec. II below) is to functionally minimize the potential so that it reproduces the prescribed energy eigenvalues to numerical accuracy. This is, in principle, a solution of the 1D quantum inverse problem for a confined system if one could fit the infinitely many levels to a local potential. Of course, this is beyond any numerical algorithm’s ability, and one has to settle for fitting only the first $`N`$ eigenvalues. Since the algorithm requires a direct solution of the Schrödinger equation (SE) for every eigenvalue, it is practically limited to fitting roughly $`N𝒪(10^3)`$ eigenvalues. In an interesting application of this algorithm, Wu and Sprung wu2  chose as their eigenvalue spectrum the complex zeros (it is implicit throughout this paper that zeros refer to the complex zeros of the zeta function) of the Riemann zeta function, which is well known in number theory edwards . By finding a local 1D potential which exactly reproduced the first $`N=500`$ “eigenvalues”, they had effectively found a 1D quantum Hamiltonian for the first 500 zeros of the Riemann zeta function. This is of course interesting in its own right because such Hamiltonians have long been conjectured to hold the key to proving the celebrated Riemann hypothesis edwards . Their reconstructed potential was found to be fractal, and the fractal dimension was estimated to be $`D=1.5`$. This result suggests the intriguing possibility that an integrable system with time-reversal symmetry may generate the Riemann zeros provided the potential is fractal. It is worth commenting that the precise meaning one can attach to such 1D potentials is not clear. In particular, taking a chaotic system (i.e., the Riemann zeros) and forcing it to be integrable will lead to a potential that is dependent on the dimension of the Hilbert space (i.e., number of levels fitted).

Motivated by the work mentioned above, Ramani et. alramani1  also investigated the reconstruction of a 1D system’s quantum Hamiltonian from its eigenvalues. Their approach to solving the inverse problem was to use a “dressing transformation” based on techniques developed in the solution of the nonlinear Korteweg-de Vries equation. Their method, being much more efficient numerically than the direct SE approach resulted in many more levels being fitted for the same numerical effort. Unfortunately, rather than directly examining the Riemann zeros, as Wu and Sprung had previously done, they investigated the fractal properties of potentials generated from energy levels whose nearest-neighbour spacing distribution (NNSD) were the Gaussian Unitary Ensemble (GUE) stockmann . Although the NNSD of the Riemann zeros are also known to belong to the GUE, Ramani et. al could only make general comments about 1D potentials generated from the GUE statistics, but nothing specific about the Riemann zeros. Nevertheless, their work suggests, using $`𝒪(10^4)`$ levels, a fractal dimension of $`D=2`$ as $`N\mathrm{}`$ for GUE reconstructed potentials. By extension, they conjectured that their more accurate estimate of the fractal dimension leads to a value of $`D=2`$ for the Riemann potential.

In a reply to this work, Wu and Sprung wu3  have pointed out that using a different spectrum (i.e., a set of energy levels whose NNSD is the GUE rather than the actual complex Riemann zeros) does not necessarily imply that the fractal dimension of the reconstructed potential should be the same. One of the key reasons behind this statement is the lack of long-range correlations (i.e., the so-called spectral rigidity $`\mathrm{\Delta }_3`$) in the generic GUE spectrum used in Ref. ramani1 . Indeed, it was argued that if the fitting spectrum does not contain proper long-range correlations, the resulting potential may appear to be more random than the Riemann zeros potential. This, of course, can result in two different estimates for the fractal dimension even if a large number of levels have been fitted. Another important concern raised in their reply is that the potential obtained by inversion of the levels need not be unique. In other words, even if Ramani et. al had used the Riemann zeros as their fitting spectrum, they may not have obtained the same potential; that is, the fractal dimension may depend on the inversion technique applied.

In light of the above, we believe that there are several questions which have been left unanswered: (i) For the same set of energy levels, are the potentials reconstructed from different inversion procedures unique ? (ii) Does the fractal dimension depend on the inversion technique ? (iii) What is the fractal dimension for the Riemann zeros potential ? (iv) Is the fractal dimension a generic quantity for reconstructed 1D potentials or does it depend on the NNSD and $`\mathrm{\Delta }_3`$ statistics of the energy level spectrum ?

In what follows, we will attempt to provide answers to these questions. We will begin by presenting a brief overview of the numerical inversion techniques used in Refs. wu2 ; ramani1 , and then apply them to reconstruct the potentials of the Riemann zeros and prime numbers. Our decision to focus on these number sequences is three-fold. First, the contradictory results $`D=1.5`$ and $`D=2`$ for the fractal dimension of the Riemann potential remains unresolved. Secondly, the Riemann zeros are known to be dual to the prime numbers in the sense that all of the primes are encoded through the complex Riemann zeros via Riemann’s formula edwards ; bhaduri . One may then speculate as to whether any additional relationships between the Riemann zeros and prime numbers can be obtained from the reconstructed fractal potentials. Finally, since the NNSD of the Riemann zeros belong to the GUE class, and the prime numbers obey a more poisson-like distribution, it is interesting to investigate the dependence of the fractal dimension on the level statistics. The reconstructed potentials obtained from the two inversion methods will then be used to address the questions of uniqueness and fractal dimension.

## II Reconstructing the potential from a set of levels

In this section, we will give a brief overview of the two inversion techniques used in this paper to reconstruct the local 1D potentials. After this presentation, we will compare the reconstructed potentials for the Riemann zeros and prime number sequence obtained from both methods.

### II.0.1 Direct SE approach

Perhaps the most intuitive approach to solving the inverse problem is to directly invoke the SE. One initially starts from some arbitrary potential $`V(x)`$ (although for reasons of convergence, it is wiser to start from the semiclassical potential as described in Ref. wu2 ) and solve the SE to obtain a set of eigenvalues $`\epsilon _n`$. One then defines a “cost function”

$$F=\underset{n}{}(\epsilon _nE_n)^2,$$
(1)
where the $`E_n`$ are the exact eigenvalues one wishes to reproduce. The potential $`V(x)`$ is now adjusted in such a way as to minimize the cost function. Mathematically, this results in the following functional equation:

$$\frac{\delta F}{\delta V(x)}=2\underset{n}{}(\epsilon _nE_n)\varphi _n^2(x)=0,$$
(2)
where $`\varphi _n(x)`$ is the $`n`$-th normalized wavefunction. The functional minimization can be performed using a conjugate gradient method and the Numerov technique can be used to solve the discrete 1D SE. We have used this approach to successfully reproduce all of the numerical data for the Riemann zeros reported in Ref. wu2 . Unfortunately, even on modern PC/Workstations, this direct approach involves a large computational undertaking. We have found it impractical to go beyond $`N=2000`$ levels when using this technique.

### II.0.2 Method of the dressing transformation

The basis for the dressing transformation can be found in the soliton theory of nonlinear wave equations. The details of this approach are presented in Ref. ramani1 , and we give here only the essential equations required for its implementation. The starting point is a given set of levels $`ϵ_0,ϵ_1,\mathrm{},ϵ_{N1},ϵ_N`$ which have been shifted so that one has exactly $`N1`$ negative eigenvalues with the last eigenvalue satisfying $`ϵ_N=0`$. In this way, the potential is constructed from the “top down”. This does not impose any restrictions other than the requirement that the number of levels to be fitted is fixed in advance. An initial potential $`V(x)=0`$ is then used as input and the first-order differential equation (diffeq)

$$f^{}f^2+V(x)=ϵ$$
(3)
is evaluated numerically with $`f(0)=0`$ and $`ϵ=ϵ_{N1}`$. Once the function $`f`$ has been found at all points, a new potential $`W(x)`$ is constructed from the diffeq:

$$W(x)=2ϵ+2f^2V(x).$$
(4)
This new potential is then substituted back into Eq. (3) (i.e., $`W(x)V(x)`$) with the eigenvalue now being $`ϵ=ϵ_{N2}`$. The choice $`f(0)=0`$ is enforced throughout the calculation and ensures that the final potential $`V(x)`$ is even. This iterative procedure is continued until all of the eigenvalues have been exhausted. The final $`W(x)`$ is the desired reconstructed potential which reproduces exactly all of the $`N`$ eigenvalues. The integration of Eq. (3) can be performed using any standard method, but we have used a fourth-order Runge-Kutta (RK) method with step size $`h=1\times 10^5`$ to ensure high accuracy. Nevertheless, it should be clear that the dressing transformation is numerically a far superior method since it does not require a large parameter search in the construction of the potential. Consequently, it is easy to fit orders of magnitude more levels than in the direct SE approach.

### II.0.3 Uniqueness of the reconstructed potentials

Having provided the basic background for numerically solving the inverse problem, we are now in a position to address the question of uniqueness. To this end, we first reconstruct the potential for the Riemann zeros with $`N=500`$. This number was chosen because it represents the maximum number of Riemann zeros fitted in Ref. wu2 . Figure 1 displays the potential obtained from the direct SE approach (solid curve) and the dressing transformation (dashed). In the inset to the figure, we zoom-in on the the region $`x[0,5]`$ to illustrate the remarkable agreement between the two methods, which continues until one approaches energies near the last fitted level $`(x35)`$. Obviously above this threshold the potentials will not agree since the SE approach uses the semiclassical profile as its zeroth-order potential whereas the dressing transformation takes $`V(x)=0`$. However, as the number of levels is increased, this threshold is moved to higher and higher energies, so that as $`N\mathrm{}`$, the two potentials should agree over all space. The equivalence between the two inversion techniques has been checked for up to $`N=2000`$ (which computationally speaking, represents the upper limit of the SE approach), but we see absolutely no reason why this equivalence should not be true for arbitrarily large $`N`$. Nevertheless, to ensure that this agreement is not in some way fortuitous, we have also reconstructed the potential for the prime number sequence with $`N=500`$. The results of the two approaches for the primes are shown in Fig. 2. Again, the potentials are indistinguishable except for regions near the last fitted level.

It also proves instructive to compare in more detail the reconstructed potentials for the Riemann zeros and prime numbers. To facilitate this comparison, we present in Fig. 3 details of the potentials for $`N=2000,5000,10000`$ levels around the minimum (one could look at any other region on the curves and obtain similar graphs). Figure 3 (a)-(c) corresponds to the Riemann zeros whereas Fig. 3 (d)-(f) corresponds to the prime potential. Several interesting observations can be made at this point. First, we note that for a given number of fitted levels, the Riemann zeros potential looks less “noisy” than the prime potential. In fact, by $`N=10000`$, the prime potential takes on the characteristics of (roughly speaking) white noise whereas the Riemann zero potential appears to have far more structure. In particular, the Riemann potential contains large jumps at very specific spatial positions which appear to be quite robust (i.e., the positions of the jumps do not change as $`N`$ increases). In contrast, any local structure present in the prime potential is completely washed out by $`N=10000`$. This suggests that the prime potential may have a larger fractal dimension (i.e., close to $`D=2`$ of white-noise) than the Riemann zeros. Indeed, these observations lead us to consider more closely the fractal potentials generated by Ramani et. al ramani1 .

As we mentioned earlier, the potentials in Ref. ramani1  were obtained from a generic GUE distribution without accounting for the $`\mathrm{\Delta }_3`$ statistics. A comparison of our Fig. 3 (a)-(c) with their Fig. 1 (a)-(c) clearly illustrates that a generic GUE distribution cannot capture, even qualitatively, the details of the Riemann zeros potential. We are therefore inclined to agree with the assertion in Ref. wu3  that long-range correlations in the level spectrum play an important role in determining the fractal properties of the potential. This is further emphasized by noting that for large $`N`$, the prime potential \[see Fig. 3 (f)\] looks very similar to the large $`N`$ GUE fractal potential in Ref. ramani1 , which was conjectured to have $`D=2`$ as $`N\mathrm{}`$. It is perhaps then not so ambitious to suggest that the prime potential also has fractal dimension $`D=2`$. We cannot however make such claims for the Riemann zeros, which evidently require a more careful quantitative analysis.

## III Fractal dimension

We have seen in the previous section that for the same number of fitted levels, the oscillatory behaviour of the Riemann zeros and prime potentials are quite different. What significance then (if any) do these rapid oscillations have on the fractal dimension ? Before we answer that question, it is perhaps useful to give a short account of the method used by the authors in Refs. wu2 ; ramani1  to estimate the fractal dimension of their reconstructed potentials.

The standard procedure for numerically estimating the fractal dimension of curves such as the ones we have generated is to use the box-counting technique feder . Simply put, the box-counting technique involves choosing a region on the curve (over which the character of the oscillations does not change) and normalizing the axes so that the region is a square of side length equal to unity. The square region is then divided into $`n^2`$ cells of side length $`\mathrm{\Delta }x=1/n`$ and the number of cells $`P`$ that contain a portion of the curve is counted. A plot of $`\mathrm{log}(P\mathrm{\Delta }x)`$ versus. $`\mathrm{log}(\mathrm{\Delta }x)`$ is then constructed and the region over which the resulting graph has a constant slope $`1D`$ determines the box-dimension $`D`$; the box-dimension is usually simply referred to as the fractal dimension. Although it is not the only method for estimating the fractal dimension of a curve, we will make use of the box-counting technique to allow for a sensible comparison with the results in Refs. wu2 ; ramani1

In Fig. 4, we display the results of the box-counting method as applied to the Riemann potential for $`N=500,5000,10000`$. In keeping with Ref. wu2 , we have focused on the range $`x[0,10]`$ for the analysis. What is clear from the figure is that the slope, $`m`$, of the linear region in the graph is essentially unchanged for $`N>500`$. In fact, the slopes for $`N=5000`$ and $`N=10000`$ are virtually identical. Given that the fractal dimension is given by $`D=1+|m|`$, we conclude that to two significant figures, $`D=1.5`$ for the the Riemann zeros potential. This result is consistent with the findings of Wu and Sprung wu2  who likewise obtained $`D=1.5`$. We see no evidence at all for $`D=2`$ as was suggested by the work of Ramani et. al ramani1 .

Figure 5 shows the results of the box-counting method for the prime number potential. In contrast to the Riemann zeros, there is a noticeable increase in the slope of the curve (i.e., larger fractal dimension) as one increases the number of fitted levels from $`N=500`$ to $`N=5000`$. However, for $`10000<N<40000`$ (not shown), the slope is found to be practically unchanged. We are therefore led to conclude that the fractal dimension of the primes is $`D=1.8`$. We are, however, open to the possibility that as $`N\mathrm{}`$ the fractal dimension of the prime potential approaches $`D=2`$ based on the qualitative behaviour of the potential seen in Fig. 3 (d)-(f). Indeed, a similar slow convergence of the fractal dimension was also observed in Ref. ramani1  for the GUE statistics. However, numerically investigating this convergence would require possibly millions of levels, which is well beyond the capabilities of our current computational facilities.

## IV Summary and Discussion

We have numerically examined the reconstruction of a 1D local potential from a set of energy levels using two distinct inversion procedures. Applying these methods to the Riemann zeros and the prime number sequence we find that both inversion techniques produce the same quantum potential when applied to the same set of energy levels. This provides evidence that the inversion problem in 1D results in a unique potential. Thus, the answer to the first question raised in the introduction is affirmative. Based on this finding, it follows that the fractal dimension cannot depend on the method of inversion. We have also examined the issue of the fractal dimension of the Riemann zeros for up to $`N=40000`$ levels and estimate that $`D=1.5`$. This result is in complete agreement with Wu and Sprung’s wu2  earlier estimate (albeit using only $`N=500`$ levels), and illustrates that the potentials investigated by Ramani et. al ramani1  have little bearing on the Riemann zeros potential. In particular, we demonstrated that it is insufficient to simply account for the NNSD of the energy levels to capture the local details responsible for the fractal properties of the reconstructed potential; the Riemann zeros display local structure that is not reproduced by a generic GUE level spectrum. Therefore, long-range correlations appear to play a pivotal role in determining the fractal dimension of the potential.

The prime number potential has also been investigated, and found to have fractal dimension $`D=1.8`$. This result is somewhat surprising since the NNSD of the primes is almost poisson-like, whereas the Riemann zeros obey the (chaotic) GUE statistics. Unlike the Riemann potential, our calculation of the prime potential’s fractal dimension suggests that there may be a very slow convergence (as $`N\mathrm{}`$) of the fractal dimension to $`D=2`$. This possibility has also been noted in Ref. ramani1  where the fractal properties of white noise-like potentials corresponding to the GUE statistics were studied.

###### Acknowledgements.

We would like to thank Rajat Bhaduri and Donald Sprung for very fruitful discussions. BVZ acknowledges Prof. Dr. R. K. Bhaduri for financial support through the NSERC of Canada. DAWH acknowledges support from the Marsden Fund of the Royal Society of New Zealand and a University of Otago Research Grant.


text:
# 1 Introduction

## 1 Introduction

The correlated random walk, whose diffusive analog is described by hyperbolic diffusion equation, is possibly the simplest mathematical model allowing one to incorporate a form of momentum in addition to random or diffusive motion. The correlated random walk differs from the ordinary random walk in that the probabilistic element used at each step is the probability of continuing to move in a given direction rather than the probability of moving in a given direction independent of the direction of the immediately preceding step. Thus the process remains a Markov process.

The solution to a variety of forms of the hyperbolic diffusion equation has been given by Goldstein .

There is considerable literature  on physical processes in the field of thermophysics leading to a mathematical formulation in terms of a hyperbolic heat conduction equation. The physics behind such application implies that the signal propagation speed is finite rather than infinite as in the case of ordinary diffusion.

Until recently the thermal phenomena studied with the help of hyperbolic heat transfer equation do not consider the quantum aspects of heat transfer . All the results are obtained in classical approximations.

In this paper the quantum limit of heat transfer will be considered and quantum heat transport equation will be developed. The study was motivated by the advent of the ultra – short duration laser pulses, particularly those in attosecond ($`10^{18}`$ and zeptosecond  domains. Considering that time in which electron orbits the atomic nucleus is of the order of femtoseconds — with attosecond pulses it is possible penetrate the quantum structure of the path of electrons.

## 2 Overview of research

In present day materials sciences laser technology is a fundamental tool in such processes as surface deformation, materials hardening, welding, etc. Recently it has been shown possible to study accompanying phenomena on an atomic level. It is well recognized that on deposition of the laser pulse, electromagnetic energy is absorbed by a photon – electron interaction, the magnitude of which is described by the fine structure constant $`\alpha =1/137`$. The response of material to the thermal perturbation is governed by the relaxation time $`\tau `$. These two parameters fully describe the dynamics of energy transfer of the ultra – short (femtosecond and shorter) laser pulses on a nanometer scale. Conditions in the processes induced by ultra – short laser pulses are in stark contrast with those of longer laser pulse duration, the distinction being that at longer pulse duration ($``$ 1 ps) excited particles and their surroundings have had sufficient time to approach thermal equilibrium. On the other hand at temporal resolution $``$ 1 fs it is possible to resolve the dynamics of nonstationary transport phenomena.

For the transport phenomena of the nonrelativistic electrons there is a stochastic process that is ultimately connected with its propagation, namely the Brownian motion. As is well known Nelson  succeed in deriving the Schrödinger equation from the assumption that quantum particles follow continous trajectories in a chaotic background. The derivation of the usual Schrödinger equation follows only if the diffusion coefficient D associated with quantum Brownian motion takes the value $`D=\mathrm{}/2m`$ as assumed by Nelson.

In this paper the diffusion process of the quantum particles in context of the thermal energy will be studied. In the paper  the heat transport in a thin metal film (Au) was investigated with the help of the hyperbolic heat conduction equation (HHC). It was shown that when the memory of the hot electron gas is taken into account then the HHC is the dominant equation for heat transfer. The hyperbolic heat conduction equation for heat transfer in an electron gas can be written in the form :

$$\frac{1}{c^2}\frac{^2T}{t^2}+\frac{1}{c^2\tau }\frac{T}{t}=\frac{\alpha ^2}{3}^2T.$$
(1)
In equation (1) $`c`$ is the light velocity in vacuum, $`\alpha =e^2\mathrm{}^1c^1=1/137`$ is the fine structure constant and $`\tau `$ is the relaxation time.

In the present paper the quantum limit of the equation (1) will be obtained. When a high energy laser pulse hits the thin metal film it is possible that very short relaxation time are obtained. Since the relaxation time $`\tau `$ strongly depends on the temperature the relaxation times are of the order of 1 fs for $`T>10^3`$ K. In this case the mean free path of the hot electron is the same order as the thermal de Broglie wavelength. The temperature region for which mean free path $`\lambda `$ is of the order of $`\lambda \text{ -}_B`$ (de Broglie wavelength) we will define as the quantum limit of the heat transport and the master equation — the quantum heat transport equation (QHT).

## 3 Quantum limit of the heat transport

The classical notion of the field was born from a description of material media: the propagation of certain perturbations in such a medium led to the motion of a wave. A field, or a wave, gave initially, therefore, a picture of the collective motions of the medium and their propagation. For example the fluctuations of the electronic density in a solid state plasma are described by “plasma” waves. The classical formalism of the theory of fields is the same whether it is applied to “true” fields, which are thought of as existing in themselves, even in the absence of any material substratum, such as the electromagnetic field, or to a “phenomenological field” produced by the collective motions of a propagating medium. The concept of a frequency, for example, plays the same fundamental role in the analysis of these fields. The quantum synthesis of the concepts of frequency and energy is thus universal, i. e. they do not depend upon the specific nature of the phenomenon being considered and characterize the quantum theory in all its generality. The Planck – Einstein relation

$$E=\mathrm{}\omega $$

holds equally well, therefore, in quantum description of “phenomenological” fields. In other words, one must have a discretization of the energy into packets or quanta, $`E=\mathrm{}\omega `$, for any waves: acoustic or thermal, with a vibration frequency $`\omega `$, exactly as for an electromagnetic wave.

In  we developed the new hyperbolic heat transport equation which generalizes the Fourier heat transport equation for the rapid thermal processes. The hyperbolic heat conduction equation for the fermionic system can be written in the form:

$$\frac{1}{\left(\frac{1}{3}v_F^2\right)}\frac{^2T}{t^2}+\frac{1}{\tau \left(\frac{1}{3}v_F^2\right)}\frac{T}{t}=^2T,$$
(2)
where $`T`$ — denotes the temperature, $`\tau `$ — the relaxation time for the thermal disturbance of the fermionic system and $`v_F`$ is the Fermi velocity.

In the subsequent we develop the new formulation of the HHC considering the details of the fermionic system.

To start with we recapitulate the well known results from the simplest fermionic model of the matter: the Fermi gas of electrons in metals.

For the electron gas in metals the Fermi energy has the form :

$$E_F^e=(3\pi ^2)^{2/3}\frac{n^{2/3}\mathrm{}^2}{2m_e},$$
(3)
where $`n`$ — density and $`m_e`$ — electron mass. Considering that

$$n^{1/3}a_B\frac{\mathrm{}^2}{m_ee^2}$$
(4)
and $`a_B=`$Bohr radius, one obtains

$$E_F^e\frac{n^{2/3}\mathrm{}^2}{m_e}\frac{\mathrm{}^2}{ma^2}\alpha ^2m_ec^2,$$
(5)
where $`c=`$light velocity and $`\alpha =1/137`$ is the fine structure constant. For the Fermi momentum $`p_F`$, we have

$$p_F^e\frac{\mathrm{}}{a_B}\alpha m_ec$$
(6)
and for Fermi velocity $`v_F`$,

$$v_F^e\frac{p_F}{m_e}\alpha c.$$
(7)
Considering formula (7), equation (2) can be written as

$$\frac{1}{c^2}\frac{^2T}{t^2}+\frac{1}{c^2\tau }\frac{T}{t}=\frac{\alpha ^2}{3}^2T.$$
(8)
As it is seen from (8) the HHC equation is the relativistic equation as it takes into account the finite velocity of light. In order to derive the Fourier law from equation (8) we are forced to break the special theory of relativity and put in equation (8) $`c\mathrm{}`$, $`\tau 0`$. In addition it was demonstrated from HHC in a natural way, that in electron gas the heat propagation velocity $`v_hv_F`$ in the accordance with the results of the laser pump probe experiments.

In the following the procedure for the discretization of temperature $`T(\stackrel{}{r},t)`$ in hot fermion gas will be developed. First of all we introduce the reduced de Broglie wavelength

$$\lambda \text{ -}_B^e=\frac{\mathrm{}}{m_ev_h^e}v_h^e=\frac{1}{\sqrt{3}}\alpha c$$
(9)
and mean free path $`\lambda ^e`$,

$$\lambda ^e=v_h^e\tau ^e.$$
(10)
Considering formulae (9), (10) we obtain HHC for electron

$$\frac{\lambda \text{ -}_B^e}{v_h^e}\frac{^2T^e}{t^2}+\frac{\lambda \text{ -}_B^e}{\lambda ^e}\frac{T^e}{t}=\frac{\mathrm{}}{m_e}^2T^e.$$
(11)
Equation (11) is the hyperbolic partial differential equation which is the master equation for heat propagation in Fermi electron gas. In the following we will study the quantum limit of heat transport in the fermionic systems. We define the quantum heat transport limit as follows:

$$\lambda ^e=\lambda \text{ -}_B^e.$$
(12)
In that case equation (11) has the form:

$$\tau ^e\frac{^2T^e}{t^2}+\frac{T^e}{t}=\frac{\mathrm{}}{m_e}^2T^e,$$
(13)
where

$$\tau ^e=\frac{\mathrm{}}{m_e(v_h^e)^2}.$$
(14)
Equation (13) defines the master equation for quantum heat transport (QHT). Having the relaxation time $`\tau ^e`$ one can define the “pulsations” $`\omega _h^e`$:

$$\omega _h^e=(\tau ^e)^1$$
(15)
or

$$\omega _h^e=\frac{m_e(v_h^e)^2}{\mathrm{}}$$

, i.e.

$$\omega _h^e\mathrm{}=m_e(v_h^e)^2=\frac{m_e\alpha ^2}{3}c^2.$$
(16)
Formula (16) defines the Planck-Einstein relation for heat quanta $`E_h^e`$

$$E_h^e=\omega _h^e\mathrm{}=m_e(v_h^e)^2.$$
(17)
The heat quantum with energy $`E_h=\mathrm{}\omega `$ can be named as the heaton in complete analogy to the phonon, magnon, roton and etc. For $`\tau ^e0`$ equations (13), (16) are the quantum Fourier transport equations (QFT) with quantum diffusion coefficients $`D^e`$:

$$\frac{T^e}{t}=D^e^2T^e,D^e=\frac{\mathrm{}}{m_e}.$$
(18)
For finite $`\tau ^e`$ for $`\mathrm{\Delta }t<\tau ^e`$ equation (13) can be written as follows:

$$\frac{1}{(v_h^e)^2}\frac{^2T^e}{t^2}=^2T^e.$$
(19)
Equation (19) is the wave equation for quantum heat transport (QHT).

It is interesting to observe that equation (19) is the wave equation which describes the propagation of the de Broglie thermal wave. For the wave length $`\lambda \text{ -}`$ of the wave which is the solution of equation (19) one obtains:

$$\lambda \text{ -}=\frac{v_h}{\omega }=\frac{v_h\mathrm{}}{\omega \mathrm{}}=\frac{\mathrm{}}{p_e},p_e=m_ev_h,$$
(20)
i.e. $`\lambda \text{ -}`$ is equal the reduced de Broglie wave length (9).

On the other hand the quantum Fourier equation (18) resembles the free Schrödinger equation. The replacement $`tit/2`$ turns the quantum diffusion equation into the Schrödinger equation. Both are parabolic and require the same boundary and initial conditions in order to be “well posed”.

The heaton energies for electron gas can be calculated from formula (17). For electron gas we obtain from formulae (9), (17) for $`m_e=0.51\mathrm{MeV}/c^2`$, $`v_h=(1/\sqrt{3})\alpha c`$

$$E_h^e=9\mathrm{eV},$$
(21)
which is of the order of the Rydberg energy. The numerical values of the relaxation time $`\tau `$ and the quantum diffusion coefficient $`D`$ can be calculated from formulae (14) and (18).

$$\tau =10^{17}\mathrm{s},D=\mathrm{1.2\hspace{0.17em}10}^4\mathrm{m}^2\mathrm{s}^1.$$

Concluding, one can say that for the temporal resolution $`\mathrm{\Delta }t`$ of the order or shorter of the relaxation time $`\tau 10^{17}\mathrm{s}`$ the heat transport phenomena are adequately described by the quantum heat transport equation (QHT), formula (13). It seems that the advent of lasers with attosecond laser pulses open quite new possibility for the study theses discrete thermal phenomena .


text:
# Study of the 𝜂𝜋 and 𝜂'𝜋 spectra and interpetation of possible exotic 𝐽^{𝑃𝐶}=1⁻⁺ mesons

## Abstract

We discuss a coupled channel analysis of the $`\eta \pi `$ and $`\eta ^{}\pi `$ systems produced in $`\pi ^{}p`$ interactions at 18 GeV/$`c`$. We show that known $`Q\overline{Q}`$ resonances, together with residual soft meson-meson rescattering, saturate the spectra including the exotic $`J^{PC}=1^+`$ channel. A possibility of a narrow exotic resonance at a mass near 1.6 GeV/$`c^2`$ cannot, however, be ruled out.

Gluons play a central role in the structure of matter as they mediate the strong interactions leading to quark and gluon confinement. The quantitative description of the dynamics of low energy gluons, however, is still unavailable and the direct manifestation of gluons in the current spectroscopy of hadrons, for the most part, seems to be elusive. Any theoretical approach is complicated by the strong character of the interactions of soft gluons in QCD. Furthermore, the phenomenology of gluonic excitations is challenging because gluons do not carry electromagnetic or weak charge to which probes can directly couple. The spectroscopy of hadrons with excited gluonic degrees of freedom remains as the only avenue to obtain information about the nature of the soft gluonic modes. In this context the spectroscopy of exotic mesons is of fundamental importance. Exotic mesons have quantum numbers that cannot be attributed to valence quarks alone. Lattice gauge simulations indicate the gluonic field confining quarks forms flux tubes and the quantum numbers of the excited flux tubes couple with those of the quarks leading to exotic quantum numbers. The spin ($`J`$), parity ($`P`$) and charge conjugation ($`C`$) quantum numbers of the lightest exotic isovector should be $`J^{PC}=1^+`$ and have mass around 1.9 GeV/$`c^2`$ l1 ; l2 ; l3 . It is possible, however, that corrections due to extrapolation of the lattice mass predictions to the physical light quark masses could introduce a 100 - 200 MeV/$`c^2`$ downward shift astt . Based on the large $`N_C`$ expansion it has been shown that exotic mesons ought to have hadronic decay widths comparable to the other mesonic resonances cohen  i.e. to be of the order of $`\mathrm{\Gamma }=`$100 - 200 MeV/$`c^2`$. Lattice simulations also indicate that at least one exotic meson multiplet should exist below the threshold energy for string breaking and quark production colin .

Finding a $`J^{PC}=1^+`$ wave is not equivalent to establishing the existence of a QCD exotic. A partial wave can have many poles as a function of complex energy in the unphysical domain of the complex energy plane (Riemann sheets). With an incomplete knowledge of the underlying strong interaction dynamics it is impossible to unambiguously determine positions of distant poles. Similarly limited experimental data are unable to isolate and interpret such poles. Only poles relatively close to the real axis can potentially be unambiguously identified, however, even in this case theoretical input is needed to discriminate between a preexisting bound state and a dynamical resonance.

The current experimental situation appears to be puzzling. Recently three exotic $`J^{PC}=1^+`$ meson candidates have been reported. The $`\pi _1(1400)`$, which has an unusually light mass, $`M=1370\pm 16_{30}^{+50}\text{ MeV}/c^2`$, and a large width, $`\mathrm{\Gamma }=385\pm 40_{105}^{+65}\text{ MeV}/c^2`$, was reported by the E852 collaboration in $`\pi ^{}p\pi _1^{}(1400)p\pi ^{}\eta p`$ eta  at 18 GeV/$`c`$. This state was confirmed in the $`\eta \pi `$ channels in $`\overline{p}d`$ annihilations to $`\pi ^{}\pi ^0\eta p`$ by the Crystal Barrel collaboration cb . Two other $`\pi _1(1600)`$ exotics were also reported by the E852 collaboration, one with $`M=1597\pm 10_{10}^{+45}\text{ MeV}/c^2`$ and a large width of $`\mathrm{\Gamma }=340\pm 40\pm 50\text{ MeV}/c^2`$ observed in the $`\eta ^{}\pi ^{}`$ channel eprim , and the other with $`M=1593\pm 8_{47}^{+29}\text{ MeV}/c^2`$ and a normal width of $`\mathrm{\Gamma }=168\pm 20_{12}^{+150}\text{ MeV}/c^2`$ decaying into $`\rho ^0\pi ^{}`$ rho . The $`\eta \pi `$ and $`\eta ^{}\pi `$ channels are well suited for $`J^{PC}=1^+`$ exotic meson searches, since a neutral $`P`$-wave resonance has exotic quantum numbers.

It should be noted that the above resonance parameters come from a Breit-Wigner (BW) parameterization thus assuming a resonance interpretation of the data. But a BW parameterization is only applicable if amplitudes have singularities near the real axis, which in the case of a $`300400\text{ MeV}/c^2`$ wide resonances is not reasonable. Furthermore, in the original analysis of the E852 data, only selected $`J_z`$ states of the $`\eta \pi `$ and $`\eta ^{}\pi `$ systems were considered in the mass dependent analysis that gave the resonance parameters. If an exotic $`J=1`$ resonance is produced, however, all $`2J+1`$ spin components should be included in the analysis. Even though these spin projections can have different production characteristics they should all display similar resonance behavior as a function of the invariant mass of the decay products. In a recent analysis of the $`\eta \pi ^0`$ E852 data alex  it was shown that an exotic $`P`$-wave is indeed present, and very similar to the one found in the $`\eta \pi ^{}`$ mode, but when all information was used, including all spin components and production modes, no self-consistent set of BW parameters describing the observed $`P`$-wave could be found.

Within QCD, the quark structure of a typical narrow meson resonance, $`e.g`$ the $`\rho (770)`$ or the $`a_2(1320)`$, is not much different from a bound state. The quark wave functions are compact and the small width arises from coupling to a few open channels. In contrast, broad structures like the $`\pi _1(1400)`$ or the $`\pi _1(1600)`$, as observed in the $`\eta \pi `$ and $`\eta ^{}\pi `$ channels, respectively, are likely to be of a different origin. In this paper we present a new analysis of the E852 data and explain the dynamical origin of the $`P`$-wave enhancement seen in the $`\eta \pi `$ and $`\eta ^{}\pi `$ channels. In particular, we argue that these might be very similar to the $`\sigma `$ meson which is used to parameterize the low energy $`S`$-wave, isoscalar $`\pi \pi `$ spectrum. It is generaly accepted that the $`\sigma `$ meson does not originate from a quark bound state but is instead a dynamical resonance arising from residual, low energy $`\pi \pi `$ interactions julich ; oset .

The low energy interactions of the pseudoscalar meson nonet can be constrained using standard low energy expansion methods, e.g. effective range theory. Chiral symmetry provides a natural scale for such an expansion, $`f100\text{ MeV}/c`$ gasser . In particular, for a two-body potential acting between two-meson channels, $`i,j`$ in $`L`$-th partial wave one has,

$$V_{ij}^L(p,q)=p,i|V_L|q,j=(pq)^L\left[c_{0,ij}+c_{2,ij}\frac{p^2+q^2}{f^2}\mathrm{}\right]$$
(1)
with, $`p`$ and $`q`$ denoting the relative momenta between mesons. It is expected that the coefficients, $`c_i`$ are $`O(1)`$, and that these coefficients determine the low momentum expansion of the scattering amplitude. The real part of the scattering amplitude, however, arises from iteration of the potential leading to important non-analytical contributions from integration over low momentum components of the virtual states. These non-analytical contributions have, in many cases, been found to be sufficient to extend the range of applicability of the low momentum expansion to regions well above the momentum scale set by $`f`$, even up to $`1\text{ GeV}/c`$ oset ; tony ; tony1 . Furthermore, it is expected that higher terms in the expansion can be approximated by pole terms corresponding to genuine QCD bound states. The approach can also be verified a posteriori by studying how sensitive are the parameters of these bound sates and intensity of the coherent backgrounds introduced by $`V_L`$ to the details of the approximations. In particular, we have studied the $`N/D`$ approach with $`V_L`$ determining the numerator and the poles determining the CDD terms CDD . We also examined the Lippmann-Schwinger formalism with and without relativistic kinematics. Finally, in all cases we studied various methods for removing the large momentum contributions from loop integrals, e.g. with or without explicit subtraction of power divergences. Below we present results which we have found to be robust with respect the various implementations.

The distribution of events, $`N(M_X,\mathrm{\Omega },t)`$, at fixed energy in the reaction $`\pi ^{}pX^{}p`$ with $`X=\eta \pi ^{}`$ or $`\eta ^{}\pi ^{}`$, is proportional to $`|_{LM}a_{LM}(t,M_X)Y_{LM}(\mathrm{\Omega })|`$. The partial wave amplitudes, $`a_{LM}`$, determine the production strengths of a two-pseudoscalar meson system, with spin $`J=L`$, as a function of their invariant mass ($`M_X`$) and the momentum transferred squared, $`t`$, to the target. The peripheral nature of the high-$`s`$, low-$`t`$ production implies that the $`M_X`$ and $`t`$ dependence is expected to factorize. Verification of the assumption can be found in alex . The dependence on the solid angle, $`\mathrm{\Omega }`$, is defined with respect to the Gottfried-Jackson frame. For a mass dependent analysis we bin the data in $`M_X`$ and integrate over all $`t`$ (the limited statistics prevents us from considering more than two $`t`$-bins). The number of events, $`N(M_X,\mathrm{\Omega })`$, in a given mass bin is decomposed into moments defined as $`H_{LM}(M_X)=𝑑\mathrm{\Omega }N(M_X,\mathrm{\Omega })`$ which are corrected for detector acceptance. The partial wave amplitudes are calculated from two-meson interactions as described above and fitted to the moments. In addition, it is necessary to specify the production amplitudes, $`P_i^L`$, whose mass dependence can, in principle, be constrained by duality. Instead, we use the requirement that they should have a expansion similar to the potential,

$$P_i^L(q)=q,ip|V_L|\pi ^{}p=A_i^Lq^L\left(1+d_{2,i}\frac{q^2}{f^2}\mathrm{}\right)$$
(2)
with the relative strengths satisfying, $`A_i^L/A_i^L^{}=O(1)`$. In the mass range from threshold to $`M_X<2\text{ GeV}/c^2`$ we find that only the $`L=S,P,D`$ waves are needed to describe the data. The partial waves are further classified by the magnitude of the spin projection, $`M=0,1,\mathrm{}L`$ onto the beam axis and parity under reflection in the production plane. We find that only amplitudes with $`M<1`$ are relevant, leading to a set of seven partial waves, $`S,P_0,P_\pm ,D_0,D_\pm `$. Finally the positive and negative parity waves do not interfere and arise from $`t`$-channel natural and unnatural parity exchange, respectively.

For completeness, in addition to the $`\eta \pi `$ and $`\eta ^{}\pi `$ channels, we also include the isosvector $`K\overline{K}`$ channel. In the $`S`$-wave the $`K\overline{K}`$ attraction near threshold produces a cusp at $`M_X1\text{GeV}/c^2`$ which corresponds to the $`a_0(980)`$ meson oset . A QCD resonance is expected at $`M_X1.31.4\text{ GeV}/c^2`$ and could be associated with the $`a_0(1400)`$ meson. This is similar to the scalar-isoscalar channel where the low energy $`\pi \pi `$ attraction results in the initial growth of the phase shift from threshold to $`700\text{ MeV}/c^2`$ which is commonly attributed to an effective $`\sigma `$ meson. The opening of the attractive $`K\overline{K}`$ channel responsible for the $`f_0(980)`$ meson leads to the rapid phase increase at $`1\text{ GeV}/c^2`$. Finally further growth of the phase in the $`1.31.5\text{ GeV}/c^2`$ can be attributed to presence of QCD resonances. In the $`S`$-wave we thus include all three coupled two-meson channels and a single pole term with a mass expected around $`1.31.4\text{ GeV}/c^2`$ julich ; oset . The $`D`$ wave is dominated by the narrow $`a_2(1320)`$ meson which we introduce as a pole term. The soft interactions in the $`D`$-wave correspond to $`O(p^4/f^4)`$ in the chiral expansion and thus cannot be unambiguously disentangled from pole contributions, also assumed to be $`O(p^4/f^4)`$. This introduces a theoretical uncertainty in the coherent $`D`$-wave background. Thus the only constraint imposed on the $`D`$-wave is that coefficients are of $`O(1)`$. The ability to correctly describe the known spectrum in the $`S`$ and $`D`$-wave in an important check on both the theoretical method and data selection.

The $`P`$-wave is of the prime interest here. We have studied cases with and without a pole term. The soft interactions in the $`P`$-wave channel are less constrained than, for example, the $`S`$-wave since they require mapping between the physical, $`\eta ,\eta ^{}`$ and the $`SU(3)`$ flavor, $`\eta _0,\eta _8`$ states marco . In general it is predicted however that both $`\eta \pi `$ and $`\eta ^{}\pi `$ interactions are attractive with the former weaker by $`\mathrm{sin}^2(\theta )0.01`$ where $`\theta `$ is the $`\eta _0\eta _8`$ mixing angle. With these general assumptions we have found that the predictions for the spectrum of the $`P`$-wave is very insensitive to the details of the channel interactions.

The results of the combined analysis of the $`\eta \pi ^{}`$ and $`\eta ^{}\pi `$ channels are shown in figures 1 and 2. The $`\eta \pi `$ channel is dominated by the $`a_2(1320)`$ resonance. The contribution from the coherent $`D`$-wave background is small leading to a pole mass within $`10\%`$ of the resonance mass listed by the PDG. The $`S`$-wave production is very weak and does not constrain parameters of the $`S`$-wave pole, the $`a_0(1400)`$ not the dynamically generated, $`a_0(980)`$. The discrepancy between data and theory in the $`\eta \pi `$ intensity at $`M_X1\text{ GeV}/c^2`$ is within systematical errors of the $`S`$-wave parameters.

We did find, however, that the $`H_{21}`$ moment which is sensitive to the interference between the unnatural exchange, $`D_{}`$ wave and the $`S`$ wave, has nontrivial structure at $`M_X=1.4\text{ GeV}/c^2`$ which could be related to the the $`a_0(1400)`$ pole. The $`P`$ wave intensity in the $`\eta \pi `$ channel is weak, $`5\%`$ of the $`D`$ wave, and leads to a broad structure consistent with the absence of nearby resonance in this channel. As discussed in alex  the peak-like structure at $`M_X=1.3\text{GeV}/c^2`$ is most likely due to the dominant $`a_2(1320)`$ leaking into the $`P`$-wave due to imperfect acceptance corrections. The phase of the natural exchange $`P_+`$ wave is found to be very weakly mass dependent, $`|\delta \varphi _P|<10^0`$ over the mass range considered, and as shown in figure 2(d) consistent with the data. The variations in the $`\varphi _D\varphi _P`$ phase difference at $`M_X>1.6\text{ GeV}/c^2`$ is due entirely to the weak coherent $`D`$-wave background. Similarly the sharp variation of this phase difference at low mass can be fully accounted for by the $`D_+`$ coherent background. Even though the channel interaction in the $`D`$ wave are not well constrained it seems clear that both phase and magnitude of the $`P`$-intensity in the $`\eta \pi `$ channel can be well accounted for by $`\eta \pi `$ rescattering. In contrast, in the $`\eta ^{}\pi `$ channel the contribution from the $`P`$ wave is significant and consistent with the theoretical expectations discussed above. As shown in figure 1b, the $`P`$ is as large as the $`D`$-wave while the $`S`$-wave is small. Since the $`D`$ and $`P`$ waves are comparable there is no danger of one leaking into another and the partial wave analysis is quite unambiguous. It is seen from figure 2 that the soft channel interactions describe both the intensity and the phase motion of the $`P`$ wave very well.

We have found, however, that there is a structure near $`M_X=1.6\text{ GeV}/c^2`$ in several moments which could be accounted for by a pole term in the $`P`$-wave. Since the bulk of the $`P`$-wave intensity comes from soft re-scattering, the resonance associated with the pole term becomes narrow. If parameterized as a BW resonance it has $`\mathrm{\Gamma }200\text{MeV}/c^2`$ and is thus consistent with the narrow resonance observed in the $`\rho \pi `$ channel.

In conclusion, we find that the broad $`P`$ wave in the $`\eta \pi `$ and $`\eta ^{}\pi `$ spectra can be accounted for by low energy re-scattering effects. Comparing this to the $`\pi \pi `$ channel, we interpret this enhancement as the equivalent of the $`\sigma `$ meson, i.e. arising from correlated meson-meson state and not from a QCD bound state. The applicability of low momentum expansion used here is well justified since, after subtracting the threshold energies, the relevant momentum range is comparable to the $`\pi \pi `$ case and the two analyses are completely consistent. In the non-exotic $`P`$-wave $`\pi \pi `$ channel there are, in addition to residual meson-meson interactions, narrow resonances like the $`\rho `$ and $`\varphi `$ mesons that are believed to be true QCD bound states. We have found that the $`\eta ^{}\pi `$ data does not rule out a narrow exotic $`P`$-wave resonance which might be the $`\pi _1(1600)`$ that is found to decay into $`\rho \pi `$ but from no final conclusions can be drawn at this time.

This work was supported in part by the U. S. Department of Energy.