format: text

text:
# Discovery and Photometric Observation of the Optical Counterpart in a Possible Galactic Halo X-ray Transient, XTE J1118+480

## 1. Introduction

Soft X-ray transients (SXTs), or X-ray novae (e.g. Chen et al. 1997) are binary systems which exhibit luminous X-ray and optical outburst (Tanaka, Shibazaki 1996). They uniquely provide the most compelling evidence for the existence of steller mass black holes using radial velocity studies, giving mass functions exceeding the maximum mass of a stable neutron star ($`3M_{}`$), and we know eight such black hole candidates (van Paradijs, McClintock 1995; Bailyn et al. 1998; Orosz et al. 1998). Their outburst light curves often have an common feature, that is, after a rise of a few days, the X-ray intensity comes to the maximum, which was typically followed by the exponential decay with an e-folding time of $`40\mathrm{d}`$ (Chen et al. 1997). At the maximum, the X-ray luminosity reaches $`10^{3839}\mathrm{erg}\mathrm{s}^1`$ and the optical to X-ray flux ratio of $`500`$ is typical in SXTs (Tanaka, Shibazaki 1996).

Many astronomers have recently believed that the compact object of SXTs in quiescence is surrounded by the accretion disk whose inner part is the geometrically thick and optically thin advection dominated accretion flow (ADAF) and the outer part is the geometrically thin and optically thick disk which becomes thermally unstable when the disk becomes too hot for hydrogen to remain neutral (Narayan, Yi 1995; Shakura, Sunyaev 1973; Osaki 1974). This model of the outburst mechanism is called the disk instability model and satisfactorily explains the outburst cycle of a few tens of years and the outburst duration by the viscous diffusion time scale of the accretion disk in the cool and hot states, respectively (Mineshige 1996).

Almost all observed SXTs are distributed on the galactic disk and few had been discovered in the galactic halo (Chen et al. 1997; Bradt et al. 2000). White, van Paradijs (1996) have studied the galactic distribution of BHCs low-mass X-ray binaries and found an rms value for the distance from the galactic disk to be $`0.4`$ kpc. We can easily understand these biased distribution because massive stars, which are responsible for producing neutron stars and black holes, have been generated in the disk rather than the halo.

A new SXT, XTE J1118+480 whose galactic latitude is high ($`62^{}`$), has been discovered at an intensity of 39 mCrab with the All-Sky Monitor (ASM; Levine et al. 1996) on the Rossi X-Ray Timing Explorer (RXTE) in 2000 March 29 (Remillard et al. 2000). The X-ray spectrum just after the discovery is similar to Cyg X-1 in its hard state (Remillard et al. 2000). The hard X-ray spectrum is well characterized by a power law with a photon index of 2.1 and the source is visible up to 120 keV which implies that it is a possible black hole X-ray transient (Wilson, McCollough 2000). In this paper, we report the discovery of the optical counterpart of XTE J1118+480 and the optical short-time variability, and discuss the binary nature and a distance estimation.

## 2. Discovery and Observations

Following communication of the X-ray discovery of XTE J1118+480 to the Variable Star NETwork (VSNET, http://www.kusastro.kyoto-u.ac.jp/vsnet) on 2000 March 30, we started CCD time-series observations, which immediately revealed the presence of the optical counterpart at 12.92 mag, within the error circle of $`6^{}`$ proposed by RXTE/ASM (Uemura et al. 2000a). The position of the optical counterpart is R.A. = $`11^\mathrm{h}\mathrm{\hspace{0.17em}18}^\mathrm{m}\mathrm{\hspace{0.17em}10}^\mathrm{s}.85`$, Decl. = $`+48^{}\mathrm{\hspace{0.17em}02}^{}\mathrm{\hspace{0.17em}12}^{\prime \prime }.9`$ (equinox 2000.0; accuracy $`0.2^{\prime \prime }`$). A star of 18.8 mag in the USNO A1.0 and A2.0 catalogues (USNO A1.0 1350.08089912 = USNO A2.0 1350.07924726) lies within $`2^{\prime \prime }`$ of this position. No other object can be seen on images of the second DSS images available from the Space Telescope Science Institute within $`10^{\prime \prime }`$, indicating that this star is the quiescent optical counterpart.

Figure 1 gives the finding chart of XTE J1118+480. The left and right panels show the DSS image and outburst image on March 30.83 (UT), respectively. After detection of March outburst in optical and X-ray wavelengths, the reanalysis of previous data revealed the presence of another X-ray outburst and simultaneous optical activity in 2000 January (Remillard et al. 2000; Uemura et al. 2000a; Wren, McKay 2000).

The CCD photometric observations were made using the following cameras and telescopes; unfiltered ST-7 and 25-cm Schmidt-Cassegrain (Kyoto), unfiltered SXL8 and 33-cm Newtonian (Conder Brow Observatory), CB245 with clear filter and 44-cm Newtonian (California), unfiltered ST-7 and 28-cm Schmidt-Cassegrain (Ceccano), and unfiltered ST-9e and 28-cm Schmidt-Cassegrain (Nayoro), whose integration times are $`30\mathrm{s}`$ , $`30\mathrm{s}`$, $`16\mathrm{s}`$, $`40\mathrm{s}`$, and $`25\mathrm{s}`$, respectively. After correcting for the standard de-biasing and flat fielding, we processed object frames with the PSF and aperture photometry packages. We performed differential photometry relative to the comparison star, C, shown in figure 1, whose constancy was confirmed by check star of GSC3451.938. Photographic observations were performed using an unfiltered T-Max400 film with 10-cm camera at Nagano and Aichi.

## 3. Results

We observed no activity brighter than 15 mag from HJD 2449660 (1994 November) to 2451547 (1999 December) with 17 negative observations. Figure 2 shows the light curve of XTE J1118+480 after it became active in 2000 January. The abscissa denotes time in heliocentric julian date. The filled circles, filled diamonds, open circles, and open triangle denote the photographic magnitude, nightly averaged unfiltered CCD magnitude which almost equals to $`R_c`$-magnitude, the $`V`$-magnitude, and upper-limit reported by ROTSE collaboration, respectively (Uemura et al. 2000a; Wren, McKay 2000). The error of each photometric observation is typically less than 0.1 mag. We also depict the X-ray flux observed with RXTE/ASM shown in gray crosses with errorbars (http://xte.mit.edu/ASM\_lc.html). The discovery date corresponds to the first of the filled diamonds on HJD 2451634.

In figure 2, we can see two distinct X-ray outbursts reported by Remillard et al. (2000) and two optical outbursts which correlate with X-ray intensity. At the beginning of the second outburst, figure 2 clearly shows the optical precursor which began to rise about 10 days prior to the X-ray outburst. We also see a fast optical flare around HJD 2451572 which anti-correlates with the X-ray intensity. The flare occurred just before the first X-ray outburst ended and lasted only for a few days. A possible optical flare is also seen around HJD 2451627.

Figure 3 shows the light curve of XTE J1118+480 from HJD 2451634 to 2451647. In figure 3, we can see the sinusoidal periodic modulation upon the plateau of 12.89 mag and no general tendency of rise nor decay was observed on the limited data available. The data on HJD 2451643 were special in that they showed a rise over $`0.2\mathrm{d}`$, which was not seen in the other data. Although we can see a possible hump feature during this rise, we have excluded this segment from the following period analysis. We performed a period analysis using the PDM method (Stellingwerf 1978).

The resultant period-theta diagram is shown in figure 4. As can be seen in figure 4, the best candidate of the period is $`0.17078\pm 0.00004\mathrm{d}`$. This period is also confirmed with the subsequent observations by VSNET collaboration team (Uemura et al. 2000b). We see a remarkable dip around HJD 2451643.0 while the all other data well fits with the single peak averaged light curve. The amplitude of 0.055 mag is not changed during our time-series CCD observations. We noted that the scatter seemed large for a star of this brightness even if subtracting the periodic modulation, indicating the presence of fluctuations on a short (seconds) time scale.

## 4. Discussion

XTE J1118+480 shows some features different from those of typical X-ray transients, that is, the multiple peaks within three months observed both in X-ray and optical, an extremely low optical to X-ray flux ratio, and the high galactic latitude of $`62^{}`$. For the first outburst shown in figure 3, the light curve profile around the peak is relatively sharp compared with the second one and no clear time lag is seen between the X-ray and optical outburst. We can interpret this correlation through the occurrence of an “inside-out” type outburst in accretion disk which will generate the rapid rise in X-ray and simultaneous optical brightening, however this picture cannot predict the fast optical flare at the end of outburst (Smak 1984). On the other hand, we suggest the “outside-in” type outburst for the second outburst to explain the optical precursor. To observe these different types of outburst indicates that they may be caused by the different instability faculty such as normal and super outbursts of SU UMa-type dwarf novae (Warner 1995).

From the X-ray flux of 39 mCrab reported in Remillard et al. (2000) and the simultaneous magnitude of optical counterpart of $``$ 13 mag, the optical to X-ray flux ratio is calculated as $`5`$ which is by two orders of magnitude smaller than the typical value of $`500`$. The X-ray intensity particularly shows a low value rather than optical, and this suggests an idea that it may be a high inclination system (Garcia et al. 2000). The HST observation revealed the continuum slope somewhat flatter, which suggests the intrinsic X-ray flux is relatively low and hence removes the suggestion of the high inclination system (Haswell et al. 2000). A high inclination system generally provides a high probability of occurrence of eclipses by secondary, whereas we have not detected them.

It is almost certain that the periodicity appeared in figure 3 reflects the orbital motion of the underlying binary, through the reprocessing the X-rays by the secondary, or possibly superhump modulations as seen in SU UMa-type dwarf novae (O’Donoghue, Charles 1996). We therefore suggest the orbital period of $`0.17078`$ d, which is relatively short compared with other low mass X-ray binaries (Ritter, Kolb 1998).

Since the quiescent magnitude is 18.8 mag, if we assume a K dwarf secondary star and no other optical source in quiescence, the distance is at least 1.5 kpc, that is, 1.3 kpc above the galactic plane. The color of the quiescent counterpart in the USNO catalog ($`br=+0.6`$) is much bluer than that of a K or M dwarf, implying the substatial contribution from the accretion disk. The above distance estimate should be thus regarded as a lower limit, implying that XTE J1118+480 is a galactic halo X-ray transient. We can generally expect that such a galactic halo object is older than that in disk, and this is consistent with the short orbital period of 0.17078 d which implies a more evolved binary system. It should be noted, however, if we assume a short orbital period of 0.17078 d, the secondary may be a smaller star than in other typical SXTs. This case may be comparable to the system GRO J0422+32, whose orbital parameter and secondary type are 0.212 d and M2V, respectively (Orosz, Bailyn 1995; Filippenko et al. 1995). An M2V star secondary leads to the lower limit of the distance of 0.55 kpc. Optical or infrared spectroscopic observations are definitely needed to unambiguously determine the orbital period and the type of the secondary.

On the point of the high galactic latitude SXTs, 3U 0042+32 is similar to XTE J1118+480 (Ricketts, Cooke 1977). The outburst in 1977 of 3U 0042+32 was detected by Uhuru which classified it as a high galactic latitude ($`30^{}`$) X-ray transient. During the active phase in 1977 Feburuary, 3U 0042+32 experienced at least four distinct outbursts whose interval is estimated as 11.6 d (Watson, Ricketts 1978). Compared with XTE J1118+480, the outburst duration seems to be quite shorter (decay e-folding time $`=2.8`$ d). Another noteworthy, possibly related, high galactic latitude X-ray binary is MS 1603+2600 = UW CrB (Morris et al. 1990). This object is a persistent source whose orbital period is 111 min. Although Hakala et al. (1998) suggested this X-ray binary may correspond to a quiescent state of SXTs, the discovery of a short-period transient system XTE J1118+480 may suggest that these two systems comprise a new variety of halo X-ray binaries.

## 5. Summary

We discovered the optical counterpart of the high galactic latitude soft X-ray transient, XTE J1118+480, whose position is R.A. = $`11^\mathrm{h}\mathrm{\hspace{0.17em}18}^\mathrm{m}\mathrm{\hspace{0.17em}10}^\mathrm{s}.85`$, Decl. = $`+48^{}\mathrm{\hspace{0.17em}02}^{}\mathrm{\hspace{0.17em}12}^{\prime \prime }.9`$. There are 18.8 mag star within the $`2^{\prime \prime }`$ error circle in USNO catalogue and we consider it as an optical counterpart in quiescence. We revealed the two distinct optical outbursts before discovery. From our time-series CCD photometry, we detected a 0.17078 d periodicity which we consider as the orbital period. If we assume the secondary star of K dwarf, XTE J1118+480 is the possible first firmly identified black hole candidate X-ray transient in the galactic halo.

## References

Bailyn C. D., Jain R. K., Coppi P., Orosz J. A. 1998, ApJ 499, 367

Bradt H., Levine A., Remillard R., Donald A. S. astro-ph/0003438

Chen W., Shrader C. R., Livio M. 1997, ApJ 491, 312

Filippenko A. V., Matheson T., Ho L. C. 1995, ApJ 455, 614

Garcia M., Brown W., Pahre M., McClintock J. 2000, IAU Circ. 7392

Hakala P. J., Chaytor D. H., Vihu O., Piirola V., Morris S. L., Muhli P. 1998, ApJ 333, 540

Haswel C. A., Hynes R. I., King A. R. 2000, IAU Circ. 7407

Levine A. M., Bradt H., Cui W., Jernigan J. G., Morgan E. H., Remillard R., Shirey R. E., Smith D. A. 1996, ApJ 469, L33

Mineshige S. 1996, PASJ 48, 93

Morris S. L., Liebert J., Stocke J. T., Gioia I. M., Schild R. E., Wolter A. 1990, ApJ 365, 686

Narayan R., Yi I. 1995, ApJ 452, 710

O’Donoghue D., Charles P. A. 1996, MNRAS. 282, 191

Orosz J. A., Bailyn C. D. 1995 ApJ 446, L59

Orosz J. A., Jain R. K., Bailyn C. D., McClintock J. E., Remillard R. A. 1998, ApJ 499, 375

Osaki Y. 1974, PASJ 26, 429

Remillard R., Morgan E., Smith D., Smith E. 2000, IAU Circ. 7389

Ricketts M. J., Cooke B. A. 1977, IAU Circ. 3039

Ritter H., Kolb U. 1998, A&AS 129, 83

Shakura N. I., Sunyaev R. A. 1973, A&A 24, 337

Smak J. 1984, Acta Astron. 34, 161

Stellingwerf R. F. 1978, ApJ 224, 953

Tanaka Y., Shibazaki N. 1996, ARA&A 34, 607

Uemura M., Kato T., Yamaoka H. 2000a, IAU Circ. 7390

Uemura M., Kato T., Matsumoto K., Honkawa M., Cook L. M., Martin B., Masi G., Oksanen A., Moilanen M., Novak R., Sano Y., Ueda Y. 2000b, IAU Circ. 7418

van Paradijs J., McClintock J. E. 1995, in X-ray Binaries, ed. W.H.G. Lewin J., van Paradijs E. P. J. van den Heuvel (Cambridge: Cambridge Univ. Press), 58

White N., van Paradijs J. 1996, ApJ 473L, 25

Warner B. 1995, Cataclysmic Variable Stars, p126–215 (Cambridge Univ. Press, Cambridge)

Watson M. G., Ricketts M. J. 1978, MNRAS. 183, 35P

Wilson C. A., McClollough M. L. 2000, IAU Circ. 7390

Wren J., McKay T. 2000, IAU Circ. 7394


text:
# 1 Introduction

## 1 Introduction

Astrophysical systems in the Universe are characterized by the gravitation. The structure formed through this long-range force is quite different from those formed through the other short-range forces. If the system does not strongly depend on the initial conditions of the Universe, we can apply statistical mechanics for describing such self-gravitating systems (SGS). However, we cannot directly apply the standard Boltzmann statistical mechanics for SGS since the long-range nature of gravity strongly violates the extensive property of the system which is the premise of statistical mechanics. Actually, the total energy increases much faster than the particle number $`N`$, the partition function $`Z`$ often becomes complex, reflecting the fact that there is no absolute stable state in SGS.

In order to seek for workable statistical mechanics of SGS, we try an approach based on the new definition of entropy whose extensivity is violated from the beginning; Tsallis statistical mechanics.

We formulate the count-in-cell method for the large scale galaxy distributions in this new statistical mechanics. First we calculate the expression of the composite entropy and the generalized Euler relation in this new statistical mechanics. These are applied to the data of CfA II South redshift survey. The parameter q becomes negative, which represents the instability of gravity.

## 2 Tsallis Statistical Mechanics

The ordinary Boltzmann statistical mechanics is characterized by the entropy $`S=\underset{i}{}p_i\mathrm{ln}p_i`$. The distribution function $`p_{N,E}=\mathrm{exp}\left[\left(E\mu N\right)/T\right]/\mathrm{\Xi }`$ maximizes this entropy with constraints of the probability conservation, the energy conservation, and the particle number conservation. This statistical mechanics is originally aimed to describe the multi-fractal and chaos structures. It is characterized by the entropy of the form  $`S_q\left[p\right]=\left(\underset{i}{}p_i^q1\right)/(1q)`$, where $`q`$ is a real parameter. Tsallis distribution function is obtained so that it maximizes this entropy with the same constraints. The solution has a power law tail:

$$p_{N,E}=\frac{1}{\mathrm{\Xi }_q}\left\{1\frac{1q}{\stackrel{~}{T}}\left(E\overline{E}\mu \left(N\overline{N}\right)\right)\right\}^{\frac{1}{1q}},$$
(1)
where $`\stackrel{~}{T}CT`$, and $`C=\underset{i}{}p_i^q`$. Partition function is defined as

$$\mathrm{\Xi }_q=\underset{N,E}{}\left\{1\frac{1q}{\stackrel{~}{T}}\left(E\overline{E}\mu \left(N\overline{N}\right)\right)\right\}^{\frac{1}{1q}}.$$
(2)
Note that $`p_{N,E}`$ reduces to the ordinary Boltzmann form for $`q1`$. For the consistent formulation, it is important to notice that the observable expectation value is calculated by the escort distribution: $`P_i=p_i^q/C`$, $`Q_q=\underset{i}{}P_iQ_i`$. The averaged quantities $`\overline{E},\overline{N}`$ in Eqs.(1) and (2) should be understood in this sense.

## 3 Galaxy Distribution in Count-in-Cell Method

There are many approaches to describe large scale structure of the Universe. In this paper, we will concentrate on distribution of galaxies. One of the well-known method to describe the distribution of galaxies is the two-point correlation function. However, when we solve equations of motions for the two-point correlation function, we need three-point correlation function. Such higher order essentialness is called BBGKY chain. Since we have to cut the BBGKY chain, we have to apply some kind of approximations.

On the other hand, count-in-cell method is often used to describe distribution of galaxies. In this method, we use an analytic formula for probability $`f(N)`$ of finding $`N`$ galaxies in a randomly positioned volume $`V`$. As we don’t have to use any approximation in count-in-cell method, we can apply it even for clustered system. Saslaw and Hamilton, and S. Inagaki introduced the virial parameter b=(gravitational correlation energy)/(kinetic energy of random motion), which measures the deviation from the dynamical-equilibrium. Then they found that in thermal-equilibrium, their theoretical investigation of $`f(N)`$ is consistent with he N-body simulations or catalogues of observations. Strictly speaking, we should not apply thermal-equilibrium theory for expanding Universe. However their consistency let us further study thermal-equilibrium statistical description. In evolution of the Universe, dynamical-equilibrium $`b=1`$ must be realized before the Universe reach thermal-equilibrium. Therefore in this paper we consider the dynamical-equilibrium $`b=1`$ case. We believe that there exist adequate fitting parameter other than the virial parameter $`b`$. That is the reason why we consider non-extensive statistics which gives us parameter $`q`$.

Supposing the galaxies distribute according to Tsallis statistics, we consider a system described by equilibrium thermodynamics. That is the grand canonical ensemble, characterized by the given temperature and the chemical potential.

The probability to find no galaxy in the volume $`V`$ is

$$f\left(0\right)\underset{E}{}P_{0,E}=P_{0,0}=\frac{\left(p_{0,0}\right)^q}{C},$$
(3)
where

$$p_{0,0}=\mathrm{\Xi }_q^1\left[1+\frac{1q}{\stackrel{~}{T}}\left(\overline{E}\mu \overline{N}\right)\right]^{\frac{1}{1q}}.$$
(4)
The partition function can be decomposed as

$$\mathrm{\Xi }_q=\underset{N,E}{}\left(1\frac{1q}{\stackrel{~}{T}}\left(E\overline{E}\mu \left(N\overline{N}\right)\right)\right)^{1+q/\left(1q\right)}=\left(\mathrm{\Xi }_q\right)^qC$$
(5)
and therefore $`\mathrm{\Xi }_q=C^{1/\left(1q\right)}=\left(1+\left(1q\right)S\right)^{1/\left(1q\right)}`$. Thus we obtain

$$f\left(0\right)=\left(1+\left(1q\right)S\right)^{\frac{1}{1q}}\left(1+\frac{1q}{\stackrel{~}{T}}\left(\overline{E}\mu \overline{N}\right)\right)^{\frac{q}{1q}}.$$
(6)
In order to reduce the above expression, we need to calculate the composite entropy and the Euler relation in Tsallis statistics. When we compose two systems A and B, the distribution function is given by $`p_{ij}(A,B)=p_i(A)p_j(B)`$ and the composed entropy is $`S_{A+B}=S_A+S_B+\left(1q\right)S_AS_B.`$ Sequentially using this composition $`S_{N+1}=s+S_N+\left(1q\right)sS_N`$, (where $`s`$ means entropy of 1 system), we obtain the total entropy of N identical systems as

$$S_N=\frac{\left(1+\left(1q\right)s\right)^N1}{1q}.$$
(7)
Generalized Euler relation is given by the following arguments. The variables $`E,V,N`$ are the natural arguments of the entropy: $`S_N=S(E,V,N)`$. Differenciating $`S_{\alpha N}=S(\alpha E,\alpha V,\alpha N)`$ with respect to $`\alpha `$ and setting $`\alpha 1`$, we obtain

$$\frac{S_{\alpha N}}{\alpha }|_{\alpha 1}=\frac{S(\alpha E,\alpha V,\alpha N)}{\alpha }|_{\alpha 1}=\frac{E}{T}+\frac{pV}{T}\frac{\mu N}{T},$$
(8)
where we used $`S/V=p/T,S/E=1/T,S/N=\mu /T`$, which are guaranteed by the Legendre structure of the Tsallis thermodynamics. In the above, we can put any other values for $`\alpha `$ as well. In general, the $`\alpha `$-dependence is inherited from $`S`$ to $`T`$ on the right hand side of Eq.(8). Non-extensivity of $`S`$ necessarily accompanies the non-intensivity of the Legendre-conjugate variable $`T`$. Actually, the fact that the variables $`E,V,N`$ are extensive and $`p,\mu `$ are intensive in Eq.(8) renders the $`\alpha `$-dependence of the temperature:

$$T\left(\alpha \right)=T_1\left(1+\left(1q\right)s\right)^{1\alpha N}$$
(9)
where $`T_1`$ is the temperature of one particle system. Note that $`\stackrel{~}{T}CT`$ defined just after Eq.(1) is $`\alpha `$-independent! Probably this quantity $`\stackrel{~}{T}`$ should be related with the physical temperature as defined from the velocity dispersion of the system. However at present, we do not have idea on the meaning of the quantity $`T`$111 This situation is similar to the argument on the physical probability distributions at the end of the section two. Among $`P_i`$ and $`p_i`$, related with each other by the factor $`C`$, we have chosen the escort distribution $`P_i`$ as the physical distribution. . This point will be further discussed in the last section.

Then we obtain the Euler relation in non-extensive statistical mechanics

$$\frac{\left(1+\left(1q\right)S\right)\mathrm{ln}\left[1+\left(1q\right)S\right]}{1q}=\frac{E+pV\mu N}{T}.$$
(10)
This temperature on the right hand side is $`T\left(\alpha =1\right)`$ in Eq.(9), though the explicit form of the temperature does not appear in the final expression of $`f\left(N\right)`$. Using this Euler relation, we can now eliminate $`\overline{E}\mu \overline{N}`$ in $`f\left(0\right)`$ and Eq.(6) reduces to

$$f\left(0\right)=\left(1+\left(1q\right)S\right)^{\frac{1}{1q}}\left(1\frac{pV}{\stackrel{~}{T}}\left(1q\right)+\mathrm{ln}\left[\left(1+\left(1q\right)S\right)\right]\right)^{\frac{q}{1q}}.$$
(11)
Since we consider the dynamical-equilibrium system, the system is fully virialized $`b=1`$, and therefore the pressure $`p`$ must be $`0`$, we obtain:

$$f\left(0\right)=\left(1+\left(1q\right)S\right)^{\frac{1}{1q}}\left(1+\mathrm{ln}\left[\left(1+\left(1q\right)S\right)\right]\right)^{\frac{q}{1q}}.$$
(12)
Moreover using Eq.(7), we finally obtain

$$f(0)=\left\{1+\left(1q\right)s\right\}^{\frac{N}{1q}}\left[1+N\mathrm{ln}\left\{1+\left(1q\right)s\right\}\right]^{\frac{q}{1q}}.$$
(13)
Our parameters are $`s`$ (the unit of non-extensive entropy per galaxy) and $`q`$ (the Tsallis statistical parameter).

Probability of finding $`N`$ galaxies in the volume V is generated from $`f\left(0\right)`$:222 Note that this expression is slightly different from that by Saslow and Hamilton .

$$f\left(N\right)\frac{(n)^N}{N!}\frac{d^N}{dn^N}f(0)$$
(14)
where $`n`$ is the galaxy number density. This is because the void probability $`f\left(0\right)`$ contains all the information of the whole correlation functions:

$$f\left(0\right)=1P\left\{X_1\right\}+P\left\{X_1X_2\right\}P\left\{X_1X_2X_3\right\}\pm \mathrm{},$$
(15)
where $`P\left\{X_1X_2\right\}`$ is the probability that there are galaxies in $`dV_1`$ at $`X_1`$ and $`dV_2`$ at $`X_2`$. The above expression guarantees that the probability is properly normalized: $`\underset{N=0}{\overset{\mathrm{}}{}}f\left(N\right)=1`$.

## 4 Comparison with observations

We use the data of CfA II South redshift observations which includes 4392 galaxies . We have to reduce the data to the uniform sample. First we restrict the data to the galaxies whose absolute luminosity is brighter than the magnitude -19.1 and the distance $`4,0008,000[\mathrm{km}/\mathrm{sec}]`$. The distance is measured by the cosmic redshift. We further exclude the edge of the observation region. We applied the K-correction for compensating the reddening. Finally the data is reduced to 870 galaxies. The number density is $`n=4.44\times 10^9[\left(\mathrm{km}/\mathrm{sec}\right)^3]`$.

We first fit the void probability $`f\left(0\right)`$ by varying the parameters q and s. The best fit is realized by $`q=5.66847`$ and $`s=0.164142`$ (Fig.1). We fix these values and do not change them hereafter in this paper. With these parameters, general probability $`f\left(N\right)`$ is given. In Figs.1-2, we compared our calculations and the CfA data.

We have checked the normalization of probabilities $`\underset{N=0}{\overset{\mathrm{}}{}}f\left(N\right)=1`$ is realized.

The negative value of $`q`$ we obtained may not be so surprising. In reference , the value of $`q`$ appears to decrease monotonically from 1 to $`\mathrm{}`$ for one-dimensional logistic maps. Their model reveals unstable onset-to-chaos attractor. For negative $`q`$, the entropy functional loses its convexity and the distribution becomes unstable. Thus the intrinsic instability of SGS is faithfully represented in this formalism and this fit.

## 5 Conclusions and Discussions

We have constructed the non-extensive statistical mechanics based on the non-extensive entropy. Especially calculating the entropy of composite systems and deriving the generalized Euler relation in thermodynamics, we could evaluate the void probability function $`f\left(0\right)`$ and the probability to find $`N`$ galaxies $`f\left(N\right)`$. This result was applied to the CfAII South galaxy observations and we have obtained negative parameter $`q`$. This is thought to be another representation of the intrinsic instability of SGS. It will be also interesting to notice the fact that the multi-fractal scaling is observed in this CfAII South data within the scale-region from $`500[\mathrm{km}/\mathrm{sec}]`$ to $`3000[\mathrm{km}/\mathrm{sec}]`$. We would like to clarify possible connection between the non-extensive distributions and the multi-fractal nature in the context of gravity.

On the way we derive $`f\left(N\right)`$, we encountered “scale ($`\alpha `$) dependent temperature $`T`$”. If we put the values we obtained $`q=5.66847`$ and $`s=0.164142`$ to Eq.(9), we can explicitly plot the scale dependence of the temperature. It turns out to reduce with increasing scale and abruptly drops to zero at about $`r600[\mathrm{km}/\mathrm{sec}]`$ or, assuming the cosmic expansion speed as $`H=72[\mathrm{km}/\mathrm{sec}/\mathrm{Mpc}]`$, at about $`8.3[\mathrm{Mpc}]`$, which is almost the scale that the galaxy correlation function becomes unity. On the other hand, it is apparent that the galaxies do have peculiar velocity of order $`1000[\mathrm{km}/\mathrm{sec}]`$ at this scale. Therefore, at least, the quantity $`T`$ cannot be interpreted as the ordinary temperature as defined from the velocity dispersions. One of our next task will be to elucidate the meaning of $`T`$.

In relation with the astronomical velocity distributions, the authors claim that the velocity distribution of the clusters of galaxies can be well fitted by the Tsallis distribution with the parameter $`q=0.23_{0.05}^{+0.07}`$, which is apparently different from our negative value for galaxy distributions. Further study on the velocity distributions in various scales (galaxies, clusters, super-clusters) would reveal the origin and evolution of the large-scale structure of Universe. We would like to report these results in the near future.

## Acknowledgment

I. J. would like to thank Prof. T. Yokobori and Prof. A. T. Yokobori, Jr for their hospitality. All of us would like to thank S. Abe, M. Hotta and K. Sasaki for useful discussions and valuable comments.


text:
# Roughness of fracture surfaces

## Abstract

We study the roughness of fracture surfaces of three dimensional samples through numerical simulations of a model for quasi-static cracks known as Born Model. We find for the roughness exponent a value $`\zeta 0.5`$ measured for “small length scales” in microfracturing experiments. Our simulations confirm that at small length scales the fracture can be considered as quasi-static. The isotropy of the roughness exponent on the crack surface is also showed. Finally, considering the crack front, we compute the roughness exponents of longitudinal and transverse fluctuations of the crack line ($`\zeta _{}\zeta _{}0.5`$). They result in agreement with experimental data, and support the possible application of the model of line depinning in the case of long-range interactions.

A large amount of studies has been devoted to the problem of material strength and to the study of fractures in disordered media . In this paper we focus our attention on the self-affine properties of the fracture surface . By self-affinity one means that the surface coordinate $`z`$ in the direction perpendicular to the crack or fracture (x-y) plane has the following scaling properties:

$$\begin{array}{ccc}\hfill z(\lambda x,y)& \hfill =& \hfill \lambda ^{\zeta _x}z(x,y)\\ \hfill z(x,\lambda y)& \hfill =& \hfill \lambda ^{\zeta _y}z(x,y)\end{array}$$
(1)
where the two $`\zeta `$ exponents are known as roughness exponents, and the $`\widehat{y}`$ direction is chosen along the direction of propagation of the crack. Even if these two directions seem to play a different role in the morphological description of the fracture surface, experimental measurements showed that these two directions have similar scaling properties for very different materials . A unique roughness exponent is then generally considered (therefore $`\zeta _x=\zeta _y\zeta `$), and it has been claimed that it has a universal value of $`\zeta =0.8`$ . This behaviour has been confirmed for a large variety of experimental situations . However, more extended studies have also shown that in some experimental conditions, for fracture surfaces of metallic materials, one deals with a different value of $`\zeta =0.5`$ . These two values of the roughness exponent are connected to the length scale at which the crack is examined. In particular, at small length scales one observes a roughness exponent $`\zeta =0.5`$, whereas at large length scales the larger value $`\zeta =0.8`$ is found. These results have recently been connected with the velocity of the crack front , and interpreted as a quasi-static regime and a dynamic one: the value of $`\zeta =0.5`$ should be expected in the quasi-static regime. This connection comes from the suggestion that the crack surface can be thought as the trace left by the crack front  so that the problem of the surface roughness can be mapped on the problem of the evolution of a line moving in a random medium. Ertaş and Kardar  introduce a couple of non-linear Langevin equations to describe this evolution and their model can be usefully considered to evaluate the statistical properties of the evolution of a crack surface line. In this case, these equations describe longitudinal and transverse fluctuations with respect to the fracture plane containing the line velocity and the pulling force . Values obtained are close to the roughness measured at large length scales. At small length scales however, the crack front behaves as a moving line undergoing a depinning transition, and in this case results from Ertaş and Kardar model , should be revised  considering long-range interactions, to give  $`\zeta _{}`$ and $`\zeta _{}`$ equal to $`0.5`$ respectively for longitudinal and transverse fluctuations. Another interesting model has been proposed by Roux et al.  where the fracture surface is expected to be a minimal energy surface. In two dimension this problem maps directly in the random directed polymer problem: the polymer with the minimum energy is the collection of the “weakest” monomers in the medium that form a directed path. This corresponds to the surface crack of a fuse model and for brittle fractures gives a roughness exponent $`\zeta =2/3`$ in 2 dimensions. Similar arguments hold for $`d=3`$ and are discussed for a scalar model of fracture in . In this case on has $`\zeta =0.42\pm 0.01`$, differently from experimental results .

In this paper we present a numerical study of a fracture propagation for a model of quasi-static fractures to show that in such a regime, the roughness exponent is in agreement with the one found at small length scales. Comparisons with theoretical models are discussed at the end of this work. To model the fracture propagation we use a mesoscopic model known as Born Model (BM), describing the sample through a discrete collection of sites connected by springs . The statistical properties of the two dimensional BM have been previously considered in . In particular for $`d=2`$  a value of roughness of the fracture surface of $`\zeta 0.64(3)`$ in agreement with other measurements  has been found. In the BM the elastic energy of the sample under load is given by the energy of deformation for the springs connecting the sites of the sample. The elastic potential energy consists of two different terms, describing, respectively, a central force and a non-central force contribution:

$$V=\frac{1}{2}\underset{i,j}{}V_{i,j}=\frac{1}{2}\underset{i,j}{}\left\{(\alpha \beta )[(𝐮_i𝐮_j)\widehat{𝐫}_{i,j}]^2+\beta [𝐮_i𝐮_j]^2\right\}$$
(2)
where $`𝐮_i`$ is the displacement vector for site $`i`$, $`\widehat{𝐫}_{i,j}`$ is the unit vector between $`i`$ and $`j`$, $`\alpha `$ and $`\beta `$ are force constants tuning the effects of central and non-central force contributions, and the sum is over the nearest neighbour sites connected by a non-broken spring. By imposing the condition $`_𝐮V=0`$ one obtains a series of equations for the fields $`𝐮_i`$. Solving them one obtains the equilibrium positions of the springs in the sample. It has to be noticed that since one is interested in the equilibrium position one has to consider only the ratio of the two parameters $`\beta /\alpha `$ (hereafter we will consider $`\alpha =1`$ and a varying $`\beta `$). Simulations show that varying $`\beta `$ between $`1`$ and $`0`$ corresponds to varying the Poisson coefficient between $`0`$ and $`1/2`$, as expected from the theory of elasticity. At this point with a probability proportional to $`\sqrt{V_{i,j}}`$, which represents a generalized elongation, one selects a spring to remove on the fracture boundary, with the result of obtaining a connected fracture ; by doing that the boundary conditions of the system change and one has to compute a new equilibrium position. Breaking a new bond after the complete relaxation of the lattice, results in a slow velocity for the fracturing process, which mimics a quasi-static process.

The elastic springs can be arranged in different kinds of networks: however in two dimensions one has to consider a triangular lattice, for a square lattice does not correctly describe the response of the system to the applied stress. In three dimensions as well, one has to consider a network with the correct response, which results in a more complex arrangement with respect to the case of the simple cube. In our case we chose a sample described by a Face Centered Cubic (FCC) structure and we applied a mode I loading on two opposite faces by fixing their displacement field. We then applied periodic boundary conditions on a second couple of opposite faces and on one of the two remaining faces we put a starting notch (see Fig.1).

Starting from this setup we realized different simulations by stopping the algorithm when the sample is divided in two parts. At this point we started considering the surface of the fracture and we analyzed its statistical properties. We performed simulations for different values of the $`\beta `$ parameter, for 20 different samples of $`32\times 32\times 16`$ cells (each cell contains 4 sites for a total number of $`2^{16}`$ sites) for each value of $`\beta `$, for which we obtained all the relevant results. Further simulations on $`40\times 40\times 20`$ and $`50\times 50\times 25`$ FCC cells lattices were performed to verify the generality of the results. Simulations lasted from a minimum of 18-hours of CPU time for each $`32\times 32\times 16`$ cells lattice, up to more than 180 hours for a $`50\times 50\times 25`$ cells lattice, on a Digital alpha-station (500 MHz).

An example of a final fracture surface is shown in Fig.2, whereas a typical broken sample is showed in Fig.3: different colors show damaged (with at least one broken bond) and undamaged sites, and the structure of the FCC lattice.

To compute the roughness exponent, we considered different cuts of the fracture surface, some of them along the $`\widehat{x}`$ direction and some of them along the $`\widehat{y}`$ direction which is the direction of propagation of the crack. In principle we did not consider the fracture to be isotropic, but we tried to recover the roughness exponents $`\zeta _x`$ and $`\zeta _y`$. To measure the roughness of the surface, we followed the procedure described in , by introducing the two spatial correlation functions

$$\begin{array}{ccc}\hfill C_x(\rho )& \hfill =& \hfill [z(x_i+\rho ,y_i)z(x_i,y_i)]^2\\ \hfill C_y(\rho )& \hfill =& \hfill [z(x_i,y_i+\rho )z(x_i,y_i)]^2\end{array}$$
(3)
where the average $`\mathrm{}`$ is taken over the different $`x`$ and $`y`$ of the sites on the surface, and over different realization of the surfaces. Then we considered the power spectra $`\stackrel{~}{C}_{x,y}(k)`$ of the profile, that is to say we studied the Fourier transform of the previous introduced correlation functions. In this way the boundary effects are considered only in the large $`k`$ modes . For self-affine profiles these power spectra are expected to scale as

$$\stackrel{~}{C}_{x,y}(k)k^{(1+2\zeta _{x,y})}$$
(4)
Fits for the Fourier transform are shown in Fig.4. Results show that the value of $`\zeta _x`$ is equal within the error bars to the value of $`\zeta _y`$. The two directions on the fracture surface show the same statistical properties as expected for a large variety of materials; the surface can then be described by a unique roughness exponent $`\zeta `$. Moreover, the value of $`\zeta `$ does not seem to depend on the value of $`\beta `$, and is in complete agreement with the value expected for fractures in a quasi-static regime . All the results are summarized in Tab.I.

To test such a measure (as suggested by Ref.) we also studied the scaling behaviour of the surface width in direct space, along cuts perpendicular to the direction of the crack propagation: the same behaviour for $`\zeta `$ is recovered (see Fig.5). In this case results are quite striking, since a small deviation from the value of 0.5 leads to slightly displaced curves. This observation allows us to be enough confident in these results even if they extend for about one decade.

An interesting analisys is the computation of the roughness of fluctuations of the crack front, compared with experiments and with the theory of line depinning by Ertaş and Kardar. It has to be noticed that the definition of the crack front has some sort of ambiguity, because during the fracturing process more than one fracturing plane can develop, but just one among these will belong to the final fracture surface. Experimentally, this exponent can be found arresting the fracture during propagation, and injecting indian ink into the cracks under moderate vacuum. Samples are then dried and the process of fracturing is continued until complete separation is reached . From this point of view, one has to look at the border of the fracture belonging to what will be the final surface. Also, this corresponds to considering the fracture as the trace left from the crack front: this therefore belongs to the final surface.

Following this idea, we measured the roughness for the fluctuations of the crack front along the direction parallel to the line velocity ($`\zeta _{}`$) and for those perpendicular to it ($`\zeta _{}`$), during the crack evolution. In the same way as the experimental one, we mark the sites reached through the fracturing process until a certain timestep. The process then continues up to complete separation. We then look at the marked part of the final surface, to recover the crack front. Final results are from Fourier transform of an average at subsequent steps of the autocorrelation function in the steady state. In either cases we again obtained for both a value close to $`0.5`$ (see Tab.II), in good agreement with experiments.

As regards the roughness exponent we can conclude from our simulations that its value is the one characteristic of the quasi-static evolution of cracks. This result is confirmed by quantitative analysis  and is to be compared with results from the different approaches. Our conclusions seem to be different from the conclusions of the directed polymers approach as presented with their scalar problem in . The value of $`\zeta =0.42\pm 0.01`$ found with this model could be related to the different physics of the fuse networks. The use of this model in fact, comes from the assumption that in two dimensions a fracture can be described by means of a scalar model. This is not stated in three dimensions, where a description like the one that comes from the theory of elasticity can be obtained only through a vectorial model. This could also explain the difference between our results, which agree with experimental values, and those for the fuse network in .

The analysis of the roughness of longitudinal and transverse fluctuations of the crack front, give results still in agreement with experiments and seems to confirm the mapping of the crack front to a line undergoing a pinning-depinning transition. Our result is also to be compared with the one from . In this paper in fact it is stated that an explanation of the roughness in terms of a quasi-static fracturing process seems unlikely. This seems to suggest a different conclusion, as in our simulations elastic waves are cut-off through relaxation of the lattice after each bond-breaking. However, fluctuations are still enclosed in the stocastic process for the fracture.

In conclusion we presented numerical simulations of the fracture of a three dimensional sample. Our result supports the idea that the fracture roughness exponent is related to the different length scale at which the sample is analysed and then to the different dynamics of the crack. In particular for short length scales where the fracture can be identified as quasi-static the roughness exponent is $`\zeta =0.5`$. We also show that for elastic fractures one can expect isotropic behaviour in the developing of the surface: our results show no dependance from the direction on crack surface. As regards the crack front, our results agree with experimental measurements, and support the mapping to line depinning in the case of long-range interactions.

We whish to thank A.Petri and R.C.Ball for useful discussions. We also acknowledge support of EU contract N. ERBFMRXCT980183


text:
# A new approach to texture analysis of quark mass matrices

## 1 Mass matrices in the triangular basis

Texture analysis of quark mass matrices has been a subject of interest for over two decades, and the advent of new data has rendered many popular textures into oblivion. A general and systematic approach based on the observed mass and mixing spectrum becomes both welcome and necessary. What we summarize here is one such approach using triangular matrices.

It is recently observed that, as a result of the chiral nature of the electroweak force and the observed hierarchical structure among quark masses and mixing angles, the ten physical parameters of the quark mass matrices can be most simply encoded in hierarchical matrices of upper triangular form. For example, in the basis where $`M^U`$ is diagonal, $`M^D`$ is accurately given by,

$$M^U=\left(\begin{array}{ccc}m_u& & \\ & m_c& \\ & & m_t\end{array}\right)M^D=\left(\begin{array}{ccc}m_d/V_{ud}^{}& m_sV_{us}& m_bV_{ub}\\ 0& m_sV_{cs}& m_bV_{cb}\\ 0& 0& m_bV_{tb}\end{array}\right)(1+𝒪(\lambda ^4)),$$
(1)
where $`\lambda =|V_{us}|=0.22`$. Note that in this form, the unphysical right-handed rotations have been eliminated to a good approximation.

Eq.(1) is actually one of ten triangular pairs in the minimal-parameter basis (m.p.b.), where each matrix element consists of a simple product of a quark mass and certain CKM elements, and the CP-violating weak-phase is a linear combination of the phases of certain matrix elements. These ten pairs are simply related by weak basis transformation. Given any quark mass matrices, one can read off their physical content after converting them into one of the ten triangular pairs. Alternatively, one can start from upper triangular matrices and obtain mass matrices in other form, e.g. hermitian, and analyze texture zeros in the new basis.

## 2 Texture zeros of hermitian mass matrices

Based on the weak-scale quark mass relations $`m_u:m_c:m_t\lambda ^8:\lambda ^4:1`$ and $`m_d:m_s:m_b\lambda ^4:\lambda ^2:1`$, and on the hierarchical CKM mixings $`V_{us}=\lambda `$, $`V_{cb}\lambda ^2`$, $`V_{ub}\lambda ^4`$, and $`V_{td}\lambda ^3`$, we can write the properly normalized Yukawa matrices in the general hierarchical triangular form,

$$T^U=\left(\begin{array}{ccc}a_U\lambda ^8& b_U\lambda ^6& c_U\lambda ^4\\ 0& d_U\lambda ^4& e_U\lambda ^2\\ 0& 0& 1\end{array}\right),T^D=\left(\begin{array}{ccc}a_D\lambda ^4& b_D\lambda ^3& c_D\lambda ^3\\ 0& d_D\lambda ^2& e_D\lambda ^2\\ 0& 0& 1\end{array}\right).$$
(2)
Note that the diagonal elements are essentially the quark masses. The left-handed (LH) rotations are directly related to the off-diagonal elements, whose coefficients are of order one or much smaller to avoid fine-tuning in getting the CKM mixing. This direct correspondence to masses and LH rotations is a unique feature of upper triangular matrices, and this feature makes the triangular matrices especially useful in analyzing quark mass matrices. As an example, we now turn to the analysis of texture zeros of hermitian quark mass matrices.

Starting from the triangular matrices of Eq. (2), we can easily write down their corresponding hermitian form:

$`Y^U`$ $`=`$ $`\left(\begin{array}{ccc}(a_U+c_Uc_U^{}+b_Ub_U^{}/d_U)\lambda ^8& (b_U+c_Ue_U^{})\lambda ^6& c_U\lambda ^4\\ (b_U^{}+c_U^{}e_U)\lambda ^6& (d_U+e_Ue_U^{})\lambda ^4& e_U\lambda ^2\\ c_U^{}\lambda ^4& e_U^{}\lambda ^2& 1\end{array}\right)(1+𝒪(\lambda ^4))`$ (6)
$`Y^D`$ $`=`$ $`\left(\begin{array}{ccc}(a_D+b_Db_D^{}/d_D)\lambda ^4& b_D\lambda ^3& c_D\lambda ^3\\ b_D^{}\lambda ^3& d_D\lambda ^2& e_D\lambda ^2\\ c_D^{}\lambda ^3& e_D^{}\lambda ^2& 1\end{array}\right)(1+𝒪(\lambda ^2)).`$ (10)
It is seen that diagonal zeros in the hermitian matrices imply definite relations between the diagonal and off-diagonal triangular parameters (i.e. masses and LH rotation angles). As a result of the different mass hierarchies in the up and down quark sector, there is a clear asymmetry between Eqs. (6) and (10) regarding their texture zeros. For example, the $`(2,2)`$ element can be zero for $`Y^U`$ but not for $`Y^D`$; both the $`(1,1)`$ and $`(1,2)`$ elements can vanish with $`Y^U`$ but not with $`Y^D`$.

Hermitian mass matrix textures can then be analyzed as follows. We first list all possible texture pairs ($`M^U`$, $`M^D`$) directly from Eqs. (6) and (10). Each pair is then transformed into one of the ten triangular forms in the m.p.b. so that we can read-off possible relations between quark masses and mixing. Finally, the viability of each texture pair is tested by confronting its predictions with data. In this way, we find no viable hermitian pairs with six zeros, including, e.g., the Fritzsch texture. We are thus lead to consider hermitian mass matrices with five texture zeros.

## 3 Hermitian mass matrices with 5 texture zeros

Following the procedure outlined above, we identify five candidates for viable hermitian pairs with five texture zeros, as was first obtained by Ramond, Roberts, and Ross (RRR) . We now examine these five pairs analytically using the m.p.b. triangular matrices, and confront their predictions with data. The running quark mass values are taken from , and we use a recent update on $`V_{ub}/V_{cb}`$ and $`V_{td}/V_{ts}`$: $`\left|V_{ub}/V_{cb}\right|_{\mathrm{exp}}=0.093\pm 0.014`$, and $`0.15<\left|V_{td}/V_{ts}\right|_{\mathrm{exp}}<0.24`$. The results of this analysis  can be summarized as follows.

RRR patterns 1, 2, and 4 lead to the same predictions: $`\left|V_{td}/V_{ts}\right|\sqrt{m_d/m_s}=0.224\pm 0.022`$ and $`\left|V_{ub}/V_{cb}\right|\sqrt{m_u/m_c}=0.059\pm 0.006`$. Note that the numbers are independent of the scale at which the texture is valid, and that the latter prediction is on the low side and is disfavored by data.

The 5th RRR pattern also gives rise to two relations, each with two solutions depending on the sign of $`d_U`$. For the first relation,

$$\left|V_{us}\right|\left|\sqrt{\frac{m_d}{m_s}}e^{i\delta }\sqrt{\frac{m_u}{m_c}}\sqrt{\frac{V_{cb}^2}{(m_c/m_t)\pm V_{cb}^2}}\right|,$$
(11)
where the relative phase is free as CP violation depends on additional phases in the mass matrices. Thus Eq. (11) is valid with a properly chosen phase. For the second relation,

$$\left|\frac{V_{ub}}{V_{cb}}\right|\sqrt{\frac{m_u}{m_c}\left(\frac{m_c}{m_tV_{cb}^2}\pm 1\right)}=\{\begin{array}{c}0.107\pm 0.012(d_U>0,5\mathrm{a})\\ 0.068\pm 0.011(d_U<0,5\mathrm{b})\end{array}$$
(12)
where the numbers are given at $`M_Z`$ with $`V_{cb}=0.040`$. If the texture is valid at a much higher scale like $`M_{\mathrm{GUT}}`$, the weak scale predictions for the quark mixing will decrease by only a few percent. Comparing to data, we see that these predictions are marginally acceptable.

Finally, assuming the 3rd RRR texture to be valid at the weak scale, two predictions result: $`\left|V_{ub}\right|\sqrt{m_u/m_t}=0.0036\pm 0.0004`$, and $`\left|V_{us}/V_{cs}\right|\sqrt{m_d/m_s}=0.224\pm 0.022`$. Assuming the texture valid at the GUT scale will decrease the prediction for $`V_{ub}`$ by a few percent. These predictions are in excellent agreement with data.

In summary, the 3rd RRR five-texture-zero pattern is currently most favorable.

## 4 Remarks on parallel textures

There are some recent studies of 4-zero hermitian hierarchical mass matrices with parallel textures between $`M^U`$ and $`M^D`$ (i.e. $`M^U`$ and $`M^D`$ have the same texture zeros in the same locations). In particular, the generalized Fritzsch texture with zeros at (1,1), (1,3) and (3,1) has been a focus of interest. Using triangular matrices, one can easily see that generalized Fritzsch texture does not change the predictions for $`V_{ub}/V_{cb}`$ and $`V_{td}/V_{ts}`$ from its 5-zero hermitian counterparts (i.e. the 1st, 2nd, and 4th RRR patterns), and thus not favorable with data. In fact, no parallel hermitian hierarchical textures with 4 zeros is found to be favorable by current data. In other words, the viable hermitian pairs display an asymmetry between the up and down textures. This asymmetry could serve as a useful guidline in building realistic models for quark-lepton masses.

## Acknowledgments

I would like to thank T.K. Kuo, S. Mansour, and S.-H. Chiu for very enjoyable collaboration.


text:
# Spontaneous pattern formation in driven nonlinear lattices

\[

## Abstract

We demonstrate the spontaneous formation of spatial patterns in a damped, ac-driven cubic Klein-Gordon lattice. These patterns are composed of arrays of intrinsic localized modes characteristic for nonlinear lattices. We analyze the modulation instability leading to this spontaneous pattern formation. Our calculation of the modulational instability is applicable in one and two-dimensional lattices, however in the analyses of the emerging patterns we concentrate particularly on the two-dimensional case.

\]

Complex spatial patterns are often observed in systems driven away from equilibrium . Typically, the patterns emerge when relatively simple systems are driven into unstable states that will deform dramatically in response to small perturbations. As the patterns are arising from an instability, the pattern-forming behavior is likely to be extremely sensitive to small changes in system parameters. The description of deterministic pattern forming systems is typically accomplished in the form of partial differential equations such as the Navier-Stokes equations for fluids and reaction-diffusion equations for chemical systems. These phenomena have primarily been studied exclusively in continuum systems such as hydrodynamical, optical, chemical systems and liquid crystals, although more recently pattern formation of similar type have also been reported in periodically vibrated granular media .

Complementing the development of the theoretical understanding of pattern formation in continuum systems, the localized mode forming ability of discrete lattices has also received significant recent attention. There is now a fairly complete understanding of the existence and stability of localized structures (often referred to as intrinsic localized modes (ILM’s) or discrete breathers) in a variety of nonlinear lattices, undriven as well as driven. It is fair to claim that these collective patterns are well understood while the process of their creation and interaction remains relatively unexplored.

In the present communication we study the pattern forming abilities of a damped and periodically driven nonlinear lattice. Specifically, we demonstrate how a driven nonlinear (cubic) Klein-Gordon lattice forms a variety of patterns via modulational instabilities. We analyse the modulational instabilities and show how these relate to the length-scale of the patterns that are formed. Normally, the spatial extend (characteristic length-scale) of ILM’s is directly related to the frequency of the ILM’s . However, in our case the (generally different) length-scale emerging from the instability may lead to a length-scale competion, the results of which we will explore.

Model and stability of homogeneous solution. First, we study the Klein-Gordon lattice

$`\ddot{x}_n+\gamma \dot{x}_n+\omega _0^2x_n=\mathrm{\Delta }_nx_n+\lambda x_n^3+ϵ\mathrm{cos}\omega t,`$ (1)
where $`\gamma `$ is the damping parameter, $`\omega _0`$ the natural frequency of the oscillators, $`\lambda =\pm 1`$ the nonlinearity parameter and finally $`ϵ`$ is the amplitude of the ac-drive at frequency $`\omega `$. In one dimension the nearest neighbor coupling is $`\mathrm{\Delta }_nx_n=x_{n+1}2x_n+x_{n1}`$. The amplitude $`A_0`$ (and phase $`\delta _0`$) of the spatially homogeneous solution $`x_n=A_0\mathrm{cos}(\omega t+\delta _0)`$ of Eq.(1) can, within the rotating wave approximation, be shown to satisfy

$$A_0^2\left(\gamma ^2\omega ^2+(\omega ^2\omega _0^2+\frac{3}{4}\lambda A_0^2)^2\right)=ϵ^2.$$
(2)
The amplitude $`A_0`$ of the response to a driving amplitude $`ϵ`$ is shown for the soft  ($`\lambda =1`$) potential in Fig. 1.

For $`\omega <\omega _0`$ (dashed line) three solutions are possible, while for $`\omega >\omega _0`$ (solid line) a single solution is possible. A similar picture is valid for the hard ($`\lambda =1`$) potential except that the multiple solutions then appear in the case $`\omega >\omega _0`$.

Analyzing the stability of the homogeneous solution with respect to spatial perturbations, we introduce $`x_n=y+z_n`$ into Eq.(1). Assuming periodic boundary conditions, we may expand $`z_n`$ in its Fourier components $`z_n=_k\mathrm{exp}(ikn)\xi _k(t)`$, where the mode amplitude $`\xi _k(t)`$ is then governed by

$`\ddot{\xi }_k+\gamma \dot{\xi }_k+\omega _k^2\xi _k`$ $`=`$ $`{\displaystyle \frac{3}{2}}\lambda A_0^2\left[1+\mathrm{cos}(2\omega t+2\delta _0)\right]\xi _k`$ (3)
with $`\omega _k^2=\omega _0^2+4\mathrm{sin}^2(k/2)`$ denoting the linear dispersion relation of the system.

Finally, the transformation $`\xi _k(t)=\zeta _k(\omega t+\delta _0)\mathrm{exp}(\frac{\gamma }{2}(\omega t+\delta ))\zeta _k(\tau )\mathrm{exp}(\frac{\gamma }{2}\tau )`$ reduces Eq.(3) to a standard Mathieu equation

$`{\displaystyle \frac{d^2\zeta _k}{d\tau ^2}}+a\zeta _k2q\mathrm{cos}(2\tau )\zeta _k=0,`$ (4)
where

$$a=\frac{1}{4\omega ^2}\left(4\omega _k^26\lambda A_0^2\gamma ^2\right)\text{and}q=\frac{3\lambda A_0^2}{4\omega ^2}.$$
(5)
As is well-known  the Mathieu equations exhibit parametric resonances when $`\sqrt{a}i`$, where $`i=1,2,3,\mathrm{}`$. The width of the resonance regions depends on the ratio $`q/a`$ (see, e.g. ). In the framework of Eq.(4) the extent of the primary resonance $`a1`$ can easily  be estimated to be $`(a1)^2<q^2`$. However, in the presence of the damping $`\gamma `$ the resonance condition for Eq.(3) becomes

$$q^2>\frac{\gamma ^2}{\omega ^2}+(a1)^2.$$
(6)
With $`a`$ and $`q`$ defined in Eqs.(5), given $`\lambda ,\gamma ,\omega ,\omega _0`$, and $`ϵ`$, this translates into an instability band of certain wavenumbers $`k`$.

Figure 2 shows this instability band as given by Eq.(6) for parameters corresponding to the solid curve in Fig. 1. The shaded region is the band of wavenumbers that are unstable according to Eq.(6) and the dashed line indicates the most unstable wavenumber, i.e.  $`a=1`$. The effect of the damping clearly is to pinch off the instability region at a finite driving $`ϵ>0`$. Similar, using Eq.(6), instability regions can be determined for the solutions indicated by the dashed curve in Fig.1.

We have verified the presence and location of the instability band by direct simulations of Eq.(1) starting from the initial condition $`x_n(t=0)=A_0\mathrm{cos}(\delta _0)+\eta _n,`$ where, $`\eta _n`$, represents a small ($`|\eta _n|A_0`$) spatially random perturbation. This initial condition injects energy into all wavenumbers and in the presence of an unstable region of wavenumbers the dynamics will enhance the energy content in this region and thereby identify the unstable region.

The above analysis is easily extended into the case of two spatial dimensions, the only required change being that the dispersion relation now is $`\omega _\stackrel{}{k}^2=\omega _0^2+4\mathrm{sin}^2(k_x/2)+4\mathrm{sin}^2(k_y/2)`$, where the wavevector is $`\stackrel{}{k}=(k_x,k_y)`$. The instability in this case appears on an annulus in the wavevector plane, with a radius given by $`a=1`$ (see, Eq.(5)) and a width determined by Eq.(6).

Pattern formation. Numerical simulations allow us not only to verify the predicted instability band, but also to follow the full nonlinear development and saturation initiated by the instability. In particular, in regions of parameter space we obtain the spontaneous formation of patterns of distinct spatial geometry. Although we have observed this phenomenon in one as well as in two dimensions, in the present communication we focus on the two-dimensional system, where the pattern formation is very rich.

Although the dynamics show different features according to the specific region of parameter space, it is possible to trace a typical behavior as follows: Initializing the system in the spatially homogeneous state described above with a small amount of randomness added, the instability sets in after a certain number of cycles of the ac-drive, depending on the strength of the parametric resonance, i.e. on the value of $`\sqrt{a}1,2,3,\mathrm{}`$. Thereafter the system usually evolves through a sequence of different patterns (rhombi, stripes, etc.), composed of localized regions of high amplitude oscillations, before reaching a final configuration that may or may not result in a structure of definite symmetry.

Due to the sensitive response to very small changes of the parameters, determining stability regions for the different pattern geometries is a difficult task. However, in the case of a hard potential ($`\lambda =1`$), Fig. 3 shows a limited area of ($`ϵ`$,$`\omega `$)-space in which distinct spatial patterns emerge and remain stable.

This diagram is constructed following the full dynamics of the system for thousands of cycles. As our study concentrates on patterns arising from instabilities of the homogeneous solutions, we do not discuss possible hysteretic behavior as is sometimes observed in similar systems (see, e.g. Ref.).

Figure 4 shows representative examples of the spontaneously emerging patterns corresponding to the points marked in Fig. 3.

The patterns consist of localized regions of high amplitude oscillations (ILM’s) residing on a background that oscillates at the frequency, $`\omega `$, of the ac-drive. In all the considered cases we have observed the natural result that patterns are only energetically sustained when the ILM, $`\omega _{ILM}`$, and driving, $`\omega `$, frequencies are commensurate, i.e. $`\omega _{ILM}=n\omega `$, where $`n`$ is an integer. For the patterns displayed in Fig. 4, $`n=2`$. Further, the motion of the ILM’s is out of phase, i.e. , at the points in time where the background oscillation reaches its maximal excursion the ILM’s obtain their minimal amplitude such that at these points the state is completely homogeneous.

At a fixed driving $`ϵ`$, for increasing values of the frequency $`\omega `$, as in Fig. 4, we can observe the following behavior: Due to the presence of the damping $`\gamma `$, at sufficiently small $`\omega `$ the spatially homogeneous solution is stable towards all possible spatial modulations such that the flat state is sustained. However, for values near point (a) the system becomes unstable with respect to certain spatial modulations and spatial patterns in the form of large stripes emerge (Fig. 4(a)). Increasing the frequency, these stripes become thinner and denser and begin to show an increasingly clearer modulation (Fig. 4(b)). The characteristic length scale of these patterns is set by the size of the unstable $`\stackrel{}{k}`$-vector according to the above analysis. The nonlinear character of the system results in a transition towards a more isotropic geometry (rhombic) as the driving frequency is increased further (Fig. 4(c)). As in the case of the stripes, stronger localization of the ILM’s arranged in the rhombic pattern (Fig. 4(d)) is observed for even larger driving frequencies. The angle between the sides of the rhombus unit cell varies but for the hard potential it is always close to $`\pi /2`$. For values of $`ϵ`$ larger than those displayed in Fig. 3, the final mesoscopic patterns of the system dynamics are spatially disordered much like the phenomena observed in granular media . It is important to realize that the length of the unstable wavevector determines the length-scales of the final patterns, while the symmetry of the patterns is determined by the nonlinear character of the system.

As a result of the periodic boundary conditions, the length scale of the emerging patterns must be commensurate with the system size. However, we have observed that the existence of a band (variability in length and angle) of unstable $`\stackrel{}{k}`$-vectors (see Fig. 2 and related discussion) allows continuous accommodation of this constraint except for the discontinuous changes in length scales occurring when it is energetically favorable for the system to add (or subtract) an additional stripe (or row of ILM’s).

For the soft potential ($`\lambda =1`$) the variation of the amplitude $`ϵ`$, and the frequency, $`\omega `$ of the ac-drive is particularly problematic as the dynamics in this case can lead to the development of catastrophic instabilities as one or several oscillators overcome the finite barrier in the quartic potential.In all the cases we have been able to simulate, the early stage time evolution of the system is characterized by the formation of ILM’s regularly arranged in a square pattern. This spatial configuration, which is sustained for up to hundreds of cycles, seems always to suffer from a weak instability and eventually deforms into a rhombic pattern, as shown in Fig. 5.

Contrasting with the case of the hard potential, here in the soft potential the angle between the sides of the rhombus unit cell is always close to $`2\pi /3`$ (so almost hexagonal).

We now analyze this pattern more closely. The pattern shown in Fig. 5 consists of ILM’s spontaneously organized into a regular rhombic pattern. The ILM’s are spatially localized and perform harmonic temporal oscillations (at the frequency, $`\omega `$, of the ac-drive) and are therefore objects that are well described  in the literature. A particular feature of these ILM’s is that they reside on a background that oscillates at the same frequency. To expose the dynamics of the ILM’s more closely, we show in Fig. 6 a time sequence along a one-dimensional cut of the two-dimensional system (the cut is indicated by the line in Fig. 5). It should be noted that in Fig. 6 we have removed the oscillations of the background in order to expose the ILM dynamics most clearly.

Although we have focused on the ILM dynamics in the case of a soft potential, the features are analogous in the case of the hard potential and only the symmetry of the patterns is different.

In summary we have studied the modulational instability in the damped and ac-driven cubic nonlinear Klein-Gordon lattice. The analysis applies to one as well as two spatial dimensions. We have further demonstrated how these instabilities lead to a variety of mesoscopic patterns of intrinsic local modes. In the case of a hard potential we characterized the patterns in a stability diagram and in the case of the soft potential we showed that the dynamics always result in a rhombic pattern with an angle close to $`2\pi /3`$. These rhombic patterns were never observed in the case of the hard potential. This difference in the shapes of the rhombic patterns in the hard and soft cases can be understood by exploiting the analogy between the changes in the steady states of a dissipative systems and phase transitions in systems at thermodynamic equilibrium (see. e.g. Ref. ). In terms of this analogy the appearance of spatially periodic structures in a nonequilibrium system can be connected to a perturbation of the translational symmetry of thermodynamic states in equilibrium. A study based on this philosophy was pursued in Ref. with the results that the angle defining rhombic patterns is given by the coefficient of the cubic term in the system, which is precisely our observation. From the present analysis it appears that experimental studies of the pattern forming abilities of discrete systems present an excellent opportunity to study ILM’s and their mutual interactions and mesoscopic patterning. For example, optical systems  and spin systems  appear good candidates for such studies.

Research at Los Alamos National Laboratory is performed under the auspice of the US DoE.