text:
# Criticality in non-linear one-dimensional maps: RG universal map and non-extensive entropy

## 1 Introduction

The proposal, in 1988, by Tsallis  of a non-extensive generalization of statistical mechanics unleashed an unprecedented discussion  on the foundations of this most basic branch of physics. This is an exceptional event, since its inception, with Boltzmann and Gibbs (BG), in the long and heretofore unchallenged history of statistical mechanics. Many recent studies  have offered experimental and numerical evidences indicating limiting conditions for the domain of validity of the canonical BG theory. But, so far there have been no firm developments to tilt the balance towards universal acceptance of these new ideas, nor a definite understanding of the deep-lying physical reasons (believed in some cases to be a breakdown in the chain of increasing randomness from non-ergodicity to completely developed chaoticity) for the failure of the BG statistics and the competence of the non-extensive generalization. However, here we present exact analytical results for a class of critical points in non-linear maps from which a rigorous corroboration and understanding of the non-extensive theory can be drawn  \- .

The renowned period-doubling and intermittency routes to chaos are now classical examples of how the use as starting points of simple non-linear discrete maps have often led to important developments in the theory of non-linear dynamical systems . The discovery of the universal properties associated to these mechanisms, comparable with those of conventional critical phenomena in statistical physics, triggered, over two decades ago, an explosion of activity in the field and now both routes are well understood and verified experimentally. Here we might add yet other universal aspects to their properties (explained in more detail below) the non-extensivity marked by the power-law growth or decay of the sensitivity to initial conditions, and the entropy extrema associated to each of the fixed points that define the RG universality classes.

The well-known period-doubling and intermittency transitions are based on the pitchfork and the tangent bifurcations. In particular, the one-dimensional logistic map and its generalization to non-linearity of order $`𝔷>1`$, $`f_\mu (x)=1\mu \left|x\right|^𝔷`$,$`1x1`$, exhibit these two different types of bifurcations an infinite number of times as the control parameter $`\mu `$ reaches critical values. At these critical points universal scaling laws and indexes hold, independently of the details of the map, and are therefore shared by all maps of the same order $`𝔷`$. As with other statistical mechanical systems an explanation of universality is provided via the renormalization group (RG) method, which has been successfully applied to this type of iterated map. The RG doubling transformation, consisting of functional composition and rescaling, was first devised to study the cascade of period-doubling transitions and its accumulation point, and was later applied to the intermittency transition. As pointed out  the latter case is one of the rare examples where the RG equations can be solved exactly. Here we indicate that this analytic solution holds generally not only for tangent bifurcations (as known time ago ) but also for pitchfork bifurcations (as noticed recently ). We also point out that the fixed-point map can be determined as an extremum of the non-extensive entropy, thus providing evidence for a possible connection between extremal properties of entropy expressions and those for the fixed points obtained from the usual RG recursive method .

As required, at each of the map critical points the Lyapunov $`\lambda _1`$ exponent vanishes, and the sensitivity to initial conditions $`\xi _t`$ for large iteration time $`t`$ ceases to obey exponential behavior, exhibiting instead power-law behavior . As the means to describe the dynamics at such critical points, and based on the non-extensive entropy of Tsallis $`S_Q`$  (see Eq. (16) below), the $`q`$-exponential expression

$$\xi _t=\mathrm{exp}_q(\lambda _qt)[1(q1)\lambda _qt]^{1/(q1)}$$
(1)
containing a $`q`$-generalized Lyapunov exponent $`\lambda _q`$ was proposed , together with generalizations for the rate of entropy production (referred here as the Kolmogorov-Sinai (KS) entropy ) $`K_Q`$ (see Eq. (15) below) and for the Pesin identity $`\lambda _q=`$ $`K_Q`$ . The standard expressions are recovered when $`q,Q1`$. The use of the two indexes $`q`$ and $`Q`$ in our application will become clear below. Furthermore, recent studies -, focusing mainly at the onset of chaos (the accumulation point of the period-doubling bifurcations) of the $`𝔷`$-logistic map have revealed a series of connections between the Tsallis entropic index $`q`$ and the non-linearity $`𝔷`$, the fractal dimension $`d_f`$ of the chaotic attractor, and the end points of its multifractal singularity spectrum $`f(\alpha )`$. Here we identify $`q`$ and $`\lambda _q`$ in terms of the tangent and pitchfork critical map parameters.

The properties of dynamical systems at critical points can also be probed by considering the deterministic diffusive processes that take place in non-linear maps with discrete translational symmetry. These maps consist of a sequence of cells containing each elements of a basic nonlinear map. One such construction  makes use of the tangent bifurcation, and, interestingly, anomalous diffusional behavior has been observed close to the intermittency transition , . Continuous-time random-walk theory concepts have been applied ,  and, characteristically, the associated waiting-time distribution function $`\psi (t)`$ has been shown to have non-exponential power-law decay. The validity of the non-extensive approach to criticality has been recently checked by examination of the anomalous diffusion at the intermittency transition , . The key results in these studies center on the analytical form predicted for $`\psi (t)`$ which bears a close relationship to the power-law Eq. (1). Here we point out that the analytical forms for this and related distributions actually correspond to those obtained by the repeated application of the RG transformation. As mentioned, if the fixed point solution can also be determined via an entropy optimization method the expressions for these distributions are indicative of the same entropy extremum.

## 2 Outline of results

The main points in the following analysis are:

1) The connection between the RG fixed-point map expression for the iterate $`x_t`$ associated to trajectories $`x_t`$ with $`x_00`$ at $`t=0`$ and the expression in Eq. (1) for the sensitivity of initial conditions $`\xi _t`$. This leads to simple identifications of the non-extensive parameters as $`q=2z^1`$ and $`\lambda _qa`$. We corroborate the earlier known results $`q=3/2`$ and $`q=5/3`$ for the tangent and pitchfork bifurcations for the logistic map $`𝔷=2`$, respectively .

2) The connection between the RG fixed-point map recursion relation for the iterate $`x_t`$ and the expression obtained from $`S_Q`$ for the temporal evolution of the number of cells $`W_t`$ occupied by an ensemble of initial conditions spread over a finite-length interval $`[0,x_0]`$. (See the definition of $`S_Q`$ in Eq. (16), and the explanatory comments at the end of this section about the relationship $`q=2Q^1`$ between the indexes $`q`$ and $`Q`$). This leads to the identifications $`Q=z`$ and $`K_qa`$, confirming $`q=2z^1`$ and giving support to the validity of the generalized Pesin identity $`\lambda _q=`$ $`K_Q`$ at the critical transitions. See Ref.  for a recent proof of the validity of this generalized identity at the edge of chaos of the logistic map $`𝔷=2`$.

3) The general property of the non-extensive entropy function which decreases monotonically as the RG transformation flows away from the fixed point where it attains its maximum value.

4) The maximum entropy property of the distribution $`\mathrm{\Psi }(\tau )`$, $`\mathrm{\Psi }(\tau )=_0^\tau 𝑑t\psi (t)`$, of cell residence time intervals $`[0,\tau ]`$ in the critical diffusive processes.

For clarity of presentation it is convenient to advance the following property of the $`q`$-exponential function. For this function the ordinary derivative $`d\mathrm{exp}(x)/dx=\mathrm{exp}(x)`$ and inverse $`[\mathrm{exp}(x)]^1=\mathrm{exp}(x)`$ properties become, respectively, $`d\mathrm{exp}_q(x)/dx=[\mathrm{exp}_q(x)]^q`$ and $`[\mathrm{exp}_q(x)]^1=\mathrm{exp}_{2q}(x)`$, and these combine into the identity

$$\mathrm{exp}_q(x)\left[\mathrm{exp}_Q(\frac{x}{Q})\right]^Q,$$
(2)
where $`q=2Q^1`$. Thus, in a given problem where the $`q`$-exponential function applies there appear related values of the $`q`$-index, that become the same (and the above identity becomes trivial) as $`q1`$.

## 3 Critical dynamics from fixed-point map

So, we recall ,  the solution of Hu and Rudnick to the Feigenbaum RG recursion relation obtained for the case of the tangent bifurcation. For the transition to periodicity of order $`n`$ consider the $`n`$-th composition $`f^{(n)}`$ of the original map $`f`$ in the neighborhood of one of the $`n`$ points tangent to the line with unit slope and shift the origin of coordinates to that point. Then, one obtains

$$f^{(n)}(x)=x+u\left|x\right|^z+o(\left|x\right|^z),$$
(3)
where $`u>0`$ is the leading expansion coefficient (the right-hand-side in Eq. (3) may depend on $`n`$ only through the omitted terms). The RG fixed-point map $`x^{}=f^{}(x)`$ was found to be

$$x^{}=x\mathrm{exp}_z(ux^{z1})=x[1(z1)ux^{z1}]^{1/(z1)},$$
(4)
as it satisfies $`f^{}(f^{}(x))=\alpha ^1f^{}(\alpha x)`$ with $`\alpha =2^{1/(z1)}`$ and has a power-series expansion in $`x`$ that coincides with Eq. (3) in the two lowest-order terms. (We have used $`x^{z1}\left|x\right|^{z1}\mathrm{s}gn(x)`$). Eq. (4) can be rewritten as $`x^{(z1)}=x^{(z1)}(z1)u`$, a translation of the power law $`y=`$ $`x^{(z1)}`$, the property that led Hu and Rudnick in their derivation. The effect of a perturbation of the form $`x^p`$ transforms this into

$$x^{(z1)}+\epsilon ^{}x^p=x^{(z1)}+\epsilon x^p(z1)u,$$
(5)
or, to lowest order in $`\epsilon ^{}=\epsilon `$,

$$x^{}=x\mathrm{exp}_z(ux^{z1})\frac{\epsilon x^{zp}}{z1}\left[\mathrm{exp}_z(ux^{z1})\right]^z$$

$$+\frac{\epsilon x^{zp}}{z1}\left[\mathrm{exp}_z(ux^{z1})\right]^{zp}.$$
(6)
The recurrence relation in Eq. (5) remains invariant under iteration and rescaling by $`\alpha `$ when $`\epsilon ^{}=\epsilon `$ is multiplied by the factor $`\gamma =2^{(pz+1)/(z1)}`$. As it is known , , the leading relevant perturbation about this critical point is that which corresponds to a shift from tangency, for this case the lowest order term in $`x`$ in Eq. (6) must be of order one, and so we fix $`p=z`$ for the perturbation exponent (in the unormalized eigenfunction).

Fittingly, it can be readily concluded that the RG recursion relation in the neighborhood of the pitchfork bifurcations $`df^{(2^{k1})}(x)/dx|_{x=0}=`$ $`1`$ conforms with the above scheme and that its solution is also applicable to these transitions provided the sign of $`u`$ is changed for $`x>0`$. It is interesting to point out that at the time the RG method was first applied to the tangent bifurcation it was thought that the simple closed form obtained for the fixed point map, the $`q`$-exponential in Eq. (4), would most likely prove not to be fruitful for the study of period doubling , . Not only is this solution applicable to each period doubling transition with the minor modification of the indicated sign change, but also, as it has been recently demonstrated , the time evolution at their accumulation point, the edge of chaos, is expressable in terms of similar $`q`$-exponentials. Below we comment on the caveats that appear in relating $`z`$ above to the $`𝔷`$ in the initial map e.g. $`f_\mu (x)=1\mu \left|x\right|^𝔷`$.

When the number of iterations (after an initial $`x_0`$) is large and the continuous time $`t`$ approximation is taken ($`x^{}=x_t`$, $`x_0=x`$ and $`u=at`$ in Eq. (4)) the fixed-point map can be written as

$$x_t=x_0[1(z1)ax_0^{z1}t]^{1/(z1)},$$
(7)
with properties $`dx_t/dt=ax_t^z`$and$`dx_t/dx_0=(x_t/x_0)^z`$. Thus, the fixed-point and perturbation relations can be used to describe the time evolution of intervals containing iterates, an initial interval $`[0,x_0<x_+]`$ of length $`\mathrm{\Delta }x_0`$ evolves into $`[0,x_t<x_+]`$ with length $`\mathrm{\Delta }x_t`$ where $`x_t`$ is the image of $`x_0`$ at time $`t`$. Accordingly, Eq. (5) with ($`p=z`$) becomes

$$\mathrm{\Delta }x_t^{(z1)}+\epsilon _t\mathrm{\Delta }x_t^z=\mathrm{\Delta }x_0^{(z1)}+\epsilon _0\mathrm{\Delta }x_0^z(z1)at.$$
(8)
We should like to emphasize that the expression for time evolution Eq. (7) has the same universal form as that for the fixed-point map Eq. (4). In the neighborhood of the points of tangency $`f^{(n)}(0)`$ the iterates follow monotonic paths with time shape set by the map shape itself, i.e. long time dynamics follow the static solution of Hu and Rudnick. It should also be kept in mind that Eq. (8) applies only to small perturbations as this corresponds to the perturbation recursion relation obtained to linear order around the RG fixed point. We also comment that our use here of the continuous time approximation (valid for large iteration time $`t`$) in deriving Eq. (4) is only for expediency reasons and has no effect in the generality of our results as this equation is equally obtained from the scaling property

$$x_m=f^{^{(m)}}(x)=m^{1/z1}f^{}(m^{1/z1}x),m=1,2\mathrm{},$$
(9)
that in turn is obtained from $`f^{}(f^{}(x))=\alpha ^1f^{}(\alpha x)`$ with $`\alpha =2^{1/z1}`$ .

## 4 Sensitivity to initial conditions

We have reached a point at which a straightforward proof of the validity of Tsallis’ Eq. (1) can be obtained. This is done by direct evaluation of the sensitivity to initial conditions from its definition $`\xi _tlim_{\mathrm{\Delta }x_00}\left(\mathrm{\Delta }x_t/\mathrm{\Delta }x_0\right)`$. From Eq. (7) and $`dx_t/dx_0=(x_t/x_0)^z`$ one obtains

$$\xi _t(x_0)=[1(z1)ax_0^{z1}t]^{z/(z1)},$$
(10)
and comparison with Eq. (1) yields $`q=2z^1`$ and $`\lambda _q(x_0)=zax_0^{z1}`$.

We can also confirm the values of $`q`$ obtained numerically by Tsallis et al  for both the tangent and pitchfork bifurcations when $`𝔷=2`$ . The composition $`f_\mu ^{(n)}`$ for the $`𝔷`$-logistic map $`f_\mu (x)=1\mu \left|x\right|^𝔷`$,$`𝔷>1`$, at the tangent bifurcations is $`f_\mu ^{(n)}(x)=x+ux^2+o(x^2)`$,$`u>0`$, and from $`z=2`$ one has $`q=3/2`$. For the pitchfork bifurcations one has instead $`f_\mu ^{(n)}(x)=x+ux^3+o(x^3)`$, because $`d^2f_\mu ^{(2^k)}/dx^2=0`$ at these transitions and $`u<0`$ is now the coefficient associated to $`d^3f_\mu ^{(2^k)}/dx^3<0`$. In this case we have $`z=3`$ in $`q=2z^1`$ and one obtains $`q=5/3`$. Notably, these specific results for the index $`q`$ are valid for all $`𝔷>1`$ and therefore define the existence of only two universality classes for unimodal maps, one for the tangent and the other one for the pitchfork bifurcations .

Further, if one uses the fact  that the invariant distribution of $`f^{(n)}`$ in Eq. (3) is of the form $`\rho (x)x^{(z1)}`$, the average $`\overline{\lambda }_q`$ of $`\lambda _q(x_0)`$ over $`x_0`$ yields the constant $`\overline{\lambda }_qa`$. In relation to this we observe that this average is compatible with

$$\overline{\lambda }_q=𝑑x\rho (x)\mathrm{ln}_q\left|\frac{df^{}(x)}{dx}\right|,$$
(11)
the $`q`$-extension of the customary  expression for $`\lambda _1`$ as the average of $`\mathrm{ln}\left|df^{}(x)/dx\right|`$ over $`\rho (x)`$. By taking the $`q`$-logarithm (where $`\mathrm{ln}_qy(y^{1q}1)/(1q)`$ is the inverse of $`\mathrm{exp}_q(y)`$) of both sides of

$$\xi _t(x_0)=[1(q1)\lambda _q(x_0)t]^{1/(q1)},$$
(12)
and by recalling the definition of $`\xi _t(x_0)`$ one obtains

$$\lambda _q(x_0)=t^1\mathrm{ln}_q\left|\frac{df^{(t)}(x_0)}{dx_0}\right|,$$
(13)
where $`\mathrm{ln}_q\left|df^{(t)}(x_0)/dx_0\right|t`$ as $`\lambda _q(x_0)`$ does not depend on $`t`$. (Notice that the limit $`t\mathrm{}`$ of Eq. (13) is the $`q`$-extension of the regular definition of $`\lambda _1(x_0)`$ ). Next, use of the fixed-point relation $`f^{(t)}(x_0)=t^{1/1z}f^{}(t^{1/z1}x_0)`$ and $`t=(x/x_0)^{z1}`$ into (13) leads to $`\lambda _q(x_0)=(x/x_0)^{(z1)}\mathrm{ln}_q\left|df^{}(x)/dx\right|`$ and Eq. (11).

## 5 Rate of entropy change

Even though our basic results have been derived, it is useful to extend our analysis and reveal more details about the dynamics at these critical points. Our following arguments focus on $`x_t`$ rather than $`dx_t/dx_0`$ and imply $`Q=z`$, thus confirming $`q=2z^1`$. They also indicate that the rate of entropy change is $`K_Qa`$, suggesting the validity of the $`q`$-generalization of the Pesin $`\overline{\lambda }_q=`$ $`K_Q`$ at these critical states. Therefore, let’s consider a cell partition ($`N\mathrm{}`$ cells of length $`l0`$) in the ’entire’ phase space $`[0,x_+]`$. Consider also the time evolution of occupied cells $`W_t`$ by an ensemble of iterates with an extended initial occupation $`W_0`$ that spans $`[0,x_0<x_+]`$. We write the cell distribution of iterates $`p_i`$ as

$$p_i=\frac{1+\delta _i}{W},\underset{i}{\overset{W}{}}\delta _i=0,$$
(14)
as this is a convenient form when considering normalized perturbations to the uniform distribution. It is a general expression when no further conditions are imposed on the $`\delta _i`$ but below we require small $`\delta _i`$. When Eq. (14) is introduced into the KS $`q`$-entropy difference

$$K_Qt=S_Q(t)S_Q(0),t\text{large,}$$
(15)
where

$$S_Q=\frac{1_i^Wp_i^Q}{Q1}$$
(16)
is the Tsallis entropy, we obtain

$$W_t^Q\underset{i}{\overset{W_t}{}}(1+\delta _{i,t})^Q=W_0^Q\underset{i}{\overset{W_0}{}}(1+\delta _{i,0})^Q(Q1)K_Qt.$$
(17)
To lowest order in $`\delta _i`$ this reduces to

$$W_t^{(Q1)}+\frac{1}{2}Q(Q1)\delta _t^2W_t^Q=$$

$$W_0^{(Q1)}+\frac{1}{2}Q(Q1)\delta _0^2W_0^Q(Q1)K_Qt,$$
(18)
with $`\delta ^2_i^W`$ $`\delta _i^2`$.

Remarkably, Eqs. (18) and (8) can be seen to be equivalent if $`Q=z`$ and if $`W_{t,0}`$ $`\mathrm{\Delta }x_{t,0}`$, $`K_Qa`$ and $`\frac{1}{2}Q(Q1)\delta _{t,0}^2\epsilon _{t,0}`$. The entropic index $`Q`$ is simply given by the non-linearity $`z`$ while $`K_Q`$ is the expansion coefficient $`a`$ (and through the proviso that $`\overline{\lambda }_q=`$ $`K_Q`$ this would imply $`\overline{\lambda }_qa`$). Our comparison between Eqs. (18) and (8) is consistent in that each is the lowest-order expression valid for small perturbations, the first to the uniform cell distribution of iterates and the second to the fixed-point map solution. We can draw the following conclusion from the association of Eqs. (18) and (8): The entropy $`S_Q^{}`$ for the fixed-point map is maximum, as $`S_Q^{}S_Q=Q\delta ^2/2W^q>0`$, and the distribution becomes uniform, $`p_i=W^1`$. The RG flow is away from the fixed-point since $`\gamma >1`$ and the entropy decreases as the RG transformation is applied. The proposed equivalence between Eqs. (18) and (8) is justified by making the following connection with Eq. (1). For a uniform distribution $`W_{t,0}`$ $`\mathrm{\Delta }x_{t,0}`$, as the ratio of the number of cells $`W_t`$ occupied at time $`t`$ to the distance $`\mathrm{\Delta }x_t`$ is the length $`l`$, and in the limit $`N\mathrm{}`$ and $`l0`$, so that $`\mathrm{\Delta }x_00`$, with say $`\mathrm{\Delta }x_0/l=b`$, one obtains $`W_t=b\xi _t`$. This explicit link between the occupation numbers $`W_t`$ and $`\xi _t`$ is in line with the all-important connection between the loss of information embodied by the KS entropy and the sensitivity to initial conditions in non-critical regimes , and with the parallel connection suggested by Tsallis et al ,  for the analogous quantities at critical points.

It is worth the mention that the use of $`q`$-exponential (and the inverse $`q`$-logarithm) functions in the derivation of the properties of a given system will in general introduce pairs of conjugate indexes such as $`q=2Q^1`$ above. As pointed out, taking derivatives and inverse operations lead to expressions involving one or the other of these simply related values of the entropic index (that become trivially the same $`q=Q=1`$ for the usual exponential and logarithmic functions). Consequently, while some theoretical features are expressed through $`q`$ some others appear expressed via $`Q`$, but in all cases their occurrence and meaning is apparent.

## 6 Deterministic diffusion

We comment now on the deterministic diffusion process that takes place for the variable $`x_t`$ in the map $`x_{t+1}=f(x_t)`$, where $`f(m+x)=m+f(x)`$, $`m`$ integer and $`f(x)=f(x)`$. Geisel and Thomae  found that when the unit-cell map $`f(0x1/2)`$ is chosen to have the (intermittency transition) form $`x_{t+1}=x_t+ax_t^z`$ within $`[0,x_+1/2]`$ anomalous diffusion, of the dispersive type, is generated when $`z>2`$ since then the mean-square displacement $`<\mathrm{\Delta }x^2(t)>t^{1/(z1)}`$. This behavior follows  after the distribution of cell residence times $`\psi (t)`$ is found to be

$$\psi (t)=ax_+^{z1}[1+(z1)ax_+^{z1}t]^{z/(z1)}.$$
(19)
Grigolini and co-workers  recognized the relationship between the expressions for $`\psi (t)`$ and $`\xi _t`$. They also analyzed, and verified, the validity of the non-extensive scheme involving the $`q`$-generalized KS entropy its equality with the $`q`$-generalized Lyapunov exponent , . Here we can recap these results as follows: The properties of the diffusion process driven by the intermittency transition are universal since they stem from the RG fixed-point map expression. Namely, $`\psi (t)`$ is, besides normalization, the derivative of the inverse fixed-point map $`dx_0/dx_t`$ with $`x_t=x_+`$ and is therefore inversely proportional to $`\xi _t=`$ $`dx_t/dx_0`$. The distribution $`\mathrm{\Psi }(\tau )=_0^\tau 𝑑t\psi (t)<\mathrm{\Delta }x^2(\tau )>`$ is in turn proportional to the inverse of the fixed-point map $`x_0(x_\tau )`$. Also,

$$\mathrm{\Psi }(\tau )=A[1+(z1)ax_+^{z1}\tau ]^{1/(z1)},$$
(20)
has the form of the distribution obtained by optimization of the non-extensive entropy with index $`Q=z`$ for fixed $`ax_+^{z1}\tau `$. The related entropy index $`q=2Q^1=2z^1`$ applies to the residence time distribution $`\psi (t)`$.

## 7 Summary and discussion

Thus, for the unimodal maps analyzed we have indicated that the exact static solution of the RG equations for the tangent and pitchfork bifurcations also describe the dynamics of iterates. Further, these fixed-point expressions are the same $`q`$-exponential expressions describing the temporal evolution of ensembles of iterates prescribed by the non-extensive formalism for non-linear maps at criticality. Therefore the RG fixed-point map for the universality classes determined by $`z`$ can be constructed via Tsallis entropic arguments. The equivalence rests on the derivation of Eq. (1) exclusively from RG procedures. The study of the pitchfork and period-doubling transitions has been expanded via detailed derivation of their $`q`$-generalized Lyapunov exponents $`\overline{\lambda }_q`$ and interpretation of the different types of sensitivity $`\xi _t`$ . Likewise, the properties of the intricate trajectories at the edge of chaos of the logistic map $`𝔷=2`$ have been analytically obtained leading too to the determination of $`\overline{\lambda }_q`$ and interpretation of the dynamics at the strange attractor . Further, the validity of the $`q`$-generalized Pesin identity at this important state has been recently rigorously proved . More generally, we addressed here the properties of the fixed-point map with the purpose of finding additional evidence  for an existing relationship between the variational properties of entropy expressions and the RG approach when applied to systems with scale invariance properties. We found that the non-extensive entropy function decreases monotonically as the RG transformation flows away from the fixed point where it attains its maximum value.

The ergodic hypothesis lies at the foundation of statistical mechanics implying that trajectories in phase space cover uniformly the entire pertinent regions. But is this hypothesis always correct? Already many years ago the answer to this question has been probed by the study of simple dynamical systems with only a few degrees of freedom . Besides their uncomplicated description these systems display extremely convoluted motion in phase space that is neither regular nor simply ergodic, and the mechanism by which ergodicity emerges in these and more complex deterministic systems has been effectively explored by studying the sensitivity to initial conditions and the associated Lyapunov exponents . The distinction between periodic and chaotic motion is signaled, respectively, by the long-time exponential approach or departure of trajectories with close initial positions. Here we have analyzed the borderline critical states in one-dimensional non-linear maps at which the exponential sensitivity law stops working and find that the universal dynamical behavior at these states precisely follow the predictions of the generalized non-extensive theory.

The currently developed non-extensive generalization ,  of BG statistical mechanics entails a rare examination of the domain of validity of this most fundamental base of physics. The suggested physical circumstances for which BG statistics fails to be applicable are thought to be associated to states that lack the full degree of chaotic irregular dynamics that probes phase space thoroughly and that is necessary for true equilibrium. Such anomalous states are signalled by the vanishing of the largest Lyapunov exponent and exhibit non-ergodicity or unusual power-law mixing. One particular case in point is that of critical states in nonlinear unimodal maps. Here impeded or incomplete mixing in phase space (the interval $`1x1`$) arises from the special ’tangency’ shape of the map at the pitchfork and tangent transitions that produces the monotonic trajectories given by Eq. (7). This has the effect of confining or expelling trajectories causing anomalous phase-space sampling, in contrast to the thorough coverage in generic states with $`\lambda _1>0`$. By construction the dynamics at the intermittency transitions, describe a purely nonextensive regime, since the map studied here, Eq. (3), does not consider access of trajectories to an adjacent or neighboring chaotic region, as in the setting of Refs. ,  or as in trajectories in conservative maps with weakly developed chaotic regions . Hence there is no reappearance of trajectories from chaotic regions that would cause the relaxation from the nonextensive regime with vanishing ordinary Lyapunov exponent to an extensive regime with a positive one at some crossover iteration time $`\tau `$.

The links we have exhibited through exact analytical results between the various properties of one-dimensional non-linear maps at critical conditions provide a clear-cut corroboration, with a universal attribute, in the RG sense, of the validity of non-extensive statistics at such states.

Acknowledgments. I would like to thank C. Tsallis and F. Baldovin for useful discussions and comments. I gratefully acknowledge the hospitality of the Centro Brasileiro de Pesquisas Fisicas where this work was carried out and the financial support given by the CNPq processo 300894/01-5 (Brazil). This work was also partially supported by CONACyT grant P-40530-F (Mexico).


format: text

text:
# 1 Introduction

## 1 Introduction

In the conventional hydrodynamic models for analyses of relativistic heavy-ion collisions, both chemical and thermal equilibrium are assumed and the chemical freeze-out temperature $`T_{\mathrm{ch}}`$ is the same as the thermal freeze-out temperature $`T_{\mathrm{th}}`$ which is so chosen as to reproduce the slope of $`p_t`$ (or $`m_t`$) spectra. For analyses of the $`p_t`$ slope at the RHIC energy by using the conventional hydrodynamics, see, e.g., Refs. \[1–5\]. On the other hand, various analyses based on thermal models indicate the following picture of time evolution of hot and dense hadronic matter : the system first undergoes the chemical freeze-out at $`T_{\mathrm{ch}}`$ where the observed particle ratios are fixed, continues to cool down and finally goes through the thermal freeze-out where the momentum distribution is fixed.111The observed momentum distribution can be affected by feeding from resonance decay after thermal freeze-out. We take account of this picture in a hydrodynamic model and discuss how the early chemical freeze-out affects the space-time evolution of fluids. We also analyse the particle spectrum at the RHIC energy. In the present work, we assume $`\mu _B=0`$ for simplicity. Although we neglect finite baryon effects, our hydrodynamic simulations are performed in a genuine three dimensional space. As a future work, we will include finite $`\mu _\mathrm{B}`$ in our model and analyse particle ratios and spectra at various collision energies .

## 2 Equations of State

We construct three models equation of state (EOS) to discuss the space-time evolution of fluids. These models describe the first order phase transition between the QGP phase and the hadron phase at $`T_\mathrm{c}=170`$ MeV. We suppose the QGP phase is composed of massless u, d, s quarks and gluons and that it is common to three models. For the hadron phase, we choose three different models EOS as follows.

The first model is a conventional resonance gas model in which complete chemical equilibrium is always assumed (model CE). We include hadrons up to the mass of $`\mathrm{\Delta }(1232)`$ for simplicity.

The second model is the simplest one which describes the picture of chemical freeze-out (model CFO). Below $`T_{\mathrm{ch}}`$, we assume the numbers of all hadrons $`N_i`$ included in the EOS are fixed and that the particle number densities obey $`_\mu n_i^\mu =0`$. We introduce a chemical potential $`\mu _i(T)`$ associated with each species so that $`N_i`$ becomes a conserved quantity. From the conservation of entropy $`_\mu s^\mu =0`$, the ratio of the particle number density to the entropy density obeys

$`{\displaystyle \frac{n_i(T<T_{\mathrm{ch}})}{s(T<T_{\mathrm{ch}})}}={\displaystyle \frac{n_i(T_{\mathrm{ch}})}{s(T_{\mathrm{ch}})}}`$ (1)
for all hadrons along the adiabatic path. Here we assume $`T_{\mathrm{ch}}=170`$ MeV, which is consistent with a recent analysis based on a thermal model . From Eq. (1), we obtain a chemical potential as a function of temperature for each hadron.

Fig. 1: (Left) Pressure as functions of energy density for three models (Right) Temperature as functions of energy density.

The third model represents a more realistic EOS than the second one. The observed particle numbers are always composed of the direct particles and the contribution from resonance decays, e.g., $`\overline{N}_\pi =N_\pi +_{i\pi }d_{i\pi }N_i`$. So some processes with large cross sections can be equilibrated even below $`T_{\mathrm{ch}}`$ as long as the equality $`\overline{n}_i(T<T_{\mathrm{ch}})/s(T<T_{\mathrm{ch}})=\overline{n}_i(T_{\mathrm{ch}})/s(T_{\mathrm{ch}})`$ is kept instead of Eq. (1). We assume $`\pi `$, $`K`$, $`\eta `$, $`N`$, $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ are “stable” particles and that all chemical potentials can be represented by chemical potentials associated with these stable particles, e.g., $`\mu _\rho =2\mu _\pi `$, $`\mu _K^{}=\mu _\pi +\mu _K`$, $`\mu _\mathrm{\Delta }=\mu _\pi +\mu _N`$ and so on . Thus the third model describes the partial chemical equilibrium (model PCE) even below $`T_{\mathrm{ch}}`$.

One can easily evaluate EOS’s for these models by using chemical potentials as obtained above. We represent pressure and temperature as functions of energy density for three models in Fig. 1. We find pressure as a function of energy density is similar to each other. On the other hand, temperature as a function of energy density in the models CFO or PCE is deviated from the model CE. See also Fig. 1 (right). Since the resonance population of the models CFO or PCE is larger than that of the model CE due to the chemical freeze-out, the energy density at a fixed temperature in the hadron phase is also large in those models. Conversely, temperature in the models CFO or PCE at a fixed energy density is smaller than in the model CE.

## 3 Radial and Elliptic Flow

By using the EOS’s in the previous section, we perform the hydrodynamic simulation with appropriate initial conditions and compare the hydrodynamic evolution with each other. We choose initial conditions so as to roughly reproduce pseudorapidity distribution in Au+Au 130$`A`$ GeV central (0-6 %) and semi-central (35-45 %) collisions obtained by the PHOBOS Collaboration . At the initial time $`\tau _0=0.6`$ fm/$`c`$, the energy density (temperature) for central collisions at the centre of the fluid, where it has a maximum value, is 30.8 GeV/fm3 (349.5 MeV). For further details of our hydrodynamic model, especially initial longitudinal and transverse profiles, see Ref. .

Fig. 2: (Left) Time evolution of temperature for three models. (Right) Thermal freeze-out temperature dependence of average radial flow for three models.

Figure 2 (left) shows the time evolution of temperature at the centre of the fluid for three models. From EOS’s represented in Fig. 1 (left), we expect the time evolution of energy density is similar to each other. However, when we transpose the energy density to the temperature by using EOS’s in Fig. 1 (right), the difference between three models appears in the time evolution of temperature. We find temperature in the models CFO or PCE cools down more rapidly than in the conventional model CE. We next show in Fig. 2 (right) the thermal freeze-out temperature dependence of average radial flow at vanishing space-time rapidity $`\eta _s=0`$ and compare our results with the ones from thermal models. Radial flows in the models CFO and PCE are considerably reduced. This is understood qualitatively by the same reason as the above discussion on the time evolution of temperature. These results seem to be consistent with the one obtained by Burward-Hoy , while these are completely inconsistent with the model II in Ref. . Since the results from thermal models are different from each other and we do not take into account the baryonic chemical potential, more detailed analyses based on both thermal and hydrodynamical models are needed to understand the thermal freeze-out.

Fig. 3: Thermal freeze-out temperature dependence of transverse momentum spectrum for negative charged hadrons (left) and $`v_2(p_t)`$ for charged hadrons (right). Solid, dashed and dotted lines correspond to, respectively, $`T_{\mathrm{th}}`$ = 140, 120 and 100 MeV

Finally we show in Fig. 3 the $`p_t`$ spectrum and $`v_2(p_t)`$ only for the model PCE. Left figure represents the transverse momentum distribution of negative charged hadrons for various thermal freeze-out temperatures. Experimental data is observed by the STAR Collaboration . Although the $`p_t`$ slope usually depends on $`T_{\mathrm{th}}`$ in the conventional EOS where the chemical equilibrium is assumed, our result for the model PCE is found to be almost independent of $`T_{\mathrm{th}}`$ due to the reduction of radial flow. Figure 3 (right) shows the transverse momentum dependence of $`v_2`$ for charged hadrons in minimum bias events. Our result with $`T_{\mathrm{th}}=140`$ MeV gives the best fit to the data obtained by the STAR Collaboration  below 1 GeV/$`c`$. The discrepancy of $`T_{\mathrm{th}}`$ between our result and the one in Ref.  may be due to the fact that the resonance decays and the chemical potential of pions are not taken into account in Ref. . We consistently reproduce experimental data of both $`p_t`$ spectrum and $`v_2(p_t)`$ below 1 GeV/$`c`$ by choosing $`T_{\mathrm{th}}=140`$ MeV, so contribution from hard processes can become dominated above 1 GeV/$`c`$. Note that the concavity of $`v_2(p_t)`$ in small $`p_t`$ region ($`p_t<0.3`$ GeV/$`c`$) is due to the resonance decay after thermal freeze-out .

## 4 Summary and Discussion

We have investigated the effect of early chemical freeze-out on radial and elliptic flow at the RHIC energy. We found the system cools down more rapidly than the conventional model due to the large population of resonances. We also found that the resultant radial flow is considerably reduced when we compare with each other at a fixed $`T_{\mathrm{th}}`$. Although the space-time evolution of temperature in the model PCE is considerably different from the conventional model, we consistently reproduced $`p_t`$ spectrum and $`v_2(p_t)`$ below 1 GeV/$`c`$ by choosing the thermal freeze-out temperature $`T_{\mathrm{th}}=140`$ MeV. This indicates the contribution from hard processes can become dominated above 1 GeV/$`c`$. It should be noted that the average radial flow results in $`<v_r>_{\eta _\mathrm{s}=0}0.35c`$.

Since the two pion correlation function is believed to be sensitive to the size and the lifetime of the fluid, it is interesting to see the effect of early chemical freeze-out on the HBT radii. Early chemical freeze-out in hydrodynamics will be one of the possible explanations to interpret the “HBT puzzle”  at the RHIC energy . In addition to the HBT radii, detailed analyses of particle ratio and spectra with finite $`\mu _\mathrm{B}`$ are also our future works.


text:
# Possible Origin of Lognormal Distributions in Gamma-Ray Bursts

## 1 INTRODUCTION

Among the statistical properties of the observed quantities, lognormal distributions are frequently seen in Gamma-Ray Bursts (GRBs). The lognormal distribution may be defined as the distribution of a random variable $`x`$ whose logarithm is normally distributed,

$`f(x)dx=\{\begin{array}{cc}\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left[\frac{(\mathrm{log}x\mu )^2}{2\sigma ^2}\right]d\mathrm{log}x,\hfill & \mathrm{if}x>0,\hfill \\ 0,\hfill & \mathrm{if}x0,\hfill \end{array}`$ (3)
where $`f(x)`$ is the probability density function for $`x`$, and $`\mu `$ and $`\sigma ^2`$ are the sample mean and the variance of $`\mathrm{log}x`$ (e.g., Crow & Shimizu 1988). This distribution is unimodal and positively skew. McBreen, Hurley, Long & Metcalfe (1994) pointed out that the total duration of the long and short bursts and the time interval between pulses are consistent with the lognormal distributions. Li & Fenimore (1996) showed that the pulse fluence and the pulse interval distributions within each burst are consistent with lognormal distributions. Nakar & Piran (2001) found that the pulse duration also have a lognormal distribution. The break energy distribution is also lognormal (Preece et al. 2000; see below).

Lloyd, Fryer & Ramirez-Ruiz (2001) suggested that $``$ 10% of GRBs might have the redshift larger than 6 so that the redshift distribution might be wide. Therefore it is quite strange that the observed break energy distribution and the duration distribution are lognormal since it does not seem that the observed lognormal distribution reflects the redshift distribution of the GRBs. Even if the break energy distribution is lognormal at the source, the observed break energy should be smaller than the intrinsic one by a factor of $`(1+z)`$ while the observed duration should be longer than the intrinsic one by a factor of $`(1+z)`$. These factors change by order of unity between $`z=0`$ and $`z6`$.

In this Letter, we will consider a possible origin of the observed lognormal distributions in GRBs from the viewpoint of the central limit theorem.111The lognormal distribution of the time interval between pulses and the pulse fluence might be reproduced by a fine tuning of model parameters (Spada, Panaitescu & M$`\stackrel{´}{\mathrm{e}}`$sz$`\stackrel{´}{\mathrm{a}}`$ros 2000).

## 2 BREAK ENERGY DISTRIBUTIONS

Figure 1 shows the histogram of the break energy $`E_b`$ taken from the electronic edition of Preece et al. (2000).222 It is not known whether the paucity of the soft and hard bursts is real or not, because harder bursts have fewer photons (Cohen, Piran & Narayan 1998; Lloyd & Petrosian 1999; but see Brainerd et al. 1999) and there may exist relatively many soft bursts with low luminosities, so called X-ray rich GRBs (or X-ray flashes or Fast X-ray transients) (Strohmayer et al. 1998; Heise et al. 2001; Kippen et al. 2001). Here we assume that the selection effect is small. The $`\chi ^2`$ test of all data gives the probability of $`1.4\times 10^{185}`$ (the reduced $`\chi ^2`$ is 16.5 with 66 degrees of freedom) that the data was taken from the lognormal distribution. Therefore the null hypothesis that the break energy distribution is lognormal fails. However, if we exclude the data in the high and low energy ends, the fit becomes good as shown in Figure 1. The $`\chi ^2`$ test of the data between $`70.8`$ keV and $`708`$ keV gives the probability of 0.497 (the reduced $`\chi ^2`$ is 0.963 with 17 degrees of freedom) that the data was taken from the lognormal distribution, with $`\mu =2.38\pm 0.004`$ ($`E_b238`$ keV) and $`\sigma =0.240\pm 0.004`$ ($`1\sigma `$ width is between $`137`$ keV and $`413`$ keV). The improvement of the lognormal fit to the break energy distribution excluding the high and low energy ends may suggest that the soft and hard bursts originate from a different class of GRBs or emission mechanisms. Anyway, hereafter we will assume that the observed break energy distribution is lognormal.

Now let us assume that the intrinsic break energy distribution is lognormal. Next we numerically calculate the observed break energy distribution assuming that the redshift distribution has the form,

$`f(z)dz=\{\begin{array}{cc}A(1+z)^{a1}(1+z_0)^{b1}dz,\hfill & \mathrm{if}0<z<z_0,\hfill \\ A(1+z_0)^{a1}(1+z)^{b1}dz,\hfill & \mathrm{if}z_0<z,\hfill \end{array}`$ (6)
where $`b<0`$. This redshift distribution rises as $`(1+z)^{a1}`$ to a redshift of $`z_0`$ and then declines as $`(1+z)^{b1}`$, and it is similar to that in Figure 8 of Lloyd, Fryer & Ramirez-Ruiz (2001). To mimic Figure 1, we generate 155 bursts with each burst having 35 spectra for each realization. We take the mean and the variance of the intrinsic break energy distribution so that the mean and the variance of the observed break energy are close to that in Figure 1. For each realization, we make a $`\chi ^2`$ fitting to obtain the probability that the data is taken from the lognormal distribution. As in Figure 1, we use the data between 70.8 keV and 708 keV for the $`\chi ^2`$ test. From this simulation, we have found that the average probability can reach $`0.5`$ even when the variance of the redshift distribution is comparable to the observed one $`\sigma _z^2(0.240)^2`$, although it is slightly stronger condition to preserve the lognormal form than $`\sigma _z^2<(0.240)^2`$. At first glance it is strange that the simulations do not reflect the redshift distribution contrary to the argument in Section 1. In the next section we will show that it is not strange but natural due to the central limit theorem.

Figure 2 shows the histogram of the observed break energy for one experimental realization with $`(z_0,a,b)=(3,3,3)`$. The $`\chi ^2`$ test of the data between $`70.8`$ keV and $`708`$ keV gives the probability of $`0.128`$ (the reduced $`\chi ^2`$ is $`1.39`$ with $`17`$ degrees of freedom) that the data was taken from the lognormal distribution, with $`\mu =2.37\pm 0.004`$ ($`E_b235`$ keV) and $`\sigma =0.256\pm 0.004`$ ($`1\sigma `$ width is between $`130`$ keV and $`423`$ keV). It is interesting to note the excess of soft bursts relative to the lognormal fit as in Figure 1. The average redshifts of these soft bursts are relatively high.

## 3 LOGNORMAL DISTRIBUTIONS

The standard model of the GRB emission is the optically thin synchrotron shock model (e.g., Piran 1999). A similar discussion in the following will be applied to the inverse Compton model. Let us consider a slow (rapid) shell with a Lorentz factor $`\gamma _s`$ ($`\gamma _r`$), a mass $`m_s`$ ($`m_r`$) and a width $`l_s`$ ($`l_r`$). When a separation between two shells is $`L`$, the collision takes place at a radius $`R_s2L\gamma _s^2`$. At the collision, the forward and the reverse shock are formed. Here we consider the reverse shock propagating into the rapid shell. The discussion for the forward shock is similar. We assume that a fraction of electrons $`\zeta _e`$ is accelerated in the shock to a power law distribution of Lorentz factor $`\gamma _e`$, $`N(\gamma _e)d\gamma _e\gamma _e^pd\gamma _e`$ for $`\gamma _e\gamma _{min}[(p2)/(p1)](ϵ_eu^{}/\zeta _en^{}m_ec^2)`$, where $`n^{}`$ and $`u^{}`$ are the number density and the internal energy density in the local frame, respectively, $`p2`$, and we assume that a fraction $`ϵ_e`$ of the internal energy goes into the electrons. We also assume that a fraction $`ϵ_B`$ of the internal energy goes into the magnetic field, $`B^2=8\pi ϵ_Bu^{}`$. The local frame quantities, $`u^{}`$ and $`n^{}`$, can be calculated using the shock jump conditions (Blandford & Mckee 1976; Sari & Piran 1995). We assume that the unshocked shells are cold and the shocked shells are extremely hot. If the Lorentz factor of the shocked region is $`\gamma `$, the relative Lorentz factor of the unshocked and the shocked region is given by $`\gamma _{rel}(\gamma _r/\gamma +\gamma /\gamma _r)/2\gamma _r/2\gamma `$, so that $`u^{}=(\gamma _{rel}1)n^{}m_pc^2\gamma _{rel}n^{}m_pc^2`$. The number densities of the unshocked and the shocked region are given by $`n_r^{}=m_r/4\pi m_pR_s^2l_r\gamma _rm_r/16\pi m_pL^2\gamma _s^2l_r\gamma _r`$ and $`n^{}=(4\gamma _{rel}+3)n_r^{}4\gamma _{rel}n_r^{}`$, respectively. Thus, the characteristic synchrotron energy is given by

$`E_b={\displaystyle \frac{\mathrm{}q_eB\gamma \gamma _{min}^2}{m_ec(1+z)}}260\left[{\displaystyle \frac{p2}{p1}}\right]^2ϵ_e^2ϵ_B^{1/2}\zeta _e^2m_{r,28}^{1/2}L_{10}^1l_{r,9}^{1/2}\gamma _{s,2}^2\gamma _{r,2}^{1/2}\gamma _{rel}^2(1+z)^1\mathrm{keV},`$ (7)
where we assume that the source is at a redshift $`z`$. Note that the relative Lorentz factor of the unshocked and shocked region $`\gamma _{rel}`$ depends on the relative Lorentz factor of the rapid and slow shell $`\gamma _{rs}(\gamma _r/\gamma _s+\gamma _s/\gamma _r)/2\gamma _r/2\gamma _s`$ and the ratio between the number densities in these shells $`fn_s^{}/n_r^{}=m_sl_r\gamma _r/m_rl_s\gamma _s`$ (Sari & Piran 1995). For the ultrarelativistic shock case $`\gamma _{rs}^2f`$, $`\gamma _{rel}=\gamma _{rs}^{1/2}f^{1/4}/\sqrt{2}=(m_sl_r\gamma _r^3/m_rl_s\gamma _s^3)^{1/4}/2`$.

Equation (7) shows that the break energy is written in the form of a product of many variables. For such a variable made from the product of many variables, the lognormal distribution may have a very simple origin, that is, the central limit theorem (Crow & Shimizu 1988; Montroll & Shlesinger 1982). Let a variable $`q`$ be written in the form of a product of variables,

$`q=x_1x_2\mathrm{}x_n.`$ (8)
Then,

$`\mathrm{log}q=\mathrm{log}x_1+\mathrm{log}x_2+\mathrm{}+\mathrm{log}x_n.`$ (9)
When the individual distributions of $`\mathrm{log}x_i`$ satisfy certain weak conditions that include the existence of second moments, the central limit theorem is applicable to the variable $`\mathrm{log}q`$, so that the distribution function of $`\mathrm{log}q`$ tends to the normal distribution as $`n`$ tends to infinity.

As an example, we numerically generated random variables $`x_i`$ $`(i=1,2,\mathrm{})`$ whose logarithms are uniformly distributed between 0 and 1. Figure 3 shows the histogram of the product of these three variables, $`q=x_1x_2x_3`$, for $`10^4`$ experimental realizations. The distribution of $`q`$ agrees with the lognormal distribution quite well. It is surprising that the $`\chi ^2`$ test gives the probability of $`0.483`$ (the reduced $`\chi ^2`$ is 1.00 with 278 degrees of freedom) that the distribution of $`q`$ is taken from the lognormal distribution. This example shows that the lognormal distributions may be achieved by a relatively small number of variables (Yonetoku & Murakami 2001). Note that when the number of the variables is two, i.e., $`q=x_1x_2`$, the probability that the distribution is taken from the lognormal distribution was only $`1.64\times 10^5`$, so that a product of only one more variable may make a distribution lognormal.

Therefore the lognormal distribution of the break energy may be a natural result from the central limit theorem. We may say, “Astrophysically, not $`n=\mathrm{}`$ but $`n=3`$ gives the lognormal distribution !!”. The effect of the redshift is just to add one variable in equation (7) if the redshift in the observed data is randomly chosen.

## 4 PULSE FLUENCE/DURATION/INTERVAL DISTRIBUTIONS

Let us consider the lognormal distributions in other quantities related to GRBs. When the rapid shell catches up the slow one in the internal shock, using the conservation of the energy and the momentum, the Lorentz factor of the merged shell $`\gamma _m`$ and the internal energy $`E_{int}`$ produced by the collision are given by $`\gamma _m[(m_r\gamma _r+m_s\gamma _s)/(m_r/\gamma _r+m_s/\gamma _s)]^{1/2}`$ and $`E_{int}=m_r\gamma _r+m_s\gamma _s(m_r+m_s)\gamma _m`$, respectively (e.g., Piran 1999). If we assume that a fraction $`ϵ_e`$ of the internal energy goes into the electrons and a fraction $`ϵ_w`$ of the energy radiated by the electrons is within the gamma-ray band, the observed energy is given by

$`E_{obs}=ϵ_wϵ_eE_{int}(1+z)^1ϵ_wϵ_em_r\gamma _r(1+z)^1.`$ (10)
Equation (10) shows that the observed energy, which is proportional to the pulse fluence, is written in the form of a product of five variables,333 Although the fraction $`ϵ_e`$ of the internal energy that goes into electrons may be a fundamental constant, there is still a large dispersion of order unity in $`ϵ_e`$ deduced from the the afterglow observations (Panaitescu & Kumar 2001). $`ϵ_w`$, $`ϵ_e`$, $`m_r`$, $`[\gamma _r\gamma _m+(m_s/m_r)(\gamma _s\gamma _m)]\gamma _r`$ and $`(1+z)^1`$. Therefore the lognormal distribution of the pulse fluence may be a natural result from the central limit theorem.

The pulse duration is determined by three time scales: the hydrodynamic time scale, the cooling time scale, and the angular spreading time scale (Kobayashi, Piran & Sari 1997; Katz 1997; Fenimore, Madras & Sergei 1996). The cooling time scale is usually much shorter than the other two time scales in the internal shocks (Sari, Narayan & Piran 1996). The hydrodynamic time scale $`l/c`$ and the angular spreading time scale determine the rise and the decay time of the pulse, respectively. Since most observed pulses rise more quickly than they decay (Norris et al. 1996), we assume that the pulse duration is mainly determined by the angular spreading time $`R_s/2c\gamma _m^2`$. Then, the pulse duration $`\delta t`$ is given by

$$\delta t(L/c)(\gamma _s^2/\gamma _m^2)(1+z).$$
(11)
On the other hand, the interval between pulses $`\mathrm{\Delta }t`$ is determined by the separation between shells,

$`\mathrm{\Delta }t(L/c)(1+z),`$ (12)
since all shells are moving towards us with almost the speed of light (Kobayashi, Piran & Sari 1997; Nakar & Piran 2001).

Equation (12) shows that the pulse interval $`\mathrm{\Delta }t`$ reflects the separation between shells $`L`$, while equation (11) shows that the pulse duration $`\delta t`$ is multiplied by one more factor $`(\gamma _s/\gamma _m)^2`$ other than $`(L/c)(1+z)`$. Therefore, if we consider that the distribution of a product of variables tends to the lognormal distribution as the number of the multiplied variables increases, the distribution of the pulse duration $`\delta t`$ may be closer to the lognormal distribution than that of the pulse interval $`\mathrm{\Delta }t`$. In fact, Nakar & Piran (2001) argued that the pulse duration $`\delta t`$ has the lognormal distribution while the pulse interval $`\mathrm{\Delta }t`$ does not, as noticed by Li & Fenimore (1996). The pulse interval has an excess of long intervals relative to the lognormal distribution. This may suggest the existence of a different distribution, i.e., quiescent times (long periods with no activity) (Nakar & Piran 2001; Ramirez-Ruiz & Merloni 2001).444 The correlation between the pulse interval $`\delta t`$ and the pulse duration $`\mathrm{\Delta }t`$ is broken when the separation between shells is too large to collide each other before the external shock. However, almost all quiescent times in Nakar & Piran (2001) are smaller than the limit $`\delta t100\mathrm{sec}E_{52}^{1/3}n_1^{1/3}\gamma _{100}^{8/3}`$ so that shells may collide. But the central limit theorem may be also responsible for the lognormal distribution of the pulse duration $`\delta t`$.

## 5 DISCUSSIONS

We considered the possible origin of the lognormal distributions in the break energy, the pulse fluence and the pulse duration as a result of the central limit theorem. Astrophysically the lognormal distribution may be achieved by a product of only a few variables. The effect of the redshift is just to add one variable to the product so that the redshift distribution is hidden.

We have no idea about the origin of the lognormal distributions in the pulse interval555 Note that the distribution of the pulse interval is lognormal if we exclude the quiescent times (Nakar & Piran 2001). $`\mathrm{\Delta }t`$ and the total duration $`\mathrm{\Delta }T`$. However the viewing angle may be one factor to be multiplied to the pulse interval $`\mathrm{\Delta }t`$. Recently we suggested that the luminosity-lag relation could be explained by the variation in the viewing angle $`\theta _v`$ from the axis of the jet (Ioka & Nakamura 2001; Nakamura 2000). The duration of the pulse from the jet also depends on the viewing angle, and according to Figure 2 of Ioka & Nakamura (2001) we have

$`\mathrm{\Delta }t(L/c)(1+z)(1+\gamma ^2\theta _v^2),`$ (13)
when $`\theta _v\mathrm{\Delta }\theta `$ where $`\gamma `$ is the Lorentz factor of the jet and $`\mathrm{\Delta }\theta `$ is the opening half-angle of the jet. The multiplied factor $`(1+\gamma ^2\theta _v^2)`$ may be responsible for the lognormal distribution of the pulse interval $`\mathrm{\Delta }t`$. Note that the total duration $`\mathrm{\Delta }T`$ is equal to the lifetime of the central engine and thus does not depend on $`\theta _v`$.

We are grateful T. Murakami, D. Yonetoku and S. Kobayashi for useful comments. We thank the Yukawa Institute for Theoretical Physics at Kyoto University, where this work was initiated during the YITP-W-01-06 on “GRB2001”. This work was supported in part by Grant-in-Aid for Scientific Research Fellowship of the Japanese Ministry of Education, Science, Sports and Culture, No.00660 (KI) and by Grant-in-Aid of Scientific Research of the Ministry of Education, Culture, and Sports, No.11640274 (TN) and 09NP0801 (TN).


text:
# THE MILLENNIUM ARECIBO 21-CM ABSORPTION LINE SURVEY: EXTOLLING THE ANALYSIS PROCEDURES

## 1 OBTAINING THE OPACITY AND EXPECTED HI PROFILES

One can derive spin temperature $`T_s`$ from emission-absorption observations because the opacity $`\tau (\nu )`$ depends on $`T_s`$ differently than the brightness temperature $`T_B(\nu )`$ does. Naively, one makes one ON-source and one or more OFF-source observations. One combines the off-source observations to produce the expected profile $`T_{exp}(\nu )`$, which is the $`T_B(\nu )`$ that would be seen at the source position if the source were magically turned off. Both $`\tau (\nu )`$ and $`T_{exp}(\nu )`$ are affected by HI structure on the sky; these uncertainties are a major item.

For each source we observed a grid of 16 OFF-source positions, plus the ON-source position. This allows us to expand the observed brightness temperature $`T_B(\nu )`$ in a second-order two-dimensional Taylor series around the source position and to derive the following quantities, plus their uncertainties:

1. The opacity profile $`\tau (\nu )`$.
2. The expected profile $`T_{exp}(\nu )`$.
3. The first-derivative profiles $`\frac{T_B(\nu )}{\alpha }`$ and $`\frac{T_B(\nu )}{\delta }`$ ($`\alpha `$ is right ascension and $`\delta `$ is declination).
4. The second-derivative profiles $`\frac{^2T_B(\nu )}{\alpha ^2}`$, $`\frac{^2T_B(\nu )}{\alpha \delta }`$, and $`\frac{^2T_B(\nu )}{\delta ^2}`$.

Figure 1 exhibits results for 3C75. In the top panel, the solid line profile extending above zero is $`T_{exp}(\nu )`$. The dashed curves show the observed emission from each individual WNM Gaussian, including absorption by the CNM, and the heavy dashed curve shows the totality of WNM fitted Gaussians. The dotted curves show the intrinsic emission from each individual CNM component, and the heavy dots the observed emission of all CNM components including absorption from the CNM components that lie in front. The solid line near mid-profile-height is the residuals of the data from the fit. The solid line profile extending below zero is the (negative of the) uncertainties in $`T_{exp}(\nu )`$. The annotation shows the properties of the WNM Gaussian components.

In the bottom panel, the solid line profile extending above zero is the opacity profile $`[1e^{\tau (\nu )}]`$. The dotted curves show each individual CNM opacity Gaussian and the heavy dots show the opacity sum of all the CNM Gaussians ($`[1e^{{\scriptscriptstyle \tau (\nu )}}]`$. The solid line near mid-profile-height is the residuals of the data from the fit; the dashed line shows the contribution to opacity from each WNM component assuming that its temperature is $`T_{kmax}`$. The solid line profile extending below zero is the (negative of the) uncertainty in the opacity profile. The annotation shows the properties of the CNM Gaussian components.

## 2 THE GAUSSIAN ANALYSIS PROCEDURE

We choose 3C75 as an example because both the opacity and expected profiles are fairly complicated, yet the fits are unambiguous. There are many sources with less complicated profiles that are fit with fewer components. Of the 79 sources, 58 have two or WNM fewer components.

In the fitting process we first fit the opacity profile with the fewest number of CNM Gaussian components required to make the residuals, on the middle line, comparable with the uncertainties, the negative-going profile at the bottom. For 3C75 this requires three CNM Gaussians, one for each of the major peaks. They are unambiguous: one could certainly not fit fewer. One could include one or more additional components, but this is not required in view of the uncertainties. In the main papers we show that this arbitrariness does not affect derived spin temperatures.

We next perform a least squares fit on $`T_{exp}(\nu )`$. The parameters include one spin temperature for each CNM component; emission from an arbitrary number of WNM components; the ordering of the CNM components along the line of sight; and the fraction of WNM that lies behind the CNM components. These last two are required because the CNM components absorb emission from WNM that lies behind them, and each CNM component absorbs emission from the CNM components that lie behind it.

It is important to note that the CNM emission, shown by the dotted lines on the top panel, is almost always narrower than the total emission. This highlights the requirement for WNM emission that is both broader than the CNM and, also, produces no detectable absorption. This width difference means that:

* The spin temperatures depend mainly on the $`T_{exp}(\nu )`$ intensity within the velocity ranges covered by the dotted-line CNM opacity Gaussians.
* The properties of the WNM Gaussians are determined mainly by the $`T_{exp}(\nu )`$ profile shape lying outside the velocity ranges covered by the dotted-line CNM opacity Gaussians.

In this case of 3C75, we derive upper limits to WNM kinetic temperature $`T_{kmax}=3720`$ and 2095 K; both are in the thermally unstable range. Looking at the $`T_{exp}(\nu )`$ profile, it is clear that there is no way to obtain wider WNM Gaussians; the WNM widths are constrained by the portions of the emission profile that lie outside the dotted CNM emission. One can obtain narrower WNM Gaussians in two ways:

1. Include a very broad, weak Gaussian to sit underneath the main emission profile. This would eliminate the systematic residuals in the profile wings and produce a better fit. We included such a component for many sources; indeed, for most sources with two WNM components, one of the components is usually such a broad, weak pedestal.
2. Fit a larger number of narrow WNM Gaussians to $`T_{exp}(\nu )`$. One can often replace a single-Gaussian fit by more narrower Gaussians. However, this violates the philosophy of using the minimum number of free parameters. Moreover, with more narrower WNM components, the WNM temperatures would decrease and their associated opacities would increase; given the uncertainties in $`\tau (\nu )`$ there is some room for this, but not much.

In this case, then, we conclude that WNM lying in the unstable temperature range is undeniable. Similar conclusions apply in most cases. Some sources have less unambiguous fitting results, but we always attempt to be “reasonable” and this philosophy leads to the statistics that we quote in the abstract and introduction.

This work was supported in part by NSF grants AST-9530590 and AST-0097417, and by the NAIC.


text:
# Strong-field terahertz-optical mixing in excitons

## Abstract

Driving a double-quantum-well excitonic intersubband resonance with a terahertz (THz) electric field of frequency $`\omega _{THz}`$ generated terahertz optical sidebands $`\omega =\omega _{THz}+\omega _{NIR}`$ on a weak NIR probe. At high THz intensities, the intersubband dipole energy which coupled two excitons was comparable to the THz photon energy. In this strong-field regime the sideband intensity displayed a non-monotonic dependence on the THz field strength. The oscillating refractive index which gives rise to the sidebands may be understood by the formation of Floquet states, which oscillate with the same periodicity as the driving THz field.

Electromagnetically-induced coherent quantum effects in semiconductor systems include Rabi oscillations in donors and quantum dots, and electromagnetically-induced transparency in quantum wells. In these effects, to induce coherence the driving field must overcome dephasing processes. Thus, the Rabi frequency $`\mu E`$ must be greater than the dephasing rates $`\gamma `$: $`\mu E>\mathrm{}\gamma `$, where $`\mu `$ is the dipole moment and $`E`$ is the electric field strength.

If the field strength is increased still further, a strong-field regime is encountered as the Rabi energy becomes comparable to the photon energy: $`\mu E\mathrm{}\omega `$. Previous studies into this regime have looked at atoms in microwave cavities and strong laser fields. Since atomic excited states are near the continuum, strong-field effects manifest themselves as nonmonotonic multi-photon ionization rates as a function of field strength.

However, the atoms are ionized, meaning the bound states do not exist after the microwave or laser field is turned on. Multi-photon ionization is an inescapable result of driving excited atomic levels which get closer and closer together.

There is a rich body of theoretical predictions on the effect of strong fields on the bound states of a quantum system. In quantum wells, upper states get further and further apart. Thus the strong-field condition can be satisfied by resonantly driving the lowest subbands of a quantum well at THz frequencies. Yet the quantum well is far deeper than the Rabi energy, so the bound states still survive.

We describe experiments in which a near-infrared (NIR) probe laser beam is mixed with a terahertz (THz) pump beam in a gated, asymmetric double quantum well (DQW). The THz field couples to an excitonic intersubband excitation while the NIR field couples to an excitonic interband excitation. Applying a DC voltage to the gates brings the intersubband transition into resonance with the THz field, and the interband transition into resonance with the NIR field. When these resonance conditions are met the NIR probe is modulated resulting in the emission of optical sidebands which appear at frequencies $`\omega _{sideband}=\omega _{NIR}+n\omega _{THz}`$ where $`\omega _{NIR}`$ ($`\omega _{THz}`$) is the frequency of the NIR (THz) beam and $`n=\pm 1,2,\mathrm{}`$. Strong-field effects manifest themselves as non-monotonic sideband intensities as a function of THz field strength.

The sample consisted of an active region, gates, and a distributed Bragg reflector (DBR). The active region consisted 5 periods of DQW, each consisting of a 120 Å GaAs QW and a 100 Å GaAs QW separated by a 25 Å Al0.2Ga0.8As tunnel barrier. Each period was separated by a 300 Å Al0.3Ga0.7As barrier. The active region was sandwiched between two gates, each consisting of a Si delta-doped 70 Å QW with carrier density $`1\times 10^{12}`$ cm-2. The gates were separated from the active region by 3000 Å Al0.3Ga0.7As barriers. Since the gate QWs are much narrower than the active DQWs, the gate QWs are transparent to both the THz and NIR beams. The gated DQWs were grown on top of a DBR which consisted of 15 periods of 689 Å AlAs and 606 Å Al0.3Ga0.7As. It had a low-temperature passband nearly centered on the low-temperature bandgap of the DQW, making it about 95% reflective for the NIR probe beam and sidebands.

We etched a mesa and annealed NiGeAu ohmic contacts to the gate QWs. The sample was then cleaved into a 1 mm wide strip, 8 mm long. We then cleaved a 400 $`\mu `$m-thick wafer of crystal sapphire substrate material into a $`1\times 5`$ mm strip, oriented with the optically slow axis along the long dimension. The sapphire was mechanically pressed against the surface of the sample with a beryllium copper clip as shown in Fig. 1. At cryogenic temperatures sapphire is index-matched to GaAs and transparent at THz wavelengths. At the same time, sapphire is transparent to the NIR probe. Pressing the sapphire against the sample forms a rectangular dielectric waveguide with half of the waveguide defined by the sample substrate and the other half defined by the sapphire. The epilayer containing the active region lies in the middle of the heavily over-moded waveguide. This results in much higher fields at the active region than would be possible without the sapphire, where the active region would be at the edge of a dielectric waveguide.

The sample was cooled to 21 K in a closed-cycle He cryostat. The 1 kW THz beam from a free-electron laser is polarized in the growth direction, propagated in the QW plane and focused by a 90o off-axis parabolic mirror (F/1) into the edge of the waveguide. The maximum THz electric field strength at the focus is estimated to be between $`520`$ kV/cm, where our uncertainty comes from uncertainty of the exact size of the mode in the waveguide. The THz power can be continuously varied by a pair of wire-grid polarizers.

NIR light from a continuous-wave Ti:sapphire laser is chopped by an acousto-optic modulator into 25 $`\mu `$s pulses which overlap the 5 $`\mu `$s FEL pulses at maximum repetition rate of 1.5 Hz. The vertically-polarized NIR beam propagates normal to the THz beam, and was focused (F/10) at a power density of $``$ 5 W/cm2 to the same small interaction volume in the sample. The reflected beam, sidebands, and photoluminescence (PL) are analyzed by a second polarizer, dispersed by a 0.85 m double-monochromator, and detected by a photomultiplier tube.

Previous work  showed that sideband generation is enhanced when the NIR field resonantly couples the vacuum state $`|0`$ with an exciton $`|1`$, and the THz field resonantly couples with the intersubband transition between excitons $`|1`$ and $`|2`$. The exciton states can be labelled as EμHHνX, indicating an exciton consisting of an electron in conduction subband $`\mu `$ and a heavy hole in valence subband $`\nu `$. Changing the DC electric field tunes the excitonic intersubband transition into resonance with the THz field, as shown in Fig. 2. The sideband resonance we focused on is E2HH2X-E2HH3X, a peak assignment made by comparing low-field sideband spectroscopy with a nonlinear susceptibility calculation for excitons.

Fig. 3 shows the dependence on THz field strength of sideband generation at various THz frequencies. Each curve was taken at a gate bias near the E2HH2X-E2HH3X resonance. The most striking feature is the non-monotonic behavior in the strong-field regime when the Rabi energy is comparable to the photon energy. Clearly at lower frequencies the power dependence rolls over at a lower field strength.

By sitting at the peak NIR frequency for E2HH2X-E2HH3X and varying the DC gate voltage we took THz power dependence scans at various THz detunings $`E_2E_1\mathrm{}\omega _{THz}`$ where $`E_2`$ ($`E_1`$) is the energy of the upper (lower) exciton state. The results are shown in Fig. 4. As the detuning is varied, so does the shape of the THz field dependence.

The power dependence cannot be explained by a nonlinear susceptibility, which is inherently a low-field theory because it relies on a Taylor expansion about the field strength . Saturation effects in a nonlinear susceptibility can only come about when contributions from virtual transitions initiated from excited states $`|1`$ or $`|2`$ are comparable to those initiated from the vacuum $`|0`$. However the populations of excited states $`|1`$ or $`|2`$ are never significant compared with $`|0`$ in our undoped sample and weak NIR beam. Therefore a nonlinear susceptibility can only predict a linear (quadratic) dependence of sideband intensity on THz power (field strength).

The following simple, phenomenological, three-state model captures qualitative features of the dependence of sidebands on THz power. The Hamiltonian $`H_0`$ has eigenstates $`\varphi _0(x,z)`$, $`\varphi _1(x,z)`$, and $`\varphi _2(x,z)`$ with eigenenergies $`E_0=0`$, $`E_1`$, and $`E_2`$, respectively. Near the experimental resonance, the effect of the DC voltage is to tune $`\varphi _2`$ while the other $`\varphi _1`$ remains relatively unchanged. The two exciton states $`\varphi _1(x,z)`$ and $`\varphi _2(x,z)`$ are coupled by the $`z`$-dipole operator for the DC and THz fields. The two upper states $`\varphi _1(x,z)`$ and $`\varphi _2(x,z)`$ are coupled to the ground state by the $`x`$-dipole operator for the NIR field.

Therefore, the x and z matrices are

$`z=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& z_{12}\\ 0& z_{21}& z_{22}\end{array}\right)x=\left(\begin{array}{ccc}0& x_{01}& x_{02}\\ x_{10}& 0& 0\\ x_{20}& 0& 0\end{array}\right)`$ (7)
All the nonzero terms $`x_{\alpha \beta }=\varphi _\alpha |x|\varphi _\beta `$ and $`z_{\alpha \beta }=\varphi _\alpha |z|\varphi _\beta `$ are set equal to unity. The Hamiltonian is

$$H=H_0+zE_\omega \mathrm{cos}\omega t+x\lambda E_\mathrm{\Omega }\mathrm{cos}\mathrm{\Omega }t$$

where $`E_\omega `$ represents the strong THz electric field strength and $`E_\mathrm{\Omega }`$ represents the weak NIR electric field. We solved the Hamiltonian nonperturbatively within a Floquet formalism for $`\lambda =0`$, while the weak probe was added later using time-dependent perturbation theory.

The solutions to the Schrodinger equation for the time-periodic Hamiltonian have the form

$`\phi _i(z,t)=e^{iϵ_it/\mathrm{}}u_i(z,t),u(z,t)=u(z,t+{\displaystyle \frac{2\pi }{\omega }})`$ (8)
where $`u_i(z,t)`$ has the same periodicity as the driving frequency and and $`ϵ_i`$ is called the quasienergy. The states $`\phi _i(z,t)`$ are called Floquet states and are mathematically analogous to Bloch states for a spatially-periodic Hamiltonian. Meanwhile the quasienergies are analogous to the crystal momenta of Bloch theory.

The spatial-dependence of the wavefunction $`\phi _i(z,t)`$ can be expanded in terms of the original eigenstates $`\varphi _\alpha (z)`$, and the time-dependence can be Fourier-expanded

$$\phi _i(z,t)=e^{iϵ_it/\mathrm{}}\underset{\alpha ,n}{}c_{\alpha ,n}^ie^{in\omega t}\varphi _\alpha (z)$$
(9)
The coefficients $`c_{\alpha ,n}^i`$ are what will determine the part of the index of refraction which will oscillate at particular multiples of $`\omega `$. Solving for these coefficients is the key to understanding the power dependence of the sideband. These Floquet coefficients are closely related to the photon-assisted tunnelling sidebands which appear in the irradiated current-voltage curves in superconducting weak-link junctions and coupled quantum wells.

The time-dependent Schrodinger equation can be cast in terms of the Floquet operator $`F=Hi\mathrm{}\frac{}{t}`$, which can be written as a matrix of the form

$$F=\left(\begin{array}{ccccc}\mathrm{}& [\mu ]& 0& 0& 0\\ [\mu ]& [H_0]\mathrm{}\omega & [\mu ]& 0& 0\\ 0& [\mu ]& [H_0]& [\mu ]& 0\\ 0& 0& [\mu ]& [H_0]+\mathrm{}\omega & [\mu ]\\ 0& 0& 0& [\mu ]& \mathrm{}\end{array}\right)$$
(10)
where $`[H_0]`$ and $`[\mu ]`$ are matrix representations of the operators $`H_0`$ and $`\mu `$ in the $`|\varphi _\alpha `$ basis. In practice the Floquet matrix must be truncated up to $`\pm N`$ photons. $`N`$ can be made arbitrarily large for an arbitrarily precise result at high field strengths. The solution to the Schrodinger equation for a time-periodic Hamiltonian reduces to finding the eigenvalues and eigenvectors of the Floquet matrix (10). Other than the $`\pm N`$ truncation the solution involves no perturbation or rotating-wave approximations. Here we used $`N=16`$ photons.

The Floquet solutions are labeled $`\phi _i(x,z,t)`$ with $`i=0,1,2`$, and expanded as in (9). Their expansions are explicitly expressed for clarity:

$`\phi _0(x,z,t)`$ $`=`$ $`\varphi _0(x,z)`$
$`\phi _1(x,z,t)`$ $`=`$ $`e^{iϵ_1t/\mathrm{}}{\displaystyle \underset{n}{}}e^{in\omega t}[c_{1,n}^1\varphi _1(x,z)+c_{2,n}^1\varphi _2(x,z)]`$
$`\phi _2(x,z,t)`$ $`=`$ $`e^{iϵ_2t/\mathrm{}}{\displaystyle \underset{n}{}}e^{in\omega t}[c_{1,n}^2\varphi _1(x,z)+c_{2,n}^2\varphi _2(x,z)]`$
The ground state $`\varphi _0(x,z)`$ is not coupled by the the strong field to the 2-d subspace spanned by $`|\varphi _1`$ and $`|\varphi _2`$. Therefore the coefficients for the Floquet state $`\phi _0(x,z,t)`$ vanish except $`c_{0,0}^0=1`$. In other words, the ground vacuum state remains untouched by the intense THz field.

The dipole response of the Floquet states to the weak NIR probe is calculated perturbatively. The state of the system with both driving fields $`\psi (x,z,t)`$ can be expanded in terms of the Floquet states $`\phi _i(x,z,t)`$ to first order. Discarding second order terms and anti-resonant contributions, the dipole expectation is given by

$$\psi |x|\psi =E_\mathrm{\Omega }e^{i\mathrm{\Omega }t}\underset{i=1}{\overset{2}{}}\frac{\phi _0|x|\phi _i\phi _i|x|\phi _0}{ϵ_i\mathrm{}\mathrm{\Omega }}$$
(11)
The form of (11) is the same as the linear susceptibility of an undriven system, except the states $`\phi _i`$ are oscillating Floquet states instead of stationary states. Explicitly expanding the Floquet states in the numerator gives us an expression for the polarization

$$x(t)=E_\mathrm{\Omega }e^{i\mathrm{\Omega }t}\underset{i=1}{\overset{2}{}}\underset{n,\alpha }{}_{m,\beta }\frac{c_{n,\alpha }^ic_{m,\beta }^ie^{i(nm)\omega t}x_{0\beta }x_{\alpha 0}}{ϵ_i\mathrm{}\mathrm{\Omega }}$$
(12)
The $`+1\omega `$ sideband is given by the Fourier component of $`x(t)`$ oscillating at the frequency $`\mathrm{\Omega }+1\omega `$. Also, for the resonant condition illustrated in Fig. 2, $`\mathrm{}\omega ϵ_1`$, so we keep only the $`i=1`$ term in the sum (12). This condition is satisfied when $`nm=1`$. Thus we obtain an expression for the sideband polarization

$$x_{sideband}(t)=\frac{E_\mathrm{\Omega }e^{i(\mathrm{\Omega }+\omega )t}}{ϵ_1\mathrm{}\mathrm{\Omega }}\underset{n,\alpha ,\beta }{}c_{n+1,\alpha }^1c_{n,\beta }^1x_{0\beta }x_{\alpha 0}$$
(13)
The intensity of the sideband is proportional to $`x_{sideband}^2`$. This is plotted vs. field strength for various detunings in Fig. 5. The detuning parameter $`d`$ is the level spacing normalized by the photon energy of the strong field, $`d=\frac{E2E1}{\mathrm{}\omega }`$.

Our model captures the rollover of the resonant power dependence up to field strengths of around $`\frac{\mu E}{\mathrm{}\omega }1.5`$. Given a calculated excitonic intersubband dipole of $`\frac{\mu }{e}=`$ 12 nm, the strong-field condition is met at a THz field of 5 kV/cm. However, experimentally the sideband never completely disappears, a striking prediction of the theory. It is unlikely that a 3-state simplification is entirely valid because there is a spectrum of other exciton states as well as an electron-hole continuum that may provide significant off-resonant contributions.

To summarize, we drove an quantum-well excitonic intersubband transition with an intense THz laser field in a regime where the Rabi energy was comparable to the photon energy. The nonlinear mixing between the THz pump and a weak NIR probe displayed a nonmonotonic dependence on THz power, which could not be explained by a conventional nonlinear susceptibility. Rather, the excitons were dressed by the THz field, a process which we described in a simple three-level model in which the interaction between excitons and the THz field is solved non-perturbatively within the Floquet formalism.

This research is funded by NSF-DMR 0070083.