text:
# 1 Introduction

## 1 Introduction

Algebraic models are popular in several branches of physics . In ideal cases, algebraic models are able to describe in good approximation several properties of a physical system using a limited number of appropriate (usually collective) degrees of freedom. In addition to providing successful predictions by themselves, exactly soluble algebraic models can also serve as a useful testground for more sophisticated microscopic theories, in which heavy numerical work is inevitable.

An algebraic approach to the structure of metal clusters  has been introduced recently , taking advantage of the three-dimensional $`q`$-deformed harmonic oscillator (3D $`q`$-HO) , which is constructed using the techniques of quantum algebras (quantum groups) \[see  and references therein\]. The symmetry of this oscillator is uq(3)$``$soq(3) , while its derivation involves use of irreducible tensor operators under soq(3) . In terms of the 3D $`q`$-HO a good description of the magic numbers of alkali clusters (up to 1500 atoms, which is the limit of validity for theories based on the filling of electronic shells ), as well as of Al clusters, has been obtained . In addition, it has been proved that supershells occur naturally in this model , which is characterized by only one parameter (the deformation parameter $`q=e^\tau `$, with $`\tau `$ being real) in addition to the overall scale.

An alternative approach to shell structure coming from mathematical physics has been introduced long before the discovery of metal clusters  by Balian and Bloch , in the framework of the theory of periodic orbits and classical quantization conditions . In this approach the valence electrons of the metal cluster are supposed to move in a spherical cavity (representing the effects of the mean field). Supershells then occur from the superposition of closed classical stationary orbits  of the electrons, the superposition of triangular and square orbits being the simplest example . The theory of Balian and Bloch leads to very specific predictions about the behaviour of magic numbers, which can serve as a test for alternative approaches. For example, the plot of $`N_i^{1/3}`$ vs. $`i`$, where $`N_i`$ are the magic numbers of clusters of a specific alkali metal and $`i`$ their index ($`i=1`$, 2, 3, …), should be a straight line with a slope of 0.61 .

In this paper we confront the results provided by the 3D $`q`$-HO to the restrictions imposed by the theory of Balian and Bloch. It turns out that agreement between the predictions of the two theories can be established by introducing in the Hamiltonian of the 3D $`q`$-HO an additional term, which is characterized by the same symmetry as the original Hamiltonian. Furthermore, it can be seen that this additional term occurs naturally through a variational procedure, analogous to the one used for obtaining the Morse oscillator  from the usual harmonic oscillator.

In Section 2 of this paper a brief account of the 3D $`q`$-HO is given, while in Section 3 the application of the periodic orbit theory of Balian and Bloch to metal clusters is briefly described. The predictions of the 3D $`q`$-HO are compared to the restrictions imposed by the theory of Balian and Bloch in Section 4, while in Section 5 a modified Hamiltonian for the 3D $`q`$-HO is introduced, allowing for full agreement with the theory of Balian and Bloch. Numerical details concerning the modified 3D $`q`$-HO Hamiltonian are given in Section 6, together with a study of supershells in the framework of this Hamiltonian. In Section 7 a variational method leading from the usual harmonic oscillator to the Morse oscillator is introduced, while in Section 8 this method is applied for deriving the modified Hamiltonian introduced in Section 5 from the original 3D $`q`$-HO Hamiltonian. Finally in Section 9 a discussion of the present results and plans for future work are given.

## 2 The 3-dimensional $`q`$-deformed harmonic oscillator (3D $`q`$-HO)

The space of the 3-dimensional $`q`$-deformed harmonic oscillator consists of the completely symmetric irreducible representations of the quantum algebra uq(3) . In this space a deformed angular momentum algebra, soq(3), can be defined . The Hamiltonian of the 3D $`q`$-HO is defined so that it satisfies the following requirements:

a) It is an soq(3) scalar, i.e. the energy is simultaneously measurable with the $`q`$-deformed angular momentum related to the algebra soq(3) and its $`z`$-projection.

b) It conserves the number of bosons, in terms of which the quantum algebras uq(3) and soq(3) are realized.

c) In the limit $`q1`$ it is in agreement with the Hamiltonian of the usual 3-dimensional harmonic oscillator.

It has been proved  that a Hamiltonian of the 3D $`q`$-HO satisfying the above requirements takes the form

$$H_q=\mathrm{}\omega _0\left\{[N]q^{N+1}\frac{q(qq^1)}{[2]}C_q^{(2)}\right\},$$
(1)
where $`N`$ is the number operator and $`C_q^{(2)}`$ is the second order Casimir operator of the algebra soq(3), while

$$[x]=\frac{q^xq^x}{qq^1}$$
(2)
is the definition of $`q`$-numbers and $`q`$-operators. In the framework of quantum algebras in general, the deformation parameter $`q`$ can either be real ($`q=e^\tau `$, with $`\tau `$ being real) or a phase factor ($`q=e^{i\tau }`$, with $`\tau `$ being real). In the present case, however, the parameter $`q`$ is restricted to obtain real values only, since this is required in the derivation leading to the Hamiltonian of Eq. (1). (See the comment following Eq. (9) in Ref.  for more details.)

The energy eigenvalues of the 3D $`q`$-HO are given by

$$E_q(n,l)=\mathrm{}\omega _0\left\{[n]q^{n+1}\frac{q(qq^1)}{[2]}[l][l+1]\right\}=\mathrm{}\omega _0e_q(n,l),$$
(3)
where $`n`$ is the number of vibrational quanta and $`l`$ is the eigenvalue of the angular momentum, obtaining the values $`l=n,n2,\mathrm{},0`$ or 1, and

$$e_q(n,l)=[n]q^{n+1}\frac{q(qq^1)}{[2]}[l][l+1].$$
(4)
In the limit of $`q1`$ one obtains $`\mathrm{lim}_{q1}E_q(n,l)=\mathrm{}\omega _0n`$, which coincides with the classical result.

For small values of the deformation parameter $`\tau `$ (where $`q=e^\tau `$) one can expand Eq. (3) in powers of $`\tau `$ obtaining

$$E_q(n,l)=\mathrm{}\omega _0n\mathrm{}\omega _0\tau \left\{l(l+1)n(n+1)\right\}$$

$$\mathrm{}\omega _0\tau ^2\left\{l(l+1)\frac{1}{3}n(n+1)(2n+1)\right\}+𝒪(\tau ^3).$$
(5)
The last expression to leading order bears great similarity to the modified harmonic oscillator suggested by Nilsson  (with the spin-orbit term omitted). Comparisons between the predictions of the 3D $`q`$-HO and Nilsson’s modified oscillator for the magic numbers and supershells of metal clusters have been given in . One of the main differences between the two models is that the term $`\tau n(n+1)`$ in the eigenvalues of the 3D $`q`$-HO given in Eq. (5) is occuring as a consequence of the overall symmetry, while the corresponding term $`\mu ^{}n(n+3)/2`$ in the eigenvalues of Nilsson’s modified oscillator (see, for example, Eq. (3) of Ref. ) is put in “by hand”. It should also be noticed that the use of $`q`$-deformations does not mean that an additional parameter is introduced in the theory. In both the 3D $`q`$-HO and Nilsson’s modified oscillator only one parameter appears (in addition to the overall scale). In the 3D $`q`$-HO case this is the deformation parameter $`\tau `$, while in Nilsson’s modified oscillator it is the parameter $`\mu ^{}`$ (see, for example, Eq. (3) of Ref. ).

## 3 The theory of Balian and Bloch

The theory of Balian and Bloch  is an example of periodic orbit theory , developed as a semiclassical bridge between quantum mechanics and classical mechanics. In the approach of Balian and Bloch, shell effects in metal clusters are studied by considering the valence electrons moving on straight lines within a smooth sphere and being reflected on the inner surface of the sphere. Periodic orbits in this case correspond to various polygons with three or more corners, the triangular and the square orbits being the simplest examples .

In the theory of Balian and Bloch a cluster of $`N`$ alkali atoms can be studied  by considering the $`N`$ valence electrons moving in a sphere with radius

$$R=r_SN^{1/3},$$
(6)
where $`r_S`$ is the Wigner–Seitz radius. The lengths of the triangular and the square orbits are then respectively

$$L_3=3\sqrt{3}r_SN^{1/3},$$
(7)
$$L_4=4\sqrt{2}r_SN^{1/3}.$$
(8)
Assuming that electrons move with the Fermi velocity $`v_F`$, the beating pattern created by the superposition of the two orbits leads to supershell structure. It turns out that the magic numbers (shell closures) $`N_i`$ are related to the index $`i`$ counting their number ($`i=1`$, 2, 3, …) by

$$N_i^{1/3}=\frac{h}{mv_Fr_S}\frac{2}{3\sqrt{3}+4\sqrt{2}}i=0.605i,$$
(9)
where $`m`$ is the electron mass and the fact that $`v_Fr_S`$ is constant for all metals has been used. Therefore the plot of $`N_i^{1/3}`$ vs. $`i`$ should be a straight line with a slope of 0.61 . Furthermore, a phase shift of a half unit of the index $`i`$ should occur in the node region of the beat pattern, i.e., when passing from a supershell to the next one . In addition, within each supershell the shell closures should appear at equidistant positions (i.e., they should exhibit a periodicity) when plotted vs. $`N^{1/3}`$ , where $`N`$ is the number of particles (valence electrons in the present case).

Comparisons to experimental sets of magic numbers have shown  that a slope of 0.61 is obtained in the cases of Na and Li clusters, revealing that the triangular and square orbits represent a good approximation in these cases. For Al clusters, though, a slope of 0.32 has been obtained , indicating that in this case, if the theory is applicable, more complicated orbits enter .

## 4 Comparing the predictions of the 3D $`q`$-HO and of the Balian–Bloch theory

As we have seen in the previous section, a hallmark of the theory of Balian and Bloch is that the plot of $`N_i^{1/3}`$ vs. $`i`$, where $`N_i`$ are the magic numbers and $`i`$ the index counting them ($`i=1`$, 2, 3, …), should be a straight line having a slope of 0.61 in the case of alkali metals, while in the case of Al clusters the slope should be 0.32 .

In order to compare the predictions of the 3D $`q`$-HO to the restrictions imposed by the theory of Balian and Bloch, we plot in Fig. 1(a) (line labelled by $`ϵ=0.0`$) the magic numbers obtained from the 3D $`q`$-HO for $`\tau =0.038`$, the parameter value found appropriate in Ref.  for reproducing the magic numbers of alkali clusters (up to 1500 atoms, which is the limit of validity of theories based on the filling of electronic shells ). Numbers considered as magic, listed in Table 1, correspond to energy gaps larger than $`\delta =0.38`$, as in Ref. , with $`\mathrm{}\omega _0=1`$. In the same figure, a straight line with a slope of 0.61 appears. It is clear that the 3D $`q`$-HO magic numbers follow the straight line up to $`i=14`$ quite well, while beyond this point the predictions of the 3D $`q`$-HO are clearly lower than the straight line, indicating that “too many” magic numbers are produced by the model in this region.

The same conclusion is arrived at by looking at Fig. 1(b), where the predictions of the 3D $`q`$-HO for $`\tau =0.050`$, the parameter value found in Ref.  appropriate for reproducing the magic numbers of Al clusters, are reported (line labelled by $`ϵ=0.0`$). Again magic numbers, listed in Table 1, correspond to energy gaps larger than $`\delta =0.38`$ (with $`\mathrm{}\omega _0=1`$). Small magic numbers (below 186, i.e., below $`i=9`$) are not shown, since it is known that small magic numbers in Al clusters cannot be explained by models based on the filling of electronic shells, because of the symmetry breaking caused by the ionic lattice , while for large magic numbers this problem does not exist. In the plot a straight line with a slope of 0.32, which is expected to be appropriate for Al clusters , as mentioned in the previous section, is also seen. It is clear that the predictions of the 3D $`q`$-HO follow a line parallel to the one with slope 0.32 roughly up to $`i=23`$, while beyond this point the slope is gradually reduced, indicating that the model predicts “too many” magic numbers.

## 5 A modified Hamiltonian for the 3D $`q`$-HO

The discrepancy between the 3D $`q`$-HO and the theory of Balian and Bloch can be lifted by considering the Hamiltonian

$$H_q^{}=H_qϵH_q^2,$$
(10)
with eigenvalues

$$E_q^{}(n,l)=E_q(n,l)ϵE_q^2(n,l),$$
(11)
where $`ϵ`$ is a small real positive constant. Justification for this choice will be given in Sections 7 and 8 through a variational procedure. For the moment the following comments suffice:

a) It is clear that $`H_q^{}`$ is a function of $`H_q`$, which is by construction an soq(3) scalar, as mentioned in Section 2. Therefore $`H_q^{}`$ is also an soq(3) scalar.

b) The energy eigenvalues of the new Hamiltonian can be written in the form

$$E_q^{}(n,l)=\mathrm{}\omega _0\left(1ϵ\mathrm{}\omega _0e_q(n,l)\right)e_q(n,l)=\mathrm{}\omega (n,l)e_q(n,l),$$
(12)
where

$$\omega (n,l)=\omega _0\left(1ϵ\mathrm{}\omega _0e_q(n,l)\right)$$
(13)
is a variable frequency, depending on the quantum numbers $`n`$, $`l`$, and on the small parameter $`ϵ`$. We shall call this oscillator the Variable Frequency Oscillator (VFO) corresponding to the 3D $`q`$-HO, a term for which justification will be provided in Sections 7 and 8.

The magic numbers provided by the VFO for a few appropriate values of $`ϵ`$ in the case of $`\tau =0.038`$, which is relevant for alkali clusters , are shown in Fig. 1(a) and listed in Table 1. Details of the calculation will be given in Section 6. Once more magic numbers are separated by gaps larger than $`\delta =0.38`$, while $`\mathrm{}\omega _0=1`$. It is clear that the predictions of the VFO roughly follow the straight line with a slope of 0.61 even for large values of $`i`$, thus overcoming the difficulties faced by the 3D $`q`$-HO.

A similar picture is obtained for $`\tau =0.050`$ (and $`\delta =0.38`$, with $`\mathrm{}\omega _0=1`$), which is appropriate for Al clusters . For a few appropriate values of $`ϵ`$, shown in Fig. 1(b) and listed in Table 1, the predictions of the VFO roughly follow a straight line with a slope of 0.32, as they should, according to the previous section.

We therefore conclude that the addition of the second term in the Hamiltonian of the 3D $`q`$-HO, leading to the VFO, makes the predictions of the 3D $`q`$-HO compatible with the predictions of the theory of Balian and Bloch.

## 6 Numerical details

In this section the calculations leading to the results reported in Section 5 will be described. Throughout this paper we put $`\mathrm{}\omega _0=1`$ for simplicity.

An important difference between the 3D $`q`$-HO and the VFO of Eq. (10) lies in the way truncations of the spectrum are made. The following comments apply:

a) In the 3D $`q`$-HO the level with $`l=n`$ always lies lowest in energy within each shell, the level with $`l=n2`$ lies immediately above it, and so on. Therefore stopping the level scheme at the $`l=n`$ level of a given shell and taking into account all levels with lower $`n`$ (i.e., all levels of the shells lying below the given one), one makes sure that all levels up to the given level have been included .

b) In the case of the VFO of Eq. (10) the following picture occurs: For a given (small) value of $`ϵ`$ the first several shells exhibit the same behaviour as in the case of the 3D $`q`$-HO, i.e., the level with $`l=n`$ lies lowest in energy within each shell, the level with $`l=n2`$ lies immediately above it, and so on. As the shell number is increasing, however, an inversion of the order of the levels occurs, with the levels with $`l=0`$ (for $`n`$ being even) or $`l=1`$ (for $`n`$ being odd) lying lowest in energy within the shell. It is then clear that beyond this inversion point truncation at a given shell should be made at the level with $`l=0`$ (for $`n`$ being even) or at $`l=1`$ (for $`n`$ being odd).

The magic numbers obtained for the values of $`\tau `$ and $`ϵ`$ used in the previous section are given in Table 1. In each case the maximum value of $`n`$, $`n_{max}`$, included in the calculation is indicated. Care has been taken that no inversion of the order of the levels, of the type described in comment b) above, occurs for the values of $`n`$ included in the calculation. Therefore in all cases truncation is made at the level with $`n=n_{max}`$ and $`l=n_{max}`$. The total number of levels up to the truncation point, $`N_{max}`$, is also shown in Table 1. The following remarks are now in place:

a) In the cases considered here and up to the truncation point, for a given value of $`\tau `$ the order of the levels is not modified as $`ϵ`$ is changed. The only modification occuring is that the spectrum gets “squeezed” as $`ϵ`$ increases.

b) As a result of a), the magic numbers reported in Table 1 present the following feature. For each value of $`\tau `$ and for $`ϵ=0.0`$ the largest number of magic numbers appears. As $`ϵ`$ increases, some of the magic numbers cease to be magic any more, since the “squeezing” of the spectrum brings the levels closer to each other. For each value of $`ϵ`$ the magic numbers occuring are a subset of the magic numbers occuring for lower values of $`ϵ`$ (with the same $`\tau `$). No new magic numbers appear, within the limits considered here, as $`ϵ`$ increases.

It is interesting to examine at this point what the influence of the additional term to the appearance of supershells is. For this purpose we are going to use the procedure employed by Nishioka et al. . For a given number of particles $`N`$ the single particle energies $`E_j(n,l)`$ of the $`N`$ occupied states are summed up

$$E(N)=\underset{j=1}{\overset{N}{}}E_j(n,l).$$
(14)
This sum is then divided into two parts: A smooth average part $`E_{av}`$ and a shell part $`E_{shell}`$, which will exhibit the supeshell structure

$$E(N)=E_{av}(N)+E_{shell}(N).$$
(15)
For the average part of the total energy a Liquid Drop Model expansion is used

$$E_{av}(N)=a_1N^{1/3}+a_2N^{2/3}+a_3N+a_4N^{4/3}+a_5N^{5/3}+a_6N^2.$$
(16)
The parameters of the fits occuring in the cases considered here are shown in Table 2, together with the number of levels, $`N_{max}`$, included in the fit and the rms deviation $`\sigma `$. We remark that for a given value of $`\tau `$ the parameters change smoothly for the nonzero values of $`ϵ`$, while the case with $`ϵ=0.0`$ is characterized by quite different values of the parameters, but also by a higher rms deviation $`\sigma `$. The addition of the second term in Eq. (10) improves the agreement of the average part of the total energy to the Liquid Drop Model expansion, thus resulting in lower rms deviations $`\sigma `$.

The procedure of the calculation was as follows: First the summations described by Eq. (14), resulting in the total energy $`E(N)`$ for each particle number $`N`$, have been performed. Subsequently, in order to reduce the size of the calculation approximately by a factor of 10, the average $`E(N)`$ was calculated every 11 points (i.e., for $`N=6`$, 17, 28, …) up to the cutoff point which is reported in Table 2 as $`N_{max}`$. These averaged values of $`E(N)`$ were subsequently fitted by the expansion of Eq. (16), resulting in the determination of $`E_{av}(N)`$ at these points. Finally $`E_{shell}(N)`$ has been obtained at these points as the difference $`E(N)E_{av}(N)`$ and plotted in Figs. 2 and 3.

The shell energy, $`E_{shell}`$, is plotted vs. the particle number, $`N`$, for $`\tau =0.038`$, which is appropriate for alkali clusters , in Fig. 2. In Fig. 2(a) the predictions of the original 3D $`q`$-HO (with $`ϵ=0`$) are shown, while in Fig. 2(b) the results of the VFO with $`ϵ=0.006`$ are depicted. The VFO with $`ϵ=0.007`$, 0.008 gives results which look almost identical with Fig. 2(b) and therefore are not shown for brevity. The similarities between Figs. 2(a) and 2(b) are clear. Not only the supershell appears in both cases around $`N=1000`$, as it is expected for Na clusters , but in addition even the maxima and minima of the shell energy appear at the same particle numbers and have roughly the same magnitude. Even the local maxima and minima present striking similarities. These results corroborate the remarks made above, namely that the addition of the second term to the Hamiltonian of Eq. (10) does not influence the order of the energy levels, the main effect of the second term being the gradual “squeezing” of the spectrum as energy increases. Of course this conclusion is valid only within the region of particle numbers studied and for small values of $`ϵ`$, like the ones used here.

Similar results are obtained in Fig. 3 for the case of $`\tau =0.050`$, which is appropriate for Al clusters . Besides the $`ϵ=0`$ case, shown in Fig. 3(a), the results corresponding to $`ϵ=0.0050`$ are shown in Fig. 3(b), since the cases with $`ϵ=0.0053`$, 0.0055 provide results almost identical with the ones shown in Fig. 3(b). In all cases there is some evidence for a supershell below $`N=1000`$, although its appearance is not as clear as in the case of Fig. 2. The appearance of a supershell in this region is in agreement with the results of more sophisticated calculations, as, for example, spherical jellium model predictions in Local Density Approximation , but it is not in good agreement with experiment, where no evidence for supershell in Al clusters exists in this region . An advantage of the VFO in comparison to earlier calculations  is that at least it can reproduce the slope of 0.32 in Fig. 1(b), something which is not occuring in spherical jellium model calculations, although it occurs experimentally .

In the theory of Balian and Bloch, as mentioned above, a phase shift by a half unit of the running index $`i`$ should be observed in the plot of the magic numbers $`N_i^{1/3}`$ vs. $`i`$ when passing from a supershell to the next . In Fig. 1(a) such a shift is seen quite clearly around $`N_i=1000`$ in the cases of $`ϵ=0.006`$, 0.007, 0.008, while no clear shift of this type is seen in Fig. 1(b).

Furthermore, the theory of Balian and Bloch, as mentioned above, predicts that within each supershell the minima of the shell energy, $`E_{shell}`$, should appear at equidistant positions (i.e. they should exhibit a periodicity) when plotted vs. $`N^{1/3}`$. The change in the periodicity when passing from the first supershell to the second one is clear in Figs. 2 and 3 (although in these cases, for reasons of clarity, $`E_{shell}`$ is plotted vs. $`N`$ and not vs. $`N^{1/3}`$). The fact that the predictions of the VFO corresponding to the 3D $`q`$-HO approximately exhibit the right periodicity features does not come as a surprise, since the 3D $`q`$-HO is known to show this feature , while, as we have seen above, the addition of the second term in Eq. (10) does not influence the position of the minima.

In conclusion, the VFO corresponding to the 3D $`q`$-HO is able to reproduce the right slope in the $`N_i^{1/3}`$ vs. $`i`$ plot in both the alkali and Al clusters. In addition it predicts correctly the first supershell in alkali clusters, while in Al clusters its prediction for a supershell is in rough agreement with results of spherical jellium models but not with experiment.

## 7 A variational method

In nuclear physics it is well known that nuclear spectra can be described very accurately in terms of the Variable Moment of Inertia (VMI) model . In this model the energy levels are given by

$$E(J)=\frac{J(J+1)}{2\mathrm{\Theta }(J)}+\frac{1}{2}C\{\mathrm{\Theta }(J)\mathrm{\Theta }(0)\}^2,$$
(17)
where $`J`$ is the angular momentum and $`\mathrm{\Theta }(J)`$ is the moment of inertia, which is supposed to be a function of the angular momentum. $`C`$ and $`\mathrm{\Theta }(0)`$ are free parameters, the latter representing the ground state moment of inertia. It is clear that the VMI formula is a generalization of the rigid rotator formula

$$E(J)=\frac{J(J+1)}{2\mathrm{\Theta }},$$
(18)
in which the moment of inertia is assumed to be constant. The rigid rotator formula is known to fail beyond the first few levels of a rotational nucleus, since the experimental levels appear “squeezed” in comparison to the rigid rotator predictions. This difficulty is overcome in the framework of the VMI model by determining the moment of inertia for each value of the angular momentum $`J`$ through a minimization of the energy with respect to the moment of inertia for given angular momentum

$$\frac{E(J)}{\mathrm{\Theta }(J)}|_J=0.$$
(19)
This variational condition leads to a cubic equation for $`\mathrm{\Theta }(J)`$, which turns out to have only one real solution , corresponding to the appropriate value of the moment of inertia for the given value of the angular momentum. The second term in Eq. (17) is justified by the well known fact that many perturbing potentials near their origin can be approximated by a harmonic oscillator potential.

Following the same reasoning, it is interesting to examine what happens to the usual harmonic oscillator if, by analogy, one allows the angular frequency to be a function of the quantum number $`n`$. The energy will then read

$$E(n)=\mathrm{}\omega (n)\left(n+\frac{1}{2}\right)+\frac{1}{2}C\{\omega (n)\omega (0)\}^2,$$
(20)
where $`C`$ and $`\omega (0)`$ are free parameters, the latter corresponding to the ground state angular frequency. The variational principle in this case should correspond to the minimization of the energy with respect to the angular frequency for constant value of the quantum number $`n`$

$$\frac{E(n)}{\omega (n)}|_n=0.$$
(21)
It is clear that this condition leads to

$$\mathrm{}\left(n+\frac{1}{2}\right)+C\{\omega (n)\omega (0)\}=0\omega (n)=\omega (0)\frac{\mathrm{}}{C}\left(n+\frac{1}{2}\right).$$
(22)
Substituting this result in Eq. (20) we obtain

$$E(n)=\mathrm{}\omega (0)\left(n+\frac{1}{2}\right)\frac{1}{2}\frac{\mathrm{}^2}{C}\left(n+\frac{1}{2}\right)^2,$$
(23)
which is reminiscent of the spectrum of the Morse potential .

Indeed, solving the Schrödinger equation for the Morse potential

$$V(x)=D(1e^{\alpha x})^2,$$
(24)
one obtains the energy spectrum

$$E(n)=\mathrm{}\omega \left\{\left(n+\frac{1}{2}\right)x_e\left(n+\frac{1}{2}\right)^2\right\},$$
(25)
where

$$x_e=\frac{1}{2}\frac{\mathrm{}\alpha }{\sqrt{2mD}},$$
(26)
and

$$\omega =\alpha \sqrt{\frac{2D}{m}}.$$
(27)
We therefore conclude that by allowing the angular frequency of the simple harmonic oscillator to vary with the quantum number $`n`$, we obtain the spectrum of the Morse oscillator. Again the second term in Eq. (20) is in agreement to the fact that most perturbing potentials near their origin can be approximated by the harmonic oscillator potential. We shall refer to the oscillator of Eq. (20) as the Variable Frequency Oscillator (VFO).

## 8 Derivation of the modified Hamiltonian for the 3D $`q`$-HO
through a variational method

The idea leading to the VFO of the previous section can be appropriately generalized in the case of the 3D $`q`$-HO. In this case we consider the energy expression

$$E_q^{}(n,l)=\mathrm{}\omega (n,l)\left\{[n]q^{n+1}\frac{q(qq^1)}{[2]}[l][l+1]\right\}+\frac{1}{2}C\{\omega (n,l)\omega (0,0)\}^2$$

$$=\mathrm{}\omega (n,l)e_q(n,l)+\frac{1}{2}C\{\omega (n,l)\omega (0,0)\}^2,$$
(28)
where the angular frequency $`\omega (n,l)`$ depends on the quantum numbers $`n`$ and $`l`$, while $`C`$ and $`\omega (0,0)`$ are real positive constants, the latter corresponding to the ground state angular frequency, since the ground state of the 3D $`q`$-HO is characterized by $`n=0`$ and $`l=0`$. It is thus clear that $`\omega (0,0)`$ corresponds to $`\omega _0`$ appearing in Eq. (3), i.e., $`\omega (0,0)\omega _0`$. The variational condition in the present case should read

$$\frac{E_q^{}(n,l)}{\omega (n,l)}|_{n,l}=0.$$
(29)
In other words, the energy is minimized with respect to the angular frequency for constant values of the quantum numbers $`n`$ and $`l`$. The variational condition leads to the equation

$$\mathrm{}e_q(n,l)+C\{\omega (n,l)\omega (0,0)\}=0\omega (n,l)=\omega (0,0)\frac{\mathrm{}}{C}e_q(n,l)$$

$$=\omega (0,0)\frac{\mathrm{}}{C}\left\{[n]q^{n+1}\frac{q(qq^1)}{[2]}[l][l+1]\right\}.$$
(30)
Substituting this result in Eq. (28) one then obtains

$$E_q^{}(n,l)=\mathrm{}\omega (0,0)e_q(n,l)\frac{1}{2}\frac{\mathrm{}^2}{C}e_q^2(n,l),$$
(31)
which is the same as Eq. (11), with

$$ϵ=\frac{1}{2C\omega _0^2},$$
(32)
since $`\omega _0\omega (0,0)`$, as mentioned above.

## 9 Discussion

In this paper we have attempted a comparison of the predictions for the shell structure of metal clusters of the 3D $`q`$-HO model to the ones of the periodic orbit theory of Balian and Bloch. It turns out that the predictions of the 3D $`q`$-HO for the magic numbers of metal clusters can be made compatible with the predictions of the theory of Balian and Bloch by adding to the 3D $`q`$-HO a symmetry-preserving correction term reminiscent of the anharmonicity term in the spectrum of the Morse potential, while this addition does not influence the predictions for the supershells. This extended expression for the 3D $`q`$-HO can be justified through a variational method, similar to the one used in the Variable Moment of Inertia (VMI) model of nuclear physics, leading to the concept of the Variable Frequency Oscillator (VFO), which gives promise of wider applicability.

Acknowledgements

One of the authors (PPR) acknowledges support from the Bulgarian Ministry of Science and Education under contracts $`\mathrm{\Phi }`$-415 and $`\mathrm{\Phi }`$-547.

Figure captions

Fig. 1 Cubic roots of the magic numbers $`N_i`$ plotted vs. the running index $`i`$ counting them. The magic numbers are the ones listed in Table 1. (a) $`\tau =0.038`$, (b) $`\tau =0.050`$.

Fig. 2 Shell part ($`E_{shell}`$) of the total energy \[in units of $`\mathrm{}\omega _0`$, see Eqs. (3) and (11)\] vs. the number of particles $`N`$, in the case of $`\tau =0.038`$, (a) for the 3D $`q`$-HO, (b) for a corresponding VFO. The values of the dimensionless parameters $`\tau `$ and $`ϵ`$ are listed in Table 2, together with the details of the calculation. See Section 6 for further discussion.

Fig. 3 Same as Fig. 2, but for $`\tau =0.050`$.


format: text

text:
# 1 Introduction

## 1 Introduction

The evaluation of Casimir energies  associated with confined zero-point fluctuations of quantum fields depends on several factors, the more relevant ones are the nature of the quantum field, the type of space-time manifold and its dimensionality, and the type of boundary conditions imposed on the fields. All these factors taken together can lead to relatively simple dispersion relations, e.g.: a massless scalar field under Dirichlet b.c., or to rather complex ones, e.g.: a massive fermion field under MIT b.c.. Recently a very simple, but rather powerful way of calculating regularized Casimir energies suitable for non-trivial frequency spectra was introduced in the literature . This technique relies on well-known theorems of complex analysis, to wit: the Cauchy integral formula and the Mittag-Leffler expansion theorem in one of its simplest versions. The method was successfully employed in the evaluation of the Casimir energy of a massive fermionic field confined in a $`d+1`$-dimensional slab-bag  and also in the evaluation of the Casimir energy of confined bosonic and fermionic fields in the presence of pertinent b.c. and a uniform magnetic field . Here we will sketch the application of this method to the evaluation of another example of a sum over a non-trivial frequency spectrum, namely: the one associated with the regularized Casimir energy of a kappa-deformed scalar field.

## 2 A sum formula

Consider a quantum field in a $`d+1`$-dimensional flat space-time under b.c. imposed on two hyperplanes of area $`L^{d1}`$. The distance between the hyperplanes is $`\mathrm{}`$, hence there is a restriction on one of the spatial dimensions, the $`OX_d`$ axis. Suppose $`L\mathrm{}`$. At the one-loop level, the (unregularized) Casimir energy is given by

$$E_0\left(d\right)=\alpha \left(d\right)\frac{L^{d1}}{2}\underset{n}{}\frac{d^{d1}}{\left(2\pi \right)^{d1}}\mathrm{\Omega }_n,$$
(1)
where $`\alpha \left(d\right)`$ is a dimensionless factor that takes into account the internal degrees of freedom of the field and

$$\mathrm{\Omega }_n=\sqrt{p_{}^2+\frac{\lambda _n^2}{\mathrm{}^2}+m^2},$$
(2)
where $`p_{}^2=p_1^2+p_2^2+\mathrm{}+p_{d1}^2`$, $`m`$ is the mass of an elementary excitation of the field and $`\lambda _n`$ is the $`n`$th real root of the transcendental equation determined by the b. c.. It is shown in Ref.  that Cauchy integral formula and Mittag-Leffler expansion theorem carefully applied allow us to recast the sum above into the form

$$E_0\left(d\right)=\alpha \left(d\right)\frac{L^{d1}}{2}\frac{d^dp}{\left(2\pi \right)^d}\mathrm{log}\left[1+\frac{K_1\left(z\right)}{K_2\left(z\right)}\right],$$
(3)
where $`z`$ is defined by

$$q^2+\mathrm{\Omega }_n^2=\frac{z^2+\lambda _n^2}{\mathrm{}^2},$$
(4)
with $`q`$ as an auxiliary complex momentum. The functions $`K_1\left(z\right)`$ and $`K_2\left(z\right)`$ must be constructed from the boundary conditions. This last relation gives the regularized Casimir energy. For details see , and see also  for a simple example. It can be shown  that this method can be extended to the case where the sum we must evaluate is performed over a function of $`\mathrm{\Omega }_n`$, i.e.:

$$E_0\left(d\right)=\alpha \left(d\right)\frac{L^{d1}}{2}\underset{n}{}\frac{d^{d1}}{\left(2\pi \right)^{d1}}f\left(\mathrm{\Omega }_n\right).$$
(5)
This is precisely what happens when we consider a kappa-deformed scalar theory. Kappa-deformed theories are theories whose associated symmetry groups depart from the usual symmetry groups of standard quantum field theories, e.g.: the Poincaré group. This departure is measured by some convenient parameter, say $`\kappa `$.

## 3 The Casimir effect for kappa-deformed scalar field theory

The kappa-deformed Klein-Gordon equation in $`3+1`$ dimensions is given by (see  and references therein)

$$\left[\stackrel{}{}^2\left(2\kappa \mathrm{sinh}\left(\frac{P_0}{2\kappa }\right)\right)^2m^2\right]\varphi (\stackrel{}{x},x_0)=0,$$
(6)
where $`\kappa `$ is a parameter with dimensions of a mass. The on shell condition reads

$$\stackrel{}{P}^{\mathrm{\hspace{0.17em}2}}\left(2\kappa \mathrm{sinh}\left(\frac{P_0}{2\kappa }\right)\right)^2=m^2,$$
(7)
Upon diagonalizing simultaneously the operators $`\stackrel{}{P}`$ and $`P_0`$ we obtain a new frequency spectrum which is given by

$$\omega \left(\stackrel{}{p}\right)=2\kappa \mathrm{sinh}^1\left(\frac{1}{2\kappa }\sqrt{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+m^2}\right)$$
(8)
Suppose Dirichlet boundary conditions are imposed on two parallel surfaces perpendicular to the $`OX_3`$-axis. Then

$$E_0\left(\mathrm{}\right)=\frac{L^2}{\eta }\frac{d^2p_{}}{\left(2\pi \right)^2}\underset{n}{}\mathrm{sinh}^1\left(\eta \mathrm{\Omega }_n\right),$$
(9)
where for convenience we have defined $`\eta :=1/2\kappa `$. By applying the techniques described in  to this more complex case we obtain

$$E_0(\mathrm{})=\frac{L^2}{\eta }\frac{d^2p_{}}{(2\pi )^3}𝑑q\underset{n}{}\frac{2q^2}{q^2+\mathrm{\Omega }_n^2}\frac{\mathrm{sin}^1(\eta q)}{q},$$
(10)
where as before $`q`$ plays the role of an auxiliary momentum. The pole structure of the above equation is the same as in the absence of deformation, moreover, we can consider the boundary conditions as being the same. In our case this will mean Dirichlet boundary conditions applied to the undeformed case.

It follows from Eq. (10) that

$$E_0(\mathrm{})=L^2\frac{d^2p_{}}{(2\pi )^3}_0^{\frac{1}{\eta }}\frac{dq}{\sqrt{1\eta ^2q^2}}\mathrm{log}\left[1+\frac{K_1\left(z\right)}{K_2\left(z\right)}\right]$$
(11)
where $`K_1\left(z\right)`$ and $`K_2\left(z\right)`$ are two functions that can be read out from the boundary conditions . For Dirichlet b.c. we have $`K_1\left(z\right)=e^z/2`$ and $`K_2\left(z\right)=e^z/2`$, with $`z=\mathrm{}\sqrt{\stackrel{}{p}_{}^{\mathrm{\hspace{0.17em}2}}+q^2}`$ . Taking these two functions into Eq. (11) we obtain

$$E_0\left(\mathrm{}\right)=\frac{L^2}{8\pi ^3}d^2p_{}_0^{\frac{1}{\eta }}\frac{dq}{\sqrt{1\eta ^2q^2}}\mathrm{log}\left[1e^{2\mathrm{}\sqrt{\stackrel{}{p}_{}^{\mathrm{\hspace{0.17em}2}}+q^2}}\right],$$
(12)
which is in agreement with . If we take the limit $`\eta 0`$ after some manipulations we obtain

$$E_0\left(\mathrm{}\right)=\frac{L^2}{2\pi ^2}\mathrm{}^3_0^{\mathrm{}}𝑑zz^2\mathrm{log}\left[1e^{2z}\right],$$
(13)
which leads to the original Casimir energy for this situation .

## 4 Final remarks

In this brief note we described an alternative method for summing over non-trivial frequency spectra and applied the main result to the problem of evaluating the Casimir energy of a kappa-deformed scalar theory under Dirichlet boundary conditions. The results obtained are consistent with the ones obtained by other methods and hold for quantum electrodynamics between parallel plates. This method can be also applied to other non-trivial examples. A more detailed account of the extension of the technique employed in the example discussed here is in preparation and the results will be published elswhere.

## Acknowledgments

D. Pinheiro acknowledges the financial support of the Instituto de Física, UFRJ, graduate program; F. C. Santos and A. C. Tort acknowledge the financial support of XXIV ENFPC Organizing Comitee that allowed them to be part of the XXIV Brazilian National Meeting on Particles and Fields held in Caxambu MG, Brazil, September 30 – October 4, 2003.


text:
# Andrade, Omori and Time-to-failure Laws from Thermal Noise in Material Rupture

## Abstract

Using the simplest possible ingredients of a rupture model with thermal fluctuations, we provide an analytical theory of three ubiquitous empirical observations obtained in creep (constant applied stress) experiments: the initial Andrade-like and Omori-like $`1/t`$ decay of the rate of deformation and of fiber ruptures and the $`1/(t_ct)`$ critical time-to-failure behavior of acoustic emissions just prior to the macroscopic rupture. The lifetime of the material is controlled by a thermally activated Arrhenius nucleation process, describing the cross-over between these two regimes. Our results give further credit to the idea proposed by Ciliberto et al. that the tiny thermal fluctuations may actually play an essential role in macroscopic deformation and rupture processes at room temperature. We discover a new re-entrant effect of the lifetime as a function of quenched disorder amplitude.

preprint: APS preprint

Constant stress (so-called “creep”) experiments constitute a standard testing procedure in material sciences. The typical response to the sudden application of a constant stress is that the strain rate as well as the acoustic emission rate first jump to high values followed by slow universal power law decays, respectively called the Andrade law Andrade  for the strain rate and the Omori law Omori  for the acoustic rate. Then, after a long decay whose duration may vary within extraordinary large bounds (see below), the rates rebound and accelerate (while the applied stress remains constant) by following a power law acceleration resulting in a finite-time singularity (the rupture of the sample). The two regimes of decelerating followed by accelerating rates and the lifetime of the structure are the result of a subtle interplay between the pre-existing micro-heterogeneity of the material and the self-organized evolving deformation and damage due to dislocation motion and/or micro-cracking. Up to now, there are no theory encompassing all these regimes. Here, we propose a simple mechanism that provides an explanation of all these observations, which is based on the recent proposal cilietal1 ; Politi  that thermal noise is strongly renormalized by quenched heterogeneities. Based on the analysis of a simple fiber bundle rupture model, Refs.cilietal1 ; Politi  showed that the average lifetime of the fiber-bundle takes an Arrhenius form with an effective temperature renormalized from the bare temperature $`T`$ to a value strongly amplified by the presence of the frozen disorder in the rupture thresholds $`f_c(i)`$, in agreement with experiments and numerical simulations. This result suggests that the usual assumption of neglecting the role of thermal fluctuations in material rupture processes at room temperature may actually be incorrect (see early  for early discussions): due to frozen heterogeneities, tiny thermal fluctuations can be amplified many times, thus actually controlling the time-dependent aspects of failure. Our purpose is to extend the analysis of this model by showing that it is able to reproduce all the empirical observations mentioned above in creep (constant applied stress) experiments. This shows that the simplest possible ingredients of a rupture model together with thermal fluctuations can render essentially all of the richness of creep experiments.

The democratic fiber-bundle model (DFBM) with thermal noise cilietal1 ; Politi  can be seen as a mean field treatment of rupture. A macroscopic constant load $`F=Nf_0`$ is applied at time $`t=0`$ to an initially undamaged system made of a very large number $`N`$ of parallel elastic fibers (the results derived below are obtained in the thermodynamic limit $`N\mathrm{}`$). At all times, $`F`$ is shared democratically among all $`(1\mathrm{\Phi }(t))N`$ surviving fibers, where $`\mathrm{\Phi }(t)`$ is the fraction of broken fibers at time $`t`$. The externally applied force per surviving fiber is thus

$$f_a=\frac{f_0}{1\mathrm{\Phi }(t)}.$$
(1)
The strength of each fiber $`i`$ is characterized by a critical value $`f_c(i)`$ drawn for a distribution $`P_d(f)`$, centered on the mean equal to $`1`$ and with variance $`T_d`$. Putting the mean strength to $`1`$ sets the force scale. The heterogeneous strengths are constant and determined once for all, corresponding to a frozen disorder, which is “read” in a certain organized way as the rupture develops. Microscopic thermal fluctuations are taken into account by assuming that a fiber with load $`f_a`$ and threshold $`f_c(i)>f_a`$ has a non-zero probability $`G(f_c(i)f_a)`$ to rupture per unit time governed by the rate with which a thermal fluctuation can activate a microscopic force $`\mathrm{\Delta }f_if_c(i)f_a`$ to pass the rupture threshold $`f_c(i)`$: $`G(f_c(i)f_a)=\frac{\gamma }{2}\text{erfc}\left(\frac{f_c(i)f_a}{\sqrt{2T}}\right)`$, where $`\text{erfc}(x)`$ is the complementary error function, $`T`$ is the variance of the thermal force fluctuations $`\mathrm{\Delta }f_i`$ and $`\gamma `$ is a microscopic constant rate fixing the time scale of the thermal activation process. This expression amounts to introducing a zero-mean normal distribution of thermal fluctuation forces $`\mathrm{\Delta }f_i`$ with variance $`T`$ and with correlation time proportional to $`1/\gamma `$.

We first follow Politi  and introduce the distribution $`Q(f,t)`$ of the rupture thresholds of the unbroken fibers at time $`t`$. Obviously, $`Q(f,t)=0`$ for $`f<f_a`$, since all these fibers are already broken. Politi et al. Politi  have shown that $`Q(f,t)`$ can be approximated with a very high accuracy in the limit $`N\mathrm{}`$ by the initial distribution $`P_d(f)`$ of rupture strengths truncated at a lower value $`f_s(t)`$,

$$Q(f,t)=P_d(f),\mathrm{for}f>f_s(t)$$
(2)
and $`0`$ otherwise, where $`f_s(t)`$ is determined by the self-consistent equation

$$\mathrm{\Phi }(t)=_{\mathrm{}}^{f_s(t)}𝑑fP_d(f)$$
(3)
expressing that all fibers whose strengths are below $`f_s(t)`$ have failed at some time before $`t`$. This approximation for $`Q(f,t)`$ with (3) amounts to view the time-dependent rupture as a “front” propagating and “eating” the distribution $`P_d(f)`$ from the weakest towards the strongest fibers. We also have by definition $`\mathrm{\Phi }(t)=1_{\mathrm{}}^+\mathrm{}𝑑fQ(f,t)`$. Taking the time derivative of $`\mathrm{\Phi }(t)`$ and replacing $`\dot{Q}(f,t)`$ by $`Q(f,t)G(ff_a)`$ expressing that the rate of breaking is controlled the thermally activated rupture process acting on each fiber independently, we get $`\dot{\mathrm{\Phi }}=_{\mathrm{}}^+\mathrm{}𝑑fQ(f,t)G(ff_a)`$. Putting (2) in this equation and taking for $`P_d(f)`$ a normal distribution centered on $`1`$ with variance $`T_d`$ as in cilietal1 ; Politi  yields

$$\dot{\mathrm{\Phi }}=\frac{\gamma }{2}_{f_s}^{\mathrm{}}\frac{1}{\sqrt{2\pi T_d}}\mathrm{exp}\left[\frac{(1f)^2}{2T_d}\right]\text{erfc}\left(\frac{ff_a}{\sqrt{2T}}\right)𝑑f.$$
(4)
Making explicit $`P_d(f)`$ in (3) gives

$$\mathrm{\Phi }=\frac{1}{2}\left[\text{erf}\left(\frac{f_s1}{\sqrt{2T_d}}\right)+1\right],f_s=1+\sqrt{2T_d}\text{irf }(2\mathrm{\Phi }1),$$
(5)
where $`y=\text{irf }(z)`$ is the inverse function to the error function $`z=\text{erf }(y)`$. Putting $`f_s`$ in (4) gives

$$\dot{\mathrm{\Phi }}=R(\varphi )\frac{\gamma }{2}_\mathrm{\Phi }^1\text{erfc}\left[L(\mathrm{\Phi },z)\right]𝑑z,$$
(6)
$$L(\mathrm{\Phi },z)=\frac{1}{\sqrt{2T}}\left(1\frac{f_0}{1\mathrm{\Phi }}\right)+\mu \text{irf }(2z1),$$
(7)
with $`\mu =\sqrt{T_d/T}`$. The solution of equation (6) with (7) provides in principle all the information on the fraction $`\mathrm{\Phi }(t)`$ of broken fibers. This equation (7) is valid as long as the approximation (2) holds (see below).

The parameter $`\mu `$ quantifies the relative importance of the thermal fluctuations compared with the quenched heterogeneities. The relevant regime for applications to macroscopic ruptures at room temperature is $`\mu >1`$ and often $`\mu 1`$, that is, thermal fluctuations are tiny contributions to the applied macroscopic mechanical forces. Indeed, assuming that the energy barrier to rupture a fiber corresponds to the Griffith energy $`gc^2`$ necessary for nucleating a crack of half-length $`c`$ in the solid with surface energy $`g`$, we obtain $`\mu 1.5410^3`$ for $`c=1`$ micron and $`\mu 1.54`$ for $`c=1`$ nanometer, using $`g=1050`$ erg/cm2 for most solids. Thus, even for the smallest microcracks, thermal fluctuations are very small in relative value.

It turns out that this regime $`\mu 1`$ allows for a very convenient approximation of $`R(\mathrm{\Phi })`$ obtained by linearizing $`L(\mathrm{\Phi },z)`$ with respect to $`z`$. Then, the integral over $`z`$ in (6) can be calculated explicitly to yield

$$\dot{\mathrm{\Phi }}=R(\mathrm{\Phi })=\frac{\gamma T}{4\pi \mu D(\mathrm{\Phi })U(\mathrm{\Phi })}e^{U(\mathrm{\Phi })/T},$$
(8)
where

$$U(\mathrm{\Phi })=TL^2(\mathrm{\Phi },\mathrm{\Phi })=\frac{1}{2}\left[f_s(\mathrm{\Phi })f_a(\mathrm{\Phi })\right]^2,$$
(9)
and $`D(\mathrm{\Phi })=(1/\sqrt{2\pi T_d})df_s(\mathrm{\Phi })/d\mathrm{\Phi }=\mathrm{exp}\{\text{irf }^2(2\mathrm{\Phi }1)\}`$. Notice that $`U(\mathrm{\Phi })`$ in (9) has a clear physical interpretation. It is the energy barrier between the actual force $`f_a(\mathrm{\Phi })`$ (1) and the front value $`f_s(\mathrm{\Phi })`$ (5) of the distribution $`Q(f,t)`$ in (2). Equation (8) is valid as long as $`U(\mathrm{\Phi })T`$, which implies $`\mathrm{\Phi }<\mathrm{\Phi }_c`$, where $`\mathrm{\Phi }_c`$ is such that the force $`f_a`$ per surviving fiber reaches the average strength $`1`$: $`f_a(\mathrm{\Phi }_c)=1`$ yielding $`\mathrm{\Phi }_c=1f_0`$. Such fractions $`\mathrm{\Phi }\mathrm{\Phi }_c`$ correspond to the ultimate regime of explosive failure. We have checked by direct numerical calculations that (8) provides an exceedingly precise approximation of equation (6) as long as $`\mu 1`$ and not too close to $`\mathrm{\Phi }_c`$ (in practice to within a few percent).

To go further, we need to distinguish between two regimes, $`\mathrm{\Phi }<\mathrm{\Phi }^{}`$ and $`\mathrm{\Phi }>\mathrm{\Phi }^{}`$, where $`\mathrm{\Phi }^{}`$, solution of $`dR(\mathrm{\Phi })/d\mathrm{\Phi }=0`$, corresponds to the minimum failure rate. For $`\mu 1`$, $`\mathrm{\Phi }^{}`$ is actually independent of the temperature $`T`$ and is the root of the equation

$$D(\mathrm{\Phi }^{})(1\mathrm{\Phi }^{})^2=\alpha ,\alpha =\frac{f_0}{\sqrt{2\pi T_d}},$$
(10)
where $`\alpha `$ is a important physical parameter quantifying the strength of the disorder relative to the applied force. The explicit approximate solution of equation (10) is

$$\mathrm{\Phi }^{}(\alpha )=\{\begin{array}{cc}\frac{20\pi 4(4\pi )\alpha }{242\pi +8(4+\pi )\alpha },\hfill & \alpha <3/2,\hfill \\ \frac{1}{2}\text{erfc }\left(\sqrt{\mathrm{ln}\alpha }\right)\frac{\alpha \sqrt{\pi \mathrm{ln}\alpha }}{1+\alpha \sqrt{\pi \mathrm{ln}\alpha }},\hfill & \alpha >3/2.\hfill \end{array}$$
(11)
For example, this gives $`\mathrm{\Phi }^{}(\alpha =2)=0.089`$ compared with the exact value $`0.092`$.

It follows from (8) that the time to reach some $`\mathrm{\Phi }`$ is given by $`\gamma Tt=4\pi \mu _0^\mathrm{\Phi }D(z)U(z)e^{U(z)/T}𝑑z`$. For $`0<\mathrm{\Phi }<\mathrm{\Phi }^{}(\alpha )`$, due to the exponential factor $`e^{U/T}`$, the main contribution to the last integral comes from a small neighborhood of the upper integration limit. This yields

$$\gamma t4\pi \mu \frac{D(\mathrm{\Phi })U(\mathrm{\Phi })}{A(\mathrm{\Phi })}e^{U(\mathrm{\Phi })/T},A(\mathrm{\Phi })=\frac{dU(\mathrm{\Phi })}{d\mathrm{\Phi }},$$
(12)
for $`\mathrm{\Phi }<\mathrm{\Phi }^{}`$. This approximation is correct under the assumption that $`e^{U(\mathrm{\Phi })/T}`$ is rapidly increasing with $`\mathrm{\Phi }`$, i.e., if $`\left|A(\mathrm{\Phi })\right|T`$. The absolute value $`|\mathrm{}|`$ stresses that this condition applies also for $`\mathrm{\Phi }>\mathrm{\Phi }^{}(\alpha )`$. The same reasoning in this case gives a similar approximation

$$\gamma (t_ct)4\pi \mu \frac{D(\mathrm{\Phi })U(\mathrm{\Phi })}{\left|A(\mathrm{\Phi })\right|}e^{U(\mathrm{\Phi })/T},\mathrm{\Phi }>\mathrm{\Phi }^{},$$
(13)
where $`t_ct`$ is the time to complete rupture. The condition $`\left|A(\mathrm{\Phi })\right|T`$ shows that both relations (12) and (13) do not work in the vicinity of the minimum rate of fiber failures given by the solution of (10), for which $`A(\mathrm{\Phi }^{})=0`$. For $`0<\mathrm{\Phi }<\mathrm{\Phi }^{}`$, combining (8) and (12), we obtain $`\dot{U}=T/t`$ for $`t<t^{}`$, where $`t^{}`$ is defined by $`\mathrm{\Phi }(t^{})=\mathrm{\Phi }^{}`$. This gives

$$U\left[\mathrm{\Phi }(t)\right]=T\mathrm{ln}\gamma t(t<t^{})$$
(14)
valid for $`\gamma t1`$ (this condition simply means that the thermal fluctuations have had time to contribute several independent jolts). The constant of integration gives the $`\mathrm{ln}\gamma `$ contribution determined from matching with the initial stage. Replacing the lhs of (9) by (14) and putting $`f_af_0`$ (for $`\mathrm{\Phi }`$ small) gives, in view of equation (8), the fraction rate

$$\dot{\mathrm{\Phi }}\frac{1}{4\pi \mu t\mathrm{ln}\gamma t}\mathrm{exp}\left[\frac{1}{2T_d}\left(1f_0\sqrt{2T\mathrm{ln}\gamma t}\right)^2\right].$$
(15)
Expression (15) is one of our main results: the failure rate $`\dot{\mathrm{\Phi }}`$ of fibers decreases after application of the load proportionally to $`1/t`$, up to logarithm corrections. This $`1/t`$ decay lasts as long as $`\mathrm{\Phi }`$ remains smaller than $`\mathrm{\Phi }^{}`$. This $`1/t`$ law is known in seismology as the Omori law. It is also ubiquitous in creep experiments with exponents that are often close to or smaller than our prediction $`1`$. For intermediate times such that $`\gamma t<e^{1/2T}`$, $`\dot{\mathrm{\Phi }}\frac{1}{t\mathrm{ln}\gamma t}e^{(1f_0)\sqrt{2T\mathrm{ln}\gamma t}/T_d}`$, which gives an apparent exponent $`1/t^p`$ with $`p<1`$. For $`\mathrm{ln}\gamma t(1f_0)^2/2T`$, $`p1+(T/T_d)`$ which is close to but slightly larger than $`1`$. Numerical simulations confirm these predictions accurately. See for instance figure 2 of Politi  which our theory explains quantitatively. Exact numerical integration and our analytical approximation coincide everywhere, excluding a time interval corresponding to a very small vicinity of the stationary point $`\mathrm{\Phi }^{}`$. Note that Andrade’s law Andrade  also derives from the deformation rate being proportional to $`df_a(t)/dt=f_0\dot{\mathrm{\Phi }}/(1\mathrm{\Phi }(t))^2\dot{\mathrm{\Phi }}`$ as $`\mathrm{\Phi }(t)`$ varies much more slowly than $`\dot{\mathrm{\Phi }}`$.

Let us now turn to the description of the second regime $`\mathrm{\Phi }(t)>\mathrm{\Phi }^{}`$, relevant to obtain the failure rate up to global failure. Combining (8) and (13), we obtain the expression

$$U\left[\mathrm{\Phi }(t)\right]=T\mathrm{ln}[\gamma (t_ct)](\mathrm{\Phi }^{}<\mathrm{\Phi }<\mathrm{\Phi }_c),$$
(16)
which is analogous to (14). The regime $`\mathrm{\Phi }^{}<\mathrm{\Phi }(t)<\mathrm{\Phi }_c`$ is strongly influenced by thermal fluctuations, so that the disorder term can be neglected to obtain, in view of (16) and (9), $`\mathrm{\Phi }_c\mathrm{\Phi }(t)=f_0\frac{\sqrt{2T\mathrm{ln}\gamma (t_ct)}}{1\sqrt{2T\mathrm{ln}\gamma (t_ct)}}`$. Differentiating both sides of this expression with respect to $`t`$ yields the failure rate

$$\dot{\mathrm{\Phi }}(t)=C(t)/(t_ct),$$
(17)
where $`C(t)=f_0T/[c(1c)^2]`$ with $`c=\sqrt{2T\mathrm{ln}[\gamma (t_ct)]}`$. This is the second important result of our analysis (see also eq.(B11) in cilietal1 ), which shows that, for $`\mathrm{\Phi }>\mathrm{\Phi }^{}`$, the failure rate accelerates towards a finite-time singularity approximately as $`1/(t_ct)`$. Such a behavior has been documented extensively in experiments on rupture of heterogeneous material samsor . Our analysis provides a novel mechanism for the ubiquitous time-to-failure regime observed in heterogeneous material. Strong quenched heterogeneity has been shown to play an essential role in controlling the critical nature of the rupture process Leungetal  and in the existence of a time-to-failure power law such as (17). Here, we confirm that the heterogeneity is essential to renormalize the thermal fluctuations cilietal1 ; Politi . While the philosophy is similar, the mechanism is different and novel. As for the Omori law, the logarithmic corrections in (17) may give an apparent exponent of the power law, slightly smaller than $`1`$, as observed in experiments. Our numerical tests show that expression (17) provides an approximation which coincide almost everywhere with the exact solution inside the interval $`t^{}<t<t_c`$.

There is a simple physical interpretation of the transition between the two above mentioned rate behaviors (15) and (17). To explain the first (rate decaying) regime, consider the degenerate case of spontaneous fracture ($`f_0=0`$, $`\mathrm{\Phi }^{}=1`$) for which $`U(\mathrm{\Phi })=f_s^2/2`$. As time increases, $`f_s`$ grows, the remaining fibers are stronger and the failure rate decays together with the rate of change of the energy barrier. The second regime can be qualitatively understood by taking the limit of zero disorder ($`T_d=0`$, $`\mathrm{\Phi }^{}=0`$), leading to $`U(\mathrm{\Phi })=(1f_a)^2/2`$. The force $`f_a`$ per remaining fiber grows with time, the fibers break more and more easily and the failure rate grows to give the fracture in finite-time. In the intermediate case $`0<\mathrm{\Phi }^{}<\mathrm{\Phi }_c`$, due to the competition between the growth of $`f_s`$ and $`f_a`$, the two regimes co-exist. At early times, the growth of $`f_s`$ dominates giving the Omori and Andrade laws, followed by the growth of $`f_a`$ in the second regime $`\mathrm{\Phi }>\mathrm{\Phi }^{}`$ giving the power law finite-time singularity.

Lastly, we turn to the behavior for $`\mathrm{\Phi }\mathrm{\Phi }^{}`$, which turns out to provide the dominant contribution for the total time for rupture, as shown in cilietal1 ; Politi . Indeed, the fiber bundle spends most of its time in the vicinity of the stationary point $`\mathrm{\Phi }^{}`$, corresponding to the minimum failure rate. In this case, $`U`$ can be expanded as

$$U(\mathrm{\Phi })=U(\mathrm{\Phi }^{})B(\mathrm{\Phi }^{})(\mathrm{\Phi }\mathrm{\Phi }^{})^2,B(\mathrm{\Phi })=\frac{1}{2}\frac{d^2U(\mathrm{\Phi })}{d\mathrm{\Phi }^2},$$

and equation (8) becomes

$$\dot{\mathrm{\Phi }}R(\mathrm{\Phi }^{})\mathrm{exp}\left[\frac{B(\mathrm{\Phi }^{})}{T}\left(\mathrm{\Phi }\mathrm{\Phi }^{}\right)^2\right].$$
(18)
The solution of this equation is

$$\text{erfi }\left(\sqrt{\frac{B(\mathrm{\Phi }^{})}{T}}(\mathrm{\Phi }\mathrm{\Phi }^{})\right)=2\sqrt{\frac{B(\mathrm{\Phi }^{})}{\pi T}}R(\mathrm{\Phi }^{})(tt^{}),$$
(19)
where $`\text{erfi }(z)=\frac{1}{i}\text{erf }(iz)`$ is the imaginary error function. Using its asymptotics $`\text{erfi }(z)\frac{1}{\sqrt{\pi }z}e^{z^2}`$ for large $`z`$ together with (19), equation (18) becomes

$$\frac{d\mathrm{\Psi }^2}{dt}\frac{T}{B(\mathrm{\Phi }^{})(tt^{})},\mathrm{\Psi }=\mathrm{\Phi }\mathrm{\Phi }^{},$$

whose solution yields the fracture rate

$$\dot{\mathrm{\Phi }}\frac{\sqrt{T}\text{sign }(tt^{})}{2|tt^{}|\sqrt{B(\mathrm{\Phi }^{})\mathrm{ln}(\gamma |tt^{}|)}}(\gamma |tt^{}|1).$$

Expression (18) allows us additionally to calculate the total lifetime $`t_c`$ of the fiber bundle:

$$\gamma t_c=\frac{1}{R\left(\mathrm{\Phi }^{}\right)}_{\mathrm{}}^{\mathrm{}}\mathrm{exp}\left[\frac{B(\mathrm{\Phi }^{})}{T}\left(\mathrm{\Phi }\mathrm{\Phi }^{}\right)^2\right]𝑑\mathrm{\Phi }.$$
(20)
The calculation of this integral with the use of (8) gives

$$\gamma t_c4\pi \sqrt{\pi }\frac{D(\mathrm{\Phi }^{})U(\mathrm{\Phi }^{})}{\sqrt{B(\mathrm{\Phi }^{})}}\frac{\sqrt{T_d}}{T}\mathrm{exp}\left[\frac{U(\mathrm{\Phi }^{})}{T}\right].$$
(21)
This expression, together with (11), recovers the main result of cilietal1 ; Politi , while improving on the prefactors to the main Arrhenius-type dependence.

Using (9) and (10), $`U(\mathrm{\Phi }^{})`$ can be written explicitly

$$U(\mathrm{\Phi }^{})=\left[\frac{\mathrm{\Phi }_c\mathrm{\Phi }^{}}{1\mathrm{\Phi }^{}}\pm \sqrt{\mathrm{ln}\left(\frac{\alpha }{(1\mathrm{\Phi }^{})^2}\right)}\right]^2$$
(22)
where the sign $`+`$ (resp. $``$) corresponds to the case $`\mathrm{\Phi }^{}>1/2`$ (resp. $`\mathrm{\Phi }^{}<1/2`$). As shown in Fig. 1, $`U(\mathrm{\Phi }^{})`$ is a non-monotone function of $`T_d`$. Due to the above mentioned competition between quenched disorder and the growth of the actual force $`f_a`$, $`U(\mathrm{\Phi }^{})`$ decreases as long as $`T_d<T_d^{}`$ and then increases with increasing $`T_d`$ beyond $`T_d^{}`$. The first regime $`T_d<T_d^{}`$ corresponds to the effect discovered in Refs. cilietal1 ; Politi  and mentioned above of the renormalization of thermal fluctuations by quenched disorder, and consequently of decreasing strength by increasing the disorder. Since a larger $`U(\mathrm{\Phi }^{})`$ corresponds to a large lifetime through (21), we uncover the new effect of a strengthening of the fiber system by increasing the disorder beyond a certain threshold. All our formulas have been checked by direct numerical integration with excellent agreements. We expect that extensions of the DFBM to non-mean field power law interactions Kunetal  will not change our results qualitatively but may modify the Omori’s and time-to-failure exponents.

We acknowledge useful exchanges with S. Ciliberto and A. Politi.


text:
# Conduction in Carbon Nanotubes Through Metastable Resonant States

(November 12, 2003)

## Abstract

We report here on electrical measurements on individual multi-walled carbon nanotubes (MWNTs) that show that the presence or movement of impurities or defects in the carbon nanotube can radically change its low temperature transport characteristics. The low temperature conductance can either decrease monotonically with decreasing temperature, or show a sudden increase at very low temperatures, sometimes in the same sample at different times. This unusual behavior of the temperature dependence of the conductance is correlated with large variations in the differential conductance as a function of the dc voltage across the wire. The effect is well described as arising from quantum interference of conduction channels corresponding to direct transmission through the nanotube and resonant transmission through a discrete electron state, the so-called Fano resonance.

There is tremendous interest in the transport properties of carbon nanotubes due to their potential for use in future nanodevices nextgeneration , and from their role as canonical models of one-dimensional electron transport luttinger . While individual single-walled carbon nanotubes are expected to show either semiconducting or metallic behavior depending on their chirality physicalproperties , the presence of impurities, defects and interactions is expected to modify this behavior backscattering ; defects . However, measuring the intrinsic transport properties of single or multiwalled carbon nanotubes is complicated by the experimental problem of making electrical contact to the nanotube. Quite often, the contact resistance between the metallic electrodes and nanotube is much higher than the resistance of the nanotube itself, so that the transport properties of the device are determined in large part by the properties of the metal-nanotube contacts. In spite of this problem, transport measurements on high contact-resistance carbon nanotube devices by a number of groups have elucidated the wide variety of physical problems that can be studied in these systems, including the Coulomb blockade set , the Kondo effect kondo , and Luttinger liquids luttinger .

In devices with lower contact resistances, the intrinsic properties of carbon nanotubes can be directly measured. For example, by using different techniques to obtain low resistance contacts, observation of quantization of the conductance in ballistic carbon nanotubes has been reported ballistic , and there are also suggestions of superconductivity in nanotube rope devices made with low resistance contacts superconductivity . However, much interesting physics still remains to be explored. In this Letter, we report on transport measurements on multiwalled carbon nanotubes (MWNTs) with low resistance contacts. The low-temperature differential conductance $`G(V_{dc})`$ of these devices as a function of the dc voltage bias $`V_{dc}`$ across them is highly asymmetric, showing large, reproducible fluctuations that can be as large as 10% of the total conductance. This conductance ‘fingerprint’ can change on thermal cycling to just 2 K, indicating that it might be associated with metastable impurities or mechanical instabilities. The fluctuations in $`G(V_{dc})`$ are reflected in the zero-bias ($`V_{dc}`$=0) conductance as a function of temperature $`G(T)`$, which can either decrease or increase with decreasing temperatures at low temperatures, sometimes showing both dependences in a single sample upon thermal cycling. The sharp structure in $`G(V_{dc})`$ can be described as arising from Fano resonances through resonant states in the device. The position of these peaks and dips frequently change on thermal cycling, indicating that the resonant states may arise from metastable impurities or defects in the device.

Our devices consist of an isolated MWCNT with four Au/Ti electrodes (Fig. 1). The arc-grown MWNTs, typically 2-5 $`\mu `$m long and 25-50 nm wide, were spun onto an oxidized Si substrate from a suspension in dimethylformamide. After locating the nanotubes with respect to fiducial marks, and patterning electrodes by electron-beam lithography, thin film electrodes (5 nm Ti/50 nm Au) were deposited to make contact to the nanotube. Prior to deposition, a short-time oxygen plasma etch was used to improve the metal-nanotube contact resistance. Without this cleaning, the contact resistances were typically in the range of a few k$`\mathrm{\Omega }`$ to a few M$`\mathrm{\Omega }`$, while this process reduced the contact resistance to a few hundred ohms. The presence of four electrodes enables us to make four-terminal resistance measurements on the MWCNTs, eliminating effects of the contact resistance between the MWCNT and the electrodes. In some devices, two external contacts were made to each electrode, as shown in Fig. 1, enabling direct measurements of the electrode-nanotube contact resistance; in the other devices, the contact resistance was inferred from differences between four-terminal and two-terminal measurements. The typical distance between the voltage probes was 300 nm. The samples were measured in a 3He refrigerator and a dilution refrigerator. A home-made ac resistance bridge was used to measure the four-terminal differential resistance $`dV/dI`$ as a function dc current $`I_{dc}`$ using the probe configuration shown in Fig. 1, with ac excitation in the range of 0.25-1 nA to avoid heating effects. In this paper, we shall plot our measurements in terms of the differential conductance $`G`$=1/($`dV/dI`$) as a function of the dc voltage $`V_{dc}`$, obtained by numerically integrating the measured $`dV/dI`$ vs $`I_{dc}`$ curves.

Figure 2 shows $`G(V_{dc})`$ for four different samples with contact resistances ranging from 1 k$`\mathrm{\Omega }`$ to less than 100 $`\mathrm{\Omega }`$, at two different temperatures each, and demonstrates the sample-specific behavior that can be seen in our devices. $`G(V_{dc})`$ is highly asymmetric and shows large fluctuations, with the pattern of the fluctuations being different for different devices. In general, the structure in $`G(V_{dc})`$ becomes sharper as the temperature is lowered: the peaks increase in conductance, while the valleys decrease in conductance. This gives rise to characteristically different temperature dependences for nominally identical samples: for example, if a peak is observed at $`V_{dc}`$=0, $`G(T)`$ will increase with decreasing temperature, while if a valley is observed, it will decrease. (Of course, it is not necessary that either a peak or a valley occur at $`V_{dc}=0`$). This is demonstrated in the insets to the panels in Fig. 2, which show the corresponding $`G(T)`$ for each device. As demonstrated in Fig. 2, the structure in $`G(V_{dc})`$ appears to increase on average as the resistance of the contacts decreases. Indeed, large fluctuations in $`G(V_{dc})`$ are not seen in our samples with contact resistances in the k$`\mathrm{\Omega }`$ to M$`\mathrm{\Omega }`$ range; one sees instead a largely symmetric curve with reduced conductance near $`V_{dc}`$=0, and correspondingly, a monotonically decreasing $`G(T)`$.

$`G(V_{dc})`$ varies not only from sample to sample, but can also change in a single sample as a function of thermal cycling. Figures 3(a) and 3(b) show $`G(V_{dc})`$ for Sample 3 (see the caption of Fig. 2), for two different cooldowns of the dilution refrigerator. As can be seen, there is a radical change after warming to room temperature: while the most prominent feature at low temperature in Fig. 3(a) is a peak, this feature has changed into a dip in Fig. 3(b). As we noted before, this change is also reflected in $`G(T)`$, as shown in the inset to the panels in Fig. 3. Indeed, we found that it was not necessary to warm the samples all the way to room temperature; warming to temperatures on the order of a few Kelvin changed the behavior of the sample, although the change was not as dramatic as shown in Figs. 3(a) and 3(b).

Sample-specific fluctuations of the conductance are well known in mesoscopic systems. In the case of disordered metals and semiconductors, they are associated with quantum interference of electron waves that are scattered by impurities, defects or grain boundaries in the sample ucf . Each scattering event introduces a finite but time-independent shift in the phase of the electron. The phases and phase shifts of the electrons can be modulated by external parameters such as a gate voltage or magnetic field. So long as time-dependent scattering processes that destroy the phase of the electron (processes such as electron-electron or electron-phonon scattering) are negligible, the interference of electrons manifests itself as aperiodic or periodic fluctuations of the conductance as a function of the external parameter ucfexp . Furthermore, since scattering from each impurity or defect introduces a phase shift, the movement of the impurity by even a very small distance (equivalent to 1/$`k_F`$, where $`k_F`$ is the Fermi wave vector) will change the phase shift, resulting in a corresponding change in the conductance pattern noise . Time-dependent conductance changes have been observed in disordered metals, and this mechanism has been shown to be a source of $`1/f`$ noise in metals at low temperatures noiseexp .

A similar interference mechanism can also exist in relatively clean carbon nanotubes with one or a few number of impurities or defects. In this case, the interference can be between electron waves that are directly transmitted, and those that are transmitted via a resonant state. This interference between directly transmitted channels and resonant channels gives rise to a Fano resonance, well known in atomic scattering Fano . For simplicity, we consider the case of a single directly transmitted channel and a single resonant state. Transmission through the resonant channel is described by an amplitude $`t_r(ϵ)=z_r\mathrm{\Gamma }/(2(ϵϵ_0)+i\mathrm{\Gamma })`$ Clerk . Without the factor $`z_r`$, this gives the usual expression for resonant transmission through a localized state of energy $`ϵ_0`$ (measured with respect to the Fermi energy $`ϵ_F`$) and intrinsic energy width $`\mathrm{\Gamma }`$, $`T_r=|t_r|^2=\mathrm{\Gamma }^2/(4(ϵϵ_0)^2+\mathrm{\Gamma }^2)`$, with a transmission of $`T_r=1`$ on resonance ($`ϵ=ϵ_0`$). The transmission amplitude of the direct channel does not depend on $`ϵ`$, and can be described by an expression of the form $`t_d=\sqrt{T_d}e^{i\alpha _d}`$, where $`T_d`$ is the transmission of the direct channel, and $`\alpha _d`$ describes the phase difference between the direct and resonant path. The interference between the two paths is taken into account by taking the sum of the transmission amplitudes to calculate the total transmission coefficient, $`T_t=|t_r+t_d|^2`$. The resulting conductance can be expressed in the Fano form

$`G(ϵ)={\displaystyle \frac{2e^2}{h}}T_t={\displaystyle \frac{2e^2}{h}}T_d{\displaystyle \frac{|2(ϵϵ_0)+q\mathrm{\Gamma }|^2}{4(ϵϵ_0)^2+\mathrm{\Gamma }^2}}`$ (1)
where $`q=i+z_re^{i\alpha _d}/\sqrt{T_d}`$ is the complex Fano parameter Clerk . Note that far from resonance ($`ϵ\mathrm{\Gamma }`$), the conductance reduces to $`G=(2e^2/h)T_d`$.

A number of characteristics distinguish the resulting Fano resonance from other resonances that might occur in carbon nanotubes. First, depending on the phase difference between the resonant and non-resonant transmission channels, the Fano resonance can give rise to a peak or dip in the conductance (or something in between). Second, the Fano lineshape can be asymmetric about the transmission maximum or minimum. Finally, the position of the resonance is determined by the energy of the resonant state, and does not necessarily occur at the Fermi energy (zero bias). These distinguishing characteristics, which can be clearly seen in our data, rule out other possible mechanisms (such as the Kondo effect kondo ) for the structure we observe in the differential resistance of our devices.

Two groups have recently reported observing Fano resonances in carbon nanotube devices. Kim et al. Kim  measured the conductance of crossed MWNTs, and observed a Fano resonance in two of the devices. Noting that a Fano resonance was never observed in devices without crossed MWNTs, they associated the presence of the Fano resonance with the MWNT cross, although the mechanism by which a discrete electron level is created was not discussed. Yi et al. Yi  measured the conductance of crosses consisting of metal electrodes patterned across MWNT bundles. In two of these devices, they observed non-monotonic behavior of the conductance near zero voltage bias, which they ascribed to a Fano resonance arising from interference between a Kondo resonance and non-resonant channels. However, the nature of the localized state giving rise to the Kondo resonance was not made clear. Furthermore, this interpretation is suspect in our opinion, because a Kondo resonance typically arises at the Fermi energy, and the resonances observed by Yi et al. were typically observed at voltage biases of 0.3-0.6 mV. In both papers, the metastable behavior we observe was not reported.

Fig. 4(a) shows the measured Fano resonances for the lower temperatures shown in Figs. 3(a) and 3(b), along with fits to the Fano function given above. To obtain these curves, the background conductance outside the region $`V_{dc}`$ = \[-0.1 mV,-0.02 mV\] for curve 1, and $`V_{dc}`$ = \[-0.05 mV, 0.02 mV\] for curve 2 has been fit, and subtracted from the experimental data. After subtraction, the conductance far from resonance would be 0. However, we arbitrarily introduce an offset of $`4e^2/h`$, corresponding to the conductance of two channels of a single-walled nanotube that coherently interfere with the resonant state, with the assumption that other channels that contribute to the background do not interfere with the resonant state. It can be seen that data are well described by the Fano equation. For curve 1, $`t_d`$ is in phase with $`t_r`$ on resonance, resulting in a maximum of conductance, while for curve 2, $`t_d`$ is out of phase with $`t_r`$ on resonance, resulting in a minimum of conductance. This change in the phase of the non-resonant channel and the resonant channel is a result of the annealing process, which presumably causes a change in the position of the impurity or defect that gives rise to the resonant state.

Fig. 4(b) shows the fitted value of $`|q|`$ and $`\mathrm{\Gamma }`$ as a function of temperature for the conductance dip shown in Fig. 4(a). Both $`|q|`$ and $`\mathrm{\Gamma }`$ increase with increasing temperature. Although the values we obtain for $`\mathrm{\Gamma }`$ are comparable to those obtained by Kim et al., the temperature dependence is different. Kim et al. observed a linear temperature dependence, which they ascribed to thermal broadening of the linewidth, even though $`\mathrm{\Gamma }`$ was less than $`k_BT`$. Our temperature dependence is not linear, which is not surprising, since we are in the regime $`\mathrm{\Gamma }k_BT`$. Fitting to a power law gives a dependence $`\mathrm{\Gamma }T^{1/5}`$, shown as the solid line in Fig. 4(b); the origin of this power law is not clear to us. In contrast, while $`|q|`$ also increases with temperature at higher temperatures, it appears to saturate below 100 mK. The phase of $`|q|`$ remains essentially constant over the fitted temperature range.

What is the possible origin of the resonant states in our samples? Although we are not certain, the metastability we observe suggests impurities or structural defects cause variations in the potential seen by the electrons in the nanotube, which in turn lead to localized electronic states with well-defined resonant energies. Movement of the impurities or defects can give rise changes in the conductance through interference effects. Although the changes we observed in our samples typically occurred on the scale of days at low temperatures, we have also observed changes on the scale of hours. If there are many such impurities moving on a sufficiently rapid time scale, this may be one mechanism for the large $`1/f`$ noise in carbon nanotubes reported by some groups tubenoise .

###### Acknowledgements.

We thank the group of R. P. H. Chang for providing us the nanotubes used in these experiments, and D. E. Prober for a critical reading of the manuscript. Funding for this work was provided by a NASA/MSFC Phase II SBIR, Contract No. NAS8-02102, through a subcontract from Lytec, LLC.


text:
# A Discrete-Event Analytic Technique for Surface Growth Problems

## I INTRODUCTION

Broad interest in surface growth problems and the interface motion stems from many applications that these problems have in practically all sciences. Some well-known examples are crystallization problems, epitaxial growth, wetting phenomena, electrophoretic deposition of polymer chains, and growth of bacterial colonies. The mainstream studies of physical mechanisms that lead to growth and roughening, and interface properties, focus on the so-called universal properties, i.e., large-scale properties of growing surfaces as determined by universal scaling exponents. A central point here is a coarse-grained description of the interface, which leads to an effective growth equation, typical for a universality class. It is safe to say that universal properties of time-evolving surfaces are well understood BS95 . There are many situations, however, where non-universal properties, i.e., those pertaining to the microscopic structure of the interface, are of importance in understanding the physical system. One case is the density of local minima of a virtual-time interface, evolving in parallel discrete-event simulations KTNR00 ; KNGTR03 ; KNR03 . Another case is the movement of random walkers on a growing surface Ch02 . Another case is the density of maxima after solidification and processing of a surface, where this density is related to the frictional force SDG01 . Even for a one-dimensional interface, as in the growth of a two-dimensional crystal SJB00 , a general description of the microscopic structure is challenging, both for analytics and for simulations, due to a variety of length and time scales. However, for a strongly non-equilibrium interface (e.g., in a cold system) a net picture of the interface gains and losses can be simplified. This leads to a variety of simulation models such as, e.g., ballistic deposition, Eden or solid-on-solid models.

Here we introduce a discrete-event analytic technique that provides a means for calculating some non-universal properties of interfaces. The key idea is the construction of a probability distribution for events that take place on the surface. This new approach gives a mean-field like approximation of averages that are determined in simulations. In Sec. II, we explain the main concept by examples of applications to selected non-equilibrium one-dimensional model interfaces. The method is summarized in Sec. III.

## II THE METHOD

We focus on a model in which a one dimensional (1D) surface grows by the deposition of local height increments $`\eta _k`$ at $`L`$ lattice sites, with periodic boundary condition $`h_{L+1}=h_1`$ ($`k`$ enumerates the lattice sites, $`h_k`$ is the local height). The deposited $`\eta _k`$ is a real positive number that can take on continuous values in the interval $`[0;\mathrm{})`$. Our discussion is not specific to the probability distribution from which $`\eta _k`$ is sampled. We consider a rough surface at saturation (Fig. 1). Such a surface can be generated, e.g., in the single-step solid-on-solid models xxx  or in random deposition with surface relaxation Fam86  or by depositing random $`\eta _k`$ at local surface minimum KTNR00 ; KNR03 .

In the latter case, the growth is simulated by the nearest-neighbor rule: $`h_k(t+1)=h_k(t)+\eta _k(t)`$, if $`h_k(t)\mathrm{min}\{h_{k1}(t),h_{k+1}(t)\}`$; otherwise, $`h_k(t+1)=h_k(t)`$ ($`t`$ is the time index). The above deposition rule produces surfaces characterized by the universal roughness exponent KTNR00  $`\alpha =1/2`$. Our objective is to show that it is possible to construct a theoretical probability distribution of deposition events on the surface generated by the above rule. This distribution can then be used to compute measurable quantities at saturation, e.g., the mean density of local minima or maxima. Also, it can serve as the departure point for the construction of other distributions that describe more complex growth processes. Our derivation makes the following two simplifying assumptions. First, we neglect correlations between nearest-neighbor local slopes, which depend on the type of deposition, i.e., the distribution from which $`\eta _k`$ is sampled. Second, we neglect temporal correlations among the groups of the surface configuration classes, as explained later. Because of the above simplifications, our theoretical results for the averages are a mean-field like approximation to the averages measured in simulations.

The adopted growth rule provides the principle for the classification of local lattice-site configurations of the interface. There are only four groups of these configurations (Fig. 2). Each group corresponds to one of the four mutually exclusive discrete events A, B, C and D that take place at site $`k`$ during the deposition attempt at $`t`$: “A” denotes an event when the deposition rule is satisfied from both sides, i.e., $`h_k(t)h_{k1}(t)`$ and $`h_k(t)h_{k+1}(t)`$; “B” denotes an event when the rule is not satisfied from the right; “C” denotes an event when the rule is not satisfied from the left; “D” denotes an event when the rule is not satisfied from either side. The above events can be mapped in unique way onto elementary local site configurations A, B, C and D, presented in Fig. 2. For a closed chain of sites, the set of events when all sites are simultaneously in the A-configuration is of measure zero, and not all $`L`$ sites can have the same elementary configuration. Therefore, in the set of $`L`$ sites there must be at least one site with configuration A. We assign the index $`k=1`$ to one of the sites in the A-configuration and enumerate the other sites accordingly, progressing to the right. Its right neighbor (site $`k=2`$) can be only either in configuration C or D. Similarly, its left neighbor (site $`k=L`$) can be only either in B or D. If site $`k=2`$ is in C then site $`k=3`$ can be only either in C or D. If site $`k=2`$ is in D then site $`k=3`$ can be only either in B or A. If site $`k=3`$ is in B then its right neighbor can be either in B or A, etc. Starting from the site $`k=1`$ and progressing to the right towards $`k=L`$, applying these elementary neighbor rules, we can construct all possible configuration equivalency classes of the entire surface generated by the adopted growth rule. These can be categorized into groups (called $`p`$-groups), based on the number $`p`$ of the deposition events on the surface at $`t`$, i.e., here, the number of local minima in the surface configuration (coded by “A”) at $`t`$. The probability distribution $`f(p;L)`$ of the deposition events is obtained as the quotient of the multiplicity $`M(p)`$ of the $`p`$-group and the total number $`M`$ of configuration classes. (We omit $`t`$ in the notation since the analysis concerns steady-state evolution.) The mean density $`u(L)`$ of local minima is given by the generally valid formula:

$$u(L)=\underset{p}{}u(p)f(p;L),$$
(1)
where the summation is over all $`p`$-groups of the admissible surface-configuration classes, $`u(p)`$ is the density characteristic for each group, and $`f(p;L)`$ is the frequency of the occurrence of the $`p`$-group at saturation, i.e., the probability of generating a surface with $`p`$ local minima.

For example, the binary tree for the construction of possible surface-configuration classes for $`L=5`$ is shown in Fig. 3. Looking along its branches, starting from the leading A at the fixed $`k=1`$ position, one identifies $`M=8`$ configuration classes of the entire surface: (1) ACCCD; (2) ACCDB; (3) ACDBB; (4) ACDAD; (5) ADBBB; (6) ADBAD; (7) ADACD; and (8) ADADB. Each surface configuration represents a class of infinitely many topologically equivalent deformations (Fig. 4) since the deposited random height increment is a real positive number. There are only two $`p`$-groups. In the first, there are four classes with one letter A: $`M(1)=4`$, $`f(1;5)=1/2`$, and $`u(1)=1/5`$. In the second, there are four classes with two letters A: $`M(2)=4`$, $`f(2;5)=1/2`$, and $`u(2)=2/5`$. By Eq. (1), $`u(5)=3/10`$.

To find $`f(p;L)`$, one can exploit the recurrent structure of the binary tree for general $`L`$ in counting classes that contain the A-configuration at exactly $`p`$ sites KNR03 . Counting gives: $`M=2^{L2}`$, $`M(p)=(L1)!/((2p1)!(L2p)!)`$, and $`p=1,2,\mathrm{},[L/2]`$ ($`[L/2]=L/2`$ for even $`L`$; $`[L/2]=(L1)/2`$ for odd $`L`$). Thus,

$$f(p;L)=\frac{1}{2^{L2}}\left(\genfrac{}{}{0pt}{}{L1}{2p1}\right).$$
(2)
In deriving $`f(p;L)`$ the underlying assumption is that at saturation any class of the entire surface configurations is generated with probability $`1/M`$. The theoretical $`f(p;50)`$ is compared with the simulated distributions in Fig. 5. We considered three deposition types: the Poisson type, where $`\eta _k`$ was sampled from the Poisson distribution; Gaussian type, where $`\eta _k`$ was sampled from the Gaussian; and, the uniform deposition with $`\eta _k`$ being a uniform random deviate. In all cases, Eq. (2) is a very close approximation to time-averaged distributions measured in simulations.

Having Eq. (2), all moments of $`f(p;L)`$ can be computed exactly KNR03 , e.g., for $`L4`$ its variance is

$$\sigma ^2(L)=\underset{m}{}\left(mp(L)\right)^2f(m;L)=\frac{L1}{16}.$$
(3)
For $`L3`$ the mean number of local minima is

$$p(L)=\underset{m}{}mf(m;L)=\frac{L+1}{4}.$$
(4)
The mean density of local minima, presented in Fig. 6, is simply obtained from Eq. (4) as $`u(L)=p(L)/L`$.

Suppose, the adopted growth rule is slightly modified so the deposition at a local minimum (site A) happens with probability $`q_A`$. Then, the new probability distribution $`F(n;L,q_A)`$ is the product of the probability $`f(n_A;L)`$ of drawing a surface that has exactly $`n_A`$ local minima, given by Eq. (2), and the probability $`B(n;n_A,q_A)`$ that $`n`$ depositions occur on this surface. The latter is given by binomial distribution

$$B(n;n_A,q_A)=\left(\genfrac{}{}{0pt}{}{n_A}{n}\right)q_A^n\overline{q}_A^{n_An},$$
(5)
where first factor gives the number of ways in which $`n`$ depositions can be distributed among $`n_A`$ sites, and $`\overline{q}_A=1q_A`$ is the complementary probability of the no-deposition event at site A. In the old rule: $`B(n;n_A,q_A=1)=\delta _{n,n_A}`$ so, trivially, $`F(n;L;q_A)=f(n_A;L)`$.

Consider next a more complicated rule, where at the $`t`$-th deposition attempt, the deposition of $`\eta _k(t)`$ is not exclusively confined to local minimum but may also happen at other sites B, C or D with the corresponding probabilities $`q_B`$, $`q_C`$ and $`q_D`$, respectively. We assume now $`q_A=1`$. The probability $`F(n;L)F(n;L;n_B,n_C,n_D;q_B,q_C,q_D)`$ of having $`n=n_A+n_B+n_C+n_D`$ deposition sites in the interface is the product:

$`F(n;L)`$ $`=`$ $`f(n_A;L)B(n_B;N_B,q_B)`$ (6)
$`\times `$ $`B(n_C;N_C,q_C)B(n_D;N_D,q_D).`$
In Eq. (6), $`f(n_A;L)`$ is the probability of generating a surface with $`n_A`$ local minima, each of which is the deposition site; the binomial $`B(n_B;N_B,q_B)`$ is the probability of having $`n_B`$ depositions at $`N_B`$ sites in the local B-configuration on the generated surface; similarly, $`B(n_C;N_C,q_C)`$ and $`B(n_D;N_D,q_D)`$ are probabilities that $`n_C`$ and $`n_D`$ depositions occur at $`N_C`$ and $`N_D`$ sites in the local C- and D-configurations, respectively. For any surface, $`N_A+N_B+N_C+N_D=L`$. For surfaces growing on a closed ring of $`L`$ sites, the number $`N_D`$ of local maxima matches the number $`N_A`$ of local minima. This gives $`N_B+N_C=L2n_A`$. By Eq. (1), the mean number $`n(L)`$ of deposition sites in the interface is

$$n(L)=\underset{n_A=1}{\overset{\left[\frac{L}{2}\right]}{}}\underset{n_B=0}{\overset{N_B}{}}\underset{n_C=0}{\overset{N_C}{}}\underset{n_D=0}{\overset{n_A}{}}nF(n;L).$$
(7)
The mean density of these sites is $`U(L)=n(L)/L`$.

In Eq. (7), the summation over $`n_D`$ is easily performed: $`_{n_D=0}^{n_A}(n_A+n_B+n_C+n_D)B(n_D;N_D,q_D)=n_A(1+q_D)+n_B+n_C`$. This gives

$$n(L)=\underset{n_A=1}{\overset{\left[\frac{L}{2}\right]}{}}f(n_A;L)\left(n_A\left(1+q_D\right)+S_{BC}\right),$$
(8)
where $`S_{BC}`$ is the mean number of deposition events in the set of sites other than A or D. The computation of $`S_{BC}`$, though elementary, is a bit tricky since one must distinguish between depositions in the set of B-sites and depositions in the set of C-sites to avoid double counting:

$`S_{BC}`$ $`=`$ $`{\displaystyle \underset{n_B=0}{\overset{N_B}{}}}B(n_B;N_B,q_B){\displaystyle \underset{n_C=0}{\overset{N_C}{}}}B(n_C;N_C,q_C)`$
$`\times `$ $`(n_B+n_C)=(L2n_A)q_Bq_C.`$
Substitution to Eq. (8) gives for $`L3`$:

$$n(L)=\frac{1+q_D+2\overline{Q}}{4}L+\frac{1+q_D2\overline{Q}}{4},$$
(9)
where $`\overline{Q}=1Q`$, $`Q=\overline{q}_B+\overline{q}_C\overline{q}_B\overline{q}_C`$.

One application of Eq. (9) is the case when $`\overline{q}_B=\overline{q}_C=q/2`$ and $`\overline{q}_D=q`$, where the probability $`q=\sqrt{2/N}`$ is parametrized by an integer $`N2`$. Here, the mean density of lattice sites where the deposition takes place is

$$U(L)=\left(1\frac{q}{2}\right)\left(1\frac{q}{4}\frac{L1}{L}\right).$$
(10)
Predictions of Eq. (10), compared to simulations with the Poisson deposition in Fig. 7, are consistent with the fact that the approximation $`f(n_A;L)`$ gives slight overestimates of the density of local minima at the high-$`L`$ end (Fig. 6). Nevertheless, the analytical $`U(L)`$ reflects well the overall tendencies in the data. As $`q0`$, $`U(L)`$ approaches $`1`$. In this limit-case each attempted deposition event at each lattice site becomes certain. Such a case corresponds to the random deposition growth rule. The corresponding interface is entirely uncorrelated and never saturates. This is the asymptotic limit to which the curves in Fig. 7 converge when $`N\mathrm{}`$.

## III OUTLOOK AND CONCLUSION

In summary, we have introduced a new approach for the study of non-universal properties of surfaces. We showed the explicit connection between the interface morphology and the event statistics on the interface at saturation. The key concept here is the ability to build distinct equivalency classes of the entire surface configurations from its local-site configurations. In this construction, the growth rule defines the set of events that can be mapped in a unique way on the set of elementary local-site configurations.

Small differences between simulation results and the approximate analytical results come mainly from neglecting temporal correlations among $`p`$-groups of surface-configuration classes. These correlations are intrinsically present in the computation of averages over time series in simulations but they are absent in our derivation. They depend on the probability distribution from which the deposited random height increments are sampled. The proper handling of temporal correlations is essential in the microscopic description of the growth phase (before saturation), which is a challenging problem still open. Some insight to this issue may be found in recent works of Rikvold and Kolesik RK , who used the concept of equivalency classes in connecting the interface microstructure with its mobility for Ising and solid-on-solid models with various stochastic dynamics.

Although the focus of the study presented here is on obtaining analytical results in the mean-field spirit, it is interesting to notice the sensitivity of simulation results to the probability distribution used to select random height increments. Negligible variations between the derived analytical probability distribution and simulated distributions with Gaussian, Poissonian and uniform increments demonstrate that the theoretical distribution presents well the statistics of events on the surface. Temporal correlations between interfaces generated via random deposition at local minima deserve further analytical studies.

In conclusion, we observe that our approach could be applied to a variety of other growth rules and models. The current applications to 1D problems are illustrative examples. The main advantage of the approach is that it enables one to compute analytically quantities that can be only estimated qualitatively by continuum methods.

###### Acknowledgements.

This work is supported by NSF grants DMR-0113049 and DMR-01200310, and by the ERC Center for Computational Sciences at MSU. It used resources of the National Energy Research Scientific Computing Center, supported by the Office of Science of the US Department of Energy under contract No. DE-AC03-76SF00098.