text:
# Percolation and localization in the random fuse model

## 1 Introduction

Understanding the scaling properties of fracture in disordered media represents an intriguing theoretical problem with important implications for practical applications breakdown . Experiments have shown that in several materials under different loading conditions, the fracture surface is rough and can be described by self-affine scaling man  with universal exponents bouch . Scaling is also observed in acoustic emission experiments, where the distribution of pulses decays as a power law over several decades. Experimental observations have been reported for several materials such as wood ciliberto , cellular glass strauven , concrete ae  and paper paper , but universality does not seem to hold.The experimental observation of power law behavior suggests an interpretation in terms of critical phenomena and scaling theories, but a complete theoretical explanation has not been found. The statistical properties of fracture in disordered media are often studied with lattice models, describing the medium as a discrete set of elastic bonds with random failure thresholds breakdown . These numerical simulations are used in estimating the roughness of the fracture surface, which is found to be self-affine hansen91b ; bat-98 ; rai-98 ; sep-00 , and the power law distribution of avalanche precursors hansen ; zrvs ; zvs ; gcalda ; alava . While the results agree qualitatively with experiments a quantitative comparison is not always satisfactory.

Apart from the comparison with experiments, an important theoretical issue is to understand the origin of the scaling behavior observed in the numerical simulations of lattice model. A very well studied model is the random fuse model (RFM), where a lattice of fuses with random threshold are subject to an increasing voltage deArcangelis85 ; breakdown ; zrvs ; kahng88 ; delaplace ; deArcangelis89 ; nukalajpamg . A resistor network represents a scalar analog of an elastic medium and is thus relatively simple to analyze, while retaining some important characteristic features of the problem.

Simulations of the RFM show that the type of behavior at macroscopic fracture is significantly influenced by the amount of disorder kahng88 . When the disorder is narrowly distributed, materials breakdown without significant precursors. As the disorder increases, substantial damage is accumulated prior to failure and the dynamics resembles percolation sornette . Indeed, in the limit of infinite disorder, the damage accumulation process can exactly be mapped onto a percolation problem guyon88 . It is still debated, however, if percolation is applicable to the case of strong but non-infinite disorder. Recently, some evidence has been provided in this direction hansen003  suggesting that the critical exponent, $`\nu `$, of the correlation length in the RFM with strong disorder is same as that of uncorrelated percolation (i.e. $`\nu =4/3`$), and that the fuse model is in the same universality class of percolation. Close to failure damage would then localize and the resulting crack roughness would ensue from a gradient percolation mechanism hansen003 .

In the present work, we propose a different interpretation of damage localization in the RFM: while in the initial stages damage accumulates as in a percolation process, the final crack nucleates abruptly due to current enhancement, yielding the observed localization profiles. As a consequence of this, the damage localization length follows closely the scaling of crack width. In addition, we test percolation scaling by simulating the RFM with triangular and diamond (square lattice with 45 degrees inclined bonds to the bus bars) lattice topologies deArcangelis85 ; breakdown  and different disorder distributions (uniform and power law). The numerical results allows us to exclude that $`\nu =4/3`$, but could still be compatible with another value of $`\nu `$, although devations from scaling are observed for large sizes.

Finally, we show that the cumulative damage probability distribution at failure can be collapsed for different lattice sizes and follows the normal distribution. This suggests that damage is the sum of uncorrelated variables and thus no long-range correlation are present as would be expected in critical phenomena. This fact, together with the abrupt localization occurring right at failure, suggests the validity of the first-order transition hypothesis discussed in Ref. zrvs .

The paper is organized as follows: in section II we define the model and in section III we discuss the role of disorder. The percolation hypothesis is tested in section IV, while section V is devoted to damage localization. The scaling of the number of broken bonds is discussed in section VI. Section VII presents the analysis of the failure probability distribution and section VIII summarizes briefly our conclusions.

## 2 Model

In the random thresholds fuse model, the lattice is initially fully intact with bonds having the same conductance, but the bond breaking thresholds, $`t`$, are randomly distributed based on a thresholds probability distribution, $`p(t)`$. The burning of a fuse occurs irreversibly, whenever the electrical current in the fuse exceeds the breaking threshold current value, $`t`$, of the fuse. Periodic boundary conditions are imposed in the horizontal direction to simulate an infinite system and a constant voltage difference, $`V`$, is applied between the top and the bottom of lattice system bus bars.

Numerically, a unit voltage difference, $`V=1`$, is set between the bus bars and the Kirchhoff equations are solved to determine the current flowing in each of the fuses. Subsequently, for each fuse $`j`$, the ratio between the current $`i_j`$ and the breaking threshold $`t_j`$ is evaluated, and the bond $`j_c`$ having the largest value, $`\text{max}_j\frac{i_j}{t_j}`$, is irreversibly removed (burnt). The current is redistributed instantaneously after a fuse is burnt implying that the current relaxation in the lattice system is much faster than the breaking of a fuse. Each time a fuse is burnt, it is necessary to re-calculate the current redistribution in the lattice to determine the subsequent breaking of a bond. The process of breaking of a bond, one at a time, is repeated until the lattice system falls apart. In this work, we consider a uniform probability distribution, which is constant between 0 and 1, and a power law distribution $`p(t)t^{1+\beta }`$, $`0t1`$ with $`\beta =\frac{1}{20}`$.

Numerical simulation of fracture using large fuse networks is often hampered due to the high computational cost associated with solving a new large set of linear equations every time a new lattice bond is broken. The authors have developed rank-1 sparse Cholesky factorization updating algorithm for simulating fracture using discrete lattice systems nukalajpamg . In comparison with the Fourier accelerated iterative schemes used for modeling lattice breakdown batrouni , this algorithm significantly reduced the computational time required for solving large lattice systems. Using this numerical algorithm, we were able to investigate damage evolution in larger lattice systems (e.g., $`L=1024`$), which to the authors knowledge, is so far the largest lattice system used in studying damage evolution using initially fully intact discrete lattice systems. However, in this paper, we consider results up to $`L=512`$ due to insufficient number of available sample configurations for $`L=1024`$, which will be considered in a future publication.

Using the algorithm presented in nukalajpamg , we have performed numerical simulations on two-dimensional triangular and diamond (square lattice inclined at 45 degrees between the bus bars) lattice networks. For many lattice system sizes, the number of sample configurations, $`N_{config}`$, used are excessively large to reduce the statistical error in the numerical results (see Tables 1 and 2). Each numerical simulation was performed on a single processor of Eagle (184 nodes with four 375 MHz Power3-II processors) supercomputer at the Oak Ridge National Laboratory. The statistically independent $`N_{config}`$ number of configurations were simulated simultaneously on number of processors available for computation.

## 3 The role of disorder

The disorder distribution has clearly an important effect on the fracture behavior and a classification has been proposed in terms of a scale-invariant spectrum hansen91 ; hansen001 , which for the uniform thresholds distribution $`(0t1)`$ in terms of intensive variables, $`\alpha _t`$ and $`f_t(\alpha _t)`$, is given by hansen91 ; hansen001

$$f_t(\alpha _t)=2\alpha _t\text{for}0\alpha _t2$$
(1)
where

$`\alpha _t`$ $`=`$ $`{\displaystyle \frac{\text{log}t}{\text{log}L}}`$ (2)
$`f_t(\alpha _t)`$ $`=`$ $`{\displaystyle \frac{\text{log}L^2tp(t)}{\text{log}L}}`$ (3)
Hence, for the uniform distribution between 0 and 1, the two control parameters $`\varphi _0`$, and $`\varphi _{\mathrm{}}`$ that characterize the thresholds distribution $`p(t)`$ close to zero and infinity, respectively, are given by $`\varphi _0=1`$, and $`\varphi _{\mathrm{}}=\mathrm{}`$ hansen91 ; hansen001 , where $`\varphi _{0/\mathrm{}}`$ are defined by

$$\varphi _{0/\mathrm{}}=\underset{t0/\mathrm{}}{lim}\left(\frac{\mathrm{log}(tp(t))}{\mathrm{log}(t)}\right)$$
(4)
For the power law thresholds distributions as the one we are using, the scale-invariant thresholds spectrum is given by

$$f_t(\alpha _t)=2\beta \alpha _t\text{where}0\alpha _t\frac{2}{\beta }$$
(5)
and the two control parameters are $`\varphi _0=\beta `$ and $`\varphi _{\mathrm{}}=\mathrm{}`$ hansen91 ; hansen001 . According to Ref. hansen91 ; hansen001 , based on the values of these two control parameters, $`\varphi _0`$ and $`\varphi _{\mathrm{}}`$, both the uniform thresholds distribution and power-law thresholds distributions (as long as $`\beta <2`$ hansen001 ) belong to the same scaling regime characterized by diffusive damage and localization (see regime B of Figure 18 in reference hansen001 ). According to this analysis, if the exponent $`\nu `$ were to be universal and is in the same universality class as that of uncorrelated percolation as conjectured in hansen003 , then we expect to find $`\nu =4/3`$ for uniform and power law threshold distributions.

## 4 Damage and percolation

In the case of strong disorder, in the initial stages of damage evolution, bond breaking events occur in an uncorrelated manner and thus resemble percolation. As the damage starts to accumulate, some degree of correlation can be expected because of the current enhancement present at the crack tips. A natural question to ask concerns the relevance of these correlations as failure is approached. If correlations are irrelevant one should observe percolation scaling up to failure, as in the case of infinite disorder. On the other hand, in the weak disorder case, the current enhancement at the crack tips is so strong that a spanning crack is nucleated soon after a few bonds (or even a single bond) are broken kahng88 . The interesting situation corresponds to the diffuse damage and localization regime, where a substantial amount of damage is accumulated prior to failure. In this regime one should test if damage follows percolation scaling up to failure.

Hansen and Schmittbuhl hansen003  have considered broad threshold distributions $`p(t)t^{1+\beta }`$, $`0t1`$ with two different $`\beta =\frac{1}{10}`$ and $`\beta =\frac{1}{3}`$ values, to see if percolation scaling was observed for $`\beta >0`$. Based on the similarities with percolation, they hansen003  suggested the following finite size scaling law for the fraction of broken bonds, given by

$`p_fp_c`$ $``$ $`L^{\frac{1}{\nu }}`$ (6)
In the Eq. (6), $`p_f`$ and $`p_c`$ represent the fracture thresholds in a lattice system size of $`L`$ and infinity, respectively. As the system size $`L\mathrm{}`$, the broken bonds at failure $`p_fp_c`$. The correlation critical exponent $`\nu `$ was found in Ref. hansen003  to be consistent with the percolation value $`\nu =4/3`$. An additional test is provided by the damage standard deviation at failure $`\mathrm{\Delta }_f`$ ramstad  which should scale as

$$\mathrm{\Delta }_fL^{\frac{1}{\nu }}$$
(7)
Here we test these scaling laws for a wider finite size range than Ref. hansen003 , wherein simulations with sizes up to $`L=60`$ are reported. The fraction of broken bonds for each of the lattice system sizes is obtained by dividing the number of broken bonds with the total number of bonds, $`N_{el}`$, present in the fully intact lattice system. For triangular lattice topology, $`N_{el}=(3L+1)(L+1)`$, and for diamond lattice topology, $`N_{el}=2L(L+1)`$. The lattice system sizes considered in this work are $`L=\{8,16,24,32,64,128,256,512\}`$. However, since corrections to the scaling laws are strongest for small lattice systems, in the following, we use lattice sizes $`L16`$ for obtaining the scaling exponents. Table 1 presents mean and standard deviations in the broken bond density (fraction of broken bonds) at the peak load and at failure for various lattice system sizes in both the triangular and diamond lattice systems for uniform thresholds distribution. In order to compare diamond and triangular lattice topologies, we find it more convenient to use $`N_{el}`$ rather than $`L`$ as a finite size parameter. This is because the two lattices have a different dependence of the real lattice size (i.e., $`N_{el}`$) on the linear size $`L`$. We plot in Fig. 1, the mean fraction of broken bonds at failure $`p_f`$ as a function of $`N_{el}^{3/8}`$ for diamond and triangular lattices, which in principle should obey Eq. 6 as well. While to accept the percolation hypothesis one should observe a linear regime, a net curvature is apparent in the data especially for sizes $`L>100`$. A similar curvature is found in the 50% survival probability $`p_s`$.

We repeated the same analysis using the power law threshold distribution with $`\beta =1/20`$ in triangular lattices of sizes $`L=8,16,24,32,64,128,256,512`$. Table 2 presents mean and standard deviations in the broken bond density (fraction of broken bonds) at the peak load and at failure for power law thresholds distribution. The result is reported in Fig. 2, and is once again in contrast with percolation scaling (i.e. $`\nu 4/3`$)

Thus we can conclude that in the random fuse model with a uniform and power law distributions $`\nu `$ is not equal to $`4/3`$. On the basis of the results, however, we can not exclude the possibility that Eq. (6) is valid with a different value of $`\nu `$. This would correspond to some sort of correlated percolation, which, being a second order critical phenomenon, is expected to be universal with respect to the lattice structure. A way to test this idea is to plot the data as in Figs. 1-2 using a different value of $`\nu `$. In Figs. 3a and 3b we report similar plots for the triangular and diamond lattice topologies with uniform thresholds distribution, where it can be seen that a straight line is obtained only for large values of $`\nu `$. In addition, the $`\nu `$ values for which a straight line is obtained in Figs. 3a and 3b are quite different from one another. For example, a direct fit of the data for $`L16`$ yields

$`p_f0.0816`$ $`=`$ $`0.42N_{el}^{0.19},`$ (8)
for triangular lattices, and

$`p_f0.0751`$ $`=`$ $`0.57N_{el}^{0.25}.`$ (9)
for diamond lattices. Hence, Eqs. (8) and (9) for triangular and diamond lattice topologies, estimate the scaling exponent $`\nu `$ to be equal to $`2.63`$ and $`2.0`$, respectively.

An additional test is provided by plotting the standard deviation of bonds at failure, which in case of percolation should follow Eq. 7. For small lattice sizes the data follow a power law with an exponent $`\nu 1.56`$, but deviations occur at large sizes. In the framework proposed in Ref. hansen003  one could interpret the data in Fig. 4 as an indication of an initial correlated percolation process with $`\nu =1.56`$, which then crosses over to a localized state where scaling is obscured by the presence of a damage concentration profile. We will show in the next section, however, that there is apparently no localization at peak load (i.e. the maximum current before catastrophic failure). On the other hand, the damage standard deviation at peak load follows closely the behavior at failure as shown in the inset of Fig. 4. Thus the damage profile does not appear to be responsible from the deviation from scaling. A different explanation would be that there is no scaling just because the lattice fails abruptly far from a (correlated) percolation critical point.

## 5 Damage localization

As discussed in the previous section, damage is not described by percolation critical scaling up to the largest sizes. Fracture is abrupt and damage localizes. Here we clarify when and how localization takes place. In particular, we will consider the damage accumulated up to the peak load (i.e. the highest current that the lattice can bear without breaking) and after the peak load up to failure.

In Fig. 5 we display the damage profile $`p(y)`$ at failure and at peak load in a single simulation, for the cases of uniform and power law disorders. For uniformly distributed thresholds, localization appears clearly at failure, while at peak load the damage spreads homogeneously though the lattice. For power law disorder, it is difficult to assess from these curves the extent of localization. To obtain a quantitative description of the localization process it is necessary to average over different realizations.

Averaging the profiles is a delicate task since localization does not necessarily take place in the center, but can in principle occur anywhere along the length of the lattice. Thus one can not perform a simple average because this would yield a flat average profile irrespective of the individual profile shapes in a single realization. The authors of Ref. mikko04  proposed to first shift the profiles so that they are centered around the maximum and then average. This method emphasizes the noise too much, yielding a spurious cusp in the center. Another possibility is to shift instead by the center of mass of the damage or, to avoid any effects due to shifting, one can use the Fourier method.

We first consider the damage accumulated up to peak load shifting the data by the center of mass method. The result displayed in Fig. 6 for the uniform disorder distribution clearly shows that there is no localization at peak load. In the case of power law disorder the profile as presented in Fig. 7 is not completely flat for small scales, but flattens more and more as the size is increased. We tend thus to attribute the apparent profile to size effects. These data imply that the localization profiles already observed in the literature hansen003 ; hansen-other ; mikko04  reflect mostly the final breakdown event, after the final spanning crack is nucleated. To quantify the corresponding damage localization, we can then average the profiles $`\mathrm{\Delta }p(y)`$ obtained considering only the damage accumulated between the peak load and failure. In addition, this procedure considerably reduces the background noise.

The center of mass shifted averaged profiles for various system sizes are reported in Figs. 8a and 8b for uniform and power law disorder cases and these profiles show that the profile shapes decay exponentially at large system sizes. The damage peak $`\mathrm{\Delta }p(0)`$ scales as $`L^{0.3}`$ for the uniform distribution and is roughly constant for the power law distribution. We have tried different ways to collapse these profiles to extract a localization length. The simple linear scaling $`\mathrm{\Delta }p(y,L)/\mathrm{\Delta }p(0)=f((yL/2)/L)`$, proposed in Ref. hansen003 , does not yield a very good collapse (see Fig. 9a). A perfect collapse is instead obtained using the form

$$\mathrm{\Delta }p(y,L)/\mathrm{\Delta }p(0)=f(|yL/2|/\xi ),$$
(10)
where $`\mathrm{\Delta }p(0)=L^{0.3}`$ and $`\xi L^\alpha `$, with $`\alpha =0.8`$ (see Fig. 9b). The situation for power law disored is similar, but scaling is less precise. As shown in Fig. 8b, the profiles are decaying again exponentially. We could not, however, perform a reliable collapse: a linear scaling seems appropriate for system sizes $`L<256`$ but fails for larger sizes.

To obtain additional confirmation of these results, we perform a Fourier analysis of the $`\mathrm{\Delta }p`$ profiles, thus avoiding any possible bias due to the shifting. We first compute the magnitude of the Fourier transform for each realization and then average over disorder. From Eq. 10 we would expect the power spectrum of the profile to follow

$$|\stackrel{~}{p}(k)|^2/|\stackrel{~}{p}(1)|^2=\stackrel{~}{f}(kL^{0.8}).$$
(11)
This result applies to an infinite system and finite size deviations and other problems of the discrete Fourier transform are expected to affect the data. Nonetheless, as shown in Fig. 10, we can collapse reasonably well the curves using the same exponent as in real space.

It is also possible to estimate the localization length directly, independently of the profile averaging. This is done by computing the width of the damage cloud as $`W((y_b\overline{y_b})^2)^{1/2}`$, where $`y_b`$ is the $`y`$ coordinate of a broken bond and the average is taken over different realizations. We have first measured $`W`$ at peak load and at failure, and obtained a result $`WL`$, consistent with earlier results hansen91  (see Fig. 11). This result is expected, since for a uniformly distributed damage such as at the peak load, $`WL/\sqrt{12}0.288L`$, and is in excellent agreement with the numerical data. This is consistent with the fact that there is no localization at peak-load and at failure the damage cloud is dominated by the uniform fluctuations already present at peak load. Any additional scaling due to localization is obscured. To uncover it we can restrict the average only to the bonds broken in the last failure event (post-peak damage), obtaining $`WL^{0.81}`$ which is consistent with the data collapse of the profiles. It is interesting to notice that the width of the final crack scales as well with a similar exponent. This supports the idea that localization is produced by catastrophic failure which yields at the same time the crack and the damage profile.

## 6 Scaling of the damage density

It has been noted in the previous section that the final breakdown event is very different from the initial precursors. Thus, we consider the scaling of the number of broken bonds at the peak load, $`n_p`$, that excludes the last catastrophic event. In Fig. 12 we plot $`n_p`$ as a function of the lattice size $`N_{el}`$ for triangular and diamond lattices. The data displays a reasonable power law behavior $`n_pN_{el}^b`$, with $`b=0.93`$ and $`b=0.91`$ for triangular and diamond lattice, respectively as previously shown in Ref. delaplace . Thus the difference between the two lattices is marginal and may be attributed to the results obtained from the smaller lattice sizes, where corrections to the fractal scaling may exist. By plotting $`n_p/N_{el}^b`$ vs $`N_{el}`$ (see the inset of Fig. 12) we show that some systematic deviations appear. The data could be equally well fit by a linear law times a logarithmic correction $`n_pN_{el}/\mathrm{log}(N_{el})`$ as suggested in Ref. mikko04  (see Fig. 13). Both these fits imply that in the limit of large lattices the fraction of broken bonds prior to fracture vanishes (i.e. $`p_c=0`$ in the thermodynamic limit).

Figure 14 presents the scaling for the number of broken bonds, $`(n_fn_p)`$, after crossing the peak load. Once again, the scaling exponent for $`(n_fn_p)`$ is similar for both the triangular and diamond lattice topologies, and is equal to $`0.72`$ and $`0.69`$ respectively. This exponent is consistent with the parameters estimated from the profile $`(n_fn_p)<p(0)>(3L+1)\xi L^{1.5}`$, The behavior of power law disorder appears to be simpler. The number of broken bonds at peak load scales linearly with $`N_{el}`$, as shown in Fig. 15. A similar result holds asymptotically for $`(n_fn_p)`$ (see Fig. 15), which is again consistent with the profile scaling.

## 7 The failure probability distribution

The cumulative probability distribution for the fraction of broken bonds at failure (also termed as cumulative failure probability distribution), defined as the probability $`\mathrm{\Pi }_f(p_b,L)`$ that a system of size $`L`$ fails when the fraction of broken bonds equals $`p_b=\frac{n_b}{N_{el}}`$, where $`n_b`$ is the number of broken bonds, is plotted in Fig. 16 for varying triangular lattice system sizes. In Ref. hansen003  a data collapse of a different, but, related distribution (i.e., survival probability) was attempted using percolation scaling. As it is evident from the failure of such a scaling (i.e., Eq. (6)), this collapse is poor.

On the other hand, we obtain a very good collapse by simply plotting the distribution in terms of $`\overline{p}_f\frac{(n_b\mu _{n_f})}{\sigma _{n_f}}=\frac{(p_bp_f)}{\mathrm{\Delta }_{p_f}}`$, where $`\mu _{n_f}`$ and $`\sigma _{n_f}`$ denote the mean and standard deviation of the number of broken bonds at failure, and $`p_f`$ and $`\mathrm{\Delta }_{p_f}`$ denote the mean and standard deviation of fraction of broken bonds at failure. Fig. 17 shows that $`\mathrm{\Pi }_f(p,L)`$ may be expressed in a universal scaling form such that $`\mathrm{\Pi }_f(p,L)=\mathrm{\Pi }_f(\overline{p}_f)`$ for both triangular and diamond lattice topologies of different system sizes $`L`$. A similar collapse can be performed for the power law disorder distribution (see Fig. 18). The excellent collapse of the data in Figs. 17 and 18 suggests that the cumulative failure probability distribution, $`\mathrm{\Pi }_f(p_b,L)=\mathrm{\Pi }_f(\overline{p}_f)`$, may be universal in the sense that it is independent of lattice topology and disorder distribution.We have also checked that the distributions at peak load obey essentially the same laws, i.e., $`\mathrm{\Pi }(\overline{p})=\mathrm{\Pi }_f(\overline{p}_f)=\mathrm{\Pi }_p(\overline{p}_p)`$, where $`\mathrm{\Pi }_p(\overline{p}_p)`$ is the probability that a system of size $`L`$ is at the peak load when the fraction of broken bonds equals $`p_b`$, and $`\overline{p}_p`$ is the corresponding reduced variable at the peak load. Finally, the collapse of the data in Fig. 19 indicates that a Gaussian distribution adequately describes $`\mathrm{\Pi }_f`$.

The fact that damage is Gaussian distributed suggests that there is no divergent correlation length at failure consistent with the conclusions of Ref. delaplace  that reported a finite correlation length at failure. Long-range correlations in the damage would imply that the central limit theorem does not hold and hence the normal distribution would not be an adequate fit to the data. The absence of long-range corralation is again in agreement with the hypothesis that fracture is analogous to a first-order transition zrvs .

## 8 Discussion

This paper presents numerical simulations on large two-dimensional triangular and diamond lattice fuse networks with uniform and power law disorder distributions. We focus our analysis on the statistical properties and localization features of the damage as a function of the lattice type, size and disorder form. The use of high statistical sampling and relatively large lattice sizes is essential to obtain reliable results.

The picture emerging from our analysis is that, for strong disorder, damage accumulates first in the system in an uncorrelated (or short range correlated) manner. This process continues up to the peak load, where no apparent sign of localization is present. Further increase of the current leads to catastrophic failure through a large avalanche event, whose size scales with the lattice size as $`L^{1.45}`$. The damage accumulated in this event is localized in a cloud surrounding the final crack. The damage profile has exponential tails and can be collapsed for different system sizes using an appropriate scaling law.

The accumulated damage density at failure $`p_f`$ appears to be far from a percolation like critical point, since we could not find a reliable universal scaling law as the system size is varied. A different possibility that we have tested is that $`p_f0`$ as $`L\mathrm{}`$, but a fit using simple scaling forms show some deviation as well. Thus we can not conclude whether $`p_f`$ decays slowly to zero, or to a non vanishing asymptotic value. This value, however, would not necessarily coincide with a critical point.

Finally, this study also presents the scaling of cumulative failure probability distribution, which is defined as the probability that a lattice system fails at a given fraction of broken bonds. Based on the numerical results presented, we show that the cumulative failure probability distribution is universal in the sense that it does not depend on the lattice topology, i.e., the distributions are identical for both triangular and diamond lattice topologies. Furthermore, a normal distribution presents an adequate fit to the data.

This research is sponsored by the Mathematical, Information and Computational Sciences Division, Office of Advanced Scientific Computing Research, U.S. Department of Energy under contract number DE-AC05-00OR22725 with UT-Battelle, LLC. We thank M. J. Alava, A. Hansen, H. J. Herrmann and S. Roux for useful remarks and discussions.

## References


format: text

text:
# WERE THE LMC GLOBULAR CLUSTERS FORMED IN A DISK?

## 1 INTRODUCTION

In a pioneering investigation of the clusters in the Large Magellanic Cloud, Freeman, Illingworfth & Oemler (1983) surmised that “\[T\]here is no old, kinematical halo population among the clusters of the LMC.” At that time radial velocities were only available for nine old clusters in the Large Cloud. Presently 13 globular clusters are known to be associated with the Large Magellanic Cloud (Schommer 1991, Schommer et al. 1992, Suntzeff 1992, van den Bergh 2000, p. 104). From their discussion of the kinematics of old clusters in the LMC Schommer et al. (1992) concluded that: “The oldest clusters still present an enigma; they do not have the kinematics of an isothermal, or slowly rotating, pressure-supported halo. These objects rotate with an amplitude comparable to that of the younger disk”. Taken at face value this conclusion would appear to indicate that the oldest clusters in the Large Cloud formed quite differently from those in M33 (Schommer et al. 1991, Chandar et al. 2002), in which the oldest (globular) clusters are observed to have halo kinematics, whereas younger (open) clusters exhibit disk-like motions. The view that late-type galaxies such as M33 (M V = -18.9) and the LMC (M V = -18.5), which have similar uminosities, had very different evolutionary histories would not fit comfortably with most current views of galaxy evolution. This prompts one to ask how secure is the conclusion that the oldest massive clusters in the Large Cloud formed in a disk, rather than in the LMC halo?

## 2 KINEMATICS OF OLDEST LMC CLUSTERS

In their Table 2 Schommer et al. (1992) provide radial velocities (corrected for an Eastward transverse velocity of 150 km s-1 for the 13 oldest LMC clusters which belong to SWB class VII (Searle, Wilkinson & Bagnuolo 1980). These data are plotted in Figure 1, which shows that the rms velocity dispersion around the LMC HI rotation curve (thin line) is 28 km s -1, compared to a dispersion of 30 km s-1 around the mean velocity (dashed line) of the SWB VII clusters. This result shows that the H I rotation curve does not provide a significantly better fit to the the kinematics of the oldest LMC clusters than does a model in which these clusters are not assumed to be supported by rotation. In this respect the old clusters of class WWB VII differ significantly from the intermediate-age clusters of SWB classes IV, V and VI which (see Figure 5 of Schommer et al. 1992) clearly exhibit a significant rotational velocity component. It is interesting to note that the 30 $`\pm `$6 km s-1 velocity dispersion of the 13 SWB VII clusters appears to be smaller, albeit only at the 2 $`\sigma `$ level, than the 53 $`\pm `$10 km s-1 velocity dispersion that Minniti et al. (2003) have recently found for 42 RR Lyrae variables in the Large Cloud. However, it is noted that the velocity dispersion of the LMC RR Lyrae variables is comparable to the 54 $`\pm `$8 km s-1 velocity dispersion that Chandar et al. obtain for the old M33 globular clusters with $`R<2.25`$ kpc, but smaller than the 82 $`\pm `$13 km s-1 dispersion that these authors find for the old M33 halo clusters with $`R>2.25`$ kpc.

## 3 CONCLUSIONS

By lumping old and intermediate-age LMC clusters together Schommer et al. (1992) found that these clusters exhibited significant rotation. However, the data plotted in Figure 1 show that the oldest (SWB class VII) clusters do not exhibit a significant rotational component. In this respect the Large Cloud cluster system resembles M33 in which Chandar et al. find that 85 $`\pm `$ 5 percent of the oldest clusters have halo kinematics. It is therefore concluded that worries about the apparent difference between the kinematics of the oldest population components in M33 and in the LMC (van den Bergh 2000, p. 83) were probably premature.

I thank Ken Freeman for pointing out to me that the number of old LMC globular clusters might not be large enough to prove the conclusion that these clusters have disk kinematics. Thanks are also due to Peter Stetson for statistical advice.


text:
# Time scales involved in market emergence

## 1 Introduction

A number of studies carried out over the past years showed that the time evolution of individual securities depends strongly on the evolution of other securities and even of the whole market. This observation leaded to the development of various models which help investors to minimize risk and to choose the optimal investment strategies. Therefore the magnitude, the temporal stability and the time-scale characteristics of the correlations are crucial factors for these models. Also from the theoretical point of view the existence and strength of correlations play an important role in development of proper models of the stock market dynamics and could help understanding the mechanisms which are responsible for the emergence of the collective signals out of noise in complex systems.

It is well-documented in literature that what dominates the market dynamics at the microscopic level is noise . The movements of stock’s price are governed by a series of buy/sell orders reaching the market essentially at random moments (although there are some long-time dependences which can be a source of the non-Gaussian tails of the distributions of the inter-transaction time intervals). Among the factors leading to this microscale randomness there is a difference in reaction times of the investors to arriving news. This can be well related to the different investment time horizons, acting through different market makers and so on. Thus, even though an important piece of news arrives on the market, different investors absorb this information and adjust their positions at distinct moments. Similarly, there exists randomness in the transaction volume, bid-ask spread etc. On short time scales, all these elements cause the price to fluctuate stochastically around its “true” value like it is in a random walk. In these circumstances, if there is some amount of persistence in the price evolution, it can be observed only after many transactions take place, i.e. on longer time scales.

The above-described randomness in the price movements is even better evident after a parallel inspection of the tick-by-tick data for two or more assets is made. Apart from the-already-mentioned difference in the reaction time of the investors for the same piece of news, there is a separate news flow regarding each of the companies under study, which can be yet another source of randomness. Therefore, we can safely assume that on very short time scales comparable with the mean time interval between consecutive trades, the correlations among the stock price fluctuations do not differ from the noise level and they are insignificant. Consequently, on such short time scales it is not justified to consider the market as a coherent whole; one rather deals with a set of the elements evolving independently from each other. Going from short time scales to the ones much longer than the mean inter-transaction interval, new effects occur. Firstly, all the investors have opportunity to react to the news, which gives the complete picture of how this piece of news affects the price. Secondly, the investors analyse price changes of other assets and correct their positions and strategies accordingly. In the price evolution of each asset there is thus information of other assets’ prices which causes the inter-stock correlations to emerge. This diffusion of information increases with increasing time, bringing about the correlations to be strenghtened either. Moreover, the longer time passed, the more investors manage to act, which leads to an even wider flow of information between the assets and the time scales . From this angle, the coupling strength reaches its maximum after all the investors can correct their positions. Macroscopically, the existence of the inter-stock correlations, both those originating from similar responses of different stocks to the same piece of news and those being the effect of a directed network of influence among the stocks , is an important aspect of the collective market formation (see also ). The most striking evidence of this phenomenon is the strong index movements and trends which cannot be observed in a completely decorrelated system.

In our paper, we would like to address the closely related question of what are the time scales at which the significant correlations emerge. i.e. at which time scale range the transition from noise to a collective behaviour takes place, and what are the factors responsible for the inter-stock coupling strength. We shall compare the properties of the inter-stock correlations at several different time scales for different groups of stocks by applying the globally-oriented correlation matrix formalism in order to look at the couplings between more than two assets. Additionally, motivated by the results from our earlier study  dealing with the statistical properties of the stock price fluctuations, we intend to compare the correlation properties of the contemporaneous and the historical data.

## 2 Results

We analyze high frequency data from the American and from the German stock market contained in the TAQ and the KKMDB databases , respectively. In both cases, the data covers the over-two-years-long period from Dec 1, 1997 to Dec 31, 1999. For each company listed on NYSE and NASDAQ markets, the TAQ database contains a record of all transactions which took place within a given time interval of trading (the same refers to the KKMDB database and the Deutsche Börse). As the transactions are made at random moments, first we need to create a time series of price values being sampled with constant frequency. Following the standard prescription, we assume that the price $`x_\alpha (t_i)`$ of an asset $`\alpha `$ at time $`t_i`$ is equal to the price of the last preceiding transaction on the corresponding stock. Then, given a time scale $`\mathrm{\Delta }t`$, we calculate a time series of normalized returns defined by

$$g_\beta (t_i)=\frac{G_\beta (t_i)G_\beta (t_i)_{t_i}}{\sigma (G_\beta )},\sigma (G_\beta )=\sqrt{G_\beta ^2(t_i)_{t_i}G_\beta (t_i)_{t_i}^2}$$
(1)
where

$$G_\beta (t_i)=\mathrm{ln}x_\beta (t_i+\mathrm{\Delta }t)\mathrm{ln}x_\beta (t_i).$$
(2)
Here $`\mathrm{}_{t_i}`$ stands for averaging over discrete time. Let us say we have a set of $`N`$ stocks and from the corresponding time series of length $`T`$ we construct an $`N\times T`$ matrix $`𝐌`$ and finally we calculate an $`N\times N`$ correlation matrix $`𝐂`$ according to the formula

$$𝐂=(1/T)\mathrm{𝐌𝐌}^\mathrm{T}.$$
(3)
Each element $`C_{\alpha ,\beta }`$ of the correlation matrix is simply the correlation coefficient for the pair of stocks $`\alpha `$ and $`\beta `$. Afterwards, the correlation matrix can be diagonalized in order to obtain the spectrum of its eigenvalues $`\lambda _k`$ ($`k=1,\mathrm{},N`$. We repeat this procedure for several distinct time scales ranging from 1 min (or even 1 s) up to two trading days (780 min in New York and 1020 min in Frankfurt). Typically for the stock market data, the correlation matrix develops at least one eigenvalue which is repelled from the rest of the spectrum and which corresponds to a collective behaviour of group of stocks or a whole market . It is convenient to confront the eigenvalue spectrum with the universal predictions of the Random Matrix Theory as all the apparent discrepancies can be related to the existence of market-specific information.

We start from an investigation of the coupling strength’s dependence on the time scale $`\mathrm{\Delta }t`$ for pairs of stocks. It is convenient to quantify this dependence in terms of the correlation coefficient $`C_{\alpha ,\beta }(\mathrm{\Delta }t)`$. In Figure 1 we show this quantity for two different but typical pairs of the DJI stocks: Alcoa (AA) $``$ Exxon (XON) and Chevron (CHV) $``$ Exxon. This Figure shows that $`C_{\alpha ,\beta }`$ is definitely not invariant under a change of the time scale. Essentially, the magnitude of the correlation coefficient increases while going from the smaller intra-hour or minute to the larger daily $`\mathrm{\Delta }t`$’s. As CHV and XON belong to the same market sector (energy), while AA does not, for each $`\mathrm{\Delta }t`$ the correlations in the former case are much more significant then in the latter one. Nevertheless, both plots demonstrate only small correlations at the shortest analyzed time scale of 1 min while for larger $`\mathrm{\Delta }t`$ the correlation coefficient significantly increases. For both pairs of stocks the picture is qualitatively similar up to $`\mathrm{\Delta }t20`$ min but then we observe a difference in the behaviour of the correlations for larger $`\mathrm{\Delta }t`$’s: $`C_{\mathrm{CHV},\mathrm{XON}}`$ still goes up with increasing the time scale reaching as high as 0.65 for daily returns, while $`C_{\mathrm{AA},\mathrm{XON}}`$ almost saturates below 0.20. However, a trace of saturation is also identifiable for the CHV-XON pair.

Due to the fact that such an increase of correlation magnitude is characteristic for all pairs of stocks, one may expect that a similar effect can also be observed if one looks at a more global measure of correlations e.g. the largest eigenvalue of the correlation matrix. Thus we select two sets of 30 stocks listed in the DJIA and DAX indices and we evaluate the corresponding function $`\lambda _1(\mathrm{\Delta }t)`$; the results are displayed in Figure 2. Indeed, for both DJIA (Fig. 2(a)) and DAX (Fig. 2(b)) this effect is strong, although the detailed behaviour of the largest eigenvalue is not market invariant. In DJIA, $`\lambda _1`$ increases for $`\mathrm{\Delta }t`$ up to 30 min; for longer time scales there is no further increase but rather a saturation of $`\lambda _1`$ can be seen. In contrast, the largest eigenvalue for DAX gradually rises up to a daily time scale ($`\mathrm{\Delta }t=510`$ min); this increase, however, is rather slow for $`\mathrm{\Delta }t30`$ min. These results show how the markets change their behaviour from decorrelated and completely noisy dynamics to the collective one. This observation can be compared with earlier analysis of a market-sector formation while going from short to long time scales presented in ref. . Interestingly, for all time scales the DAX stocks seem to be more strongly coupled than their DJI counterparts; for $`\mathrm{\Delta }t=1`$ min $`\lambda _1`$ for DJI only moderately differs from the RMT prediction ($`\lambda _1^{\mathrm{RMT}}(\mathrm{\Delta }t=1min)1.0`$. The innately more correlated nature of the German market, which leaves its fingerprints especially for $`\mathrm{\Delta }t>30`$ min, has already been pointed out earlier (see eg. ) for daily returns. Magnitude of this effect is substantially time-dependent, however (compare with the results from 1990-2001 in ).

These results go in parallel with the so-called Epps effect , named after the first researcher who demonstrated that when going from the daily to the intra-hour time scales the inter-stock correlations decay. Althoug his analysis was entirely based on the stock market data, the similar results were obtained later for the currency exchange markets as well (see  for recent results).

Price correlations among different stocks at high frequencies are caused by market makers who quickly react to important news and to changes of prices of other securities. As a security’s price can be determined only in a transaction, statistically a piece of news influences those securities first which are traded more frequently than others. This of course implies that also correlations are more likely to occur earlier for the most active securities. This can be seen for the currency exchange rates where the correlations are significant already on the time scales shorter by an order of magnitude than in the case of the stock markets (e.g. ). In order to compare the size of correlations between stocks of different transaction frequencies (being positively correlated with market capitalization of the corresponding companies), we select a few distinct sets of 30 stocks in such a way that the companies within each set are characterized by similar capitalization (from $`10^8`$$ to $`10^{11}`$$). For each set of stocks we evaluate $`\lambda _1(\mathrm{\Delta }t)`$ and compare it across the sets (Figure 3). The so-quantified correlation range shows a systematic and monotoneous dependence on the capitalization: for a given $`\mathrm{\Delta }t`$, the larger the company, the stronger the average coupling with its same-size counterparts. $`\lambda _1(\mathrm{\Delta }t)`$ for the DJI market being a mixture of $`10^910^{11}`$$ companies is also presented (diamonds in Fig. 3). In the case of the smallest firms considered (triangles down), $`\lambda _1`$ is essentially at the RMT level for all the time scales up to $`\mathrm{\Delta }t=30`$ min. It is interesting to note that a saturation level occurs only for the largest companies worth at least $`10^{11}`$$ each (circles) and for the DJI stocks; all the other groups of companies display increase of $`\lambda _1`$ even for the largest time scales analyzed. It may be hypothesized that the saturation could manifest itself on time scales much longer than two days $``$ the smaller the companies, the later the saturation. From this point of view it is clear why $`\lambda _1`$ for the DJI stocks saturates earlier than for the DAX stocks (Fig.2): as pointed out in , the average number of transactions for the DAX companies is significantly smaller than for the DJI ones.

One of the possible sources of the analyzed effect can be the lack of synchronicity in transaction moments for different securities  and the associated nonsyncronicity of their price determination. We choose a pair of the most correlated stocks in DJIA: Citigroup (C) and General Electric (GE) and by following ref.  we remove from the original data all transactions which didn’t take place simultaneously for both stocks. (By “simultaneous” we understand transactions which were made within the same second.) Obviously, this strongly reduces the average number of transactions per business day (e.g. for GE: from 4240 to 440) but nevertheless there is still more than one simultaneous transaction per minute. Next we proceed in the usual way by creating a time series of $`\mathrm{\Delta }t`$-returns and then by calculating the correlation coefficient $`C_{\mathrm{C},\mathrm{GE}}(\mathrm{\Delta }t)`$. Figure 4(a) displays the functional dependence of $`C_{\mathrm{C},\mathrm{GE}}`$ on the time scale for the original data comprising all the transactions (circles) and for the synchronous data only (squares). For $`\mathrm{\Delta }t<60`$ min, values of the correlation coefficient for the synchronous data are elevated in respect to the nonsynchronous data with the difference increasing with decreasing $`\mathrm{\Delta }t`$. For the shortest 1 min time scale, the elimination of nonsynchronous transactions almost doubles $`C_{\mathrm{C},\mathrm{GE}}`$. In Fig. 4(b) analogous calculation is carried out by using data from the two most active NASDAQ stocks: Dell Computer (DELL) and Intel Corp. (INTC). As the average nonsynchronous transaction number per a business day exceeds 20,000 with over 10,000 simultaneous transactions per day, here the correlations are detectable even at a-few-seconds time scales. The two analyzed data sets differ from each other only for $`\mathrm{\Delta }t<30`$ s; this difference is, however, impressive: for $`\mathrm{\Delta }t=1`$ s $`C_{\mathrm{DELL},\mathrm{INTC}}`$ is about an order of magnitude larger for the synchronous data than in the other case.

A straightforward generalization of this procedure towards more than two degrees of freedom (stocks) requires some care, though. Selection of precisely synchronous transactions leads to a drastic reduction of the time resolution and of the range of available time scales, therefore this proves inefficient for more than two stocks. In order to overcome this problem, we weaken our definition of synchronicity by introducing a tolerance parameter $`\tau `$; now transactions are considered to be synchronous if they are made within an interval $`(t_\alpha \tau ,t_\alpha +\tau )`$, where $`t_\alpha `$ stands for the transaction time for a reference stock $`\alpha `$ being the most active stock in a set of the stocks under study. The functional dependence of $`\lambda _1`$ on $`\mathrm{\Delta }t`$ for a set of the 10 most frequently traded stocks and for several values of $`\tau >0`$ is presented in Figure 5. It is clear from the figure that the more synchronous transactions are considered ($`\tau 0`$), the stronger short-time-scale couplings among the stocks occur.

Nonsynchronicity of trading being related to the microscopic randomness of the stock market dynamics cannot alone account for the observed size of the correlations decay and thus there must be other influential factors here, like for example the possible existence of lagged correlations amongst the assets  or the differences in the time horizon of trading strategies of individual market agents . Since these factors cannot be directly implemented in the correlation matrix formalism, we shall not discuss this issue in the present paper and instead we refer reader to the literature (see also e.g. ).

Owing to the fact that the high frequency data from the American stock market has been extensively studied over past decades, we can compare outcomes from our study with the ones from other studies. In his original paper , Epps investigated data for the stocks of the automobile sector, recorded during a few months of 1971. In Figure 6 we confront the Epps’ historically distant results with the more contemporary ones from 1998-1999 for an exemplary pair of stocks (Ford Motor - F, General Motors \- GM), studied also in ref. . The phenomenon of a decrease of the correlation coefficient with decreasing $`\mathrm{\Delta }t`$ is evident in both cases, but at present it is strongly shifted towards the shorter time scales. In contrast, for the daily time scales $`C_{\mathrm{F},\mathrm{GM}}`$ assumes comparable magnitudes for both 1971 and 1998-99. In another study Andersen et al. (ref. ), who analyzed data from the time period 1993-1998, showed that the average value of the correlation coefficient calculated for 30 DJI stocks assumes 0.12 at the time scale of 5 min. In our case, this is equivalent to the averaged correlation matrix element, which for $`\mathrm{\Delta }t=5`$ min equals 0.19. Although this comparison is not fully decisive because of only one time scale investigated, it suggests that the short-time-scale couplings among the DJI stocks were stronger at the end of 1990’s than they used to be on average during this decade; again, this is in the spirit of the above conclusions.

## 3 Conclusions

For a summary, we study the process of the emergence of a collective market out of noise by investigating a time-scale dependence of the magnitude of the inter-stock correlations for the companies listed in the American and German markets. We observe that the correlations’ magnitude quantified in terms of the Pearson’s correlation coefficient increases with increasing the time scale from a statistically insignificant level at the time scales comparable with the average inter-transaction interval to high values at the long hourly or daily time scales. By applying the correlation matrix formalism to the high-frequency data, in a simple way we generalize these results to more than two degrees of freedom. We show that such behaviour of the inter-stock correlations can be observed also globally for the market as a whole. Our results convince us that one of the most important factors that determines the time scales at which the collective behaviour of assets occur is the trading frequency: for a given time scale, the most active stocks are also among the most correlated ones. In contrast, the stocks of small companies which are characterized by a low number of transactions, present non-significant correlations up to daily time scales. We next demonstrate that the synchronous data with a supressed level of randomness of transaction moments reveals stronger couplings than the original data. Finally, our results provide us with the indication that nowadays the collective market emerges at significantly shorter time scales than it used to do in the more distant history, i.e. for the large American companies the market shows the trace of a weak collectivity already at the minute or even the intra-minute time scales compared with the intra-hour scales previously.

This result recalls a congenial effect of a faster convergence of the stock returns distributions towards a Gaussian in the contemporary data if compared with the historical data. As documented in ref. , in the same period 1998-99, the scaling law with the exponent $`\alpha 3.0`$ breaks already at the time scales of tens of minutes, while earlier it was still hold even at the time scales of several days . We interpreted this phenomenon as being a direct consequence of a faster information processing and a faster loss of memory in the evolution of the market when going from past to present. Now, while looking at correlations among the stocks at various time scales, we receive a further firm support for such conclusions. Due to the fact that nowadays the emergence of the collective dynamics of stocks can be observed at time scales much shorter then before, and keeping in mind that one of its fundamental governing factors is an asset’s trading frequency, we shall underline that a time flow in the stock market (and possibly in other financial markets as well) is not constant over long periods but instead, owing to the technological progress and an enlarged flow of the arriving information, the market time tends to accelerate: effectively, one day in 1980 might not be completely equivalent to one day in 2000. A straightforward and far-reaching consequence of this fact is that potential models of the financial dynamics which do not take this acceleration into account and which assume that properties of the market dynamics are time-invariant might be not fully adequate. It is also interesting to notice that such an acceleration of the financial market evolution as viewed from the linear time scale is consistent, at least qualitatively, with a scenario priovided by the log-periodicity effect, especially the one that refers to the last 200 years .


text:
# The Extragalactic Origin of the Arcturus Group

## 1 Introduction

The discovery of substructure in the distribution of stars in the outer regions of galaxies like the Milky Way and M31 constitutes indisputable evidence that accretion events and mergers have played a role in building up their stellar halos (Ibata et al. 1994, Helmi et al. 1999, Majewski et al. 2003). Although the role of such events in the formation of the Galactic disk components is not yet fully understood, models of the thick disk component typically envision the tidal heating of an extant thin disk by a merging satellite (Quinn & Goodman 1986). Recent numerical simulations of disk galaxy formation also suggest that disrupted satellites might have contributed a significant fraction of the old stars in the disk of the Galaxy (Abadi et al. 2003a,b). As with the stellar halo, long-lived traces of such events may be preserved in the form of substructure in the orbital distribution of stars in the disk (Helmi et al 2003).

Identifying such debris ought to be easier in samples where contamination by the thin disk component is minimized, such as samples of metal poor/old stars, or else samples collected in situ away from the disk or far from the Galactic center. The discovery of a “ring” of stars of almost certain extragalactic origin in the outer disk of the Milky Way by the SDSS (Newberg et al 2002, Yanny et al 2003), as well as recent report of a galaxy being cannibalized in the Milky Way disk (Martin et al 2003) provide examples that this process indeed seems to be at work in the Galaxy (Helmi et al 2003, Rocha-Pinto et al 2003).

Gilmore, Wyse & Norris (2002) have also recently argued for substructure in the Galactic thick disk. They postulate the existence of a distinct population of stars, with kinematics intermediate between the canonical thick disk and the stellar halo (mean rotational velocity of $`100\mathrm{km}\mathrm{s}^1`$) to explain the results of their spectroscopic survey of F–G stars located 0.5–5 kpc away from the Galactic plane.

Substructure in the disk has been noted since the early days of Galactic structure research (see, e.g., Eggen 1998 and references therein), and has usually been interpreted as signature of the gradual dissolution of loosely bound star clusters in the clumpy potential of the disk. One example is the group of stars associated by Eggen (1971) with the bright star Arcturus, whose rotation speed lags the Local Standard of Rest (LSR) by roughly $`100\mathrm{km}\mathrm{s}^1`$. We discuss below that evidence for a population of stars with rotation speeds lagging the LSR by $`100\mathrm{km}\mathrm{s}^1`$ is also present in the metal-poor star catalog of Beers et al (2000; hereafter B00), and in the compilation of Gratton et al (2003, hereafter GCCLB).

In this Letter we reassess the status of Eggen’s Arcturus group by using modern distance, proper motions, radial velocities, and metallicity determinations available in the literature. We argue that the Arcturus group; the structure in the B00 and GCCLB compilations; as well as the kinematically-odd population of stars discovered by Gilmore, Wyse & Norris (2002), are all part of the same, dynamically coherent group of stars that might have originated in the disruption of a fairly massive satellite early during the Galaxy’s assembly.

## 2 The Arcturus group

### 2.1 Eggen’s stars

Over the years, Eggen compiled a list of stars kinematically similar to Arcturus from catalogs of nearby stars with large proper motions (Eggen 1971a; 1987; 1996; 1998). Distances to some of these stars, however, were unavailable at the time, and Eggen would assign them distances so that they had negligible velocity dispersion in $`V`$ ($`\sigma _V0.5`$ $`\mathrm{km}\mathrm{s}^1`$), 111We use $`U`$, $`V`$, and $`W`$ to denote radial, tangential and vertical velocities in Galactic coordinates, measured relative to the LSR. $`U`$ is positive outwards, $`V`$ in the direction of the Sun’s rotation, and $`W`$ towards the North Galactic Pole. We assume that the Sun’s ($`U`$,$`V`$,$`W`$) velocity is ($`9`$,$`12`$,$`7`$) $`\mathrm{km}\mathrm{s}^1`$. as required by his interpretation. This somewhat arbitrary step has made Eggen’s interpretation controversial but, remarkably, his choice led most candidates to lie along relatively narrow isochrones, enabling Eggen to derive ages for these stars. He concluded that most Arcturus group members are old ($`\tau >10`$ Gyrs) and metal poor (\[Fe/H\]$`0.6`$; Griffin & Lynas-Gray, 1999).

The shaded histogram (labelled “Eggen”) in Fig. 1a shows the distribution of the vertical component of the angular momentum, $`J_z`$, for stars listed by Eggen (1996, 1998) as members of the Arcturus group. We have included here only stars with available proper motions and Hipparcos parallaxes. This effectively restricts Eggen’s sample to $`46`$ stars within a sphere of radius $`300`$ pc centered on the Sun. The dispersion in $`J_z`$ (corresponding to a dispersion in $`V`$ of $`\sigma _V50\mathrm{km}\mathrm{s}^1`$) is much larger than demanded by Eggen’s interpretation of this association as a dissolving star cluster. Still, many of the candidates are indeed on orbits that lag the Sun’s rotation speed by $`100\mathrm{km}\mathrm{s}^1`$. Furthermore, many of these stars are on disk-like orbits; $`70\%`$ of them have vertical speeds that do not exceed $`50\mathrm{km}\mathrm{s}^1`$ and do not climb farther than $`1`$ kpc out of the Galactic plane. These stars show no obvious net motion in the radial direction.

### 2.2 Metal-deficient stars in the solar neighborhood

B00 provide a catalog of metal-deficient stars compiled with the specific aim of minimizing kinematic selection biases. This is the largest sample of metal-poor (\[Fe/H\] $`0.6`$) stars with available abundances, distances, and radial velocities. Proper motions are also available for many of these stars, making it a good sample to search for substructure in the orbits of metal-deficient stars (see, e.g., Chiba & Beers 2000; Brook et al. 2003).

The top dot-dashed histogram in Fig. 1a shows the angular momentum distribution for all stars in the B00 sample. The full sample is dominated by the canonical “thick disk” component, whose $`J_z`$ distribution peaks at $`1300\mathrm{km}\mathrm{s}^1`$ kpc, and lags the LSR rotation speed by $`50`$$`60`$ $`\mathrm{km}\mathrm{s}^1`$. The halo component contributes a significant number of counter-rotating stars, and is fairly prominent in this sample biased towards metal-poor stars.

The B00 distribution shows a number of prominent features that persist even when the sample is restricted to stars with disk-like kinematics ($`|W|50\mathrm{km}\mathrm{s}^1`$, see the lower dot-dashed histogram in Fig. 1a). We shall focus on two of these features, also seen in the distribution of rotation-speeds shown in Fig. 1b ($`V_{\mathrm{rot}}=J_z/R`$, where $`R`$ is the polar radial distance measured in the meridional plane).

The first feature, coincident with the peak velocity expected for the canonical thick disk, likely signals the “Hercules” stream (Eggen 1971b). This stream has $`V60\mathrm{km}\mathrm{s}^1`$ and $`U35`$ $`\mathrm{km}\mathrm{s}^1`$ (Skuljan et al. 1999), and is thought to originate from dynamical resonances induced by the Galactic bar (Dehnen 2000, Fux 2001). Intriguingly, the second feature, at $`J_z950`$ $`\mathrm{km}\mathrm{s}^1`$ kpc, has an angular momentum similar to that of Arcturus. The same feature is also seen in the GCCLB compilation of nearby metal-poor stars (see solid histogram in Figure 1a).

The statistical significance of this feature is difficult to assess, since the sample is subject to a number of subtle and complex biases and selection effects. We can nevertheless estimate it roughly by comparing the actual number of stars with velocities similar to that of Arcturus, with that expected for a random sample drawn from the standard Galactic components: thin disk, thick disk and stellar halo. This is shown in Fig. 1b, where the $`V_{\mathrm{rot}}`$ distribution has been decomposed into three Gaussians. From left to right, the halo, thick disk, and thin disk components are assumed to have ($`\overline{V}_{\mathrm{rot}}`$, $`\sigma `$) equal to ($`0`$,$`110`$), ($`160`$,$`50`$) and ($`220`$,$`25`$) $`\mathrm{km}\mathrm{s}^1`$, respectively.

The halo distribution is normalized to match the number of counter-rotating stars, whereas the relative contributions of the thick and thin disks are chosen to match the $`V_{\mathrm{rot}}`$ distribution of the remaining co-rotating stars. The fit obtained is shown by the thick solid curve in Fig. 1b. Most random realizations are unable to account for the “excess” of stars labelled “Arcturus group”; only $`16`$ out of $`1000`$ realizations give a number of stars comparable to what is observed. Such odds seem compelling enough to warrant exploring further the provenance of stars with such peculiar kinematics.

### 2.3 The origin of the Arcturus group

Could the Arcturus group be the result of perturbations to the orbits of disk stars induced by the Galactic bar? As discussed above, detailed modeling of such perturbations have linked stellar groupings in phase space to orbital resonances in a rotating barred potential (see, e.g., Dehnen 2000, Fux 2001). However, these perturbations typically induce small velocity changes ($`20`$$`50`$ $`\mathrm{km}\mathrm{s}^1`$), and therefore are unlikely to be responsible for the $`110`$$`130`$ $`\mathrm{km}\mathrm{s}^1`$ lag of stars in the Arcturus group. Moreover, such perturbations tend to produce a net mean radial velocity with respect to the LSR, which is not observed neither for the B00 candidates nor for stars in the original Eggen’s compilation (see Fig. 2).

#### 2.3.1 The group as tidal debris

A compelling alternative is that Arcturus is part of the tidal debris of a disrupted satellite, akin to the interpretation of the SDSS “ring” proposed by Helmi et al (2003). Intriguingly, the angular momentum of Arcturus is roughly half of that of a circular orbit at the solar circle, reminiscent of the dynamics of the SDSS “ring”. A further similarity with the outer Galactic “ring” is the lack of net radial motion, which suggests that these structures have become recognizable as stars of similar energy crowd at the apocenter of their orbits. As discussed below, an extragalactic origin for the Arcturus group would explain its sizeable velocity dispersion, as well as its vertical extension and distinct metal abundance.

We can test the plausibility of this scenario by following the disruption of a small satellite with orbit similar to that of Arcturus in the potential of the Galaxy. To this aim, we have integrated Arcturus’ orbit $`8`$ Gyrs back in time and use those orbital parameters to release a $`10^8\mathrm{M}_{}`$ satellite (modelled with $`10^5`$ particles after a King sphere with core and tidal radius $`233`$ pc and $`1.57`$ kpc, respectively). The numerical simulations include the self-gravity of the satellite using a quadrupole expansion of the internal potential (Zaritsky & White 1988), but are otherwise run in a fixed Galactic potential, of the form given by Johnston et al. (1996). We shall focus on particles (“stars”) located in a small volume around the Solar neighborhood, 5 to 8 Gyr after the satellite is disrupted.

The open circles in Fig. 2 show the velocities of particles located within 500 pc of the “Sun” (assumed to be at 8.5 kpc from the center of the Galaxy, and on the Galactic plane), which should be compared to stream candidates taken from the Eggen (solid black circles) and GCCLB (solid gray/red circles) samples. The good agreement between the mean velocities of the candidate stars and particles in the simulation indicate that the dynamics of the stream is consistent with that expected for the debris of a shredded satellite. Clues to the original mass of the satellite are preserved in the velocity dispersion of the group, but will depend on the number of streams expected to run through the solar neighbourhood at present. This is very hard to ascertain given current uncertainties in the velocities.

#### 2.3.2 Chemical properties of the Arcturus group

Supporting evidence for the “shredded-satellite” interpretation may be found in the chemical properties of stars in the Arcturus group. The Beers et al (2000) catalog lists \[Fe/H\] estimates for all stars, and further examination shows that stars with angular momentum consistent with the Arcturus group span a wide range of metallicities, roughly $`2.5<`$\[Fe/H\]$`<0.5`$. Does this broad metallicity distribution argue against Eggen’s suggestion that these stars are part of a disrupting cluster born out of the same molecular cloud ? It is indeed clear that most of the B00 stars in the appropriate $`J_z`$ range are metal-poor, but we would expect to find a significant contamination from halo and metal-poor thick disk stars in this range of $`J_z`$ (c.f. Fig 1b).

On the other hand, the tightness of the \[$`\alpha `$/Fe\] relation for the group candidates in Fig. 3, compared with the non-group candidates, does suggest a common brief star formation history for the stars in this range of $`J_z`$. This is best shown using the GCCLB dataset, which includes abundance estimates for elements other than Fe. This dataset may thus be used to search for potential features in the chemical properties of group candidate stars. In Figure 3, open symbols are used for all stars in the GCCLB sample, whilst filled circles are used for Arcturus group candidates; i.e., stars with $`J_z`$ in the range $`(700,1100)`$ $`\mathrm{km}\mathrm{s}^1`$ kpc, and $`|W|<50`$ $`\mathrm{km}\mathrm{s}^1`$. Candidate stars of the group define a tight sequence in the \[$`\alpha `$/Fe\]–\[Fe/H\] plane, much tighter than expected for stars of similar \[Fe/H\] selected at random from the full GCCLB sample.

Interestingly, the group candidates appear to follow, with minimal scatter, a well defined enrichment pattern roughly consistent with simple self-enrichment models, where the overabundance in \[$`\alpha `$/Fe\] (relative to solar) suggests a gas consumption timescale short compared with that of supernova type Ia ignition. The curve in Fig. 3 shows the predictions of Matteucci & Francois (1989) for a one-zone model222The curve corresponds to \[Si/Fe\] vs \[Fe/H\] in their model 1. Results are similar for other $`\alpha `$ elements. that self-enriches to \[Fe/H\]$`=1`$ in about $`1`$ Gyr, before the Fe-rich yield of SNIa start to play a role, leading to a sudden change in the slope of the curve. The good agreement between the model prediction and the candidate stars should be taken as a demonstration of the plausibility of our interpretation rather than as an endorsement of this particular model; other choices of model parameters may lead to similarly good fits, although all are likely to require a short timescale for enrichment.

We conclude that the coherence in the dynamical and chemical properties of stars in the group suggest that they originate from a single satellite system that was dragged and shredded into the disk of the Galaxy early during its assembly.

## 3 Summary and Discussion

The peculiar angular momentum of the Arcturus group, together with its moderate velocity dispersion and singular metal abundance, all point to a possible extragalactic origin for the group. Age estimates are available for some stars in the group (Fuhrmann 2000, 2003), and suggest that it is perhaps as old as $`10`$-$`12`$ Gyr. This would imply that the group is the debris from an ancient accretion event.

This accretion event has probably contributed a significant number of metal-deficient stars to the solar neighborhood, as shown by the prominence of the group in samples that favor metal-deficient stars. The presence of this dynamically-coherent group may explain some of the variance in the results of dynamical studies of the thick disk. Indeed, some of these studies have concluded that the bulk of the thick disk lags the LSR by only $`30`$$`50\mathrm{km}\mathrm{s}^1`$ (Beers & Sommer-Larsen 1995, Ojha et al 1996, Norris 1999), whereas others have favored a lag of $`100\mathrm{km}\mathrm{s}^1`$ instead (Wyse & Gilmore 1986, Chen et al 2000, Fuhrmann 2000, 2003), possibly reflecting the different contribution of stars in the group to each sample.

There are a number of inquiries that may further validate (or refute) our interpretation. For example, confirming the link between the Arcturus group and the debris identified by Gilmore, Wyse & Norris outside the plane would help to uphold the group’s extragalactic origin. Most of the stars in the Arcturus group are located close to their orbital apocenter, where density enhancements are expected: we may have detected a “shell” similar to those observed around elliptical galaxies (Helmi et al 2003). Other distinct features of this shell interpretation would include a rapid decline in the importance of the group at increasingly large Galactocentric distance (shells are “sharp”) and, of course, the presence of other shells at different Galactocentric radii. It may also be possible to find other members of the progenitor of the Arcturus group by searching for stars on disk-like eccentric orbits with metal abundances consistent with the “sequence” outlined by the group in Figure 3.

Finally, it is oddly gratifying to think of stars visible to the naked eye, such as Arcturus, as silent night-sky witnesses of the merging history of the Milky Way.

We thank Brad Gibson and the hospitality of Swinburne University of Technology, where this work was conceived. JFN is supported by several grants from Canada’s NSERC and from the Canadian Foundation for Innovation.


text:
# Trapping and Steering on Lattice Strings: Virtual Slow Waves, Directional and Non-propagating Excitations

## I Introduction

Can a local excitation (source) in classical field theories be invisible to observers outside the region of excitation? This question has recently received renewed interest.

Berry et al. Berry et al. (1998) described a peculiar excitation case for the one-dimensional wave-equation of a perfectly elastic string under tension. They show that the response of the string can be made to be confined to a bounded region by carefully choosing a forced excitation of oscillatory type. This means that the excitation will not propagate away along the string. Denardo gives a simple and intuitive explanation by using a wave interference argument Denardo (1998). Gbur et al Gbur et al. (1999) discuss conditions of finite string length and dissipation.

Other recent work investigated non-propagating excitations include Marengo and Ziolkowski Marengo and Ziolkowski (1999, 2000a, 2000b) who discuss the generalization of non-propagating conditions of D’Alembertian ($`\mathrm{}\stackrel{\text{def}}{=}^2c^2\frac{^2}{t^2}`$) operators and its temporally reduced version the Helmholtz operator ($`^2+k^2`$) on various related classical scalar and vector fields. Marengo, Devaney and Ziolkowski Marengo et al. (1999) give the condition for time-dependent but not necessarily time-harmonic non-radiating sources and for selective directional radiation for the inhomogeneous wave equation in three spatial dimensions. Marengo and Ziolkowski Marengo and Ziolkowski (2000c) generalize these conditions to more general scalar and vector field dynamics. Marengo, Devaney and Ziolkowski Marengo et al. (2000) also give examples in one and three spatial dimension for the time-harmonic case. Hoenders and Ferwerda Hoenders and Ferwerda (2001) discuss the relationship of non-radiating and radiating parts of the case of the reduced Helmholtz equation, which can be derived from the string equation by assuming general oscillatory time solutions (see Berry et al. (1998)). Denardo and Miller Denardo and Miller (2003) discuss the related case of leakage from an imperfect non-propagating excitation on a string. Gbur Gbur (2003) provides a comprehensive recent review of this topic and the reader is referred to this review for more detailed historical context. Of the earlier work the following contributions are particularly relevent for the discussion here: Schott Schott (1933, 1937) gave the condition for non-radiation of a spherical shell on a circular orbit. Bohm and Weinberg Bohm and Weinstein (1948) extended this result to more general spherical charge distributions and Goedecke Goedecke (1964) showed how an asymmetrical charge distribution with spin is non-radiating. All of these works are concerned with the case of spatially moving sources. Finally it is noteworthy, that non-radiating sources play an important role in inverse problems and have been investigated in a one-dimensional electrodynamic situation by Habashy, Chow and Dudley Habashy et al. (1990).

In this paper our purpose is to describe this phenomena in the case of a lattice string in one dimensions by discretizing D’Alembert’s solution. This approach is used extensively to simulate vibrating strings and air tubes of musical instruments. See Smith (1997) and references therein.

This leads to explicit dynamical constructions of previously reported non-propagating excitations. Its simplicity allows for additional insight into the mechanism that allows for the local confinements and the conditions under which they occur. I will show how the basic mechanisms that provide a time-harmonic stationary non-propagating excitation in one dimension as studied by Berry et al. and Gbur et al. Berry et al. (1998); Gbur et al. (1999) allows for a much wider class of excitations. For instance can such an excitation be relieved from the time-harmonic assumption beyond one period allowing for non-propagating sources that are short-lived. Directional excitations can easily be achieved using very simple bidirectional excitation patterns. These are explicit constructions of such waves in one spatial dimension whose general condition of existence in the three-dimensional case has been derived by Marengo, Devaney and Ziolkowski Marengo et al. (1999). Wave propagation can be virtually slowed down. In general I will show that non-propagating excitations can be extended to steered excitation regions with basic physical restrictions imposed by the underlying field dynamics.

First I will give a quick derivation of the simple lattice model from the wave equation as can also be found in Smith (1997). Then I will give a new argument and construction of the Berry et al. type non-propagating excitation purely based on discrete string dynamics. This will then be compared to the original approach. Then I will extend the discussion to examples of additional types of non-propagating waves, including directional and slowed waves. Finally I will discuss very general constraints on such “steered” localized excitations.

## II Lattice String Model

The lattice string model can easily be derived from the wave-equation by discretizing the D’Alembert solution. Hence the continuous case will be discussed first.

### II.1 Continuous Wave Solutions

The one-dimensional homogeneous wave equation of the perfectly elastic string under tension is:

$$\mu \frac{^2y}{t^2}K\frac{^2y}{x^2}=0$$
(1)
where $`c^2=K/\mu `$ is derived from mass density $`\mu `$ and tension $`K`$. The D’Alembert solution of the homogeneous “free field” case has the well known form (Kreyszig, 1999, p. 596, eq. (4)):

$$y(x,t)=w^+(xct)+w^{}(x+ct)$$
(2)
Hence the solution of the general of the homogeneous wave-equation are two propagating waves whose content is restricted by initial and boundary conditions. As wave-equation is linear we have a connection between initial conditions and external driving forces. Driving forces can be seen as infinitesimal time frames that act on the wave dynamics by imposing a new initial condition at each point in time. Hence we need to consider the initial value problem to gain insight into both processes at once.

At a given time frame $`t_i`$ let the following initial conditions hold:

$`y(x,t_i)`$ $`=f(x,t_i)`$ (3)
$`y_t(x,t_i)`$ $`=g(x,t_i)`$ (4)
Equation (3) with (2) gives a particular solution $`v^+`$:

$`v^+(xct_i)+v^{}(x+ct_i)`$ $`=f(x,t_i)`$ (5)
Taking the first temporal derivative of (2) and satisfying equation (4) we get:

$`cv_t^+(xct_i)+cv_t^{}(x+ct_i)`$ $`=g(x,t_i)`$ (6)
Integrating with respect to $`x`$ we get (Kreyszig, 1999, eq. (10) p. 596):

$`\begin{array}{c}\hfill cv^+(xct_i)+cv^{}(x+ct_i)=k(x_0)+{\displaystyle _{x_0}^x}g(s)𝑑s,\\ \hfill k(x_0)=cv^+(x_0)+cv^{}(x_0)\end{array}`$ (7)
From equations (5) and (7) we can solve for the traveling wave components:

$`\begin{array}{cc}\hfill v^+(xct_i)={\displaystyle \frac{1}{2}}f(x,t_i){\displaystyle \frac{1}{2c}}{\displaystyle _{x_0}^{xct_i}}& g(s)ds\hfill \\ & {\displaystyle \frac{1}{2}}k(x_0)\hfill \end{array}`$ (8)
$`\begin{array}{cc}\hfill v^{}(x+ct_i)={\displaystyle \frac{1}{2}}f(x,t_i)+{\displaystyle \frac{1}{2c}}{\displaystyle _{x_0}^{x+ct_i}}& g(s)ds\hfill \\ & +{\displaystyle \frac{1}{2}}k(x_0)\hfill \end{array}`$ (9)
We see that forced displacement $`f()`$ splits evenly between left and right traveling waves and the integrated forced velocity $`g()`$ splits with a sign inversion.

For our current discussion I will share the assumption of no initial velocity of Berry et al. Berry et al. (1998) and hence the integral over $`g()`$ will vanish.

For the infinite string this is already the complete solution for any twice differentiable function of free solutions and external forced displacements.

### II.2 Discrete Wave Solutions

To arrive at lattice equations we discretize the solution of the wave-equation (2) in time via the substitution $`tTn`$ where $`T`$ is the discrete time-step and $`n`$ is the discrete time index. This automatically corresponds to a discretization in space as well, because in finite time $`T`$ a wave will travel $`X=cT`$ distance according to (2). The spatial index will be called $`m`$. The free-field discrete D’Alembert solution:

$$\begin{array}{cc}\hfill y(mX,nT)=& w^+(mXcnT)+w^{}(mX+cnT)\hfill \end{array}$$
(10)
In general we can always express all discrete equations in terms of finite time steps or finite spatial lengths. We chose a temporal expression and substitute $`X=cT`$ and suppress shared terms in $`cT`$ to arrive at the index version of the discrete D’Alembert solution Smith (1997):

$$y(m,n)=w^+(mn)+w^{}(m+n)$$
(11)
By equations (8) and (9) we see that at an instance $`m_i,n_i`$ the discrete contribution of external forced displacements splits evenly between the traveling waves and we arrive at the discrete field equations including external forced displacements:

$`W^+(m_in_i)=w^+(m_in_i)+{\displaystyle \frac{1}{2}}f(m_i,n_i)`$ (12)
$`W^{}(m_in_i)=w^{}(m_in_i)+{\displaystyle \frac{1}{2}}f(m_i,n_i)`$ (13)
## III Non-propagating Excitation

Next we will construct the non-propagating excitation from the lattice string dynamics directly.

For simplicity and without loss of generality, we will assume a region aligning with the discretization domain throughout. We want to construct an excitation which is confined to a length $`LxL`$. For now we will assume that the string should otherwise stay at rest. This implies that there are no incoming waves into the region $`\mathrm{\Omega }=[L,L]`$ from the outside. We are interested in a non-trivial excitation within the region.

First we consider the contributions to the position $`L`$. As there are no incoming external waves we get:

$`w^+(L+n)=0`$ (14)
We do expect non-trivial wave $`w^{}(Ln)`$ to reach the boundary but we require the total outgoing wave to vanish we have:

$`W^{}(Ln)=w^{}(Ln)+{\displaystyle \frac{1}{2}}f(L,n)=0`$ (15)
The necessary external forced displacement contribution to for cancellation needs to be:

$`{\displaystyle \frac{1}{2}}f(L,n)=w^{}(Ln)`$ (16)
The complete incoming wave (12) will see the same forced contribution (16) and with equation (14) we get:

$`W^+(L+n)={\displaystyle \frac{1}{2}}f(L,n)=w^{}(Ln)`$ (17)
Hence the matched forced displacement leads to a reflection with sign inversion at the region boundary at $`L`$.

Following the same line of argument at point $`L`$ we get the related condition:

$`W^{}(Ln)={\displaystyle \frac{1}{2}}f(L,n)=w^+(L+n)`$ (18)
With these two conditions we can study the permissible form of excitations. First we assume an initial forced displacement impulse from a position $`p`$ in the interior of the domain $`\mathrm{\Omega }\mathrm{\Omega }=(L,L)`$. Hence $`L<p<L`$ and $`f(p,0)=a_p`$ with $`a_p`$.

It will take half the impulse $`L+p`$ steps to reach the left boundary and the other half $`Lp`$ steps to reach the right one.

At each boundary the respective condition (17,18) needs to be satisfied and we get:

$`f(L,L+p)`$ $`=f(p,0)`$ (19)
$`f(L,Lp)`$ $`=f(p,0)`$ (20)
The impulse will then reflect back and create periodic matching conditions.

$`f(L,L+p+4L\nu )`$ $`=f(p,0)`$ (21)
$`f(L,Lp+(2\nu 1)2L)`$ $`=f(p,0)`$ (22)
$`f(L,Lp+4L\nu )`$ $`=f(p,0)`$ (23)
$`f(L,L+p+(2\nu 1)2L)`$ $`=f(p,0)`$ (24)
with $`\nu =1,2,\mathrm{}`$.

Hence we see that a single impulse will necessitate an infinite periodic series of forced external displacements at the boundaries to trap the impulse inside as each “annihilation” of a half-pulse reaching the boundary leads to a “creation” of a reflected one.

The required impulse response of a boundary forced function $`f(\pm L,)`$ can easily be observed from equations (2124) to be spatially periodic in $`4L`$ with an initial phase factor dictated by the starting position $`p`$. Additionally the functional shape of the impulse responses $`f(\pm L,)`$ is completely defined for all time steps as $`f(\pm L,)=0`$ for all times that equations (2124) don’t apply.

A condition for stopping a non-propagating excitation can be derived from the fact that a impulse will return to its initial position every $`4L`$ time steps. Additionally it is easy to see that the traveling impulses will occupy the same spatial position every odd multiple of $`2L`$ with a sign inversion. Hence an impulsive forced displacement $`f((1)^{\mu 1}p,4L\mu )=(1)^{\mu 1}a_p`$ with $`\mu =1,2,\mathrm{}`$ will cancel an initial impulse $`f(p,0)=a_p`$. From this we can immediately deduce the following property:

###### Theorem 1

The shortest possible single impulse finite non-propagating excitation takes $`2L`$ time-steps.

and more generally:

###### Theorem 2

The time of any single impulse excitation finite non-propagating excitation has to be $`2\mu L,\mu `$.

More importantly we observe the property: Non-propagating excitations can be finite in duration.

This is an extension beyond Berry et al. Berry et al. (1998) which assumes infinitely periodic temporal progressions in their derivations.

The general solution for discrete non-propagating wave functions can be derived by observing that any initial “phase” $`p_i`$ is orthogonal to other phases $`p_j`$ for $`i,j\mathrm{\Omega }\mathrm{\Omega }=(L,L)`$, i.e. $`f(p_i,0),f(p_j,0)=0`$ for $`ij`$. Within a $`2L`$ period $`f(\pm L,)`$ is well-defined by $`_if(p_i,)`$. Interestingly though this provides the only restriction to the forced boundary functions. This can be seen by Theorem 1. After $`2L`$ each $`p_i`$ will find constructive interference and can be annihilated or rescaled to an arbitrary other value $`a_i(2L)`$. Hence any arbitrary succession of $`2L2`$ force distributions with a $`2L`$ termination is permissible. Hence periodicity is not necessary.

The time harmonic case can be derived if the initial force distribution within the domain is not modified over time. Then a configuration will repeat after traveling left and right, being reflected at the domain boundary twice, traversing the length of the region twice. Hence the lowest permissible wave-length is $`4L`$. By reflecting twice the wave will have gone through a $`2\pi `$ phase shift, but we note that the periodicity condition is also satisfied if any number of additional $`2\pi `$ shifts have been accumulated. Hence we get for permissible wave-numbers:

$`k={\displaystyle \frac{2\pi n}{4L}},\text{where}n=1,2,\mathrm{}.`$ (25)
or

$`kL={\displaystyle \frac{n\pi }{2}}`$ (26)
By allowing only even $`n`$ we get the Berry at al. condition Berry et al. (1998) for an even square distribution. The odd $`n`$ situation corresponds to the odd-harmonic out-of-phase construction proposed by Denardo Denardo (1998).

Many of these properties can be seen visually in the numerical simulation depicted in Figure 1.

It is interesting to observe that two synchronous point-sources oscillating with the above phase condition will not be completely non-propagating. They will only be non-propagating after waves created at the wave onset have escaped. This is a refinement of the argument put forward by Denardo Denardo (1998) and can intuitively be described as non-interference of the first trap period. Hence the first pairs of pulses will have half-Amplitude components escaping in either direction but every subsequent period will be trapped. This behavior, which could be called imperfect trapping or trapping with transient radiation, is depicted in Figure 2. Sources presented by Berry et al. and Denardo Berry et al. (1998); Denardo (1998) do not display this behavior because the force is assumed to be oscillatory at all times and hence has no onset moment.

Non-propagating excitations can be used as generic building blocks for other unusual excitation induced behavior on the string. In particular I will next describe how to construct an uni-directional emitter, and a virtually slowed propagation. In fact a non-propagating excitation can be seen as virtually stopping a wave at a particular position.

## IV Directional Excitations

A one-sided open trap immediately suggests another unusual excitation type, namely the directional excitation. The string is to be excited in such a way that a traveling wave in only one direction results.

We start with a one-sided open trap. This is a trap that uses a reflection condition (17) and (18) only on one side of an initial excitation. Evidently the wave then can only travel in the opposite direction. For the discussion we will describe a right-sided propagator (i.e. a propagator traveling with increasing negative index). The trapping condition then reads:

$`f(m+1+p,n+p)=f(m,n1)`$ (27)
Hence the trapping excitation point is a $`p`$ time-step lagging negative copy of the original excitation. The emitted wave will have the form

$`{\displaystyle \frac{1}{2}}f(m+1,n+2p){\displaystyle \frac{1}{2}}f(m+1,n)`$ (28)
The emitting wave will show self-interference at a phase of $`2p`$ time-steps, as can be seen in the simulation depicted in Figure 3. In general the self-interference phase can be chosen by the distance $`p`$ between the wave creation point and the trapping point. It is worth noting that it is possible to eliminate interference by trapping the lagging contribution and hence create a wave non-interference directional wave left of the trapping region.

## V Virtual Slow Waves

Virtual slow waves can be achieved by alternating directional wave propagation with trapping. The slowness of the wave propagation can be controlled by the number and and duration times of the traps along a propagation. The propagation characteristics of the dynamic operator has not changed at all, hence we call the this state “virtually slow” as opposed to the case where the field itself induces a change in wave propagation speed. This also means that within a slowed or “steered” region the wave propagation is the one prescribed by the dynamic operator $`(\frac{}{x}+c\frac{}{t})(\frac{}{x}c\frac{}{t})`$ on the string $`y(x,t)`$.

The amount of time spend in traps determines the overall slowness. One example of slow wave consists of an immediate alteration between one stage of trapping and one step of one-sided propagation is illustrated in Figure 4. The effective propagation speed of the wave can easily be read from the diagram to give $`c_{\text{eff}}=c\frac{X}{3T}=c/3`$. As is evident from Theorem 1, a unit $`L=1`$ trap will last $`2`$ time-steps and will not propagate spatially and one step of free propagation will last one time-step and and make one spatial step, hence resulting in a spatial to temporal ratio of $`1:3`$.

The trapping relations are:

$`\begin{array}{cc}\hfill f(m2\nu ,n+1+6\nu )& =f(m+1\nu ,n+6\nu )\hfill \\ & =f(m,n1)\hfill \end{array}`$ (29)
$`\begin{array}{cc}\hfill f(m3\nu ,n+4+6\nu )& =f(m\nu ,n+3+6\nu )\hfill \\ & =f(m,n1)\hfill \end{array}`$ (30)
with $`\nu =0,2,4,\mathrm{}`$.

## VI Steering

The generalized interpretation of the excitation interaction lead to the general dynamical confinement of waves by external excitation. For instance following very similar arguments as for virtual slow waves a construction is possible which gives a slowed “cone of influence” by successively widening the trap boundaries at a speed slower than the the wave speed $`c`$. By this argument it is sufficient for the trap boundaries’ change to be less than $`c`$ for it to be trapping the wave. This is not a necessary condition by the following counter-example: Let the trap width be $`L`$ and change rapidly by some slope $`dL>c`$ to some new constant width $`L_2`$ at which it becomes constant. Obviously the wave will then be able to reach the new boundary even though a local change of the boundary exceeded the dynamical speed $`c`$. The necessary condition can be seen from our previous construction. At a trap boundary a wave is reflected and will propagate in the opposite direction of the domain following the linear characteristic $`c`$. Only if this characteristics intersects with the dynamic trapping boundary will there be another externally forced reflection as illustrated in Figure 5. These may in fact have regions where no trapping is necessary and possible.

## VII Interaction with Background Fields

It is important to note that while we assumed that the incoming wave vanishes, see equation (14), the outgoing wave condition (15) does not change if there is in fact an incoming wave. The “reflection wave” (17) and (18) can be rewritten for a non-zero incoming field without affecting the trapping:

$`\begin{array}{cc}\hfill W^+(L+n)& =w^+(L+n)+{\displaystyle \frac{1}{2}}f(L,n)\hfill \\ \hfill {\displaystyle \frac{1}{2}}f(L,n)& =w^{}(Ln)\hfill \end{array}`$ (31)
and

$`\begin{array}{cc}\hfill W^{}(Ln)& =w^{}(Ln)+{\displaystyle \frac{1}{2}}f(L,n)\hfill \\ \hfill {\displaystyle \frac{1}{2}}f(L,n)& =w^+(Ln)\hfill \end{array}`$ (32)
These conditions are “absorbing” in the sense that an external field entering the trapping region will not leave it.

The “non-interacting” property of a trap defined by the periodic matching conditions (2124) can be seen by assuming a non-zero incoming wave at one point of the trap boundary $`\delta \mathrm{\Omega }`$. Then the total wave entering the trapping region the sum of the wave created by the trapping condition and the incoming wave value $`\frac{1}{2}f(\delta \mathrm{\Omega }^1,)+w^\pm (\delta \mathrm{\Omega }^1,)`$, where $`\delta \mathrm{\Omega }^1`$ denotes the first trap boundary reached. When reaching the second trapping boundary $`\delta \mathrm{\Omega }^2`$ the now outgoing wave will see a matching force $`f(\delta \mathrm{\Omega }^2,)=\frac{1}{2}f(\delta \mathrm{\Omega }^1,)`$ leaving an outgoing wave contribution $`w^\pm (\delta \mathrm{\Omega }^2,)=w^\pm (\delta \mathrm{\Omega }^1,)`$ to escape the trapping region $`\mathrm{\Omega }`$.

In order to achieve selective radiation, only part of the content of a trapped region are trapped at the boundary as can be achieved by using a reduced force at the trapping boundary or by selectively omitting certain phases in the trapping force pattern.

### VII.1 Relationship of Traps to Non-Radiating Sources

Marengo and Ziolkowski Marengo and Ziolkowski (1999) present ideas very much related to ideas presented here and in Berry et al.Berry et al. (1998).

However, they arrive at a definition of non-radiating (NR) sources that is not obviously similar to the traps presented here. In particular they define NR sources as being non-interacting. While Marengo and Ziolkowski (1999) note that a central property of NR sources is that they store non-trivial field energy, traps described here can not only store, but accumulate and selectively radiate waves.

The difference can be understood by observing that for example Berry et al. assume a simple time-harmonic driver (Berry et al., 1998, eq. (3)) throughout their discussion:

$`f(x,t)=\text{Re}\left\{f(x)e^{i\omega t}\right\}`$ (33)
By our earlier discussion we see that the temporal progression of the boundary has to match the content of the interior domain. Hence once the boundary is defined to be oscillatory the interior of the domain needs to be spatially harmonic as derived in Berry et al. (1998); Marengo and Ziolkowski (1999) and has been rederived here. Hence a NR source as noted in literature, with exception of the general orthogonality formulation for time-varying sources given by Marengo, Devaney and Ziolkowski Marengo et al. (1999), can be thought of as a time-oscillatory trap.

The arguments made here use a formalism that is discrete in nature. However, the discreteness of the arguments are not necessarily restrictive. The continuous case can be imagined with the discrete time-step made small ($`T0`$) or alternatively, discrete pulses can be substituted with narrow distributions of compact support. In neither case are the results of interest derived here altered.

As has already been derived in Berry et al. (1998); Gbur et al. (1999) the critical condition for non-propagating waves lie at the boundary of the domain-range that the wave ought not to leave. In the discrete case it is easy to see how this insight can be used and generalized. In fact, the boundaries of the confining domain need not be static, nor need the condition be used in a two-sided fashion.

## VIII Conclusion

In summary, this paper presented constructions of a broad class of non-propagating sources on a string lattice model using trapping conditions. In particular this includes numerical demonstrations of finite-duration non-propagating excitations, directional excitations, as well as virtually slowed waves. These examples help explain the extension of non-propagating sources beyond the time-periodic case and include treatment of onset, annihilation and spatial steering. These properties ought to be observable in experiments well-described by the wave equation. This equation often arises in problems in acoustics, elasticity, optics and electromagneticsm. And hence the results presented here apply to these domains of application. While here I discussed the forward problem, these results also relate to the inverse problem of finding source contributions from the one-dimensional field state as occur for example in acoustical, optical and electromagnetic detection problems.

###### Acknowledgements.

I am grateful for reprints provided by Bruce Denardo and Greg Gbur, who also brought relevant references to my attention. I’m also thankful for helpful comments and pointers to relevant literature by an anonymous referee.