text:
# The Random Fuse Network as a Dipolar Magnet

## Abstract

We introduce an approximate mapping between the random fuse network (RFN) and a random field dipolar Ising model (RFDIM). The state of the network damage is associated with a metastable spin configuration. A mean-field treatment, numerical solutions, and heuristic arguments support the broad validity of the approximation and yield a generic phase diagram. At low disorder, the growth of a single unstable ‘crack’ leads to an abrupt global failure. Beyond a critical disorder, the conducting network sustains significant damage before the coalescence of cracks results in global failure.

Heterogeneities are present in real materials, in a wide range of magnitudes and scales, and generate fluctuations in local toughness. The fracture mechanics then results from an interplay between disorder and stress fluctuations. In attempting to elucidate the role of the disorder, a number of classifications and characterizations of rupture have been proposed . In particular, the behavior of the damage as the applied force approaches its critical value (global rupture)—partial versus no damage, macroscopic abrupt versus continuously increasing damage—allowed to distinguish various types of fracture. Most studies were carried out either numerically or in the framework of extreme statistics (see Refs.  and references therein), while few satisfactory analytical approaches exist. In parallel, a number of authors introduced the so-called random fuse network (RFN), less involved than fracture models because of its scalar nature but still capturing the main aspects of rupture mechanics. Thus the RFN often is seen as a toy model for the investigation of mechanical rupture. The mechanical bonds are replaced by fuses burning if the current (or equivalently the electrical field) exceeds a random threshold. As the total applied current flowing through the network is increased, more and more fuses burn, thereby decreasing the effective conductivity of the network until a global failure at which no current can go through anymore. The similarity of the RFN with a driven random field Ising model at zero temperature has been recently noted, but no systematic correspondence has been established.

In this letter, we present a mapping, valid at least at the initial stages of failure, between the RFN and a driven random field Ising model with dipolar coupling. This mapping associates the configuration of the burnt fuses at a given value of the applied electric field with a zero temperature metastable state of the random field dipolar Ising model (RFDIM). The mapping between fracture and a classical model of statistical mechanics, is supported by mean-field and heuristic arguments and by a preliminary numerical study, from which a generic phase diagram with a disorder induced transition emerges.

The RFN consists of a $`d`$-dimensional lattice of fuses, each with a conductivity $`\sigma _i(𝐱)=\sigma _\mathrm{𝟎}`$ initially, where $`𝐱=a𝐦`$, $`𝐦𝐙^d`$ and $`a`$ is the lattice spacing (we will assume the continuum limit $`a0`$), and $`i=1,\mathrm{},d`$ corresponds to the orientation of the fuse. If the electric field along $`i`$ at $`𝐱`$, $`E_i(𝐱)`$, exceeds a random threshold $`_i(𝐱)`$, the fuse burns and $`\sigma _i(𝐱)=0`$. We consider a scenario in which the electric field in the direction $`i=1`$ is slowly increased from zero and $`E_i=0`$ for $`i1`$ at all times.

Assuming locality and analyticity, we can write a general equation for the time evolution of the conductivity field as

$$\eta _t\sigma _i(𝐱,t)=f_i(\{\sigma _j(𝐱,t)\},\{_j\sigma _k(𝐱,t)\},\mathrm{})\theta (\mathrm{\Delta }_i(𝐱,t)),$$
(1)
where $`\mathrm{\Delta }_i=_i(𝐱)^2E_i(𝐱,t)^2`$ and the Heaviside function ensures that the conductivity of the fuse remains $`\sigma _0`$ as long as the electric field it sustains does not exceed its threshold. The simplest choice for the function $`f_i`$ within the constraints of positivity $`(\sigma _i(𝐱,t)0)`$ and monotonicity $`(_t\sigma _i(𝐱,t)0)`$, $`f_i=\mathrm{constant}\times \sigma _i(𝐱,t)`$, yields

$$\eta _t\sigma _i(𝐱,t)=\sigma _i(𝐱,t)\theta (\mathrm{\Delta }_i(𝐱,t))$$
(2)
after absorption of the constant in a redefinition of $`\eta `$. Equation (2) incorporates two time scales: the relaxation time $`\eta `$ and the characteristic time $`\tau E_1(dE_1/dt)^1`$. We focus on the limit $`\eta /\tau 0`$, in which the conductivity relaxes instantaneously from $`\sigma _0`$ to 0 as the local field crosses the local threshold. To capture the sequential dynamics associated with the $`\eta 0`$ retardation effect, we discretize  Eq. (2) in time as

$$\sigma _i(𝐱,t+\delta t)=\sigma _i(𝐱,t)\theta (\mathrm{\Delta }_i(𝐱,t)).$$
(3)
Formally, this corresponds to setting $`\eta =\delta t`$ and discretizing Eq. (2) in the usual fashion. As expected, Eq. (3) prescribes an irreversible dynamics such that $`\sigma _i(𝐱)=0`$ at time $`t`$ implies $`\sigma _i(𝐱)=0`$ for all later times.

For a given configuration of the conductivities, Kirchoff’s laws may be recast  into the integral equation

$$E_i(𝐱)=E_0\delta _{i1}+d^dy\underset{j}{}G_{ij}(𝐱𝐲)\delta \sigma _j(𝐲)E_j(𝐲),$$
(4)
where $`E_0`$ is the magnitude of the external field applied along the first axis, $`\delta \sigma _i(𝐱)=\sigma _i(𝐱)\sigma _0`$, $`G_{ij}(𝐱)=_i_jg(𝐱)`$ the dipolar tensor \[17a\], and the Green’s function $`g`$ is defined by $`\sigma _0^2g(𝐱)=\delta (𝐱)`$. To lowest non-trivial order, Kirchoff’s laws become

$$E_i(𝐱)=E_0\left[\delta _{i1}+d^dyG_{i1}(𝐱𝐲)\delta \sigma _1(𝐲)\right].$$
(5)
As we shall see, this simple approximation captures much of the physics of network damage. First, due to the presence of the dipolar kernel $`G_{ij}`$, Eq. (5) reproduces the expected shielding of the current ‘upstream’ and ‘downstream’ to a burnt region, as well as the lateral amplification . Second, the $`k`$th term of the expansion contains a product $`_{i=1,k}\delta \sigma _1(y_i)`$ which vanishes with high probability, especially at the early stages of the damage. In addition, this approximation satisfies the necessary property that the field at the tip of a crack increases with crack growth. For one isolated burnt fuse, the field on the next parallel link in the direction perpendicular to the applied field is enhanced by a factor $`\alpha >1`$. The approximation (5) gives $`\alpha =1G_{11}(0,1)2.14`$ on the $`2d`$ square lattice (and $`\alpha =2`$ in the continuum limit \[17b\]); the exact result on the square lattice is $`\alpha =4/\pi `$. When $`n`$ such parallel fuses are broken, the field at the tip of the crack is known to scale as $`n^{1/2(d1)}`$, while in the approximation of Eq. (5)

$$E_{\mathrm{tip}}=E_1(0,\rho =n+1)=E_0(1+\underset{1<\rho <n+1}{}1/\rho ^d)$$
(6)
(where $`\rho ^2=x_2^2+\mathrm{}+x_d^2`$ is the transverse square distance) converges to a finite value as $`n`$ increases, corresponding formally to the $`d=\mathrm{}`$ case. For uniform disorder distributed according to $`p()=\theta (^2_m^2)\theta (_m^2+w^2)/w`$ with $`w<2`$ or a binary disorder $`E_{\mathrm{tip}}`$ can be larger than the largest threshold and we expect the same results as if $`E_{\mathrm{tip}}n^{1/2(d1)}`$. For a disorder distribution allowing with a high probability very large values of thresholds (such as distributions with fat tails), the scaling of $`E_{\mathrm{tip}}`$ with the size of the crack becomes relevant and we cannot use the approximation (5). Here, we present results for the uniform disorder while similar results are obtained for binary disorder.

Squaring Eq. (5) we obtain, to lowest nontrivial order, the distance $`\mathrm{\Delta }_i`$ from the local threshold in terms of the damage $`\delta \sigma _i`$, as

$`\mathrm{\Delta }_1(𝐱,t)`$ $``$ $`_1(𝐱)^2E_0^22E_0^2{\displaystyle d^dyG_{11}(𝐱𝐲)\delta \sigma _1(𝐲,t)}`$ (7)
$`\mathrm{\Delta }_j(𝐱,t)`$ $``$ $`_j(𝐱)^2E_0^2\left[{\displaystyle d^dyG_{j1}(𝐱𝐲)\delta \sigma _1(𝐲,t)}\right]^2,`$ (8)
(with $`j1`$) which, along with the update rule of Eq. (3) completely determines the state of the network under an applied field $`E_0`$.

In terms of the ‘Ising spins’ defined by $`s(𝐱,t)=2\delta \sigma _1(𝐱,t)/\sigma _0+1`$, Eq. (3) simply stipulates that each spin $`s(𝐱,t)`$ will lie along a local field $`\mathrm{\Delta }_1(𝐱,t)`$. Eq. (7) shows that this local field $`\mathrm{\Delta }_1`$ results from a uniform applied component, a random component, and a spin-spin interaction mediated by the component $`G_{11}`$ of the dipolar kernel. Therefore, as the applied field is switched from $`0`$ to $`E_0`$, spins flip in avalanches until they reach the state $`s(𝐱)=\mathrm{sgn}[\mathrm{\Delta }_1(𝐱,t=\mathrm{})]`$ which minimizes the Hamiltonian

$`={\displaystyle \frac{1}{2}}{\displaystyle d^dxd^dys(𝐱)J(𝐱𝐲)s(𝐲)}{\displaystyle d^dxh(𝐱)s(𝐱)},`$ (9)
where $`J(𝐱)=E_0^2\sigma _0G_{11}(𝐱)`$, $`h=H+_1(𝐱)^2`$, and $`H=E_0^2[1\sigma _0\stackrel{~}{G}_{11}(\mathrm{𝟎})]=E_0^2(1+1/d)`$. As the driving field $`|H|`$ is increased, more and more spins flip downward, and the state of the RFN follows the metastable (nonequilibrium) spin state of $``$, connected by its history to the initial condition $`s(𝐱)=+1`$ for all $`𝐱`$.

While similarities between the RFN and the random field Ising model (RFIM) have been pointed out in the past, here a precise mapping is elucidated. We stress that, at odds with previous studies on the RFIM in which the spin-spin coupling is taken to be short-ranged, constant, and isotropic, the coupling $`J(𝐱)`$ is rather peculiar as its magnitude is proportional to the applied field $`H`$ and its modulation is dipolar. Even in the absence of disorder, the dipolar interaction in Ising spin models leads to a rich behavior  not yet fully understood. In particular, a striped phase was identified  with alternating up and down spins. A similar scenario occurs in the RFN at very low (narrow) disorder. All the fuses remain intact until $`E_0`$ reaches the lowest threshold and the corresponding fuse burns. The dipolar coupling then relieves its longitudinal neighbors while further stressing its transverse neighbors, who burn in turn, and a stripe of burnt fuses develops, eventually spanning the whole system (the network then stops conducting and the formation of additional stripes is consequently prohibited). At high (wide) disorder, a resilient fuse in the passage of a burnt stripe will stop the growth of the crack, and new cracks will nucleate elsewhere. This suggest a transition, upon increase of the disorder, from an abrupt regime characterized by the growth of a single macroscopic unstable crack to a continuous regime in which the system is significantly damaged before its global failure.

For a qualitative understanding of the role of disorder, we first revert to a simpler version of the problem in which the dipolar kernel is replaced by a uniform infinite-range coupling $`J(𝐱)=E_0^2J/N`$, where $`J>0`$ measures the strength of the interaction and $`N`$ is the total number of spins (equivalently, fuses). This ‘democratic’ model is mean-field like in that the evolution of the local field (Eq. (7)) acting on a spin $`s(𝐱)`$ may be described in terms of a single degree of freedom, the fraction $`n_{}`$ of down spins (burnt fuses),

$$\mathrm{\Delta }_1(𝐱)=_1(𝐱)^2E_0^22JE_0^2n_{},$$
(10)
and depends on the position only through the local random threshold. Upon increase of the applied field a number of spins flip, thus increasing $`n_{}`$ which in turn causes more flips, and so on until this avalanching terminates at $`n_{}(E_0)`$, the metastable fraction of downward spins at an applied field $`E_0`$. We follow the graphical scheme of Ref.  to obtain the resulting phase diagram. For the uniform disorder of width $`w`$, clearly $`n_{}=0`$ as long as $`E_0^2<_m^2`$, the lowest threshold. If $`w<2J_m^2`$, $`n_{}`$ then jumps to 1 abruptly; if $`w>2J_m^2`$, on the other hand, $`n_{}=(E_0^2_m^2)/(w2JE_0^2)`$ increases continuously with $`E_0`$ up to $`n_{}=1`$ at $`E_0^2=(_m^2+w)/(1+2J)`$. As expected, a critical width $`w_02J_m^2`$ of the disorder separates two regimes: at low disorder, a huge avalanche is leading to global failures, while for large disorder, tiny avalanches can be stable, until cracks coalesce. We note the broad validity of this picture, at least in the mean field. In contrast with the ‘democratic fiber bundle model’  for which it was shown  that abrupt failures are but an artifact of a large discontinuity in the threshold distribution, here the phase diagram extends to distributions that are not continuous, e.g., to any uniform distribution on a support $`[_m^2,_m^2+w]`$. These mean-field results are also recovered from the homogeneous saddle point of the partition function

$$Z=\underset{\{s(𝐱)=\pm 1\}}{}e^\beta $$
(11)
This saddle point corresponds to the stability condition for a microcrack: if the electric field at the tip of the crack is larger than the largest threshold, the crack is unstable. This implies the existence of the critical width $`w_0`$ found above. It is quite remarkable that this result—which was already found by probabilistic methods—appears here naturally as the uniform saddle point or mean field.

Because of the form of the coupling in this simplified democratic model, we can carry out a similar investigation without recourse to the linear approximation of Eq. (5), but using Eq. (4) directly. Assuming $`J<1/2`$ (without which the problem is ill-defined), we find not only a similar phase diagram but also, surprisingly, the same value for the critical disorder $`w_0=2J_m^2`$. At least in this mean field version of the original problem, the linear approximation is legitimate. This also suggests that when Eq. (5) is valid, the early stages of the damage dominate the whole breakdown process.

Evidently, a mean field phase diagram may be significantly modified by fluctuations governed by a dipolar coupling. On the one hand the dipolar kernel decays as $`1/r^d`$, on the other hand its angular dependence invests the model with a so-called nonmonotonicity: burning a fuse relieves the current on a fraction of neighboring fuses. Monotonicity is a simplifying feature of many driven systems, and its lack clearly introduces complications and a possibly richer phenomenology. Furthermore, while a mean field treatment predicts the evolution of damage, it fails to capture the random fluctuations of the breakdown field $`E_b`$ as well as finite-size effects. These two features are however present in the approximation (5). In particular, it is easily shown that the maximum of the electric field scales as $`\mathrm{ln}L`$, as a consequence of the long range of the dipolar tensor. Investigation of finite size effects in the framework of this approximation present a complementary route to most of the approaches to finite size effects, which rely on extreme statistics and neglect correlations.

For a better grasp of the full problem we obtained numerical solutions of the RFN (using Eqs. (3,4)) and the RFDIM (using the approximation Eq. (7)). We found (Fig. 1) that the RFDIM reproduces quite well the behavior of the RFN. In both cases, we observe a critical value $`w_0`$ of the disorder beyond which there is a large window of damage that preceeds the global failure.

Before concluding, we propose a first attempt to describe the stability of cracks in the spirit of the Imry-Ma argument. In three dimensions, the creation of a spherical crack of radius $`R`$ modifies the total energy (cf. Eq. (9)) by

$$\mathrm{\Delta }2\overline{h}R^3wR^{3/2}+gR^3$$
(12)
where the last term comes from the dipolar interaction ($`g>0`$) and merely renormalizes the average of the random field. If the applied field $`E_0`$ is small, the average $`\overline{h}`$ is positive and the formation of a large crack is prohibited. In this case, the typical size of a crack will be $`R_0[w/(2\overline{h}+g)]^{2/3}`$. On the other hand, if $`2\overline{h}+g<0`$ it is favorable for a crack to grow indefinitely. A more precise argument based on an oblate (‘penny-shaped’) crack perpendicular to the applied field leaves the above conclusions unchanged. As $`w`$ is increased, the typical size $`R_0`$ of damaged regions becomes larger, up to a critical value $`w_0`$ at which the latter percolate. This occurs when $`R_0`$ becomes comparable to the typical distance between nucleation sites, $`a[w/(H_m^2)]^{1/3}`$ for an uniform disorder. As long as $`w<w_0`$, breakdown results from a single crack that spans the system when $`2\overline{h}+g`$ becomes negative upon increase of the applied field $`H`$. When $`w>w_0`$, localized cracks grow and percolate before $`2\overline{h}+g`$ changes sign, which allows significant damage while the network is still conducting. This heuristic argument also substantifies the broad validity of the approximation of Eq. (5) a priori only valid at the initial stages of the damage: the breakdown is controlled either by an instability or by the coalescence of (small) cracks.

In summary, we have presented an approximate mapping of the RFN to the RFDIM, and argued for its broad validity. The state of the RFN—in particular the extent of damage (number of burn fuses)—is mapped into a metastable spin state. From a mean-field investigation, numerical solutions, and heuristic arguments, a generic picture emerges, characterized by a critical value of the disorder. At low disorder, a macrosocopic crack grows until it spans the system. At large disorder, threshold fluctuations stabilize (micro)cracks and a significant precursor damage develops before the latter coalesce into a percolating structure. Our investigation were carried out for a bounded disorder. By contrast, a slowly decaying threshold distribution may modify the character of the transition or suppress it altogether, an interesting question for future study.

###### Acknowledgements.

We thank T. Garel, M. Kardar, Y.-P. Pellegrini, D. Rittel, S. Zapperi for interesting discussions. RdS is supported by a Young Scientist Grant from the Swiss National Science Foundation, the NSF through Grant No. DMR-98-05833, and the Harvard University Society of Fellows.


format: text

text:
# 1 𝐵_𝑐(𝑛_𝑠) calculated using the data from Refs. (). The dashed and dotted straight lines are fits to the data from Refs. (), respectively.

Comment on “Low-Density Spin Susceptibility and Effective Mass of Mobile Electrons in Si Inversion Layers”

In a recent Letter pudalov02 , the effective $`g`$-factor and the effective mass, $`m`$, have been studied by Pudalov et al. in a dilute 2D electron system in silicon. By analyzing Shubnikov-de Haas oscillations in superimposed parallel and perpendicular magnetic fields, the authors reproduce the strong increase of $`gm`$ with decreasing electron density previously reported in Ref. kravchenko00 . However, they contrast their data with the data obtained by our group shashkin01  and claim that the spin susceptibility (or the product $`gm`$) “increases gradually with decreasing density”, which “does not support the occurrence of spontaneous spin polarization and divergence of $`gm`$ at $`n_s=n_c`$” (here $`n_s`$ is the electron density and $`n_c`$ is the critical density for the metal-insulator transition, MIT). The purpose of this Comment is to show that all available experimental data, including those of Pudalov et al., are consistent with each other and are in favor of a spontaneous spin polarization in this 2D system and a divergence of $`gm`$ at a finite electron density.

Spin polarization and spin susceptibility were recently studied by measuring the (parallel) magnetic field $`B_c`$, required to fully polarize the electrons’ spins, using scaling of magnetoresistance shashkin01  and magnetoconductivity vitkalov01 . To compare the data of the above two groups with the newer data of Pudalov et al. pudalov02 , in Fig. 1 we plot all three sets of data (two sets in the inset for better visibility). To convert $`gm`$ from Ref. pudalov02  into $`B_c`$, we use the condition for the full spin polarization: $`g\mu _BB_c/2=\pi \mathrm{}^2n_s/2m`$, where $`\mu _B`$ is the Bohr magneton. (The factor of 2 on the right side of the equation reflects the valley degeneracy.) The agreement between all three sets of data is remarkable, especially if one takes into account that different groups used different methods, different samples, and different field/spin-polarization ranges. $`B_c`$ is a linear function of $`n_s`$ and extrapolates to zero at a finite electron density which we will designate $`n_\chi `$. The linear fit of the data from Ref. pudalov02  yields $`n_\chi =8\times 10^{10}`$ cm-2, which is identical with ours suggesting that $`n_\chi `$ is sample-independent.

Contrary to the claim made in Ref. pudalov02  of a gradual increase of $`\chi `$ with decreasing $`n_s`$, the linear dependence of $`B_cn_s/gm`$ (Fig. 1) points to the critical behavior of the spin susceptibility: $`\chi n_s/(n_sn_\chi )`$. The divergence of $`\chi `$ should occur at the sample-independent electron density $`n_\chi `$, which in the samples studied in Ref. shashkin01  coincides with the critical density $`n_c`$ for the MIT indicated by arrow. In more disordered samples, however, $`n_c`$ may be noticeably higher than $`n_\chi `$. This is the case for the samples from Ref. pudalov02  with $`n_c1\times 10^{11}`$ cm-2 being well above the expected ferromagnetic transition point $`n_\chi 8\times 10^{10}`$ cm-2. This explains why no divergence of $`gm`$ is seen by Pudalov et al. at electron densities down to $`n_c`$ in their samples. Note that even in the least disordered samples $`gm`$ is still expected to be finite near $`n_\chi `$, as it normally occurs for any ferromagnetic transition due to non-zero temperature, inhomogeneous broadening, etc.

Of course, for the spin susceptibility to diverge at $`n_s=n_\chi `$, the extrapolation of $`B_c(n_s)`$ to zero must be valid. To verify its validity, accurate data at lower densities, lower temperatures, and on much less disordered samples are needed. We emphasize that, in contrast to their claim, the method used by Pudalov et al. pudalov02  certainly cannot be applied “down to and across the 2D MIT” because at $`n_s10^{11}`$ cm-2, (i) the amplitude of oscillations is too large (and even diverges as $`T0`$) pudalov94 , which is inconsistent with the Lifshitz-Kosevich formula they use, and (ii) there are too few oscillations pudalov94  to study the beating pattern. We also note that the Lifshitz-Kosevich formula was deduced for the case of weak electron-electron interactions, and its application to strongly-correlated system is not justified.

S. V. Kravchenko
Physics Department
Northeastern University
Boston, Massachusetts 02115

A. A. Shashkin and V. T. Dolgopolov
Institute of Solid State Physics
Chernogolovka, Moscow District 142432
Russia

PACS numbers: 71.30.+h, 73.40.Qv, 73.40.Hm


text:
# Exact solutions of effective mass Schrödinger equations

## Abstract

We outline a general method of obtaining exact solutions of Schrödinger equations with a position dependent effective mass. Exact solutions of several potentials including the shape invariant potentials have also been obtained.
Keywords: Effective mass, position dependent, shape invariant

PACS numbers: 03.65.-w, 03.65.Ca, 03.65.Ge

1. Apart from being an interesting topic itself Schrödinger equations with a position dependent effective mass have found wide applications in the study of electronic properties of semiconductors , quantum dots , quantum liquids  etc. Although exact solutions of are difficult to obtain, for certain potentials the effective mass Schrödinger equation can be solved . Recently supersymmetric techniques have also been applied to obtain a few exactly solvable potentials . In this article we shall use a simple method ( which depends on the equivalence of second order differential operators ) to obtain exact solutions of Schrödinger equation with a position dependent mass. In particular we shall obtain exact solutions of a number of hitherto unknown potentials. To demonstrate the simplicity of the method we shall also rederive the solutions of some of the potentials studied earlier .

2. When the mass depends on the position the kinetic energy can be defined in several ways. Here we shall be following Lévy-Leblond  and in this case the Schrödinger equation is given by

$$\frac{d}{dx}(\frac{1}{2m(x)}\frac{d\psi (x)}{dx})+V(x)\psi (x)=E\psi (x)$$
(1)
The wave function $`\psi (x)`$ should be continuous at the mass discontinuity and the derivative of the wave function should satisfy the following condition:

$$\frac{1}{m(x)}\frac{d\psi (x)}{dx}|_{}=\frac{1}{m(x)}\frac{d\psi (x)}{dx}|_+$$
(2)
We now perform the transformation

$$\psi (x)=[2m(x)]^{\frac{1}{4}}\varphi (x)$$
(3)
on equation (1) and obtain

$$\frac{1}{2m(x)}\varphi ^{\prime \prime }(x)+\frac{1}{4}(\frac{m^{}(x)}{m^2(x)})\varphi ^{}(x)+[\frac{4m^{\prime \prime }(x)7m_{}^{}{}_{}{}^{2}(x)}{32m^3(x)}]\varphi (x)+V(x)\varphi (x)=E\varphi (x)$$
(4)
where the prime indicates differenciation with respect to $`x`$. Next we make a change of the independent variable defined by

$$\overline{x}=^x\sqrt{2m(y)}𝑑y$$
(5)
Using (5) in equation (4) we get

$$\frac{d^2\varphi (\overline{x})}{d\overline{x}^2}+\mathrm{\Omega }(\overline{x})\varphi (\overline{x})=E\varphi (\overline{x})$$
(6)
where for the sake of convenience we have used $`\varphi (x)|_{x=\overline{x}}=\varphi (\overline{x})`$ and $`\mathrm{\Omega }(\overline{x})`$ is defined by

$$\mathrm{\Omega }(\overline{x})=V(x=\overline{x})+[\frac{4m^{\prime \prime }(x)7m_{}^{}{}_{}{}^{2}(x)}{32m^3(x)}]_{x=\overline{x}}=V(x=\overline{x})+V_1(x=\overline{x})$$
(7)
It is important to note that the change of variable in (5) may not always be invertible or at least not easily invertible. But this does not really pose a problem as far as solvability of (1) is concerned. This is because $`\overline{x}`$ as a function of $`x`$ is explicitly known from (5) and if we choose $`V(x)`$ such that

$$V(x)=V_2(\overline{x})V_1(\overline{x})$$
(8)
where $`V_2(\overline{x})`$ is a solvable potential then the spectrum of (6) will be known and this in turn will give us the spectrum of (1). The corresponding wave functions can be obtained using (3). In the next section we illustrate the method with a few examples.

3. In order to deal with specific potentials it is necessary to prescribe the mass function $`m(x)`$. First we consider the choice used in :

$$m(x)=(\frac{\alpha +x^2}{1+x^2})^2,m(0)=\alpha ^2,m(\pm \mathrm{})=1$$
(9)
Then form (5) we get

$$\overline{x}=\sqrt{2}[x+(\alpha 1)tan^1x],\mathrm{}<\overline{x}<\mathrm{}$$
(10)
and using (10) in (7) we find

$$V_1(x)=\frac{(\alpha 1)}{2(\alpha +x^2)^4}[3x^4+(2\alpha 4)x^2+\alpha ]$$
(11)
As we mentioned in the last section it is now necessary to choose $`V_2(\overline{x})`$ to be a solvable potential. Let us first consider the harmonic oscillator:

$$V_2(\overline{x})=\frac{1}{4}\overline{x}^2$$
(12)
so that

$$E_n=(n+\frac{1}{2})$$
(13)
Then the original potential $`V(x)`$ is given by

$$V(x)=\frac{1}{2}[x+(\alpha 1)tan^1x]^2+\frac{(\alpha 1)}{2(\alpha +x^2)^4}[3x^4+(42\alpha )x^2\alpha ]$$
(14)
The potential (14) has harmonic oscillator spectrum given by (13) and using (3) and (5) the wave functions are found to be

$$\psi _n(x)=N_n\sqrt{\frac{\alpha +x^2}{1+x^2}}exp(\overline{x}^2/2)H_n(2^{\frac{1}{4}}\overline{x})$$
(15)
where $`\overline{x}`$ is given by (10). We note that the potential (14) is a shape invariant potential . Next we choose $`V_2(\overline{x})`$ as the Morse potential

$$V_2(\overline{x})=\lambda ^2[1exp(x)]^2,E_n=[2\lambda (n+\frac{1}{2})(n+\frac{1}{2})^2],0n[\lambda 1]$$
(16)
Then for $`V(x)`$ we get the following potential which has the same spectrum (16) as the Morse potential:

$$V(x)=\lambda ^2\left\{1exp[(x+(\alpha 1)tan^1x)]\right\}^2+\frac{(\alpha 1)}{2(\alpha +x^2)^4}[3x^4+(42\alpha )x^2\alpha ]$$
(17)
As another example let us consider the soliton potential

$$V_2(\overline{x})=\lambda (\lambda +1)sech^2\overline{x}$$
(18)
with energy

$$E_n=(\lambda n)^2,n<\lambda $$
(19)
The initial potential is then given by

$$V(x)=\lambda (\lambda +1)sech^2[\sqrt{2}(x+(\alpha 1)tan^1x)]+\frac{(\alpha 1)}{2(\alpha +x^2)^4}[3x^4+(42\alpha )x^2\alpha ]$$
(20)
with energy given by (19).

Finally we consider an important class of potential,namely,quasi exactly solvable potential  for which a part of the spectrum can be determined analytically. A typical representative potential of this class is given by

$$V(\overline{x})=\overline{x}^6(8j+3)\overline{x}^2,j=0,\frac{1}{2},1,\mathrm{}$$
(21)
and for (21) $`(2j+1)`$ levels can be determined exactly. For example for $`j=\frac{1}{2}`$ two levels $`E_{}`$ can be found:

$$E_{}=2\sqrt{2}$$
(22)
The quasi exactly solvable potential $`V(x)`$ corresponding to (1) is given by

$$V(x)=[x+(\alpha 1)tan^1x]^67[x+(\alpha 1)tan^1x]^2+\frac{(\alpha 1)}{2(\alpha +x^2)^4}[3x^4+(42\alpha )x^2\alpha ]$$
(23)
and its eigenvalues are given by (22). We note that wave functions corresponding to the potentials (17),(20) and (23) can be obtained using (3) from those of (16),(18) and (21) respectively. It is thus clear that for every potential for which spectral properties of the constant mass Schrödinger equation are known we can obtain a corresponding potential for the effective mass Schrödinger equation with identical spectral properties.

4. The problem of isospectrality of Hamiltonians (when the mass is constant) has been thoroughly investigated . When the mass depends on the space coordinate has also been investigated using supersymmetry . Here we shall demonstrate that two Hamiltonians can indeed be isospectral even when both the mass and the potential are different. To show this we choose a mass of the form

$$m(x)=(\frac{\alpha +x^2}{1+x^2})^4$$
(24)
Then from (3) we get

$$\overline{\overline{x}}=2\sqrt{2}[2x+\frac{(\alpha 1)^2x}{(1+x^2)}+(\alpha 1)(\alpha 3)tan^1x],\mathrm{}<\overline{\overline{x}}<\mathrm{}$$
(25)
while (7) gives us

$$V_1(x)=\frac{(\alpha 1)(1+x^2)^2}{(\alpha +x^2)^6}[3x^4+(5\alpha 7)x^2+\alpha ]$$
(26)
Now for $`V_2(\overline{\overline{x}})`$ we take

$$V_2(\overline{\overline{x}})=\frac{\overline{\overline{x}}^2}{4}$$
(27)
so that its spectrum is given by (13). In this case $`V(x)`$ is given by

$$V(x)=\frac{1}{8}[2x+\frac{(\alpha 1)^2x}{(1+x^2)}+(\alpha 1)(\alpha 3)tan^1x]^2+\frac{(\alpha 1)(1+x^2)^2}{(\alpha +x^2)^6}[3x^4+(75\alpha )x^2\alpha ]$$
(28)
Since (27) is a harmonic oscillator potential it follows that the effective mass Schrödinger equations for the potentials (14) and (28) with masses given by (9) and (24) respectively are exactly isospectral.

5. In this article we have described a general method of obtaining exact solutions of effective mass Schrödinger equation. In particular the method has been applied to obtain solutions of several potentials which have so far been unknown. It has also been shown that Schrödinger equations with different effective masses as well as different potentials can be isospectral.


text:
# Is the X-ray variable source in M82 due to gravitational lensing ?

## 1 Introduction

Observations of the most famous starburst galaxy M82 with the High-Resolution Camera on board the Chandra X-Ray Observatory, show that there are nine sources in the central $`1^{}\times 1^{}`$ region, but no source was detected at the galactic center(Matsumoto et al., 2001). Comparing the observations on 1999 October 28 and on 2000 January 20, the authors found an extremely large time variability of the source CXOM82 J095550.2+694047, which is $`9^{\prime \prime }`$ away from the galactic center. They concluded that this source is the origin of the hard X-ray time variability of M82 detected with ASCA. Assuming a spectral shape obtained by the ASCA observation, its luminosity in the 0.5 - 10 KeV band changed from $`1.2\times 10^{40}`$ erg s-1 on 1999 October 28 to $`8.7\times 10^{40}`$ erg s-1 on 2000 January 20.

It is difficult to explain such a short-term variability in terms of a supernova remnant. If the spectral shape of the source is described by an absorbed thermal bremstrahlung model, the observations show that the probability of such a bright source existing in the $`1^{}\times 1^{}`$ field is $`0.3\%`$. Consequently, this source is probably not a background AGN.

CXOM82 J095550.2+694047 may be a medium-massive black hole (Matsumoto et al., 2001), and the possibilities that an X-ray binary source whose jet is strongly beamed towards us cannot be excluded. In this letter, we explore an alternative possibility that the time variability of the source is originated from the microlening effect by MACHO along the line of sight.

## 2 Microlensing explanation

The source CXOM82 J095550.2+694047 is $`9^{\prime \prime }`$ away from the center of M82 (its distance from us is 3.9 Mpc). It may be located in a star cluster with typical size of 2.34 - 10 pc (de Grijs et al., 2001)and typical velocity dispersion of member stars of 15 km/s.

If a foreground star is traversing across the line of sight, the apparent magnitude of the source would be magnified according to microlensing theory. The duration of the crossing time is $`T=2a_E/v`$, in which $`a_E`$ is the Einstein radius defined as

$$a_E=\left(\frac{4GM}{c^2}\frac{D_{ds}D_d}{D_s}\right)^{1/2},$$
(1)
where $`D_d`$, $`D_s`$ and $`D_{ds}`$ are the angular diameter distances to the lens, to the source and from the lens to the source, respectively, $`v`$ is the relative velocity perpendicular to the line of sight between the source and the lens.

The mass of microlens is related to the crossing time as (Wu, 1996)

$$M=\frac{c^2}{32G}\frac{D_s}{D_dD_{ds}}\frac{v^2T^2}{0.4919},$$
(2)
where we have chosen the total maximum magnification to be $`\mu _T=7.0`$, since the luminosity of the source in the 0.5 - 10 KeV band changed from $`1.2\times 10^{40}`$ erg s-1 on 1999 October 28 to $`8.7\times 10^{40}`$ erg s-1 on 2000 January 20.

A microlens of mass of $`M_{}`$ could lie at any point between the source and the observer, and different positions would lead to different crossing time. For the source CXOM82 J095550.2+694047, the observed time duration must exceed 84 days. For simplicity, we can take $`T=2\times 84`$ days for an approximate estimate. The microlens mass versus its distance is plotted in Fig.1 for three kinds of lensing models: the microlens is located in our Galaxy with a typical velocity of 200km/s, in M82 with a typical velocity of 100km/s and in the star cluster of the source with a velocity of 15km/s.

At three different positions of a microlens in the Galaxy, the M82 and the star cluster of the source, the same microlens mass would lead to different crossing time durations, as is shown in Fig.2. We set $`D_s=3.9`$Mpc in all the panels. Because of the symmetry between $`D_d`$ and $`D_{ds}`$ in equation (2), the first and second panels both in Fig.1 and Fig.2 look very similar.

## 3 Discussions and conclusions

It can be seen clearly from Fig.1 and Fig.2 that the microlens object, if any, should be located in M82 or the Galactic halo. However, we can exclude the possibility of the self-lensing by a star in the star cluster where CXOM82 J095550.2+694047 may be harbored if the typical size of a star cluster is smaller than $`10`$ pc.

The probability that a source is gravitationally lensed is described by the optical depth which amounts to $`10^6`$. Indeed, the event like the hard X-ray variable source CXOM82 J095550.2+694047 seen in M82 is very rare if it is due to microlensing. A conclusive resolution needs a detailed sample of the light curve of CXOM82 J095550.2+694047. Whether or not the light curve demonstrates the time-symmetry around the maximum intensity will be a crucial test for the lensing hypothesis.

I’m grateful to Prof. Xiang-Ping Wu for reading through the manucript and giving constructive suggestions. This work was supported by the National Science Foundation of China, under grant 19725311 and 19873014.


text:
# Subjet Multiplicity in Quark and Gluon Jets at DØ

## I Introduction

The Tevatron proton-antiproton collider is a rich environment for studying high energy physics. The dominant process is jet production, described in Quantum Chromodynamics (QCD) by scattering of the elementary quark and gluon constituents of the incoming hadron beams. In leading order (LO) QCD, there are two partons in the initial and final states of the elementary process. A jet is associated with the energy and momentum of each final state parton. Experimentally, however, a jet is a cluster of energy in the calorimeter. Understanding jet structure is the motivation for the present analysis. QCD predicts that gluons radiate more than quarks. Asymptotically, the ratio of objects within gluon jets to quark jets is expected to be in the ratio of their color charges $`C_A/C_F=9/4`$.

## II The $`k_T`$ Jet Algorithm

We define jets in the DØ detector  with the $`k_{}`$ algorithm . The jet algorithm starts with a list of energy preclusters, formed from calorimeter cells or from particles in a Monte Carlo event generator. The preclusters are separated by $`\mathrm{\Delta }=\sqrt{\mathrm{\Delta }\eta ^2+\mathrm{\Delta }\varphi ^2}>0.2`$, where $`\eta `$ and $`\varphi `$ are the pseudorapidity and azimuthal angle of the preclusters. The steps of the jet algorithm are:

1. For each object $`i`$ in the list, define $`d_{ii}=p_{T,i}^2`$, where $`p_T`$ is the energy transverse to the beam. For each pair $`(i,j)`$ of objects, also define $`d_{ij}=min(p_{T,i}^2,p_{T,j}^2)\frac{\mathrm{\Delta }_{ij}^2}{D^2}`$, where $`D`$ is a parameter of the jet algorithm.

2. If the minimum of all possible $`d_{ii}`$ and $`d_{ij}`$ is a $`d_{ij}`$, then replace objects $`i`$ and $`j`$ by their 4-vector sum and go to step 1. Else, the minimum is a $`d_{ii}`$ so remove object $`i`$ from the list and define it to be a jet.

3. If any objects are left in the list, go to step 1.

The algorithm produces a list of jets, each separated by $`\mathrm{\Delta }>D`$. For this analysis, $`D=0.5`$.

The subjet multiplicity is a natural observable of a $`k_{}`$ jet. Subjets are defined by rerunning the $`k_{}`$ algorithm starting with a list of preclusters in a jet. Pairs of objects with the smallest $`d_{ij}`$ are merged successively until all remaining $`d_{ij}>y_{cut}p_T^2(jet)`$. The resolved objects are called subjets, and the number of subjets within the jet is the subjet multiplicity $`M`$. The analysis in this article uses a single resolution parameter $`y_{cut}=10^3`$.

## III Jet Selection

In LO QCD, the fraction of final state jets which are gluons decreases with $`xp_T/\sqrt{s}`$, the momentum fraction of initial state partons within the proton. For fixed $`p_T`$, the gluon jet fraction decreases when $`\sqrt{s}`$ is decreased from 1800 GeV to 630 GeV. We define gluon and quark enriched jet samples with identical cuts in events at $`\sqrt{s}=1800`$ and 630 GeV to reduce experimental biases and systematic effects. Of the two highest $`p_T`$ jets in the event, we select jets with $`55<p_T<100`$ GeV and $`|\eta |<0.5`$.

## IV Quark and Gluon Subjet Multiplicity

There is a simple method to extract a measurement of quark and gluon jets on a statistical basis, using the tools described in the previous sections. $`M`$ is the subjet multiplicity in a mixed sample of quark and gluon jets. It may be written as a linear combination of subjet multiplicity in gluon and quark jets:

$$M=fM_g+(1f)M_q$$
(1)
The coefficients are the fractions of gluon and quark jets in the sample, $`f`$ and $`(1f)`$, respectively. Consider Eq. (1) for two similar samples of jets at $`\sqrt{s}=1800`$ and 630 GeV, assuming $`M_g`$ and $`M_q`$ are independent of $`\sqrt{s}`$. The solutions are

$$M_q=\frac{f^{1800}M^{630}f^{630}M^{1800}}{f^{1800}f^{630}}$$
(2)
$$M_g=\frac{\left(1f^{630}\right)M^{1800}\left(1f^{1800}\right)M^{630}}{f^{1800}f^{630}}$$
(3)
where $`M^{1800}`$ and $`M^{630}`$ are the experimental measurements in the mixed jet samples at $`\sqrt{s}=1800`$ and 630 GeV, and $`f^{1800}`$ and $`f^{630}`$ are the gluon jet fractions in the two samples. The method relies on knowledge of the two gluon jet fractions.

## V Results

The HERWIG 5.9 Monte Carlo event generator provides an estimate of the gluon jet fractions. The method is tested using the detector simulation and CTEQ4M PDF. We tag every selected jet in the detector as either quark or gluon by the identity of the nearer (in $`\eta \times \varphi `$ space) final state parton in the QCD 2-to-2 hard scatter. Fig. 1 shows that gluon jets in the detector simulation have more subjets than quark jets. The tagged subjet multiplicity distributions are similar at the two center of mass energies, verifying the assumptions in § IV.

We count tagged gluon jets and find $`f^{1800}=0.59\pm 0.02`$ and $`f^{630}=0.33\pm 0.03`$, where the uncertainties are estimated from different gluon PDF’s. The nominal gluon jet fractions and the Monte Carlo measurements at $`\sqrt{s}=1800`$ and 630 GeV are used in Eqs. (2-3). The extracted quark and gluon jet distributions in Fig. 1 agree with the tagged distributions and demonstrate closure of the method.

Figure 2 shows the raw subjet multiplicity in DØ data at $`\sqrt{s}=1800`$ GeV is higher than at $`\sqrt{s}=630`$ GeV. This is consistent with the prediction that there are more gluon jets at $`\sqrt{s}=1800`$ GeV compared to $`\sqrt{s}=630`$ GeV, and gluons radiate more than quarks. The combination of the distributions in Fig. 2 and the gluon jet fractions gives the raw subjet multiplicity distributions in quark and gluon jets, according to Eqs. (2-3).

The quark and gluon raw subjet multiplicity distributions need separate corrections for various detector-dependent effects. These are derived from Monte Carlo, which describes the raw DØ data well. Each Monte Carlo jet in the detector simulation is matched (within $`\mathrm{\Delta }<0.5`$) to a jet reconstructed from particles without the detector simulation. We tag detector jets as either quark or gluon, and study the subjet multiplicity in particle jets $`M^{ptcl}`$ vs. that in detector jets $`M^{det}`$. The correction unsmears $`M^{det}`$ to give $`M^{ptcl}`$, in bins of $`M^{det}`$. Figure 3 shows the corrected subjet multiplicity is clearly larger for gluon jets compared to quark jets.

The gluon jet fractions are the largest source of systematic error. We vary the gluon jet fractions by the uncertainties in an anti-correlated fashion at the two values of $`\sqrt{s}`$ to measure the effect on $`R`$. The systematic errors listed in Table I are added in quadrature to obtain the total uncertainty in the corrected ratio $`R=\frac{M_g1}{M_q1}=1.84\pm 0.15(\mathrm{stat})_{0.18}^{+0.22}(\mathrm{sys})`$.

## VI Conclusion

We extract the $`y_{cut}=10^3`$ subjet multiplicity in quark and gluon jets from measurements of mixed jet samples at $`\sqrt{s}=1800`$ and 630 GeV. On a statistical level, gluon jets have more subjets than quark jets. We measure the ratio of additional subjets in gluon jets to quark jets $`R=1.84\pm 0.15(\mathrm{stat})_{0.18}^{+0.22}(\mathrm{sys})`$. The ratio is well described by the HERWIG parton shower Monte Carlo, and is only slightly smaller than the naive QCD prediction 9/4.

## Acknowledgements

We thank the staffs at Fermilab and collaborating institutions for their contributions to this work, and acknowledge support from the Department of Energy and National Science Foundation (U.S.A.), Commissariat à L’Energie Atomique (France), State Committee for Science and Technology and Ministry for Atomic Energy (Russia), CAPES and CNPq (Brazil), Departments of Atomic Energy and Science and Education (India), Colciencias (Colombia), CONACyT (Mexico), Ministry of Education and KOSEF (Korea), and CONICET and UBACyT (Argentina).


text:
# Measurement of the quadratic slope parameter in the 𝐾_𝐿→3𝜋⁰ decay Dalitz plot

## Introduction

The $`K3\pi `$ decays Dalitz plot distributions can be expanded in powers of Dalitz plot variables $`u`$ and $`v`$ :

$$\left|M(u,v)\right|^21+gu+jv+hu^2+kv^2$$
(1)
$$u=\frac{\left(s_3s_0\right)}{m_{\pi ^+}^2},v=\frac{\left(s_1s_2\right)}{m_{\pi ^+}^2}$$

$$s_0=\frac{s_1+s_2+s_3}{3},s_i=\left(P_KP_i\right)^2,i=1,2,3$$

where $`P_K`$ and $`P_i`$ are the four-momenta of the decaying kaon and $`i`$th pion ($`i=3`$ for the “odd charge” pion), respectively.

In the case of the $`K_L3\pi ^0`$ decay the expression (1) reduces to :

$$\left|M_{000}(R^2,\theta )\right|^21+hR^2,$$
(2)
$$R^2=\left(u^2+v^2/3\right),\theta =\mathrm{arctan}\left(v/\sqrt{3}u\right)$$

because of the identical final state particles. A positive/negative value of the quadratic slope parameter $`h`$ would mean that asymmetric/symmetric final states are favoured.

Based only on isospin symmetry and some general assumptions it is possible to define a relation

$$\rho h+3k\frac{g^2}{4\mathrm{cos}^2\beta }$$
(3)
between the quadratic $`h,k`$ and linear $`g`$ Dalitz plot slope parameters and the final state interaction phases $`\beta `$ which should be approximately the same for all $`K3\pi `$ decays . For $`K_L3\pi ^0`$ decays the parameter $`\rho `$ is equal to $`2\times h`$ ($`kh/3`$) and is not influenced by radiative corrections since all the final state particles are neutral. Thus a high statistics measurement of the quadratic slope parameter in the $`K_L3\pi ^0`$ decay Dalitz plot could establish strong constraints to the Dalitz plot slope parameters and f.s.i. phases in other $`K3\pi `$ decays. Also, it could be useful in the evaluation of radiative corrections applied to $`K^\pm \pi ^\pm \pi ^\pm \pi ^{}`$ and $`K_L\pi ^+\pi ^{}\pi ^0`$ decays. On the other hand, by combining the $`K_L3\pi ^0`$ decay Dalitz plot slope parameter $`h`$ with the linear and quadratic slope parameters in the $`K^\pm \pi ^\pm \pi ^\pm \pi ^{}`$ decay it is possible to probe the validity of the $`\mathrm{\Delta }I=1/2`$ rule.

The $`K_L3\pi ^0`$ Dalitz plot slope parameter $`h`$ has been estimated in the framework of ChPT ,,. The analysis includes a nonzero $`\mathrm{\Delta }I=3/2`$ amplitude in the quadratic term of the Dalitz plot expansion (1) by taking into account the next-to-leading $`O(p^4)`$ corrections. The phenomenological coupling constants in the chiral Lagrangian are evaluated using a fit to the recent experimental data in both $`K2\pi `$ and $`K3\pi `$ decays.

In this paper we present a measurement based upon a data set collected with the NA48 detector in 1998. The detector acceptance has been estimated using a detailed Monte Carlo (MC) simulation.

## Detector setup and data taking

The NA48 experiment is dedicated to the measurement of direct CP violation in $`K_{L,S}^02\pi `$ decays. The $`K_L`$ beam is produced by $`450`$ GeV/c protons from the CERN SPS hitting a $`2`$ mm diameter and $`400`$ mm long beryllium target at a production angle of $`2.4`$ mrad. Passing through a set of collimators, a $`\pm 0.15`$ mrad divergent $`K_L`$ beam enters the decay volume $`126`$ m downstream of the target. The evacuated $`89`$ m long decay volume is followed by a helium tank which contains four drift chambers of the charged particle spectrometer. After the tank a scintillator hodoscope, a liquid krypton electromagnetic calorimeter, a hadron calorimeter and muon veto counters are placed.

The liquid krypton calorimeter (LKr) is designed to measure the energy, position and timing of electromagnetic showers . The $`127`$ cm long detector consists of 13212 readout cells in a projective tower geometry which points to the middle of the decay volume. Each cell of $`2\times 2`$ cm2 cross section is made of copper-beryllium ribbons which are extended longitudinally in a $`\pm 48`$ mrad accordion structure. The readout cells are contained in a cryostat filled with about $`20`$ t liquid krypton at $`121`$ K. The initial current induced on the copper-beryllium electrodes is measured by $`80`$ ns FWHM and digitized by $`40`$ MHz FADCs. For the 1998 data taking period, the energy and position resolutions of the calorimeter were determined to be:

$$\frac{\sigma (E)}{E}=\frac{(3.2\pm 0.2)\%}{\sqrt{E}}\frac{(0.09\pm 0.01)}{E}(0.42\pm 0.05)\%$$
(4)
$$\sigma (x)\sigma (y)\frac{0.4\mathrm{cm}}{\sqrt{E}}0.05\mathrm{cm}$$
(5)
with $`E`$ measured in GeV.

The time resolution was better than $`300`$ ps for photons with energies above $`20`$ GeV.

The energy non-linearity in the calorimeter response was found to be less than $`1\times 10^3`$ between $`6`$ GeV and $`100`$ GeV .

The LKr detector contains the so-called neutral hodoscope which is made of a 4 mm thick plane of scintillating fibres placed near the maximum of the e.m. shower development. The neutral hodoscope is used in the estimation of the main trigger inefficiency and for an additional measurement of the event time.

A description of the whole experimental setup can be found elsewhere .

The $`K_L3\pi ^06\gamma `$ data sample was acquired by three different triggers. The first one is the neutral $`2\pi ^0`$ trigger (NUT $`2\pi ^0`$) in which each $`25`$ ns the calibrated analogue sums of $`2\times 8`$ LKr cells are used to construct $`64`$-channel $`x`$ and $`y`$ projections of the energy deposited in the calorimeter . Based on these projections the total deposited energy $`E_{sum}`$, the first and second moments of the $`x`$ and $`y`$ energy distributions, and the time of the energy peaks are computed. The centre-of-gravity COG of the event and the longitudinal vertex position are reconstructed as:

$$COG^{trig}=\frac{\sqrt{M_{1,x}^2+M_{1,y}^2}}{E_{sum}}$$
(6)
$$z_{vertex}^{trig}=z_{LKr}\frac{\sqrt{E_{sum}\left(M_{2,x}+M_{2,y}\right)\left(M_{1,x}^2+M_{1,y}^2\right)}}{m_K}$$
(7)
where $`M_{1,x}`$ , $`M_{1,y}`$ and $`M_{2,x}`$, $`M_{2,y}`$ are the first and second moments of the energy deposition, $`z_{LKr}`$ is the longitudinal position of the calorimeter and $`m_K`$ is the nominal kaon mass.

The proper decay time $`\tau ^{trig}`$ from the beginning of the decay region (just after the final collimator) in units of the $`K_S`$ lifetime $`\tau _S`$ is derived taking into account calibration constants. The events are accepted by the NUT $`2\pi ^0`$ trigger if $`E_{sum}>50`$ GeV, $`COG^{trig}<15`$ cm and $`\tau ^{trig}/\tau _S<5`$. In addition, less than $`6`$ energy peaks within $`13.5`$ ns in both the $`x`$ and $`y`$ projections are required.

The second trigger (NUT $`3\pi ^0`$) uses the same hardware chain as NUT $`2\pi ^0`$ trigger but is specially set-up to select $`3\pi ^0`$ events by applying no condition on the number of the energy peaks.

The third trigger (Nhodo) is based on the information from the neutral hodoscope. It requires a coincidence of signals from the upper and lower or the left and right parts of the calorimeter. In order to reduce the trigger output rates the Nhodo trigger and the NUT $`3\pi ^0`$ one were properly down-scaled.

## Data analysis

The $`K_L3\pi ^06\gamma `$ events were selected from all events which met at least one of the trigger requirements. The selection criteria required 6 or more LKr clusters satisfying the following conditions:

* The energy of the cluster had to be between $`3`$ GeV and $`100`$ GeV;
* To avoid energy losses in the LKr, the distance from the cluster to the closest dead calorimeter cell was required to be greater than $`2`$ cm and the cluster had to be more than $`5`$ cm from the beam pipe and the outer edge of the calorimeter.

All possible combinations of $`6`$ clusters which passed these requirements were considered. In addition, the following further selection criteria were applied on each combination:

* The distance between all cluster pairs had to be greater than $`10`$ cm in order to avoid overlapping of the clusters;
* All $`6`$ clusters had to lie within $`5`$ ns from the average cluster time;
* The sum of the clusters energies had to exceed $`60`$ GeV which is sufficiently far from the NUT $`2\pi ^0`$ and NUT $`3\pi ^0`$ trigger threshold of $`50`$ GeV;
* The radial position of the center-of-gravity

$$COG=\frac{\sqrt{\left(_{i=1}^6E_ix_i\right)^2+\left(_{i=1}^6E_iy_i\right)^2}}{_{i=1}^6E_i},$$

where $`E_i`$, $`x_i`$, $`y_i`$ are the energies and positions of the six selected clusters, had to be less than $`10`$ cm;
* No additional clusters with an energy $`>1.5`$ GeV were allowed within $`\pm 3`$ ns from the event time to minimize possible accidental effects.

The longitudinal vertex positions were reconstructed analogous to (7):

$$z_{vertex}=z_{LKr}\frac{\sqrt{_{i=1}^6_{j>i}^6E_iE_j\left[\left(x_ix_j\right)^2+\left(y_iy_j\right)^2\right]}}{m_K}$$
(8)
where $`E_i`$ and $`x_i`$, $`y_i`$ are the energies and positions of the six selected clusters.

For each combination the invariant masses $`m_1`$, $`m_2`$ and $`m_3`$ of all $`15`$ possible photon pairings were computed using $`z_{vertex}`$. By applying a $`\chi ^2`$ criteria, the combination and pairing most compatible with the hypothesis that $`m_1`$, $`m_2`$ and $`m_3`$ are equal to the nominal $`\pi ^0`$ mass were picked.

To ensure the purity of the sample, events which had one or more reconstructed tracks in the spectrometer or more than one hit per plane in the third and fourth drift chambers were rejected.

To improve the resolution on the Dalitz plot variables $`R^2`$ and $`\mathrm{tan}(\theta )`$ and to reduce differences in the energy resolutions and non-linearities between data and MC which may bias the measurement of the quadratic slope parameter $`h`$, all the events were passed through a kinematical fitting procedure with constraints. Assuming the neutral kaon coming from the $`K_L`$ target and a $`K3\pi ^06\gamma `$ decay the cluster energies and positions are adjusted to minimize

$$\chi _{fit}^2=\underset{i=1}{\overset{6}{}}\frac{\left(E_iE_i^{fit}\right)^2}{\sigma \left(E_i\right)^2}+\underset{i=1}{\overset{6}{}}\frac{\left(x_ix_i^{fit}\right)^2}{\sigma \left(x_i\right)^2}+\underset{i=1}{\overset{6}{}}\frac{\left(y_iy_i^{fit}\right)^2}{\sigma \left(y_i\right)^2}$$
(9)
where $`\sigma \left(E_i\right)`$, $`\sigma \left(x_i\right)`$, $`\sigma \left(y_i\right)`$ are given by (4),(5) and $`E_i^{fit}`$, $`x_i^{fit}`$, $`y_i^{fit}`$ are the adjusted energies and positions of the photons as functions of the $`15`$ decay parameters (the kaon energy, the vertex position, the three Euler angles of the decay plane in the kaon rest frame, the Dalitz plot variables $`R^2`$ and $`\theta `$ and the 2 angles of the $`\gamma `$ directions for each of the decaying $`\pi ^0`$s). While the effect of the kinematical fitting on the $`E_K`$ and $`z_{vertex}`$ resolutions was small, the $`R^2`$ and $`\theta `$ resolutions were improved significantly. In particular the RMS widths of $`(R^2R_{true}^2)`$ and $`(\theta \theta _{true})`$ from $`8.8\times 10^2`$ and $`48`$ mrad became $`3.0\times 10^2`$ and $`31`$ mrad, respectively.

Finally a cut $`\chi _{fit}^2<8.5`$ was applied to exclude the tails in the $`\chi _{fit}^2`$ distribution which differ for data and MC.

The numbers of data and MC events which passed the selection criteria are presented in Table 1.

To estimate the quadratic slope parameter $`h`$, the $`R^2`$ distribution for the data events was corrected for the detector acceptance. For this purpose about $`8\times 10^8`$ Monte Carlo events were generated with $`h=0`$. Then the quadratic slope parameter $`h`$ was evaluated by a linear fit to the ratio of the data and MC $`R^2`$ distributions (Fig. 1).

The fit was done separately for each of the triggers using the corresponding MC samples within appropriate decay regions.

An important feature of the analysis was that the ratio of the data and MC $`R^2`$ distributions was fitted up to $`1.9`$ (Fig. 1). This allowed us to significantly improve the stability of the result by excluding the part of the Dalitz plot most affected by the energy resolution and non-linearities. Since the shape of the detector acceptance on $`R^2`$ depends mainly on the kaon energy, the whole analysis was done in kaon energy bins. Combining the measured values of the Dalitz plot slope parameter $`h`$ presented in Figure 2

and Table 1 we get a overall result of $`(6.1\pm 0.9)\times 10^3`$.

## Systematics

Contributions of the detector acceptance estimation, the $`\chi _{fit}^2`$ cut, the trigger conditions and the differences in the energy resolution, non-linearities and energy scale between the data and MC to the systematic error were considered.

In order to study possible uncertainties in the acceptance we divided the data and MC samples upon the decay length $`\tau /\tau _s`$, the longitudinal vertex position $`z_{vertex}`$ and the centre-of-gravity COG. For each considered trigger sample the maximum deviation between the average value of $`h`$ from Table 1 and the values obtained averaging over the dependence of a second variable ($`\tau /\tau _s`$, $`z_{vertex}`$, COG) was taken as systematic error (Table 2).

An important check of the reliability of the MC was done by comparing the data and MC Dalitz plot variable $`\theta `$ (Fig. 3).

A linear fit to the ratio of the data and MC $`\theta `$ distributions clearly points to a very good agreement between experimental data and the simulation.

The second considered source of systematics was the difference between the data and MC cluster energy resolution. Using the measured and kinematically adjusted cluster energies we established a strong control over the resolutions seen in the data and the MC. The differences in the resolutions were adjusted by an additional smearing to the cluster energies in the MC. The maximum change in $`h`$ after these adjustments was found to be $`\pm 0.20\times 10^3`$ and was taken as a systematic error. A similar approach was applied in the study of the systematics due to the energy non-linearities leading to a $`\pm 0.17\times 10^3`$ systematic error on the $`h`$. To take into account small non-gaussian tails in the cluster energy resolution seen in the data, a simplified parameterization was used to add the tails to the MC events. The effect of the tails was found to be negligible. Changing the energy scale in the MC by $`\pm 1\times 10^3`$ we got a shift in $`h`$ of $`\pm 0.17\times 10^3`$.

The influence of the trigger conditions was studied by subdividing the data and MC in $`\tau /\tau _s`$ bins. Although that for the NUT $`2\pi ^0`$ and Nhodo trigger samples no systematic effect was found, the result for NUT $`3\pi ^0`$ trigger sample was affected by a small NUT inefficient region near the $`\tau ^{trig}/\tau _s`$ cut. The data and MC were processed excluding this region and the change in the result of $`0.1\times 10^3`$ gave us a systematic error due to the trigger.

The stability of the result upon the $`\chi _{fit}^2`$ cut was checked by loosening it to $`17`$. No statistically significant effect on $`h`$ was established. Furthermore, the result was found to be stable upon different cuts on the clusters energies and positions in LKr.

No influence of a possible background contamination from $`\gamma `$ conversions and Dalitz decays was observed.

All the sources of systematic error were added quadratically and yield a total systematic error of $`0.5\times 10^3`$.

In addition, a second independent analysis has been performed yielding a consistent result.

## Conclusions

Our result for the Dalitz plot slope parameter in the $`K_L3\pi ^0`$ decay is $`h=(6.1\pm 0.9_{stat}\pm 0.5_{syst})\times 10^3`$, where the errors are statistical and systematic, respectively. This result is in an agreement with the value of $`(12\pm 4)\times 10^3`$ predicted in . The previous measurement of $`\left(3.3\pm 1.1_{stat}\pm 0.7_{syst}\right)\times 10^3`$ performed by the E731 experiment , is comparable with the result presented in this paper.

## Acknowledgements

We would like to thank the technical staff of the participating laboratories, universities and affiliated computing centres for their support and co-operation.