text:
# I Introduction

## I Introduction

The interaction of traveling solitary pulses (which, in particular but not necessarily, may be solitons, that we realize here as pulses in nonlinear media maintaining a steady shape while propagating) with various local imperfections or pinned dynamical excitations is a problem of fundamental importance. For direct experimental observation, the most straightforward case is the propagation of a pulse in an optical fiber (linear or nonlinear), in which a strong localized nonlinear defect represents an inserted piece of a dispersion-shifted fiber (DSF) , that has finite nonlinearity and negligible dispersion. This configuration, which may be quite useful for optical telecommunications , is described by the following version of the nonlinear Schrödinger (NLS) equation,

$$iu_z\frac{1}{2}\beta _2u_{\tau \tau }+\frac{i}{6}\beta _3u_{\tau \tau \tau }+\gamma \left|u\right|^2u+\mathrm{\Gamma }\delta \left(z\right)\left|u\right|^2u=0,$$
(1)
where $`u(z,\tau )`$ is a local amplitude of the electromagnetic wave, $`z`$ is the propagation coordinate,

$$\tau tz/V_0$$
(2)
is the so-called reduced time, $`V_0`$ being the group velocity of the carrier wave, $`\beta _2`$ and $`\beta _3`$ are coefficients of the second-order and third-order group-velocity dispersion (GVD) in the fiber, and $`\gamma `$ is the nonlinearity coefficient of the host (system) fiber. The strength of the localized nonlinear defect is $`\mathrm{\Gamma }\gamma _{\mathrm{DSF}}L`$, where $`\gamma _{\mathrm{DSF}}`$ and $`L`$ are the nonlinearity coefficient and actual length of the above-mentioned finite-length DSF inserted into the host fiber at the point $`z=0`$; this insertion may be represented by the delta-function in Eq. (1), as the DSF’s dispersion is negligible . Note that the nonlinear coefficient $`\gamma `$ is always positive, as it is induced by the Kerr effect in the optical fiber, that always has the sign corresponding to self-focusing . On the other hand, the GVD coefficient $`\beta _2`$ may be either positive or negative, which corresponds, respectively, to normal and anomalous dispersion . The third-order dispersion (TOD) coefficient $`\beta _3`$ is frequently neglected, unless $`\beta _2`$ is very small or if the pulse is very narrow in the $`\tau `$-domain. However, TOD will play an essential role in the present work. Note that $`\beta _3`$ is usually positive in optical fibers .

A qualitatively different version of the model containing the local defect pertains to the case when it describes the interaction of the propagating modes (pulses) with a localized defect in the form of a pinned breather . These are spatially localized, time-periodic dynamical states, which are ubiquitous in nonlinear physical systems ranging from quasi-one dimensional polymers  and charge-density-wave materials (e.g., metal-halogen electronic chains ) to Josephson ladders . Scattering of a moving solitary pulse on a breather is an important process in many physical, chemical and biological systems which combine charge, spin and energy localization with transport of these quantities by pulses. Detailed understanding of this type of the scattering will not only yield valuable information on the dynamical properties of breathers, such as their mobility and stability, but also reveal the dynamical response of materials in which the breathers are excited. In particular, the interaction of a narrow optical pulse with a pre-existing electron-phonon breather in a lattice chain is a process of considerable fundamental and technological interest, since it may enhance optical nonlinearities of the material, and thus the efficiency of the second- and third-harmonic generation. Indeed, there is some evidence of such an effect in conjugated polymers such as polyenes .

Another relevant example of the breathers is provided by the rotational dynamics of certain chemical groups in a molecular chain . Analysis of the pulse scattering on them, in conjunction with recent developments in ultrafast (femtosecond) spectroscopy  and chemistry , will enable an efficient use of the rich photophysics of functional optical and electronic materials. Damage tracks in certain mica minerals and sputtering on crystal surfaces have also been attributed to moving breathers . In addition to mixed-valence transition-metal complexes , there is experimental evidence for localized breather-like states in antiferromagnetic chains , as well as in the above-mentioned Josephson junction arrays and ladders .

The one-dimensional (1D) model introduced and considered in this work is a first step in a systematic study of the pulse-breather scattering, the eventual objective being to elucidate the interaction of pulses and breathers in 2D and 3D nonlinear dynamical lattices. However, the formulation and consideration of the corresponding multidimensional models is a complicated issue, which is beyond the scope of this work.

The simplest generalization of Eq. (1) for the case when the local nonlinearity is induced by a small-size breather oscillating at a frequency $`\omega `$ is

$$iu_z\frac{1}{2}\beta _2u_{\tau \tau }+\frac{i}{6}\beta _3u_{\tau \tau \tau }+\gamma \left|u\right|^2u+\mathrm{\Gamma }\delta \left(z\right)\mathrm{cos}^2\left(\omega t\right)\left|u\right|^2u=0.$$
(3)
Equation (3) implies that oscillations of the breather modulate the strength of the corresponding local nonlinear defect, but do not change its sign, which is a natural assumption in the case of the cubic nonlinearity. Note that, despite the difference between the time $`t`$ and the reduced time $`\tau `$, see Eq. (2), it makes no difference what definition of time is used in the argument of $`\mathrm{cos}^2`$ in Eq. (3), as $`\tau t`$ at $`z=0`$.

Both equations (1) and (3) conserve the norm of the wave field (which has the physical meaning of energy in the application to fiber optics ),

$$E=_{\mathrm{}}^+\mathrm{}\left|u(\tau )\right|^2𝑑\tau .$$
(4)
Additionally, the model (1) including the static defect conserves the field momentum,

$$P=i_{\mathrm{}}^+\mathrm{}u_\tau u^{}𝑑\tau ,$$
(5)
the asterisk standing for the complex conjugation.

The simplest and perhaps most fundamental version of both models (1) and (3) is the one with all the nonlinearity concentrated at the point $`z=0`$, while the bulk nonlinearity is negligible, i.e., $`\gamma =0`$. It appears, rather surprisingly, that, unlike the more technically involved case of the fully nonlinear model with $`\gamma >0`$, which was studied in a part in Ref. , in the case $`\gamma =0`$ the interaction of the moving pulse with the nonlinear defect has not been yet considered. Therefore, the detailed study of the pulse-defect interaction in the model with $`\gamma =0`$ is the first objective of the present work. Then, we will also consider the full nonlinear model with $`\gamma 0`$, concluding that, as a matter of fact, the most interesting results can be obtained just in the model with the linear host medium, $`\gamma =0`$.

In the linear medium, it is natural to take an incident pulse in the form of a Gaussian, which is an eigenmode of the linear system. The problem of the interaction of a Gaussian with the delta-like nonlinear defect can be partially solved in an analytical form, which is done below in section 2. In fact, the analytical result can be obtained in an advanced form for the case of the static defect with $`\omega =0`$, which corresponds to Eq. (1) and, as it was explained above, has direct application to fiber optics. A basic result produced by the analytical consideration is that the initial Gaussian pulse, passing through the nonlinear defect, generates an infinite series of strongly overlapped Gaussian pulses. At this stage of the analysis, TOD becomes a crucially important ingredient of the model: we show that it lends each Gaussian pulses its own velocity, which will bring about eventual separation of the pulses. In section 2, the velocity generated by TOD is found in an analytical form, which yields an exact result for an initial stage of the evolution. The splitting of the solution into an array of pulses is considered in detail in section 3.

In section 4 we display results of the fully numerical solution for the case when the incident Gaussian pulse interacts with the dynamical defect \[the breather, see Eq. (3)\]. In this case ($`\omega 0`$), the eventual pattern is much less regular, consisting of separated pulses with random shapes, rather than Gaussian-like ones. However, if the frequency $`\omega `$ is very large, one should expect that the dynamical defect described by the last term in Eq. (3) may be replaced by its averaged static counterpart. In accordance with this expectation, the numerical computations demonstrate return to a regular pattern for very large values of $`\omega `$.

In section 5 we present a solution for the most general model, which includes the bulk nonlinearity at $`z0`$, i.e., $`\gamma 0`$ in Eq. ( 3). In this case, it is natural to take the incident pulse as a soliton of the corresponding homogeneous NLS equation, rather than a Gaussian pulse, and the interaction may only be simulated numerically. Results of simulations of the interaction of the soliton with the defect turn out to be quite different from those in the case when the host medium was linear: the soliton is found to be much more resistant to the action of the local nonlinear defect than the Gaussian pulse in the linear medium. Unless the defect is extremely strong, the pulse remains virtually intact, the TOD term separating it from small-amplitude wave packets generated by the local defect. However, in the absence of the TOD term, the effect of the perturbation may accumulate and destroy the soliton. Thus, the third-order dispersion plays a crucially important role in both versions of the present model, i.e., for the linear and nonlinear host medium.

If the defect is very strong, it may split the soliton into two. The results obtained for the model with the nonlinear host medium are not sensitive to the value of the defect’s intrinsic frequency $`\omega `$, being essentially the same for the static defect and its dynamical counterpart. Moreover, the soliton is found to be stable against the action of the local defect irrespective of the relative sign of the bulk and local nonlinearities. These results, obtained for the case when the input signal is a soliton, rather than a linear pulse, is another manifestation of the well-known general principle, according to which solitons are very robust eigenmodes of nonlinear media.

## II Analytical consideration of the interaction between a Gaussian pulse and a nonlinear defect.

### A The pulse

We start the consideration with the simplest version of the model, viz., Eqs. (1) or (3) with $`\gamma =\beta _3=0`$. At $`z0`$, we thus have a linear Schrödinger equation, which gives rise to the well-known exact solution in the form of a Gaussian pulse (which is termed a coherent state in quantum mechanics),

$$u(z,\tau )=\frac{A_0}{\sqrt{12i\beta _2cz}}\mathrm{exp}\left[\frac{c\left(1+2i\beta _2cz\right)}{1+\left(2\beta _2cz\right)^2}\tau ^2\right],$$
(6)
where $`c>0`$ and $`A_0`$ are arbitrary real constants. The constant $`c`$ determines the width of the Gaussian pulse and its chirp (i.e., the imaginary part of the coefficient in front of $`\tau ^2`$ in the argument of the exponential, which shows a slope of the local frequency across the pulse ).

Note that the solution (6) has zero velocity in the present reference frame \[which is defined by Eq. (2)\]. An exact solution for a moving pulse can be generated from Eq. (6) by the action of a boost (Galilean transform)

$$u(z,t)u(z,\tau sz)\mathrm{exp}\left(\frac{is^2}{2\beta _2}z\frac{is}{\beta _2}\tau \right),$$
(7)
where $`s`$ is a real velocity-shift parameter. It is important to notice that, for the boosted pulse, the values of its velocity, momentum (5) and energy (4) are related, irrespective of the particular form of the pulse, in a simple way,

$$s=\beta _2\frac{P_{\mathrm{pulse}}}{E_{\mathrm{pulse}}}.$$
(8)
### B Pulse acceleration by the third-order dispersion and numerical verification

Before proceeding to the consideration of the passage of the pulse (6) through the nonlinear defect, we need to understand how a free pulse will move under the action of the TOD term in Eqs. (1) or ( 3). The consideration of this issue will help to understand how an array of pulses splits in the presence of the TOD term.To this end, a solution is sought for by means of the Fourier transform, which yields the following integral representation for it:

$$u(z,\tau )=\frac{A_0}{2\sqrt{\pi c}}_{\mathrm{}}^+\mathrm{}\mathrm{exp}\left[\frac{\omega ^2}{4c}i\omega \tau +i\left(\frac{1}{2}\beta _2\omega ^2\frac{1}{6}\beta _3\omega ^3\right)z\right]𝑑\omega $$
(9)
\[setting $`z=0`$, this expression goes over into Eq. (6) taken at $`z=0`$\]. In the limit of large values of $`z`$, the integral (9) is dominated by a contribution from a vicinity of the stationary-phase point, which is

$$\omega _0\frac{\tau }{\beta _2z}\frac{i\tau }{2c\beta _2^2z^2}+\frac{\beta _3\tau ^2}{2\beta _2^3z^2}.$$
(10)
In particular, a contribution of the stationary-phase point (10) to the phase of the pulse is

$$\varphi (z,\tau )=\frac{\tau ^2}{2\beta _2z}\frac{\beta _3\tau ^3}{6\beta _2^3z^2}.$$
(11)
Direct simulations show that, under the action of TOD, the initial Gaussian pulse is gradually destroyed, generating a long “tail”, while a part of the wave packet may still be interpreted as a surviving pulse. Using the expression (11), one can find a relation between the momentum \[see Eq. (5)\] and energy of the pulse-like part of the wave packet,

$$P_{\mathrm{pulse}}=\frac{\beta _3c}{2\beta _2}E_{\mathrm{pulse}}.$$
(12)
Note that, as the net momentum of the wave field must be conserved, a “recoil” momentum $`P_{\mathrm{pulse}}`$ is carried away by the above-mentioned tail. The comparison of Eqs. (8) and (12) yields an analytical prediction for the asymptotic value of the velocity shift acquired by the surviving pulse under the action of the TOD term in the limit $`z\mathrm{}`$:

$$s=(1/2)\beta _3c.$$
(13)
To this end, a natural definition of the position of the wave-packet’s center of mass is adopted,

$$\tau _cE^1_{\mathrm{}}^+\mathrm{}\tau \left|u(\tau )\right|^2𝑑\tau ,$$
(14)
where $`E`$ is the net energy defined by Eq. (4), and the derivative $`d\tau _c/dz`$ may be regarded as a velocity of the pulse.

In the case of the free propagation, one can derive an exact relation

$$\frac{d\tau _c}{dz}=\beta _2\frac{P}{E}\frac{1}{2}\beta _3E^1_{\mathrm{}}^+\mathrm{}\left|u_\tau \right|^2𝑑\tau ,$$
(15)
where $`P`$ and $`E`$ are the net momentum and energy of the wave field. Comparing this to Eq. (8), we conclude that, in the absence of TOD, the velocity (15) is identical to the boost parameter (in particular, $`d\tau _c/dz=0`$ if $`P=0`$); however, this identity is broken by TOD, that is why acceleration of the pulse by TOD is observed. If the TOD term is treated as a small perturbation, one can analytically calculate the second term on the right-hand side of Eq. (15), using the expression (6), which is an exact solution in the absence of TOD. As a result, we find the velocity of the pulse in the case when the net field momentum vanishes, $`P=0`$, which is true for the initial Gaussian pulse (6):

$$\frac{d\tau _c}{dz}=\frac{1}{3}\beta _3c$$
(16)
\[note that this expression does not contains neither $`z`$ nor $`\beta _2`$, despite the fact the solution (6), used for the calculation of the right-hand side of Eq. (16), does depend on $`z`$ and $`\beta _2`$\]. The expression (16) shows that the velocity lent to the pulse by TOD linearly depends on the pulse’s parameter $`c`$, hence initially overlapping pulses with different values of $`c`$ \[see Eq. (21) below\] are expected to separate under the action of TOD.

The analytical result (21) was checked against direct numerical simulations of the linear version of Eq. (1). For instance, in the case $`\beta _2=1`$, $`\beta _3=0.1`$, and $`c=1`$, it was found that the pulse was indeed moving at a constant velocity, the value of which exactly coincided with that given by Eq. (16), in the interval $`0<z<25`$. At larger values of $`z`$, the absolute value of the velocity decreases, which can be explained by the fact that the TOD term essentially alters the shape of the pulse at that later stage of the evolution.

### C Passage of the pulse through the nonlinear defect

Our next aim is to consider transformation of the Gaussian pulse passing the nonlinear defect. Obviously, in an infinitesimal vicinity of the defect (as $`|z|0`$), only the first and last terms should be kept in Eqs. ( 1) and (3), which yields a simplified equation,

$$\frac{u}{z}=i\mathrm{\Gamma }\delta \left(z\right)\mathrm{cos}^2\left(\omega \tau \right)\left|u\right|^2u.$$
(17)
To solve Eq. (17), we represent the solution as $`u(z)a(z)\mathrm{exp}\left[i\varphi (z)\right]`$, with real amplitude $`a`$ and phase $`\varphi `$. Substituting this into Eq. (17), one immediately finds that $`a/z=0`$, and

$`{\displaystyle \frac{\varphi }{z}}=\mathrm{\Gamma }\delta \left(z\right)\mathrm{cos}^2\left(\omega \tau \right)a^2.`$

A solution to the latter equation is obvious,

$$\varphi (z=+0,\tau )\varphi (z=0,\tau )=\mathrm{\Gamma }a^2\left(\tau \right)\mathrm{cos}^2\left(\omega \tau \right),$$
(18)
where we take into regard that $`a`$ may be a function of $`\tau `$.

Thus, we take the input pulse at the point $`z=0`$ in the general form \[cf. the expression (6)\],

$$u(z=0)=A_0\mathrm{exp}\left[\left(c_0+ib_0\right)\tau ^2\right],$$
(19)
where $`c_0>0`$ determines the initial width of the pulse, and $`b_0`$ is its initial chirp. The substitution of the expression (19) into Eq. (18) yields the form of the pulse appearing after the passage of the nonlinear defect:

$$u(z=+0,\tau )=A_0\mathrm{exp}\left[\left(c_0+ib_0\right)\tau ^2+i\mathrm{\Gamma }A_0^2\mathrm{exp}\left(2c_0\tau ^2\right)\mathrm{cos}^2\left(\omega \tau \right)\right].$$
(20)
Further analytical consideration for the general case of the dynamical defect (breather), with $`\omega 0`$, is extremely cumbersome. Therefore, in the rest of this section and in the next one, we focus on the static case, $`\omega =0`$. The aim will be to realize a result of the further evolution of the transformed pulse (20), governed by the linear equation (1) Including the TOD term) at $`z>0`$. To this end, we notice that the expression (20) can be expanded into an infinite series:

$$u(z=+0)=A_0\mathrm{exp}\left(ib_0\tau ^2\right)\underset{n=0}{\overset{+\mathrm{}}{}}\left(n!\right)^1\left(i\mathrm{\Gamma }A_0^2\right)^n\mathrm{exp}\left[\left(1+2n\right)c_0\tau ^2\right].$$
(21)
Comparing it to the exact fundamental-pulse solution (6), we conclude that the expression (21), if considered as an initial condition to the linear equation (1) with $`\gamma =\beta _3=0`$, gives rise to a superposition of an infinite number of Gaussians with the values of the width constant $`c_n=\left(1+2n\right)c_0`$, and with the common initial value $`b_0`$ of the chirp. The evolution of the wave packet (21) can be presented in a relatively simple form in the case when the initial chirp is absent, $`b_0=0`$ (and TOD is neglected, $`\beta _3=0`$):

$$u(z,\tau )=A_0\mathrm{exp}\left(ib_0\tau ^2\right)\underset{n=0}{\overset{+\mathrm{}}{}}\frac{\left(i\mathrm{\Gamma }A_0^2\right)^n}{n!\sqrt{12i\beta _2c_0\left(1+2n\right)z}}\mathrm{exp}\left[\frac{c_0\left(1+2i\beta _2c_0\left(1+2n\right)z\right)}{1+\left(1+2n\right)^2\left(2\beta _2c_0z\right)^2}\tau ^2\right]$$
(22)
\[in the case $`\omega 0`$, Eq. (20) shows that the Gaussian in each term of the initial series (21) with $`n0`$ is additionally multiplied by $`\left[\mathrm{cos}(\omega t)\right]^{2n}`$, which will make the subsequent result much more complex than that given by Eq. ( 22)\].

Thus, Eq. (22) gives an exact solution to the nonlinear model equation (1) in the case when the incident pulse has no chirp and $`\gamma =\beta _3=0`$. However, while taking the input pulse to be chirpless, and disregarding the bulk nonlinearity ($`\gamma =0`$) are quite acceptable assumptions, the TOD term may not be neglected, as, without this term, the exact solution (22) remains strongly degenerate. Indeed, centers of all the Gaussian pulses, the superposition of which constitutes this solution, exactly coincide, staying at $`\tau =0`$ \[note that this degeneracy is not lifted if the input pulse has nonzero chirp, nor if the bulk nonlinearity ($`\gamma 0`$) is added\]. On the other hand, Eq. ( 16) shows that the TOD term will lend each pulse its own velocity, depending on the initial width parameter $`c`$ of the pulse. Obviously, this will eventually split the superposition (22) of the overlapping Gaussians into an array of separating pulses.

However, the analytical prediction (16) for the TOD-induced velocity shift of each pulse is only valid for finite values of $`z`$. Therefore, to find an actual shape of the evolving wave train, it is necessary to solve numerically the linear Schrödinger equation with the TOD term,

$$iu_z\frac{1}{2}\beta _2u_{\tau \tau }+\frac{i}{6}\beta _3u_{\tau \tau \tau }=0,$$
(23)
with the initial condition in the form (20). Results generated by Eq. (23) are presented in the next section.

## III Generation of the wave train by the third-order dispersion in the case of the static nonlinear defect

As Eq. (23) is linear and has constant coefficients, it can be solved by Fourier transform. Thus, the numerical part of the solution amounts to the computation of the Fourier transform

$`u_0(\omega ){\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{exp}(i\omega \tau )u(z=+0,\tau )𝑑\tau `$

for the initial configuration (20), and subsequent computation of the inverse Fourier transform

$`u(z,\tau )(2\pi )^1{\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{exp}(i\omega \tau )u(z,\omega )𝑑\tau ,`$

where, as it immediately follows from Eq. (23),

$`u(z,\omega )=u_0(\omega )\mathrm{exp}\left[i\left({\displaystyle \frac{1}{2}}\beta _2\omega ^2z{\displaystyle \frac{1}{6}}\beta _3\omega ^3z\right)\right].`$

In Fig. 1 we display the profiles of $`|u(\tau )|`$, computed at the points $`z=10,\mathrm{\hspace{0.17em}20},\mathrm{\hspace{0.17em}30},`$ and $`40`$ for a case when the static nonlinear defect is weak, having $`\mathrm{\Gamma }=0.1`$. The figure demonstrates that, in accordance with the analysis presented in the previous section, the initial wave packet tends to split into an array of regular Gaussian-like pulses. The splitting actually takes place for larger values of the nonlinear-defect’s strength, as is shown in Fig. 2, which pertains to $`\mathrm{\Gamma }=1`$. For the same case, the eventual shape of the pulse array is shown in more detail in Fig. 3.

The result for a still stronger nonlinear defect, with $`\mathrm{\Gamma }=10`$, is displayed in Fig. 4. In this case, the splitting into pulses takes a violent character, with appearance of huge gradients and formation of a very sharp front, the latter effect being accounted for by the interplay of the strong local nonlinearity and TOD. These results may be explained by the fact that, as follows from Eq. (21), the number $`n_{\mathrm{max}}`$ of the largest-amplitude pulse in the series grows $`\mathrm{\Gamma }`$ with increase of $`\mathrm{\Gamma }`$, hence the width $`W`$ of an individual pulse decreases $`1/\sqrt{\mathrm{\Gamma }}`$, which gives rise to the large gradients revealed by the subsequent evolution. We also note that, although the TOD coefficient $`\beta _3`$ remains small, the size of the TOD term in Eqs. (1) and (3) grows as $`1/W^3`$ with the decrease of the pulse’s width. This explains the formation of the abrupt front under the action of the asymmetric TOD term.

## IV Generation of the pulse array by a localized breather

An example of the splitting of the initial Gaussian as a result of its interaction with the dynamical defect (with the same strength $`\mathrm{\Gamma }=1`$ as in the case shown in Fig. 2, and the frequency $`\omega =1/2`$) is shown in Fig. 5. The splitting is displayed in more detail by blowups collected in Fig. 6.

Comparison of Figs. 2 and 5 shows that the static defect and the dynamical one with a moderate value of the frequency produce quite similar results. Taking larger values of the frequency strongly changes the situation: as is shown in Fig. 7 and in the accompanying blowup (Fig. 8), the same value of the defect strength as in the cases displayed in Figs. 2 and 5, i.e., $`\mathrm{\Gamma }=1`$, but combined with $`\omega =10`$, gives rise to an essentially more disordered pattern. Note that this pattern is disordered in a way essentially different from that observed as a result of the action of a strong static defect ($`\omega =0`$) , cf. Fig. 4. In particular, strong asymmetry of the pattern and sharp fronts are not found in the present case, and the local gradients are not as huge as in Fig. 4. The absence of those features in the present case is easy to understand, as they may only be generated by a large value of $`\mathrm{\Gamma }`$, as explained above.

Keeping to increase $`\omega `$ at a fixed value of $`\mathrm{\Gamma }`$, we have concluded that the situation shown in Figs. 7 and 8 is rather similar to that observed at $`\mathrm{\Gamma }=1`$ and $`\omega =100`$ (not shown here). On the other hand, it is obvious that, if the frequency is extremely large, one should be able to replace $`\mathrm{cos}^2(\omega t)`$ in Eq. (3) by its average value $`1/2`$, thus reverting to Eq. (1) for the static nonlinear defect. To verify this argument, in Fig. 9 we present the results of the numerical computations for $`\mathrm{\Gamma }=1`$ and $`\omega =1000`$, which are supported by the blowup shown in Fig. 10. From these pictures, it is evident that, in the case of extremely large $`\omega `$, the situation is indeed nearly the same as in the case of the static nonlinear defect, cf. Figs. 2 and 3.

## V Transformation of a soliton by the local defect in the nonlinear host medium

In all the cases considered above, the input pulse was taken in the form of a Gaussian, since the host medium in which this pulse propagated before the collision with the nonlinear defect was linear, the Gaussian being its eigenmode. The situation is completely different in the case when the host waveguide is itself nonlinear, as in that case the incoming pulse must be a $`\mathrm{sech}`$ soliton . The action of the localized nonlinearity on the soliton is described by the same general expression (18) as above; however, the result cannot be interpreted in such a straightforward manner (even if $`\omega =0`$) as it was done above for the case of the input Gaussian pulse, see Eq. (21). In fact, in the case of the nonlinear host medium all the analysis following the use of expression (18) for the defect-induced phase change of the input soliton can only be performed numerically.

As is commonly known, in the homogeneous part of the nonlinear model (1) or (3) (at $`z0`$), the soliton may only exist if the second-order dispersion is anomalous while the nonlinearity is self-focusing, or vice versa , i.e., only if $`\beta _2\gamma <0`$ . In the simulations, we fixed $`\beta _2+1`$, and took different negative values of $`\gamma `$, see below. The input soliton was always taken as a fundamental one with a fixed width corresponding to these values of the parameters:

$$u\left(z=0,\tau \right)=|\gamma |^{1/2}\mathrm{sech}\tau .$$
(24)
If the TOD term is neglected ($`\beta _3=0`$), the nonlinear defect produces a strong perturbation around the soliton. However, the TOD term helps to separate the soliton and the perturbation, as is illustrated by Fig. 11 for the case $`\beta _3=0.3`$. This figure shows a large $`\tau `$-domain, in order to demonstrate the evolution over a large propagation distance; as is seen, the perturbation spreads out indefinitely, while the soliton remains essentially intact.

The comparison of Fig. 11 with Fig. 2, which shows the interaction of the Gaussian pulse in the linear host medium with the local nonlinear defect that has the same strength, $`\mathrm{\Gamma }=1`$, suggests a conclusion that the soliton in the nonlinear host medium is much more resistant to the action of the nonlinear defect than the Gaussian pulse in the linear host medium. In fact, this conclusion is strongly supported by many other simulations with different values of the parameters.

Another difference of the present case from that for the Gaussian pulse in the linear host medium is that the result of the action of the nonlinear defect on the soliton is less sensitive to the defect’s intrinsic frequency. For instance, if the defect is dynamical with $`\omega =1`$, while the other parameters take the same values as in the case shown in Fig. 11, the evolution (not shown here) seems nearly the same as in Fig. 11, with a difference that the generated perturbation takes a somewhat larger portion of energy from the soliton (as a result, the soliton reappears after the collision with the amplitude $`0.8`$, to be compared with the amplitude $`1`$ in Fig. 11). Nevertheless, in this case too, the perturbation remains rather small and it does not strongly affect the soliton.

With further increase of the dynamical defect’s frequency, the perturbation again becomes smaller, and the general picture is reverting to that corresponding to the static defect with $`\omega =0`$. This trend can be easily explained by the self-averaging of the dynamical defect similar to that observed above in the case of the linear host medium. However, unlike that case, where the self-averaging was evident only for extremely large values of the frequency, $`\omega 1000`$ (see Fig. 9), in the present case simulations (not shown here) clearly show that $`\omega =10`$ is sufficient for the self-averaging to manifest itself.

The above results were obtained for the cases with $`\mathrm{\Gamma }>0`$ and $`\gamma <0`$, which implies that the localized nonlinearity is self-defocusing if the bulk nonlinearity is self-focusing, or vice versa. It is also interesting to consider the case when both nonlinearities have the same sign, i.e., $`\mathrm{\Gamma }<0`$. An example is displayed in Fig. 12 for $`\mathrm{\Gamma }=1`$ (this example pertains to the dynamical defect with $`\omega =1`$, but the results are virtually the same for the static defect with $`\omega =0`$). In this case, the action of the defect produces a weaker effect, as the soliton reappears after the interaction with the amplitude $`1`$, to be compared with the above-mentioned value $`0.8`$ of the soliton’s amplitude found for the same values of parameters but $`\mathrm{\Gamma }=+1`$. This difference is quite natural, as the local defect which has the same sign of the nonlinearity as in the host medium is acting to compress the soliton stronger, while the local defect with the opposite sign of $`\mathrm{\Gamma }`$ was acting to stretch the soliton, which produces a more destructive effect, helping some radiation (wave packets) to escape from the weakened soliton pulse.

Taking a much stronger nonlinear defect (with $`\mathrm{\Gamma }=11`$, see Fig. 13), we conclude that, naturally, the soliton loses a larger part of its energy to the generation of the perturbation separating from it. Nevertheless, the soliton survives even in this case (cf. Fig. 7, which shows a complete and fast chaotization of the wave field in the case of the collision of the Gaussian pulse with a strong defect, having $`\mathrm{\Gamma }=10`$, in the linear host medium). In the cases when $`\mathrm{\Gamma }`$ is large but negative (not shown here), when the nonlinearity of the defect has the same sign as in the host medium, the soliton’s losses are smaller than in the case displayed in Fig. 13, which is quite similar to what is shown in Fig. 12 for a moderately strong defect with $`\mathrm{\Gamma }<0`$.

## VI Conclusion

We have introduced a model describing collision of a pulse in a linear or nonlinear waveguide with a strong nonlinear local defect, that may be either static or breather-like. The model with the static defect should directly apply to an optical pulse in a long fiber-optic link with an inserted section of a nonlinear (dispersion-shifted) fiber of an arbitrary length. On the other hand, a local breather, which gives rise to the nonlinear defect affecting the propagation of narrow optical pulses, may be realized in a molecular chain with electron-phonon coupling.

In the case when the host waveguide is linear, the pulse was naturally taken as a Gaussian. A result of its interaction with the nonlinear defect was found analytically, amounting to its transformation into a “lump” consisting of an infinite number of overlapping Gaussian pulses. Further evolution of the lump in the linear medium is generated by the corresponding Fourier transform. An important ingredient of the medium is the third-order dispersion, that splits the lump into an array of individual pulses; the velocity shift lent to an initially Gaussian pulse by the TOD term was found in an analytical form. The influence of the intrinsic frequency, in the case when the defect is a breather, was also investigated.

The full numerical solution for the model in which the host medium is nonlinear, and the input pulse is taken as a fundamental soliton shaped by this medium, produces a result quite different from that for the linear host medium: the soliton is much more resistant to the action of the local nonlinearity. If the local defect is not very strong, the soliton remains essentially intact, the third-order-dispersion separating the soliton and small wave packets generated by the collision. Even if the local nonlinearity is very strong, the soliton survives, losing a limited part of its energy.

Beyond straightforward application of the results reported above to fiber-optic links, the model introduced in this work, and current developments in experimental techniques, such as ultrafast spectroscopy  of optical and electronic materials , the results presented here may help to understand the dynamics of scattering of pulses on static and breather-like defects in many physical systems, e.g., conjugated polymers , Josephson ladders , and coupled electron-vibron lattice systems .

## Acknowledgements

B.A.M. appreciates hospitality of the Center for Nonlinear Studies at the Los Alamos National Laboratory and of the Department of Physics at the University of Athens (Greece). Work at Los Alamos was supported by the U.S. Department of Energy.

## Figure Captions

Fig. 1. The evolution of the field $`|u(\tau )|`$ after the interaction of the input Gaussian pulse with the nonlinear defect, in the case $`\omega =\gamma =0`$, $`\beta _2=1`$, $`\beta _3=0.1`$, and $`\mathrm{\Gamma }=0.1`$.

Fig. 2. The same as in Fig. 1 for a stronger nonlinear defect, with $`\mathrm{\Gamma }=1`$.

Fig. 3. Blowups of the field pattern from the last panel in Fig. 2, clearly showing the formation of an array of separated pulses.

Fig. 4. The same as in Fig. 3 for a very strong static nonlinear defect, with $`\mathrm{\Gamma }=10`$.

Fig. 5. The evolution of the field $`|u(\tau )|`$ after the interaction of the input Gaussian with the nonlinear dynamical defect, in the case $`\omega =1/2`$ , $`\gamma =0`$, $`\beta _2=1`$, $`\beta _3=0.1`$, and $`\mathrm{\Gamma }=1`$. Except for $`\omega `$, these parameters are the same as in the case shown in Fig. 2.

Fig. 6. Blowups of segments of the panels from Fig. 5, showing the gradual splitting of the wave packet into an array of pulses.

Fig. 7. The same as in Fig. 5, but with $`\omega =10`$.

Fig. 8. The same as in Fig. 6, but with $`\omega =10`$.

Fig. 9. The same as in Fig. 7, but with $`\omega =1000`$.

Fig. 10. The same as in Fig. 6, but with $`\omega =1000`$.

Fig. 11. The evolution of the field $`|u(\tau )|`$ over a very large propagation distance in a large temporal domain after the interaction of the input soliton (24) with the nonlinear static defect, in the case $`\gamma =\beta _2=1`$, $`\beta _3=0.3`$, and $`\mathrm{\Gamma }=1`$, $`\omega =0`$. It is evident that the local defect initially creates a large perturbation around the soliton; however, the perturbation is diffused away under the combined action of the second- and third-order dispersions, and a robust soliton reappears.

Fig. 12. The same as in Fig. 11, but with $`\omega =1`$ and $`\mathrm{\Gamma }=1`$.

Fig. 13. The same as in Fig. 12, but with $`\beta _3=0.5`$ and $`\mathrm{\Gamma }=+11`$.


format: text

text:
# References

hep-th/0204143

Field Theory of Tachyon Matter

Ashoke Sen

Harish-Chandra Research Institute

Chhatnag Road, Jhusi, Allahabad 211019, INDIA

and

Department of Physics, Penn State University

University Park, PA 16802, USA

E-mail: asen@thwgs.cern.ch, sen@mri.ernet.in

Abstract

We propose a field theory for describing the tachyon on a brane-antibrane system near the minimum of the potential. This field theory realizes two known properties of the tachyon effective action: 1) absence of plane-wave solutions around the minimum, and 2) exponential fall off of the pressure at late time as the tachyon field evolves from any spatially homogeneous initial configuration towards the minimum of the potential. Classical solutions in this field theory include non-relativistic matter with arbitrary spatial distribution of energy.

Dynamics of the tachyon on D-brane anti-D-brane system or unstable D-branes in superstring theory and bosonic string theory have been investigated during the last few years. Of the various known properties of the classical tachyon effective action,111Although we shall use the language of tachyon effective action obtained by integrating out the other open string modes, we should note that many of these results are derived using the technology of full string field theory, containing all the infinite number of fields describing the dynamics of the open string. two specific properties which specify the behaviour of the tachyon effective action around the minimum of the potential are as follows:

1. Since the minimum of the potential describes a configuration where there are no D-branes, around this minimum there are no physical open string excitations. Translated to a property of the tree level tachyon effective action, this implies that there are no plane wave solutions to the linearized equations of motion around the minimum of the tachyon potential.
2. If we let the tachyon roll beginning with any spatially homogeneous initial configuration, it evolves asymptotically towards its minimum instead of oscillating about the minimum. The total energy density is conserved during the evolution but the pressure evolves to zero. The information about the tachyon effective action near the minimum of the potential is encoded in the late time behaviour of the pressure. This is given by222Since we are interested in studying tree level properties of the tachyon potential, we shall ignore the backreaction of gravity and other closed string fields on the evolution of the tachyon. Once the form of the action is derived by this method, one can of course couple it to background gravitational field and study the evolution of the coupled system. These effects have been studied in refs.. For some earlier attempts at the study of time dependent solutions for the tachyon field see .

$$p=Ke^{\alpha x^0},$$
(1)
where $`K`$ is an irrelevant normalization constant which can be changed by shifting the origin of the time coordinate $`x^0`$, and $`\alpha `$ is given (in the $`\alpha ^{}=1`$ unit) by

$`\alpha `$ $`=`$ $`1\text{for bosonic string theory}`$ (2)
$`=`$ $`\sqrt{2}\text{for superstring theory}.`$

Several field theory models have been proposed for realising the first property. These models are of two kind.333Here we are stating the results in terms of properly redefined field $`\varphi `$ for which the two derivative term has the standard form $`\frac{1}{2}\eta ^{\mu \nu }_\mu \varphi _\nu \varphi `$. In the first kind the tachyon potential is non-singular at the minimum of the potential but the absence of physical plane-wave solutions has its origin in the structure of the kinetic term (which typically involves higher derivative terms). This possibility seems to be realized in the explicit analysis of the string field theory based on $``$-product interaction. In the second kind of field theory models the tachyon potential is singular at the minimum. In particular the second derivative of the potential blows up and hence the tachyon mass becomes infinite. Analysis based on boundary string field theory seems to support this view. Of course these results are not mutually incompatible, since an effective field theory of the first kind may be turned into one of second kind and vice versa by a suitable field redefinition that involves (infinite number of) derivatives of the tachyon field.

A field theory model that realises the evolution of the system to a pressureless gas was proposed in  following earlier proposals in . The action was taken to be

$$S=d^{p+1}xV(T)\sqrt{detA},$$
(3)
where

$$A_{\mu \nu }=\eta _{\mu \nu }+_\mu T_\nu T.$$
(4)
However in this proposal the tachyon potential $`V(T)`$ was left undetermined. The only property of $`V(T)`$ that was determined is that the minimum of $`V(T)`$ where $`V(T)=0`$ should be at $`T=\mathrm{}`$.

In this paper we shall show that by taking $`V(T)e^{\alpha T/2}`$ for large $`T`$ in the action (3), we can satisfy both the requirements given above: the absence of plane-wave solutions around the tachyon vacuum, and the exponential fall off of the pressure for large $`x^0`$. The constant of proportionality can be absorbed into a constant shift of $`T`$. Thus the proposed tachyon effective action for large $`T`$ on an unstable D-$`p`$-brane system is

$$S=d^{p+1}xe^{\alpha T/2}\sqrt{detA}=d^{p+1}xe^{\alpha T/2}\sqrt{1+\eta ^{\mu \nu }_\mu T_\nu T},$$
(5)
where $`A_{\mu \nu }`$ is given by eq.(4). We expect this to be a valid approximation when the second and higher order derivatives of $`T`$ are small. As we shall see, the tachyon matter does satisfy this condition.

On a brane-antibrane system where $`T`$ is a complex field, the tachyon potential is given by $`e^{\alpha |T|/2}`$. In this case the action (5) can be regarded as the restriction of the full action to the case $`\mathrm{}(T)=0`$. We shall carry out our analysis for this restricted field configuration. The imaginary part of the tachyon, which can be regarded as the Goldstone mode associated to the broken $`U(1)`$ phase symmetry, is absorbed by the U(1) gauge field due to Higgs mechanism, and is expected to decouple from the dynamics near the minimum of the potential.

We shall first verify that the action (5) produces the correct large $`x^0`$ behaviour of the pressure. For this analysis we can restrict to spatially homogeneous, time dependent field configurations. As discussed in , for such configurations the conserved energy density is given by:

$$T_{00}=e^{\alpha T/2}(1(_0T)^2)^{1/2}.$$
(6)
Since $`T_{00}`$ is conserved, we see that for any given $`T_{00}`$, as $`T\mathrm{}`$, $`_0T1`$. In particular, for large $`x^0`$ the solution has the form:444Note that since for this configuration $`_0^nT`$ for $`n>1`$ falls off exponentially, the contribution from possible higher derivative corrections like $`(_\mu ^\mu T)^2`$ inside the square root on the right hand side of (5) will be suppressed.

$$T=x^0+Ce^{\alpha x^0}+𝒪(e^{2\alpha x^0}).$$
(7)
In order to see that (7) gives the correct form of the solution we simply need to note that the leading contribution to $`T_{00}`$ computed from this configuration remains constant in time:

$$T_{00}\frac{1}{\sqrt{2\alpha C}}.$$
(8)
The pressure associated with this configuration is given by:

$$p=e^{\alpha T/2}(1(_0T)^2)^{1/2}\sqrt{2\alpha C}e^{\alpha x^0}.$$
(9)
This is in precise agreement with (1).

Next we shall demonstrate the absence of plane-wave solutions. First let us expand the action in powers of derivatives of $`T`$ up to terms containing at most two powers of derivatives. This gives

$$S=d^{p+1}xe^{\alpha T/2}(1+\frac{1}{2}\eta ^{\mu \nu }_\mu T_\nu T+\mathrm{}),$$
(10)
where $`\mathrm{}`$ denotes terms with higher powers of derivatives. Defining

$$\varphi =e^{\alpha T/4},$$
(11)
we get

$$S=d^{p+1}x\frac{16}{\alpha ^2}(\frac{1}{2}\eta ^{\mu \nu }_\mu \varphi _\nu \varphi \frac{\alpha ^2}{16}\varphi ^2+\mathrm{}).$$
(12)
The minimum of the potential is at $`\varphi =0`$. If we ignore the higher derivative terms then the quantization of the action seems to lead to a scalar particle of mass $`\alpha /2\sqrt{2}`$. This is reflected in the fact that the equations of motion without the higher derivative terms contain a plane-wave solution of the form:

$$\varphi =ae^{ik_\mu x^\mu },$$
(13)
with $`\eta ^{\mu \nu }k_\mu k_\nu =\alpha ^2/8`$ and $`a`$ an arbitrary constant.

We shall now show that when we take into account the effect of the higher derivative terms, (13) ceases to be a solution of the equations of motion for any $`k_\mu `$.555Note that (13) describes a complex solution. Since the full equations of motion are homogeneous of degree one, but not linear, even if we had a solution of the form (13), it is not clear how we could construct a real solution since we cannot superpose the solutions. We shall show however that even this complex solution does not exist. To show this we note that the full action (5) written in terms of $`\varphi `$ defined in (11) takes the form:

$$S=d^{p+1}x\varphi ^2\sqrt{1+\frac{16}{\alpha ^2}\varphi ^2\eta ^{\mu \nu }_\mu \varphi _\nu \varphi }.$$
(14)
The action is homogeneous of degree 2 in $`\varphi `$. The equations of motion derived from this action is given by:

$$\eta ^{\mu \nu }_\mu \left(\frac{_\nu \varphi }{\sqrt{1+\frac{16}{\alpha ^2}\varphi ^2\eta ^{\mu \nu }_\mu \varphi _\nu \varphi }}\right)+\frac{\alpha ^2}{8}\frac{\varphi +\frac{8}{\alpha ^2}\varphi ^1\eta ^{\mu \nu }_\mu \varphi _\nu \varphi }{\sqrt{1+\frac{16}{\alpha ^2}\varphi ^2\eta ^{\mu \nu }_\mu \varphi _\nu \varphi }}=0.$$
(15)
Substituting (13) into (15) we get:

$$\frac{a}{\sqrt{1\frac{16}{\alpha ^2}\eta ^{\mu \nu }k_\mu k_\nu }}=0.$$
(16)
Clearly this equation has no non-trivial solution for finite values of $`k_\mu `$. This establishes the absence of plane-wave solutions.666Note again that since for the configuration (13) $`_\mu _\nu T`$ vanishes, contribution to the equations of motion from possible higher derivative corrections like $`(_\mu ^\mu T)^2`$ inside the square root on the right hand size of (5) vanishes. (Note that if we expand the left hand side of this equation in powers of $`k_\mu `$ and keep up to quadratic terms, we reproduce the mass-shell condition $`\eta ^{\mu \nu }k_\mu k_\nu =\alpha ^2/8`$ derived earlier.)

Absence of plane wave solutions does not imply absence of other classical solutions however. We have already seen the existence of solutions with constant (but arbitrary) energy density. As emphasized in  (see also ) the correct excitations of the system can be found by working in the hamiltonian formalism. Defining the momentum conjugate to $`T`$ as:

$$\mathrm{\Pi }(x)=\frac{\delta S}{\delta (_0T(x))},$$
(17)
we can construct the Hamiltonian $`H`$ following :

$$H=d^px,=T_{00}=\sqrt{\mathrm{\Pi }^2+e^{\alpha T}}\sqrt{1+_iT_iT}.$$
(18)
For large $`T`$ we can ignore the $`e^{\alpha T}`$ term, and the equations of motion take the form:

$$_0\mathrm{\Pi }(x)=\frac{\delta H}{\delta T(x)}=_j\left(\mathrm{\Pi }(x)\frac{_jT}{\sqrt{1+_iT_iT}}\right).$$
(19)
$$_0T(x)=\frac{\delta H}{\delta \mathrm{\Pi }(x)}=\sqrt{1+_iT_iT},$$
(20)
for $`\mathrm{\Pi }(x)>0`$. Thus we can get a solution to the equations of motion by taking $`_iT=_0\mathrm{\Pi }=0`$, $`_0T=1`$. This gives

$$\mathrm{\Pi }(x)=f(\stackrel{}{x}),T(x)=x^0,$$
(21)
where $`f(\stackrel{}{x})`$ is any arbitrary function of the spatial coordinates. The energy density associated with such a solution is proportinal to $`f(\stackrel{}{x})`$. This shows that the system admits classical solutions with energy density which is time independent but has arbitrary dependence on the spatial coordinates. Classically this energy density can be as low as we like. The ability to create configurations with arbitrarily low energy density can be traced to the scale invariance:

$$T(x)\lambda ^1T(\lambda x),\mathrm{\Pi }(x)\lambda ^{p+1}\mathrm{\Pi }(\lambda x),$$
(22)
under which the Poisson brackets and the equations of motion remain unchanged, and $`(x)\lambda ^{p+1}(\lambda x)`$.

Since the original theory is Lorentz invariant, it is clear that given a localized density of tachyon matter we should be able to boost it to an arbitrary velocity. Intuitively we also expect that is should be possible to construct configurations with different local velocities at different points. This can be seen to be the case by rewriting the equations of motion in a slightly different form. We define:

$$u_\mu =_\mu T,ϵ(x)=\mathrm{\Pi }(x)/_0T(x).$$
(23)
The equations of motion then take the suggestive form:

$$\eta ^{\mu \nu }u_\mu u_\nu =1,_\mu (ϵ(x)u^\mu )=0.$$
(24)
The energy momentum tensor can be computed following . Expressed in terms of these new variables, $`T_{\mu \nu }`$ take the form:

$$T_{\mu \nu }=ϵ(x)u_\mu u_\nu .$$
(25)
These are precisely the equations govering gradient flow of non-interacting dust, with $`u_\mu `$ interpreted as the local velocity vector. These equations are expected to be valid as long as $`_\mu u_\nu `$ is small in magnitude (in string units).

Although in our analysis we have ignored the coupling of the tachyon to various massless fields on the D-brane world-volume, these can be easily incorporated following . We can generalize (5) to curved background space-time by replacing $`\eta _{\mu \nu }`$ in eq.(4) by the closed string metric $`g_{\mu \nu }`$. If we also have anti-symmetric tensor field background then $`\eta _{\mu \nu }`$ should be replaced by $`g_{\mu \nu }+B_{\mu \nu }`$. The dilaton $`\varphi `$ couples via an overall factor of $`e^\varphi `$ multiplying the lagrangian density. For non-BPS D-$`p`$-brane of type II string theories, another quantity of interest is the coupling of the tachyon to Ramond-Ramond (RR) field $`C^{(p)}`$. In particular it is known that the D-$`p`$-brane world-volume theory has a coupling:

$$d^{p+1}xf(T)𝑑TC^{(p)},$$
(26)
where $`f(T)`$ is some function of the tachyon field $`T`$. We can determine the behaviour of $`f(T)`$ for large $`T`$ if we can find the source of the RR field that is generated by the rolling tachyon field. This can be computed from the boundary state associated with the solution. Although we have not performed a systematic analysis of the problem along the lines of , general arguments based on symmetry and other considerations lead us to guess that (in $`\alpha ^{}=1`$ unit) the source for the RR $`p`$-form field is proportional to:

$$\mathrm{sin}(\stackrel{~}{\lambda }\pi )\left[\frac{e^{x^0/\sqrt{2}}}{1+\mathrm{sin}^2(\stackrel{~}{\lambda }\pi )e^{\sqrt{2}x^0}}\frac{e^{x^0/\sqrt{2}}}{1+\mathrm{sin}^2(\stackrel{~}{\lambda }\pi )e^{\sqrt{2}x^0}}\right],$$
(27)
where $`\stackrel{~}{\lambda }`$ is the parameter labelling the total energy density of the system:

$$T_{00}=\frac{𝒯_p}{2}(1+\mathrm{cos}(2\stackrel{~}{\lambda }\pi )),$$
(28)
$`𝒯_p`$ being the tension of the non-BPS D-$`p`$-brane. Thus for large $`x^0`$ the source is proportional to $`e^{x^0/\sqrt{2}}`$. Since $`Tx^0`$ for large $`x^0`$, this indicates that for large $`T`$, the coupling of the RR $`p`$-form field to the tachyon has the form:

$$d^{p+1}xe^{T/\sqrt{2}}𝑑TC^{(p)},$$
(29)
up to an overall normalization constant.

There are many issues which call for further investigation. As we have seen, the effective field theory describing tachyon matter has classical solutions with localized energy densities. It will be interesting to construct the two dimensional conformal field theories associated with these solutions along the lines of refs. where the conformal field theories associated with spatially homogeneous energy densities were constructed. (Of course, one possibility is to take an arbitrary spatial distribution of coincident D0-$`\overline{\mathrm{D}}0`$ pair and construct the rolling tachyon solution on each such pair.) Another issue of importance is the effect of quantum corrections on such background.

Acknowledgement: I would like to thank L. Kofman, J. Maldacena and T. Padmanabhan for useful discussions. This work was supported in part by a grant from the Eberly College of Science of the Penn State University.


text:
# Coronal heating distribution due to low-frequency wave-driven turbulence

## 1 Introduction

The mechanism by which the solar corona is heated is a fundamental problem in astrophysics and plasma physics, and one that remains incompletely explained. An acceptable mechanism transports energy to the corona, heating it to temperatures greatly in excess of chromospheric values. This mechanism accelerates the solar wind and establishes the boundary conditions for the entire plasma heliosphere. There are several constraints on an acceptable coronal heating model. Spectrometers (Kohl et al., 1995, 1997) show that proton temperatures in polar coronal holes rise to several $`10^6`$ K at heliocentric distances of 1.5 to 2 solar radii ($`R_s`$) and the bulk wind speed reaches 200 km/s. Remote sensing observations (Grall et al., 1996) indicate wind speeds up to 600 km/sec at 4 $`R_s`$. Fluid models of the solar wind (Habbal et al., 1995; McKenzie et al., 1995; Evje & Leer, 1998) support the conclusion that wind acceleration at such low altitudes requires heating low in the corona. In the fluid models an *ad hoc* heating function $`Q(r)`$ is assumed, which represents the energy deposition per unit volume as $`QQ_0\mathrm{exp}\left[\frac{rr_0}{L}\right]`$, where $`r`$ is heliocentric distance, $`r_0=R_s7\times 10^5`$ km, $`L`$ is a prescribed dissipation lengthscale, and $`Q_0`$ is chosen so that the energy flux has the required value of about $`5\times 10^5\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ (Withbroe & Noyes, 1977). Despite the fact that there is no direct observational evidence for the exponential heat function, it has been long known to provide (Holzer & Axford, 1970; Kopp & Orrall, 1976; Hammer, 1982) rapid acceleration of the wind, and with $`L0.2`$$`0.5`$$`R_s`$ it accounts for many observed features of the corona (McKenzie et al., 1995; Habbal et al., 1995; Axford & McKenzie, 1997). Although some possible sources of heating are proposed, in those models there is no physical mechanism for explaining this heating function and, being an ad hoc quantity, the exact explicit exponential function expression is perhaps irrelevant and the main point to be remarked is the requirement of strong heating per unit volume in the first solar radius of the corona. Yet other two-fluid (electrons and protons) models (Tu & Marsch, 1997; Marsch & Tu, 1997) have included some mechanisms for the heating, based on the requirement of high frequency waves launched from the Sun, but those high frequencies (of the order of KHz) waves have yet remained unobserved in the coronal base where the energy is injected. In another previous approach, Hollweg, 1986; Hollweg & Johnson, 1988; Isenberg, 1990 conjectured that waves are launched from the Sun at low frequencies but magnetohydrodynamic (MHD) turbulent-like cascades are assumed to transfer the energy to higher frequencies, where ion cyclotron mechanisms can be invoked to dissipate it. This idea, recently employed by Li et al., 1999 and reviewed by Cranmer, 2000 involve phenomenological turbulent rates (being called of the Kolmogorov and Kraichnan type) but two important problems are faced. First, a cascade in wave frequencies is assumed, which implies a turbulent cascade proceeding in the parallel (to the mean field) wavenumber direction and a linear dispersion relation to persist between parallel wavenumber and frequency. MHD simulation studies show on the other hand that in the presence of a strong mean magnetic field, like in the corona, the parallel wavenumber cascade is almost completely supressed (Shebalin et al., 1983; Oughton et al., 1994; Matthaeus et al., 1998) in favor of the perpendicular cascade. This is precisely the kind of anisotropy which is well documented (Armstrong et al., 1990) for the lower coronal density.

Second, the necessity of counterpropagating fluctuations to sustain the turbulence (Kraichnan, 1965; Dobrowolny et al., 1980) is not considered in the models and a WKB approach is implicitly assumed. In those models, the turbulent cascade was not the central focus of study, but we believe that a useful perspective can be gained into the heating problem by further study of this issue. In this present work, instead of assuming a parallel wave cascade, we describe *ab initio* the conditions under which a strong (perpendicular) turbulent cascade can actually be developed and sustained by low frequency driving at the base. Reflections from density and mean field gradients provide the necessary counterpropagating fluctuations to sustain the turbulence and this effect is fully taken into account in our model. We will employ both direct numerical simulation and turbulent phenomenologies to obtain the heating rate distribution.

We show that the heating distribution profile is related to the coronal background magnetic field and density profile, i.e., the heating and the density profile are not independent quantities. On the one hand, we find that the heating per unit volume is significant mostly within the first solar radius, with a dissipation length scale (defined as the length in which the heating decays by an order of magnitude) of a fraction of solar radius. On the other hand the heating per unit mass varies much less dramatically and is significant through the complete domain. For a complete model, appropriate kinetic processes will need to be identified to absorb the small perpendicular scale cascaded energy. Candidates for small scale damping are plentiful, and include oblique kinetic Alfvén waves, collisionless reconnection and associated heating at current sheets and nonlinear beam instabilities (Leamon et al., 1998). The kinetic physics of heating is beyond the scope of the present paper. The main purpose of this work is to show that a heating distribution emerges naturally from a consistent MHD turbulent cascade driven by low frequency fluctuations at the coronal base. This result may bring us closer to understand the puzzle of coronal heating and the origin of the solar wind.

## 2 Model

Most coronal heating models derive their energy supply from either wave power or quasistatic field line motions, originating in the photosphere or chromosphere. These differ mainly in their characteristic timescales. A general scenario for wave production in open magnetic regions has been proposed (McKenzie et al., 1995) in which activity in the network regions, driven by convective motions, provides a source of upward traveling Alfvén waves. Higher frequency waves are easily transported upward, whereas low-frequency waves experience greater reflection and transmit less readily into the corona (Parker, 1965; Heinemann & Olbert, 1980; An et al., 1990; Zhou & Matthaeus, 1990; Velli, 1993). However even in the limit of sharp gradients or discontinuities, there are indications (Hollweg, 1981) that a substantial energy flux can reach the coronal base (Kudoh & Shibata, 1999).

Here we examine the following scenario (Matthaeus et al., 1999): MHD-scale ( low-frequency ) fluctuations are excited below the base of the corona, on transverse length scales that are characteristic of the chromospheric network, and propagate upwards along the mean magnetic field $`𝐁`$. A portion of this wave flux enters the corona. Reflections from large-scale inhomogeneities within the lower corona produce counter-propagating fluctuations (Heinemann & Olbert, 1980; Zhou & Matthaeus, 1990) a situation that permits strong nonlinear MHD couplings (Dobrowolny et al., 1980). The couplings drive a perpendicular (to $`𝐁`$) turbulent cascade, producing small scale reconnection events that couple to kinetic processes at small scales and oblique wavevectors (Leamon et al., 1998).

In our model, we simulate the propagation of low-frequency MHD fluctuations in a medium whose properties are chosen to represent an open magnetic field region (like a polar coronal hole) in the lower corona. We consider a volume section which extends in the radial direction a distance $`L_s`$ of order $`2R_s`$ above the surface. A sketch of such a region is shown in Fig. 1.

Our goal is to compute the dynamics of transverse fluctuations of magnetic field $`𝐛`$ and velocity $`𝐯`$ as influenced by specified large-scale inhomogeneities (mean magnetic field and density gradients) as well as by local nonlinear couplings. To do this, in a medium where there is a strong mean magnetic field, we use the reduced MHD (RMHD) approximation (Strauss, 1976; Montgomery, 1982; Zank & Matthaeus, 1992). This approach is appropriate for fluctuations in a nearly incompressible plasma, where the plasma beta (ratio of thermal to magnetic pressure) is low (Zank & Matthaeus, 1992). We employ the Elsasser variables $`𝐳_\pm `$ = $`𝐯`$ $`\pm `$ $`𝐛`$ (where the magnetic field is written in units of $`𝐛/\sqrt{4\pi \rho }`$, with $`\rho `$=density) to represent downward (+) and upward (-) propagating fluctuations. Fluctuations depend on the radial coordinate $`r`$, the transverse spatial coordinates (including an areal expansion factor $`A(r)`$) and on time $`t`$. Mean quantities like density $`\rho `$ and magnetic field $`𝐁`$ are assumed to depend only on $`r`$. Consistency with the RMHD assumptions (Zank & Matthaeus, 1992) requires $`_r_{}`$, $`V_Az_\pm `$, and $`_{}𝐳_\pm =0`$, where $`𝐕_𝐀=𝐁/\sqrt{4\pi \rho }`$ is the Alfven speed which gives the propagation velocity of fluctuations in the linear limit.

In the presence of a non-uniform $`V_A`$, the nonlinear RMHD equations can be written (Zhou & Matthaeus, 1990) as

$`{\displaystyle \frac{𝐳_{}}{t}}+V_A{\displaystyle \frac{𝐳_{}}{r}}`$ $`=`$ $`R_1𝐳_++R_2𝐳_{}𝐳_+_{}𝐳_{}+\eta _{}^2𝐳_{}`$
$`{\displaystyle \frac{𝐳_+}{t}}V_A{\displaystyle \frac{𝐳_+}{r}}`$ $`=`$ $`R_1𝐳_{}R_2𝐳_+𝐳_{}_{}𝐳_++\eta _{}^2𝐳_+`$ (1)
where $`\eta `$ is the resistivity (assumed equal to the viscosity), and $`R_1(r),R_2(r)`$ are specified reflection rates due to the inhomogeneous $`V_A(r)`$. For a symmetric radial geometry ($`A(r)=r^2`$), $`B=B_0(r_0/r)^2`$ (with $`B_0`$ the magnetic field at $`r_0=R_s`$), the reflection rates are $`R_1(r)=\frac{1}{2}\frac{dV_A}{dr}`$ and $`R_2(r)=\frac{1}{2}\frac{dV_A}{dr}+\frac{V_A}{r}`$. In the case of a general geometry with cross section area $`A(r)`$, where $`A(r)=A_0B_0/B(r)`$, the reflection rates are $`R_1(r)=\frac{1}{2}\frac{dV_A}{dr}`$ and $`R_2(r)=\frac{1}{2}\frac{dV_A}{dr}+\frac{1}{2}\frac{V_A}{A}\frac{dA}{dr}`$. The previous expressions are valid for a section with transverse dimensions small compared to the longitudinal size and close to the polar region in the Sun (dependence of the mean field $`B`$ with the latitudinal angle should be assumed to maintain the free divergence condition, $`𝐁=0`$).

A Chebyshev-Fourier representation is chosen to solve the equations (1) using a direct numerical simulation of the pseudospectral type. Transverse coordinates are periodic, while non periodic boundary conditions are imposed on the bottom and top sections of the volume. In order to model the input of velocity fluctuations at the base of the coronal domain, we refer to observations of non-thermal velocity amplitudes in the upper transition region of the Sun. These amplitudes are reported to range from 20 to 55 km/sec (Chae et al., 1998; Doyle et al., 1998; Hassler et al., 1990). The system of nonlinear equations (1) is driven by imposing a time dependent oscillatory pattern of velocity and magnetic field at the bottom surface, with a transverse lengthscale corresponding to supergranule diameters of 13-32 Mm at the coronal base (Hagenaar et al., 1997) and the typical inter-network spacing of 30 Mm (Axford & McKenzie, 1997). Transverse dimensions of the simulation box are chosen to include about 4 supergranules. In the bottom panel of Fig. 2 the boundary conditions imposed on the form of waves at the bottom of the domain are shown. The top panel of Fig. 2 shows a cross section of the volume above in the coronal simulation (corresponding to one of the simulations to be described below). Numerical resolution of $`128\times 128\times 33`$ Fourier+Chebyshev grid points are employed in this simulation, for a macroscopic Reynolds number of order $`600`$. The fluctuating fields are highly structured on the transverse directions as shown by concentrations of current density in the form of sheet-like small scale structures.

Driving is at a low frequency of $`f=0.6/t_A`$, where $`t_A=R_s/V_{A_0}`$ is a vertical Alfvén crossing time based on the Alfvén velocity at the base and a distance of a solar radius. This correspond to driving wave periods of the order of 1000 sec.

Open boundary conditions are assumed at the top surface. If reflections were absent, this would lead to a wavetrain of upward propagating fluctuations that would escape continuously from the upper boundary without depositing any energy in the system. The system would dynamically relax toward a maximal cross helicity state (Dobrowolny et al., 1980; Grappin et al., 1983; Ting et al., 1986) in which downward propagating fluctuations vanish and there would be no incompressible MHD nonlinearities and hence no cascade. However, the first two terms on the right hand sides of Eqs. (1) are associated with gradients of the wave speed, and these lead to partial reflection of the excited waves. This results in a state with both upward and downward type fluctuations and thus, sustained turbulence. We have presented in a previous numerical study (Dmitruk et al. 2001a, ) the conditions under which an efficient MHD turbulent perpendicular cascade can be sustained; what is required is occurrence of reflections and non-propagating structures (controlled by the type of boundary conditions imposed). The efficiency of this mechanism on dissipating the injected energy has been demonstrated by both phenomenologies (Matthaeus et al., 1999; Dmitruk et al. 2001b, ) as well as direct numerical simulations (Dmitruk et al. 2001a, ; Oughton et al., 2001). Efficiencies (ratio of turbulent dissipation rate by the energy injection rate) averaging between 10 and 40 % can be obtained. Even very weak reflections are capable of sustaining the turbulence and provide effective dissipation.

In the next section we obtain the heating distribution profiles resulting from this turbulent state.

## 3 Heating profiles

We performed numerical simulations employing several different mean magnetic field and density radial profiles. The values at the bottom boundary (Axford & McKenzie, 1997; Feldman et al., 1997) have been fixed to $`B_0=9`$ Gauss, and a number density of $`n_0=3.2\times 10^8\mathrm{cm}^3`$, which gives an Alfven speed at the base of $`1100`$ km/s. With the assumed amplitude of fluctuations of the order of $`2055`$ km/s, this gives an average energy input flux per unit area, $`F_{A_0}\rho _0V_{A_0}<z_{}^2z_+^2>`$ of order $`5\times 10^5\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ corresponding to the Withbroe & Noyes, 1977 value. These values have to be taken as order of magnitude quantities but a different choice would not invalidate the general argument we want to present here. In our first simulation we consider a symmetric radial mean magnetic field profile $`B(r)=B_0(R_s/r)^2`$. The mean density follows from the hydrostatic equilibrium approximation (An et al., 1990; Velli, 1993; Nakariakov et al., 2000): $`dp/dr`$ = -$`m_ing`$, with $`p`$ = thermal pressure, $`n`$ = number density, $`m_i`$ the proton mass and $`g`$ = gravity. In an isothermal corona, this equation takes the form $`dn/dr=\alpha _TnR_s/r^2`$ where $`\alpha _T=\mu GM/R_gTR_s`$ ($`\mu `$ mean molecular weight, $`G`$ gravitational constant, $`M`$ solar mass, $`R_g`$ gas constant and $`T`$ temperature). The quantity $`\alpha _T12`$ for the Sun and a fully ionized hydrogen coronal plasma at $`T=10^6`$ K. The number density distribution is then $`n(r)`$ = $`n_0`$ exp\[$`\alpha _T(R_s/r1)`$\] and the density is $`\rho (r)=m_in(r)`$. We consider two such density profiles, shown in the top left panel of Fig. 3 with continuous and dotted lines, corresponding to values of $`\alpha _T=12`$ and $`\alpha _T=6`$ (representing a hotter corona). This gives two different density profiles, with the hotter corona having the shallower profile. The corresponding Alfven speed profiles $`V_A(r)=B(r)/\sqrt{4\pi \rho (r)}`$ are shown on the top right panel of Fig. 3. When the system reaches a statistically steady state the turbulent dissipation rate $`ϵ`$ (equal to the cascade rate) is computed from the fluctuating fields $`𝐳_{}`$, $`𝐳_+`$ from which the heating per unit volume is obtained as $`Q=\rho ϵ`$.

This quantity is plotted as a function of $`r`$ in the bottom left panel of Fig. 3. It is apparent that the dissipation per unit volume is peaked near the model coronal base. It is important to note that $`Q(r)`$ decreases faster for the model with the steeper density profile. The radial distance $`\mathrm{\Delta }r`$ at which $`Q`$ decreases by an order of magnitude with respect to the coronal base value $`Q_0`$ is $`0.36R_s`$ and $`1.1R_s`$ for the two profiles in Fig. 3. This would be the analog of the dissipation length in the ad-hoc exponential heating distribution profiles. The bottom right panel on Fig. 3 shows the heating per unit mass, a useful quantity to consider the heating in a per particle basis. Over our entire domain this quantity maintains a relatively constant value. Extended heating per unit mass fits the apparent requirement of a hot (non adiabatic) solar wind and has also emerged in empirical models of radiation loss and other heating models (Mullan and Cheng, 1994; Gibson, 1973).

A third simulation has been performed employing a super radial geometry, where the mean magnetic field profile has been taken as

$$B=1.5\left[(f_{max}1)\left(\frac{R_s}{r}\right)^{3.5}+\left(\frac{R_s}{r}\right)^2\right]\mathrm{Gauss}$$
(2)
Profiles of this type have been considered in Hollweg, 2000 and similar super expansion profiles have been used by McKenzie et al., 1995; Evje & Leer, 1998; Li et al., 1999 in their solar wind models. We have taken $`f_{max}=6`$ which gives a boundary value of $`B_0=9`$ Gauss. An observationally based number density is assumed in this case (Feldman et al., 1997),

$$n(r)=3.2\times 10^8\left(\frac{R_s}{r}\right)^{15.6}+2.5\times 10^6\left(\frac{R_s}{r}\right)^{3.76}\mathrm{cm}^3$$
(3)
The density profile and resulting Alfven speed profile are shown in the top panels of Fig. 3 (dashed line). The obtained heating per unit volume and heating per unit mass are shown in the bottom panels of Fig. 3 (dashed line). The result confirms the tendency shown for the earlier cases: this profile, with the fastest decreasing density, also has the fastest decreasing heating per unit volume. The decay length for this case is $`\mathrm{\Delta }r=0.2R_s`$. The heating per unit mass remains at a relatively constant value.

The relevant results to be extracted from these simulations for different Alfven speed profiles is that the heating per unit volume has significant values only within the first solar radius of the lower corona and that the radial profile of the heating is controlled by the radial profile of the density. As it will be more clear from the phenomenological analysis that follows in the next section, it is actually through the Alfven speed profile and its gradient that this behavior is manifested.

Whether the radial profile of the heating per unit volume is exactly exponential or not is, we believe, irrelevant and we conjecture that a steep heating profile, like the ones we obtained here, would work as well to explain the fast acceleration of the solar wind.

Fig. 3 contains the essence of the present results, that a reasonable and efficient heating can be obtained from the MHD perpendicular turbulent cascade driven by low frequency waves injected from the bottom. The properties of this heating are controlled by both the amount of energy injected into the system and the background profiles of magnetic field and mean density. To further illustrate this idea and, in some limiting case, to obtain an analytical expression for the heating profile in the low region of the volume considered here, we present in the following section a phenomenological approach to the previous model Eqs (1).

## 4 Heating profile from a phenomenological model

A nonlinear phenomenology (Dmitruk et al. 2001b, ) simplifies the study of the wave-driven RMHD heating model. We modify Eqs. (1) as follows: the independent variables, now designated as $`Z_\pm (r)`$, are treated as one dimensional wave amplitudes with the same linear transport terms (LHS and reflection terms) as in Eqs. (1), but with the nonlinear and dissipative terms on the RHS replaced by the nonlinear models $`Z_{}|Z_\pm |/2\lambda _{}(r)`$. This model entirely suppresses the transverse structure, with the strength of the nonlinear cascade effects represented through a single correlation lengthscale $`\lambda _{}(r)`$. Previous investigations (Dmitruk et al. 2001b, ) have shown that this phenomenology portrays many of the same physical effects as wave-driven models based upon the full RMHD equations (Dmitruk et al. 2001a, ; Oughton et al., 2001). We expect this model to correspond to high Reynolds number fully developed RMHD turbulence since the amplitudes $`Z_{}`$ and $`Z_+`$ are strongly coupled in the adopted forms for the nonlinear terms. We extracted a heating function $`Q(r)`$ from a numerical solution to this phenomenological model, using the same parameters as those employed in the full nonlinear RMHD simulations and assuming a correlation lengthscale $`\lambda _{}(r)`$ linearly increasing with radius, with a value at the coronal base set to correspond to the inter-network spacing of 30 Mm. The derived $`Q(r)`$ is illustrated in Fig. 4 as the thin line. The direct numerical solution of the previous section is shown as the thick line. This corresponds to the simpler hydrostatic density and $`1/r^2`$ mean magnetic field profile (similar results are obtained for the super expansion geometry). Despite some differences on the far region, where $`dV_A/dr`$ approaches $`0`$ and the heating is very low, the agreement is pretty good, considering the extreme simplification involved in the phenomenological treatment. Of course, the phenomenology would not be useful for studies of the turbulent cascades and spectra, but it gives the right order of magnitude for the averaged quantities such as the dissipation (heating), when this quantity is not extremely low. The interesting issue of how the correlation length and the strength of the cascade depend on the external parameters warrants further study comparing the phenomenology with the full non linear direct numerical simulation.

In Fig. 5 we explore the phenomenological heat function $`Q`$ obtained by varying the correlation lengthscale at the coronal base, with $`\lambda _{}=1`$ corresponding to 30 Mm. The limit $`\lambda _{}0`$ corresponds to stronger local turbulence since the implied eddy turnover time $`\lambda _{}/z_{}`$ becomes much smaller than the wave transit time $`R_s/V_A`$. Fig. 5 shows that $`Q`$ increases as $`\lambda _{}0`$, while maintaining the fast radially decreasing behavior.

In this limit, an analytical asymptotic expression can be obtained. Although limited in its validity, this analysis serves to illustrate better the properties of the heating in this model. We begin by formulating a flux balance equation, based upon Eq.(1). This takes the form $`dF/dr=A(r)Q`$, where the net upward energy flux $`F=\rho AV_A(Z_{}^2Z_+^2)`$ and $`A`$ is the cross section area. Note that in geometrical optics (or “WKB” theory for noninteracting waves) (Jacques, 1977) the flux balance formally requires that $`Q0`$. When fluctuations are damped, $`Q0`$ (Hollweg, 1986). Here the crucial relation that connects the heating to the Alfvén speed profile emerges by identifying $`Q`$ with a physically correct turbulent dissipation rate per unit volume at distance $`r`$, namely $`Q=\rho ϵ`$, where the dissipation rate per unit mass $`ϵ`$ is expressed in the phenomenological model as $`ϵ(Z_+^2|Z_{}|+Z_{}^2|Z_+|)/\lambda _{}`$. We now proceed to form an asymptotic strong turbulence limit as follows. Letting $`\lambda _{}0`$, one can conclude (Dmitruk et al. 2001b, ) that the downward flux is small, $`|Z_+||Z_{}|`$, and that in a steady state this limit is characterized by $`|Z_+|/\lambda _{}=|dV_A/dr|`$. Thus $`ϵZ_{}^2|Z_+|/\lambda _{}`$ and we obtain the asymptotic result that the dissipation per unit length can be related to the Alfvén speed profile and the energy flux: $`AQ\rho AZ_{}^2|dV_A/dr|V_A^1|dV_A/dr|F`$. Thus, the asymptotic flux balance equation can now be written as $`dF/drV_A^1|dV_A/dr|F`$, whose solution is

$$F\{\begin{array}{cc}F_0V_{A_0}/V_A\hfill & \text{if }r<r_m\hfill \\ F_0(V_{A_0}/V_{A_m})^2V_A/V_{A_0}\hfill & \text{if }r>r_m\hfill \end{array}$$
(4)
where $`F_0`$ is the net flux at the base and $`r_m`$ is the radial distance at which $`V_A`$ reaches its maximum value $`V_{A_m}`$.

This allows us finally to obtain an explicit expression for the heating deposition

$$Q(r)\{\begin{array}{cc}F_{A_0}(A_0/A)|dV_A/dr|(V_{A_0}/V_A^2)\hfill & \text{if }r<r_m\hfill \\ F_{A_0}(A_0/A)|dV_A/dr|(V_{A_0}/V_{A_m}^2)\hfill & \text{if }r>r_m\hfill \end{array}$$
(5)
where $`F_{A_0}=F_0/A_0`$ is the energy flux per unit area at the base (i.e, the Withbroe & Noyes, 1977 value).

For the particular case of the $`A=r^2`$ geometry and the isothermal atmosphere used before, with $`r_m3R_s`$, this relation predicts the top dot dashed curve of Fig. 5. As seen in this figure, the phenomenological solutions for different $`\lambda _{}`$ values approach this asymptotic case as $`\lambda _{}`$ tends to 0. Already for $`\lambda _{}=1/10`$ the asymptotic solution would be a good approximation. This corresponds to correlation length scales of order 3000 km. Although much less than the typical inter network distances, observations do not rule out correlation scales of this magnitude. Eq. (5) also shows that in this limiting case the value of the heating goes to $`0`$ at points where $`dV_A/dr=0`$, a limitation which is only apparent at corresponding very low values of $`Q`$ in the numerical solution (a similar result applies for the super expansion profile presented before, where the point at which $`dV_A/dr=0`$ is $`r_m1.6R_s`$).

For this asymptotic limit, eq. (5) explicitly shows that the spatial distribution of turbulent heating is directly connected to the radial profile of the Alfvén speed and its absolute value depends on the energy input flux. This illustrates the connection already pointed out through the direct numerical simulations results.

## 5 Conclusions

We have presented a model for MHD turbulence driven by low-frequency waves imposed at the lower boundary. We applied this model to an open region in the lower corona, adopting reasonable profiles for magnetic field and density. Both direct numerical simulations and phenomenological studies have been performed and compared. Results show that turbulence can be sustained, due to reflection of upward traveling waves by Alfvén speed gradients, producing a perpendicular turbulent cascade and efficient heating. The energy dissipation per unit volume is significant mostly in the lower part of the model corona with a dissipation length scale of a fraction of solar radius. The obtained heating profile $`Q`$ is similar to profiles commonly used in solar wind acceleration models. However the more general result emerges, for strong turbulence, that the heat function $`Q(r)`$ should be related to the Alfvén speed profile $`V_A(r)`$. This shows that the heating distribution is determined by the large-scale magnetic field and density profiles. In this way one sees that the scale-height of the heating per unit volume is directly related to the density scale-height, thus explaining the spatial confinement of $`Q`$ as a consequence of the rapid density decrease in the lower corona. However, the deposition of heat per unit mass $`Q/\rho `$ is much more extended, perhaps well beyond the lower coronal region considered in the present model. These conclusions do not rule out other models, such as direct cyclotron absorption of high-frequency waves (Axford & McKenzie, 1997; Tu & Marsch, 1997) but they imply that it is possible to provide the proper kind of heating per unit volume using a low-frequency wave-driven MHD turbulence mechanism.

Research supported by NASA (NAG5-7164, NAG5-8134 SEC Theory Program), NSF (ATM-9713595, ATM-0296113), and UK PPARC (PPA/G/S/1999/00059).

Fig. 1 – Sketch of the coronal region considered in this model

Fig. 2 – In the bottom panel the wave forcing applied at the bottom boundary is illustrated with the imposed transverse fluctuating field $`𝐳_{}`$ (arrows) overlaying the current density (grayscale). In the top panel a cross section at a fraction $`0.2R_s`$ above the base is shown with transverse fluctuating field (arrows) over the current density (grayscale) obtained in the numerical simulation. Note the small scale structure formation.

Fig. 3 – The panel at the top left shows three different number density profiles assumed for the simulations. The continuous and dotted line correspond to the hydrostatic model for different values of the density scale height. The dashed line correspond to an observationally based density (see text). The top right panel is the corresponding Alfven speed profile, with $`1/r^2`$ geometry for the hydrostatic density and super expansion geometry for the observational density. The bottom left panel shows the result of the heating per unit volume from the simulations. The bottom right panel is the heating per unit mass.

Fig. 4 – Comparison between the heating per unit volume obtained with the direct numerical simulation of Eqs.(1) (thick line) and the numerical solution of the phenomenological model (thin line), for the case of an hydrostatic density profile and $`1/r^2`$ geometry.

Fig. 5 – Different solutions for the heating per unit volume of the phenomenological model (thin lines) for values of the correlation length $`\lambda _{}=1,0.5,0.1,0.02`$ (in units of the inter network length) which approach the asymptotic solution when $`\lambda _{}0`$ (Eq.(5) in the text) shown as the dot dashed line.


text:
# An alternative search for the electron capture of 123Te

## I Introduction

In the past decades the investigation of $`\beta `$-decay has played a major role in understanding weak interactions and their structure. Even though being a little bit out of fashion, there are still interesting problems to be investigated. Among them is the second forbidden unique electron capture (EC) of 123Te occuring to the ground state of 123Sb with a transition energy of $`53.3\pm 0.2`$ keV a . Two positive evidences obtained for this decay mode show a discrepancy by six orders of magnitude, namely $`(1.24\pm 0.10)10^{13}`$ yrs b  and $`(2.4\pm 0.9)10^{19}`$ yrs c . Even with the lower value already in slight disagreement with other bounds on the lifetime of K shell capture in 123Te d  it is quoted quite often in literature e . However, still only a lower bound of about $`10^{13}`$ yrs is used toi . One of the main problems of the past measurements was associated with the fact that observations relied on the outside detection of K X-rays following the deexcitation of 123Sb. This line is at 26.1 keV, very close to the Te K X-ray line at 27.3 keV which might be excited by other processes like cosmic rays and radioactive background.
In the measurement reported here, a new technique for searches of rare nuclear decays using CdZnTe semiconductor detectors has been used for the first time. These offers several advantages. The natural abundance of 123Te is 0.9 %, hence it is intrinsic to any detector containing Te. Due to this fact, the total cascade after K-capture is contained within the detector resulting in a peak at 30.5 keV. Cd(Zn)Te detectors are well known devices in X-ray and $`\gamma `$-ray astronomy and offer good energy resolution. The study was performed within a pilot project for the COBRA-experiment COBRA , which plans to use a large amount of CdZnTe detectors to explore double beta decay.

## II Experimental setup

The measurements were performed with a $`10\times 10\times 5`$ mm3 CdZnTe detector provided by eV-Products111eV-PRODUCTS, Saxonburg, PA, USA, www.evproducts.com. The detector is encapsulated within an aluminium-tube. It is working with the coplanar grid technology, focussing on the readout of electrons CPG . This avoids asymmetric peaks normally seen because of the insufficient collection of holes due to trapping within CdZnTe.

The energy resolution was measured with various sources, among them 241Am and 57Co, resulting in lines at 59.5, 122.1 and 136.5 keV respectively. An FWHM-energy resolution of $`\mathrm{\Delta }E=4.7`$ keV is determined for 30.5 keV by extrapolation.
To be sensitive for rare decays, a special shielding had to be constructed to reduce background events due to environmental radioactivity and cosmic rays. A schematic plot of the setup is shown in Fig.1. In this way, a $`50\times 50\times 50`$ cm box made of copper and lead was designed. The inner $`20\times 20\times 20`$ cm consisted of electrolytic copper, the outer part was made of spectroscopy lead. In a first step, the copper was cleaned in an ultrasonic bath. All lead bricks were etched shortly before installation with HCl and HNO3 to clean the surfaces, the copper was electropolished. The complete shielding was surrounded by an air-tight box of aluminium, which was flushed with nitrogen. This prevents air and with it 222Rn from penetrating into the apparatus. To suppress signals due to cosmic ray muons, an active veto consisting of 19 plastic scintillators was used. They covered all four sides of the box in a single layer with about 95 % detection efficiency as well as the top side in two layers resulting in about 99.5 % efficiency. All scintillators were fed into discriminators and a common OR was used to form the veto signal. In that way after the firing of any of the scintillators the data aquisition was blocked for 20 $`\mu `$s resulting in a dead-time of about $`1.6`$ %. Only accepted signals were fed into a multichannel analyser card in a PC, serving as DAQ-system. The shielding depth was about 5 meter water equivalent (mwe).

## III Results

Single runs were taken with 30 minutes life-time each. 392 runs were obtained of which 2 had to be removed resulting in 195 hours of good data. The obtained spectrum of interest is shown in Fig.2. No signal is visible at 30.5 keV. Fitting the region from $`17.5`$ to $`45.5`$ keV with a $`1/E+const`$ function results in a $`\chi ^2/\text{dof}`$ of $`11.7/11`$ within the $`2\mathrm{\Delta }`$E-region around the peak position, a fit including a Gaussian peak at the right position and with the correct energy resolution gives a $`\chi ^2/\text{dof}`$ of $`11.8/11`$ while the gaussian content ist consistent with zero. Excluding the $`2\mathrm{\Delta }`$E-region around $`30.5`$ keV (i.e. 25 to 36 keV) from the fitpoints, the resulting $`\chi ^2/\text{dof}`$ of the pure background fit is still $`12.0/11`$. Therefore, the background-only hypothesis is well within allowed parameters.

Using the fit without the peak region and adding a Gaussian more than 31.2 events can be ruled out (95 % CL). Knowing the mass of $`2.89`$ g the upper bound can be converted into a lower limit on the half-life for the electron capture of 123Te :

$$T_{1/2}>3.210^{16}\text{yrs}95\%\text{CL}$$
(1)
This clearly rules out possible observations with half-lives of the order of $`10^{13}`$ yrs and is in strong favor of the results of Alessandrello et al. c . Longer half-lives are also supported from theoretical considerations f ; g .

## IV Summary

A novel technique for studying rare nuclear decays using CdZnTe semiconductor detectors was used to investigate the second forbidden unique electron capture of 123Te . It takes advantage of the fact that source and detector are identical. The main motivation for this search was to discriminate between two published results which differ by six orders of magnitude. No evidence for this process was found and a half-life limit larger than $`3.210^{16}`$ yrs (95 % CL) is obtained. This clearly rules out the solutions obtaining $`10^{13}`$ yrs. The measurement was performed during tests for the planned COBRA-experiment to search for double beta decay. A proof for the quoted half-life of $`2.4\pm 0.910^{19}`$ yrs can only be done with a setup using a larger amount of CdZnTe detectors.

## V Acknowledgements

We thank H. Kiel and C. Gößling for valuable discussions and support. We also acknowledge the support of the mechanical workshop of the University of Dortmund during the construction of the test-setup. Finally we want to thank eV-Products for generously lending the detector that was used for our measurements.