text:
# Instability and Spatiotemporal Dynamics of Alternans in Paced Cardiac Tissue

\[

## Abstract

We derive an equation that governs the spatiotemporal dynamics of small amplitude alternans in paced cardiac tissue. We show that a pattern-forming linear instability leads to the spontaneous formation of stationary or traveling waves whose nodes divide the tissue into regions with opposite phase of oscillation of action potential duration. This instability is important because it creates dynamically an heterogeneous electrical substrate for inducing fibrillation if the tissue size exceeds a fraction of the pattern wavelength. We compute this wavelength analytically as a function of three basic length scales characterizing dispersion and inter-cellular electrical coupling.

\]

It is well established that the duration of cardiac excitation can oscillate from beat to beat at sufficiently short pacing interval . Pioneering studies by Nolasco and Dahlen  and Guevara et al.  have demonstrated that the generic sequence LSLS… of long and short action potential duration ($`APD`$), known as alternans, is a direct consequence of the restitution relationship

$$APD^{n+1}=f(DI^n)$$
(1)
between the $`APD`$ generated by the $`n^{th}+1`$ stimulus, $`APD^{n+1}`$, and the diastolic time interval $`DI^n`$ during which the tissue recovers its resting properties after the end of the previous ($`n^{th}`$) action potential. If we denote the interval between the $`n^{th}`$ and $`n^{th}+1`$ stimulus by $`T^n`$, we must have $`DI^n=T^nAPD^n`$. Then, for a fixed period: $`T^n=\tau `$ for all $`n`$, Eq. 1 yields the map $`APD^n=f(\tau APD^{n1})`$ that has a period doubling instability if the slope $`f^{}`$ of the restitution curve evaluated at its fixed point exceeds unity, which generically occurs as the period is decreased.

Over the last decade, the study of alternans , and their control , has become a main focus of research because of the potentially crucial link of this dynamical instability with cardiac fibrillation . However, there is presently no simple analytical understanding of how the bifurcation to alternans is manifested *spatiotemporally* in paced cardiac tissue. Analytical progress to date is limited to the one-dimensional circulation of electrical impulse in a ring of tissue .

In this letter, we derive an equation that governs the spatiotemporal dynamics of alternans close to the onset of instability. This enables us to obtain a quantitative analytical understanding of the formation of recently observed complex patterns of $`APD`$ oscillations that can promote fibrillation . A crucial feature of these patterns is that the $`APD`$ oscillates with opposite phases in two (or more) spatially extended regions of tissue, i.e. with a sequence LSLS… in one region and SLSL… in the other. These “discordant alternans” have been observed experimentally in both two-dimensional  and linear strands  of cardiac tissue, as well as in ionic model simulations . Moreover, they have been shown to lead to the formation of conduction blocks  as well as to the onset of spiral wave formation and fibrillation . We show here that discordant alternans result from a *pattern-forming* linear instability that has interesting similarities with classic instabilities leading to the spontaneous formation of spatially periodic patterns in nature (such as Rayleigh-Bénard convection, Taylor-Couette flow, etc ), but also presents some unique features.

We consider a one-dimensional (1-d) homogeneous cable of length $`L`$ paced at period $`\tau `$ from one end ($`x=0`$). Close to the onset of instability, we can expand the $`APD`$ and the period in the form

$`APD^n(x)`$ $``$ $`APD_c+a(x,t)e^{i\pi n},`$ (2)
$`T^n(x)`$ $``$ $`\tau _c\delta \tau +b(x,t)e^{i\pi n},`$ (3)
where $`APD_c`$ and $`\tau _c`$ are the $`APD`$ and the period evaluated at the bifurcation point of the map ($`f^{}=1`$), $`x`$ measures the position along the cable, and $`\delta \tau \tau _c\tau \tau _c`$. In this range of period, $`a`$ and $`b`$ vary slowly from beat to beat, which allows us to treat the time, $`tn\tau `$, as a continuous variable; the fast beat-to-beat oscillations are contained in the exponential factor $`e^{i\pi n}`$.

A relation between $`a`$ and $`b`$ can first be derived by noting that $`T^n(x)`$ is the difference of arrival time of two subsequent action potentials at $`x`$ , or

$$T^n(x)=\tau +_0^x\frac{dx^{}}{c(DI^n(x^{}))}_0^x\frac{dx^{}}{c(DI^{n1}(x^{}))},$$
(4)
where the first (second) integral on the right-hand-side is the time required for the leading front of the action potential to travel from the paced end to $`x`$ at the $`n^{th}`$ ($`n^{th}1`$) stimulus; concomitantly, $`c(DI)`$ is the standard dispersion curve that relates the propagation speed of this front with the local diastolic interval. The dispersion curve is typically steeply increasing at small $`DI`$ and flat at large $`DI`$. Substituting Eqs. 2-3 in Eq. 4 with $`DI^n(x)=T^n(x)APD^n(x)`$, and expanding to linear order in $`a`$ and $`b`$, we obtain $`b(x)_0^xa(x^{})𝑑x^{}/\mathrm{\Lambda }`$ where we have defined $`\mathrm{\Lambda }c^2/(2c^{})`$, with $`c`$ and $`c^{}dc/dDI`$ evaluated at the bifurcation, and assumed that $`\mathrm{\Lambda }`$ is much larger than the scale over which $`a`$ varies.

Next, in order to derive an evolution equation for the amplitude $`a(x,t)`$, we first neglect the influence of the electrical coupling between cells on the $`APD`$. This allows us to assume that the restitution relationship (1), and hence the second iteration of the map

$$APD^{n+2}=f(T^{n+1}f(T^nAPD^n)),$$
(5)
continues to be valid even when the $`APD`$ is non-spatially uniform; we shall soon see why it is crucial to relax this assumption. We substitute Eqs. 2-3 in Eq. 5, and expand the right-hand-side keeping only the dominant linear and weakly nonlinear terms. Furthermore, we use the aforementioned fact that, close to onset, $`a`$ varies slowly from beat to beat, and we therefore expand the left-hand-side as $`APD^{n+2}APD^n+2a/ne^{i\pi n}`$, where $`a/n=\tau a/t`$. Finally, we use the integral relation between $`a`$ and $`b`$ derived earlier. After equating both sides of Eq. 5, we obtain

$$\tau _ta=\sigma aga^3_0^x\frac{dx^{}}{\mathrm{\Lambda }}a(x^{}),$$
(6)
where $`\sigma f^{\prime \prime }(\tau \tau _c)/2`$, $`gf^{\prime \prime 2}/4f^{\prime \prime \prime }/6`$, and all derivatives are evaluated at the bifurcation point.

In order to test this evolution equation, we simulate the standard cable equation

$$_tV=D_x^2V\left(I_{\mathrm{ion}}+I_{ext}\right)/C_m,$$
(7)
for the membrane current $`I_{\mathrm{ion}}`$ given by the Noble model  with time in units of millisecond (ms), $`D=2.5\times 10^4`$ cm2/ms, $`C_m=12\mu `$F/cm2, dx=0.01 cm, dt=0.05 ms, and $`I_{ext}`$ modeling a sequence of stimuli applied at $`x=0`$ at the pacing interval $`\tau `$. To determine the parameters of the amplitude equation, we compute the restitution and dispersion curves by pacing Eq. 7 in a short cable and by using two subsequent stimuli spaced by different intervals to vary $`DI`$; we also use $`V=40`$ mV as threshold of the transmembrane voltage to define the $`APD`$. We impose zero gradient boundary conditions on $`V`$ and $`a`$ at the two ends of the cable in all our simulations.

Fig. 1(a) shows that Eq. 6 produces discordant alternans consistent with the picture that restitution and dispersion suffice to produce this state . However, the magnitude of the spatial gradient of $`a`$ at the node increases with time, and can be shown to diverge in a finite time. This yields an unphysical spatial discontinuity of $`APD`$, which also occurs if the system of coupled maps (Eqs. 1 and 4), from which the amplitude equation is derived, is solved numerically. Note that the analogous system of coupled maps for a 1-d pulse circulating in a ring produces a smooth $`APD`$ modulation starting from a smooth initial condition , which highlights the much stronger role of dispersion at producing heterogeneity of $`APD`$ during pacing than circulation.

This discontinuity, which is absent in the simulation of the cable equation in Fig. 1(b), can be cured by adding spatial derivative terms to the amplitude equation. Since the cable is paced at one end, the underlying basic state (i.e. traveling pulses) is not invariant under parity symmetry. Hence, in addition to $`_x^2a`$, a term proportional to $`_xa`$ must generally be included, which yields the final form of our amplitude equation

$$\tau _ta=\sigma aga^3_0^x\frac{dx^{}}{\mathrm{\Lambda }}a(x^{})w_xa+\xi ^2_x^2a.$$
(8)
To calculate the new lengthscales $`w`$ and $`\xi `$, we must determine how the electrical coupling between cells modifies the restitution relationship (Eq. 1). For complex electrophysiological models such as Noble, this generally needs to be done numerically using a procedure that will be discussed elsewhere. For a simple two-variable ionic model, which we study below, the analytical expressions

$`w`$ $`=`$ $`2D/c,`$ (9)
$`\xi `$ $`=`$ $`(D\times APD_c)^{1/2},`$ (10)
can be derived by interpreting Eq. 7 as a diffusion equation with a source $`I_{ion}`$, and expressing $`V`$ as space-time integral of $`GI_{ion}`$, where $`G`$ is the standard Green’s function of the diffusion equation. This integral can then be evaluated analytically because the action potential shape is simply triangular for this model, and used to calculate $`w`$ and $`\xi `$. Eq. 10 has the simple physical interpretation that $`V`$ diffuses a length $`\xi `$ in the time interval of one $`APD`$. Therefore, the repolarization of a given cell is influenced by other cells within a length $`\xi `$ of cable. In addition, repolarization of this cell is influenced unequally by its left and right neighboring cells because these cells are activated at different times by the propagating wavefront. Clearly, this asymmetry must vanish in the limit $`c\mathrm{}`$ where all cells are activated simultaneously consistent with Eq. 9. Fig. 1(b) shows that our regularized amplitude equation (8) now produces a smoothly varying stationary profile of $`a`$ that agrees well with the simulation of the cable-Noble equation, where $`a`$ is obtained from the $`APD`$ using Eq. 2.

The genesis of discordant alternans can be understood by computing the linear stability spectrum of the spatially homogeneous state ($`a=0`$). We have calculated this spectrum both numerically for different $`L`$, and analytically for the large $`L`$ limit. The main result is that the wave pattern can emerge from the amplification of *either* a unique finite wavelength mode, which yields a stationary pattern, *or* from a discrete set of complex modes that approach a continuum in the limit $`L\mathrm{}`$, and yields a traveling pattern. There is indeed experimental evidence for both stationary  and traveling  waves.

We treat here explicitly the large $`L`$ limit since it provides the basis to understand finite-$`L`$ patterns. In this limit, we can analyze stability by differentiating Eq. 8 with respect to $`x`$ and letting $`a(x,t)e^{ikx+\mathrm{\Omega }t}`$, with both $`\mathrm{\Omega }\mathrm{\Omega }_r+i\mathrm{\Omega }_i`$ and $`kk_r+ik_i`$ complex, which yields at once the eigenvalue equation

$$\mathrm{\Omega }=\sigma \xi ^2k^2i\left[wk1/(\mathrm{\Lambda }k)\right].$$
(11)
When dispersion is weak, a unique real mode with $`k_i=0`$ grows faster than the other complex modes and yields a stationary pattern. Its wavelength $`\lambda =2\pi /k_r`$, is determined by the condition $`\mathrm{\Omega }_i=0`$, which yields

$$\lambda =2\pi (w\mathrm{\Lambda })^{1/2},$$
(12)
in good agreement with the wavelength observed in simulations of the cable-Noble equation (Table I). Note that $`\mathrm{cos}k_rx`$ is an exact eigenvector of Eq. 8 linearized around $`a=0`$ that satisfies $`_xa=0`$ at the two cable ends when $`L`$ is an integer multiple of $`\lambda /2`$. The threshold of instability occurs when $`\mathrm{\Omega }_r=0`$, or for a period $`\tau _{th}`$ defined by $`\sigma _{th}=f^{\prime \prime }(\tau _{th}\tau _c)/2=\xi ^2/(w\mathrm{\Lambda })`$.

In the opposite limit where dispersion is strong, complex modes that grow exponentially at large $`x`$ are the most unstable. It is simple to deduce from Eq. 11 that each $`k`$-mode travels towards the pacing end of the cable, but a wave-packet constructed from a linear superposition of these modes has a group velocity that makes the packet move away from the pacing end. This is the signature of a convective instability  where perturbations are transported as they grow, similarly to Taylor-Couette vortices developing in an axial flow . In such a situation, patterns are only transient unless the group velocity vanishes, or $`\mathrm{\Omega }_i/k_r=0`$, and hence they grow at a fixed point in space. Moreover, the fastest growing wavelength that dominates at large time must correspond to a maximum of $`\mathrm{\Omega }_r`$, which yields the additional condition $`\mathrm{\Omega }_r/k_r=0`$. These two conditions are equivalent to requiring that $`d\mathrm{\Omega }/dk=0`$. From this, we deduce that the threshold of absolute instability occurs when $`\sigma _{th}=(3/2)(\xi /2\mathrm{\Lambda })^{2/3}`$, with a pattern of wavelength

$$\lambda =(4\pi /\sqrt{3})(2\xi ^2\mathrm{\Lambda })^{1/3},$$
(13)
which travels with phase velocity $`\mathrm{\Omega }_i\lambda /(2\pi )`$ where $`\mathrm{\Omega }_i=(3\sqrt{3}/2)(\xi /2\mathrm{\Lambda })^{2/3}`$. It can also be deduced that traveling waves are favored over a stationary pattern when $`\mathrm{\Lambda }=2c^2/c^{}\xi ^4/w^3`$, and hence for large enough dispersion ($`c^{}`$), and vice versa in the opposite limit.

Nonlinearities are mainly responsible for saturating the amplitude of the linear modes, but both the wavelength and evolution of the pattern can generally vary with distance from onset. There is a close analogy between the amplitude equation (8) and the real Ginzburg-Landau equation that has been extensively studied in the context of phase transitions and front propagation . The dynamics is richer here because the integral term originating from dispersion causes a non-local interaction of the fronts separating two out-of-phase oscillating regions with the pacing end of the cable.

As a final test of our theory, we study a two-variable model where all the parameters of the amplitude equation can be calculated analytically including $`w`$ and $`\xi `$ given by Eqs. 9-10, and for which our theory predicts traveling waves. In this model, $`I_{\mathrm{ion}}/C_m`$ is the sum of a slow outward current, $`I_o/C_m=\tau _0^1(S+(1S)V/V_c)`$, and a fast inward current, $`I_i/C_m=\tau _a^1hS`$, where inactivation of the latter is controlled by

$$_th=(1Sh)/(\tau ^{}(1S)+S\tau ^+)$$
(14)
In addition, $`V`$ is dimensionless and $`S(1+\mathrm{tanh}((VV_c)/ϵ))/2`$. We choose $`V_c=0.1`$, $`\tau _0=150`$ ms, $`\tau _a=6`$ ms, $`\tau ^{}=60`$ ms, $`\tau ^+=12`$ ms, $`ϵ=0.005`$, and we simulate Eq. 7 with the same $`D`$ as before, dx=0.01 cm, dt=0.02 ms, and $`V=0.1`$ as threshold for the $`APD`$.

As predicted, we observe traveling waves in this model with a wavelength that agrees well with Eq. 13. Fig. 2 illustrates sustained and transient wave patterns above (left) and below (right) the onset of absolute stability, respectively. Furthermore, Fig. 3 shows that the results of simulations of the cable and amplitude equations are in good agreement over a wide range of $`L`$ and $`\tau `$, including the fact that the onset of instability occurs at a shorter period in a larger cable (and, hence, for a slope of the restitution curve larger than unity). Fig. 3 is qualitatively similar for stationary waves, but this overstabilization is smaller because dispersion is weaker.

The $`APD`$ oscillations at a fixed $`x`$ are periodic when the wave is stationary and quasiperiodic when it travels. For both cases, we find here that $`\lambda `$ is independent of $`L`$, in contrast to the oscillations produced by a pulse circulating in a ring where it is known that $`\lambda 2L/(2i+1)`$ for weak dispersion  (with $`i`$ integer and $`L=`$ ring perimeter). As will be discussed elsewhere, our theory applied to the ring shows that the bifurcation to alternans is finite dimensional with $`i=0`$ being the most unstable mode, in agreement with the fact that it is this mode that is generically selected in experiments  and ionic model simulations . In addition, it shows that the gradient term ($`w_xa`$) can lead to quasiperiodicity even in the absence of dispersion ($`c^{}=0`$).

Our results demonstrate that the formation of discordant alternans is crucially affected by the effect of electrical coupling (diffusion) on repolarization, in addition to restitution and dispersion. Dispersion is responsible for the formation of nodes and spatial gradients of $`APD`$ that steepen with time. Diffusion, in turn, tends to spread the $`APD`$ spatially, and also induces a drift of the pattern away from the pacing site that is induced by the more subtle gradient term ($`w_xa`$) in the amplitude equation. When dispersion is sufficiently weak, drift balances dispersion and produce a stationary pattern. In the opposite limit, the tendency for dispersion to form steep gradients of $`APD`$ is balanced by the spreading effect of diffusion. Nodes then travel, cyclically disappearing (appearing) at the pacing (opposite) end of the cable.

In conclusion, we have derived a simple evolution equation that describes the universal spatiotemporal dynamics of small amplitude alternans in paced cardiac tissue. Moreover, we have shown that discordant wave patterns , which are linked to fibrillation , result from a finite wavelength linear instability. Hence, their formation requires a minimum tissue size $`L_{min}\lambda /4`$, required for at least one node to form. The value of $`L_{min}`$ that we measure in simulations of reaction-diffusion models are actually close to $`\lambda /4`$ with $`\lambda `$ predicted by Eq. 12 and Eq. 13, respectively (Table I). This equation can be extended to higher dimensions as well as to include heterogeneities. Our preliminary results in a 2-d homogeneous tissue are very similar to 1-d. Finally, this equation is readily applicable to model a non-constant pacing interval and provides a theoretical basis to study the control of alternans in spatially extended tissue.

This research was supported in part by NIH SCOR in Sudden Cardiac Death P50-HL52319.


format: text

text:
# Crumpling a Thin Sheet

## Abstract

Crumpled sheets have a surprisingly large resistance to further compression. We have studied the crumpling of thin sheets of Mylar under different loading conditions. When placed under a fixed compressive force, the size of a crumpled material decreases logarithmically in time for periods up to three weeks. We also find hysteretic behavior when measuring the compression as a function of applied force. By using a pre-treating protocol, we control this hysteresis and find reproducible scaling behavior for the size of the crumpled material as a function of the applied force.

If you, the reader, were to rip a page from this journal and crumple it, squeezing it with your hands into a ball as hard as you can, the resulting object is still more than $`75\%`$ air. What gives this crumpled sheet its surprising strength and how does the ultimate size of the sheet depend on the forces applied?

The energy stored in crumpled sheets has been investigated in the limit where the thickness, $`\delta `$, of the sheet is much less than the lateral dimension, $`L`$ (e.g., the diameter of a circular sheet) 1 ; 2 ; 3 ; 4 ; 5 . Analysis of these objects has focused on two characteristic structures caused by the crumpling: singular, conical points and the curved ridges, which store most of the energy, connecting them. As the sheet is crumpled farther, the ridges must collapse so as to remain within the confining container; and an increasing number of smaller ridges must form. In the limit where the diameter of the container, $`D`$, is much smaller than the initial size of the sheet but still much larger than the sheet thickness ($`\delta DL`$) a scaling relation is predicted between the energy stored and the size of the resulting crumple: $`ED^\beta `$. This prediction assumes that the forces resisting crumpling are conservative. One knows however, that frictional forces are involved as the sheet rubs against itself and against the constraining walls. How do these forces affect the resulting energy balance? In addition plastic flow can occur in regions of high curvature effectively constraining the geometry of the energy bearing structures.

The notion of scaling properties describing a crushed sheet goes back to the initial experiments of Kantor et al. 5  Further investigations of such sheets were made by Gomes and collaborators 6 ; 7 . More broadly, scaling behavior of force with compression in tenuous structures was reported for random fibrous material by Baudequin et al. 8  and for colloidal aggregates by several groups 9 . Here we report distinctive forms of scaling not seen in those previous studies.

To study how the external dimensions of a crumpled sheet depend on the confining force, we have placed a large circular ($`L=34cm`$) sheet of thin ($`\delta =12.5\mu m`$) aluminized Mylar 10  under a weighted piston in a plastic cylindrical cell of diameter $`D=10.2cm`$. The cylinder is surrounded by a box to protect it from air currents in the room and is mounted on an optical table to minimize the effects of vibration. We measure compression by measuring the height, $`h`$, of the piston above the base. No special care is taken in the initial crumpling that allows the large sheet to be placed within the cylinder.

After a mass, $`M`$, is placed on the piston, the piston settles to a new height held up by the crumpled sheet below it. However, this behavior is not as simple as one might have expected: $`h`$ continues to decrease long after the mass is initially introduced. As shown in Fig. 1(a), over seven decades in time, $`2\times 10^1sec<t<2\times 10^6sec`$, the entire duration of the experiment, the Mylar height decreases logarithmically:

$$h=ab\mathrm{log}(t/sec)$$
(1)
where $`a`$ and $`b`$ are constants. Even after three weeks the piston did not reach its asymptotic height. Similar relaxation occurred in other tenuous materials, as shown in the insets to Fig. 1(a). Absorbent tissue paper 11  showed a $`logt`$ dependence like that of the Mylar. Cotton balls 12  showed a slightly smaller relaxation over several weeks, in which the logarithmic dependence extended only up to $`10^4sec`$. Although this logarithmic behavior is robust and is found each time the Mylar is compressed, the constants $`a`$ and $`b`$ fluctuate between runs with no monotonic evolution of those constants as the sheet is crumpled more times. Fig. 1(b) shows a weak correlation between $`a`$ and $`b`$ for the same $`200g`$ load on the piston.

Albuquerque and Gomes 7  studied stress relaxation of crumpled aluminum foil under fixed strain and found stretched-exponential behavior quite different from the logarithmic dependence of strain that we measured at fixed stress. Aluminum foil is much more maleable than the Mylar we have used and presumably has more plastic flow along the ridges and vertices. This may be the origin of the different relaxation kinetics. Menon 13  has found similar $`\mathrm{log}t`$ behavior as our results in stress relaxation.

Relaxation continues even if the compression is stopped. Fig. 2(a) shows the time dependence of the height of a Mylar sheet after a $`200g`$ mass has been place on it. After $`500sec`$. the piston is lifted very slightly and fixed in place so that it can no longer compress the material. After $`1000sec`$, the piston is released and quickly sinks to a new position comparable to where it would have been had it never been restrained; thereafter it relaxes at nearly the same logarithmic rate as it did initially.

What is responsible for the slow decrease in height over these extended periods? External vibrations affect the compression. Fig. 2(b) shows $`h(t)`$ both before and after the container is vibrated laterally with an acceleration of $`2.5m/sec^2`$ at $`7.5Hz`$. Initially, the container was at rest and the logarithmic slope, $`b`$, was determined. After approximately $`200sec`$, vibrations were begun and $`b`$ nearly doubled. This change, while substantial, is not so large as to indicate that the slope in the quiescent state is caused by the ambient vibrations in the room which are less than $`10^2m/sec^2`$.

This logarithmic decay of $`h`$ suggests activated relaxation (with an effective temperature $`T_{\mathrm{𝑒𝑓𝑓}}`$ introduced to account for external vibrations 14 ) of the sheet between metastable minima separated by energy barriers, $`\mathrm{\Delta }E`$: $`dh/dt=\zeta \mathrm{exp}[\mathrm{\Delta }E/(k_BT_{\mathrm{𝑒𝑓𝑓}})]`$. As $`h`$ decreases, the crumpled sheet becomes progressively stiffer suggesting that the barriers are increasing. To lowest order (i.e., $`\mathrm{\Delta }EE_0Ch`$) we find:

$$dh/dt=\zeta \mathrm{exp}[(E_0Ch)/(k_BT_{\mathrm{𝑒𝑓𝑓}})]$$
(2)
which has a solution of the form of Eq. 1 where $`a=d_1+E_0/C(k_BT_{\mathrm{𝑒𝑓𝑓}}/C)\mathrm{ln}(C\zeta /k_BT_{\mathrm{𝑒𝑓𝑓}})`$ and $`b=k_BT_{\mathrm{𝑒𝑓𝑓}}/C`$ where $`d_1`$ is a constant. As $`T_{\mathrm{𝑒𝑓𝑓}}`$ increases (e.g., due to external vibrations) the logarithmic slope, $`b`$, also increases as observed in the experiments. If the precise initial crumpling determines $`E_0`$ and $`C`$, then $`a`$ and $`b`$ should be correlated between runs: $`a=d_1+b(\mathrm{ln}b+d_2)`$ where $`d_2=E_0/k_BT_{\mathrm{𝑒𝑓𝑓}}\mathrm{ln}\zeta `$. The curve in Fig. 1(b) shows the fit with $`d_1=6.2`$ and $`d_2=9.7`$. We note however, that $`d_1`$ and $`d_2`$ may also vary between trials.

So far, we have reported the time dependence of the height, $`h`$, under a fixed weight. We are also interested in determining the effect of the piston mass, $`M`$, on the height. There are two problems which interfere with obtaining such data in a straight-forward manner. First, as we have emphasized, $`h`$ varies logarithmically with time and it is thus impractical to find the final height of a crumpled sheet under a given load. Second, there is large hysteresis so that different loading procedures produce different values of $`h`$. Such hysteresis is already suggested by the variation between runs of the constants, $`a`$ and $`b`$, describing the relaxation. Despite these problems we can obtain reproducible results by following two protocols for handling the crumpled sheets.

In order to avoid the problems created by the slow relaxation, we measure the height after a fixed interval of time (e.g., $`100sec`$) and use this value, instead of the unattainable asymptotic ($`t=\mathrm{}`$) value, for comparison between different loadings. Such a procedure means that we neglect relaxation processes slower than $`100sec`$. To obtain meaningful results we confine ourselves to comparing the piston height for relatively large changes in the mass, $`\mathrm{\Delta }M`$ so that $`\mathrm{\Delta }h_{\mathrm{\Delta }M}`$ (the change in height measured at $`t=100sec`$ due to the increased piston mass) is much greater than $`\mathrm{\Delta }h_{\mathrm{\Delta }t}`$ (the change in height due to waiting extra time after the mass has been added, e.g., between $`100sec`$ and $`1day`$ ): $`\mathrm{\Delta }h_{\mathrm{\Delta }M}\mathrm{\Delta }h_{\mathrm{\Delta }t}`$. In what follows we measure $`h_{100}(M)`$, the height at $`t=100sec`$ after a mass, $`M`$, has been added to the piston.

Dealing with the second problem of hysteresis requires a protocol for preparing a sheet prior to measurement. Hysteresis is evident if we start with a small initial mass and measure (as prescribed above) the height after subsequent increments $`\mathrm{\Delta }M`$ have been added until a large final mass is applied. If we then remove this load and repeat the measurement with the initial mass, we find that the newly measured height is smaller than that initially measured for the same mass.

In order to avoid this hysteresis, we have followed the protocol of Baudequin et al. for studying compressed glass wool 8 . We pre-train a sheet by first crumpling it under a mass $`M_1`$ for a period of one day. The mass is then removed and we measure $`h_{100}(M)`$ for $`M<M_1`$ (as shown for $`M_1=2.6kg`$ by the open O’s in Fig. 3(a)). This curve is reproducible (shown by the X’s) as long as $`M_1`$ is not exceeded. However, if we exceed $`M_1`$ and place a larger mass on the piston, we no longer obtain the same curve. With this larger mass $`M_2=4.5kg`$ we obtain a new curve (shown by the solid triangles) that is itself reproducible as long as $`M_2`$ is not exceeded. We continue this measurement to obtain a family of curves each of which corresponds to a pre-training mass $`M_i`$. For each curve, considered separately, we find a scaling relation:

$$h_{100}(M,M_i)=h_i+cM^\alpha $$
(3)
where $`h_i`$ is the asymptotic value of the height for the pre-treating mass $`M_i`$. Such a family of curves is shown in Fig. 3(b). We find $`\alpha =0.53\pm 0.04`$ for all pre-treating masses. This is smaller than the value $`\alpha =0.67`$ found by Baudequin et al. 8  for fibrous glass wool. The inset to Fig. 3(b) shows a scaling relation between the asymptotic height, $`h_i`$, and $`M_i`$: $`h_iM_i^\gamma `$ with $`\gamma 0.8`$. Despite the problems associated with the logarithmic relaxation and hysteresis, Fig. 3(b) shows that it is possible, under well-defined conditions, to obtain reproducible data on crumpled sheets.

We can predict the scaling relation between $`M`$ and $`h`$ if we assume a small volume fraction and a uniform network of ridges between nearly flat facets of characteristic size $`X`$ that depends on compression. The energy of each facet is proportional to $`\kappa (X/\delta )^{1/3}`$, where $`\kappa `$ is the bending modulus of the sheet and $`\delta `$ is its thickness 1 ; 2 ; 3 ; 4 ; 5 . Facets are oriented at random, so there is roughly one facet per volume $`X^3`$. Thus the energy contained in the volume $`V`$ is $`EV/X^3\kappa (X/\delta )^{1/3}`$. The space occupied by one facet is large enough to contain a stack of $`X/\delta `$ facets so that $`\varphi \delta /X`$ is a volume fraction. Macroscopically, $`\varphi =(\pi \delta (L/2)^2)/(\pi h(D/2)^2)`$, where $`h`$ is the height of the container, $`D`$ is its diameter and $`L`$ is the diameter of the sheet. $`E`$ can be re-expressed:

$$E=\kappa V/\delta ^3(X/\delta )^{8/3}=(\kappa /\delta ^3)(\pi h(D/2)^2)\varphi ^{8/3}h^{5/3}.$$
(4)
On the other hand, the energy $`E`$ supplied by the applied force $`Mg`$ is given by $`E=_{\mathrm{}}^hM(h)g𝑑hMgh`$. As noted previously, we have assumed here that all the work goes into ridge energy. We conclude that $`Mghh^{5/3}`$, so that $`hM^{3/8}`$. This predicted value of $`\alpha =0.375`$ is to be compared with the measured $`\alpha =0.54`$. A more careful estimate 15  yields $`h0.31cm`$ when $`M=0.1kg`$. The predicted height, shown by the dashed line of Fig. 3(b), is a factor 1.5 to 3 lower than those measured for the smallest masses, and approaches the measured heights for the largest masses.

The computer simulations of Lobkovsky et al. 1  indicated that the asymptotic scaling worked well when the volume fraction was no more that a few tenths of a percent. However, the experiments are well outside this regime, where the volume fraction is large (over one half). The size of the observed facets is far from uniform, and this nonuniformity could alter the stored energy. Further, we neglect friction which is expected to strengthen the crumpled sheet by a factor of order unity. Plastic deformation is expected to weaken it by at most a factor of order unity. Considering these differences between the experiment and the theoretical model, the agreement between the two exponents is surprisingly good.

The data presented here suggests that the answer to our question in the first paragraph, about how the ultimate size of a sheet depends on the forces applied, has a complicated answer. The $`\mathrm{log}t`$ dependence found in the relaxation suggests that there is dissipation, either due to friction or to plastic flow in regions of large curvature. It is unclear to what extent an analysis based only on conservative forces, which does not consider such effects, can fully describe the energetics of crumpling in real materials. Despite that caveat, our data nevertheless do support a scaling picture of the crumpling process: after a sheet has been appropriately pre-treated, the height approaches an asymptotic value as a power-law.

The results suggest several extensions such as studying how the exponent $`\alpha `$ depends on size and thickness in different materials. The most pressing question is the origin of the logarithmic relaxation. We do not know if the behavior occurs because of plastic flow in the material or because of frictional sliding at the contacts between sheets. Neither effect is included in the theoretic or simulation studies of crumpling. In order to determine the role of plastic flow in regions of large curvature, one could measure the time-dependent relaxation in rubber sheets where such plasticity would be minimized.

This research is supported by NSF DMR-0089081 and MRSEC DMR-9808595. We thank A. Marshall for help in preparing this manuscript.


text:
# Quantum dynamics in two and three-dimensional quasiperiodic tilings

\[

## Abstract

We investigate the properties of electronic states in two and three-dimensional quasiperiodic structures : the generalized Rauzy tilings. Exact diagonalizations, limited to clusters with a few thousands sites, suggest that eigenstates are critical and more extended at the band edges than at the band center. These trends are clearly confirmed when we compute the spreading of energy-filtered wavepackets, using a new algorithm which allows to treat systems of about one million sites. The present approach to quantum dynamics, which gives also access to the low frequency conductivity, opens new perspectives in the analyzis of two and three-dimensional models.

\]

Since the discovery of quasicrystals , the nature of the eigenstates and of the energy spectrum, as well as the propagation of electrons in quasiperiodic potentials, have been the subject of many theoretical studies. One-dimensional ($`1D`$) quasiperiodic Hamiltonians have been widely investigated by means of numerical  and perturbative approaches . Their spectrum consists in a set of zero width bands associated to critical eigenstates responsible of a sub-ballistic propagation of the electrons. In higher dimensions, the situation is more complex and, given the geometrical complexity of the structures, analytical treatments are difficult. Most of the studies have focussed on topologically trivial systems whose properties can be easily extracted from the $`1D`$ case (Fibonacci quasilattices , labyrinth tiling ). It has been shown that if the propagation is always sub-ballistic, the spectrum could be either absolutely continuous, singular continuous or any mixture. Other works based on exact diagonalizations, have also been performed on Penrose-like $`2D`$ and $`3D`$ systems. In $`2D`$, there is some evidence that eigenstates are critical  and responsible of an anomalous diffusion . In $`3D`$ systems, the nature of eigenstates remains unclear . In this context, it is of great interest to determine whether critical states are generic of true quasiperiodic systems especially in the $`3D`$ case, and to study the physical properties induced by such states.

In this letter we study in a tight-binding approach, the electronic properties of the of the $`2D`$ and $`3D`$ generalized Rauzy tilings (GRT) . We discuss the nature of the eigenstates by computing their localization degree for small systems ($`5500`$ sites). This analyzis suggests that the eigenstates are critical and rather more extended at the band edges than at the band center. To go beyond these exact diagonalizations, we introduce a powerful algorithm that allows us to study the quantum dynamics for very large systems (up to $`10^6`$ sites here). Sub-ballistic propagation laws are clearly established and the energy dependence of the exponents confirms that states are more mobile (almost ballistic) close to the band edges than in the center. From this point of view, these results provides evidence of critical states in the $`2D`$ and $`3D`$ GRT. Finally, we discuss the low frequency conductivity which can be directly related to the quantum diffusion .

The GRT are codimension one quasiperiodic structures that can easily be built in any dimension $`D`$ by the standard cut and project method . Here, we focus on approximant structures that are obtained for rational cut directions. These tilings have a complex topological structure with sites of coordination number ranging from $`D+1`$ to $`2D+1`$. as displayed in Fig. 1 for $`D=2`$.

This geometrical complexity implies that they cannot be separated into a combination of $`1D`$ systems as in the Fibonacci quasilattices  or as in the labyrinth tiling . We thus expect that their spectral and dynamical properties are representative of other $`2D`$ or $`3D`$ Penrose-like quasiperiodic structures.

We consider noninteracting electrons described by the following tight-binding Hamiltonian :

$$H=\underset{i,j}{}t|ij|\text{,}$$
(1)
where $`t=1`$ is the hopping energy between nearest neighbor sites, and where $`|i`$ denotes an orbital localized on the site $`i`$. As discussed in , the Hamiltonian matrix with periodic boundary conditions is straightforwardly written as a band (Toeplitz-like) matrix. This is due to a natural indexing of the sites known as the conumbering which classify them according to their local environment. The study of the spectrum of $`H`$ has already been presented in Ref. , and we focus here on the nature of the eigenstates. Therefore, we have computed by exact numerical diagonalization, the participation number $`P`$ of each normalized eigenstate $`|\psi `$ defined by : $`P(\psi )=\left(_i|i|\psi |^4\right)^1`$. For a given system size, this quantity measures the localization degree of the state $`|\psi `$ considered. We display in Fig. (2), for the $`2D`$ and $`3D`$ GRT, the participation ratio $`p=P/N`$ as a function of the energy ($`N`$ is the number of sites per unit cell). Qualitatively, we observe in both systems that, despite strong fluctuations, the most extended eigenstates are globally located near the spectrum edges and the more localized ones in the central part.

Another interesting information is given by the scaling of the mean participation number (averaged over the whole spectrum) versus the number of sites $`N`$. Indeed, when the system size increases, the mean participation number converges towards a finite value if the states are exponentially localized and goes to infinity otherwise. As shown in Fig. 3, this mean participation number behaves as $`N^\gamma `$ with $`\gamma _{2D}0.965\pm 0.001`$ and $`\gamma _{3D}=0.975\pm 0.003`$, indicating a weakly critical behaviour of the eigenstates. Note that this is only an average behaviour that conceals strong fluctuations from one energy to another as it has already been shown in Penrose-like tilings .

In the $`2D`$ GRT, the exponent is greater than the one obtained for the octagonal tiling  (codimension 2) where $`\gamma =0.87\pm 0.05`$. This difference can be understood by invoking the codimension of the structures that gives an indication on the geometrical complexity in terms of local environment. Intuitively, we indeed expect low codimension structures to be weakly quasiperiodic, and to have weakly critical states. The recent level statistics analyzis done in the $`2D`$ GRT  also corroborates this assumption. Finally, up to the numerical precision, we have $`\gamma _{2D}<\gamma _{3D}`$. It means that for a fixed codimension (here 1), the localization of the eigenstates increases when the dimension is lowered, as it is the case in disordered systems. To confirm the trends observed with diagonalizations, we have studied the spreading of energy-filtered states on much larger systems of about one million sites.

The diffusion of wave packets can be measured by the mean square spreading :

$`X^2(E,t)`$ $`=`$ $`(\widehat{X}(t)\widehat{X}(0))^2_E`$ (2)
$`=`$ $`A^{}(t)A(t)_E\text{,}`$ (3)
where $`A(t)=[\widehat{X},\text{e}^{\text{i}Ht}]`$. $`\widehat{X}(t)`$ is the Heisenberg representation of the position operator along the $`x`$ direction and $`\mathrm{}_E`$ denotes the average value over the eigenstates of energy $`E`$. We stress that $`X^2(E,t)`$ evolves with time even though it is a mean value over eigenstates.

$`X^2(E,t)`$ is related to $`\sigma (E_\mathrm{F},\omega )`$ the conductivity at frequency $`\omega `$ and Fermi energy $`E_\mathrm{F}`$ . In the limit of small frequencies one has :

$$\sigma (E_\mathrm{F},\omega )_\mathrm{\Delta }=e^2\underset{0}{\overset{\mathrm{}}{}}\text{e}^{\text{i}\omega t}n(E_\mathrm{F})\frac{\text{d}^2X^2(E_\mathrm{F},t)}{\text{d}t^2}_\mathrm{\Delta }dt$$
(4)
where $`n(E_\mathrm{F})`$ is the density of states at the Fermi energy, $`e`$ is the electron charge, and $`\mathrm{}_\mathrm{\Delta }`$ means an average over a range $`\mathrm{\Delta }`$ of values of $`E_\mathrm{F}`$. Note that at zero frequency, the integral in (4) is equal to the diffusivity $`D(E_\mathrm{F})=\text{lim}_t\mathrm{}\text{d}X^2(E_\mathrm{F},t)/\text{d}t`$ and one recovers the Einstein formula $`\sigma (E_\mathrm{F},0)=e^2n(E_\mathrm{F})D(E_\mathrm{F})`$.

To investigate the quantum diffusion, we consider an approximant of the GRT with a sufficiently large unit cell so that the eigenstates are Bloch waves characterized by a wave vector $`𝐤`$. As for the exact diagonalizations we restrict our study to the $`𝐤=0`$ subspace, i.e. to the eigenstates which are invariant by translation from one unit cell to another. Note that for an approximant the operator $`A(t)`$ is invariant by translation of a unit cell and thus lets invariant the subspace of $`𝐤=0`$ states. We calculate $`X_0^2(E,t)`$ defined by :

$`X_0^2(E,t)`$ $`=`$ $`A^{}(t)A(t)_{E,𝐤=0}`$ (5)
$`X_0^2(E,t)`$ $`=`$ $`{\displaystyle \frac{\text{Tr}\left\{A^{}(t)\delta (EH)A(t)\right\}}{\text{Tr}\left\{\delta (EH)\right\}}}`$ (6)
Actually, $`X_0^2(E,t)`$ is a good approximation of $`X^2(E,t)`$ provided that $`t`$ is lower than the time required to significantly couple orbitals separated by a distance $`l_x/2`$ during the time evolution. In this study, we have always considered wave packet evolutions for such times. The operator trace Tr is restricted to the subspace $`𝐤=0`$. An efficient trick to evaluate the traces is to use states $`|\psi `$ having a constant presence probability on each site $`|i|\psi |^2=1/N`$ and a random phase on each orbital of the unit cell. We thus approximate $`\text{Tr}\{B\}`$ by $`N\psi |B|\psi `$. In fact, we should perform an average over several random phases configurations but we have checked that, thanks to the great number of sites in the unit cell, we only need to consider one random phase state $`|\psi `$ to get a good estimate of the traces. Moreover, to eliminate fluctuations, we have performed a convolution of the $`\delta (EH)`$ by a Gaussian operator of typical width $`\mathrm{\Delta }=W/100`$ where $`W`$ is the total bandwidth . This is equivalent to compute the spreading of energy filtered wavepackets. These parts of the calculation are rather standard and in fact the crucial step is to evaluate $`A(t)|\psi `$.

We start from the expansion of the evolution operator $`\text{e}^{\text{i}Ht}`$ using a basis of Chebyshev polynomial $`P_i(E)`$ of the first kind. These polynomials are associated to the weight $`\rho (E)=1/\left(\pi \sqrt{4b^2(Ea)^2}\right)`$. This expansion converges provided the spectrum of $`H`$ (computed here using the standard recursion method ) is included in the interval $`[a2b,a+2b]`$ . The decomposition of the evolution operator is written as :

$$\text{e}^{\text{i}Ht}=\underset{i=0}{\overset{\mathrm{}}{}}c_i(t)P_i(H)\text{,}$$
(7)
with $`c_p(t)=\text{i}^ph_pJ_p(2bt)`$ where $`h_0=1`$, $`h_{p>0}=1/2`$, and where $`J_p(x)`$ is the Bessel function of order $`p`$ . At large $`p`$ and fixed $`x`$, one has : $`J_p(x)1/\sqrt{2\pi p}(x\text{e}/2p)^p`$ so that the sum in Eq. (7) converges quickly and can be truncated. The polynomials $`P_i`$ obey the following relation :

$$bP_{i+1}(H)=(Ha)P_i(H)bP_{i1}(H)\text{,}$$
(8)
with $`P_0(H)=1`$ and $`P_1(H)=(Ha)/2b`$. Then, introducing $`|\varphi _i=P_i(H)|\psi `$ and $`|\psi _i=[\widehat{X},P_i(H)]|\psi `$ and using Eqs. (7) and (8), one readily obtains :

$`b|\varphi _{i+1}`$ $`=`$ $`(Ha)|\varphi _ib|\varphi _{i1}`$ (9)
$`b|\psi _{i+1}`$ $`=`$ $`(Ha)|\psi _ib|\psi _{i1}+[\widehat{X},H]|\varphi _i\text{.}`$ (10)
Solving this latter system allows one to compute efficiently :

$$A(t)|\psi =\underset{i=0}{\overset{\mathrm{}}{}}c_i(t)|\psi _i\text{,}$$
(11)
and thus $`X_0^2(E,t)`$ with an arbitrary accuracy given by the order of the truncation.

All the results discussed here have been obtained in the $`22^{nd}`$ approximant (755476 sites) for the $`2D`$ GRT and in the $`21^{st}`$ approximant ($`1055026`$ sites) for the $`3D`$ GRT. We have computed, for both tilings and for several energy-filtered wave packets, the time dependence of $`R_0^2(E,t)=X_0^2(E,t)+Y_0^2(E,t)+Z_0^2(E,t)`$ where $`x,y`$ and $`z`$ are three orthogonal directions ($`Z_0=0`$ in $`2D`$).

As displayed in Fig. 4, one has $`R_0^2(E,t)t^{2\beta (E)}`$ at large times with $`0<\beta (E)1`$. The observation of such a sub-ballistic propagation for large approximant is consistent with the existence of critical eigenstates. The energy dependence of the diffusion exponent $`\beta `$ shown in Fig. 5 also clearly confirms the tendency observed with exact diagonalizations since the diffusion exponent is, in both systems ($`2D`$ and $`3D`$), larger at the band edges than at the spectrum center. In the $`3D`$ GRT, we have even obtained an almost ballistic propagation $`(\beta 1)`$. In addition, the shape of the density of states is nearly free- electron like at band edges which is also consistent with the less localized eigenstates at these energies. If this general trend is similar to what is observed in the octagonal tiling  and in the Penrose tiling , it is the opposite of what occurs in the $`3D`$ icosahedral Amman-Kramer tiling . Note also that in disordered systems, states at band edges are more localized than in the band center. The parameters that determine the nature of band edges states are obviously not yet understood.

We would like to emphasize that the scaling laws are valid for $`R_0^2(E,t)`$ varying by $`34`$ orders of magnitude in $`2D`$ and $`23`$ orders of magnitude in $`3D`$. In both cases the number of sites $`N`$ contained in a volume $`R_0^2(E,t)`$ in $`2D`$ and $`R_0^3(E,t)`$ in $`3D`$ varies by $`34`$ orders of magnitude. This is much better than with exact diagonalization for which the typical variation of $`N`$ is $`1`$ order of magnitude.

The power-law behaviour obtained for $`R_0^2(E,t)`$ allow us to explicitely calculate the low frequency conductivity. Indeed, assuming $`X^2(E_\mathrm{F},t)C(E_\mathrm{F})t^{2\beta (E_\mathrm{F})}`$, one obtains, using Eq. (4) :

$$\sigma (E_\mathrm{F},\omega )e^2n(E_\mathrm{F})C(E_\mathrm{F})\mathrm{\Gamma }\left(2\beta +1\right)(\text{i}/\omega )^{2\beta 1}\text{,}$$
(12)
where $`\mathrm{\Gamma }`$ is the Euler Gamma function. For the GRT, the diffusion exponent $`\beta (E_\mathrm{F})`$ is always greater than $`1/2`$ (superdiffusive propagation) so that although the system is non periodic, its dc conductivity $`\sigma (E_\mathrm{F},\omega =0)`$ is infinite for all $`E_\mathrm{F}`$. Note that experiments are qualitatively explained if $`\beta <1/2`$  but such values of $`\beta `$ probably require a strong quasiperiodic potential.

To conclude, we would like to point that the calculation of quantum diffusion presented here, implies a computational load proportionnal to the number $`N`$ of sites of the system whereas exact diagonalizations implies a load which is of order $`N^3`$ or $`N^2\mathrm{ln}N`$ for the best algorithms. This is the main reason of the great efficiency of the present method and of its ability to treat large systems. Note that in the context of large scale electronic structure calculations, order $`N`$ algorithms raise much interest and can also be much more efficient than exact diagonalizations . As a consequence, this approach should open new perspectives in the investigation of $`2D`$ and $`3D`$ models.

We thank Claude Aslangul for a careful reading of the manuscript.


text:
# Unconventional Superconductivity in UPd2Al3 from Realistic Selfconsistent Calculations

\[

## Abstract

Realistic selfconsistent calculations of unconventional superconductivity in a heavy-fermion material are reported. Our calculations for UPd2Al3 start from accurate energy band dispersions that are computed within the local spin-density functional theory and provide Fermi surfaces in agreement with experiment. Using physically motivated, realistic pairing potentials it is shown that the superconducting gap has two lines of nodes around the $`z`$-axis, thus exhibiting $`d`$-wave symmetry in the $`A_{1g}`$ representation of the $`D_{6h}`$ point group. Our results suggest that in a superconductor with gap nodes, the prevailing gap symmetry is dictated by the constraint that nodes must be as far as possible from high-density areas.

\]

Superconductivity (SC) in heavy-fermion (HF) materials exhibits a fascinating complexity of phenomena . The nature of the SC state is mostly anticipated to be unconventional, i.e., an additional symmetry is broken in the SC state, and the order parameter is therefore not the conventional spin singlet $`s`$-wave type . The identification of the order parameter symmetry is an essential issue in current studies of unconventional SC. A common approach to address the symmetry of the order parameter is through investigations of the SC gap, which is expected to have the same symmetry as the order parameter. While progress has been made for high-$`T_c`$ superconductors (HTSC), the identification of the gap symmetry remains an open question for most HF materials. Unfortunately, the eventual identification of the gap symmetry in a HF compound is not sufficient to identify the pairing mechanism. Additional elements are necessary and in the case of HF materials, it was suggested that spin-fluctuations - either in the itinerant  or localized limit as magnetic excitons  \- may replace the phonons as mediators of the pairing.

In order to identify the gap symmetry and clarify its relationship to the pairing mechanism, selfconsistent solutions of the SC gap equations with physical momentum dependent pairing potentials and accurate band structures are unavoidable. Such type of calculations have been achieved only recently for HTSC and ruthenates where usually two-dimensionality is assumed and a tight-binding dispersion of one or two electron bands is considered . In the case of HF’s the situation is much more complex. In addition to the three-dimensionality, we must deal with several anisotropic bands that cross the Fermi level ($`E_F`$) producing numerous highly anisotropic Fermi surface (FS) sheets. Due to the complicated structure of the bands in HF materials, simple tight-binding fits are impossible to obtain. As a consequence, no realistic selfconsistent calculations of SC have been achieved so far in HF materials.

Here we report for the first time a direct combination of relativistic band-structure calculations with realistic selfconsistent calculations of SC in a HF material. We focus on UPd2Al3 which is a fascinating HF superconductor having a moderately large specific heat coefficient $`\gamma =140`$ mJ/mol K2, and a relatively high critical temperature $`T_c=2`$ K . It orders antiferromagnetically below $`T_N=14.3`$ K with an ordered moment of 0.85 $`\mu _B`$/U-atom , which is large compared to the moments of other HF superconductors as, e.g., UPt3, which has a moment of only 0.03 $`\mu _B`$. A further anomalous feature is that the antiferromagnetic (AFM) order coexists with SC below 2 K. Coexistence of magnetism and SC was observed for other materials, e.g., containing $`4f`$ elements, but in those cases the magnetism is due to the localized $`4f`$ electrons far from $`E_F`$, whereas the SC is carried by itinerant electrons at $`E_F`$. Conversely, for UPd2Al3 most of the recent studies reveal that the SC, magnetic order, as well as HF behavior all involve the uranium $`5f`$ states . Moreover, some striking similarities with cuprate HTSC emerged recently. Neutron scattering experiments  have reported a distinct resonant mode that develops below $`T_c`$ as in HTSC. Tunneling spectroscopy experiments  in the SC state reported peculiar dip-hump features in the tunneling spectrum, similar to those observed in HTSC.

We solve selfconsistently the simplified gap equation starting from energy dispersions which we computed within the framework of density-functional theory in the local spin-density approximation (LSDA). In our computational approach we adopt several assumptions. First, we approximate the multiband gap equation  by an effective singleband equation. This is justified because the bands at $`E_F`$ have mainly $`5f`$ character, and there is thus no need to discriminate different band characters. The pairing potential is therefore independent of the band character. In the multiband Hamiltonian  the fermion operators for the different band characters become thereby identical, and an effective singleband Hamiltonian, with the energy dispersion $`ϵ_𝐤`$ replaced by $`_nϵ_{n𝐤}`$ ($`n`$ is the band index), results. We obtain consequently one gap symmetry defined for the whole FS. This corresponds to the experimental studies that identify a global gap symmetry and not one for every FS sheet . Second, since it is experimentally firmly established that UPd2Al3 is a spin singlet superconductor , we explore in our calculations only singlet states. Third, we adopt an adiabatic approach in our investigation of the SC order parameter. This means that we neglect the renormalization of the AFM bands by the SC, which is justified because the magnitude of the energy scale of SC is an order of magnitude smaller than that of the AFM order. Fourth, in UPd2Al3 there exist certainly many-particle effects that are esponsible for the high effective mass, which can cause a renormalization of the energy bands. Such processes are not taken into account in our LSDA calculations. Previous studies demonstrated a very close correspondence of the measured  and calculated  de Haas-van Alphen frequencies, and consequently the FS of UPd2Al3 is accurately described by the LSDA energy dispersions. Since the relevant physics for SC happens in the vicinity of the FS, the LSDA bands are the appropriate starting point for the study of the gap symmetry. Moreover, we focus here on a qualitative analysis of the gap symmetry avoiding quantitative aspects of SC for which the incorrect values of the effective masses in our LSDA scheme will indeed have some influence.

We computed the energy bands of UPd2Al3 with the relativistic augmented-spherical-wave method . The resulting LSDA band structure of UPd2Al3 along characteristic symmetry lines is shown in Fig. 1 for the paramagnetic and AFM state. The FS resulting from our AFM bands consists of four types of sheets, and is practically identical to that computed previously , as is also the ordered AFM moment of 0.81 $`\mu _B`$, which agrees with the experimental moment of 0.85 $`\mu _B`$ .

The simplified gap equation in the $`T=0`$ regime is

$$\mathrm{\Delta }_𝐤=\frac{1}{2}\underset{𝐤^{}}{}V(𝐤,𝐤^{})\mathrm{\Delta }_𝐤^{}(X_𝐤^{}+\mathrm{\Delta }_𝐤^{}^2)^{1/2},$$
(1)
where $`V(𝐤,𝐤^{})`$ is the pairing potential, and $`X_𝐤^{}=(_nϵ_{n𝐤})^2`$ is the superposition of all the relevant bands (dotted lines in Fig. 1). The momentum summation is over the full Brillouin zone (BZ). An unphysical property of Eq. (1) may occur when summed energies are zero, but not the individual energies. To avoid such situation we approximate $`X_𝐤_nϵ_{n𝐤}^2`$. We verified that this approximation did not affect the resulting gap symmetry. We consider two types of pairing kernels that were shown to produce unconventional order parameters. The small-$`q`$ pairing potential adopts at small wavevectors the following form in momentum space $`V(𝐤,𝐤^{})=`$ $`\frac{V}{𝐪_c^2+(𝐤𝐤^{})^2}+`$ $`\mu ^{}(𝐤𝐤^{})`$. This kernel is characterized by a smooth momentum cut-off $`𝐪_c`$ which selects the small-q processes in the attractive part . At larger wavevectors the repulsive Coulomb pseudopotential $`\mu ^{}(𝐤𝐤^{})`$ prevails. This type of attractive interaction may be due to phonons if screening with short range Hubbard-like Coulomb terms is involved . It has been claimed that this type of kernel may account for unconventional SC in HTSC  and in other materials, including HF’s . Agterberg et al.  have recently suggested for HF superconductors a multipocket FS with an attractive intrapocket and repulsive interpocket potential. Such situation results naturally from a potential dominated by small-q attractive pairing as the one considered here, if the FS contains multipockets. We also considered the case of pairing mediated by spin-fluctuations using a phenomenological Millis-Monien-Pines pairing potential  which has the form $`V(𝐤,𝐤^{})=\frac{V}{𝐪_c^2+(𝐤𝐤^{}𝐐)^2}`$ where $`𝐐=(0,0,\pi /2c)`$ in UPd2Al3 ($`c`$ is the $`c`$-axis lattice constant). This type of pairing has been suggested for virtually all unconventional SC’s, including UPd2Al3 .

To obtain selfconsistent gap solutions with such momentum dependent kernels we have used a Fast-Fourier-Transform technique. The problem is solved iteratively on a symmetric part of the BZ that we discretize with a $`128\times 128\times 128`$ momentum grid. For each point of our 3-D grid we have our LSDA energy bands as an input. Within our procedure the momentum space problem is fully resolved numerically without any simplification or bias on the resulting gap symmetry.

In UPd2Al3, NMR measurements revealed a singlet SC state  with a gap of even parity. Because UPd2Al3 and thus our LSDA bands obey the hexagonal $`D_{6h}`$ point group symmetry, the only even parity accessible gap states with nodes should transform according to the irreducible representations of $`D_{6h}`$  shown schematically in Fig. 2. As a first step in our iteration scheme, the initial gap configuration is chosen randomly and the system is totally free to converge within the iteration cycle towards the most favorable representation. In a second step we bias the system by adopting as initial gap configuration each of these representations and try to force convergence towards the chosen representation. In that way, we check the ground state character of the solution obtained in the first step. The precision of our numerical calculation is such that it allows to distinguish which representation corresponds to the lowest free energy. The lowest-in-energy representation, e.g., has a free energy that is 10% less than the next higher one. The second step in addition allows us to identify any eventual degeneracy between the various accessible gap symmetries.

The central result of our calculations is that the only accessible gap states with nodes in UPd2Al3 transform according to the fully symmetric irreducible representation $`A_{1g}`$. This is true irrespective of the precise momentum structure of the pairing kernel. We show in Fig. 3 some examples of our selfconsistently obtained gap functions in the plane $`\mathrm{L}\mathrm{A}\mathrm{L}\mathrm{L}\mathrm{M}\mathrm{L}`$ which contains the $`z`$-axis (the $`\mathrm{A}\mathrm{\Gamma }\mathrm{A}`$ axis) and in the plane obtained by a $`\pi /6`$ or $`\pi /2`$ rotation around the $`z`$-axis (the plane $`\mathrm{H}\mathrm{A}\mathrm{H}\mathrm{H}\mathrm{K}\mathrm{H}`$). All selfconsistent solutions shown in fig. 3 have two lines of nodes practically perpendicular to the $`\mathrm{A}\mathrm{\Gamma }\mathrm{A}`$ axis, thus belonging in the $`A_{1g}`$ representation. No other representation possessing nodes was accessible for all parameters in the pairing kernels that we have investigated. To study the influence of the coexisting AFM order on SC we have made calculations using both paramagnetic and AFM bands. In both paramagnetic and AFM cases, the only solutions with gap nodes are of $`A_{1g}`$ type. Note that our results with the spin-fluctuations kernel confirm the model analysis by Huth et al. .

There is overwhelming experimental evidence that indeed UPd2Al3 has a gap with the $`A_{1g}`$ node structure. Recent tunnel measurements along the $`z`$-axis showed the absence of a node in this direction . However, the presence of nodes is definitely established in UPd2Al3, and since $`A_{1g}`$ is the only even parity gap with nodes which is nodeless along the $`z`$-direction, it was deduced a posteriori that the gap symmetry must be $`A_{1g}`$ type . In addition, the observation of a spin-fluctuations peak below $`T_c`$ at $`(0,0,\pi /2c)`$  would be forbidden by the involved coherence factors unless the gap changes sign along the $`z`$-axis as in the $`A_{1g}`$ representation . Finally, measurements of the angular dependence of the critical field indicate two lines of nodes perpendicular to the $`z`$-axis  again suggesting the $`A_{1g}`$ representation. UPd2Al3 is perhaps the only HF superconductor for which such a clear picture of the nodal structure of the gap is obtained from the experiments. We consider the agreement of our results with the experiments as strong support of our calculations.

The surprising robustness of the $`A_{1g}`$ solution results from the following general rule: If a system must choose a gap symmetry with nodes because of the repulsive effective interaction at large wavevectors (short distances) it chooses the representation that has the minimum number of nodes as far as possible from the high-density areas. In fact, the system must maximize the condensation free energy and this is obtained when there is a gap in the high-density areas and the node (gapless) areas are kept minimal. Examining the paramagnetic bands of UPd2Al3 (Fig. 1b) we see that the high-density areas near the FS due to saddle points in the band dispersions are found essentially near the A point and near the H point. The $`A_{1g}`$ representation prevails because it is the only one without a node near the A point. In the case of the AFM bands (Fig. 1c), the saddle points of the bands near the FS are found essentially in the $`z=0`$ plane (near the $`\mathrm{\Gamma }`$ point, along the $`\mathrm{\Gamma }\mathrm{K}`$, $`\mathrm{K}\mathrm{M}`$ and $`\mathrm{M}\mathrm{\Gamma }`$ symmetry lines) and near the $`\mathrm{\Gamma }`$ point along $`\mathrm{\Gamma }\mathrm{A}^{}`$ ($`\mathrm{A}^{}=\mathrm{A}/2`$). The high-density areas in the $`z=0`$ plane exclude the $`B_{1g}`$, $`E_{1g}`$ and $`B_{2g}`$ representations which have a line node in the $`z=0`$ plane. From the remaining $`A_{1g}`$, $`E_{2g}`$ and $`A_{2g}`$ representations, $`A_{1g}`$ is the only one without a line of nodes parallel to the $`\mathrm{A}\mathrm{\Gamma }\mathrm{A}`$ axis which would cross the high-density areas at $`\mathrm{\Gamma }`$ and near $`\mathrm{\Gamma }`$ along $`\mathrm{\Gamma }\mathrm{A}^{}`$. Also, due to the AFM symmetry a high density occurs again at $`\mathrm{A}`$, as $`\mathrm{A}`$ is equivalent to $`\mathrm{\Gamma }`$. The node perpendicular to the $`\mathrm{\Gamma }\mathrm{A}`$ path in the $`A_{1g}`$ representation is close to $`\mathrm{A}^{}`$ and therefore does not cross the high-density areas near $`\mathrm{\Gamma }`$ and $`\mathrm{A}`$. Note that $`A_{2g}`$, $`B_{1g}`$ and $`B_{2g}`$ are also handicapped by the fact that they have more nodal areas.

To stress the generality of this argument we reconsidered the case of cuprate HTSC. Saddle points produce high-density areas in the vicinity of $`(0,\pi /a)`$ and symmetry related areas ($`a`$ is the lattice constant in the basal plane). The $`d_{x^2y^2}`$ gap symmetry corresponds to nodes along the $`(\pi /a,\pi /a)`$ direction as far as possible from the high-density $`(0,\pm \pi /a)`$ and $`(\pm \pi /a,0)`$ areas. Using both types of pairing kernels as adopted for UPd2Al3 (now $`𝐐=(\pi /a,\pi /a)`$) we were never able to obtain a $`d_{xy}`$ solution for which the nodes cross the high-density $`(0,\pi /a)`$ areas no matter the details of the interaction.

In conclusion, we have computed the first selfconsistent solutions of SC in the HF material UPd2Al3. Our calculations demonstrate that the only even parity accessible gap symmetry with nodes transforms according to the $`A_{1g}`$ irreducible representation of the $`D_{6h}`$ point group, independent of whether the pairing potential is dominated by small-q attractive processes or by spin-fluctuations. The robustness of the $`A_{1g}`$ representation is analyzed to be due to the presence of high-density regions of the relevant bands in the vicinity of the $`\mathrm{\Gamma }`$ and $`\mathrm{A}`$ points, that preclude node formation at these points. As a general rule we obtain that nodes must be as far as possible from high-density areas in the phase space.

We are grateful to H. Adrian, E. N. Economou, P. Fulde, M. Huth, M. Lang, P. Thalmeier, N. Sato and F. Steglich for illuminating discussions. We thank M. Huth for the figure of the irreducible representations of the $`D_{6h}`$ group.


text:
# CP Trajectory Diagram; — A tool for pictorial representation of CP and matter effects in neutrino oscillations —

## 1 Introduction

Probing into CP violation in the lepton sector is one of the most challenging goals in particle physics. Long baseline neutrino oscillation experiments are the prime candidates for observational means for detecting such effect. However, it has been known since sometime ago that the earth matter effect acts as a contamination to the measurement of genuine CP violating effects due to the leptonic Kobayashi-Maskawa phase in such experiments . Therefore, it is of crucial importance to achieve a complete understanding of the features of interplay between the genuine CP and the matter effects to have a realistic design of experiments for measuring the CP violating phase.

We describe in this article a new powerful tool which we call “CP trajectory diagram in bi-probability space” . It enables us to represent pictorially the three effects, the effects of (a) genuine CP violation due to the $`\mathrm{sin}\delta `$ term, (b) CP conserving $`\mathrm{cos}\delta `$ term, and (c) fake CP violation due to earth matter, separately in a single diagram. By using the CP trajectory diagram we observe that there is a two-fold ambiguity in the determination of $`\delta `$ which is related with the sign of $`\mathrm{\Delta }m_{13}^2`$. As described in detail in ref. , this is a remnant of the approximate degeneracy in the vacuum oscillation probability under the transformations ($`\delta \pi \delta `$) and ($`\mathrm{\Delta }m_{13}^2\mathrm{\Delta }m_{13}^2`$).

We then discuss principles of tuning beam energy for a given baseline distance to resolve the ambiguity, and to maximize the CP-violating effect. Finally, we point out that the ambiguity may be resolved with 50 % chance in the super-JHF experiment  with a megaton class water Cherenkov detector. It is quite contrary to the conventional belief that the sign of $`\mathrm{\Delta }m_{13}^2`$ can not be determined with such a short baseline as $`L=300`$ km. In a companion article, which is a report for the Proceedings of TAUP2001 , we discuss further physical implications of our results. In particular, we explore the possibility of an in situ simultaneous measurement of $`\delta `$ and the sign of $`\mathrm{\Delta }m_{13}^2`$ in a single experiment.

## 2 CP trajectory diagram in bi-probability space

We now introduce the CP trajectory diagram in bi-probability space spanned by $`P(\nu )P(\nu _\mu \nu _e)`$ and $`P(\overline{\nu })P(\overline{\nu }_\mu \overline{\nu }_e)`$. Suppose that we compute the oscillation probability $`P(\nu )`$ and $`P(\overline{\nu })`$ with a given set of oscillation and experimental parameters. Then, we draw a dot on the two-dimensional plane spanned by $`P(\nu )`$ and $`P(\overline{\nu })`$. When $`\delta `$ is varied we have a set of dots which forms a closed trajectory, closed because the probability must be a periodic function of $`\delta `$, a phase variable.

In fig. 1 plotted is the contours of oscillation probabilities $`P(\nu )`$ and $`P(\overline{\nu })`$ which is drawn by varying the CP violating phase $`\delta `$ from 0 to $`2\pi `$. They are averaged over the Gaussian neutrino energy distribution with peak energy of 500 MeV and width 100 MeV. As you might have guessed, these diagrams are elliptic. It is so exactly in vacuum and in a good approximation even in matter at relatively short baseline and as far as the mass hierarchy $`|\mathrm{\Delta }m_{12}^2/\mathrm{\Delta }m_{13}^2|1`$ is the case .

What does CP trajectory diagram actually represent? There is a very simple answer to this question. The lengths of major and minor axes represent the coefficients of CP violating $`\mathrm{sin}\delta `$ and CP conserving $`\mathrm{cos}\delta `$ terms, respectively, in the neutrino oscillation probability $`P(\nu _\mu \nu _e)`$. Whereas the distance between two (positive and negative $`\mathrm{\Delta }m_{13}^2`$) ellipses gives the size of the matter effect. The last point can be explicitly verified by running the same calculation with twice larger matter effects, as done in ref. . Therefore, you can see the size of these three terms just by eye when it was projected onto the CP trajectory diagram.

Fig. 1: CP trajectory in the bi-probability (given in %) plane for the baseline $`L=295`$ km and energy $`E=500`$ MeV. As indicated in the figures, the solid and the dashed lines are for $`\mathrm{\Delta }m_{13}^2>0`$ and $`\mathrm{\Delta }m_{13}^2<0`$ cases, respectively, and the dotted and the dash-dotted lines correspond to the same signs of $`\mathrm{\Delta }m_{13}^2`$ as above but with matter effect switched off. The mixing parameters are fixed as $`\mathrm{\Delta }m_{13}^2=\pm 3\times 10^3`$ eV2, $`\mathrm{sin}^22\theta _{23}=1.0`$, $`\mathrm{\Delta }m_{12}^2=5\times 10^5`$ eV2, $`\mathrm{sin}^22\theta _{12}=0.8`$, $`\mathrm{sin}^22\theta _{13}=0.05`$. We take $`\rho Y_e=1.4`$ g/cm3 where $`\rho `$ is the matter density and $`Y_e`$ is the electron fraction.

It is obvious from fig. 1 that the approximate degeneracy ($`\delta \pi \delta `$) and ($`\mathrm{\Delta }m_{13}^2\mathrm{\Delta }m_{13}^2`$) which exists in the vacuum case is lifted by the matter effect; matter helps! However, it is also clearly seen in fig. 1 that there are remaining degeneracies if the sign of $`\mathrm{\Delta }m_{13}^2`$ is not known a priori.

The above discussion assumes that the value of $`\theta _{13}`$ is known prior to the experiment. If the value is not known in advance there is an another ambiguity, the ($`\delta \theta _{13}`$) ambiguity, as discussed in detail in ref. . Therefore, there exits a combined ambiguity, ($`\delta `$ $``$ sign of $`\mathrm{\Delta }m_{13}^2`$) times ($`\delta \theta _{13}`$), which in the worst case can be 4-fold. The ambiguity was referred to as the “clover-leaf ambiguity” in TAUP2001 .

## 3 Principle of choosing beam energies for long-baseline neutrino oscillation experiments

If one wants to resolve the two-fold ambiguity ($`\delta `$ $``$ sign of $`\mathrm{\Delta }m_{13}^2`$) in a single experiment one cannot tune the experimental parameters so that the $`\mathrm{sin}\delta `$ term is maximal because then the trajectory shrinks to a straight line, closing the possibility of resolution of the ambiguity. Then, there should be a compromise. First, one can obtain the conditions to maximize lengths of major and minor axes, which lead to

$`\left({\displaystyle \frac{E}{\text{1 GeV}}}\right)_{\mathrm{cos}\delta }`$ $`=`$ $`1.13,0.47,0.29\left({\displaystyle \frac{L}{300\mathrm{km}}}\right)\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2}{3\times 10^3\mathrm{eV}^2}}\right),`$ (1)
$`\left({\displaystyle \frac{E}{\text{1 GeV}}}\right)_{\mathrm{sin}\delta }`$ $`=`$ $`0.62,0.24,0.14\left({\displaystyle \frac{L}{300\mathrm{km}}}\right)\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2}{3\times 10^3\mathrm{eV}^2}}\right).`$ (2)
It is shown in ref.  that one can compromize these two requirements in the energy region of 0.5 to 2 GeV.

## 4 Possibility of simultaneous determination of $`\delta `$ and $`\mathrm{\Delta }m_{13}^2`$ in the super-JHF experiment

We now address our final subject, probably the most important one, the possibility of simultaneous determination of $`\delta `$ and $`\mathrm{\Delta }m_{13}^2`$ in the super-JHF experiment. We assume that $`\theta _{13}`$ is known in a reasonable accuracy prior to this experiment. In fig. 2 plotted is the CP trajectory diagram in the number of appearance events plane for $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ channels. We assume a water Cherenkov detector of fiducial volume 0.9 Mton, and 4 MW of proton beam power which is planned in the JHF experiment in its phase II . We use the off-axis (OA) 3 degree beam whose neutrino energy peaks at $`E0.5`$ GeV. We refer ref.  for a detailed explanation of how the computation of number of events is done and the results for alternative choice of beams.

One can clearly see from fig. 2 that at least a half of the parameter space fulfilling the condition $`\mathrm{sin}\delta \mathrm{\Delta }m_{13}^2<0`$ does not suffer from the ambiguity problem. Then, with the 50 % chance (thanks to nature’s kind setting!) they will be able to determine the CP violating angle $`\delta `$ and the sign of $`\mathrm{\Delta }m_{13}^2`$ simultaneously in the super-JHF experiment.

Fig. 2: CP trajectory in event number plane $`N(e^{})N(e^+)`$ for OA beam 3 degree and the baseline $`L=295`$ km. The dotted circles correspond to 3 $`\sigma `$ statistical uncertainty.

## Acknowledgements

This work was supported by the Brazilian funding agency Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), and by the Grant-in-Aid for Scientific Research in Priority Areas No. 12047222, Japan Ministry of Education, Culture, Sports, Science and Technology.


text:
# Charge Echo in a Cooper-Pair Box

## Abstract

A spin-echo-type technique is applied to an artificial two-level system that utilizes charge degree of freedom in a small superconducting electrode. Gate-voltage pulses are used to produce the necessary pulse sequence in order to eliminate the inhomogeneity effect in the time-ensemble measurement and to obtain refocused echo signals. Comparison of the decay time of the observed echo signal with estimated decoherence time suggests that low-frequency energy-level fluctuations due to the $`1/f`$ charge noise dominate the dephasing in the system.

In a small-Josephson-junction circuit, we can construct an artificial two-level system, which is expected to have a long decoherence time and work as a qubit . There are two types of qubit, depending on which, either charge or phase, degree of freedom is used in the two-level system. The charge qubit, a Cooper-pair box, uses two charge states in a small superconducting electrode connected to a reservoir via a Josephson junction , while the phase qubit uses two phase states in a small superconducting loop intersected by Josephson junction(s) . Superposition of the two charge states  and that of two phase states  have been observed recently, and coherent control of the charge qubit has also been demonstrated . Being solid-state devices, these qubits are also expected to be scalable to an integrated system for quantum information processing . However, the actual decoherence time and the decoherence mechanism are not yet fully understood. Accordingly, we performed an experiment on the decoherence in a Cooper-pair box.

In this experiment, a Cooper-pair box is subjected to a sequence of gate-voltage pulses, which give rise to quantum-state control similar to those in the free-induction-decay (FID) and spin-echo experiments in the field of nuclear magnetic resonance . Although our experiment is on a single two-level system, such technique is useful for studying the effect of fluctuations, that is, temporal inhomogeneities, in the time-ensemble measurement. We observed rapid decay of the FID signal and significant recovery of the coherence signal in the echo experiment. These observations indicate that the effect of low-frequency energy-level fluctuations is dominant in the dephasing of the two-level system. Furthermore, we compare the decay of the observed “charge-echo” signal with estimated decoherence due to several possible sources, and suggest a contribution of $`1/f`$ background charge noise to the dephasing.

The Cooper-pair-box device  consists of a small “box” electrode coupled to a reservoir electrode via a Josephson junction, two gate electrodes capacitively coupled to the box, and an additional probe electrode connected to the box via a highly resistive tunnel junction (Fig. 1). All the electrodes are made of Al with a superconducting gap energy $`\mathrm{\Delta }`$ larger than the other important energy scales; thus, quasiparticle excitation can be neglected (except for the measurement process described below). Because of the charging effect, it is a good approximation of the quantum state of the device to take into account only two charge states with the lowest energies—the ones closest to the charge neutrality of the box—which differ by a single Cooper pair. For example, if the total gate-induced charge in the box, $`Q_tQ_0+C_pV_p(t)`$ ($`Q_0C_gV_g+C_bV_b`$), is close to one, the relevant charge states are $`|n=0`$ and $`|1`$, where $`n`$ is the excess number of Cooper pairs in the box and $`C_i`$ ($`i=p,g,b`$) is the capacitance between the box and each electrode voltage-biased at $`V_i`$.

On the basis of $`||0`$ and $`||1`$, the effective Hamiltonian of the two-level system is $`H=\frac{1}{2}\delta E(Q_t)\sigma _z\frac{1}{2}E_J\sigma _x`$, where $`\delta E(Q_t)4E_C(Q_t/e1)`$ is the energy difference between the two charge states, $`E_Ce^2/(2C_\mathrm{\Sigma })`$ is the single-electron charging energy of the box, $`C_\mathrm{\Sigma }`$ is the total capacitance of the box, $`E_J`$ is the Josephson energy of the junction, and $`\sigma _z`$ and $`\sigma _x`$ are Pauli matrices. Therefore, the quantum-state evolution of the Cooper-pair box is mapped on the evolution of the fictitious spin-$`\frac{1}{2}`$ state in an effective magnetic field, $`\stackrel{}{B}=(B_x,0,B_z)=(E_J,0,\delta E(Q_t))`$.

In the present device, the energy scales were determined by independent measurements as $`\mathrm{\Delta }=235\pm 10`$ $`\mu `$eV, $`E_C=122\pm 3`$ $`\mu `$eV, and $`E_J=34\pm 2`$ $`\mu `$eV . The present experiment was made at the temperature of about 30–60 mK ($`k_BT3`$–5 $`\mu `$eV). The details of the experimental setup have been published elsewhere .

Quantum-state control is realized by application on a high-speed gate-voltage pulse $`V_p(t)`$ to modulate $`Q_t`$ nonadiabatically . In terms of fictitious spin, $`B_z`$ is abruptly changed, and the spin precesses around the direction of the new magnetic field until $`B_z`$ is switched back after the pulse width $`\mathrm{\Delta }t`$, as depicted in the Bloch spheres in Fig. 1(b). We set our working point at $`Q_0=0.45e`$ where $`|\delta E(Q_0)|270`$ $`\mu `$eV $`E_J`$. A $`\pi `$-pulse is produced by a pulse height $`\mathrm{\Delta }Q_p=0.55e`$ which makes the two states degenerate (i.e., $`B_z=0`$) during the pulse (Fig. 1(b) top). We found that a nominal pulse width $`\mathrm{\Delta }t=80`$ ps is suitable for the $`\pi `$ rotation. For a $`\frac{\pi }{2}`$-pulse, we could, in principle, use a pulse with the same $`\mathrm{\Delta }Q_p`$ and with a different $`\mathrm{\Delta }t`$. However, for convenience, instead we use a pulse with the same $`\mathrm{\Delta }t`$ and a slightly detuned $`\mathrm{\Delta }Q_p`$ ($`B_z>0`$). The bottom sphere in Figure 1(b) illustrates the corresponding schematic quantum-state evolution. Though the pulse is not equivalent, e.g., to a $`\left(\frac{\pi }{2}\right)_x`$-pulse (which rotates the spin around $`x`$-axis), it creates fifty-fifty superposition of $`|0`$ and $`|1`$ and is thus sufficient for our purpose. For simplicity, in the following argument, we use ideal $`(\pi )_x`$\- and $`(\frac{\pi }{2})_x`$-pulses and neglect the deviations in the real operations, because such deviations only cause quantitative difference in the absolute signal size but not in the relative decay of the signal . The delay time between the pulses is controlled by programmed patterns in a multi-channel pulse-pattern generator and the shift in a mechanical delay line.

The quantum state after the control is measured by using a probe junction. The probe junction is voltage-biased appropriately (at $`V_b=645`$ $`\mu `$V) so that only the $`|1`$ state allows two electrons (quasiparticles) to escape from the box at rates $`\mathrm{\Gamma }_{\mathrm{qp1}}`$ and $`\mathrm{\Gamma }_{\mathrm{qp2}}`$, while the $`|0`$ state does not. The resistance of the probe junction is made to be very large ($`R_b=42.6`$ M$`\mathrm{\Omega }`$) in order to suppress Josephson coupling and to avoid too much decoherence due to the read-out. The measurement basis is $`\sigma _z`$-basis, and the quantum state is initialized very closely to $`|0`$ after the read-out. On the other hand, since the output of the single read-out is just two electrons, it is necessary to repeat the read-out process many times on identically prepared quantum states in order to obtain a detectable dc current proportional to the population of $`|1`$, i.e., $`\sigma _z+1`$. The repetition time $`T_r`$ (128 ns unless specially mentioned) is taken to be long enough ($`T_r1/\mathrm{\Gamma }_{\mathrm{qp1}}+1/\mathrm{\Gamma }_{\mathrm{qp2}}`$) for efficient measurement and good initialization for the next pulse operation but short enough for a sufficient current signal (at most $`2e/T_r`$). The probe current is measured with a time window of 20 ms. Therefore, each data point is a result of averaging over a time-ensemble of about $`10^5`$.

We first measured the free-induction decay (FID) by using two $`(\frac{\pi }{2})_x`$-pulses with a delay time $`\delta t_3`$. The first pulse prepares a superposition, and the quantum state acquires phase $`\phi =\mathrm{\Delta }E(Q_0)\delta t_3/\mathrm{}`$ during the delay time, where $`\mathrm{\Delta }E(Q_0)=\sqrt{\delta E(Q_0)^2+E_{J}^{}{}_{}{}^{2}}`$ is the energy difference between the two eigenstates and corresponds to the strength of the fictitious field. The second pulse projects the phase information on the measurement basis, i.e., $`\sigma _z`$-basis. As shown in Fig. 2(c) the signal oscillates with a period of about 15 ps, which agrees well with $`h/\mathrm{\Delta }E(Q_0)`$. Moreover, it is clear that the FID signal decays quite rapidly within a few hundred picoseconds. Compared with the previous result in which coherent oscillations between two degenerate states were observed up to about 2 ns , the FID time scale is much shorter. Since $`\frac{d(\mathrm{\Delta }E)}{dQ_0}`$ equals zero at the degeneracy point but not at the present working point and the inelastic relaxation rate of the two-level system is estimated to be much longer (see below), the short FID time can be attributed to the dephasing due to charge fluctuations. Our measurement is an average over a time-ensemble of 20 ms, so even very low-frequency fluctuations can contribute, as inhomogeneities in the ensemble, to the fast decay time in the FID experiment.

To remove such an effect, the three pulses shown in Fig. 2(b) are used for the echo experiment. The second pulse flips the possibly dephased ensemble of spins after the first delay time $`t_d`$, and during the second delay time, provided that the correlation time of the fluctuations in $`\mathrm{\Delta }E(Q_0)`$ is much longer than $`t_d`$, the spins respectively precess the same amount as before and are refocused at the end. The third pulse is again needed for the projection of the phase information onto $`\sigma _z`$.

The observed echo signal is shown in Fig. 2(d). When the position of the third pulse $`\delta t_3`$ was swept, an oscillating signal was observed only in the vicinity of $`\delta t_3=0`$ and decayed as $`|\delta t_3|`$ increases within nearly the same time scale as the FID signal did.

Decay of the echo signal as a function of the delay time $`t_d`$ is shown in Fig. 2(f). The echo amplitude $`\mathrm{\Delta }I`$ was obtained by fitting the echo signal as a function of $`\delta t_2`$ with a sinusoidal curve (Fig. 2(e)). The period of the echo signal is about half of those in Figs. 2(c) and (d), because here the difference between the first and the second delay times $`(t_d+\delta t_2)(t_d\delta t_2)=2\delta t_2`$ matters in the final phase. The slight phase offset at $`\delta t_2=0`$ is due to an offset in the delay line. In Fig. 2(f), the decay time is much longer than that of the FID signal. This implies that dephasing in the FID is mainly due to the low-frequency ($`<1/t_d`$) part of the fluctuations which is greatly cancelled out by the echo technique.

The above results lead to the question concerning the origin of the decay of the echo signal. In order to compare with theoretical estimations of a decoherence factor $`\mathrm{exp}i\phi `$, the echo amplitude data (Fig. 2(f)) and the FID data (Fig. 2(c)) are normalized by fitting with a gaussian curve and are plotted in Fig. 3 on a logarithmic scale.

The well-known $`1/f`$-like background charge noise, which is believed to be caused by charges randomly fluctuating around the device , may contribute to the decoherence. The $`1/f`$ noise measurement is usually restricted to a very low frequency below about 1 kHz, and the noise spectrum in the high-frequency range is unknown. Nevertheless, here we assume that the $`1/f`$ spectrum extends to the infinite frequency . Generally, the dephasing factor in the presence of energy-level fluctuation spectrum $`S_{\mathrm{\Delta }E}(\omega )`$, under the assumption of gaussian fluctuations, is given as

$`\mathrm{exp}i\phi (\tau )=\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\tau }𝑑t\stackrel{~}{\mathrm{\Delta }E}\right)^2\right]`$ (1)
$`=\mathrm{exp}\left[{\displaystyle \frac{1}{2\mathrm{}^2}}{\displaystyle _{\omega _m}^{\mathrm{}}}𝑑\omega S_{\mathrm{\Delta }E}(\omega )\left({\displaystyle \frac{\mathrm{sin}\left(\omega \tau /2\right)}{\omega /2}}\right)^2\right]`$ (2)
for the FID case. Here $`\stackrel{~}{\mathrm{\Delta }E}`$ is the deviation of $`\mathrm{\Delta }E`$ from its average and $`\omega _m2\pi /t_m`$ is the low-frequency cutoff due to the finite data-acquisition time $`t_m`$ (20 ms in the present work). Similarly, $`\mathrm{exp}i\phi `$ in the echo experiment becomes

$$\mathrm{exp}\left[\frac{1}{2\mathrm{}^2}_{\omega _m}^{\mathrm{}}𝑑\omega S_{\mathrm{\Delta }E}(\omega )\left(\frac{\mathrm{sin}^2\left(\omega \tau /4\right)}{\omega /4}\right)^2\right]$$
(3)
because of the cancellation of the dephasing before and after the $`\pi `$-pulse. With a $`1/f`$ spectrum, $`S_{\mathrm{\Delta }E}(\omega )\left(\frac{4E_C}{e}\right)^2S_e(\omega )\left(\frac{4E_C}{e}\right)^2\frac{\alpha }{\omega }`$, both Eqs.(1) and (2) give gaussian-type decay . For $`\alpha =(1.3\times 10^3e)^2`$, which we determined by a standard noise measurement on the present device used as a single-electron transistor, the two curves are plotted in Fig. 3, as well as another curve for echo experiment with $`\alpha =(0.3\times 10^3e)^2`$, which is a typical number given in the literature. Considering the uncertainty in the noise spectrum, it seems that the curves reasonably reproduce the significant increase of the decay time in the echo signal and the time scale of the gaussian decay qualitatively, suggesting that the $`1/f`$ noise is the main dephasing source. It is worth mentioning that, if this is true, long decoherence time required for quantum computing is only be possible after a drastic reduction of the $`1/f`$ noise. Because of the long tail to the high frequency, echo technique is not a perfect solution for suppressing the dephasing .

Dephasing due to the electromagnetic environment  is also calculated by Eq.(1) with $`S_{\mathrm{\Delta }E}(\omega )=4e^2\kappa _i^2S_{Vi}(\omega )`$, where $`S_{Vi}(\omega )\frac{\mathrm{}\omega }{2\pi }2R_{\mathrm{env}}\mathrm{coth}\left(\frac{\mathrm{}\omega }{2k_BT_{\mathrm{env}}}\right)`$ is the voltage fluctuation (Johnson-Nyquist noise) on each electrode, $`\kappa _iC_i/C_\mathrm{\Sigma }`$ is the coupling constant, and $`T_{\mathrm{env}}`$ and $`R_{\mathrm{env}}`$ are the temperature and the impedance real part of the environment, which is assumed to be ohmic. In the relevant temperature range, the dephasing gives exponential decay of $`\mathrm{exp}i\phi `$ at the rate of $`8\pi \frac{k_BT_{\mathrm{env}}}{\mathrm{}}\kappa _i^2\frac{R_{\mathrm{env}}}{R_K}`$ ($`R_K\mathrm{}/e^2`$). By assuming $`T_{\mathrm{env}}`$ and $`R_{\mathrm{env}}`$ for the pulse gate, which is connected to the pulse generator via a low-loss line, as 300 K and 50 $`\mathrm{\Omega }`$, and for the dc gate and the probe electrode, which are heavily filtered, as 100 mK and 100 $`\mathrm{\Omega }`$ ($`10^3\kappa _i1.3`$, 2.5, and 73 for $`i=p,g,b`$), we evaluated the total dephasing rate as (100 ns)-1. Note that the echo technique does not reduce the dephasing caused by white noise. Also, inelastic relaxation rate due to the environment is estimated from the standard perturbation theory as $`\frac{2\pi }{\mathrm{}^2}e^2\kappa _i^2\left(\frac{E_J}{\mathrm{\Delta }E}\right)^2S_{Vi}\left(\frac{\mathrm{\Delta }E}{\mathrm{}}\right)`$. The total relaxation rate calculated by using the above parameters yields (580 ns)-1. The total decay of the coherence due to the environment (dotted line in Fig. 3), including the dephasing and the relaxation, is much weaker than the decay of the observed echo signal. It should be mentioned, however, that in the present experiment, the contribution of residual noise in $`V_b`$ to the observed dephasing cannot be thoroughly excluded. In particular, the relatively short decay times of FID signal implies the existence of excess noise in the low-frequency range.

The coherence is also disturbed by the read-out process at the probe junction. The first electron tunnels out from the box, i.e., the system escapes from the relevant Hilbert space of the two-level system, at the rate $`\mathrm{\Gamma }_{\mathrm{qp1}}=\frac{\mathrm{\Delta }E_{\mathrm{qp1}}}{e^2R_b}W\left(\frac{\mathrm{\Delta }E_{\mathrm{qp1}}}{2\mathrm{\Delta }}\right)`$, where $`\mathrm{\Delta }E_{\mathrm{qp1}}eV_b+E_C+E_C(12Q_0/e)`$ is the energy dissipated in the quasiparticle tunneling and $`W(x)1`$ is a factor related to the BCS quasiparticle density of states . In the present device, $`\mathrm{\Gamma }_{\mathrm{qp1}}`$ is about (8 ns)-1 (dashed line in Fig. 3), which should contribute to the decay of the echo signal, but does not explain it completely.

In conclusion, we have demonstrated charge echo in a Cooper-pair box. The fact that the echo signal survives much longer than the free-induction decay signal indicates that the dephasing is mainly caused by low-frequency energy-level fluctuations. The observed dephasing behavior agrees qualitatively with the prediction of that due to $`1/f`$ background charge noise.

We thank Isaac Chuang and Jonathan Jones for suggesting the echo experiment and Daniel Esteve for valuable discussions on dephasing. This work has been supported by the CREST project of Japan Science and Technology Corporation (JST).


text:
# Polarization transfer in SIDIS for Λ production†

## Abstract

Polarization transfer from the longitudinally polarized positron to the final-state quark is considered for the underlying QCD subprocesses giving rise to the reaction $`e_{}^++pe^++\mathrm{\Lambda }_{}+X`$ at the HERMES energy.

Measuring the longitudinal $`\mathrm{\Lambda }`$ polarization in the semi-inclusive deeply inelastic scattering (SIDIS)

$`e_{}^++pe^++\mathrm{\Lambda }_{}+X,`$ (1)
that is under study in the experiment HERMES at DESY, may give important information for the determination of the twist-2 polarized fragmentation function $`G_1(z,\mu )`$. The latter can be found through the relation

$`D_{LL}F_q(x)\widehat{\sigma }\widehat{d}_{LL}G_1(z,\mu ),`$ (2)
where $`F_q(x)`$ is the unpolarized parton (quark/gluon) distribution function in the proton and the polarization transfer $`D_{LL}`$ is defined as

$$D_{LL}=\frac{\sigma _{++}+\sigma _{}\sigma _+\sigma _+}{\sigma _{++}+\sigma _{}+\sigma _++\sigma _+}\frac{\mathrm{d}\mathrm{\Delta }\sigma }{\mathrm{d}\sigma },$$
(3)
and it can be measured in the experiment. The hatted entities in (2) relate to the underlying QCD subprocess. The subscripts in $`\sigma `$’s in (3) denote the positron and $`\mathrm{\Lambda }`$ helicities.

The published HERMES measurement in the reaction (1) is $`D_{LL}=0.11\pm 0.17\pm 0.03`$exper .

As is seen from (2), it is very instructive to study $`\widehat{d}_{LL}`$ for the purpose of the data handling and analysis. The parameter can also be expressed as $`\widehat{d}_{LL}`$ = $`P_{out}^z/P_{in}^z`$ ——————————————————-
$``$ Talk given at IX Workshop on High Energy Spin Physics, Dubna, Aug. 2001
where $`P`$’s are the longitudinal components of the initial and final polarization vectors.

The elementary subprocesses under study are the elastic quark probing (Fig. 2)

$`e_{}^++qe^++q_{}^{}(\mathrm{LO})`$ (4)
and the diagrams with the gluon correction (Fig. 2)

$`e_{}^++qe^++q_{}^{}+G(\mathrm{NLO}).`$ (5)
Since we are interested mainly in the central quark production region, here we

neglect the contribution of the process $`e_{}^++Ge^++q_{}+\overline{q}`$ with the quark pair production by $`\gamma ^{}G`$ coupling. The spin transfer properties in the production of heavy $`Q\overline{Q}`$ quarks are discussed in ref. YuA .

The purpose of this study is to look for kinematic regions where the polarization transfer $`\widehat{d}_{LL}`$ gets enhanced in the subprocesses of fig. 1 and fig. 2. Being translated onto the hadronic level, these kinematic regions, if any, can serve as guidelines for the experiment HERMES. The interval of DIS variables ($`Q^2,y`$ and other) was chosen to be typical for HERMES.

The calculations are based on the following expression for the scattered quark vector polarization

$$P_{out}^z=\frac{\left(1+P_z^{in}\right)\left(|M_{++}|^2|M_+|^2\right)+\left(1P_z^{in}\right)\left(|M_+|^2|M_{}|^2\right)}{\left(1+P_z^{in}\right)\left(|M_{++}|^2+|M_+|^2\right)+\left(1P_z^{in}\right)\left(|M_+|^2+|M_{}|^2\right)},$$

which can be derived from the final-state spin-density matrix.

The properties of the elastic subprocess (4) are plotted in fig. 4. Both the differential cross section $`\mathrm{d}\sigma `$ and polarized differential cross section $`\mathrm{d}\mathrm{\Delta }\sigma `$ decrease with increasing $`Q^2`$ so that the polarization transfer $`\widehat{d}_{LL}`$ increases rapidly. It achieves the value $``$0.8 at the extreme value of $`Q^24`$ GeV2. For modest momentum transfers, say $`Q^22`$ GeV2, $`\widehat{d}_{LL}`$ has the ’normal’ values of $``$ 0.2.

Fig. 4 represents the parameter $`\widehat{d}_{LL}`$ calculated for the reaction (5) at four different values of the angle $`\varphi `$ between two characteristic planes of the reaction – leptonic plane and the plane formed by the virtual photon and the outgoing quark. As is known, the polarization properties of the reactions 2 $``$ 3 depend essentially, in general, on the rotational angle $`\varphi `$.

The distributions of $`\widehat{d}_{LL}`$ in fig. 4 are shown for four rotations of the planes with an increment equal to 90o. The lowest values of $`\widehat{d}_{LL}`$ are obtained at $`\varphi `$ = 180o when the normals to the said plains are anti-collinear. The more sizeable values of $`\widehat{d}_{LL}`$ appear when the normals make right angles (90o or 270o) – up to $``$0.4. And the most favorite kinematic region is predicted for the collinar normals ($`\varphi `$ = 0) where the polarization transfer increases up to 0.6.

So coming down to the experiment (1), one can expect the sizeable longitudinal component of the $`\mathrm{\Lambda }`$ polarization induced by the polarized positron under the kinematic configurations in which the normals to the leptonic plane and the plane ($`\gamma ^{}`$, $`\mathrm{\Lambda }`$) are close to be collinear. The numeric estimates made at the HERMES energy can certainly be affected by the higher-order corrections.