text:
# “Dispersion management” for solitons in a Korteweg-de Vries system

## 1 Introduction and formulation of the model

Classical models which give rise to solitons, such as the Korteweg–de Vries (KdV) or nonlinear Schrödinger (NLS) equations, describe uniform nonlinear waveguides. For example, the NLS equation is a basic model to describe light propagation in optical fibers and other guiding structures , while a well-known application of the KdV equation is the description of internal waves in stratified fluids, propagating in waveguides which exist naturally in the ocean, or can be created in a laboratory . In both these models, a (bright) soliton exists under a certain condition: in the NLS equation, the dispersion must be anomalous , while in the KdV equation, the soliton must have a definite polarity.

The necessity to improve stability, bit rate, and other operation characteristics for soliton streams in optical fibers, have recently attracted a great deal of attention to the technique of “dispersion management” (DM), i.e., transmission of solitons in an optical fiber which consists of periodically alternating sections with opposite signs of the dispersion, so that the average dispersion of the long communication link is nearly equal to zero . The corresponding NLS equation with a variable dispersion coefficient is no longer integrable by means of the inverse scattering transform, and no exact soliton solution is known for it. Nevertheless, very accurate numerical simulations show that, to an extremely high accuracy, this equation gives rise to solitons with an approximately Gaussian (rather than the ordinary $`\mathrm{sech}`$), pulsating shape . The existence of dispersion-managed solitons of this type is also strongly supported by analytical results produced by the variational approximation developed for the DM model in different forms , as well as by analysis of the model transformed into an integral equation  (note that the variational approximation can also be applied directly to the integral equation ).

A remarkable property of DM solitons is that they may exist, in a stable form, even when the average dispersion, which is not necessarily exactly zero, but may be much smaller than the local dispersion in the alternating fiber segments, is *normal*  (recall that the ordinary NLS solitons cannot exist in the case of the normal dispersion). It is necessary to mention that, prior to the appearance of the DM models with the above-mentioned piece-wise constant form of the dispersion profile, the NLS equation with a harmonically (sinusoidally) modulated sign-changing dispersion term had already been studied in detail, with a conclusion, based on both the variational approximation and direct simulations, that stable pulsating solitons do exist in that smoothly modulated model .

In fact, the NLS equation with the local dispersion subject to strong periodic modulation belongs to a class of periodic heterogeneous systems, in which stable pulsating solitons with nontrivial properties may be expected (in this work, we use the term “soliton” without assuming integrability of the corresponding model). In nonlinear optics, other systems belonging to this class are tandem waveguides for optical solitons supported by quadratic ($`\chi ^{(2)}`$) nonlinearity, based on alternation of $`\chi ^{(2)}`$ and linear segments , split-step fiber links, in which linear segments alternate with those dominated by the Kerr ($`\chi ^{(3)}`$) nonlinearity , layered bulk media, with the $`\chi ^{(3)}`$ coefficient varying between the layers so that an optical beam propagating across the layers has its power oscillating around a critical value leading to the wave collapse , and alternate nonlinear waveguides, composed of periodically alternating waveguiding and *anti*waveguiding segments . A remarkable feature, common to all these systems, despite their very different physical nature, is robustness of the propagation modes, and the absence of any apparent instability, even when this might naively be expected.

The identification of this class of models, essentially based on the NLS equation, makes it natural to ask whether other ordinary soliton-generating equations, if subjected to a periodic modulation of the dispersion coefficient, can give rise to the propagation of pulsating robust solitary wave (provided that a corresponding model describes a physically meaningful situation). The first candidate to be investigated is the Korteweg-de Vries (KdV) equation. In fact, the KdV equation with variable coefficients is a traditional object for the application of the perturbation theory for solitons . However, the case of periodic modulation, and the consequent possibility of the existence of a quasi-stable pulsating soliton in this case has not yet been studied.

We start with consideration of the propagation of weakly nonlinear long waves in a periodically inhomogeneous waveguide, governed by the variable-coefficient KdV equation:

$$u_t+c(ϵx)u_x+ϵ\left[r(ϵx)uu_x+s(ϵx)u_{xxx}\right]=0(ϵ^2),$$
(1)
where $`ϵ1`$ is a basic small parameter. Equations of this general form are commonly used to describe the propagation of solitary waves in inhomogeneous media, for instance, surface and internal waves in stratified fluids . Provided that the local nonlinear coefficient, $`r`$, and phase velocity, $`c`$, are *nonvanishing* functions of $`ϵx`$, one can introduce a propagation coordinate $`\chi `$ and a temporal variable $`\theta `$ as follows:

$$\chi \frac{ϵr}{6c^2}𝑑x,\theta t\frac{dx}{c}.$$
(2)
Then Eq. (1) reduces, in the same approximation at which it is valid, to a simpler form,

$$u_\chi +6uu_\theta +D(\chi )u_{\theta \theta \theta }=0,$$
(3)
where the local dispersion coefficient is now

$$D(\chi )=6s/\left(rc^2\right).$$
(4)
Two integral quantities, which are frequently called mass and momentum,

$$M=_{\mathrm{}}^+\mathrm{}u(\theta )𝑑\theta ,P=_{\mathrm{}}^+\mathrm{}u^2(\theta )𝑑\theta ,$$
(5)
are exact dynamical invariants of Eq. (3), i.e., $`dM/d\chi =dP/d\chi =0`$. Note that $`P`$ may be called, more accurately, wave-action flux.

Our objective is to consider the case when $`D(\chi )`$ is a periodic function, and, in particular, when it periodically changes its sign. For instance, as it is shown in the Appendix, Eq. (3) describes the propagation of gravity-capillary waves over a periodically-varying-bottom topography near the critical value of the Bond number; a similar derivation can be easily carried out for interfacial-capillary waves. It is also a possible model equation for the propagation of internal waves and Rossby waves against the background of a slowly varying shear flow . In this latter case, the wave modes are localized near critical layers, where the long-wave speed is within the range of the shear velocity. Generally speaking, in all these applications, one should also take into regard the fifth-order dispersion; nevertheless, we initially conjecture that the higher-order dispersion is negligible.

It is assumed that the variable dispersion coefficient $`D(\chi )`$ in Eq. (3) can be represented in the following, quite general, form:

$$D=SD_0(\chi /T)+D_1(\chi ),$$
(6)
where $`D_0`$ is a periodic function with period $`1`$ (i.e., the actual period of the first term in Eq. (6) is $`T`$) and amplitude $`1`$. The mean value of the function $`D_0`$ is *exactly* equal to zero. Following the analogy with the DM models in nonlinear optics, we will then refer to the parameter $`S`$ as the map strength. The function $`D_1`$ then represents the local average dispersion, which may also be subjected to a long-range modulation (a situation with the average dispersion slowly varying along the propagation distance is known for fiber-optic DM as well ). Since we are interested in the case of a sign-changing local dispersion, it is usually assumed that $`|D_1|S/2`$. In fact, without loss of generality, one may set $`S=T=1`$. Indeed, it is straightforward to see that Eq. (3) is mapped into itself, but with $`S=T=1`$, by a transformation to new variables,

$$\stackrel{~}{u}\frac{uT}{(ST)^{1/3}},\stackrel{~}{\chi }\frac{\chi }{T},\stackrel{~}{\theta }\frac{\theta }{(ST)^{1/3}},\stackrel{~}{D}(\tau )D_0(\tau )+S^1D_1(T\tau ).$$
(7)
In section 2, we develop an analytical approximation for solitary-wave solutions to Eq. (3) in the weakly nonlinear case, when the model can be reduced to a relatively simple integral equation. Numerically found solutions to the integral equations are also given in section 2. Results of direct simulations of the full model are displayed in section 3, demonstrating the existence of both stable and unstable solitons. In section 4, we briefly consider effects of the fifth-order dispersion, and the paper is concluded by section 5. The appendix gives a short account of the derivation of Eq. (3), including the fifth-order-dispersion term, for a particular hydrodynamic problem.

## 2 Dispersion-dominated waves

Following the approach of Ref. , where the fiber-optic DM model was transformed into an integral equation, we consider Eq. (3) with the modulated dispersion in the form of Eq. (6) with $`S=T=1`$. We assume that the initial values supplementing Eq. (3) are small, viz., $`f(\theta )u(\theta ,\chi =0)=O(ϵ)`$ (note that, as it follows from Eq. (7), this is equivalent to assuming that $`f=O(1)`$ and $`D`$ varies on a fast scale). It is then natural to introduce a slow time scale $`\tau ϵ\chi `$, redefine the initial values, so that

$$f(\theta )ϵg(\theta ),$$
(8)
and look for a solution in the form

$$u(\theta ,\tau )=ϵu^{(0)}(\theta ,\tau )+ϵ^2u^{(1)}(\theta ,\tau )+\mathrm{}.$$
(9)
Defining the “accumulated dispersion”

$$W(\chi )_0^\chi D_0(\chi )𝑑\chi ,$$
(10)
the zeroth-order solution to Eq. (3) can be found in the form

$$u^{(0)}=^1\left\{\mathrm{exp}\left[ik^3W(\chi )\right]\left\{A(\theta ,\tau )\right\}\right\},$$
(11)
where $``$ and $`^1`$ stand for the Fourier transform with respect to $`\chi `$ and its inverse. From the initial condition (8), it follows that $`A(\theta ,0)=g(\theta )`$ in the solution (11). Because the function $`D_0`$ is periodic with zero mean, the functions $`W(\chi )`$ and $`u^{(0)}`$ given by Eqs. (10) and (11) are also periodic.

At first order in $`ϵ`$, Eq. (3) is reduced to

$$u_\chi ^{(1)}+D_0u_{\theta \theta \theta }^{(1)}=\left[u_\tau ^{(0)}+6u^{(0)}u_\theta ^{(0)}+D_1(\tau )u_{\theta \theta \theta }^{(0)}\right]$$
(12)
with $`u^{(1)}(\theta ,\tau =0)=0`$. Defining

$$h(k,\chi ,\tau )\left\{u_\tau ^{(0)}+6u^{(0)}u_\theta ^{(0)}+D_1(\tau )u_{\theta \theta \theta }^{(0)}\right\},$$
(13)
Eq. (12) with the zero initial value has a solution

$$u^{(1)}=^1\left\{\mathrm{exp}\left(ik^3W\right)_0^\chi \mathrm{exp}\left[ik^3W(\sigma )\right]h(k,\sigma ,\tau )𝑑\sigma \right\}.$$
(14)
To avoid the growth of secularities in the first-order solution (14), we must impose a condition

$$_0^1\mathrm{exp}\left(ik^3W(\chi )\right)h(k,\chi ,\tau )𝑑\chi =0,$$
(15)
which can be written in terms of the Fourier transform, $`\widehat{A}(k)`$, of the function $`A(\theta )`$ introduced in Eq. (11), as

$$_0^1\left\{\widehat{A}_\tau iD_0k^3\widehat{A}+3ik_{\mathrm{}}^{\mathrm{}}\widehat{A}(\kappa )\widehat{A}(k\kappa )\mathrm{exp}\left[3ik\kappa (\kappa k)W(\chi )\right]𝑑\kappa \right\}𝑑\chi =0.$$
(16)
If we introduce a function

$$V(\kappa )_0^1\mathrm{exp}\left[i\kappa W(\chi )\right],$$
(17)
then it follows from Eq. (16) that the amplitude $`A(\theta ,\tau )`$ satisfies an evolution equation:

$$A_\tau +3^1\left\{ik_{\mathrm{}}^{\mathrm{}}V(3k\kappa (\kappa k))\widehat{A}(\kappa )\widehat{A}(k\kappa )𝑑\kappa \right\}+D_0(\tau )A_{\theta \theta \theta }=0,$$
(18)
with the initial condition that $`A(\theta ,0)=g(\theta )`$.

The kernel (17) of the integro-differential equation (18) can be easily calculated in the case of the piecewise-constant dispersion, which is similar to the standard DM scheme in the fiber-optic communications, ,

$$D_0(\chi )=\{\begin{array}{cc}1/2,\hfill & 1/4<\chi <3/4,\hfill \\ +1/2,\hfill & 3/4<\chi <5/4,\hfill \end{array}$$
(19)
which is repeated with the period $`1`$. In this case, the result is

$$V(\kappa )=\left(8/\kappa \right)\mathrm{sin}\left(\kappa /8\right).$$
(20)
For the smooth sinusoidal modulation of the dispersion, with

$$D_0(\chi )=(1/2)\mathrm{cos}\left(2\pi \chi \right),$$
(21)
the kernel can be found in the form

$$V(\kappa )=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{\left(n!\right)^2}\left(\frac{\kappa }{4}\right)^{2n}.$$
(22)
In general, $`V(\kappa )1+O(\kappa ^2)`$ as $`\kappa 0`$, while for $`\kappa 1`$, and assuming that $`D_0(\chi )`$ is sufficiently smooth, the integral can be evaluated using the stationary-phase approximation. The main contribution then comes from a vicinity of the point where $`W^{}(\chi )D_0(\chi )=0`$. Assume that $`D_0(\chi )`$ has two zeroes, $`\chi _1`$ and $`\chi _2`$, which are symmetric, such that $`\chi _2=1\chi _1`$ and $`D_0^{}(\chi _1)=D_0^{}(\chi _2)<0`$ and $`W(\chi _1)=W(\chi _2)`$. Then, the stationary-phase method yields

$$V(\kappa )\sqrt{\frac{8\pi }{\kappa D_0^{}(\chi _2)}}\mathrm{cos}\left(\kappa W(\chi _2)+\frac{\pi }{4}\right).$$
(23)
Hereafter, we focus on steady-state solitary-wave solutions to Eq. (18), following the approach of Ref. . To this end, we assume that the average dispersion $`D_1`$ is a positive constant, and look for solutions in the form $`A=A(\theta S\tau )`$ with a constant velocity $`S`$. Anticipating that $`S0`$, we set $`S\lambda ^2`$, so that the Fourier transform $`\widehat{A}(k)`$ must satisfy an equation following from Eq. (18),

$$(\lambda ^2+D_1k^2)\widehat{A}=3_{\mathrm{}}^{\mathrm{}}V(3k\kappa (\kappa k))\widehat{A}(\kappa )\widehat{A}(k\kappa )𝑑\kappa .$$
(24)
By means of the definition $`\widehat{A}(k)\lambda ^2\widehat{B}(k)`$, one can combine two free parameters $`D_1`$ and $`\lambda ^2`$ into a single one, $`D_1/\lambda ^2`$. Then, given the number $`N`$ of Fourier modes, the following iterative scheme is used to solve the discrete truncated version of Eq. (24):

$`\widehat{B}_k^{(n)}={\displaystyle \frac{3}{\lambda ^2+D_1k^2}}{\displaystyle \underset{l=N/2}{\overset{N/2}{}}}V(3kl(lk))\widehat{A}_l^{(n)}\widehat{A}_{kl}^{(n)},`$ (25)
$`\widehat{A}_k^{(n+1)}={\displaystyle \frac{\widehat{P}\{\widehat{A}^{(n)}\}}{\widehat{P}\{\widehat{B}^{(n)}\}}}\widehat{B}_k^{(n)}.`$ (26)
Here, the superscript is the iteration number, while the integer argument $`k`$ is written as the subscript attached to the amplitudes $`\widehat{B}^{(n)}`$ and $`\widehat{A}^{(n)}`$, and

$$\widehat{P}\{\widehat{A}\}\underset{j=0}{\overset{N1}{}}\left(\widehat{A}_j\right)^2,$$
(27)
is the momentum (see Eq. (5)) of the discrete field $`\widehat{A}_j`$. As for similar solutions of the NLS equation in Ref. , the shape of the steady solution converges; however, the momentum diverges, and hence the wave must be normalized at each step of the iterative method.

In Figs. 1 and 2, localized solutions obtained by means of the numerically implemented iterative scheme (25), (26) are shown, as $`D_1`$ is varied for the two forms (19) and (21) of the dispersion modulation. Note that localized solutions with positive extrema appear to be only possible for $`D_10`$. By symmetry, localized solutions with negative extrema and negative velocity are then only possible for $`D_10`$. As can be seen, a significant difference between the waves for the two dispersion-modulation form is observed only near $`D_1=0`$. For larger $`D_1`$, the mean dispersion dominates, and the waves for both types of the modulation approach the ordinary $`\mathrm{sech}^2`$ soliton of the constant-coefficient KdV equation. For clarity, only the domain $`20x20`$ is shown in these figures; in larger domains, the oscillations apparent in Fig. 2 near $`D_1=0`$ decay to zero.

## 3 Direct numerical results

A crucially important issue is stability of the stationary solitary waves obtained in the previous section by means of the perturbation expansion. This could be addressed by extending the perturbation expansion to the next order, but this approach proves unwieldy. A more straightforward approach, which is undertaken here, is to use solutions of Eq. (24) as initial conditions for the underlying equation (3), and then follow the evolution in direct simulations. With the increase of $`ϵ`$, the nonlinearity gets stronger (see Eq. (9)), leading to stabilization or destabilization of the steady dispersion-dominated solutions. To solve Eq. (3) numerically, we used a standard method combining pseudo-spectral techniques in $`\theta `$ and fourth-order Runge–Kutta integration in $`\chi `$.

In this section, we consider only the case when the average dispersion $`D_1`$ is constant. Further, results presented here are obtained for a slightly smoothed version of the piecewise-constant dispersion modulation (19). Results for other forms of the dispersion modulation, including the sinusoidal form (21), are very similar. Quantitative results, such as the exact position of the stability boundary in Fig. 3, see below, may be affected by the choice of the modulation form.

The transformation (7) demonstrates that, given a particular form of the dispersion modulation, the ratio $`D_1/S`$ and momentum $`P`$ uniquely define a dispersion-dominated solution. Thus, when investigating the stability of dispersion-dominated solutions of Eq. (3), we only have two free parameters, $`D_1/S`$ and $`P`$. An appropriate value of $`ϵ`$ used in the previous section is then $`ϵ=P^{1/2}`$. For the modulation form (19), numerically found stability boundaries of the dispersion-dominated solutions are shown in figure 3.

In classifying the solutions as stable or unstable, we have used a simple criterion that, for unstable solutions, in some region the Poincaré map of the solution’s amplitude (maximum value) will undergo sustained mean algebraic decay in the variable $`\chi `$. Since we solve Eq. (3) with periodic boundary conditions in $`\theta `$, at a late asymptotic stage of the evolution, the decay of the soliton will be followed by establishment of a nearly uniform state, in which the maximum and minimum values of the field are approximately equal and constant in $`\chi `$. Conversely, for stable soliton solutions, we expect that, asymptotically, the absolute mean of the Poincaré map of the maximum value of the field will be much larger than its minimum-value counterpart. Strictly speaking, the location of the stability boundary depends also on the computational domain size $`L`$ and grid size $`\mathrm{\Delta }\theta `$. However, we have taken care to keep $`L`$ sufficiently large and $`\mathrm{\Delta }\theta `$ sufficiently small to minimize their effect. The boundary also depends on the final propagation distance, $`\chi _f`$. All the solutions shown in figure 3 used $`\chi _f=400`$, which was found to be sufficiently large to make it possible to conclude whether a stable soliton was established, or the solution underwent decay.

Figure 3 clearly shows that there are two distinct regions of stable solitary-wave solutions. The stability boundary corresponds to some curve $`P_b=P(D_1/S)`$, with the minimum value of $`P_b`$ occurring at a critical value $`R_c`$ of $`D_1/S`$, where for the piecewise-constant modulation $`R_c0.25`$. Then, the stable region in $`0D_1/S<R_c`$ will be referred to as a “dispersion-dominated” one, while the stable region in $`D_1/S>R_c`$ will be referred to as a ordinary region. The numerical results suggest that the dispersion-dominated region decreases exponentially with the increase of $`D_1/S`$, while the ordinary region increases algebraically with the increase of $`D_1/S`$. In the dispersion-dominated region for $`D_10`$, it is apparent that localized solutions of Eq. (3) are only possible up to a finite value of the momentum. For the piecewise-constant dispersion modulation, this critical value was found to be $`P0.95`$. Thus, for the variable-dispersion KdV equation both nonlinearity and mean dispersion act to destabilize the dispersion-dominated solitons.

The growth of the ordinary stability region with the increase of $`D_1/S`$ is easily understood. In this region, the solitary-wave solutions exist due to balance between nonlinearity and mean dispersion, the effect of variable dispersion is then simply to modulate the form of the wave. In this region, solutions can be investigated using the perturbation theory for the ordinary $`\mathrm{sech}^2`$ solitons , or using methods similar to the guiding-center technique, which were developed for solitons in optical fibers .

In Figs. 4, 5 and 6 examples of solutions corresponding to particular points in Fig. 3 are shown. These three soliton solutions are, respectively, stable, unstable, and stable again. The dispersion-dominated solution in Fig. 4 demonstrates that some momentum is lost from the localized wave to higher-mode oscillatory waves. The subsequent interaction between the oscillatory waves and the solitary wave causes relatively large variability in the Poincaré map of the maximum value of the field; however, both the maximum and minimum values eventually set down to mean values not significantly different from their initial values. The decay of a localized wave is apparent in Fig. 5, where the mean of both the maximum and minimum values of the field can be seen to decrease. Eventually, at large enough propagation distance they would be expected to become approximately equal in the absolute values and opposite in the sign. At this stage, the mass and momentum contained initially in the localized wave would be completely transferred to the oscillatory waves, as is suggested by Fig. 5(b).

For large values of the average dispersion, as is shown in Fig. 6, which is an example of a soliton belonging to the ordinary region, it is apparent that the solitary-wave solutions are now very similar to the $`\mathrm{sech}^2`$ solitons of the constant-dispersion KdV equation. Here, only a small amount of the momentum from the initial wave is lost to the oscillatory waves. Because of the reduced amplitude of the oscillatory wave, the oscillations in the Poincaré map of the maximum field value are now much slower than in Fig. 4.

## 4 Effects of fifth-order dispersion

As was mentioned above, the fifth-order dispersion should, generally speaking, be added to the model of the DM type with a small average dispersion (similarly, third-order dispersion should be included in the optical DM models ). The main effect on the solitons will be to cause them to disperse. In particular, one may expect that, for small values of the average dispersion, the addition of the fifth-order dispersion to Eq. (3) will eventually cause stable solutions to become unstable. It is expected that this effect is most important when the third-order dispersion undergoes a change in sign, i.e. for $`D_1<S/2`$ in Fig. 3. For larger values of $`D_1`$, it is expected that the main effect of the fifth-order dispersion will be to enhance the generation of oscillatory waves . Therefore, here we briefly consider numerically the effect of fifth-order dispersion on a stable dispersion-dominated solution.

With the addition of fifth-order dispersion, Eq. (3) becomes

$$u_\chi +6uu_\theta +D(\chi )u_{\theta \theta \theta }+\delta ^2u_{\theta \theta \theta \theta \theta }=0,$$
(28)
where, as shown for the case of gravity-capillary waves in the Appendix, one can assume that the coefficient $`\delta ^2`$ is constant, in contrast to the strongly modulated third-order-dispersion coefficient. In Fig. 7, the effect of varying this parameter on the stable solution from Fig. 4 is shown. These numerical solutions are obtained using a simple extension of the pseudo-spectral Runge-Kutta method described in the previous section. First, the case of $`\delta ^2=10^4`$ was investigated, but it is not shown in Fig. 7, as it turns out to be practically identical to $`\delta ^2=0`$. Thus, for sufficiently small $`\delta ^2`$, stable solitary waves are established as well as in the case $`\delta ^2=0`$. However, as $`\delta ^2`$ is increased, the solitons disperse indeed, leaving oscillatory-wave radiation. Other simulations in the dispersion-dominated region confirm that the critical (maximum) value of $`\delta ^2`$, at which the solitary waves are able to form, decreases as both $`D_1/S`$ and $`P`$ increase, so that the critical value is $`\delta ^2=0`$ on the stability boundary in Fig. 3.

## 5 Conclusions

In this work we have introduced a KdV model with a periodically varying dispersion coefficient that takes both positive and negative values. Using a perturbative expansion based on integral equations and direct numerical simulations, we have shown that for the case when the dispersion undergoes a periodic change in sign, stable solitary-wave solutions are possible in a region defined by the average dispersion and the initial momentum of the system. Further, it has been shown that sufficiently weak fifth-order dispersion, which, generally speaking, should be taken into regard when the average third-order dispersion is very small, does not destroy these solitary waves.

It should be noted that only one analytical approach (based on the integral equation) to investigate the possibility of solitary-wave solutions has been presented in detail in this work. We have also tried two additional approaches. The first of these is the balance-equation method. Equation (3) can be derived from the Hamiltonian

$$H=_{\mathrm{}}^+\mathrm{}\left[\frac{1}{2}D(\chi )u_\theta ^2+u^3\right]𝑑\theta ,$$
(29)
which, unlike the mass and momentum given by Eq. (5), is not conserved when the dispersion coefficient is variable. Instead, it evolves with $`\chi `$ according to an immediate corollary of (3),

$$\frac{dH}{d\chi }+\frac{1}{2}\frac{dD}{d\chi }_{\mathrm{}}^+\mathrm{}u_\theta ^2𝑑\theta =0.$$
(30)
The balance-equation approach then demands a choice of a slowly varying ansatz for the form of the solitary wave. Using the conservation equations for the mass and momentum, and the evolution equation (30) for the Hamiltonian, evolution equations for parameters of the ansatz can be derived. The second approach is a numerical averaging procedure, involving successive Poincaré maps from direct numerical simulations of Eq. (3) in an attempt to generate a stable soliton solution by iterations. However, neither of these methods have so far proved successful for the variable-dispersion KdV equation.

The reason for solitary-wave solutions being less forthcoming from the variable-dispersion KdV equation than for the variable-dispersion NLS equation may have to do with the zero-dispersion limit of these two equations. For the KdV equation, it is a singular limit with wavebreaking occurring in finite time, whereas for the NLS equation in the zero-dispersion limit, only the expected phase modulation takes place. Thus, when nonlinearity dominates over the variable dispersion in the KdV equation, as in the limit $`D_1/S=0`$ and $`P1`$ in Fig. 3 above, one would not expect that stable solutions are possible, as the evolution generated by the nonlinearity is much faster than that under the action of the variable dispersion. Nonetheless, combining the results presented here with the well-known ones for the dispersion-managed solitons in the NLS models of long fiber-optic links consisting of alternating segments with anomalous and normal dispersion, we conjecture that stable pulsating solitary waves may be possible in a vast class of nonlinear wave equations subject to strong periodic modulation of the dispersion coefficient.

## Acknowledgment

B.A.M. appreciates the hospitality of Loughborough University (UK) and Monash University (Clayton, Australia).

## Appendix A Zero dispersion limits of long, weakly-nonlinear waves: gravity-capillary waves

A fundamental requirement for the application of Eq. (1) is that the linear dispersion relationship between the velocity $`c`$ and the wavenumber $`k`$ in the long wave limit can be written as

$$c=c_0\alpha k^2+O(k^4),$$
(31)
where $`c_0`$ is non-zero, and where for some special parameter combination, the coefficient $`\alpha =0`$. In geophysical hydrodynamics, the modulation of the coefficients of (1) can then be induced by a slowly-varying environment, e.g., flow past periodically varying topography. We briefly describe here an example where this may occur, namely for the flow of a long, weakly-nonlinear, gravity-capillary waves over periodically-varying topography when the Bond number is close to $`1/3`$. A similar derivation holds for interfacial-capillary waves.

Consider then the propagation of gravity-capillary waves in water of depth $`h`$ with a horizontal length scale $`L=h/\delta `$ and amplitude $`a=ϵh`$, such that $`\delta ,ϵ1`$. Let $`g`$ be the gravity acceleration, $`\rho `$ the density of the water, and $`\sigma `$ the coefficient of surface tension. The linear dispersion relationship is

$$c=c_0(1\alpha (hk)^2+\beta (hk)^4+O(k^6)),$$
(32)
where

$$c_0=(gh)^{1/2},\alpha =\frac{1}{2}\left(\frac{1}{3}B\right),\beta =\frac{1}{2}\left(\frac{2}{15}\frac{B}{3}\alpha ^2\right)$$
(33)
and $`B=\sigma /\rho gh^2`$ is the Bond number. Therefore, near $`B=1/3`$ the leading-order dispersion coefficient is approximately zero. This occurs for water of the depth $`h1\mathrm{c}\mathrm{m}`$, in which case one should generally also include diffusive effects. However, we will assume here that these are negligible; also we note that for the analogous situation of interfacial-capillary waves, the effect of gravity is reduced and correspondingly the water depth may not be quite so small, so that in that case it is more plausible to ignore diffusive effects.

For one-dimensional waves propagating to the right, the evolution of the free-surface displacement $`ϵ\eta `$ is then governed by the extended KdV equation

$$\frac{1}{c_0}\frac{\eta }{t}+\left(1+ϵ\frac{3\eta }{2h}\right)\frac{\eta }{x}+\delta ^2\alpha h^2\frac{^3\eta }{x^3}+\delta ^4\beta h^4\frac{^5\eta }{x^5}=O(ϵ^2,ϵ\delta ^2,\delta ^6),$$
(34)
Now assume that, in addition, the depth is a slowly varying function of $`x`$, such that $`h=h(ϵx)`$, and introduce new variables

$$\chi =ϵx,\theta =g^{1/2}\left(t\frac{dx}{c}\right).$$
(35)
Further, we put $`A\eta h^{\frac{1}{4}}`$ and note that in a slowly varying environment the wave-action flux density $`A^2`$ is conserved . Hence, the corrected form of Eq. (34) is

$$ϵ\frac{A}{\chi }ϵ\frac{3A}{2h^{7/4}}\frac{A}{\theta }\delta ^2\alpha h^{\frac{1}{2}}\frac{^3A}{\theta ^3}\delta ^4\beta h^{\frac{3}{2}}\frac{^5A}{x^5}=O(ϵ^2,ϵ\delta ^2,\delta ^6).$$
(36)
Typically $`\alpha `$ is $`O(1)`$; then, let $`ϵ=\delta ^2`$, which leads to the variable-coefficient KdV equation governing the wave propagation,

$$\frac{A}{\chi }\frac{3A}{2h^{7/4}}\frac{A}{\theta }\alpha h^{\frac{1}{2}}\frac{^3A}{\theta ^3}=0.$$
(37)
However, $`B=1/3`$ is a singular limit of this equation; in this limit we set $`\alpha =\delta ^2\gamma `$ and $`ϵ=\delta ^4`$, hence Eq. (36) becomes

$$\frac{A}{\chi }\frac{3A}{2h^{7/4}}\frac{A}{\theta }\gamma h^{1/2}\frac{^3A}{\theta ^3}\beta h^{3/2}\frac{^5A}{\theta ^5}=0,$$
(38)
where now $`\beta =1/90`$. Thus, if the topography and Bond number are such that the coefficient $`\alpha `$ is close to zero, and can change its sign, then a uniformly valid model equation including the fifth-order dispersion is

$$\frac{A}{\chi }\frac{3A}{2h^{7/4}}\frac{A}{\theta }\alpha h^{1/2}\frac{^3A}{\theta ^3}\frac{\delta ^2h^{3/2}}{90}\frac{^5A}{\theta ^5}=0.$$
(39)


format: text

text:
# Gamma Ray Bursts as Probes of the First Stars

## 1 Introduction

The high redshift frontier of observational cosmology currently stands at redshifts $`z6`$. The current redshift record is a quasar at $`z=5.8`$, and a few galaxies are known at marginally lower redshift. Beyond $`z=6`$, we have yet to identify any individual objects. We do know that hydrogen was predominantly neutral at redshifts $`z>30`$ based on the observed anisotropy of the cosmic microwave background, which would be smoothed out by Thomson scattering if the free electron density at $`z>30`$ were too great. The redshift range $`6<z<30`$ remains unknown territory. It is a very interesting territory, too, for it should include the formation of the first stars, galaxies, and quasars, and certainly includes the epoch at which hydrogen was reionized.

Searches for starlight (and other rest-frame near ultraviolet tracers) can make incremental progress into the low-redshift end of this period. However, these methods face a practical limit where the Lyman break redshifts out of the optical window to the near-infrared, at $`z7`$. At higher redshift, essentially no flux is expected in the optical window (observed wavelengths $`0.36\mu m<\lambda _{obs}<1\mu m`$). Atmospheric conditions and present detector technologies conspire to make searches at $`\lambda _{obs}>1\mu m`$ much less efficient. Future instrumentation like the Next Generation Space Telescope (NGST) promise extensions of “conventional” optical methods to the observed near-IR and thus to redshifts $`z6`$, but this may be a decade or more away. In the meantime, we expect the upcoming extension of our Large Area Lyman Alpha (LALA) survey (Rhoads et al 2000; Malhotra et al 2001) to $`z=6.6`$ to be at or near the practical limit for some years.

We would like to find tracers of $`z>6`$ objects that are accessible now. Fortunately, this is possible so long as we are willing to use something besides starlight. In practice, this means higher energy photons ($`\gamma `$ and x-rays), since lower energies still face either confusion or sensitivity issues.

Gamma ray bursts (GRBs) are an excellent candidate for detection at high redshift because the bursts and their afterglows are extremely bright at all wavelengths. Two conditions must be met for such a candidate to work well. First, there should be a reasonable expectation that the object exists at high redshift; and second, it should be detectable there.

The best argument that gamma ray bursts should occur at high redshift comes from the growing body of evidence linking GRBs to star formation activity (and hence presumably to the deaths of massive, short-lived stars): GRB host galaxy colors are characteristically blue (Fruchter et al 1999); the spatial distribution of GRBs on their hosts matches expectations for hypernova models (Bloom, Kulkarni, & Djorgovski 2000); and the emission lines of GRB host galaxies are unusually strong (Fruchter et al 2001). Structure formation models yield estimated redshifts $`z15\pm 5`$ for the first stars to form in the universe (cf. Barkana & Loeb 2001). This is supported by studies of heavy element abundances: It has proven extremely difficult to find objects with primordial (i.e., big bang nucleosynthesis) abundances at any redshift currently accessible. The immediate inference is that a substantial generation of stars must have existed at earlier redshifts to produced the ubiquitous metals. The association of GRBs with star formation then implies that the first GRBs also occurred in the redshift range $`z15\pm 5`$.

The detectibility of GRBs at $`z6`$ has been considered in detail by Lamb & Reichart (2000), who find that the bright end of the luminosity distribution would be detectable at very high redshifts (though quantitative predictions depend substantially on unknown details of the GRB luminosity function). This applies also to the X-ray and optical afterglows, for which time dilation of the most distant afterglows helps offset the increase in luminosity distance with redshift (Lamb & Reichart 2000; Ciardi & Loeb 2000).

The Lyman break will render afterglows at $`z>7`$ invisible to optical detectors, just as it does for galaxies. But the problem here is not so serious. Searches for $`z>7`$ galaxies suffer because galaxies at such high redshifts are faint, and the required combination of large solid angles and high sensitivity to find them is not yet practical at near-IR wavelengths. Because GRB afterglows outshine their host galaxies at early times, and because X-ray detectors can determine GRB locations with accuracy comparable to the current near-IR field of a 4m class telescope, an afterglow at this redshift is easier to find than are the galaxies around it. Indeed, published near-infrared afterglow observations (Rhoads & Fruchter 2001) already achieve a sensitivity sufficient to detect afterglows at $`z10`$ for several hours following a GRB (cf. figures 2,3 of Lamb & Reichart 2000). The followup program described in Rhoads & Fruchter 2001 is continuing at the NASA Infrared Telescope Facility, and we have a similar program at the National Optical Astronomy Observatory. The observed signature of a $`z>7`$ GRB would be a near-infrared afterglow exhibiting a Lyman break at $`\lambda _{obs}=0.1215(1+z)\mu m>1\mu m`$. Such breaks have been used to measure $`z=2.05`$ for GRB 000301C (Smette et al 2001) and to estimate $`z5`$ for GRB 980329 (Fruchter 1999; see also Reichart et al 1999). Their extension to longer wavelengths is straightforward. Thus, it is reasonable to expect that $`z>7`$ GRBs will be detected with current technology.

The prospect of detecting gamma ray bursts at $`z>7`$ opens two possible methods of studying star formation activity at these epochs: GRB rate evolution, which should trace star formation activity; and pair production cutoffs in the GeV spectra of bursts, which probe the total optical-ultraviolet background light produced by high redshift stars.

## 2 Burst Rate Evolution

The most basic inference from the observed burst rate is that the highest redshift where a burst has been detected $`z_{max,grb}`$ implies the onset of star formation at some redshift $`z_{max,}>z_{max,grb}`$. It is likely that in fact $`z_{max,}z_{max,grb}`$: The association of GRBs with star formation tracers requires short progenitor lifetimes ($`10^8`$ years), so the redshift difference between the first stars formed and the earliest possible hypernovae is small.

It will be possible to go further by measuring the GRB rate as a function of redshift, $`R_{grb}(z)`$, and taking it as a surrogate for the star formation rate. Such studies would require a large sample (several tens) of high redshift GRBs, together with an understanding of the selection effects that went into the sample. This method is likely to be limited by at least two systematic factors. First, uncertainties in the GRB luminosity function will introduce uncertain corrections to the inferred total GRB rate and the inferred star formation rate, since the high redshift sample will contain only bright bursts. Second, evolution in the burst progenitor population may influence the burst rate. One plausible example is that the GRB rate could depend on progenitor metallicity, which is likely to be lower in the early universe. Another is that the stellar initial mass function (IMF) may vary, thereby affecting the relation between GRB rate and star formation rate, and perhaps also the shape of the GRB luminosity function. Possible evidence for IMF variations has recently been found in at some high redshift Lyman $`\alpha `$ emitting galaxies (Malhotra et al 2001).

Overall, these complications suggest that calibration of the GRB rate as an indicator of global star formation might be possible to within a factor of a few. While higher accuracies would be desirable, the present uncertainties with more conventional star formation estimators are not much better. For example, rest ultraviolet continuum measurements are corrected by a factor of $`7`$ for dust absorption, and the uncertainty in this correction could easily be a factor of two given the range of possible dust properties.

## 3 Pair Production Cutoffs in GRB Spectra

The observed spectra of gamma ray bursts sometimes extend to very high photon energies: The EGRET experiment on the Compton Gamma Ray Observatory detected four bursts with unbroken power law tails extending to $`_\gamma >1\text{GeV}`$, and the Milagrito air shower experiment has tentatively detected one burst at $`_\gamma >1\text{TeV}`$. Photons with such high energies have mean free paths shorter than a Hubble distance due to $`\gamma +\gamma e^++e^{}`$ interactions with low energy background photons. The threshold for such pair production reactions is $`_\gamma ϵ_\gamma >m_e^2c^4=(511\text{keV})^2`$, corresponding to the requirement that each photon have the rest mass energy of an electron in their center of momentum frame. (Here $`_\gamma `$ and $`ϵ_\gamma `$ are the two photon energies measured in an arbitrary frame, and $`_\gamma ϵ_\gamma `$ by convention.) The cross section (for a head-on collision) peaks at $`_\gamma ϵ_\gamma =2m_e^2c^4`$ and falls asymptotically as $`1/(_\gamma ϵ_\gamma )`$ for $`_\gamma ϵ_\gamma 2m_e^2c^4`$.

Pair production cutoffs in the TeV gamma ray spectra of blazars due to interactions with the cosmic infrared background have been predicted (Stecker, De Jager, & Salamon 1992; MacMinn & Primack 1996; Madau & Phinney 1996; Malkan & Stecker 1998) and observed (e.g., De Jager, Stecker, & Salamon 1994; Konopelko et al 1999) for several years now. The extension of the same physics to higher redshifts and lower gamma ray energies has been explored recently by several groups (Salamon & Stecker 1998; Primack et al 2000; Oh 2000).

The observer frame gamma ray energy determines simultaneously the redshift and rest frame energies of the background photons that dominate the pair production optical depth. At low redshifts ($`z1`$), the effective absorption coefficient $`\alpha (_\gamma )`$ increases with $`_\gamma `$ and changes relatively little with redshift, so that the relevant physics is simply $`\alpha (_{cut})d=1`$, with $`d`$ the distance to the source. However, at $`z>1`$, redshift effects become important: The threshold energy $`ϵ_\gamma (z)1/(1+z)`$, and the background radiation field will also evolve with redshift. The optical depth for photons near $`_{cut}`$ is therefore dominated by absorption at high redshift, unless the source redshift is so high as to precede the creation of any substantial optical-IR background. By the time the photon reaches lower redshifts, the threshold for pair creation grows so large that the density of relevant photons is extremely low. Oh (2000) has shown that the highest energy background photons capable of producing optical depth $`\tau 1`$ over a Hubble distance have energies below the ionization threshold for hydrogen (i.e., $`ϵ_\gamma <13.6\text{eV}`$), since hydrogen absorption in stellar atmospheres, galaxies, and the intergalactic medium ensures a strong decrement in background photon number density at $`13.6\text{eV}`$.

The most robust observable consequence of the pair creation cutoff is the observer frame gamma ray energy $`_{cut}(z)`$ for which the optical depth $`\tau =1`$. Lower pair creation optical depths ($`\tau 1`$) cannot be measured reliably because of our imperfect knowledge of the intrinsic (i.e., unabsorbed) source spectrum, while at higher optical depths ($`\tau 1`$) the absorption reduces the flux below detection thresholds of present or near-future instruments. We might measure $`\tau (_\gamma )`$ with reasonable accuracy over the range $`1/2<\tau <2`$.

Detailed predictions of $`_{cut}(z)`$ differ from model to model, depending on the theoretical treatment adopted for the earliest star formation (see Primack et al 2000; Oh 2000). For example, the observer frame energy where $`\tau =1`$ for redshift $`z=6`$ is $`4\text{GeV}<_{cut}(6)<6\text{GeV}`$ for different models in Primack et al (2000), and $`10\text{GeV}<_{cut}(6)<26\text{GeV}`$ for models in Oh (2000). Therein lies the power of this method for learning about the first generations of stars, for these strong differences in predictions allow the models to be distinguished with comparative ease from even a modest data set.

Moreover, if we can observe the GeV cutoffs in spectra of a few GRBs spread over the redshift range $`6<z<z_{max,}`$, we can infer the evolution of the optical-UV background radiation over the same period with little dependence on models. This follows because the difference in pair creation optical depth between two bursts at redshifts $`z_1`$ and $`z_2`$ ($`z_1<z_2`$) is determined only by the background radiation in the range $`z_1<z<z_2`$.

## 4 Discussion

The two methods of using gamma ray bursts to probe high redshift star formation complement each other in many ways. GRB rate measurements at high redshift are technically easier. They require a GRB monitor plus rapid multiband near-infrared followup. Existing instrumentation and indeed existing observational programs are already adequate for this work. Pair creation cutoffs require one additional observation, namely, a GeV energy spectrum obtained during the GRB. This GeV spectrum will have to come from GLAST or a similar space mission.

The physical assumption behind the GRB rate evolution method is that the bursts are associated with star formation activity. Under this assumption, there will be some systematic uncertainties in converting the GRB rate to the star formation rate (see above). In contrast, the pair creation cutoff method requires only that some high redshift GRBs have GeV spectra that are sufficiently bright and sufficiently smooth for the cutoff to be observed. Beyond this, there is no requirement on the nature of the bursters, which are needed only as beacons to probe the intervening background radiation. The physics of pair creation is then well understood and probes the total background radiation produced by high redshift stars.

Thus, combining the two methods of studying high redshift star formation with GRBs may overcome the physical uncertainties of either method alone. Additional constraints from other techniques using other classes of objects (galaxies observed at infrared wavelengths, or quasars at X-ray wavelengths) will become available over the next few years, and will again have complementary strengths and weaknesses. By adding these to the GRB results, we can reasonably expect to understand star formation at $`z10`$ as well as we understand it at $`z3`$ today.


text:
# Sequential Bifurcations in Sheared Annular Electroconvection

## Abstract

A sequence of bifurcations is studied in a one-dimensional pattern forming system subject to the variation of two experimental control parameters: a dimensionless electrical forcing number $``$ and a shear Reynolds number $`\mathrm{Re}`$. The pattern is an azimuthally periodic array of traveling vortices with integer mode number $`m`$. Varying $``$ and $`\mathrm{Re}`$ permits the passage through several codimension-two points. We find that the coefficients of the nonlinear terms in a generic Landau equation for the primary bifurcation are discontinuous at the codimension-two points. Further, we map the stability boundaries in the space of the two parameters by studying the subcritical secondary bifurcations in which $`mm+1`$ when $``$ is increased at constant $`\mathrm{Re}`$.

In typical nonlinear pattern-forming systems, there is one easily adjusted quantity, the control parameter, whose variation results in a symmetry-breaking bifurcation to an ordered pattern. If a second control parameter is available, then it may be possible to tune the parameters to codimension-two points at which two modes with different symmetries are equally unstable. Near such points, neighboring unstable modes compete with one another leading to low-lying secondary bifurcations even in the weakly nonlinear regime which is accessible to perturbative theoretical treatment. There are few experimental systems whose patterns can be conveniently manipulated by the independent variation of two control parameters. Examples include Rayleigh-Bénard convection (RBC) in a rotating cylinder, convection in binary fluid mixtures, and Taylor Vortex Flow (TVF) with rotation of both the inner and outer cylinders. In most cases, the eventual two-dimensionality of the patterns results in a very complex phenomenology.

In this Letter, we study experimentally the bifurcations in a thin annular fluid film. The two-dimensional (2D) film is both sheared and driven into electroconvection. This system is unique in that it involves two easily adjusted control parameters, the applied voltage and the shear, with the bifurcations to uni-directionally traveling vortices only in one-dimension (1D). The annular geometry makes the 1D pattern effectively periodic and free of boundary effects. Most other pattern-forming systems are fluid-mechanically three-dimensional (3D) and the competing patterns often have different symmetries and dimensions. For example, in TVF there is a codimension-two point where Taylor and wavy vortices are simultaneously unstable. The former (latter) pattern is 1D (2D). The greatly increased size of the parameter space made available by the third dimension makes it more difficult to examine experimentally the instabilities to and between the competing pattern states. Further, it is important to understand which physical processes depend on the different symmetries of the competing modes and which on their different dimensions. In our system, symmetry alone determines the nature of the bifurcations.

The highly restricted geometry of our fluid film only allows the competition of azimuthally periodic 1D modes of various $`m`$-fold symmetries. As a function of voltage and shear, we find a simple sequence of primary and secondary bifurcations. The primary bifurcation is the transition from the conducting state to the convecting state, whereas the secondary bifurcations connect one convecting state to another. This scenario can be modelled by coupled Landau amplitude equations which can be calculated completely by perturbative methods from the primitive electrohydrodynamic equations. The amplitude equations reduce to a single generic Landau equation when applied to the primary bifurcation. We experimentally demonstrate that the coefficients of the nonlinear terms in this equation vary discontinuously at codimension-two points. Our experiments also reveal that the secondary bifurcations near codimension-two points are subcritical and suppressed by the shear.

The experimental system consists of an annular, freely suspended film of smectic A liquid crystal. The layered smectic A phase is isotropic and newtonian for flows in the plane of the film, and maintains a uniform thickness even while flowing. Such a film is driven to electroconvect by a voltage applied between the inner and outer edges which sets up an unstable distribution of surface charge on the film’s free surfaces. A Couette shear is applied by mechanical rotation of the inner electrode. Previous studies of this instability explored the behavior of the primary bifurcation to electroconvection, confirming that the weakly nonlinear regime of convection near onset is described by a complex amplitude equation of the Landau form for a single mode number $`m`$. The concentric electrodes have inner (outer) radii $`r_i`$ ($`r_o`$), from which we define the radius ratio $`\alpha =r_i/r_o`$. The film forms an annular sheet of width $`d=r_or_i1`$mm and thickness $`s0.1\mu `$m so that $`s/d10^4`$. The fluid is treated as incompressible, strictly two-dimensional, and of uniform electrical conductivity. We denote the fluid density, molecular viscosity, and conductivity by $`\rho `$, $`\eta `$, and $`\sigma `$. In a typical experiment, a steady circular Couette shear was imposed by rotating the inner electrode at $`\omega `$ rad/s and a dc voltage $`V`$ was applied between the inner and outer electrodes. Current-voltage characteristics were obtained for incremental and decremental voltages. The state of the film is specified by the radius ratio and three additional dimensionless parameters. The control parameter $`=ϵ_0^2V^2/\sigma \eta s^2`$, where $`ϵ_0`$ is the permittivity of free space, is a measure of the external driving force. The Prandtl-like parameter $`𝒫=ϵ_0\eta /\rho \sigma sd`$ is the ratio of the time scales of electrical and viscous dissipation processes in the film. The Reynolds number $`\mathrm{Re}=\rho \omega r_id/\eta `$ indicates the strength of the applied shear and is the second control parameter. Details of the experimental apparatus and procedure can be found in Refs. .

The primary bifurcation leads to a pattern of $`m`$ vortex pairs. Onset occurs at a critical voltage $`V=V_c`$ with the critical mode number $`m=m_c`$ becoming unstable. Linear stability predicts the values of $`m_c`$ and the critical control parameter $`_cV_c^2`$ at which the film becomes marginally unstable. In general, $`_c`$ and $`m_c`$ are functions of $`\alpha `$, $`𝒫`$, and $`\mathrm{Re}`$. However, when $`\mathrm{Re}=0`$, $`_c=_c^0`$ and $`m_c=m_c^0`$ are independent of $`𝒫`$. When $`\mathrm{Re}>0`$, it is found that $`_c>_c^0`$, for any $`\alpha `$ and $`𝒫`$, i.e. that the shear always acts to stabilize the conduction state and suppresses the onset of convection. In addition for nonzero shear, $`m_cm_c^0`$. Both these predictions of linear theory have been experimentally confirmed. Our main experimental tool is the total electric current $`I`$, which is a measure of the rate of charge transport due to both conduction and convection. The reduced Nusselt number $`n`$ is the ratio of the current due to electroconvection to that due to conduction and is defined as $`n=(I/I_{cond})1=(I/cV)1`$. The conductance of the film $`c`$ is obtained from the current-voltage characteristic of the film when it is not convecting. Similarly, the reduced control parameter $`ϵ`$, is given by $`ϵ=(/_c)1=(V/V_c)^21`$. The raw $`(V,I)`$ data can thus be transformed to measurements of $`(ϵ,n)`$. In these dimensionless terms, the primary bifurcation occurs at $`ϵ=0`$. When $`ϵ`$ is increased, the secondary bifurcations are transitions between flows with mode numbers $`mm+1`$. Upon decreasing $`ϵ`$, the transitions $`mm1`$ are observed. The mode changes can be observed directly, but rather crudely, using visual observations of slightly nonuniform films, or films seeded with small particles. They appear very cleanly, if more indirectly, in the form of jumps in $`n`$ as a function of $`ϵ`$. These observations are in qualitative agreement with an equivariant bifurcation theory. Increasing (decreasing) $`ϵ`$ results in secondary bifurcations at $`ϵ_1<ϵ_2<\mathrm{}`$ ($`ϵ_1<ϵ_2<\mathrm{}`$) corresponding to mode transitions $`mm+1\mathrm{}`$ ($`\mathrm{}m+1m`$).

In Figure 1 are illustrated two representative examples of secondary bifurcations. Figure 1(a) shows $`(ϵ,n)`$ data for a sheared film with subcritical primary and secondary bifurcations. The secondary bifurcation, which results in one additional traveling vortex pair, appears at $`ϵ_10.22`$ for increasing $`ϵ`$. When $`ϵ`$ is decreased, the removal of the vortex pair occurs at $`ϵ_10.12`$. Convection altogether ceases at $`ϵ0.1`$. Figure 1(b) shows a case where, for increasing $`ϵ`$, the primary and secondary bifurcations appear as a single strongly subcritical bifurcation at $`ϵ0`$. However, for decreasing $`ϵ`$, $`ϵ_1`$ is distinct from $`ϵ0.11`$ where convection stops. That shear selects the mode which is linearly unstable is well established in the case of the primary bifurcation. It has been shown from linear stability analysis that the marginally unstable mode decreases from $`m_c^0m_c^01m_c^02\mathrm{}`$ as $`\mathrm{Re}`$ is increased. Hence, for fixed $`\alpha `$ and $`𝒫`$, there are codimension-two points at particular Reynolds numbers, $`\mathrm{Re}_{m,m1}`$, at which both $`m`$ and $`m1`$ are marginally unstable. For $`\mathrm{Re}<\mathrm{Re}_{m,m1}`$, the primary bifurcation is to a convecting state with mode number $`m`$. The first secondary bifurcation occurs at $`ϵ=ϵ_1>0`$ to mode number $`m+1`$. For $`\mathrm{Re}>\mathrm{Re}_{m,m1}`$, the primary bifurcation is to mode $`m1`$. The first secondary bifurcation is at $`ϵ=ϵ_1>0`$ to mode $`m`$. Clearly, it follows that as $`\mathrm{Re}\mathrm{Re}_{m,m1}`$ from above, $`ϵ_10`$, resulting in the situation depicted in Fig. 1(b). Our experimental protocol was always to hold $`\mathrm{Re}`$ constant and vary $`ϵ`$. Neither experiment nor theory has addressed the question of what mode transitions occur if $`\mathrm{Re}`$ is varied for constant $`ϵ`$.

Figures 2(a)-(h) show a sequence of $`(ϵ,n)`$ data plots for $`\alpha =0.56`$, $`𝒫=76`$, and a range of closely spaced values of $`\mathrm{Re}`$. In Fig. 2(a), the data at $`\mathrm{Re}=0.124`$ shows a subcritical primary bifurcation with mode number $`m=8`$. For the range of $`ϵ`$ investigated, secondary bifurcations were not encountered. In Fig. 2(b), the data at $`\mathrm{Re}=0.142`$ shows a subcritical primary bifurcation with $`m=7`$ followed by a subcritical secondary bifurcation ($`m=78`$) at $`ϵ_10.26`$. Hence at some Reynolds number, $`0.124<\mathrm{Re}_{8,7}<0.142`$, there is nascent at $`ϵ=0`$ the mode change $`87`$. At $`\mathrm{Re}=0.160`$, the secondary bifurcation is suppressed to higher $`ϵ_10.48`$, as shown in Fig. 2(c). The sequence of events in Figs. 2(a) through (c) repeats with increasing $`\mathrm{Re}`$. In Fig. 2(f), data at $`\mathrm{Re}=0.214`$ shows a primary bifurcation to $`m=6`$ and a secondary bifurcation at $`ϵ_10.07`$ to $`m=7`$. Hence a second codimension-two point at $`0.196<\mathrm{Re}_{7,6}<0.214`$ is traversed. At $`\mathrm{Re}=0.231`$ and $`\mathrm{Re}=0.249`$, the secondary bifurcation is suppressed to $`ϵ_10.11`$ and $`ϵ_10.16`$ respectively, as is shown in Figs. 2(g) and (h). The data obtained while $`ϵ`$ is being reduced has equally interesting behavior. For example, for increasing $`ϵ`$ in Fig. 2(f) there is a $`m=67`$ bifurcation, but for decreasing $`ϵ`$ the $`m=7`$ convecting state exits directly to conduction by a single discontinuous transition, without recourse to the $`m=6`$ state. However, for the slightly higher $`\mathrm{Re}`$ in Figs. 2(g) and (h), the $`m=6`$ state is revisited while $`ϵ`$ decreases.

For nonzero shear, symmetry and perturbative analyses of the primitive equations lead to the amplitude equation

$`\tau (\dot{A}_mia_{Im}A_m)=ϵ(1+ic_{0m})A_m`$ (1)
$`g_{0m}(1+ic_{2m})|A_m|^2A_mg_{1m}(1+ic_{3m})|A_n|^2A_m`$ (2)
$`h_{0m}(1+ic_{4m})|A_m|^4A_m+\mathrm{},`$ (3)
and its counterpart under $`mn`$ symmetry. In the above, all of the coefficients are real. The steady state solutions of the real part of Eq. 3 are consistent with the experimental observations of a state $`m`$ being replaced by a state $`n=m\pm 1`$ via a discontinuous transition. To further interpret the data in Fig. 2 and better understand the behavior of the primary bifurcation as the Reynolds number is increased through $`\mathrm{Re}_{8,7}`$ and $`\mathrm{Re}_{7,6}`$, we consider the real and steady state part

$$ϵAgA^3hA^5=0,$$
(4)
of Eq. 3, where $`A`$ is real, $`gg_{0m}`$, and $`hh_{0m}`$. We have truncated at quintic order in the expansion. Since the reduced Nusselt number $`n=A^2`$, we can determine $`g`$ and $`h`$ for each primary bifurcation in the sequence in Fig. 2. Our results are plotted in Fig. 3. We find that both $`g`$ and $`h`$ are strongly discontinuous as $`\mathrm{Re}`$ is increased through the codimension-two points $`\mathrm{Re}_{8,7}`$ and $`\mathrm{Re}_{7,6}`$. The discontinuity in both $`g`$ and $`h`$ is smaller at $`\mathrm{Re}_{8,7}`$ than at $`\mathrm{Re}_{7,6}`$. This raises the possibility that these discontinuities may vanish in the limit $`m\mathrm{}`$ and so are a feature of the finite size and discrete mode structure of this system.

A minimal description of a state consists of identifying the value of its mode number $`m`$. Allowing for hysteresis, each state persists over certain overlapping ranges of $`\alpha `$, $`𝒫`$, $``$, and $`\mathrm{Re}`$. Here, we focus on fixed $`\alpha `$ and $`𝒫`$ for simplicity. In Fig. 4 we map the solutions in the subspace defined by $`(ϵ,\mathrm{Re},\alpha =0.56,𝒫=76)`$, where the reduced variable $`ϵ`$ is used for convenience. Figure 4 shows the values of $`ϵ`$ at which conduction and convection with various $`m`$ were found, for given values of $`\mathrm{Re}`$. Below (above) the lowermost dashed lines (horizontal solid line), only conduction (convection) occurs. Between these lines both conduction and convection are observed, depending hysteretically on whether $`ϵ`$ is increasing or decreasing. A tricritical point(TC) is shown at $`\mathrm{Re}0.11`$. For $`\mathrm{Re}\stackrel{<}{}0.11`$ the onset of convection is supercritical and hence non-hysteretic.

Within the convecting regime, the boundaries of the secondary bifurcations between $`m=8,7`$ and $`m=7,6`$ show upwards trends indicating that an analogous suppression occurs for the modes $`m=8`$ and $`m=7`$, even after they cease to be the primary mode of bifurcation. The suppression of the secondary modes increases more rapidly with $`\mathrm{Re}`$ than when they were the primary mode, with larger $`m`$ modes suppressed more rapidly than smaller $`m`$ modes. Inside the hysteresis loop of the primary bifurcation, we found as many as three possible distinct states, one being conduction. (See Figs. 2(b), (c), (g), and (h)). The entire portrait presented in Fig. 4 depends perhaps strongly on $`\alpha `$ and weakly on $`𝒫`$ and is likely to be more branched (contain more transitions) for larger $`\alpha `$. It should be possible to generate theoretically such a tree from the primitive equations, once they have been reduced to the form of Eq. 3.

We have shown that in our discrete-mode system, secondary bifurcations appear near codimension-two points. Using a generic Landau model, we have quantitatively analyzed the behavior of the primary bifurcations on either side of these points. We have found that the coefficients of the cubic and quintic terms of the Landau equation are strongly discontinuous at the codimension-two points. We expect that this scenario occurs over large regions of our parameter space and may well be a general property of discrete-mode pattern-forming systems such as TVF. Unlike TVF however, this system only has bifurcations between one-dimensional patterns with simple $`m`$-fold symmetries. The parameter space is characterized by branches that demarcate different stable and multistable domains. The surprizingly intricate bifurcation structure of this highly geometrically constrained instability is amenable to detailed quantitative analysis starting from the underlying electrohydrodynamic equations. As such, it is an ideal venue for experimental and theoretical studies of bifurcations in finite nonlinear pattern-forming systems.

We thank R. E. Ecke for helpful discussions and constructive comments. This research was supported by the Canadian NSERC and the U.S. DOE (W-7405-ENG-36).


text:
# Surface states and their possible role in the superconductivity of MgB2

\[

## Abstract

We report layer-Korringa-Kohn-Rostocker (LKKR) calculations for bulk and surface states (SFS) as well as the corresponding photoemission intensities of MgB2. Our theoretical results reproduce very well the recent angle resolved photoemission data by Uchiyama et al, cond-mat/0111152. At least two SFS are assigned. Consequences of SFS on the anisotropy of the upper critical fields and other properties in the superconducting state of grains in micro-powder samples are discussed.

PACS: 74.25.Jb, 74.70, 74.60.Ec, 73.-r

\]

Recently, the first angle-resolved photoemission spectroscopy (ARPES) measurements on a MgB2 single crystal were reported by Uchiyama et al. . They observed several dispersive states. Among them one has been ascribed to a surface state (SFS). The very existence of SFS and subsurface states for MgB2 as well as surface superconductivity have been predicted first in Refs. 2 and 3. However, there are significant quantitative differences between Ref. 2 and the ARPES data. The effects of SFS for the superconductivity in this challenging compound have not yet been discussed in detail. In particular, there is an unresolved puzzle with respect to the out-of plane anisotropy for the upper critical fields $`\gamma _H=H_{c2}^{ab}/H_{c2}^c`$. It is moderate for single crystals $`1.6\gamma _H3.5`$.  But large values $`6\gamma _H9`$ have been deduced indirectly from measurements on pure powder samples with small grains.  Enhanced superconductivity due to locally enhanced DOS (density of states) has been predicted  for idealized B- as well as Mg-terminated surfaces (BTS, MgTS) although in conflict with experimental data. We report on LKKR calculations for MgTS and BTS and compare them with the ARPES data. Furthermore we propose a scenario which involves essentially SFS to resolve some puzzles of general interest mentioned above.

The electronic DOS for the semi-infinite system MgB2 as well as the photoemission spectra, were calculated using a LKKR code . This code treats the photoemission process within the “one-step” approach. Dealing with the Green function formalism, it is capable to work with a complex potential, allowing the description of broadening effects due to the finite lifetime of quasi-particles. However, a set of input parameters must be provided phenomenologically in order to match experimental data as close as possible. As an input to the LKKR code, the real part of the one-electron effective complex potential was calculated self-consistently with the help of a LMTO (linear muffin-tin orbital) code. This preliminary calculation was done for the infinite MgB2 crystal with a basal plane lattice constant of $`a`$=3.085 Å and a ratio $`c/a`$=1.142. The von Barth-Hedin  correlation energy functional was used in connection with the local density approximation (LDA). The resulting Fermi energy $`E_F`$ 13.1 eV measured from the muffin-tin zero, while the muffin-tin radii were 1.711 Å and 0.988 Å for the Mg and B, respectively. A direct evaluation of the self-energy correction to be added to the LMTO potential is notoriously very difficult. Hence, we fixed it empirically. As for its imaginary part we chose different values for states below and above the Fermi level, -0.05 eV and -1 eV respectively. Its real part was chosen to be -1 eV and applied to the photo-electron final states only. Rigorously, the full complex self-energy correction is not constant in the energy intervals under consideration, but here such details can be omitted without losing much accuracy.

The geometry of the semi-infinite MgB2 system was simply modeled cutting the ideal crystal perpendicularly to the $`\widehat{c}`$ axis. Hence, no possible relaxation or reconstruction effects were taken into account. The topmost surface layer consists either of Mg atoms or B atoms. In the MgTS-case a work function $`\varphi =`$4.2 eV was used, in the BTS $`\varphi =`$6.1 eV, as suggested in Ref. . In both cases the surface potential barrier was modeled as a step potential, reflecting for occupied states and non reflecting, but refracting, for unoccupied ones. Its position $`d`$ was fixed to 0.987 Å with respect to the topmost layer, while its height is given by the sum $`E_F+\varphi `$. The energy of impinging photons at 28 eV and the linear polarization $``$ to the symmetry directions $`\mathrm{\Sigma }`$ and $`\mathrm{\Lambda }`$ were chosen in accordance with the experiment.

In Fig. 1 we show the DOS of the surface layer in

the case of a MgTS. The quasi parabolic feature with a minimum around -0.5 eV on the $`\overline{𝚪}`$ point is clearly a SFS, lying in a gap in the $`k_z`$-projected bulk band structure . It is $`s`$, $`p_z`$, $`d_{z^2}`$ in character, in a proportion roughly 3:2:1 at the $`\overline{𝚪}`$ point. Its characteristic decay length is 2.6 Å, whereas the distance between the Mg and B planes is 1.762 Å. In Fig. 1, the dispersion of bulk bands is also visible. The narrower stripes correspond to the quasi 2-dimensional $`\sigma `$ bonded px,y-bands, while the broad parabolic tape corresponds to the $`\pi `$ bonded pz-band.

SFS energies are much more sensitive to the values of the surface potential barrier positions $`d`$ than to the work functions $`\varphi `$. Placing the barrier at a larger $`d`$ value results in a downward energy shift (cf. Figs. 2 and 3). Thus, the situation depicted in Ref. 2, where the energy of the SFS of a MgTS is $``$ -2 eV at $`\overline{𝚪}`$, can be achieved with a surface potential barrier placed at $`d`$=1.32 Å above the topmost Mg layer. A different situation occurs in the BTS-case (cf. Fig. 3). Here the result of Ref. 2 is reproduced by a smaller $`d`$-value than our empirical one. This is probably due to the in-plane B-B distance which is 57% smaller than the Mg-Mg one, so that in a BTS slab calculation, the image plane position could have been chosen nearer to the surface than in the MgTS. The SFS which we found for a BTS in the considered energy range, runs parallel to the $`p_z`$ broad band near its upper edge. It is depicted by the solid line peak at a binding energy $``$ -0.8 eV in Fig. 3 calculated at $`𝐤_{}=0.7\overline{𝐌}`$. It decays towards the bulk with a characteristic length of 9.8 Å and exhibits pz symmetry. If steps are present on the real surface, yielding islands of Mg- and BTS, a single effective

(averaged) $`d`$ value could be justified to analyze the ARPES data. In view to have either a MgTS a BTS (or even both if we have a step-like surface), we calculated the ARPES spectra shown in Figs. 4 and 6.

We would like to stress that the crucial effective $`d`$ parameter must be fitted from experimental data, and the only available data  is shown in Fig. 5. Nice agreement with our calculations is obtained, if one superposes Figs. 4 and 6, suggesting that the measured sample presented both Mg- and BTS. According to our calculations the feature centered at $``$0.4 eV at $`\overline{𝚪}`$ is then to be assigned to a $`\pi `$-bonded Mg derived SFS in the case of MgTS. The parabolic feature centered at $``$2.7 eV is most likely, in our opinion, to be assigned to a $`\pi `$-bonded B $`2p_z`$-derived SFS arising from a BTS, in accordance with the calculations in Ref. 2. However, we admit that from the presented experimental data, an alternative assignment of this parabolic feature to the bulk B 2pz derived $`\pi `$-band edge cannot be completely excluded. An analogous consideration holds for the B $`\sigma `$ derived narrow bulk stripes, where a surface resonance state is found with our $`d`$ parameter value (cf. Fig. 3). Further ARPES experiments at variable photon energies should be performed to solve the full assignment problem. We emphasize that our interpretation mainly in terms of SFS does not reduce the importance of Ref. 1 which is the first experimental verification of the correct description of the electronic structure of MgB2 by the LDA. This is very important with respect to some times discussed strong many-body effects (correlation or non-adiabaticity) which would cause ) strong renormalizations of bulk and SFS as well. If strong correlation effects were present, they should

be observed first in the SFS, because the screening of Coulomb interactions is weaker near the surface rather than in the bulk. Taking into account the perfect agreement of LDA with the ARPES data, we conclude that such many-body effects should be of minor importance. With respect to surface superconductivity to be discussed below, the experimental determination of the electron-phonon (el-ph) coupling constant $`\lambda `$ from the temperature dependence of the ARPES spectra like in Be would be of great interest .

Now we address possible consequences of the SFS with respect to physical properties of various samples with emphasis on superconductivity and micro-powder samples. Since for typical samples the size in $`c`$-direction is about 10 nm only, about one quarter of the volume is affected by the presence of SFS. Since for MgTS the SFS shown in Fig. 1 is derived from orbitals with cylindrical symmetry and there is no contribution from orbitals with $`p_{x,y}`$ character, a very weak local el-ph interaction is expected. Thus it acts as a bad or non- superconducting layer in the sense of the proximity effect. Due to the interaction with the bulk states the resulting superconductivity near a MgTS will be weakened even for a clean MgTS, at variance with the view put forward in Ref. 3.

Thus the presence of such SFS might explain the absence of a third upper critical field $`H_{c3}`$  and the local suppression of superconductivity in relatively weak external magnetic field slightly above 2 Tesla . This was regarded as evidence for a small value of $`H_{c2}^c`$, i.e. a big anisotropy $`\gamma _H`$ 6 to 9. Anyhow, this would be in sharp contrast to the anisotropy caused by the Fermi velocities $`\gamma _v^\sigma `$ 6 to 7 which yields $`\gamma _H^\sigma `$4.2 to 5.5 as an upper limit  derived for the strongly interacting subsystem of $`\sigma `$-holes.  Due to the two-band character of MgB2  which involves also the weakly interacting nearly isotropic $`\pi `$-electrons the total $`\gamma _H`$ is further reduced till below 4 to 5, provided any order parameter anisotropy of the two different gaps can be neglected. A moderate $`\gamma _H3`$ is in accord with available single crystal data.

The interplay of superconducting bulk and SFS subsystems considered here contrasts the approach adopted in Ref.  where due to symmetry breaking at the surface mixing and pair breaking of $`\sigma `$\- and $`\pi `$\- states in the sense of a two-band model occurs. Notably, the absence or presence of superconducting SFS may affect the surface pinning. Bulk vortices will be attracted or repelled, respectively. This might explain the strong variations of pinning properties after applying pressure, which might affect the interfaces of the grains. Furthermore the presence of surface superconductivity on SFS will affect the tunneling data. In the context of frequently considered 3D multi-band models for MgB2 it is important to make a correct assignment of the observed superconducting gaps. In particular, in view of the ARPES interpretation we suggested above, the gap at about 5 meV observed by high-resolution (angle integrated) photoemission spectroscopy  might be ascribed to superconductivity on SFS.

Up to now two-dimensional (2D) electrons in metals were observed only in artificial layered structures. Let us consider in short the new physics which can be gained by studying generic 2D electrons at the surface of a 3D superconductor. As most interesting we consider the possibility for enhancement of the superconductivity of these 2D electrons due to increased local DOS or softening of local phonon modes at the surface. In this case in the framework of Ginzburg-Landau theory, for example, we have to introduce a 2D superconducting order parameter $`\mathrm{\Psi }_{2D}`$ in parallel to the bulk one $`\mathrm{\Psi }_{3D}`$; additionally the phenomenology could be complicated by necessity of taking into account pairing in different and weakly interacting bands. As a whole the surface enhancement of superconductivity could remind the twinning plane superconductivity observed in Sn and Nb. In this case the dependence of magnetization as a function of temperature and magnetic field can be similar to the fluctuation conductivity of the bulk. Roughly speaking, due to the proximity effect 2D conductivity can create a “frozen” 3D fluctuation. A systematic study of fluctuation magnetism in powdered samples with different grain sizes can give important hints in this direction. It will be very informative to perform investigation of fluctuation magnetism and to determine the mass anisotropy of fluctuation Cooper pairs using MgB2 microcrystals sintered under high pressure.

However it is more plausible that the influence of the surface is to reduce the local superconductivity: we have to take into account despite the possible Mg or B terminations of the crystal, also the adsorption of various impurities etc. In this case the influence of the surface can be traced out in the vortex pinning in the superconducting state. Indeed, strong asymmetric hysteresis loops were observed after the grinding of a bulk sample into fine powder. Surface pinning was reported for a powder sample pressed into a pellet without sintering.  Both possibilities, enhancement and suppression of superconductivity yield additional surface pinning. Finally, we argue that 2D electronic states can give an important contribution to the current dependent correction of the work function of the superconducting metal, the so called Bernoulli potential.

We have shown that the available ARPES data can be well described by the LDA. In particular, there is clear evidence for the presence of at least two SFS. Since the SFS depend very sensitively on the nature of the terminated layer, comparing our calculated intensities with the measured ones, we suggest that the studied single crystal contained both types of terminated layers. Therefore the investigation of especially prepared surfaces with only one predominant surface layer is of interest. In this context the theoretical study of steps between B and Mg terminated layers would be interesting. Since the SFS are generated by the large gapped region in the BZ, similar phenomena might be expected for other layered superconductors such as NbSe2 and even heavy fermion compounds URu2Si2 or borocarbides  where also no $`H_{c3}`$ has been observed. Taking into account all those possibilities mentioned above one can expect that comprehensive studies of MgB2 could generate interesting physics of the surface/interface of superconductors.

V.D.P.S. is grateful to A.N. Yaresko and S.V. Halilov who respectively provided the LMTO and layer-KKR codes used in the present study. We are indepted to S. Tajima for providing us with the ARPES and other single crystal data prior to publication. We thank S. Tajima, H. Uchiyama, A. Kordyuk, G. Fuchs, S. Shulga, M. Richter, G. Paasch, J. Fink, H. Eschrig, C. Laubschat, and H. Rosner for discussions. V.D.P.S. was supported by the Deutsche Forschungsgemeinschaft (project La-655/7-1).