text:
# Effects of small surface tension in Hele-Shaw multifinger dynamics: an analytical and numerical study

## I Introduction

The displacement of a viscous fluid by a less-viscous one in a Hele-Shaw cell, the so-called Saffman-Taylor problem Saffman58 ; Bensimon86a ; McCloud95 ; Tanveer00 ; Casademunt00 , is a prototypical pattern formation problem. Since the seminal work of Saffman and Taylor Saffman58  a considerable effort has been aimed at understanding both steady and unsteady interfacial patterns formed during this flow. The Saffman-Taylor problem is the simplest member of a wide class of interfacial pattern formation problems such as free dendritic growth, directional solidification, or chemical electro-deposition Langer87 ; Pelce88 ; Kessler88 . As such, a theoretical understanding of Hele-Shaw flow may help elucidate generic behavior common to many pattern forming systems. Despite its relatively simple formulation and the large amount of work devoted to it, however, several aspects of interfacial dynamics in Hele-Shaw flow are still poorly understood, in particular concerning the highly nonlinear and nonlocal dynamics of finger competition.

One of the main reasons for the wide interest in Hele-Shaw flow, at least from a mathematical point of view, is that explicit time-dependent solutions can be found in the case of zero surface tension Howison86 ; Dawson94 ; Entov95 ; Magdaleno00 . However, it is also known that the zero surface tension Saffman-Taylor (ST) problem is ill-posed as an initial value problem Tanveer93  and finite-time singularities appear frequently Howison85 . Nevertheless, rather large classes of zero surface tension solutions have been found which exhibit the variety of morphologies observed both in experiments and numerical simulations. Then, the question that naturally arises is to what extent smooth (nonsingular) zero surface tension solutions reproduce the dynamics of the physical problem with finite surface tension, in particular in the limit of vanishing dimensionless surface tension, $`B0`$.

It is well known that surface tension is a singular perturbation to the zero surface tension problem Tanveer93 . This singular character manifests dramatically in the classical selection problem posed by Saffman and Taylor Saffman58  and only solved three decades later Shraiman86 ; Hong86 ; Combescot86 ; Tanveer87 , where an arbitrarily small surface tension selects out a single, stable solution from the continuum of steady single-finger $`B=0`$ solutions. Another manifestation of the singular nature of surface tension which is directly relevant to the present work is its effect on the dynamics. Siegel, Tanveer and Dai Siegel96a ; Siegel96b  showed that interfacial evolution for the regularized problem (i.e., vanishingly small $`B`$) may differ signifficantly from that for the $`B=0`$ problem in order one time. Then, smooth time dependent solutions of the $`B=0`$ case do not coincide, in general, with the limiting solutions for $`B0`$. Accordingly, the study of finite $`B`$ dynamics encounters considerable difficulties since $`B=0`$ solutions cannot be naively used as a starting point for the study of the problem with finite $`B`$.

The physical content of exact zero surface tension solutions with pole-like singularities has been recently addressed in Refs. Casademunt00 ; Magdaleno98 ; Paune02a  using a dynamical systems approach. Through a detailed study it has been shown that the exact zero-surface tension phase flow, considered in a global sense, is structurally unstable. In other words, the zero surface tension phase dynamics are not topologically equivalent to the phase space flow of the physical problem, regularized by surface tension. Indeed, the zero surface tension phase flow omits the necessary saddle-point structure of multifinger fixed points, which is crucial to the physical finger competition process Paune02a . A natural extension of the well known solvability mechanism (first applied to ‘select’ a finger of width 1/2 out of a continuum of solutions in the single finger case) was proposed for multi-finger solutions in Paune02a ; this helps clarify how the introduction of surface tension modifies the global phase space structure of the flow and restores the hyperbolicity of multi-finger fixed points.

The approach of Ref. Paune02a , however, was qualitative in nature and could not quantify the extent to which zero surface tension trajectories might resemble the evolution with small surface tension. In particular it was recognized that, while some trajectories appear to be qualitatively correct for infinite time, others may have a dramatically different evolution. In particular, adding an infinitessimal surface tension could give the opposite outcome in the finger competition, that is, make the ‘losing’ finger for $`B=0`$ become the ‘winning’ finger when $`B>0`$, for sufficiently generic sets of initial conditions.

A satisfactory analytical understanding of the problem with $`B0`$ has been achieved in two regimes: the initial linear instability of the flat interface followed by the weakly non-linear regime Alvarez01 , and the asymptotic regime where surface tension selects the width of the single finger Shraiman86 ; Hong86 ; Combescot86 ; Tanveer87 . The highly non-linear intermediate regime that connects the quasi-planar interface with the asymptotic single-finger regime has mostly been studied through numerical computation. The first exhaustive numerical studies were reported by Tryggvason and Aref Tryggvason83 ; Tryggvason85 , who paid considerable attention to the influence of viscosity contrast on the problem and studied both single-finger and multi-finger configurations. Later, Casademunt and Jasnow Casademunt91 ; Casademunt94  showed that the basin of attraction of the single-finger solution depends strongly on viscosity contrast and that only when one of the two fluid viscosities is negligible it can be claimed that the single finger is the universal attractor of the problem. In the present work we will restrict to this limiting case. DeGregoria and Schwartz DeGregoria85 ; DeGregoria86  observed that well-developed fingers split when surface tension is sufficiently decreased. This tip-splitting instability is related to the fact that the Saffman-Taylor finger is linearly stable but non-linearly unstable, and the size of the perturbation that triggers the tip-splitting decreases quickly with surface tension Bensimon86b . Dai and Shelley Dai94  showed that for small $`B`$ numerical computations are extremely sensitive to the precision used in the computations. As a consequence noise level has to be controlled with care in order to ensure that the computation is sufficiently accurate. Hou et al. Hou94  developed a numerical method that deals with the numerical stiffness of the problem in an efficient manner, aiding the ability to perform long time computations. More recently Ceniceros et al. Ceniceros99 ; Ceniceros00 , using very high precision arithmetic have been able to study the effect of extremely small surface tension in the circular geometry with suction, and they have observed that surface tension can produce complex ramified patters even without the presence of noise.

An analytical treatment of this highly nonlinear and nonlocal free-boundary problem faces challenging difficulties. In particular, a perturbative study for small $`B`$ is complicated by the ill-posedness of the zero surface tension problem. Tanveer Tanveer93  was able to overcome this obstacle by embedding the zero surface tension problem in a well-posed one. In addition, this well-posed extension of the $`B=0`$ problem allowed Baker et al. Baker95  to develop a numerical method to compute the time evolution of zero surface tension dynamics in a well-posed manner. Once the $`B=0`$ problem is formulated in a well-posed way the $`B0`$ case can be studied using a perturbative approach. The main result of the asymptotic perturbative theory developed by Tanveer Tanveer93  is that the effect of surface tension may be manifest in a $`O(1)`$ time: the evolution of the same initial interface for $`B=0`$ and $`B0`$ will in general differ after a time of order one, even if the $`B=0`$ solution is smooth for all time. Siegel et al. Siegel96b  have extended the work of Ref. Tanveer93  to later stages of the evolution, and through numerical computation for very small values of $`B`$ they showed that smooth $`B=0`$ solutions are indeed significantly affected by the presence of arbitrarily small $`B`$ in order-one time, thus confirming the predictions of the perturbative theory. The zero surface tension solutions studied by Siegel et al. Siegel96a ; Siegel96b  in the channel geometry were single-finger solutions with an asymptotic width $`\lambda `$, specifically chosen to be incompatible with selection theory for vanishing surface tension. They found that the singular effect of surface tension was to widen the finger in order to reach the selected width. The surprising feature here is that the effect of surface tension is felt in order-one time, i.e., that the time lapse for which the regularized solution approaches the unperturbed one as $`B0`$ is bounded.

The present paper expands the work of Refs. Siegel96a ; Siegel96b  in the spirit of Ref. Paune02a , towards the study of multi-finger solutions. However, unlike the studies of Siegel96a ; Siegel96b  we chose zero surface tension multi-finger soltuions which are compatible with asymptotic selection theory, that is, with an asymptotic finger width $`\lambda =1/2`$—the selected value in the limit $`B0`$. In this way we isolate the intrinsic finger competition dynamics from the selection effects responsible of restoring the asymptotic width. Two different kinds of two-finger zero surface tension solutions are studied, and in both cases it is shown that surface tension acts as a singular perturbation to the dynamics in order-one time, modifying dramatically the late time configuration of the interface not only quantitatively but also qualitatively. Specifically we show that paths in phase space associated with zero and nonzero surface tension evolution, and indeed the global topological structure of the the phase spaces, may differ appreciably, even for arbitrarily small $`B`$. In physical terms, our evidence suggests that the presence of arbitrarily small surface tension can completely alter outcome of finger competition when compared with the zero surface tension evolution.

The paper is organized as follows: in Sec. II the equations describing Hele-Shaw flow are introduced, and a class of two-finger zero surface tension solutions relevant to two-finger competition is presented and briefly discussed. This class of solutions will be used as initial condition for numerical computation with $`B>0`$. In Sec. III the basic features of the asymptotic theory are recalled, and the theory is applied to the zero surface tension solutions introduced in the previous section. The numerical computations with finite (but small) $`B`$ are presented in Sec. IV. Sec. V discusses and summarizes the results obtained in previous sections.

## II Zero surface tension

In this section we present the equations which govern the interfacial dynamics in a rectilinear Hele-Shaw cell, following the formalism of Tanveer93 . We consider a class of exact, time-dependent zero surface tension solutions that are relevant to the finger competition problem, and briefly describe the solutions within this class.

Consider Hele-Shaw flow in the channel geometry, in which a fluid of negligible viscosity displaces a viscous liquid. The equations governing the interfacial evolution can be conveniently formulated by first introducing a conformal map $`z(\zeta ,t)`$ which takes the interior of the unit semicircle in the $`\zeta `$ plane into the region occupied by the viscous fluid in the complex plane $`z=x+\mathrm{i}y`$, in such a way that the arc $`\zeta =e^{\mathrm{i}s}`$ for $`s[0,\pi ]`$ is mapped to the interface and the diameter of the semi-circle is mapped to the channel walls111For interfaces symmetric with respect to the central axis of the channel our formulation also describes periodic boundary conditions, that is, an infinite array of fingers.. The mapping function $`z(\zeta ,t)`$ has the form $`z(\zeta ,t)=(2/\pi )\mathrm{ln}\zeta +\mathrm{i}+f(\zeta ,t)`$, and inside and on the unit semi-circle we require $`f(\zeta ,t)`$ to be analytic and $`z_\zeta (\zeta ,t)0`$. In addition, we require that

$$\mathrm{Im}f=0$$
(1)
on the real diameter of the semi-circle. This latter condition ensures that $`z`$ maps the diameter to the channel walls. Under suitable assumptions (see Tanveer93 ) the Schwartz reflection principle may be applied to show that $`f`$ is analytic and $`z_\zeta 0`$ for $`|\zeta |1`$.

The effective velocity field, averaged across the plate gap, is a two-dimensional potential flow satisfying Darcy’s law

$`𝐮=\phi .`$ (2)
Here $`\phi `$ is a velocity potential defined by $`\phi =(b^2/12\mu )p`$, where $`p`$ is the pressure, $`\mu `$ is the viscosity and $`b`$ is the gap width. Under the assumption of incompressibility $`(𝐮=0)`$ the potential satisfies Laplace’s equation

$`^2\phi =0.`$ (3)
Incompressibility also implies the existence of a stream function $`\psi `$. Therefore, one can define a complex velocity potential $`W(z,t)=\phi +\mathrm{i}\psi `$ which is analytic for $`z`$ in the fluid region of the channel. Its form as a function of $`\zeta `$ reads

$`W(\zeta ,t)=(2/\pi )\mathrm{ln}\zeta +\mathrm{i}+\omega (\zeta ,t)`$ (4)
where $`\omega (\zeta ,t)`$ is an analytic function inside the unit circle. The condition that there is no flow through the walls implies that $`\mathrm{Im}\omega =0`$ must hold on the real diameter of the unit semi-circle. In the absence of surface tension, $`\omega =0`$ (see Eq. (6)).

At the interface we impose the usual boundary conditions. The kinematic boundary condition states that the normal component of fluid velocity at a point on the interface equals the normal velocity of the interface at that point, and takes the form

$$\mathrm{Re}\left[\frac{z_t}{\zeta z_\zeta }\right]=\frac{1}{|z_\zeta |^2}\mathrm{Re}\left[\zeta W_\zeta \right].$$
(5)
The dynamic boundary condition specifies that the pressure jump across the interface is balanced by surface tension, and is given by

$$\mathrm{Re}\omega =\frac{B}{|z_\zeta |}\mathrm{Re}\left[1+\zeta \frac{z_{\zeta \zeta }}{z_\zeta }\right].$$
(6)
The parameter $`B`$ is the nondimensional surface tension and is defined by $`B=b^2T/12\mu Va^2`$, where $`T`$ is the surface tension, $`V`$ is the fluid velocity at infinity and $`a`$ is half the cell width. The equations given in (1)-(5) are in nondimensional form, with lengths and velocities nondimensionalized by $`a`$ and $`V`$, respectively.

When $`B=0`$ it is well known that pole singularities in $`z_\zeta `$ (i.e. in $`f_\zeta `$) present in the exterior of the unit disk are preserved under the dynamics, i.e., such singularities are neither created nor destroyed, although the location of those which are initially present will evolve with time. Exact $`B=0`$ solutions consisting of a collection of pole singularities with constant amplitude have been the focus of extensive studies (see e.g. Howison86  ). The simplest such solution leading to nontrivial finger competition consists of a pair of singularites in the upper halfplane of $`|\zeta |>1`$, located at positions that are symmetric with respect to the y-axis. A second pair of poles conjugate to the first pair is required to satisfy the symmetry restriction (1). This exact solution takes the form Casademunt00 ; Magdaleno98 ; Paune02a

$`z(\zeta ,t)={\displaystyle \frac{2}{\pi }}\mathrm{ln}\zeta +{\displaystyle \frac{1}{\pi }}(1\lambda +iϵ)\mathrm{ln}(1{\displaystyle \frac{\zeta ^2}{\zeta _s(t)^2}})`$
$`+{\displaystyle \frac{1}{\pi }}(1\lambda iϵ)\mathrm{ln}(1{\displaystyle \frac{\zeta ^2}{\overline{\zeta }_s(t)^2}})+d(t)+i`$ (7)
where $`\lambda `$ and $`ϵ`$ are real constants with $`0<\lambda <1`$ and $`ϵ0`$, and $`d(t)`$ is real. The singularity locations are given by the complex parameter $`\zeta _s(t)`$, which satisfies a simple differential equation given in Magdaleno98 . Analyticity of $`f(\zeta ,t)`$ in the unit circle implies that $`|\zeta _s(t)|>1`$. We employ the convention that $`\zeta _s(t)`$ is a complex number in the first quadrant. The amplitudes of the singularities, given here by the numbers $`1\lambda +iϵ`$ and its conjugate, are chosen so that the asymptotic form of the solution consists of one or two steadily propagating fingers of total width $`\lambda `$. The parameter $`ϵ`$ determines the nature of the finger competition for $`B=0`$.

The solution (II) describes generically two different fingers of the nonviscous phase penetrating the viscous one. In the linear regime $`|\zeta _s(t)|>>1`$ the interface consists of a single bump or finger, and as time increases a second finger may develop and grow, depending on the value of $`\mathrm{Arg}\zeta _s(0)`$.

We summarize the features of the solution (II) that are most relevant to the study of finger competition. Consider first $`ϵ=0`$. In this case the asymptotic configuration consists of one or two fingers of total width $`\lambda `$, depending on the initial condition. The singularities move toward the unit disk, with the limit as $`t\mathrm{}`$ denoted by $`\zeta _s(t)e^{i\theta }`$. When $`\theta =0`$ the asymptotic configuration is a single Saffman-Taylor finger growing in the center of the channel (this asymptotic configuration is denoted ST(R)), for $`\theta =\pi /2`$ it is a ‘side’ Saffman-Taylor finger i.e. a pair of half fingers of total width $`\lambda `$ with tips located at the cell walls (denoted ST(L)), and for $`\theta =\pi /4`$ it is a ‘double’ Saffman-Taylor finger, namely two identical fingers of width $`\lambda /2`$ with tips at $`x=0,\pm 1`$ (denoted 2ST). For any other value of $`\theta `$ the asymptotic configuration consists of two unequal steadily growing fingers. The two-finger asymptotic configuration is a consequence of the continuum of fixed points that is present in the phase portrait of the dynamical variables which govern the shape of the interface, namely $`(\mathrm{Re}\zeta _s(t),\mathrm{Im}\zeta _s(t))`$. In order to correspond to the notation of Magdaleno98  introduce the variable $`\alpha (t)=\alpha ^{}(t)+\mathrm{i}\alpha ^{\prime \prime }(t)=1/(i\zeta _s^2(t))`$. Then the planar interface corresponds to $`\alpha =0`$, the center Saffman-Taylor finger to $`\alpha =i`$, the side Saffman-Taylor finger to $`\alpha =i`$, and the double Saffman-Taylor finger to $`\alpha =1`$. Figure 1 shows the phase portrait of the dynamical system obtained from the substitution of the mapping defined by Eq. (II) into the evolution Eqs. (5,6) for $`B=0`$, using the dynamical variables $`(\alpha ^{},\alpha ^{\prime \prime })`$ . The asymptotic states with two advancing fingers corresponds in the dynamical system $`(\alpha ^{},\alpha ^{\prime \prime })`$ to a continuum of fixed points given by $`|\alpha |=1`$. Therefore, for $`ϵ=0`$ the solution (II) does not exhibit finger competition. In addition, it is important to note that the evolution of (II) with $`ϵ=0`$ is free of finite time singularities, i.e., $`z_\zeta 0`$ in the domain $`|\zeta |1`$ for all time.

For $`ϵ0`$ the continuum of fixed points is removed (see Fig. 1), as is the double Saffman-Taylor finger fixed point. Consequently, the solution to Eq. (II) exhibits ‘successful’ competition, in the sense that the asymptotic interface shape consists of a single Saffman-Taylor finger or side Saffman-Taylor finger. The price to pay is the appearance of finite time singularities for a certain subset of initial conditions, in the form of a zero of $`z_\zeta `$ impacting the unit disk (this is a generic feature of conformal map solutions $`z_\zeta `$ composed of a finite number of pole singularities—see Howison86 ). Then, only the subset of initial conditions free of finite time singularities is capable of sustaining finger competition all the way to the $`t\mathrm{}`$ outcome. Nevertheless, one may ask whether the class of $`B=0`$ solutions that are free of finite time singularities may describe, at least qualitatively, the physical finger competition for positive surface tension in the limit $`B0`$.

In the following sections we focus on the class of initial data for which the $`B=0`$ solutions are devoid of finite time singularities, and examine the $`B>0`$ dynamics. We develop a general theory for how the presence of positive surface tension affects the outcome of finger competition. This will enable us to predict the winner of finger competition, i.e., the eventual asymptotic state. Most interestingly, we find instances in which the presence of arbitrarily small surface tension leads to dramatically different outcomes in finger competition when compared with the zero surface tension evolution.

## III Asymptotic theory

Little is known about the effect of finite (but small) surface tension $`B`$ on the dynamics of zero surface tension multifinger solutions, and in particular on the class of exact solutions (II). For single finger configurations, steady state selection theory predicts that the finger cannot have an arbitrary width. Indeed, for vanishing surface tension $`B0`$ the width $`\lambda =1/2`$ is selected, asymptotically in time. Thus, it is clear that surface tension has a critical influence on single finger solutions with $`\lambda 1/2`$. The nature of this influence in the limit $`B0`$ has been investigated by Siegel, Tanveer and Dai Siegel96a ; Siegel96b , who present evidence that zero surface tension single finger solutions with $`\lambda <1/2`$ are significantly perturbed by the inclusion of an arbitrarily small amount of surface tension in order one time. The effect of surface tension is an increase of the finger width to reach the width predicted by selection theory.

Consider now the effect of small surface tension on the exact ($`B=0`$) two finger solution (II). When $`0<B1`$ the asymptotic perturbation theory developed in Refs. Siegel96a ; Tanveer93 ; Siegel96b  can be applied. This perturbation theory describes the effects of the introduction of a small amount of surface tension on initial data $`z(\zeta ,0)`$ specified in the extended complex plane, i.e., in a domain including the ‘unphysical’ region $`|\zeta |>1`$ (the extended domain is required to make the $`B=0`$ problem well-posed). The effect of finite $`B`$ is most important near isolated zeros and singularities of $`z_\zeta (\zeta ,0)`$, where a regular perturbation expansion in $`B`$ breaks down. (Away from these points the perturbation expansion is regular, at least initially.) For the class of solutions (II) we are discussing, the isolated singularities of $`z_\zeta (\zeta ,0)`$ are simple poles. The theory suggests that the introduction of finite surface tension modifies the poles ($`\zeta _s`$) by transforming them into localized clusters of $`4/3`$ singularities, but these clusters move at leading order according to the $`B=0`$ dynamics. Thus the effect of one of these clusters on the interface is approximately equivalent to that of the unperturbed ($`B=0`$) pole-like singularity that has given birth to it.

The influence of surface tension on the zeros of $`z_\zeta (\zeta ,0)`$ is more complex. Each initial zero instantly gives birth to two localized inner regions, i.e., regions where the $`B=0`$ and $`B>0`$ solutions differ by $`O(1)`$. One of the two inner regions moves, at least initially, according to the $`B=0`$ dynamics of the original zero $`\zeta _0`$ 222This region can move differently from the $`B=0`$ zero once the second inner region discussed below has impacted the unit disk. We do not consider this possibility here.. Since the particular zero surface tension solutions considered here have zeros that are either bounded away from the unit disk for all time or impact the unit disk only after long times, the inner region around $`\zeta _0(t)`$ has a negligible influence on the interface. The second inner region created around $`\zeta _0(0)`$ moves differently. The theory suggests that this inner region consists of a cluster of singularities, whose size scales like $`B^{1/3}`$. Unlike the case discussed above this second inner region moves away from the $`B=0`$ zero since, to leading order in $`B`$, it moves like a singularity of the zero surface tension problem and this speed is different form the speed of the zero $`\zeta _0(t)`$ which spawned the cluster. As this singularity cluster approaches the physical domain it may perturb the flow and the interface shape may differ significantly from that at $`B=0`$ shape. The location of this singularity cluster will be denoted by $`\zeta _d(t)`$, and following Tanveer93  we shall call it the daughter singularity. We emphasize that the dynamics of the daughter singularity cluster is determined at lowest order solely by the $`B=0`$ solution $`z^0(\zeta ,t)`$, at least until it arrives at the surroundings of the unit circle, and therefore can be simply computed once the initial locations of the zeros of $`z_\zeta (\zeta ,0)`$ are determined.

The daughter singularity evolution equation is given by (see Tanveer93 )

$$\dot{\zeta }_d(t)=q_1^0(\zeta _d(t),t);\zeta _d(0)=\zeta _0(0)$$
(8)
where $`q_1^0`$ is defined by

$$q_1^0=\frac{\zeta }{2\pi i}_{|\zeta ^{}|=1}\frac{d\zeta ^{}}{\zeta ^{}}\frac{\zeta +\zeta ^{}}{\zeta ^{}\zeta }\frac{\mathrm{Re}\left[\zeta ^{}W_\zeta ^0(\zeta ^{},t)\right]}{|z_\zeta ^0(\zeta ^{},t)|^2}$$
(9)
and the superscript 0 denotes that the function evaluations are done using the corresponding $`B=0`$ solution. The function $`q_1^0(\zeta ,t)`$ also gives the characteristic velocity of a pole or branch point singularity of $`z_\zeta (\zeta ,t)`$ located at position $`\zeta `$ in the region $`|\zeta |>1`$. The initial position $`\zeta _d(0)`$ is a consequence of the fact that each zero $`\zeta _0(0)`$ of the zero surface tension solution will give birth to a daughter singularity. From Eq. (8) it can be shown Tanveer93  that $`d|\zeta _d|/dt<0`$, so that the daughter singularity approaches the unit circle and it can impact it in a finite time $`t_d`$, the daughter singularity impact time, satisfyng $`|\zeta _d(t_d)|=1`$. In the limit $`B0`$, the daughter singularity impact time $`t_d`$ signals the time when the effects of the surface tension are felt on the physical interface. For times larger than $`t_d`$ the $`B=0`$ interface and the $`B0`$ are expected to differ significantly.

For the family of exact $`B=0`$ solutions the mapping function (II) has four pole-like singularities: $`\pm \zeta _s`$ and $`\pm \overline{\zeta }_s`$, and four zeros $`\pm \zeta _{0+}`$ and $`\pm \zeta _0`$ of $`z_\zeta `$ located at

$`\zeta _{0+}^2={\displaystyle \frac{(\lambda +iϵ)\zeta _s^2(\lambda iϵ)\overline{\zeta }_s^2}{2(12\lambda )}}+{\displaystyle \frac{\sqrt{\left[(\lambda +iϵ)\zeta _s^2+(\lambda iϵ)\overline{\zeta }_s^2\right]^2+4(12\lambda )|\zeta _s|^4}}{2(12\lambda )}}`$ (10a)
$`\zeta _0^2={\displaystyle \frac{(\lambda +iϵ)\zeta _s^2(\lambda iϵ)\overline{\zeta }_s^2}{2(12\lambda )}}{\displaystyle \frac{\sqrt{\left[(\lambda +iϵ)\zeta _s^2+(\lambda iϵ)\overline{\zeta }_s^2\right]^2+4(12\lambda )|\zeta _s|^4}}{2(12\lambda )}}.`$ (10b)
For the particular case $`\lambda =1/2`$ this solution presents only one pair of zeros $`\pm \zeta _0`$ located at

$$\zeta _0^2=\frac{|\zeta _s|^4}{2[(\lambda +iϵ)\zeta _s^2+(\lambda iϵ)\overline{\zeta }_s^2]}.$$
(11)
In the following it will be useful to define the real quantity $`\beta =(\lambda +iϵ)\zeta _s^2(\lambda iϵ)\overline{\zeta }_s^2`$ which appears in (10) and (11).

Depending on the value of $`\lambda `$ the initial data may have zeros on both the real and imaginary axes, or all the zeros may lie on a single axis. This difference has significant consequences in the finite surface tension dynamics. More specifically, when $`\lambda <1/2`$ the zeros described in (10a) and (10b) are located on both the real and imaginary axes of $`|\zeta |>1`$, namely at $`\pm |\zeta _{0+}|`$ and $`\pm i|\zeta _0|`$. The situation is different for $`\lambda >1/2`$, which is further divided into two cases, depending on whether $`\beta ^2+4(12\lambda )|\zeta _s|^2>0`$ or $`<0`$. In the former case all four singularities lie on the real axis (for $`\beta >0`$) or on the imaginary axis (for $`\beta <0`$). In the latter case the four zeros are located off the axes in conjugate pairs, i.e. at $`\pm \zeta _0`$ and $`\pm \overline{\zeta }_0`$. Finally, when $`\lambda =1/2`$ the solution (II) has only two zeros, located on the real axis at $`\pm |\zeta _s|^2/\sqrt{2\beta }`$ when $`\beta <0`$ and on the imaginary axis at $`\pm |\zeta _s|^2/\sqrt{2\beta }`$ when $`\beta >0`$. Note that for $`\lambda =1/2`$ the $`B=0`$ solution has two less zeros than for $`\lambda 1/2`$.

The initial zero locations described above have a critical bearing on whether the daughter singularity will impact the unit disk 333Although the daughter singularity is said to impact the unit disk when $`|\zeta _d(t)|=1`$, the singularities comprising the cluster do not actually impinge upon the unit disk. However, at $`t_d`$ the singularities are close enough to influence the interface shape, in the sense that $`|z_\zeta (\zeta _d,t_d)z_\zeta ^0(\zeta _d,t_d)|=O(1)`$.. Although all daughter singularities approach the unit disk, their impact may be shielded by the presence of an inner region corresponding to a pole singularity. More precisely, since $`\zeta _d`$ and $`\zeta _s`$ obey the same dynamical equation, they will move together if they get close enough to each other. However, the inner region around a pole moves to leading order like the $`B=0`$ pole, i.e., it moves exponentially slowly toward $`|\zeta |=1`$ when $`|\zeta _s|1<<1`$, and does not impinge upon the unit disk in finite time Howison86 . In this case the $`O(B^{1/3})`$ inner region around the daughter singularity will not affect the dynamics on $`|\zeta |=1`$, at least until $`t=O(\mathrm{ln}B)`$. Before this time, we expect the interface to be uninfluenced by the presence of the daughter singularity. This shielding mechanism is discussed in the context of single fingers in Siegel96b .

Knowledge of the $`t\mathrm{}`$ asymptotic state and the initial locations of zeros can be used to ascertain whether shielding can occur. The $`B=0`$ asymptotic state corresponds to $`\zeta _s^2(t\mathrm{})\pm 1`$. Thus, for $`\lambda <1/2`$, only one pair of daughter singularities may be shielded—never both— so at least one pair of daughter singularities will impinge on the unit disk. The daughter singularities will also not be shielded when $`\lambda >1/2`$ and $`\beta ^2+4(12\lambda )|\zeta _s|^2<0`$. However, for $`\lambda >1/2`$ and $`\beta ^2+4(12\lambda )|\zeta _s|^2>0`$ it is possible for all the daughter singularities to be shielded, since they lie on a single axis. The daughter singularities can also be completely shielded when $`\lambda =1/2`$. The different possibilities are schematically depicted in Fig. 2.

We have numerically computed the daughter singularity impact time $`t_d`$ for various values of $`\lambda `$ and $`ϵ`$, using initial conditions close to the planar interface, $`|\zeta _s|^2=20`$ and various values of $`\mathrm{Arg}[\zeta _s^2]`$. Figure 3 shows the phase portrait for different values of $`\lambda `$ and $`ϵ`$ with the daughter singularity impact indicated. From the plots it is immediately seen that for $`\lambda <1/2`$ at least one daughter singularity always hits the unit circle, and for $`\lambda 1/2`$ some trajectories are free from daughter singularity impact. In addition, it is observed that for fixed $`\lambda `$ a larger value of $`ϵ`$ causes the daughter singularities to hit in shorter times (or less developed fingers) than a smaller value of $`ϵ`$, and for fixed $`ϵ`$ larger $`\lambda `$ implies larger impact times. We have also checked that the daughter singularity impact occurs well before a finite time singularity, i.e., the impact of a zero of $`z_\zeta `$. Thus, the effect of surface tension is significant well before the curvature in the zero surface tension solution becomes large.

It is noted that the $`\lambda `$ dependence of the daughter singularity impact is consistent with the results of steady state selection theory Shraiman86 ; Hong86 ; Combescot86 ; Tanveer87 . According to selection theory, for small $`B`$ the possible values of $`\lambda `$ are discretized: $`\lambda `$ must satisfy the relation $`\lambda =\lambda _n(B)`$, given to leading order by

$$\lambda _n(B)=\frac{1}{2}\left\{1+(\frac{1}{8}\pi ^2C_nB)^{2/3}\right\},n=0,1,2,\mathrm{}$$
(12)
where $`n`$ parameterizes the branch of solutions. Note that $`\lambda _n>1/2`$ for all $`n`$. The steady finger shape is to leading order a Saffman-Taylor finger, with the above values of $`\lambda _n`$ substituted for the width $`\lambda `$. On the other hand for $`ϵ>0`$ the asymptotic state of (II) is a Saffman-Taylor finger of width $`\lambda `$. From Eq. (12) it is clear there exists a steady solution with width $`\lambda _n(B)`$ close to a Saffman-Taylor finger of arbitrary width $`\lambda >1/2`$. Thus the shielding of the daughter singularity, which leads to the persistence of a Saffman-Taylor solution with $`\lambda >1/2`$ over long times, is consistent with steady state selection theory 444The steady state with $`\lambda >1/2`$ is unstable to tip splitting modes, although specifying the initial value problem in the extended complex plane precludes the presence of noise needed to activate the instability.. In contrast for $`\lambda <1/2`$ there are no nearby steady solutions. Thus, a Saffman Taylor finger with $`\lambda <1/2`$ cannot persist over a long time. We see that the impact of a daughter singularity provides a mechanism for the onset of finger competition, finger widening, and selection of a width $`\lambda >1/2`$.

For $`ϵ=0`$ the scenario is similar, except there is an added class of exact $`B>0`$ solutions. Magdaleno and Casademunt Magdaleno99  have shown that two-finger solutions composed of steadily propagating but unequal fingers do exist for small nonzero $`B`$. The introduction of a small nonzero surface tension selects a discrete set of solutions from the continuum of fixed points of the $`B=0`$ phase portrait. The solutions are parameterized by the total width of the fingers $`\lambda =\lambda _1+\lambda _2`$ and the relative width $`q=\lambda _1/\lambda `$, and the introduction of finite $`B`$ discretizes the possible values of the parameters. In particular, they must satisfy a condition of the form $`\lambda =\lambda _n(B)`$ and $`q=q_{n,m}(B)`$ where $`n`$ and $`m`$ are integers. The expression for $`\lambda _n(B)`$ at lowest order is equivalent to Eq. (12), but with different coefficients $`C_n`$. The shape of these solutions are given to leading order (in the limit $`t\mathrm{}`$) by (II) with allowed value of $`\lambda _n(B)`$ substituted for the width $`\lambda `$. Again, $`\lambda _n(B)>1/2`$, and the consistency between daughter singularity impacts and steady state selection theory follows as above.

We conjecture that the outcome of interfacial shape evolution after the daughter singularity impinges is in general independent of the particular finger on which the impact first occurs i.e., independent of the point at which $`\zeta _d(t)`$ impacts on $`|\zeta |=1`$. More specifically, we surmise that impact on either the shorter (trailing) or larger (leading) finger retards the velocity of that finger, and is accompanied by the widening of the leading finger, so as to maintain a constant fluid flux at infinity. The widened leading finger then shields the trailing finger, preventing it from further growth. Thus, the finger which is leading at the time of the daughter singularity impact ‘wins’ the competition, in the sense that it will evolve for $`t\mathrm{}`$ to the ST steady finger. To examine this conjecture and study the dynamics of finger competition with finite (but small) surface tension we have numerically computed the evolution of an interface with initial conditions given by the conformal mapping Eq. (II) close to the planar interface ($`|\zeta _s(0)|^21`$). The results are reported in the next section.

## IV Numerical Results

Numerical computations have been performed for $`B>0`$, using an initial interface corresponding to the explicit $`B=0`$ solutions discussed in Sec. II. The effect of positive surface tension on this class of solutions is explored for various values of $`ϵ`$ and a variety of initial pole positions.

We employ the numerical method introduced by Hou et al. Hou94  and used in other studies of small surface tension effects in Hele-Shaw flow Siegel96a ; Siegel96b ; Ceniceros99 ; Ceniceros00 . The method is described in detail in Ref. Hou94 . It is a boundary integral method in which the interface is parameterized at equally spaced points by means of an equal-arclength variable $`\alpha `$. Thus, if $`s(\alpha ,t)`$ measures arclength along the interface then the quantity $`s_\alpha (\alpha ,t)`$ is independent of $`\alpha `$ and depends only on time. The interface is described using the tangent angle $`\theta (\alpha ,t)`$ and the interface length $`L(t)`$, and these are the dynamical variables instead of the interface $`x`$ and $`y`$ positions. The evolution equations are written in terms of $`\theta (\alpha ,t)`$ and $`L(t)`$ in such a way that the high-order terms, which are responsible of the numerical stiffness of the equations, appear linearly and with constant coefficients. This fact is exploited in the construction of an efficient numerical method, i.e., one that has no time step constraint associated with the surface tension term yet is explicit in Fourier space. We have used a linear propagator method that is second order in time, combined with a spectrally accurate spatial discretization. Results in this section are specified in terms of the scaled variables

$$\stackrel{~}{t}=\pi t;\stackrel{~}{B}=\pi ^2B;\stackrel{~}{x}=\pi x;\stackrel{~}{y}=\pi y.$$
(13)
instead of the original ones used in previous sections.

The number of discretization points is chosen so that all Fourier modes of $`\theta (\alpha ,t)`$ with amplitude greater than round-off are well resolved, and as soon as the amplitude of the highest-wavenumber mode becomes larger than the filter level the number of modes is increased, with the amplitude of the additional modes initially set to zero. The time step $`\mathrm{\Delta }t`$ is decreased until an additional decrease does not change the solution to plotting accuracy, nor lead to any significant differences in any quantities of interest. In a typical calculation 512 discretization points are initially used, and the initial time step is $`\mathrm{\Delta }t=510^4`$. For small values of surface tension numerical noise is a major problem, and the spurious growth of short-wavelength modes induced by round-off error must be controlled. To help prevent this noise-induced growth at short wavelengths spectral filtering Krasny86  is applied. Additionally, we minimize noise effects and also assess the time at which these effects become prevalent by employing extended precision calculations, as described in the next section.

Our main interest is to uncover the role of surface tension in the dynamics of finger competition. To isolate the features of finger competition from those of width selection, we will concentrate on $`B=0`$ solutions with $`\lambda =1/2`$, the value selected by surface tension in the limit $`B0`$. Since the $`B=0`$ dynamics for $`ϵ=0`$ and $`ϵ0`$ is quite different the numerical results for the two cases will be presented separately.

### IV.1 Solutions with $`ϵ=0`$

We first consider parameter values $`\lambda =1/2`$ and $`ϵ=0`$. A typical set of interfacial profiles is shown in Fig. 4. The initial data is given by the mapping function Eq. (II), with $`\lambda =1/2`$, $`ϵ=0`$, $`d(0)=0`$ and $`\zeta _s^2(0)=20\mathrm{exp}(i\pi /6)`$. With this value of $`\zeta _s^2(0)`$ the initial interface is well inside the linear regime. Evolutions are shown for different values of $`\stackrel{~}{B}`$, and the $`B=0`$ interface evolution is also plotted for comparison. In all these evolutions the filter level is set to $`10^{13}`$, although later we shall make comparisons to profiles computed at higher precision.

For the largest value of surface tension the computed $`B>0`$ and the exact $`B=0`$ solutions first differ appreciably at the seventh curve, corresponding to $`\stackrel{~}{t}3`$. At this point the velocity of the small finger (at the channel sides) begins to decrease and it is clearly left behind when compared with the small finger evolution in the $`B=0`$ solution. Eventually, the advance of the small finger is completely suppressed and the larger finger widens to attain a width close to $`1/2`$ of the channel. For a smaller value of surface tension, for instance $`\stackrel{~}{B}=0.001`$, the evolution displays qualitatively the same behavior. The $`B>0`$ interface differs appreciably from the $`B=0`$ sightly later than before (i.e., at the eighth curve) and the region where the two solutions differ most is to some extent more localized around the small finger than for larger values of $`B`$. Additionally, for this value of surface tension the effect of numerical noise is clearly exhibited in the interfacial profiles. Here the tip-splitting and side-branching activities are a clear effect of numerical noise, as can be easily checked redoing the computation with a different noise filter level.

In order to suppress or delay the branching induced by numerical noise that appears for small values of surface tension it is necessary to use higher precision arithmetic, e.g. quadruple precision (128-bit arithmetic). The filter level can then be reduced by a large amount and the outcome of spurious oscillations is substantially delayed. Figure 5 shows the effect of reducing the filter level to $`10^{27}`$. The $`B=0`$ solution is plotted, as well as the computation with double precision. For $`\stackrel{~}{B}=0.001`$ the branching is totally suppressed, at least for the times we have computed, but for smaller values of $`\stackrel{~}{B}`$ the use of quadruple precision is only able to delay the branching and not totally suppress it. The quadruple precision computation confirms the results observed with lower precision: the introduction of finite (but small) surface tension results in the suppression of the small finger. From Fig. 5 one can also see that for long times, when the interface is clearly affected by numerical noise (in the double precision curve), the noise-induced branching is restricted to the large finger, and the small finger is basically unaffected by noise. This observation suggests that the small finger shape, as well as its tip velocity and tip curvature, can be trusted even when the large finger has developed tip-splittings and side-branchings due to the spurious growth of round-off error.

Figure 6 shows the tip velocity of both fingers versus $`\stackrel{~}{t}`$ for decreasing values of surface tension. It can be seen that the velocity of the large finger is only slightly affected by surface tension, whereas the velocity of the small finger is substantially reduced by the inclusion of finite $`B`$. As $`\stackrel{~}{B}`$ is decreased the tip velocity of the small finger is more faithful to the $`B=0`$ evolution before the daughter singularity impact (shown by a cross), and clearly veers away from the $`B=0`$ velocity later in the evolution, consistent with asymptotic theory. Note that at the smallest value of $`\stackrel{~}{B}`$ the tip velocity of the large finger drastically differs from the $`B=0`$ velocity at late times. This discrepancy is a manifestation of noise effects in the neighborhood of the large finger tip. However, as previously seen, the small finger is basically unaffected by noise at the times we have plotted.

In order to further verify that the daughter singularity impact is responsible for the observed change in the small finger tip speed we follow the scheme introduced in Siegel96b . Define $`t_p`$ as the time when the computed tip velocity differs by $`p`$ from the $`B=0`$ tip velocity. According to asymptotic theory this $`t_p`$ will be a linear function of $`B^{1/3}`$ in the limit $`B0`$ as long as $`p`$ is small enough. Figure 7 shows $`\stackrel{~}{t}_p`$ versus $`\stackrel{~}{B}^{1/3}`$ for various values of $`p`$, and it can be seen that $`\stackrel{~}{t}_p`$ exhibits the predicted behavior. Moreover, we have extrapolated the $`B=0`$ value of $`\stackrel{~}{t}_p`$ using the two points of lowest $`\stackrel{~}{B}`$ and the result is very close to $`\stackrel{~}{t}_d`$, whose value is represented by a cross. We conclude that the impact of the daughter singularity is associated with the dramatic change of the $`B>0`$ solution when compared to the zero surface tension solution, reducing the velocity of the small finger and eventually suppressing it. In contrast, for the $`B=0`$ dynamics the small finger ‘survives’, propagating with the same asymptotic speed as the larger finger. Note that the average interface advances at unit velocity, and a tip velocity below $`1`$ implies that the finger is retreating in the reference frame of the average interface.

In summary then, our numerical results show that the computed interface for $`B0`$ follows the $`B=0`$ evolution for an $`O(1)`$ time interval—roughly corresponding to the daughter singularity impact time—and that at further times the velocity of the small finger decreases while the large finger widens. The small finger eventually comes to a halt and the larger (leading) finger reaches an asymptotic width slightly above $`1/2`$, the width singled out by selection theory. It is noted that for the initial condition we have studied the daughter singularity impact takes place on the tip of the small finger. Therefore, the influence of surface tension on the interface should be significant first around the impact point, that is, the small finger tip. Our numerical results show that in fact this is the case; the initial effect of the daughter singularity impact is to slow and then completely stop the growth of the small finger. Later on, as the singularity cluster centered in $`\zeta _d`$ spreads over the unit circle, the effect of surface tension is felt by the whole interface and the large finger widens to reach the selected width.

We have also studied the finite surface tension dynamics for a more general class of initial conditions. More precisely, we have studied initial conditions of the form $`\zeta _s^2(0)=20\mathrm{exp}(in\pi /12)`$ where $`n=0,\pm 1,\mathrm{},\pm 6`$, and have obtained the same qualitative results as in the case previously studied, namely that the presence of small surface tension suppresses the growth of the finger which is trailing at the time of daughter singularity impact. In order to compare the $`B=0`$ and the $`B0`$ dynamics in a compact and global way we have plotted the phase portrait for $`B=0`$ using the the tip velocities $`v_1,v_2`$ as dynamical variables. In the laboratory frame they read

$`v_1={\displaystyle \frac{1+\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)+\zeta _s^2\overline{\zeta }_s^2}{\zeta _s^2\overline{\zeta }_s^2+\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)/2}}`$ (14a)
$`v_2={\displaystyle \frac{1\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)+\zeta _s^2\overline{\zeta }_s^2}{\zeta _s^2\overline{\zeta }_s^2\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)/2}}.`$ (14b)
Now a comparison between dynamics for $`B=0`$ and $`B0`$ is straightforward since the trajectories can be plotted together and compared. In addition, the tip velocity is a useful variable because it contains geometric information; specifically the inverse of the tip velocity is equal to the width of the finger in the asymptotic ($`t\mathrm{}`$) regime. It is important to note that $`(v_1,v_2)`$ are dynamical variables for the $`B=0`$ problem, so that the plot of the zero surface tension trajectories onto the space $`(v_1,v_2)`$ is a true phase portrait. On the other hand $`(v_1,v_2)`$ are not state variables of the problem with finite surface tension, so in this case we simply obtain a projection onto the $`(v_1,v_2)`$ space of the original $`B0`$ trajectory, which is embeded in the infinite-dimensional phase space of interface configurations.

Figure 8 shows the phase portrait for $`B=0`$ together with the tip velocities obtained from the initial conditions described above for $`\stackrel{~}{B}=0.01`$. From the figure it is evident that the introduction of finite surface tension has substantially changed the global phase dynamics of the problem. Only one $`\stackrel{~}{B}=0.01`$ trajectory connects the planar interface $`(1,1)`$ and the 2ST point $`(2,2)`$, corresponding to the unsteady double Saffman-Taylor finger. Any other $`\stackrel{~}{B}=0.01`$ trajectory ends in one of the two ST finger points, ST(L) at $`(2,0)`$ and ST(R) at $`(0,2)`$. In contrast, the $`(2,2)`$ point, equivalent to the continuum of fixed points present with the $`(\alpha ^{},\alpha ^{\prime \prime })`$ or $`(\mathrm{Re}\zeta _s,\mathrm{Im}\zeta _s)`$ variables, has a finite basin of attraction for $`B=0`$. The introduction of finite surface tension has dramatically changed the zero surface tension $`(v_1,v_2)`$ trajectories, to the extent that the $`B=0`$ phase portrait and the $`B0`$ projection are not topologically equivalent. This result is not a complete surprise, since it was anticipated from the structural instability of the dynamical system governing the evolution of Eq. (II) for $`ϵ=0`$ Magdaleno98 . A more dramatic example of topological inequivalence of phase portraits will be given in the next subsection, when we consider the case $`ϵ0`$.

Although the use of the variables $`(v_1,v_2)`$ has allowed us to project the finite surface dynamics onto the zero surface tension phase portrait this projection has one major limitation: it only considers a local quantity, the tip velocity. We have also considered a projection that takes more global properties of the interface into account. Specifically, given a computed $`B0`$ solution for an initial condition of the form (II), one can use a suitable norm to define a ‘distance’ between the computed interface and the $`B=0`$ interface obtained from the mapping function Eq. (II). We choose this ‘distance’ to be the area enclosed between the two interfaces at a given time. Additionally, we define a projection of the $`B0`$ interface onto the $`B=0`$ phase space (with phase space variables $`(\mathrm{Re}\zeta _s,\mathrm{Im}\zeta _s)`$) by selecting the value of $`\zeta _s`$ that minimizes the ‘distance’ between the two interfaces, with the restriction that the position of the two mean interfaces must be the same. The latter condition ensures that the projection satisfies mass conservation.

Figure 9 shows the $`B=0`$ phase portrait and the corresponding projected evolution for surface tension $`\stackrel{~}{B}=0.01`$. Again, the plot clearly shows that the introduction of finite surface tension modifies the phase portrait of $`B=0`$. The projected trajectories are initially close to the $`B=0`$ dynamics, but for well developed fingers (corresponding to $`|\alpha |1`$) the projection departs from the $`B=0`$ trajectory towards the Saffman-Taylor fixed point, located at $`\alpha ^{}=0,\alpha ^{\prime \prime }=1`$. The projected trajectory only remains close to the corresponding $`B=0`$ trajectory when the latter evolves towards the Saffman-Taylor fixed point. More precisely, the continuum of fixed points present for $`B=0`$ has been removed by surface tension and the Saffman-Taylor fixed point is the universal attractor of the dynamics for finite surface tension.

In Fig. 10 the projection for decreasing values of $`\stackrel{~}{B}`$ is plotted, using the initial condition $`\zeta _s^2(0)=20\mathrm{exp}(i\pi /6)`$. As $`\stackrel{~}{B}`$ is decreased the projected trajectory gets closer to the $`B=0`$ trajectory, but as it approaches the point when the daughter singularity impinges the unit circle (this point is signaled by a cross) the projection departs from the $`B=0`$ trajectory and approaches the Saffman-Taylor fixed point, consistent with asymptotic theory.

### IV.2 Solutions with $`ϵ0`$

The continuum of fixed points present for $`ϵ=0`$ is absent for $`ϵ0`$, but in this case finite-time singularities in the form of zeros of $`z_\zeta `$ impinging on the unit disk do appear for some initial conditions. Therefore, we can expect that the effect of finite surface tension will be somewhat different than for $`ϵ=0`$. Firstly, the presence of surface tension should eliminate finite-time singularities, and secondly, finite $`B`$ could modify the basin of attraction for the two attractors of the $`B=0`$ dynamical system, namely the side Saffman-Taylor finger and the center Saffman-Taylor finger.

To explore this, we have performed computations with $`\lambda =1/2`$ and $`ϵ=0.1`$ with initial conditions $`\zeta _s^2(0)=20\mathrm{exp}(\mathrm{i}n\pi /12)`$ and $`n=0,\pm 1,\mathrm{},\pm 6`$. Initially we set $`\stackrel{~}{B}=0.01`$ and use a value of the noise filter level equal to $`10^{13}`$, which suffices due to the relatively large value of $`\stackrel{~}{B}`$. The easiest way to compare both dynamics, finite $`B`$ and $`B=0`$, is to plot their trajectories in velocity space. Thus, in Fig. 11 the tip velocities $`(v_1,v_2)`$ of the $`\stackrel{~}{B}=0.01`$ computation are plotted together with the tip velocities for $`B=0`$. For arbitrary $`ϵ`$ and $`\lambda `$ the tip velocities of the $`B=0`$ solution read

$`v_1={\displaystyle \frac{1+\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)+\zeta _s^2\overline{\zeta }_s^2}{\zeta _s^2\overline{\zeta }_s^2ϵ(\zeta _s^2+\overline{\zeta }_s^2)+\mathrm{i}\lambda (\zeta _s^2\overline{\zeta }_s^2)(12\lambda )}}`$ (15a)
$`v_2={\displaystyle \frac{1\mathrm{i}(\zeta _s^2\overline{\zeta }_s^2)+\zeta _s^2\overline{\zeta }_s^2}{\zeta _s^2\overline{\zeta }_s^2+ϵ(\zeta _s^2+\overline{\zeta }_s^2)\mathrm{i}\lambda (\zeta _s^2\overline{\zeta }_s^2)(12\lambda )}}.`$ (15b)
From the plot one can see that most $`\stackrel{~}{B}=0.01`$ velocity trajectories follow (at least qualitatively) their $`B=0`$ counterparts, in the sense that they end up in the same fixed point. However, the second, third and fourth trajectories (counting from the upper left trajectory in clockwise direction) differ significantly from their $`B=0`$ counterparts. The second $`\stackrel{~}{B}=0.01`$ trajectory moves apart from the $`B=0`$ solution simply because the latter develops a finite-time singularity, which is regularized by the introduction of finite surface tension. However, the third and fourth trajectories exhibit a quite surprising behavior: the computed interface with $`\stackrel{~}{B}=0.01`$ ends up in a different fixed point than the exact $`B=0`$ solution, despite the fact that the $`B=0`$ solution is smooth for all time and has the asymptotic width that would be selected by vanishing surface tension.

In order to get further insight into this behavior we have computed the evolution for decreasing values of $`\stackrel{~}{B}`$ using the specific initial pole position $`\zeta _s^2(0)=20\mathrm{exp}(\mathrm{i}\pi /6)`$, with $`\lambda =1/2`$ and $`ϵ=0.1`$. Quadruple precision has been used when it has been necessary. Figure 12 shows its evolution for four values of the surface tension parameter, together with the $`B=0`$ solution. The differences between the two interfaces for long times are readily apparent. When $`B=0`$ the finger in the central position stops growing and the side finger wins the competition, whereas for $`B>0`$ we encounter the opposite situation—namely, the central finger surpasses the side finger and wins the competition. For the smaller values of $`B`$ the finger on the sides has not quite stopped growing when the computation is stopped, although its tip speed shows a marked decrease over that for $`B=0`$ and is less than that of the central finger. The side finger tip speed is also decreasing at the final stage of the computation. The tip speed trend in the limit $`B0`$ is further illustrated in Fig. (13). This figure shows the tip speed versus time of each finger for a sequence of decreasing $`B`$. The plot suggests that upon impact of the daughter singularity the side finger velocity levels off and eventually decreases, whereas the velocity of the center finger is nearly unaffected and continues to increase. The trend is indicative of the center finger “winning” the competition in the $`B>0`$ dynamics, while the opposite occurs for $`B=0`$. Finally, it is noted that the influence of surface tension is first felt by the smaller finger, which is the recipient of the daughter singularity impact. Afterwards the leading finger begins to widen, in a manner consistent with the conjecture in Sec. III. Further remarks on this point are made in Sec. V.

The projection method described in the previous section has been also applied to this case, and the results are displayed in Fig. 14 in the particular case $`\stackrel{~}{B}=0.01`$. It can be seen that for most trajectories the projection stays close to the $`B=0`$ curves, even for long times. The daughter singularity impact still leads to $`O(1)`$ differences between the $`B=0`$ and $`B>0`$ solutions, although the impact does not produce changes in the outcome of finger competition. However, as expected some of the trajectories (namely the third and fourth as measured clockwise from the bottom) do indicate significant qualitative differences in the long time evolution. The plot provides a simple depiction of the topological inequivalence of the $`B>0`$ and $`B=0`$ dynamics 555It is noted that the $`B>0`$ interfacial profile and the projection profile do differ significantly in small scale features. This is another indication of the difficulties in using zero surface tension solutions to describe, even qualitatively, the finite surface tension dynamics..

It has been shown that the introduction of a finite $`B`$ has not changed the attractors of the problem, but it has changed their basins of attraction. Interestingly, in the $`B=0`$ case there does not exist a single separatrix trajectory between the two Saffman-Taylor attractors, but rather a finite region, corresponding to the set of trajectories ending in cusps, that acts as an effective separatrix. Since for finite surface tension there are no cusps, it can be assumed that there is a single trajectory that separates the two basins of attraction. Obviously, this trajectory will depend on the value of the surface tension parameter. More precisely, the initial condition $`\zeta _s^2(0)`$ corresponding to the separatrix trajectory will be a function of the surface tension $`B`$. To quantitatively characterize this set of initial conditions we have studied the dependence of the separatrix trajectory in a neighborhood of the planar interface fixed point as a function of $`\stackrel{~}{B}`$, using initial conditions of the form $`\zeta _s^2(0)=20\mathrm{exp}(i\theta )`$. For a given initial condition $`\zeta _s(0)`$ introduce the parameter $`\theta _{sep}(\stackrel{~}{B})`$, defined as the unique value for which the evolution is attracted toward the fixed point ST(L) when $`\theta >\theta _{sep}`$ and to the fixed point ST(R) when $`\theta <\theta _{sep}`$.

Figure 15 shows the plot of $`\theta _{sep}`$ versus $`\stackrel{~}{B}`$, and it is observed that as $`\stackrel{~}{B}`$ decreases, $`\theta _{sep}`$ saturates to a fixed value, namely $`\theta _{sep}(\stackrel{~}{B}0)=0.4843\pm 0.0009`$. It is interesting to compare this value to the position of the separatrix region for $`B=0`$, which is located between $`\theta _+^{B=0}=0.95758`$ and $`\theta _{}^{B=0}=1.04796`$. The separatrix for finite $`\stackrel{~}{B}`$ lays outside and far away from the separatrix region for $`B=0`$, even for vanishing surface tension. Our evidence therefore suggests that any $`B=0`$ trajectory located between the trajectories defined by $`\theta _{sep}(\stackrel{~}{B}0)`$ and $`\theta _+^{B=0}`$ will not describe, even qualitatively, the regularized dynamics in the limit $`\stackrel{~}{B}0`$, since the finger that will ‘win’ the competition under the $`B=0`$ dynamics will ‘lose’ under the $`B0`$ dynamics. Thus, there exists a positive measure set of initial conditions of the form (II) such that the evolution with $`B0`$ cannot be qualitatively described by its evolution under $`B=0`$ dynamics. This is a dramatic consequence of the singular nature of surface tension on the dynamics of finger competition which is not related to steady state selection, but confirms the ideas of the proposed dynamical solvability scenario in Ref. Paune02a .

## V Summary and concluding remarks

The asymptotic theory developed in Refs. Tanveer93 ; Siegel96b  predicts the existence of regions of the complex plane where the zero surface tension solution and the finite surface tension solution differ by $`O(1)`$. These regions are the daughter singularity clusters, and their influence is felt in the physical interface when they are close to the unit circle. Daughter singularities move towards the unit circle, and when their motion is not impeded by other singularities they reach the unit circle in $`O(1)`$ time. When the distance between the daughter singularity and the unit circle is $`O(B^{1/3})`$ the interface can display $`O(1)`$ discrepancies with respect the interface of the $`B=0`$ solutions. However, the asymptotic theory does not predict the nature of the discrepancies caused by daughter singularity impact.

Siegel et al. Siegel96b  showed numerically that the effect of the daughter singularity impact on the tip (in a single-finger configuration with $`\lambda <1/2`$) was a retard in the velocity of the finger accompanied by a widening of the finger. However, this provided small insight on the effect of the impact in multi-finger configurations, where finger competition could be substantially affected by the presence of finite surface tension, as suggested in Refs. Casademunt00 ; Magdaleno98 ; Paune02a .

Since the precise effect of the daughter singularity cannot be established by the asymptotic theory it is necessary to use numerical computation in order to establish the effects of daughter singularity on the dynamics of the interface. We have focused our efforts on uncovering the role of surface tension in the dynamics of two finger configurations, which is the simplest situation exhibiting nontrivial finger competition. Two different types of two-finger zero surface tension solutions have been studied. The first type ($`ϵ=0`$) does not exhibit finger competition when $`B=0`$ but rather contains asymptotic configurations consisting of two unequal steady fingers advancing with the same speed. These two-finger steady state solutions form a continuum of fixed points in the phase space of the corresponding, reduced dynamical system, which is structurally unstable. Numerical computations with small surface tension show that the introduction of a small $`B`$ removes the continuum of fixed points and triggers the competition process which was absent for $`B=0`$ by restoring the saddle-point (hyperbolic) structure of the appropriate multifinger fixed point. The second type ($`ϵ0`$) of two-finger solution we have studied exhibits finger competition for $`B=0`$, but the numerical computation with small $`B`$ has shown that the long time configuration of the computed interface may be qualitatively different from the $`B=0`$ solution for a broad set of initial conditions, in the sense that the finger that ‘wins’ the competition is not the same with and without surface tension. Thus, the presence of surface tension seemingly can change the outcome of finger competition even in configurations that are well behaved and smooth for all time and whose asymptotic width is fully compatible with the predictions of selection theory for vanishing surface tension. This unexpected result shows that surface tension is not only necessary to select the asymptotic width and to prevent cusp formation, but plays also an essential role in multifinger dynamics through a drastic reconfiguration of the phase space flow structure.

Our calculations support the conjecture that impact on either the shorter or larger finger retards the velocity of that finger, and is accompanied by the widening of the larger finger. As a consequence, in general the outcome of finger competition is independent of the particular finger on which the impact first occurs, and the finger which is leading at the time of the daughter singularity impact ‘wins’ the competition. This recipe fails only for interfacial configurations with very similar fingers, when not only the position of the finger (which finger is leading) but also the tip velocities (a trailing finger can have for a certain time a larger velocity than the leading one) at the impact time may play a role.

The main conclusion of the present work is that surface tension is essential to describe multifinger dynamics and finger competition, even when the corresponding zero surface tension evolution is well behaved and compatible with selection theory. That is, we have detected singular effects of surface tension on the dynamics of finger competition that are not directly related to steady state selection. These can be properly interpreted in the context of an extended dynamical selection scenario as described in Ref. Paune02a  where the reconfiguration of phase space flow by surface tension can be traced back to the restoring of hyperbolicity of multifinger fixed points.

###### Acknowledgements.

We acknowledge financial suport from the Dirección General de Enseñanza Superior (Spain), Project No. BXX2000-0638-C02-02 (J. C. and E. P.) and National Science Foundation grants DMS-9704746 and DMS-0104350 (M. S.). E. Pauné also ackowledges financial support from the Departament d’Universitats, Recerca i Societat de la Informació (Generalitat de Catalunya), and the kind hospitality of the Mathematics Department of the New Jersey Institute of Technology.


format: text

text:
# POWER SPECTRA OF X-RAY BINARIES

## 1 INTRODUCTION

The Fourier transform is widely used in timing analysis in astronomy. Through a discrete Fourier transform, a light curve $`x(t_k)`$ can be decomposed into sine wave components in the frequency domain

$$a(f_j)=\underset{k}{}x(t_k)e^{i2\pi f_jk\mathrm{\Delta }t},$$
(1)
where $`a(f_j),j=N/2,\mathrm{},N/21,`$ is the Fourier amplitude of sine wave components at frequency $`f_j=j\mathrm{\Delta }f`$ with $`\mathrm{\Delta }f=1/(N\mathrm{\Delta }t)`$ and $`x(t_k)`$ represents the photon counts during a time interval $`(t_k,t_k+\mathrm{\Delta }t)`$, $`k=1,2,\mathrm{},N`$, with $`t_k=(k1)\mathrm{\Delta }t`$. The Fourier spectrum can be a powerful technique in the search for periodic signals, pulsations, and quasi-periodic oscillations in the light curve of a X-ray source. The Fourier power density

$$p_f(f_j)=|a(f_j)|^2$$
(2)
is used to describe the variability amplitude at different frequency $`f_j`$. A peak with power in excess of the noise distribution power in the Fourier spectrum indicates the existence of a periodic component in the light curve. When estimating the amplitude of a periodic signal from the Fourier power density spectrum, one must assume that the shape of a signal pulse is sinusoidal, which is usually not true for variability due to real physical processes.

Even more care should be taken in interpreting the Fourier spectrum of an aperiodic process in the time domain. It is a prevalent misconception that the Fourier power spectrum is the only way to express the distribution of variability amplitude vs. time scale, and that $`p_f(f_j)`$ is a quantity describing the variability amplitude for the time scale $`1/f_j`$. The actual relation between the Fourier power and the variability process occurring in the time domain is given by Parseval’s theorem:

$$\frac{1}{N}\underset{j=2/N}{\overset{N/21}{}}|a(f_j)|^2=\underset{k=1}{\overset{N}{}}|x(t_k)|^2,$$
(3)
or equivalently,

$$\underset{j}{}p_f(f_j)\mathrm{\Delta }f=\underset{k}{}|r(t_k)|^2\mathrm{\Delta }t,$$
(4)
where $`r(t_k)=x(t_k)/\mathrm{\Delta }t`$ is the counting rate at $`t_k`$. Parseval’s theorem states that the integral of the Fourier power density over the whole frequency range is equal to the variability power of the same process in the time domain. In fact, that is just the reason why the quantity $`p_f(f_j)=|a(f_j)|^2`$ is called power density. However, Parseval’s theorem says nothing about the power density distribution over the time domain. The Fourier power at any given frequency within the limited range, $`01/\mathrm{\Delta }t`$, actually represents a sum of complex amplitudes from an infinite number of aliased frequencies, each contribution of which is reduced by a factor of sinc$`(\pi f\mathrm{\Delta }t)`$, with the sum of the squares of all the factors being unity. In contrast, the power at any given frequency in a time-based spectrum is most sensitive to the rate of change of the Fourier spectrum with frequency near the Nyquist frequency ($`f=1/(2\mathrm{\Delta }t)`$) and its aliases, since the sinc function given above is changing most rapidly for such frequencies. The rms variation vs. time scale of a time series may, in fact, differ substantially from its Fourier spectrum. Even in the simplest case, where the signal is purely sinusoidal with a frequency $`f`$, the Fourier spectrum is a $`\delta `$ function for the continuous Fourier transform or approximately a $`\delta `$ function for the discrete case. Such a function has little power density at frequencies except those near $`f`$, which, although in this case is a clear representation of what physical process is going on, is still not an accurate representation of the variability of the amplitude in a time series. Only when the time step is set equal to or much greater than the period, $`\mathrm{\Delta }t1/f`$, can the corresponding light curve (or pulse profile) of a sinusoidal signal be completely flattened to make the time variation vanish. At time scales shorter than $`1/f`$, however, a variation of intensity still exists in the light curve, and under some circumstances one would wish that power over this region should not be zero. In fact, different but mathematically equivalent representations with different bases or functional coordinates in the frequency domain do exist for certain light curves. The Fourier transform with its trigonometric basis is just one of many possible transforms between the time and frequency domains, and may not necessarily represent the true power distribution of a physical process in the time domain.

One can also argue that any observable physical process always occurs in the time domain, and this is why real variability amplitudes at different time scales are useful in understanding a time varying process, in addition to the understanding gleaned from the conventional Fourier spectrum. An algorithm to calculate the power density spectrum directly in the time domain without using the Fourier transform has been proposed (Li (2001)). We introduce the algorithm and compare the power spectra in the time domain with the Fourier spectra for different kinds of time process model in §1. We have studied the power density spectra of a sample of neutron stars and black hole binaries, in both the frequency and the time domains. The results are presented in §3. A discussion of the potential gain in the understanding the intrinsic nature of physical processes occurring near compact objects via comparison of the two kinds of power density spectrum is given in §4.

## 2 POWER SPECTRA IN THE TIME DOMAIN

The initial definition of the variation power in a light curve $`x(k)`$ is

$$P(\mathrm{\Delta }t)=\frac{\text{Var}(x)}{(\mathrm{\Delta }t)^2}=\frac{\frac{1}{N}_{k=1}^N(x(k)\overline{x})^2}{(\mathrm{\Delta }t)^2}=\frac{1}{N}\underset{k=1}{\overset{N}{}}(r(k)\overline{r})^2,$$
(5)
where the light curve $`x(k),k=1,\mathrm{},N`$, is a counting series obtained from a time history of observed photons with a time step $`\mathrm{\Delta }t`$, $`r(k)`$ is the corresponding rate series, and $`\overline{x}=_{k=1}^Nx(k)/N`$ and $`\overline{r}=_{k=1}^Nr(k)/N`$ are the average counts and count rates, respectively. The unit used for a power $`P`$ in spectral analysis in astronomy is usually rms2, where ’rms’ refers to the analysed time series $`x`$. In the case of $`x`$ being a counting series, the power $`P`$ can be expressed in units of counts2/s2. $`P(\mathrm{\Delta }t)`$ represents an integral of variability power for the region of time scale $`\mathrm{\Delta }t`$. It is easy to prove that $`P(\mathrm{\Delta }t_1)P(\mathrm{\Delta }t_2)`$ if $`\mathrm{\Delta }t_2>\mathrm{\Delta }t_1`$. Thus we can define the power density $`p(\mathrm{\Delta }t)`$ in the time domain as the rate of change of $`P(\mathrm{\Delta }t)`$ with respect to the time step $`\mathrm{\Delta }t`$

$$p(\mathrm{\Delta }t)=\frac{dP(\mathrm{\Delta }t)}{d\mathrm{\Delta }t},$$
(6)
in counts2/s3 or rms2/s. With Eq. (5) and (6), we can calculate the power density spectrum in the time domain for a light curve. In practice, the differential calculus in Eq. (6) can be performed numerically

$$p(\mathrm{\Delta }t)=\frac{P(\mathrm{\Delta }t_1)P(\mathrm{\Delta }t_2)}{\mathrm{\Delta }t_2\mathrm{\Delta }t_1},$$
(7)
with $`P(\mathrm{\Delta }t_1)`$ and $`P(\mathrm{\Delta }t_2)`$ being the two powers at the time scale $`\mathrm{\Delta }t_1`$ and $`\mathrm{\Delta }t_2`$ ($`\mathrm{\Delta }t_2>\mathrm{\Delta }t_1`$), respectively, and $`\mathrm{\Delta }t=(\mathrm{\Delta }t_1+\mathrm{\Delta }t_2)/2`$.

To detect the source signal in a power spectrum against a noise background, we need to know the noise power $`P_{noise}(\mathrm{\Delta }t)`$ of a time series consisting only of noise. For a white noise series where the $`x(k)`$ follows the Poisson distribution, the noise power is

$$P_{noise}(\mathrm{\Delta }t)=\frac{\text{Var}(x)}{(\mathrm{\Delta }t)^2}=\frac{<x>}{(\mathrm{\Delta }t)^2}=\frac{r}{\mathrm{\Delta }t},$$
(8)
where $`<x>`$ is the expectation value of $`x`$, but with units of counts2, which also the variance of the Poisson variable $`x`$, and $`r`$ is the expectation value of the counting rate which can be estimated by the global average of counting rate of the studied observation. The noise power density at $`\mathrm{\Delta }t=(\mathrm{\Delta }t_1+\mathrm{\Delta }t_2)/2`$ is

$$p_{noise}(\mathrm{\Delta }t)=\frac{P_{noise}(\mathrm{\Delta }t_1)P_{noise}(\mathrm{\Delta }t_2)}{\mathrm{\Delta }t_2\mathrm{\Delta }t_1}=\frac{r}{\mathrm{\Delta }t_1\mathrm{\Delta }t_2}.$$
(9)
The signal power density can be defined as

$$p_{signal}(\mathrm{\Delta }t)=p(\mathrm{\Delta }t)p_{noise}(\mathrm{\Delta }t)$$
(10)
and the fractional signal power density as

$$p_{signal}^{}(\mathrm{\Delta }t)=\frac{p_{signal}(\mathrm{\Delta }t)}{r^2},$$
(11)
in units of 1/s or (rms/mean)2/s, where both ’rms’ and ’mean’ refer to the time series.

To study the signal power density in the time domain over a background of noise in an observed photon series, we divide the observation into $`M`$ segments. For each segment $`i`$, the signal power density $`p_{i,signal}(\mathrm{\Delta }t)`$ is calculated by Eq. (10). The average power density of the studied observation is $`\overline{p}_s=_{i=1}^Mp_{i,signal}/M`$ and its standard deviation is $`\sigma (\overline{p}_s)=\sqrt{_{i=1}^M(p_{i,signal}\overline{p}_s)^2/(M(M1))}`$. We can use the statistical methods based on the normal distribution to make statistical inference, e.g. significance test, on $`\overline{p}_s`$. At short time steps, the number $`x`$ of counts per bin may be too small for it to behave as a normally distributed variable. But it is easy to get a large enough total number ($`M`$) of segments from a certain observation period to satisfy the condition for applying the central limit theorem as well as for using the normal statistics on the mean $`\overline{p}_s`$.

To compare the power spectrum in the time domain with the Fourier spectrum for the same process, we study three different kinds of time series.

(1) Periodic signal

We use Eq.(5) and (6) to calculate the power density spectrum in the time domain for a sinusoidal process with a period of 0.5 s. In Figure 1, a piece of the counting rate curve of the studied process is shown in the top panel and the power density distribution of time scale of the sinusoidal signal is shown by the solid line in the middle panel. From the rate curve, we derive a corresponding counting series with a time step 2 ms and use 8192-point FFTs to get the Fourier power densities $`p_f(f_j)`$ and the corresponding densities in the time area $`p_f(\mathrm{\Delta }t_j=1/f_j)=p_f(f_j)f_j^2`$. The Fourier spectrum is shown by the dashed line in the middle panel of Fig. 1. The areas under the two power density curves in Fig. 1 are the same, as dictated by Parseval’s theorem. It is obvious that the Fourier spectrum can not be interpreted as the distribution of variability amplitude vs. time scale. The power being concentrated at the sinusoid frequency $`f`$ does not mean that the intensity variation of the process exists only at the corresponding time scale $`1/f`$. From the light curves with time steps $`\mathrm{\Delta }t=0.05`$ s and 0.28 s, shown in the bottom panel of Fig. 1, we can see that variation of the intensity definitely exists at time scales $`\mathrm{\Delta }t<0.5`$ s, but almost no power does in the Fourier spectrum at those time scales. In contrast with the Fourier spectrum, the power density spectrum determined in the time domain (the solid line in the middle panel of Fig.1) gives a proper description of the variability amplitude distribution of time scales.

Figure 2 shows the power density spectrum (the solid line in the bottom panel) derived in the time domain and the corresponding Fourier spectrum (the dashed line in the bottom panel) for a periodic triangular signal (the top panel). The Fourier spectrum is obtained with 8192-point FFT for a light curve with time step 1 ms. The many peaks at short time scales in the Fourier spectrum (high frequency harmonics) do not mean that there really exist strong variations in the process at those time scales, which are just necessitated by mathematically decomposing the triangular signal into sine waves. As a result, a Fourier spectrum may overestimate power densities at short time scales (high frequency harmonics region) for a periodic signal with pulse shape being far from a sinusoid.

(2) Stochastic Shots

We now consider the power spectra of a signal $`s(t)`$ consisting of stochastic shots. Both the rise and decay fronts of the shots are an exponential with a time constant $`\tau `$ taken uniformly from the range between 5 ms and 0.2 s. The separation between two successive shots is exponentially distributed with an average separation $`d=0.5`$ s. The average signal rate is assumed as $`\overline{r}_s=300`$ cts/s. The peak height of the shot follows a uniform distribution between zero and the maximum. The top panel of Fig. 3 shows a piece of the signal counting curve with time step $`\mathrm{\Delta }t=0.01`$ s and the solid line in the bottom panel of Fig. 3 is the power density spectrum in the time domain expected for the signal $`s(t)`$. A simulated 2000 s light curve $`x(k)`$ with 1 ms time step is produced by a random sampling of the signal curve with Poisson fluctuations plus a white noise with mean rate 5000 cts/s. The middle panel of Fig. 3 shows a piece of the light curve with $`\mathrm{\Delta }t=0.01`$ s obtained from the simulated 1 ms light curve. In calculating the signal power density at a time scale $`\mathrm{\Delta }t`$, The total light curve with time step $`\mathrm{\Delta }t`$ was divided into $`M`$ data segments with $`N=100`$ bins each. If the segments number $`M<100`$ in cases of large time scales, we let $`M=100`$ and decreased the number $`N`$ of data points in each segment accordingly. For each segment of the light curve with time bin $`\mathrm{\Delta }t_1`$, the total powers at two time scales $`\mathrm{\Delta }t_1`$ and $`\mathrm{\Delta }t_2=2\mathrm{\Delta }t_1`$ through Eq. (5) were calculated. The corresponding noise powers were calculated by Eq. (8) with $`r`$ set equal to the average counting rate. The total power density $`p(\mathrm{\Delta }t)`$ and noise power density $`p_{noise}(\mathrm{\Delta }t)`$ at $`\mathrm{\Delta }t=1.5\mathrm{\Delta }t_1`$ can be calculated by Eq. (7) and (9), respectively. Finally, the signal power density in the time domain is $`p_{signal}(\mathrm{\Delta }t)=p(\mathrm{\Delta }t)p_{noise}(\mathrm{\Delta }t)`$. In the bottom panel of Fig. 3, the plus signs mark the average signal power densities at different time scales. For the same light curve with 1 ms time binning, we also calculated the Leahy density $`w(f_j)=2|a_j|^2/a_0`$ for each $`T=4.096`$ s segment, where $`a_j`$ is the Fourier amplitude at frequency $`f_j=j/T`$ determined from a 4096-point FFT. It is well known that the noise Leahy density $`w_{noise}=2`$, so the signal Leahy density can be written as $`w_{signal}(f_j)=w(f_j)2`$ and the Fourier power density of the signal is expressed by $`p_{F,signal}(f_j)=w_{signal}(f_j)a_0/T`$ (Leahy et al (1983), van der Klis (1988)). The dashed line in the bottom panel of Fig. 3 shows the average Fourier power density spectrum of the signal with respect to time scale $`\mathrm{\Delta }t_j=1/f_j`$. As the characteristic time $`\tau `$ of a shot is taken in the range between 5 ms and 0.2 s, there should exist considerable variability over this time scale range as represented by the power density spectrum in the time domain expected for the signal $`s(t)`$ (the solid line in the bottom panel of Fig. 3) and by that obtained from the simulated data including noise (pluses in the bottom panel of Fig. 3). Figure 4 represents the power spectra for another shot model with shorter characteristic time constant $`\tau `$ between 0.5 ms and 2 ms, average separation between two successive shots 3 ms, average signal rate $`\overline{r}_s=500`$ cts/s, and noise rate 5000 cts/s. ¿From Fig. 3 and Fig. 4 we can see that for a random shot series the Fourier spectrum is more or less consistent with the power spectrum in the time domain at time scales greater than the characteristic time scale of the model, but significantly underestimates the power densities at shorter time scales.

(3) Markov process

The Markov process or autoregressive process of first order can describe the character of the variability for many physical processes. A Markov process can be expressed as the following stochastic time series

$$u(k)=au(k1)+ϵ,$$

where $`ϵ`$ is a Gaussian random variable with zero mean and unit variance, and the relaxation time of the process is $`\tau =\mathrm{\Delta }t/\mathrm{log}|a|`$ with $`\mathrm{\Delta }t`$ being the time step. The observed light curve for the signal is

$$s(k)=cu(k)+r_s\mathrm{\Delta }t,$$

where $`r_s`$ is the average rate of the signal. We make a light curve of the signal $`s(k)`$ with $`\mathrm{\Delta }t=0.01`$ s, $`\tau =0.1`$ s, $`r_s=2000`$ cts/s and $`c=2.5`$. A piece of the produced signal light curve is shown in the top panel of Fig. 5. The final observed light curve $`x(k)=s(k)+n(k)`$ with $`n(k)`$ being a Poisson noise with mean rate 5000 cts/s, is shown in the middle panel of Fig. 5. In the bottom panel of Fig. 5, the solid line shows the power density distribution expected for the signal $`s(k)`$, the plus signs indicate excess power densities in the light curve $`x(k)`$ estimated by Eqs. (5)-(10) and the dashed line by FFT. From Fig. 5, we can see, similar to the case of stochastic shot process, that the proposed algorithm of evaluating power densities in the time domain is capable of extracting the power spectrum of the signal from noisy data and that the Fourier spectrum significantly underestimate the power densities at time scales around or shorter than the characteristic time scale of a stochastic process.

## 3 ACCRETING NEUTRON STARS AND BLACK HOLES

As revealed by above study based on simulations, the power density spectrum derived in the time domain can describe the real power distribution with respect to time scale for different time processes. The two kinds of power spectrum, Fourier spectrum and spectrum in the time domain, differ in ways depending on models of the processes. By comparing the two power spectra, more information about the nature of the variability of an object can be extracted. For this purpose we calculate both the power density spectra in the time domain and Fourier spectra for the publicly available data of the Proportional Counter Array (PCA) aboard the Rossi X-ray Timing Explorer (RXTE) for a sample of X-ray binaries, 7 neutron stars and 7 black hole binaries. The PCA observations of X-ray objects included in our study are listed in Table 1. For each analyzed observation, we use the version 4.1 of standard RXTE ftools to extract the PCA data. At 18 time steps between 0.001 s and 2.5 s, we make the corresponding light curves with $`2000`$ s duration in an energy band as noted in Table 1. Then we remove all ineffective data points caused by failure due to the satellite, detector or data accumulation system, and calculate the power density spectra in the time domain, being based on a similar procedure in our simulation. The corresponding Fourier spectra are constructed by using $`2^{12}`$ s ($`244\mu s`$) time resolution light curves divided into parts containing 8192 bins.

Figure 6 shows the results from 6 observations for 5 neutron star binaries: KS 1731-260, 4U 1705-44, GS 1826-24, 4U 0614+091, 4U 1608-522 in the low state and 4U 1608-522 in the high state. The two kinds of power spectrum in the studied accreting neutron stars are generally consistent with each other, at least for the continuum dominated region. The feature is rather complicated for such sources whose Fourier spectra have significant quasi-periodic oscillation (QPO) structure. The left panel of Fig. 7 is from the QPO source, Sco X-1, whose power density spectrum in the time domain also shows a strong QPO structure but the peak shifts to shorter time scales due to the steep slope of the sinc function near the Nyquist frequencies (7 - 12.5 Hz) as more and more of the 6+ Hz QPO is accommodated by the time sampling. The difference between the spectrum of time domain and the Fourier spectrum at short time scales for Cyg X-2 (see the right panel of Fig. 7) may also be caused by QPOs.

For black hole candidates, we first analyze the canonical source Cygnus X-1. On May 10, 1996 (day 131 of 1996), the All-Sky Monitor on RXTE revealed that Cyg X-1 started a transition from the normal low (hard) state to a high (soft) state. After reaching the high state, it stayed there for about 2 months before going back down to the low state (Cui et al (1997)). During this period, 11 pointing observations of Cyg X-1 were made by RXTE. We use one observation of PCA/RXTE for each of the four states of Cyg X-1. The signal power density spectra in the time domain of Cyg X-1 are shown by the plus signs in Fig. 8 and the corresponding Fourier spectra by the dots in Fig. 8. The power spectra of Cyg X-1 in different states (Fig. 8) have a common trait that the Fourier spectra are significantly lower than the corresponding power spectra of time domain in the time scale region of $`\mathrm{\Delta }t<0.1`$ s.

Besides Cyg X-1, the black hole candidates GRS1915+105, GRO J1655-40, GX 339-4, and XTE J1550-564 also demonstrate a significant excess in the power spectra in the time domain in comparison with the corresponding Fourier spectra at time scales shorter than $`0.1`$ s (see Figure 9). But the another analyzed black hole binary, GRS 1758-258, with QPO structures in its Fourier spectrum behaves differently (Fig. 10).

## 4 DISCUSSION

A Fourier power density spectrum presented in the time domain can not be interpreted as the real power density distribution of the physical process studied. In principal, the power densities in the time domain can be derived from a Fourier spectrum only when one knows every power density spectrum in the time domain for each sinusoidal function at all Fourier frequencies and adds them up with weight factors being the Fourier amplitudes. We propose here studying power spectra directly in the time domain. The definition of the power $`P(\mathrm{\Delta }t)`$ (Eq. 5) is based just on the original meaning of rms variation and the power density spectrum $`p(\mathrm{\Delta }t)`$ (Eq. 6) represents the distribution of the variability amplitude vs. time scale. The power density spectrum $`p(\mathrm{\Delta }t)`$ in the time domain obtained from an observation depends only on the intrinsic nature of the signal process and the statistical property of the observed data, as does the Fourier power spectrum in the frequency domain. Our simulation studies show that the proposed algorithm, Eqs. (7), (9) and (10), is capable of extracting power densities of the signal from noisy data (comparing the expected signal spectra, the solid lines in the bottom panels of Figs. 3 and 5 and in Fig. 4, with the spectra calculated from the noisy data, plus signs in corresponding plots). These results indicate that the technique of spectral analysis in the time domain is a useful tool in timing and worth applying in temporal analysis for different sources. From Figs. (1) - (5) one can see that the difference of the Fourier spectrum with the spectrum in the time domain is dependent on the model of time series and sensitive to the characteristic time of a stochastic process. We can then use the difference between two kinds of spectrum to study the intrinsic nature of a studied process.

The power spectra shown in Figs. (8) and (9) for black hole candidates have a common character that the Fourier spectra are significantly lower than the corresponding time-based spectra in the time scale region of $`\mathrm{\Delta }t<0.1`$ s, which is similar to the simulation results for the stochastic shot model (Fig. 3) or the Markov process (Fig. 5). The existence of stochastic shots in X-ray light curves of Cyg X-1 has been noticed for a long time. With an improved searching and superposing algorithm and PCA/RXTE data of Cyg X-1 in different states, Feng, Li & Chen (1998) find that the average shot profiles can be described by exponentials with characteristic time scales $`0.010.1`$ s. Pottschmidt et al (1998) point out that an autoregressive process of first order with a relaxation time $`\tau `$ of about 0.1 s can reproduce approximately the variability of Cyg X-1. Thus, the characteristics of time process in Cyg X-1 revealed from our power spectral analysis is consistent with that from modeling its light curves. The black hole candidates GRS 1915+105, GRO J1655-40, GX 339-4 and XTE J1550-564 also demonstrate characteristics similar to Cyg X-1, indicating that a stochastic process with a characteristic time $`0.1`$ s may be common in accreting black holes. While the absence of the broad-band noise above approximately 100 Hz in black hole candidates has been noticed before (Sunyaev & Revnivtsev (2000)), our analysis shows that the apparent absence in Fourier spectra is caused by the existence of a stochastic process with characteristic time $`0.1`$ s and by the insensitivity of the Fourier technique to detecting rapid variability in a stochastic process.

At the same time, the two kinds of power spectrum are more or less consistent for accreting neutron stars with continuum dominated Fourier spectrum, as shown in Fig. 6. Our simulation results, Figs. (3) - (5), show that for a stochastic process the two kinds of spectrum can be consistent with each other at time scales greater than the characteristic time of the process. Assuming that a significant variability of an accreting neutron star comes from a stochastic process with very short characteristic time constant $`\tau 1`$ ms can explain the consistence of the two kinds of power spectrum observed in neutron star systems. Sunyaev & Revnivtsev (2000) find that the power density spectra of accreting neutron stars with a weak magnetic field have significant broad noise component at the frequency 500-1000 Hz. They suggest that those X-ray transients which demonstrate significant noise in their X-ray flux at frequencies above $`500`$ Hz should be considered neutron stars. Most sources studied in this work (neutron star systems 4U 0614+091, 4U 1608-522, GS 1826-24, 4U 1705-44, KS 1731-260, Cyg X-2, and black hole candidates Cyg X-1, GX 339-4, GRS 1915+105, GRO J1655-40, XTE J1748-288, GRS 1758-258) are also studied by Sunyaev & Revnivtsev (2000) and our results support their claim under the condition that the spectra and noise in their statement are restricted within the Fourier framework. In contrast to the results from Sunyaev & Revnivtsev (2000), our results are based on inferring the characteristic time from the relation between two kinds of spectrum, no matter what the absolute magnitude of power spectrum is. The characteristic feature we find appears in all spectral states of the black hole candidate Cyg X-1 and the same is true with the neutron star binary 4U 1608-522.

The characteristic features in rapid variability for two different kinds of X-ray binary in our sample, which associate significant stochastic processes at the time scale $`0.1`$ s for black hole candidates and $`1`$ ms for neutron stars, are revealed only in continuum or broad noise. The QPO components behave in a more complicated fashion, but can be understood through the sensitivity of the time-based spectrum to variation at 2$`\delta `$t, i.e., at half the Fourier sampling frequency (the Nyquist frequency). It has been found that QPO structures in power spectra can be caused by different kinds of signal, e.g., modulated periodic signals and stochastic autoregressive processes with order $`2`$ (van der Klis (1986)). In the case of the neutron star binary Sco X-1 (the left plot of Fig. 7), the power density spectrum in the time domain reveals the QPO feature surprisingly well, though with a peak shifting to shorter time scales and with worse resolution.

Different processes with essentially different natures could result in almost the same Fourier power spectrum, thus distinguishing them is difficult through timing analysis only with the Fourier technique. This is a reason why we need to develop and apply alternative methods to supplement the Fourier technique in spectral analysis. Our results, though only preliminary, show that simultaneous use of both the Fourier and the time domain methods can help in probing the intrinsic nature of timing phenomena, and further, in distinguishing between different kinds of accreting compact object.

The authors thank the referee for helpful comments and suggestions and Dr. Qu Jinlu for help in data treatment. This work is supported by the Special Funds for Major State Basic Research Projects and the National Natural Science Foundation of China. The data analyzed in this work are obtained through the HEASARC online service provided by the NASA/GSFC.


text:
# Dynamic depletion in a Bose condensate via a sudden increase of the scattering length

## ACKNOWLEDGMENTS

We acknowledge the support from an RGC Earnmarked Grant CUHK4237/01P and the Chinese University of Hong Kong Direct Grant (Grant Nos: 2060148 and 2060150).