text:
# The origin of diffusion: the case of non chaotic systems

## I Introduction

The presence of randomness in many natural processes is traditionally ascribed to the lack of perfect control on observations and experiments. Then, usually,“disorder” and “irregular behaviors” are considered only as unavoidable factors which, in phenomena modellizations, can be described effectively through external sources (random numbers in computer simulations, noise in stochastic equations etc.). However randomness may also have another origin associated to local instabilities of the deterministic microscopic dynamics of a system, i.e. microscopic chaos (MC). For instance, macroscopic diffusion of a particle in a fluid can be interpreted as a consequence of chaos on microscopic scales generated by collisions of the particle with the surrounding others.

The possibility of describing the Brownian behavior in terms of deterministic motion has stimulated a large amount of work in the field of deterministic diffusion to study purely deterministic systems exhibiting asymptotically a linear growth of the mean square displacement ChaoDiff ; ChaoDiff1 ; KD . The deterministic interpretation of transport phenomena establishes a close connection between transport and chaos theory, and, at least in principle, macroscopic properties such as transport coefficients (viscosity, thermal and electrical conductivity, diffusion etc.) could be directly computed through indicators of chaos GasNik ; Morriss ; Dorfman  or by applying the standard techniques of chaos theory such as periodic orbit expansion Artuso ; Vance ; CEG95 ; Gaspard . However, considering randomness and transport in general as signatures of MC raises subtle and fundamental issues in statistical physics that have to be carefully discussed an investigated. First, there is the main problem of conceiving and realizing experiments able to prove, beyond any doubt, the existence of MC in real systems. An experiment in this direction was attempted in ref. Nature  but the conclusions were object of criticism by many authors, as we discuss later. Besides, for real systems, the definition of MC appears to be ambiguous. Indeed, real systems involve a large number of degrees of freedom and their model description entails taking the thermodynamic limit but, in this limit, the concept of Lyapunov exponent becomes metric dependent so not well established GS . Moreover, the large number of degrees of freedom poses severe limitations also in the applicability of standard techniques to distinguish between chaotic and stochastic signals GP ; CP ; GW  when data are extracted from real experiments. In fact, the noise-chaos distinction requires a extremely fine resolutions in the observations that is impossible to reach when one deals with systems with large dimensionality DCH ; DC ; Cencio .

References Nature ; Experiment  summarized an ingenious experiment devised to infer the presence of chaos on microscopic scales. In the experiment the position of a Brownian particle in a fluid was recorded at regular time intervals. The time series were processed by using the same techniques developed for detecting chaos from data analysis CP ; GW . Authors measured a positive lower bound for the Kolmogorov-Sinai entropy, claiming that this was an experimental evidence for the presence of microscopic chaos. However in the papers DCH  and DC  it was argued that a similar result holds for a non-chaotic system too, hence the experiment cannot provide a conclusive evidence neither for the existence of microscopic chaos nor for the relevance of chaotic behavior in diffusion phenomena. Such papers shown that, through the analysis performed in Refs. Nature ; Experiment , there is no chance to observe differences in the diffusive behavior between a genuine deterministic chaotic systems as the 2-d Lorentz gas with circular obstacles Dorfman  and its non chaotic variant: the wind-tree Ehrenfest model. The Ehrenfest model consists of free moving independent particles (wind) that scatter against square obstacles (trees) randomly distributed in the plane but with fixed orientation. Due to collisions, particles undergo diffusion, however their motion cannot be chaotic because a reflection by the flat planes of obstacles does not produce exponential trajectory separation; the divergence is at most algebraic leading to zero Lyapunov exponent. Such considerations can be extended to squares obstacles with random orientations and to every polygonal scatterers, so there is a whole class of models where diffusion occurs in the absence of chaotic motion. A difference between the 2-$`d`$ chaotic and non-chaotic Lorentz gas relies on the presence of periodic orbits DC .

The question that naturally arises “What is the microscopic origin of the diffusive behavior?” remains still open. In this paper, focusing on this problem, we extend the analysis of Dettman and Cohen showing that standard diffusion may also occur in one-dimensional models where every chaotic effect is surely absent. These systems, being one-dimensional, can be much easily studied than wind-tree models without spoiling the essence of the problem. Accurate results are obtained with relatively small computational efforts, and the interpretation of the results is clearer since we are able to better control and quantify the effects of the local instabilities. We shall see that, in such models, the relevant ingredients to obtain diffusion in the absence of deterministic chaos is the combined action of a “finite size instability”, quenched disorder and a quasi-periodic perturbation. Where by “finite size instability” we mean that infinitesimal perturbations are stable, while perturbations of finite size can grow algebraically. The interplay between these three features guarantees the system to be still diffusive as it would be chaotic. This outcome is another striking effect of finite size instabilities which make the behaviors of certain non-chaotic system similar, to some extent, to those generated by genuine chaotic systems (PLOK ; Cecconi ; Torcio ; LRB ). In ref. zasl  the reader can find a study of an interesting system with weak mixing properties and zero Lyapunov exponent, i.e. an intermediate case between integrability and strong mixing.

The paper is organized as follows. In section II, we first recall the general features of $`1`$-$`d`$ models showing deterministic diffusion. Then we present two non-chaotic, one-dimensional models that exhibit deterministic diffusion when properly perturbed. In section III we characterize numerically the diffusion properties of the models and their dependence on the external perturbation. Section IV shows how these models can be viewed as non-chaotic space-time billiards, and discusses in more detail the role of periodic orbits. Section V is devoted to conclusions and remarks.

## II Non-Chaotic Models for Deterministic Diffusion

Motivated by the problem of deterministic diffusion in non-chaotic Lorentz systems (i.e. with polygonal obstacles) we introduce models that, somehow, represent their one-dimensional analog. We begin considering one of the simplest chaotic model that generates deterministic diffusion: a $`1`$-$`d`$ discrete-time dynamical system on the real axis

$$x^{t+1}=[x^t]+F(x^t[x^t]),$$
(1)
where $`x^t`$ is the variable performing the diffusion (position of a point-like particle) and $`[\mathrm{}]`$ denotes the integer part. $`F(u)`$ is a map defined on the interval $`[0,1]`$ that fulfills the following properties

$`|F^{}(u)|>1`$, so the system has a positive Lyapunov exponent.

$`F(u)`$ must be larger than $`1`$ and smaller than $`0`$ for some values of $`u`$, so there exists a non vanishing probability to escape from each unit cell (a unit cell of real axis is every interval $`C_{\mathrm{}}[\mathrm{},\mathrm{}+1]`$, with $`\mathrm{}𝐙`$).

$`F_r(u)=1F_l(1u)`$, where $`F_l`$ and $`F_r`$ define the map in $`u[0,1/2[`$ and $`u[1/2,1]`$ respectively. This anti-symmetry condition with respect to $`u=1/2`$ is introduced to prevent the presence of a net drift.

The map, $`u^{t+1}=F(u^t)`$ (mod $`1`$) is assumed to be also ergodic. One simple choice of $`F`$ is

$`F(u)=\{\begin{array}{cc}2(1+a)u\text{if }u[0,1/2[\hfill & \\ 2(1+a)(u1)+1\text{if }u[1/2,1]\hfill & \end{array}`$ (4)
where $`a>0`$ is the parameter controlling the instability. This model has been so widely studied, both analytically and numerically KD ; Klag ; Klag1  that can be considered as a reference model for deterministic diffusion. The deterministic chaos which governs the evolution inside each unit cell is responsible of the irregular jumping between cells and the overall dynamics appears diffusive, in the sense that asymptotically the variance of $`x^t`$ scales linearly with time. Following the literature on deterministic diffusion we consider systems with no net drift ($`i.e.`$ $`x^t/t0`$ for $`t\mathrm{}`$). In this paper we will see that there is no need to invoke chaos to induce diffusive behavior in 1-$`d`$ maps, as already shown in DC  in the context of 2-$`d`$ billiards.

The construction of the models begins by dividing the first half of each unit cell $`C_{\mathrm{}}`$ into $`N`$ micro-intervals $`[\mathrm{}+\xi _{n1},\mathrm{}+\xi _n[`$ , $`n=1,\mathrm{},N`$, with $`\xi _0=0<\xi _1<\xi _2<\mathrm{}<\xi _{N1}<\xi _N=1/2`$. In each micro-interval the maps are defined by prescribing a linear function $`F_\mathrm{\Delta }`$ with unit slope. The map in the second half of the unit cell is determined by the anti-symmetry condition iii).

In the first model $`\{\xi _k\}_{k=1}^{N1}`$ is a uniformly distributed random sequence between $`[0,1/2]`$, so the size of the micro-intervals is random and the function is defined as

$`F_\mathrm{\Delta }(u)=\begin{array}{cc}u\xi _n+F(\xi _n)\text{if }u[\xi _{n1},\xi _n[,\hfill & \end{array}`$ (6)
where $`F(\xi _n)`$ is the chaotic map function in Eq. (4) evaluated at $`\xi _n`$.

In the second model, the size of the micro-intervals is uniform, i.e. $`|\xi _n\xi _{n1}|=1/N`$, but the “height” of the function is random. That is,

$`F_\mathrm{\Delta }(u)=\begin{array}{cc}u\stackrel{~}{\xi }_n+\lambda _n\text{if }u[\xi _{n1},\xi _n[,\hfill & \end{array}`$ (8)
where $`\stackrel{~}{\xi }=(\xi _{n1}+\xi _n)/2`$, and the heigth at the mean point, $`\lambda _n`$, is a random variable distributed uniformly in the interval $`[\delta ,1+\delta ]`$, with $`\delta 0`$. Figures 1 and 2 show examples of these two models. Let us note that the random variables $`\xi _n`$ and $`\lambda _n`$ are quenched variables, i.e. in a given realization of a model they depend on the cell $`C_{\mathrm{}}`$ (so one should properly write $`\xi _n^{(\mathrm{})}`$ and $`\lambda _n^{(\mathrm{})}`$) but they do not depend on time.

As Fig. 1 illustrates in model (6) each cell is partitioned in the same number $`N`$ of randomly chosen micro-intervals of mean size $`\mathrm{\Delta }=1/N`$ and in each of them the slope of the map (microscopic slope) is one. For $`\mathrm{\Delta }0`$ (equivalently $`N\mathrm{}`$) we recover the chaotic system (1) but the limit $`F_\mathrm{\Delta }F`$ has to be carefully interpreted. This kind of modification of the original chaotic system is somehow equivalent to replacing circular by polygonal obstacles in the Lorentz system DC , since the steps with unitary slope are the analogous of the flat boundaries of the obstacles. The discontinuities of $`F_\mathrm{\Delta }`$ produce a dispersion of trajectories in a way similar to that of a vertex of polygons that splits a narrow beam of particles hitting on it. However, at variance with the $`2`$-$`d`$ wind-tree model, both randomness and forcing are needed to attain a full diffusive behavior. Note that, thanks to the local preservation of the antisymmetry with respect to the cell center, no net drift is expected and also the Sinai-Golosov Sin ; Gol ; Quench  effect is absent. Also model (8) is non-chaotic and the randomness is introduced in the function but not in the partition of the unit cells.

Since $`F_\mathrm{\Delta }`$ has slope $`1`$ almost everywhere, the chaoticity condition i) is no more satisfied in these models. However, these models can exhibit diffusion provided condition i) is replaced by the following two requirements:

presence of quenched disorder,

presence of a quasi-periodic external perturbation.

In models (6) and (8), requirement i-a) is satisfied by the randomness in the partition of the micro-intervals, and the randomness in the linear function. To fulfill requirement i-b) we introduce a time-dependent perturbation into Eq. (1)

$$x^{t+1}=[x^t]+F_\mathrm{\Delta }(x^t[x^t])+\epsilon \mathrm{cos}(\alpha t)$$
(9)
where $`\epsilon `$ tunes the strength of the external forcing, $`\alpha `$ defines its frequency and the index $`\mathrm{\Delta }`$ indicates a specific quenched disorder realization. We underline again that the linear piecewise map $`F_\mathrm{\Delta }`$ explicitly depends on the cells $`C_{\mathrm{}}`$ due to the disorder, which varies from cell to cell.

Even though by construction the system is not unstable in the Lyapunov sense, it exhibits standard diffusion. In particular, numerical simulations show a linear growth of the mean square displacement. It is important to note that, for model (6) a uniform partition of cell domains plus the quasi-periodic perturbation is not sufficient to produce diffusive dynamics. Even in the presence of quenched randomness (i.e. random partitions or random heights) alone, the system spontaneously settles onto periodic or quasi-periodic trajectories and does not diffuse. The periodic and quasi-periodic behavior is broken down, only when a time dependent perturbation is explicitly superimposed onto the quenched randomness. Hence, we can say that, in our models, the time plays somehow the role of the missing spatial dimension with respect to the $`2`$-$`d`$ models of ref. DC . We will explore this idea in more detail in Sec. V.

## III Diffusion without chaos, numerical results.

In this section, we characterize the diffusion properties of model (6), while in the next, we discuss diffusion on model (8).

The diffusion coefficient,

$$D=\underset{t\mathrm{}}{lim}\frac{1}{2t}(x^tx^0)^2$$
(10)
has been numerically computed from the slope of the linear asymptotic behavior of the mean square displacement. The average has been performed over a large number of trajectories starting from random initial conditions uniformly distributed in cell $`[0,1]`$ and the results are shown in Fig. 3 for three different perturbation amplitudes $`\epsilon `$. In all performed simultations we present in this section we set $`\alpha =1`$; the influence of $`\alpha `$ on diffusive behaviour is discussed in the following section.

We see that diffusion is absent for $`\epsilon `$ small enough, thus, a transition between a diffusive to non-diffusive regime is expected upon decreasing $`\epsilon `$.

When the limit in expression (10) exists, coefficient $`D`$ can be expressed in terms of the velocity correlation function as in Green-Kubo formula Dorfman ; Klag

$$D=\frac{1}{2}(v^0)^2+\underset{k=1}{\overset{\mathrm{}}{}}v^0v^k$$
(11)
where $`x^{i+1}x^i=v^i`$ can be formally regarded as a discrete-time velocity. The simplest approximation for $`D`$ amounts to assuming that velocity correlations decay so fast that only the first term of the above series can be retained. In the pure chaotic system (1), this condition is often verified, and with the further assumption of uniform distribution for the quantity $`x^t[x^t]`$, one obtains the approximated expression

$$D(a)=\frac{a}{2(1+a)}$$
(12)
that works rather well in the parameter region we explored Klag . In the non-chaotic system (9) this approximation is expected to hold only in the large driving regime limit, where the stair-wise structure of $`F_\mathrm{\Delta }`$ is hidden by the perturbation effects.

The role of the external forcing on the activation of the diffusion process in the non chaotic system has been studied by looking at the behavior of the diffusion coefficient $`D`$ upon changing the perturbation amplitude $`\epsilon `$ and the “discretization level” defined by the number of the micro-intervals $`N`$ in each cell.

Results of Fig. 4 show that $`D`$ is significantly different from zero only for values $`\epsilon >\epsilon _c`$, for which the system performs normal diffusion. Below this threshold the values of $`D`$, when non zero, are so small that we cannot properly speak of diffusion. When $`\epsilon `$ increases well above $`\epsilon _c`$, $`D`$ exhibits a saturation toward the value predicted by formula (12) (horizontal line) for the chaotic system defined through (1,4). The onset of a diffusive regime as a threshold-like phenomenon, with respect to the external perturbation amplitude, is expected. It can be explained by noticing that, due to the staircase nature of the system, the perturbation has to exceed the typical discontinuity of $`F_\mathrm{\Delta }`$, to activate the local instability which is the first step toward the diffusion. Data collapsing in Figs 4 confirms this view, because it is achieved by plotting $`D`$ versus the scaling variable $`\epsilon N`$. This means that, for $`N`$ intervals on each cell, the typical discontinuity in the staircase-map is $`O(1/N)`$, then the $`\epsilon `$-threshold is $`\epsilon _c1/N`$. This behavior is robust and does not depend on the precise shape of the forcing, as we have verified by considering other kind of external perturbations. However, as will be discussed in the next section, there is some connection among the diffusive properties of the system, the periodic orbits, the parameter $`N`$, and the value of the perturbation frequency $`\alpha `$.

The activation of the instability performed by the perturbation can be studied by considering a well localized ensemble of initial conditions. We choose uniform initial distribution of $`𝒩`$ walkers with $`\sigma _0^210^{12},10^8,10^6`$ and monitored the distribution spreading under the dynamics by following the evolution of its variance

$$\sigma ^2(t)=\frac{1}{𝒩}\underset{i=1}{\overset{𝒩}{}}(x_i^tx_i^t)^2$$
(13)
Note that, in contrast with the average in Eq. (10) now the average is over an ensemble of initial conditions, uniformly distributed on a small size interval. The spreading of the initial cloud of points is better understood when we make a comparison between the chaotic (1) and non chaotic system (9). In the former, the instability makes the variance growing exponentially at the rate $`\mathrm{ln}[2(1+a)]`$ which is its Lyapunov exponent. This exponential separation lasts until the typical linear behavior of the diffusion takes place (see Fig. 5). In the non-chaotic system, instead, the behavior of $`\sigma ^2(t)`$ strongly depends on its initial value $`\sigma _0^2`$ and we have to consider two sub-cases, whether $`\sigma _01/N`$ or $`\sigma _01/N`$. When $`\sigma _01/N`$, basically one observes the same behavior of the continuous limit i.e. $`\sigma ^2(t)\sigma _0^2\mathrm{exp}(2\lambda t)`$ for times $`t<t_{}\mathrm{ln}(1/\sigma _0)`$ and $`\sigma ^2(t)t`$ for $`t>t_{}`$. Instead, if $`\sigma _01/N`$ we can distinguish three different regimes: a) $`\sigma (t)`$ remains constant for a certain time span $`[0,t_1]`$, b) the instability starts being active and the $`\sigma `$ grows exponentially as in the chaotic case, $`\sigma (t)(1/N)\mathrm{exp}[\lambda (tt_1)]`$, c) the system eventually reaches the regime of standard diffusion where $`\sigma ^2`$ behaves linearly. The crossover time $`t_1`$ between the regimes a) and b) depends on the size of the discontinuities of $`F_\mathrm{\Delta }`$ and the specific realization of the randomness and, of course, it decreases when either $`N`$ or $`\sigma _0`$ increases. We stress that the above scenario is due to the presence of finite size instabilities that in a non chaotic Lorentz gas would correspond to the defocusing mechanism of a beam by the vertices of polygonal obstacles.

The onset of standard diffusive behavior in model (9) depends on large-time features of the deterministic dynamics. Indeed the presence of periodic trajectories, standing or running note , have a strong influence on the diffusion, because the formers tend to suppress diffusion while the latter induce a “trivial” ballistic behavior. In Ref. DC , it has been conjectured that the main difference between a fully chaotic deterministic system and a non chaotic system exhibiting diffusion, as a consequence of quenched randomness, is in the different effects the periodic trajectories have on the systems. An escape rate method Dorfman ; Gaspard  can be applied to analyze these effects. This method consists on starting with an ensemble of $`W(0)`$ trajectories initially localized in the region $`_L`$ between $`L/2`$ and $`L/2`$ and computing the number $`W(t)`$ of them which never left $`_L`$ until the time $`t`$. For a pure random walk process one would expect an exponential decay

$$W(t)\mathrm{exp}(\gamma t)$$
(14)
where the rate $`\gamma `$ is proportional to the diffusion coefficient $`D`$ through the relation $`\gamma D/L^2`$. In this sense the method is an alternative way to estimate the diffusion coefficient from measurement of $`\gamma `$. However possible deviations from (14) generally indicate a non diffusive behavior related to the presence of periodic orbits. As clearly seen in fig 6, on all the range of sizes $`L`$ we have considered, we do not observe any substantial deviation from the exponential decay that can be associated to periodic orbit effects.

One can wonder whether the diffusive behavior could be merely due to a possible non zero spatial entropy per unit length of the quenched randomness. To unravel the role of static randomness as possible source of entropy and thus diffusion, following Ref. DC , we have considered a system were the same configuration of the disorder is repeated every $`M`$ cells (i.e. $`\xi _n^{(\mathrm{})}=\xi _n^{(\mathrm{}+M)}`$), so the entropy for unit length is surely zero. This point is crucial, because, someone can argue that a deterministic infinite system with spatial randomness can be interpreted as an effective stochastic system, but this is a “matter of taste”. Anyway, our system with a spatially periodic randonness is deterministic from any point of view.

Looking at the diffusion of an ensemble of walkers we obtain that:

there is weak average drift $`V`$, that vanishes approximately as $`V1/M`$ when $`M`$ goes to infinity.

there is still a diffusive behavior i.e. $`(x^tVt)^22Dt`$, with a coefficient $`D`$ very close to the value obtained by formula (12).

Let us note that this behavior is obtained by considering single realizations of the randomness: the value of the drift changes with the realization, in such a way that an average on the disorder would give a null drift. The fact that $`V1/M`$ can be easily interpreted as a self-averaging property, for large $`M`$. On the other hand if we write $`(x^t)^2=(Vt)^2+2Dt`$, we can define a crossover time $`\tau _cDM^2`$, from pure diffusive to a ballistic regime. When $`M`$ is finite but large enough the crossover time becomes very large. Therefore, we can safely conclude that, even with a vanishing (spatial) entropy density of randomness, at least for sufficiently large $`M`$, the system is able to effectively perform standard diffusion for a very long time. However we stress that looking at $`(x^tVt)^2`$ we do not observe any crossover.

## IV Space-time billiards

As mentioned previously, in the one-dimensional maps proposed in this paper, time plays the role of the second dimension in the two-dimensional billiard models (e.g., Lorentz and Ehrenfest models). Here we show this explicitly by describing how the one-dimensional time-dependent maps define billiards in space-time. For convenience, we carry out the construction for model (8) where the partition of the unit cells is uniform. However, the same construction can be done for model (6).

Let $`[\xi _{n1},\xi _n]`$ be a micro-interval of a given unit cell. Then, for any initial condition $`x_0[\xi _{n1},\xi _n]`$, the velocity of the particle at time $`t`$ is

$$h_n^t=x^{t+1}x^t=\lambda _n+\epsilon \mathrm{cos}(\alpha t),$$
(15)
where $`h_n^t`$ is the height (including the time-dependent perturbation) of the piece-wise linear segment defining the map according to Eq. (8). As shown in Fig. 7, in a space-time diagram the “world-line” of this orbit is the straight line joining the points $`(x^t,t)=(x^{t+1},t+1)`$ with slope $`h_n^t`$. Accordingly, as shown in Fig. 8, the space-time can be decomposed into discrete boxes $`[\xi _{n1},\xi _n]\times [t,t+1]`$, and the world-line of successive iterates of the map appears as a zig-zag trajectory joining different boxes. The rule of the billiard is that a particle entering box $`n`$ will leave this box with speed $`h_n^t`$. That is, in going from box $`m`$ to box $`n`$ a particle suffers a space-time scattering event $`\delta v=h_n^{t+1}h_m^t`$, where $`\delta v=v^{t+1}v^t`$ is the change in the particle velocity. In the space-time diagram the scattering rule is represented by drawing at each box a line segment with slope given by $`h_n^t`$ that according to (15) has a random component $`\lambda _n`$, and a time-periodic component.

According to this, the time-dependent, zero-Lyapunov exponent (slope one) map is equivalent to a space-time billiard with scattering rule $`h_n^t`$. Note that these models are simpler than the usual 2-$`d`$ billiards (e.g. Lorentz gas) because the scattering rule is independent of the angle of incidence. Although this is just a reformulation of the same model, the billiard construction provides useful insights into the problem and allows the connection with the two-dimensional billiard problems. When the map is time-independent ($`\epsilon =0`$) the scattering rule depends only on $`n`$ ($`h_n^p=h_n^q`$ for any integers $`p`$ and $`q`$) and as shown in Fig. 8-(a) the billiard structure is simply a copy of the billiard at $`t=0`$. In this case, despite the spatial randomness, the lack of a time dependence leads to a highly structured space-time billiard in which periodic orbits are ubiquitous. As a consequence, in this case there is no diffusion, a result already mentioned before. On the other hand, as Figure 8-(b) shows, when $`\epsilon 0`$ the billiard exhibits space-time disorder which leads to the possibility of diffusive behavior.

As discussed in Ref. DC  for polygonal two-dimensional billiards, periodic orbits play an important role. The existence of periodic orbits in the space-time billiard is intimately connected with the frequency $`\alpha `$ of the time periodic perturbation. If $`\{x^t\}`$ is a periodic orbit of period $`M`$ ($`x^{t+M}=x^t`$ for any $`t`$), the billiard must be $`q`$-periodic, that is $`h_n^{t+q}=h_n^t`$ for any $`t`$, where $`M=mq`$ and $`m`$ is a positive integer. According to Eq. (15) this will be the case provided

$$\alpha =2\pi \nu ,$$
(16)
with $`\nu =p/q`$. Figure 9 show two examples of periodic orbits with $`\epsilon =0.2`$. In panel (a), $`\alpha =2\pi (3/4)`$, and the period of the billiard is the same as the period of the orbit, $`q=M=4`$. In panel (b), $`\alpha =2\pi (2/5)`$, and the period of the orbit is twice period of the billiard $`M=2q=10`$. Thus, for rational $`\nu `$ the billiard shows a periodic structure in time, and periodic orbits are likely to exist leading to the suppression of diffusion. On the other hand, when $`\nu `$ is an irrational number the billiard structure never repeats exactly in time and strictly speaking there are no periodic orbits. Motivated by this, it is possible to related the degree of temporal disorder of the billiard to the degree of irrationality of $`\nu `$ as determined, for example, by the continued fraction approximation of $`\nu `$. In this regard, billiards with irrational frequencies $`\nu `$ which are easily approximated by rationals are “less disordered” than billiards with hard to approximate frequencies. To explore this ideas Fig. 10 shows the time evolution of the variance for different values of the perturbation frequency in the same random realization of model (8) with $`N=40`$ and $`\epsilon =0.4`$. For the perturbation frequencies, $`\alpha =(p/q)2\pi `$, we took the continued fraction expansion of the inverse golden mean

$$\nu =1/2,\mathrm{\hspace{0.17em}2}/3,\mathrm{\hspace{0.17em}3}/5,\mathrm{}F_{n1}/F_n\mathrm{}\left(\sqrt{5}1\right)/2=\gamma ^1,$$
(17)
where $`F_0=F_1=1`$, and $`F_n=F_{n1}+F_{n2}`$. The chose of the golden mean is motivated by the fact that this number is the hardest to approximate with rational. In the model the disorder configuration is repeated every $`M`$ cells, and for small $`M`$ the recurrence time of an orbit (i.e. the time to return to the same microinterval, module $`M`$) decreases. Since we are interested in the role of periodic orbits, here we consider the limit case $`M=1`$. That is, in the simulations reported in Fig. 10, the disorder configuration in the unit interval $`[0,1]`$ was copied to all the cells $`[\mathrm{},\mathrm{}+1]`$. As expected the figure shows the suppression of diffusion for rational frequencies. However, more interesting is the fact that there is a tendency for the onset of the suppression of diffusion to increase as $`\nu `$ approaches $`1/\gamma `$. This is related to the fact that for rationals with lower denominators the time period of the billiard is small, and the probability that the recurrence time of an orbit coincides with a multiple of the period of the billiard is high. The recurrence time is also related to the value of $`N`$. When $`N`$ is large, the probability for an initial condition to return to the same micro-interval (module $`M`$) is small, and periodic orbits are more scarce. This explains why for large enough $`N`$ and $`M`$, diffusion might not be suppressed even for rational frequencies with large denominators. This suggests that, for large enough $`N`$ and $`M`$, deviations from diffusion are unlikely to be observed numerically in the case of rational frequencies with large denominators.

## V Conclusions and remarks

From the above results we can conclude that two basic ingredients are relevant for diffusion, 1) a finite size instability mechanism: in the chaotic system this is given by the Lyapunov instability, in the ”stepwise” system this effect stems from the jumps 2) a mechanism to suppress periodic orbits and therefore to allow for a diffusion at large scales. In the presence of “strong chaos” (i.e. when all the periodic orbits are unstable) point 2) is automatically guaranteed thus the periodic orbits formalism Artuso ; CEG95 ; Dorfman  can be fully applied to compute the diffusion coefficient. In non chaotic systems, however, like the present one, the mechanism 2) is the result of a combined effect of quenched randomness and quasi-periodic forcing.

In summary, we have introduced and studied models exhibiting diffusion in the absence of any source of chaotic behavior. The models represent the $`1`$-$`d`$ analogue of non chaotic Lorentz gas (i.e. with polygonal obstacles) discussed by other authors in connection with the debate around experimental evidences for a distinction between chaotic and stochastic diffusion. Our results indicate that when the chaotic instability condition (positive Lyapunov exponent) is replaced by the presence of finite size instability (non positive Lyapunov exponent) we need both an external quasi-periodic perturbation and disorder for preventing the system falling into either a periodic/quasi-periodic evolution or a drift dominated behavior.

We believe that this kind of behaviour is rather general, in the sense that chaos does not seem to be a necessary condition for the validity of some statistical features. This is in agreement with the recent results of Ref. LRB , where the applicability of the generalized Gallavotti-Cohen fluctuation formula GC  has been proven for non chaotic systems too.

## VI Acknowledgments

This work has been partially supported by MIUR (Cofin. Fisica Statistica di Sistemi Classici e Quantistici). D dCN was supported by the Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. We thank Massimo Cencini and Stefano Ruffo for useful suggestions and discussions. D dCN gratefully acknowledges the hospitality of the Department of Physics of the University of Rome “La Sapienza” during the elaboration of this work.


format: text

text:
# Non-Drude Optical Conductivity of (III,Mn)V Ferromagnetic Semiconductors

## I  Introduction

Recent studies of (III,Mn)V diluted magnetic semiconductors (DMSs), motivated by the discovery of carrier-mediated ferromagnetism Ohno et al. (1992) in these materials, have uncovered a wealth of interesting phenomenology involving an interplay between collective electronic effects, broken symmetries, interactions, and disorder. The richness of the transport and optical properties of (III,Mn)V DMSs Ohno et al. (1992, 1996) arises from their strong valence-band spin-orbit coupling and from the sensitivity of their magnetic state to growth and annealing conditions, doping, and external fields Edmonds et al. ; Potashnik et al. (a). In situ tunability of material properties by applying gate voltages Ohno et al. (2000), by illuminationOiwa et al. (2002), or by applying external magnetic fields, suggests many opportunities for exploring new physics or building new devices using these materials. There is considerable evidence that the properties of these systems can be understood using a simplified effective model Konig et al. (2002) in which the low energy degrees of freedom are $`S=5/2`$ Mn2+ local moments, and holes in the host-semiconductor valence band. The Hamiltonian of this effective model should in general include terms that represent scalar and exchange interactions between the Mn ions and the valence band holes, and Coulomb interactions among valence band holes and between the holes and the Mn ions. The simplest theoretical approach to these models, one that appears to be relatively successful in explaining many bulk properties of strongly metallic DMS ferromagnets, either ignores disorder or treats it perturbatively and approximates exchange interactions using a mean-field theory Dietl et al. (2000, 2001); Abolfath et al. (2001); Konig et al. (2002) 111 A large database with theory predictions based on this approximation for different host materials, and different carrier and Mn densities can be accessed at  http://unix12.fzu/ms.. We have recently presented a theory of the intraband optical conductivity of (III,Mn)V ferromagnets that uses such an approach Sinova et al. (2002). In the present paper we report a finite-size numerical study of the optical conductivity of (Ga,Mn)As in which disorder is treated exactly, enabling us to address weakly metallic or insulating systems and to understand some of the qualitative strong-scattering effects that are omitted in the simpler theory.

Our interest in this property of DMS ferromagnets is motivated by a series of recent optical absorption experiments Nagai et al. (2001); Katsumoto et al. (2001); Hirakawa ; Singley et al. . The experiments, performed in thin film geometries in both metallic and insulating samples (as judged by the temperature dependence of the resistivity below $`T_c`$), exhibit several common features: (a) non-Drude behavior in which the conductivity increases with increasing frequency in the interval between $`0`$ meV and $`220`$ meV, (b) a broad absorption peak near $`220260`$ meV that becomes stronger as the samples are cooled, and (c) a broad featureless absorption between the peak energy and the effective band gap energy which tends to increase at the higher frequencies. The $`220`$ meV peak can be attributed to inter-valence-band transitions Sinova et al. (2002), to transitions between the semiconductor valence band states and a Mn induced impurity band, or to a combination of these contributions. Dynamical mean-field-theory studies Hwang et al. (2002), for a single-band model that neglect the spin-orbit coupling and the heavy-light degeneracy of a III-V semiconductor valence band, demonstrate that non-Drude impurity-band related peaks in the frequency-dependent conductivity occur in DMS ferromagnet models when the strength of exchange interaction coupling is comparable to the valence band width. The conductivities predicted by this model are, however, inconsistent with experiment in temperature trend and appear not to be in the same regime as experimental samples. On the other hand we Sinova et al. (2002) demonstrated in earlier work that peaks in the conductivity can also occur in the weak-coupling regime, provided that the realistic multi-band character of the valence bands is acknowledged. In this approach the $`220`$ meV peak is due to heavy-hole to light-hole inter-valence-band transitions and the gradual increase in absorption observed in the experiments Katsumoto et al. (2001) at frequencies above $`500`$ meV is ascribed to transitions to the bands split-off by strong spin-orbit interactions. Although this theory could account for many overall features of the measured optical conductivity in metallic samples, one prominent feature was at odds with experimental data, namely the relative magnitudes of the $`\omega 0`$ conductivity and the $`220`$ meV conductivity peak. The numerical calculations reported on in this paper suggest that this discrepancy is due to multiple-scattering effects that suppress the low-frequency conductivity.

We restrict our attention in this paper to the $`T=0`$ limit, which allows us to neglect scattering of valence band quasiparticles off thermal fluctuations in the Mn ion spin orientations. We assume a simple ground state in which the Mn spins have parallel spin orientations 222The combined effect of disorder and spin-orbit interactions implies that the classical ground state of a (III,Mn)V ferromagnet cannot have perfectly parallel Mn spin orientations. Disorder can also lead to frustration in the carrier-mediated Mn spin interactions, and to complex ground states in which spins are not aligned. Both effects become less important for more strongly metallic (III,Mn)V ferromagnets. See Ref. Zarand and Janko, 2002; Schliemann and MacDonald, 2002.. Sources of disorder that are known to have some importance in (III,Mn)V ferromagnets include randomness in the placement of Mn ions on cation sites of the host lattice, random placement of Mn ions at interstitial sites, and randomly located antisite defects in which group V elements are placed at group III element lattice sites. In the low-energy effective model, Mn ions at cation sites have charge $`Q=e`$ while interstitial Mn and antisite defects have $`Q=2e`$. Coulomb and exchange interactions with the Mn ions have comparable importance Sinova et al. (2002) for the scattering of valence band holes and must both be included in the theory. In the approach of Ref. Sinova et al., 2002 these sources of disorder scattering were treated perturbatively using a Born approximation, evaluating the valence band quasiparticle lifetimes using Fermi’s Golden rule Jungwirth et al. . In this paper we find the $`T=0`$ quasiparticle states for finite-size systems by solving the Schrödinger equation. Unlike the Born approximation theory, this approach does not assume metallic behavior and remains valid at low Mn concentrations when Yang and MacDonald  Mn acceptor impurity bands emerge. We find that corrections to the Born approximation theory are important at a quantitative level even for strongly metallic high Mn concentration samples, and that they explain the property that the low-frequency and dc conductivities are suppressed compared to their values in the Born approximation theory. We also find that the f-sum rule, obtained by integrating the conductivity over frequency, differs from its Born approximation value by less than 10% for typical metallic carrier densities.

This paper is organized as follows. Section II defines the model Hamiltonian we use to describe the valence band system in a (III,Mn)V ferromagnet which adds exchange and Coulomb interactions to the six band $`𝐤𝐩`$ envelope function Hamiltonian of the host (III,V) semiconductor. Section III describes technical issues associated with the evaluation of Kubo linear-response-theory formula for the finite-size system optical conductivity. In Sec. IV we discuss the numerical results we have obtained for different Mn and carrier density regimes, which are in qualitative agreement with experiment and discuss the degree to which aspects of the extreme impurity-band Hwang et al. (2002) and Golden-rule Sinova et al. (2002); Jungwirth et al.  limits are reflected in the exact results. We also emphasize the feasibility of using the f-sum rule to extract reasonably accurate measurements of the carrier concentration Sinova et al. (2002), circumventing the inherent uncertainties in Hall measurements in systems with large anomalous Hall effect coefficients.

## II Model Hamiltonian

We approximate the host semiconductor valence band by its six-band $`𝐤𝐩`$ Kohn-Luttinger model. This simplification takes advantage of the small number of holes per atom in these ferromagnets. The band electrons interact via Coulomb and exchange interactions with randomly distributed Mn spins and with randomly distributed antisite defects via Coulomb interactions alone; we do not explicitly account for the possible role of Mn interstitials Yu et al. (2002). The requirement of overall charge neutrality implies that the density of valence band holes ($`p`$) is related to the density of Mn ($`n_{Mn}`$) and the density of antisites ($`n_{AS}`$) by $`p=n_{Mn}2n_{AS}`$. In MBE grown (III,Mn)V ferromagnets, it is often a challenge to accurately determine the density of holes, and even the density of Mn ions that have been incorporated into the ferromagnetic state, although progress is being made. It appears at present that any annealing procedure Potashnik et al. (a, b) that leads to higher carrier densities, inevitably leads to higher ferromagnetic transition temperatures, and to much higher conductivities Edmonds et al. . For Ga0.95Mn0.05As for example, it is now possible to prepare samples with $`T=0`$ dc conductivities that are close to an order of magnitude larger than in the samples for which ac conductivity measurements have so far been carried out. Although we believe that much of the variabilty in physical properties can be accounted for in terms of variations in the fraction of Mn ions that are incorporated and, more importantly, variations in the carrier density, it appears likely that a part of this variation is still associated with sources of unintended disorder that are not present in our model. We expect that as unintended disorder is reduced, the properties of these ferromagnets will be more completely determined by the two parameters of our model calculations, $`p`$ and $`n_{Mn}`$ 333There is a broad consensus that variation in the density of Mn interstitials is responsible for most of the variation of carrier density and magnetically active Mn density that occurs when these materials are annealed. See for example Ref. Yu et al., 2002..

Below we will evaluate the frequency-dependent conductivity by using the Kubo formula description for the linear response of quasiparticles to an electromagnetic field. Coulomb interactions play an essential role in determining the quasiparticle states and interactions among the holes cannot be ignored. In this paper we use a Hartree mean-field approximation for the quasiparticles, which is consistent with the standard bubble diagram approximation we use for the conductivity.

The single particle part of the quasiparticle Hamiltonian may be written as the sum of host crystal band and disorder pieces $`\widehat{H}^0=\widehat{H}^h+\widehat{H}^D`$: The disorder Hamiltonian has three contributions $`\widehat{H}^D=\widehat{H}^{K.ex.}+\widehat{H}^{Mnh}+\widehat{H}^{Ash}`$. $`\widehat{H}^{Mnh}`$ describes the attractive Coulomb interaction between the ionized $`Mn^{2+}`$ acceptor and a valence band hole, $`\widehat{H}^{K.ex.}`$ describes the local kinetic-exchange interactions between valence band carriers and the Mn moments, while $`\widehat{H}^{Ash}`$ describes the repulsive Coulomb interaction between a hole carrier and ionized anti-site defects. The single-particle envelope function model Hamiltonian then reads:

$$\widehat{H}^0=\widehat{H}^h+\underset{I=1}{\overset{N_{Mn}}{}}\stackrel{}{S}_I\widehat{\stackrel{}{s}}J(\stackrel{}{r}\stackrel{}{R_I})+\underset{I=1}{\overset{N_{Mn}}{}}(\frac{e^2}{ϵ|\stackrel{}{r}\stackrel{}{R_I}|}V_0e^{|\stackrel{}{r}\stackrel{}{R_I}|^2/r_0^2})\widehat{I}+\underset{K=1}{\overset{N_{As}}{}}\frac{2e^2}{ϵ|\stackrel{}{r}\stackrel{}{R_K}|}\widehat{I}$$
(1)
where $`J(\stackrel{}{r})=(J_{pd})/((2\pi a_0^2)^{3/2})e^{r^2/2a_0^2}`$, and $`\widehat{\stackrel{}{s}}=(\widehat{s_x},\widehat{s_y},\widehat{s_z})`$ where $`\widehat{s}_{x,y,z}`$ are the $`6\times 6`$ matrices which describe hole spins in the representation of the six band envelope-function Hamiltonian, and $`\widehat{I}`$ is a $`6\times 6`$ unit matrix. In this Hamiltonian, $`I`$ labels Mn sites, $`K`$ labels As sites, and $`\stackrel{}{S}_I`$ stands for a Mn spin with quantum number $`S=5/2`$. The Mn and As positions are denoted by $`\stackrel{}{R_I}`$ and $`\stackrel{}{R_K}`$. The term proportional to $`V_0`$ is a central cell correction Bhattacharjee and a la Guillaume (2000) to the Coulomb attraction between holes and Mn acceptors that is known to be required to recover experimental values for the isolated acceptor energy. The quantities $`N_{Mn}`$ and $`N_{As}`$ are the number of Mn and As atoms, which we distributed randomly inside cubic simulation cells of side $`L`$. Charge neutrality in our simulation cell means that the number of quasiparticle levels occupied by holes is $`N_h=N_{Mn}2N_{As}`$. The values chosen for the phenomenological parameters that appear in this equation are $`J_{pd}=0.05`$eV nm3, $`a_0=0.3`$ nm, $`ϵ=10.9`$, $`V_0=2.5eV`$ and $`r_0=0.259`$ nm.

The host band part of the Hamiltonian is described via the six band Kohn-Luttinger model, with wavevectors $`\stackrel{}{k}`$ in the envelope function replaced by $`k_a=\frac{1}{i}_a`$ for $`a=x,y,z`$ to account for the inhomogeneity induced by disorder. Choosing the angular momentum quantization direction to be along $`z`$-axis, and ordering the $`j=3/2`$ and $`j=1/2`$ basis functions according to the list ($`3/2,1/2,1/2,3/2;1/2,1/2`$), the Luttinger Hamiltonian $`H^L`$ has the form Abolfath et al. (2001):

$$H^L=\left(\begin{array}{cccccc}_{hh}& c& b& 0& \frac{b}{\sqrt{2}}& c\sqrt{2}\\ c^{}& _{lh}& 0& b& \frac{b^{}\sqrt{3}}{\sqrt{2}}& d\\ b^{}& 0& _{lh}& c& d& \frac{b\sqrt{3}}{\sqrt{2}}\\ 0& b^{}& c^{}& _{hh}& c^{}\sqrt{2}& \frac{b^{}}{\sqrt{2}}\\ & & & & & \\ \frac{b^{}}{\sqrt{2}}& \frac{b\sqrt{3}}{\sqrt{2}}& d^{}& c\sqrt{2}& _{so}& 0\\ c^{}\sqrt{2}& d^{}& \frac{b^{}\sqrt{3}}{\sqrt{2}}& \frac{b}{\sqrt{2}}& 0& _{so}\end{array}\right)$$
(2)
In the matrix (2) we have highlighted the $`j=3/2`$ sector. The Kohn-Luttinger eigenenergies are measured down from the top of the valence band, i.e. they are hole energies. For completeness we list the expressions which define the quantities that appear in $`H^L`$:

$`_{hh}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}[(\gamma _1+\gamma _2)(k_x^2+k_y^2)+(\gamma _12\gamma _2)k_z^2`$
$`_{lh}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}[(\gamma _1\gamma _2)(k_x^2+k_y^2)+(\gamma _1+2\gamma _2)k_z^2`$
$`_{so}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\gamma _1(k_x^2+k_y^2+k_z^2)+\mathrm{\Delta }_{so}`$
$`b`$ $`=`$ $`{\displaystyle \frac{\sqrt{3}\mathrm{}^2}{m}}\gamma _3k_z(k_xik_y)`$
$`c`$ $`=`$ $`{\displaystyle \frac{\sqrt{3}\mathrm{}^2}{2m}}\left[\gamma _2(k_x^2k_y^2)2i\gamma _3k_xk_y\right]`$
$`d`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}\mathrm{}^2}{2m}}\gamma _2\left[2k_z^2(k_x^2+k_y^2)\right].`$ (3)
with $`\gamma _1=6.98`$, $`\gamma _2=2.06`$ and $`\gamma _3=2.93`$.

The Hartree potential due to interactions among holes must also be included in the single-particle Hamiltonian for quasiparticle states since it captures the screening of Coulomb interactions with impurities due to the build up or depletion of charge. The total one-particle Hamiltonian is $`\widehat{H}=\widehat{H}^0+\widehat{V}^H`$, where

$`V_{\stackrel{}{k}j,\stackrel{}{k^{}}j^{}}^H=\delta _{j,j^{}}{\displaystyle \underset{j\mathrm{"}}{}}{\displaystyle \underset{\stackrel{}{p}}{}}{\displaystyle \frac{4\pi e^2}{ϵ|\stackrel{}{k}\stackrel{}{k^{}}|^2}}\rho _{\stackrel{}{p}j\mathrm{"},\stackrel{}{k^{}}+\stackrel{}{p}+\stackrel{}{k}j\mathrm{"}},`$ (4)
and the density matrix

$`\rho _{\stackrel{}{k}j,\stackrel{}{k^{}}j^{}}`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}f_\alpha <\mathrm{\Psi }_{\stackrel{}{k}j}|\alpha ><\alpha |\mathrm{\Psi }_{\stackrel{}{k^{}}j^{}}>`$ (5)
$`=`$ $`{\displaystyle \underset{\alpha }{}}f_\alpha c_{\stackrel{}{k}j}^{(\alpha )}c_{\stackrel{}{k^{}}j^{}}^{(\alpha )},`$
with $`f_\alpha =1(0)`$ for an occupied (unoccupied) state. Here $`\mathrm{\Psi }_{\stackrel{}{k}j}(\widehat{r})=e^{i\widehat{k}\widehat{r}}u_{\stackrel{}{k}j}(\stackrel{}{r})`$, where $`u_{\stackrel{}{k}j}(\stackrel{}{r})`$ are the Bloch functions. The Hartree potential must be determined by solving self-consistently for the density matrix. We have diagonalized this Hamiltonian using a plane-wave representation for the envelope functions of $`\widehat{H}^0`$ and applying periodic boundary conditions in a cube of size $`L`$, which limits wavevectors to the discrete set $`\stackrel{}{k}=(n_x,n_y,n_z)(2\pi /L)`$, where $`n_x,n_y,n_z`$ are integers. For different system sizes we choose different maximum values of $`n_{x,y,z}`$, denoted by $`n_m`$, so that the maximum wavevector $`k_{max}`$ is held fixed. For the calculations presented here $`k_{max}=\pi `$ nm-1.

## III Kubo formula for the Optical conductivity

We evaluate the real part of the frequency dependent conductivity from the Kubo formula expression which relates it to the quasiparticle eigenvectors $`|\alpha >`$ and eigenvalues $`E_\alpha `$.

$`\sigma _{ab}(\omega )`$ $`=`$ $`{\displaystyle \frac{\pi e^2}{m^2L^3\omega }}{\displaystyle \underset{\alpha ,\beta }{}}(f_\alpha f_\beta )\alpha |\widehat{p}_a|\beta \beta |\widehat{p}_b|\alpha `$
$`\times \delta (\mathrm{}\omega E_\beta +E_\alpha ).`$
Matrix elements of the momentum operator, $`\widehat{p}_a=\frac{\mathrm{}}{i}_a`$ with $`a=x,y,z`$, are given by

$`<\alpha |\widehat{p}_a|\beta >`$ $`=`$ $`{\displaystyle \underset{j,j^{}}{}}{\displaystyle \underset{\stackrel{}{k},\stackrel{}{k^{}}}{}}c_{\stackrel{}{k}j}^{(\alpha )}c_{\stackrel{}{k^{}}j^{}}^{(\beta )}<\mathrm{\Psi }_{\stackrel{}{k}j}|\widehat{p}_a|\mathrm{\Psi }_{\stackrel{}{k^{}}j^{}}>`$ (7)
$`=`$ $`{\displaystyle \underset{\stackrel{}{k},j}{}}c_{\stackrel{}{k}j}^{(\alpha )}c_{\stackrel{}{k}j^{}}^{(\beta )}{\displaystyle \frac{m}{\mathrm{}}}{\displaystyle \frac{H_{j,j^{}}^L}{k_a}}.`$
At zero temperature and positive $`\omega `$ the difference $`(f_\alpha f_\beta )`$ restricts $`\alpha `$ to occupied quasiparticle states and $`\beta `$ to empty quasiparticle states. For finite-size systems the real part of the conductivity is more reliably evaluated by broadening the $`\delta `$ function Imry (1997)

$`C_R^{ab}(\omega )`$ $`=`$ $`\mathrm{Re}[{\displaystyle \frac{\mathrm{}e^2}{m^2L^3}}{\displaystyle \underset{\alpha ,\beta }{}}{\displaystyle \frac{f_\alpha f_\beta }{E_\beta E_\alpha }}<\alpha |\widehat{p}_a|\beta >`$ (8)
$`\times <\beta |\widehat{p}_b|\alpha >{\displaystyle \frac{\gamma }{(\mathrm{}\omega E_\beta +E_\alpha )^2+\gamma ^2}}],`$
with the level broadening, $`\gamma `$ , being of the order of the mean level spacing $`E_F/N_h1/L^3`$. In the thin film limit the real part of the conductivity is proportional to the absorption coefficient in transmission experiments, $`\stackrel{~}{\alpha }(\omega )=\frac{8\pi Re[\sigma (\omega )]}{c(1+n)}`$ with $`n=\sqrt{ϵ_0}`$ being the zero frequency limit index of refraction of the thin film material and $`ϵ_0`$ the dimensionless dielectric constant. This formula for $`\stackrel{~}{\alpha }(\omega )`$ is most accurate in the far infrared regime or in the thin film limit since at higher frequencies multiple scattering from the interfaces must be taken into account. This can be easily done using the standard electromagnetic description of light waves propagating through a medium, but the simple proportionality is lost.

## IV Results

### IV.1 Optical conductivity of parabolic band with disorder

It is useful for pedagogical purposes to first present frequency-dependent conductivities calculated with this approach for the case where the 6-band Luttinger model is replaced by a parabolic band model with a mass equal to the heavy hole mass of GaAs. In this way we can partially separate the effects of strong disorder on the optical conductivity from the complications associated with the heavy, light, and split-off valence bands. In the calculation all Mn spins point along the z-axis so that spin up and down single-particle Hamiltonians decouple. For each spin the Hamiltonian is identical to that of Eq. 1 except that $`\widehat{H}^L`$ is replaced by the usual kinetic term $`\mathrm{}^2^2/2m_{hh}`$ and the $`6\times 6`$ identity matrix is replaced by unity. The following parameters were used in the calculation of the conductivity: p=0.33 nm-3, $`n_{Mn}=1`$ nm-3 ($`x4.5`$%), $`L=8`$ nm, and $`n_m=4`$. The results presented below were obtained by averaging over $`N_D=10`$ realizations of the disorder potential.

Fig. 1 showss conductivities and absorption coefficients we have evaluated for minority and majority carriers. Even for this simple band model the conductivity displays clear non-Drude behavior, increasing as a function of frequency at low frequencies. We believe that the non-Drude behavior below 200 meV reflects multiple scattering effects that enhance high-angle scattering and would lead to localization if the disorder were stronger. These effects are explicitly neglected in the Boltzmann transport theory which leads to the Drude formula. The f-sum rule for the conductivity calculated by integrating these numerical results over frequency deviates from the parabolic band model value by less than $`10\%`$, demonstrating that the errors induced by constructing a smooth conductivity curve using broadened $`\delta `$\- functions are under reasonable control. The f-sum rule is more accurately satisfied for larger system sizes, for which we are able to choose smaller values of the broadening parameter $`\gamma `$. The non-Drude behavior we find here does capture one aspect of the experimental results. However the overall conductivity magnitude of this model is much smaller than in experiment. Evidently the mixing of heavy and light hole bands leads to quasiparticles with larger velocities.

### IV.2 Optical conductivity of Luttinger-Kohn Hamiltonian with disorder

Mixing of heavy- and light-hole bands and its interplay with magnetic order and scattering disorder can be investigated by solving the full Hamiltonian $`H=\widehat{H}^0+\widehat{V}^H`$. This approach treats the band structure realistically and the randomly located charged impurities without approximation.

Fig. 2 shows the frequency-dependent conductivity we have evaluated for $`p=0.33`$ nm-3, $`n_{Mn}=1`$ nm-3, $`L=6`$ nm, $`N_D=4`$ and $`n_m=3`$. For $`\mathrm{}\omega <200`$ meV we observe that the conductivity increases with frequency, unlike the Drude formula case. The increase is actually sharper than in the parabolic band model, a property that we believe reflects the additive contributions of non-Drude intra-band absorption and broadened inter-valence-band absorption. There are significant differences between these results and those of the Born approximation (see Fig.1 of Sinova  et al. Sinova et al. (2002)) calculated for the same model parameters. In our calculation the maximum value of the absorption is about 90000 cm-1 while in the calculation of Sinova  et al. it is about 110000 cm-1. We believe that this reduction is due to the enhanced backscattering effects captured by the exact calculation and also due to a transfer of spectral weight out of the intra-band contribution. The peak near $`\mathrm{}\omega 0.9`$ eV in these calculations is due to heavy-hole to split-off band transitions and is more pronounced here than in the Born approximation calculations. The presence of this peak at rather large frequencies emphasizes one weakness of the present calculation, namely that it neglects valence to conduction band transitions. It will be interesting to see whether or not future experiments confirm our prediction of a peak in the absorption due to these transitions that lies at high frequencies, but still clearly below the onset of valence-to-conduction band transitions.

The results shown in Fig. 2 are for parameters that are typical for high $`T_c`$ (Ga,Mn)As ferromagnets. The dc limits of the conductivities that we calculate are substantially higher than those of the samples for which $`\sigma (\omega )`$ values have been published at present, but comparable to or smaller than those measured for the highest conductivity (Ga,Mn)As ferromagnets currently available. Fig. 3 displays the absorption curves calculated for a more disordered system 444 For carrier density drops below $`0.2\mathrm{nm}^3`$, disorder localizes the valence band electrons Yang and MacDonald . More strongly localized valence band holes have more heavy-hole character, causing the the optical mass to increase in comparison to values that are appropriate in the metallic limit.. The parameters are $`p=0.2`$ nm-3, $`n_{Mn}=0.88`$ nm-3, $`L=6`$ nm, $`N_D=4`$ and $`n_m=3`$. The conductivity is substantially suppressed at all frequencies, but most strongly at low frequencies. It again has non-Drude behavior at low-frequencies, increasing with increasing frequency over the interval between 0 meV and 220 meV. The broad heavy-to-light absorption peak near 220 meV is still clear but the heavy-hole to split-off band peak is less visible, perhaps explaining the fact that it has not yet been identified experimentally. Hirakawa et al. report that the maximum value for absorption in the more disordered sample they study is about 20000 cm-1, which has the same order of magnitude as the value calculated for the parameters of Figure 3.

These results are qualitatively consistent with overall features of the measured optical conductivity, including the overall magnitude Nagai et al. (2001); Singley et al. ; Katsumoto et al. (2001). The accuracy of our numerical results is limited somewhat by the wavevector cutoffs that we must invoke in order to make the numerical calculations manageable. This is especially troublesome in estimating the importance of heavy-hole to split-off band transitions. The gradual increase in absorption beyond 500 meV seen in experiment may reflect these transitions. We expect that allowing more wavevectors in our calculations would shift oscillator strength toward higher energies and produce a broader peak, in better agreement with experiments. (Increasing $`n_m`$ from our current value, 3, by just one changes the dimension of the Hamiltonian matrix from 2058 to 4374, and the resulting self-consistent computation becomes quite substantial.) On the other hand, our investigation of parabolic band case, where $`n_m`$ can be increased to 4, suggests that system size dependence of low energy non-Drude part of absorption is negligible for $`n_m`$ bigger than 3.

## V Conclusions

We have presented a non-perturbative self-consistent study of the optical properties of diluted (III,Mn)V DMS which treats both disorder and interactions on equal footing. Our calculations for the diagonal ac-conductivity are in good qualitative or semi-quantitative agreement with optical absorption experiments in thin film geometries. Several of the non-Drude behaviors observed in the experiments are naturally accounted for in this approach which is able to capture multiple scattering effects that are beyond the scope of perturbative Born approximation approaches Sinova et al. (2002). Our calculations show that spin-orbit coupling of valence band holes must be included in the calculation in order to account for the qualitative features of the experiments. It is possible that features associated with heavy-hole to split-off band absorption will emerge more clearly near 0.9 eV when $`\sigma (\omega )`$ measurements are performed on the most metallic (Ga,Mn)As ferromagnets currently available.

We have used these calculations to test the efficacy of the f-sum rule we proposed in earlier work for carrier density measurements. In applying the f-sum rule to the present system a complication arises, compared to microscopic and parabolic band cases, in that the integrated conductivity depends on the character of the occupied quasiparticle states. In previous work Sinova et al. (2002) we estimated that for carrier densities in the range that is interesting for (Ga,Mn)As ferromagnets, the effective mass that should be used in the f-sum rule is $`0.24m_e`$, intermediate between heavy and light-hole masses. This estimate was based on calculations of the occupied quasiparticle states that treated disorder perturbatively in the Born approximation. We have found by integrating the optical conductivities evaluated here over frequency that this optical effective mass estimate is reliable at the $`10\%`$ level, when disorder is treated exactly. Uncertainties introduced by our finite momentum cut-off do not allow us to estimate the value or the variability of this mass more accurately. Nevertheless, this level of accuracy should be comparable with what is possible from Hall effect measurements which are complicated by the strong anomalous Hall effect in these ferromagnets.

###### Acknowledgements.

This work was supported in part by KOSEF Quantum-functional Semiconductor Research Center at Dongkuk University, the Welch Foundation, DARPA, the Grant Agency of the Czech Republic under grant 202/02/0912, and the Ministry of Education of the Czech Republic under grant OC P5.10.


text:
# Interplay between Vorticity and Chirality inside the Vortex Core in Chiral 𝑝-Wave Superconductors

## 1 INTRODUCTION

Site-selective nuclear magnetic resonance (NMR) method is a powerful tool for investigating the electronic structure inside vortex cores in the mixed state of type-II superconductors. $`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}\text{}\text{?}}`$ We have theoretically studied $`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}}`$ the site-selective nuclear spin-lattice relaxation rate $`T_1^1`$ inside a vortex core in the case of an isotropic chiral $`p`$-wave superconductivity, $`\mathrm{\backslash }sevenrm^\text{?}`$ $`𝐝=\overline{𝐳}(\overline{k}_x\pm i\overline{k}_y)`$. We found that $`T_1^1`$ was suppressed and almost vanished in the $`\overline{k}_xi\overline{k}_y`$ state owing to the interplay between the vorticity and chirality inside the vortex core (here, the magnetic field was applied in positive direction of the $`z`$ axis). $`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}}`$ In this paper, we investigate the effect of a tilt of the magnetic fields on $`T_1^1`$. We then investigate the effect of the anisotropy in the superconducting gap and the Fermi surface. Our result is expected to be experimentally observed as a sign of the chiral pairing state in a superconducting material Sr2RuO4$`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}}`$

## 2 QUASICLASSICAL THEORY

To investigate $`T_1^1`$, we utilize the quasiclassical theory of superconductivity. We consider the quasiclassical Green function

$$\widehat{g}(i\omega _n,𝐫^{},\overline{𝐤})=i\pi \left(\begin{array}{cc}g& if\\ if^{}& g\end{array}\right),$$
(1)
which is the solution of the Eilenberger equation, $`\mathrm{\backslash }sevenrm^\text{?}`$

$$i𝐯_\mathrm{F}(\overline{𝐤})\widehat{g}+[i\omega _n\widehat{\tau }_z\widehat{\mathrm{\Delta }},\widehat{g}]=0,$$
(2)
where the superconducting order parameter is $`\widehat{\mathrm{\Delta }}(𝐫^{},\overline{𝐤})=\left[(\widehat{\tau }_x+i\widehat{\tau }_y)\mathrm{\Delta }(𝐫^{},\overline{𝐤})(\widehat{\tau }_xi\widehat{\tau }_y)\mathrm{\Delta }^{}(𝐫^{},\overline{𝐤})\right]/2`$ and $`\widehat{\tau }_i`$ the Pauli matrices. $`𝐯_\mathrm{F}(\overline{𝐤})`$ is the Fermi velocity. The vector $`𝐫^{}=(x^{},y^{},z^{})=(r\mathrm{cos}\varphi ,r\mathrm{sin}\varphi ,z^{})`$ is the center of mass coordinate, where the magnetic field is applied along the $`z^{}`$ axis. The unit vector $`\overline{𝐤}=(\overline{k}_x,\overline{k}_y,\overline{k}_z)=(k_{}\mathrm{cos}\theta ,k_{}\mathrm{sin}\theta ,k_z)/\sqrt{k_{}^2+k_z^2}`$ represents the wave number of relative motion of the Cooper pairs in the crystallographic coordinate frame.

From the spin-spin correlation function, $`\mathrm{\backslash }sevenrm^\text{?}`$ we obtain the expression for $`T_1^1`$ in terms of $`\widehat{g}`$,

$$\frac{T_1^1(T)}{T_1^1(T_\mathrm{c})}=\frac{1}{4T_\mathrm{c}}_{\mathrm{}}^{\mathrm{}}\frac{\mathrm{d}\omega }{\mathrm{cosh}^2(\omega /2T)}W(\omega ,\omega ),$$
(3)
$$W(\omega ,\omega ^{})=a_{11}(\omega )a_{22}(\omega ^{})a_{12}(\omega )a_{21}(\omega ^{}),$$
(4)
where the spectral function $`\widehat{a}(\omega ,𝐫^{},\overline{𝐤})=\left(a_{ij}\right)`$ is given as

$$\widehat{a}(\omega ,𝐫^{},\overline{𝐤})=\frac{i}{2\pi }\widehat{\tau }_3\left[\widehat{g}(i\omega _n\omega i\eta ,𝐫^{},\overline{𝐤})\widehat{g}(i\omega _n\omega +i\eta ,𝐫^{},\overline{𝐤})\right],$$
(5)
the symbol $`\mathrm{}`$ represents the average over the Fermi surface, and $`\eta `$ is a small positive constant roughly representing the impurity effect. $`T`$ is the temperature and $`T_\mathrm{c}`$ the superconducting transition temperature.

## 3 EFFECT OF TILT OF MAGNETIC FIELD

We define $`\gamma `$ as the angle between the magnetic field ($``$ $`z^{}`$ axis) and the $`z`$ crystallographic axis ($``$ $`\overline{k}_z`$ axis). The isotropic cylindrical Fermi surface is assumed and the Fermi velocity is $`𝐯_\mathrm{F}(\overline{𝐤})=(v_\mathrm{F}\mathrm{cos}\theta ,v_\mathrm{F}\mathrm{sin}\theta ,0)`$ in the crystallographic coordinate frame. We assume that the layer perpendicular to the $`z`$ crystallographic axis is isotropic. In this section, we consider the isotropic chiral $`p`$-wave pairing, $`\mathrm{\backslash }sevenrm^\text{?}`$ $`𝐝=\overline{𝐳}(\overline{k}_x\pm i\overline{k}_y)=\overline{𝐳}\mathrm{exp}(\pm i\theta )`$. We assume that this pairing is unchanged by a tilt of the magnetic field.

The order parameter around a vortex is expressed as $`\mathrm{\Delta }(𝐫^{},\overline{𝐤})=|\mathrm{\Delta }(r)|\mathrm{exp}(i\varphi )\mathrm{exp}(\pm i\theta )`$. On the basis of an analysis of the so-called zero-core vortex model, $`\mathrm{\backslash }sevenrm^\text{?}`$ the matrix elements of $`\widehat{g}`$ at the vortex center $`r=0`$ are approximately obtained as $`\mathrm{\backslash }sevenrm^\text{?}`$ $`g=\sqrt{\omega _n^2+|\stackrel{~}{\mathrm{\Delta }}|^2}\omega _n^1,f=\stackrel{~}{\mathrm{\Delta }}\omega _n^1,f^{}=\stackrel{~}{\mathrm{\Delta }}^{}\omega _n^1`$, where $`\stackrel{~}{\mathrm{\Delta }}=|\mathrm{\Delta }(r\mathrm{})|\mathrm{exp}(i\varphi )\mathrm{exp}(\pm i\theta )`$ and $`(\mathrm{cos}\varphi ,\mathrm{sin}\varphi )(v_{\mathrm{F}x^{}},v_{\mathrm{F}y^{}})`$ in the plane perpendicular to the $`z^{}`$ axis, i.e., to the magnetic field. Without loss of generality, we tilt the magnetic field in the $`z`$-$`y`$ plane. In the present case, it follows that $`\mathrm{cos}\varphi =\mathrm{cos}\theta `$ and $`\mathrm{sin}\varphi =\mathrm{sin}\theta \mathrm{cos}\gamma `$.

Inserting these matrix elements of $`\widehat{g}`$ into Eq. (5) and taking the Fermi-surface average $`\mathrm{}=\mathrm{}d\theta /2\pi `$ in Eq. (4), we numerically calculate $`T_1^1(T)`$ at the vortex center. We show the result for $`\gamma =85^{}`$ in Fig. 1. Increasing the tilt angle $`\gamma `$, we observe in our numerical results that $`T_1^1(T)`$ are almost unchanged by $`\gamma `$ up to $`\gamma 50^{}`$. Above $`\gamma 50^{}`$ ($`\gamma <90^{}`$), each $`T_1^1(T)`$ in the $`\overline{k}_x\pm i\overline{k}_y`$ states gradually deviates from $`T_1^1(T)`$ of $`\gamma =0^{}`$ and the difference in $`T_1^1(T)`$ between the two chiral states gradually becomes small. At $`\gamma =90^{}`$, $`T_1^1(T)`$ of the two chiral states coincide each other. As seen in Fig. 1, however, even at a large tilt angle, the difference in $`T_1^1(T)`$ between the $`\overline{k}_x\pm i\overline{k}_y`$ states is still noticeable. At $`\gamma =85^{}`$ (the solid lines in Fig. 1), $`T_1^1(T)`$ in the $`\overline{k}_x+i\overline{k}_y`$ state is two times larger at $`T0.4T_\mathrm{c}`$ than that in the $`\overline{k}_xi\overline{k}_y`$ state.

## 4 EFFECT OF ANISOTROPY

Within the present framework based on the quasiclassical theory, we calculate $`T_1^1(T)`$ using an anisotropic gap $`𝐝=\overline{𝐳}(\mathrm{sin}k_x\pm i\mathrm{sin}k_y)`$ and an anisotropic dispersion relation $`\epsilon _k=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)2t^{}\left[\mathrm{cos}(k_x+k_y)+\mathrm{cos}(k_xk_y)\right]`$, $`(`$$`t^{}=0.47t`$ and the chemical potential $`\mu =1.2t`$ $`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}}`$$`)`$. In this section, the magnetic field is applied in positive direction of the $`z`$ crystallographic axis, namely $`\gamma =0^{}`$. This physical situation is the same as that of a calculation of $`T_1^1(T)`$ by Takigawa et al. $`\mathrm{\backslash }sevenrm^\text{?}`$

We show our result for $`T_1^1(T)`$ in Fig. 2. We find that, even in this case of the anisotropic gap and the anisotropic Ferm surface, the difference in $`T_1^1(T)`$ between the $`\mathrm{sin}k_x\pm i\mathrm{sin}k_y`$ states is noticeable and $`T_1^1(T)`$ in the $`\mathrm{sin}k_xi\mathrm{sin}k_y`$ state is suppressed in wide $`T`$ region in comparison with that in the $`\mathrm{sin}k_x+i\mathrm{sin}k_y`$ state. This result is in contrast to a corresponding theoretical result for $`T_1^1(T)`$ by Takigawa et al. $`\mathrm{\backslash }sevenrm^\text{?}`$ In their result, there is not such suppression of $`T_1^1(T)`$ as seen in our Fig. 2. While the same anisotropy is taken into account in both calculations, their result is different from ours. A reasonable origin of this difference may be as follows. The calculation of $`T_1^1(T)`$ by Takigawa et al. $`\mathrm{\backslash }sevenrm^\text{?}`$ is in the quantum limit ($`k_\mathrm{F}\xi 1`$) where the energy spectrum inside the vortex core is quantized and it dominantly determines $`T_1^1(T)`$ as pointed out by them, $`\mathrm{\backslash }sevenrm^\text{?}`$ while we base our calculation on the quasiclassical theory relevant in the opposite limit $`k_\mathrm{F}\xi 1`$ where the vortex core spectrum is continuous. Now, in the case of the material Sr2RuO4$`\mathrm{\backslash }sevenrm^{\text{?}\text{}\text{?}}`$ the coherence length is not so small ($`\xi 660`$ Å $`\mathrm{\backslash }sevenrm^\text{?}`$), namely $`k_\mathrm{F}\xi 1`$, and therefore our result based on the quasiclassical theory is relevant to this material.

## 5 SUMMARY

Within the framework of the quasiclassical theory of superconductivity, we numerically calculated the site-selective nuclear spin-lattice relaxation rate $`T_1^1`$ at the vortex center in the chiral $`p`$-wave superconductors. The case of the tilted magnetic field and the case of the anisotropic gap and Fermi surface were investigated. Our result (Fig. 1) can be experimentally observed as a sign of the chiral pairing state in Sr2RuO4 by applying the magnetic field tilted from the $`z`$ crystallographic axis by $`85^{}`$.

## ACKNOWLEDGMENTS

We thank M. Ichioka, M. Takigawa, M. Matsumoto, K. Machida, N. Schopohl, and M. Sigrist for helpful discussions.