text:
# Nonlinear lattice model of viscoelastic Mode III fracture

## I Introduction

The problem of dynamic fracture has received increasing attention in the physics community in the last decade review . The experiments of Fineberg, et. al. fineberg , showing interesting dynamical behavior for large velocity cracks have been at the center of this growing interest. From a theoretical point of view, the singularities at the crack tip associated with continuum treatments of the crack problem make the problem challenging. The presence of these singularities necessitates a treatment of the crack at the microscopic level, where the continuum treatment and its associated singularities are not applicable. One line of attack on this problem, initiated by Slepyan slepyan , has been through the study of lattice models of cracks. These lattice models are simpler than full atomistic simulations mdsim  in that the connectivity of the atoms is specified from the beginning, and so dislocations are excluded.

Much progress has been made in understanding steady-state propagation of cracks in lattice systems slepyan ; mg ; kl1 ; crack3 ; pk . A key simplifying assumption underlying much of this progress has been the assumption of piecewise-linear forces between the particles, so that the particles interact with Hookean springs, which break at some critical extension, reducing the force to zero. With these piecewise-linear interactions, the model admits an analytic solution via the Wiener-Hopf technique. This solution has been carried out both for Mode I and Mode III cracks, for both finite width and infinitely wide systems, with and without dissipation.

A general feature of these solutions is that they are problematic at both very small and very large (roughly above 0.7 of the wave speed) velocities. For small dissipation, the solutions are inconsistent at small velocities, in that bonds on the crack surface which are assumed to crack at time $`t=0`$, say, in fact are seen to crack earlier. For large dissipation, the small velocity solutions are consistent, but exhibit a backward (velocity increasing with decreasing driving) velocity/driving curve, indicating the solutions are unstable. Thus, in any case, stable solutions are not found for small velocities.

At large velocities, the analytic solutions are again inconsistent, this time due to the breaking of additional bonds not on the crack surface. The analytic methods thus are unable to tell us anything about the dynamics beyond this point.

An additional limitation of the piecewise-linear force law is that it complicates the task of constructing a linear stability treatment of the steady-state crack. This is due of course to the discontinuous nature of the force law. For these reasons, we choose in this paper to examine steady-state crack propagation in lattice models with arbitrary force-laws. In particular, we study a family of force-laws parameterized by $`\alpha `$, such that for small $`\alpha `$ the force-law is smooth and the force-law goes over to the piecewise-linear one in the limit of infinite $`\alpha `$. We have previously studied arrest  the behavior of arrested Mode III cracks with this family of force-laws, and found that the range of drivings for which arrested cracks exists narrows sharply as $`\alpha `$ is reduced from infinity. In this paper we examine moving Mode III cracks for varying $`\alpha `$. We find that the behaviour for small velocity is very sensitive to $`\alpha `$. For large velocities, we find that the effect of finite $`\alpha `$ is to convert the inconsistency to a regular linear instability, here of Hopf type. We show that at some distance into the unstable regime, the system adopts a new form of steady-state behavior which breaks the symmetry across the crack surface and has a period-2 structure. We conclude with some observations about the model deep in the unstable regime.

## II Model and General Methodology

In this paper, we study a triangular lattice model of Mode III cracking. We take the force exerted at site $`\stackrel{}{x}_1`$ by the displacement at nearest neighbor site $`\stackrel{}{x}_2`$, with (scalar) displacements $`u_1`$, $`u_2`$, respectively to be

$$f_{1,2}=(u_2u_1)\frac{1+\mathrm{tanh}(\alpha (1\sigma _{1,2}(u_2u_1)))}{1+\mathrm{tanh}(\alpha )}$$
(1)
where $`\alpha `$ is a parameter which controls the smoothness of the potential and $`\sigma `$ is positive if $`\stackrel{}{x}_2`$ lies to the left of $`\stackrel{}{x}_1`$ or is further from the crack plane $`y=0`$. This form has the feature that the force goes to zero at large positive extensions of the nearest-neighbor springs. In the limit of large $`\alpha `$, this force-law approaches a step-function, (with critical spring extension 1), a form introduced by Slepyanslepyan , and the subject of much investigationslepyan2 ; marderliu ; mg ; kl1 ; crack3 ; pk . The lattice has $`2N+2`$ rows in the $`y`$-direction, separated by a distance $`\sqrt{3}/2`$, so that the rows are labeled from $`y=(2N+1)\sqrt{3}/4`$ to $`y=(2N+1)\sqrt{3}/4`$. The displacements $`u`$ of the bottom and top rows at $`y=\pm (2N+1)\sqrt{3}/4`$ are constrained to be $`\pm \mathrm{\Delta }`$. We introduce a Kelvin-type viscosity parameterized by $`\eta `$ via an additional dissipative force

$$g_{1,2}=\eta k_{1,2}\frac{d}{dt}(u_2u_1)$$
(2)
where $`k_{1,2}`$ is an effective spring constant

$$k_{1,2}=f_{1,2}/(u_2u_1)$$
(3)
This form was chosen so as to ensure a purely dissipative force which in the limit of large $`\alpha `$ goes over to the form $`\eta \theta (1(u_2u_1))(\dot{u}_2\dot{u}_1)`$ studied in connection with dissipation in the piecewise-linear model langer ; kl1 ; crack3 ; pk .

These forces define the model. The equation of motion is simply

$$\ddot{u}_\stackrel{}{x}=\underset{\stackrel{}{x}^{}nn}{}\left(f_{\stackrel{}{x},\stackrel{}{x}^{}}+g_{\stackrel{}{x},\stackrel{}{x}^{}}\right)$$
(4)
We shall primarily be interested here in steady-state cracks, where the displacements have the Slepyan form $`u_{x,y}(t)=u_y(tx/v)`$. Furthermore, we will focus on symmetric cracks, such that $`u_y(t)=u_y(t+1/(2v))`$, which is the symmetry appropriate to steady-state cracks in the piecewise-linear model mg . Then, solving for the steady-state crack means solving for $`N`$ functions of time, characterizing the time development of a typical mass point in each row of the lattice. The equations are non-local in time, due to the coupling between different lattice points. As there is no hope of constructing analytic solutions at general $`\alpha `$, we solve this problem numerically.

To proceed, we discretize the time variable in units of some small $`dt`$. The key insight involved in constructing a numerical procedure is noting that the steady-state equations have a banded structure, just as was the case for the piecewise-linear model kl1 . This is due to the fact that mass points are coupled to nearest neighbors, which gives rise via the Slepyan ansatz to a coupling to displacements at a finite time separation. Of course, the equations here are everywhere nonlinear and we need to employ Newton’s method to find a solution. This is nonetheless computationally feasible, since the update step in Newton’s method is a linear problem with the aforementioned banded structure.

Imagine searching for a solution assumed to have some velocity $`v`$. There is one equation of motion for each displacement field $`u_y`$ at each discrete instant in time. In addition, fixing the translation invariance by some condition such as $`u_{\sqrt{3}/4}(0)=1/2`$ gives one additional equation, which we can use to solve for the driving displacement $`\mathrm{\Delta }`$. The only problem is that $`\mathrm{\Delta }`$ is involved in many different equations and so destroys the banded structure of the problem. Two different ways to circumvent this difficulty present themselves. The simpler of the two is to guess, for a given velocity, a value of $`\mathrm{\Delta }`$, and solve a modified system of equations where the equation of motion at $`t=0`$, $`y=\sqrt{3}/4`$ has been dropped, and the time translation has been fixed. The resulting banded problem can be solved economically, and the violation of the dropped equation of motion calculated. This defines a mismatch function, which has a zero at the true value of $`\mathrm{\Delta }`$, which is located by a standard zero-finding routine. A more efficient approach is to realize that solving the entire linear system can be accomplished at essentially no more expense than the fully banded modified system of the first approach. The algorithm for achieving this is presented in Appendix A. The advantage of this approach is that the Newton solver produces a value of $`\mathrm{\Delta }`$ directly, without the need of invoking a outer root-finding routine.

## III The small velocity regime

We have already seen in arrest  that smoothing the form of the force-law effects a dramatic reduction in the window of drivings for which arrested cracks exist. As the moving crack solution arises as a backward bifurcation of the arrested crack, we can expect that small velocity cracks are also extremely sensitive to the smoothness of the force-law. An indication of this can be found in our previous study kl1 ; crack3  of a continuous-$`x`$, discrete-$`y`$ model, which did not have any window of arrested cracks. There the velocity rose linearly from zero as the driving $`\mathrm{\Delta }`$ increased from the Griffith value $`\mathrm{\Delta }_G`$, with a slope inversely proportional to $`\eta `$. The study of small velocity solutions with sharp but not discontinuous force laws should be directly related to experimental data on the onset of propagating cracks in materials such as single crystal silicon silicon .

As mentioned above, we solve the steady-state problem via Newton’s method. We present in Fig. 1(a-c) a graph of $`v`$ versus

$`\mathrm{\Delta }/\mathrm{\Delta }_G`$ for various values of $`\eta `$ and $`\alpha `$, together with the results from the piecewise-linear limit. We chose to present data for $`N=3`$ so as to be able to investigate smaller velocities at reasonable computational cost; the small velocity regime is relatively insensitive to $`N`$ kl1 . A few general trends are evident from these plots. First, the effect of $`\alpha `$ is much more pronounced at small velocity. Second, smoothing the potential by decreasing $`\alpha `$ postpones the onset of the backward branch of the curve to lower velocity. This is consistent with the narrowing of the window of arrested cracks with decreasing $`\alpha `$, since the curve is forced to turn back to meet the end of the arrested crack branch. Decreasing $`\alpha `$ also decreases the amplitude of the oscillations present at small $`\eta `$. Examining the effect of varying $`\eta `$, we see that the larger $`\eta `$’s are less sensitive to $`\alpha `$. This is to be expected, since even in the piecewise-linear model $`\eta `$ reduces the extent of the backward branch which is one of the primary consequences of having smaller $`\alpha `$. This reduction is related to the increase in the size of the process zone with increasing $`\eta `$ crack3 , which smoothes out the lattice structure in a similar manner to that accomplished by smoothing the potential. Thus, in general, decreasing $`\alpha `$ and increasing $`\eta `$ both act to reduce the lattice-trapping effects which give rise to the backward branch. It should be noted that this analysis is consistent with the onset of running cracks directly above $`\mathrm{\Delta }_G`$ in a molecular dynamics simulation using Lenard-Jones potentialspaskin .

Another way to see that reasonably smooth force-laws give rise to essentially $`x`$-continuum behavior is to look at $`\eta v`$ as a function of $`\mathrm{\Delta }`$ for different $`\eta `$’s. In Fig. 2, we present

the data for the case $`\alpha =5`$. We see that the data overlap extremely well, except for the highest velocity data for each value of $`\eta `$. This is an indication that, as in the $`x`$-continuum problem kl1 , the only velocity scales are $`\eta `$ and the wave speed, $`c=\sqrt{3/2}`$. For the piecewise-linear model, on the other hand, where the $`x`$-lattice scale is important, the $`\eta v`$ scaling sander  only sets in for large $`\eta `$, where the process zone is large crack3 .

The other point of interest is the nature of the solutions for small $`\eta `$. In the piecewise-linear model, the solutions on the backward branch are inconsistent at small $`\eta `$, since the underdamping of the forward running waves leads to pre-cracking. This is of course not a problem in our fully nonlinear model, since cracking here is a reversible process. Nevertheless, it is amusing to note that the solutions on the backward branch for small $`\eta `$ do not behave in this manner. They do not crack and reheal. Rather, they exploit the smooth transition in the force law to crack in a slow, monotonic fashion. We show in Fig. 3a

the extension of the bond along the crack surface as a function of time, for $`\alpha =40`$, $`100`$, and $`200`$, along with the analogous result for the piecewise-linear model. We see that the limit of large $`\alpha `$ does not correspond to the behavior of the piecewise-linear model. Not surprisingly, then, it turns out the $`v`$ vs. $`\mathrm{\Delta }`$ curve for large $`\alpha `$ does not converge to the piecewise-linear result. Thus, for example, for the case considered in Fig. 3a, namely N=3, $`v=.1`$, $`\eta =0.25`$, the limiting value of $`\mathrm{\Delta }/\mathrm{\Delta }_G`$ for large $`\alpha `$ is 1.2140, as compared to 1.2778 for the piecewise-linear model. Whereas the bond extension surpasses the critical extension of unity in the piecewise-linear model, and then returns to unity before cracking, the bond extension in the nonlinear model is monotonic, and spends a long time near critical extension. A clearer picture of this emerges in Fig. 3b, where we present the force exerted by the bond as a function of time for the nonlinear model. We see that the bond exploits the nonlinearity of the force-law to crack in a very slow fashion.

## IV The large velocity regime

It is known marderliu ; mg  that the piecewise-linear solution is inconsistent at large velocity due to the cracking of other bonds. This has been extensively studied for the Mode III problem in an infinite square lattice and well as for Mode I in a triangular latticepk . We thus expect that for velocities larger than the critical one at which the inconsistency first appears, the solutions of our nonlinear model will diverge significantly from those of the piecewise-linear model.

To begin, we use our steady-state code to find large velocity solutions of the equations of motion. We present in Fig. 4 data

for $`v`$ versus $`\mathrm{\Delta }/\mathrm{\Delta }_G`$ at three values of $`\alpha `$, together with the analogous results from the piecewise-linear model. We see that for velocities less than about $`0.6c`$, the curves basically coincide, with the smallest $`\alpha `$ curve lying slightly below the others. However, in the range $`0.6c<v<0.7c`$, the data exhibit a marked sensitivity to $`\alpha `$. The low $`\alpha `$ curve is flat in this region, and as $`\alpha `$ is increased, monotonicity is destroyed and a local maximum and minimum are created. These features are completely absent in the piecewise-linear data, indicating that the additional bond breaking, absent by construction in the piecewise-linear solution, is responsible. Clearly, once a local maximum in $`v`$ versus $`\mathrm{\Delta }`$ sets in, as it does for large enough $`\alpha `$, the steady-state solutions are unstable for $`\mathrm{\Delta }`$’s past the maximum. We address this issue in the next section, where we perform a stability analysis of the steady-state solutions. Nevertheless, it is interesting to note that past this band of velocities, the dependence on $`\alpha `$ is again very mild, but very different from the piecewise-linear results.

## V Stability Analysis

Next, we consider the linear stability of the traveling wave state given by the aforementioned Slepyan ansatz $`u_{x,y}(t)=u_y(tv/x)`$. Specifically, we assume a solution of the form

$$u_{x,y}(t)=u_y^{(0)}(\tau )+e^{\omega x}\delta u_y(\tau )$$
(5)
where $`\tau `$ is the traveling wave coordinate $`tx/v`$ and $`u^{(0)}`$ is the previously determined solution. For the stability problem, we take $`\delta u<<1`$ and expand to linear order. Note that the last term can be written in the alternate form $`e^{\omega vt}\stackrel{~}{\delta u_y}(\tau )`$ with

$$\stackrel{~}{\delta u_y}(\tau )=\delta u_y(\tau )e^{v\omega \tau }$$
(6)
Hence, for stability we must have Re $`\omega <0`$, i.e. perturbations must decay in time at a fixed position in the frame moving with the crack tip.

We substitute this assumption into the equation of motion. This leads to the linear problem

$`{\displaystyle \frac{d^2}{d\tau ^2}}\delta u_y(\tau )`$ $`=`$ $`{\displaystyle \underset{nn}{}}G(\stackrel{}{x},\stackrel{}{x}^{})\left(\delta u_y(\tau )\delta u_y^{}(\tau ^{})e^{\omega (xx^{})}\right)`$ (7)
$`+\eta k_{\stackrel{}{x},\stackrel{}{x}^{}}\left(u_y^{(0)}(\tau )u_y^{}^{(0)}(\tau ^{})\right){\displaystyle \frac{d}{d\tau }}\left(\delta u_y(\tau )\delta u_y^{}(\tau ^{})e^{\omega (xx^{})}\right)`$
where we have defined

$$G(\stackrel{}{x},\stackrel{}{x}^{})=f^{}\left(u_y^{(0)}(\tau )u_y^{}^{(0)}(\tau ^{})\right)+\eta k_{\stackrel{}{x},\stackrel{}{x}^{}}^{}\left(u_y^{(0)}(\tau )u_y^{}^{(0)}(\tau ^{})\right)\frac{d}{d\tau }\left(u_y^{(0)}(\tau )u_y^{}^{(0)}(\tau ^{})\right)$$
(8)
and of course $`\tau ^{}=tx^{}/v`$.

To find the allowed values of $`\omega `$, we proceed as follows. First, we impose the asymptotic boundary conditions that at all $`\stackrel{}{x}`$ outside of our explicit lattice points, $`\delta u=0`$. We then pick an arbitrary normalization by fixing $`\delta u_{y=\sqrt{3}/4}(0)=1`$ and simultaneously relax the equation of motion at that same point. This procedure, implemented for some guessed value of $`\omega `$, converts our problem into an inhomogeneous, banded system of linear equations for the complex variables $`\delta u_y(\tau )`$ which can be solved by standard techniques. The missing equation then forms a complex mismatch function whose zeroes then determine the allowed values of the eigenvalue $`\omega `$.

Two other details should be noted. First, we do not impose any symmetry restrictions on our perturbation with respect to reflections about the crack plane. Also, any time we find a root with some specific value of Im $`\omega `$, an equivalent solution exists with

$`\text{Im}\omega `$ $``$ $`2\pi \text{Im}\omega `$
$`\delta u_y`$ $``$ $`\delta u_y^{}\text{(even rows)}`$
$`\delta u_y`$ $``$ $`\delta u_y^{}\text{(odd rows)}`$ (9)
By even and odd, we mean the parity of the lattice units of distance from the row $`y=\sqrt{3}/4`$. This result follows from the fact that the governing equation is real and that the $`x`$ spacing of nearest neighbor points on the rows separated by an even number is $`\pm 1`$, whereas it is $`\pm \frac{1}{2}`$ for odd rows. A corollary of this is that a root at exactly Im $`\omega =\pi `$ has an eigenfunction which is purely real on even rows and purely imaginary on odd ones.

One obvious test of our stability analysis arises from the fact that translation invariance guarantees that there is a root at $`\omega =0`$ with eigenfunction

$$\delta u_y(\tau )=\frac{d}{d\tau }u_y^{(0)}(\tau )/\frac{d}{d\tau }u_{y=\sqrt{3}/4}^{(0)}(0)$$
(10)
Actually, our discretiztion for numerical purposes of the continuous variable $`\tau `$ breaks this exact invariance and hence the mode should be at zero only in the small $`d\tau `$ limit. In the data we present below, the zero mode is found with a typical accuracy of 10<sup>-2</sup>, and the eigenfunction to the same accuracy agrees with its expected form.

In Fig. 5a,

we show the basic result that emerges from our analysis, namely that as the driving displacement is increased, there is a mode which crosses the stability threshold. In Fig. 5b, we show the associated $`Im\omega `$. At the point of instability, $`Im\omega `$ is close to but not exactly equal to $`\pi `$. Were it equal to $`\pi `$, this would be a period-doubling instability, as the perturbative displacements at neighboring $`x`$ values (at the same $`\tau `$) would alternate in sign. More general values of Im $`\omega `$ signify a Hopf bifurcation which in general has a frequency incommensurate with the original frequency associated with moving one lattice spacing. The signature of such a bifurcation should be a slowly oscillating crack tip speed, with the oscillation frequency $`\omega _{tip}=|Im\omega \pi |`$.

In Fig. 6, we show the real part of the eigenfunction for the rows $`y=\pm \sqrt{3}/4`$ as a function of the traveling wave coordinate $`\tau `$. Note that the perturbation is almost, but not exactly, antisymmetric around the crack plane. The eigenvector decays slowly downstream and very rapidly upstream of the crack tip; it is therefore a localized mode connected with the tip dynamics diverging from the pure Slepyan form.

In previous studies of piecewise-linear models (which correspond to the limit $`\alpha \mathrm{}`$), it was noted that the Slepyan solution breaks down above a critical velocity. There, the criterion for this breakdown was connected to the fact that certain bond extensions other than those along the presumed crack line go above the breaking threshold. In our current model where the force-law is analytic, there is no such criterion and the steady-state solution found by our procedure does not need to be checked for any auxiliary condition. On the other hand, we do observe a linear instability. It is therefore of interest to compare these two criteria, i.e. to investigate whether the stability criterion in the large $`\alpha `$ limit goes over smoothly to the inconsistency criterion for the ideally brittle case. To do this, we have plotted in Fig. 7 the horizontal bond extension

along the crack line for a steady-state solution at $`\alpha =50`$ right at the onset of the linear instability. Note that the bond extension is extremely close (roughly of order 1/$`\alpha `$) to what would be the breaking threshold. Thus, the results of the ideally brittle case as to where nontrivial spatio-temporal dynamics of the crack tip sets in are not an artifact of the force-law discontinuity; i.e., smoothing the spring law does not alter the basic conclusion. This is good news, as the ideal case has proven rather amenable to analytic techniques which enable calculations to be done for much larger $`N`$ (even for infinite $`N`$) than is possible for any direct numerical scheme.

One nagging question posed by our results concerns the fact that the instability is so close to but not exactly at the period-doubling point of Im $`\omega =\pi `$. We have investigated whether there is any tendency for the mode to get pinned at the period-doubling point as parameters are varied, as certain limits are reached etc.. For example, in Fig. 8 we show Im

$`\omega `$ at threshold versus $`\alpha `$. As far as we have been able to determine, this pinning does not occur and the generic bifurcation always leads to an incommensurate oscillating tip state. We will compare this prediction with direct numerical simulations of the equations of motion in the next section.

## VI Comparison to Simulation

In order to test the results of the stability analysis, and also to investigate the nature of the dynamics past the instability threshold, we have implemented a direct numerical simulation of the equations of motion. Our system is a triangular lattice, of width $`2N+2`$ whose top and bottom rows are fixed to have displacement $`\pm \mathrm{\Delta }`$. We start with a uniformly strained state, except for the first few left hand columns of mass points. The initial displacement of these points is taken to vary linearly in $`x`$ from $`\pm \mathrm{\Delta }`$ at $`x=0`$ to the that of the uniformly strained state at $`x=10`$. The equations of motion are then stepped forward in time, using a simple Euler scheme. We track the velocity of the crack by monitoring the bond extensions, and noting when they exceed unity, which we take as our criterion of cracking. Note that in this model, cracking is in principle a reversible process and so this criterion is merely a convenient way of keeping track of the dynamics and not an intrinsically important threshold.

As expected, for moderate $`\mathrm{\Delta }`$ the only bonds that “crack” are the diagonal bonds that span the midline. The time between cracking events is constant, and, to best compare with our steady-state calculations above, we fix our time step $`dt`$ in each case so that there are 6 time steps between each cracking event; this can then be directly compared to our steady-state solutions computed with time bandwidth $`n_b1/(vdt)`$ =12. We find that the velocities obtained in this way reproduce extremely well the velocities calculated by our steady-state code.

As we increase $`\mathrm{\Delta }`$ above the critical value for instability calculated in the previous section, we indeed find that the pattern of steady-state cracking breaks down. Some horizontal bonds on the crack surface begin to crack, and the average velocity of the crack falls dramatically, presumably due to the draining of energy by these extra cracking events. It is important to note that these broken horizontal bonds are always behind the crack tip, in accord with the calculations of the piecewise-linear model. It is not clear whether there is in fact a discontinuous drop in the average velocity. The shape of the drop does appear to be insensitive to the $`dt`$ chosen, so that it in any event does not appear to be a numerical artifact of finite $`dt`$.

As a test of our stability calculation, we have computed the Hopf frequency directly from the numerical simulation. We fixed $`\mathrm{\Delta }`$ to lie slightly higher than the critical value for the onset of the instability and measured, at the moment of the breaking of a northeast- southwest diagonal bond, the extension of the bond. Due to the discreteness of the time step, this is somewhat larger than unity. If we had true steady-state propagation, this value would be constant. Instead, for this $`\mathrm{\Delta }`$, we find that the value oscillates about some value, with the magnitude of the oscillation growing in time. The growth we associate with the small positive growth rate of the perturbation in the unstable state, and the frequency of oscillation with the Hopf frequency. In Fig. 9, we show the data. The period of oscillation is essentially 1, (in other words, the time between breakings of this type of bond), which corresponds to a Hopf frequency of essentially $`\pi `$. By numerically fitting a sine wave to the data, we find that the Hopf frequency is approximately 3.272, in very good agreement with the stability calculation result shown in Fig. 5b.

As we increase $`\mathrm{\Delta }`$ further, at some point the average velocity starts to increase again. Examining the solution in this region in more detail, we show in Fig. 10

the time history of three adjacent points on the crack surface. We see that the third trace is identical to the first, and differs significantly from the second. It appears then that the system has settled into a new kind of steady-state solution, more complex than the simple Slepyan form. To test this, we implemented a new steady-state algorithm, designed to allow for this period-2 type solution. We assumed that in each row there are two possible time histories for particles, which alternate. We also did not assume any symmetry across the crack surface. Taking care to organize the storage so that the banded structure is maintained, we employed our banded Newton’s method. Working for the sake of computational convenience at $`N=4`$, the algorithm succeeded in obtaining steady-state solutions where the symmetry between the two time histories (as well as up-down symmetry across the crack surface) is broken. We traced out the velocity-$`\mathrm{\Delta }`$ relation for this solution. We find that the line of solutions is disconnected from the symmetric branch and that in general there are multiple solutions for each $`\mathrm{\Delta }`$. Starting from a solution obtained from using the simulation to obtain an initial guess, we tracked the solution as we varied $`\mathrm{\Delta }`$. The velocity decreased with decreasing $`\mathrm{\Delta }`$ and then turned around and after some wandering headed off to large $`\mathrm{\Delta }`$ with the velocity higher than the original branch. Thus it is the slower branch which is physically realized. This data is presented in Fig. 11, along with the (time-averaged) velocities measured in simulation.

It is reasonable to assume that the existence of a period-doubled solution very close to the original symmetric branch is associated with the fact that the Hopf bifurcation is nearly period-doubling. To check this notion, we note that it was found in the previous section that the Hopf-frequency passed through $`\pi `$ for some particular value of $`\alpha `$. For this value, one expects the period-doubled branch to hit the main branch, corresponding to a higher co-dimension bifurcation. To address this question, we have also traced out the bifurcated solution for $`\alpha =40`$, which is very slightly below this crossing point. The amazingly baroque solution curve is presented in Fig. 12,

together with the symmetric branch. We see that the bifurcated solution curve turns up to approach the symmetric curve, before veering away. Examination of the solutions indicate that the asymmetry is extremely small in this region of close approach. Presumably, then, at the crossing value of $`\alpha `$ the bifurcated solution actually meets the symmetric branch, and the bifurcation is perfect at this point. As the meeting is not on the presumably stable branch of the bifurcated solution, the bifurcation should prove to be subcritical. A detailed study of the nature of the bifurcation should prove mathematically interesting, though not necessarily physically relevant. It should be noted that we have also uncovered another, apparently disconnected, branch of period-2 solutions. These other solutions presumably play an important role in the exact nature of the bifurcation, though again they do not seem to be important for the physics.

Increasing $`\mathrm{\Delta }`$ further, the simulation tracks the bifurcated steady-state solution until at some point yet additional bonds are broken, and the velocity falls below that of the bifurcated solution. The time-dependence of this state becomes more complicated than the simple period-2 structure of the bifurcated steady-state. It is important to point out that the additional bonds are also always located behind the crack tip, which remains on the midline. Thus the additional bond breaking is reminiscent of side-branching in dendritic growth kkl ; i.e., a phenomenon induced by the tip and which is left behind by the growing crack. As opposed to dendritic growth, in Mode III fracture the sidebranching is generated intrinsically by the tip, and is not noise-induced brener . We see, however, no tendency for the tip itself to leave the midline. This can be seen in Fig. 13, where the broken bonds are portrayed for two different values of $`\mathrm{\Delta }`$. Also, the picture shows that as the driving is increased the size of the sidebranches increases, but they are always microscopic, growing to a length of about 8 at the larger $`\mathrm{\Delta }`$. In this latter picture, the crack is moving at an average velocity of $`v/c=.789`$, which is quite fast. Furthermore, the sidebranches are not only short, but the sidebranching period is also on the lattice scale. Thus, the claim of marderliu  that these sidebranches are related to the the experimentally seen micro-branching appears to us doubtful, especially in light of our preliminary work on Mode I cracking mode1new , where the tip dynamics is very different. This issue is worthy of further exploration.

###### Acknowledgements.

The work of DAK is supported in part by the Israel Science Foundation. DAK thanks the Lawrence Berkeley National Laboratory, where this work was initiated, for its hospitality. HL is supported in part by the US NSF under grant DMR98-5735.

## VII Appendix - Solving the almost banded problem

In this appendix, we outline how to efficiently solve the linear problem associated with the Newton’s method treatment of the nonlinear steady-state problem. It is important to set up the problem correctly to achieve an almost-banded structure. We order the variables with the $`y`$ index running more quickly, since we want $`N2L+1`$, the length of the system in the $`x`$ direction. We order the equations similarly, replacing the equation of motion for the site $`x=0`$, $`y=\sqrt{3}/4`$ (where $`x`$ runs from $`L`$ to $`L`$) by the constraint equation $`u=u0`$, where the constant $`u0`$ can be chosen arbitrarily. The equation of motion for the site $`(0,\sqrt{3}/4)`$ is taken to be the last equation. We also note that equations of motion have to be taken for the columns $`L1xL1`$, since the number of linear modes which diverge as $`x\mathrm{}`$ is one greater than the number which diverge for large positive $`x`$ (see kl1 ; crack3 ) due to the presence of the third derivative operator.

The problem then has the structure

$$\left(\begin{array}{cc}_b& V\\ W^T& 0\end{array}\right)\left(\begin{array}{c}\delta U\\ \delta \mathrm{\Delta }\end{array}\right)=\left(\begin{array}{c}A\\ b\end{array}\right)$$
(11)
Here, $`_b`$ is a banded square matrix of size $`N(2L+1)`$ with lower bandwidth $`N(n_b+2)`$ and upper bandwidth $`N(n_b+1)`$ (with $`n_b`$ an integer such that the velocity is $`v=1/(n_bdt)`$). The replacement of the equation for the distinguished site $`(0,\sqrt{3}/4)`$ is necessary to ensure that $``$ is not singular. The length $`N(2L+1)`$ vector $`V`$ contains the derivatives of the equations of motion with respect to $`\mathrm{\Delta }`$ and the vector $`W`$ contains the equation of motion for the distinguished site (which does not depend on $`\mathrm{\Delta }`$). The shift in the displacements $`\delta U`$, and the souurce $`A`$ are also vectors of length $`N(2L+1)`$, and $`\delta \mathrm{\Delta }`$ and $`b`$ are numbers. We can explicitly solve this system in terms of $`_b^1`$, which is easy to compute due to its banded structure. We find

$`\delta \mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{W^T_b^1Ab}{W^T_b^1V}}`$
$`\delta U`$ $`=`$ $`_b^1(A(\delta \mathrm{\Delta })V)`$ (12)
Hence, our entire system can be solved with no more effort than would be required for a fully banded problem.


format: text

text:
# The Discovery of Broad P Cygni X-ray Lines from Circinus X-1 with the CHANDRA High Energy Transmission Grating Spectrometer

## 1 Introduction

The most sensitive indicators of mass loss from cosmic sources are the resonance absorption lines of abundant ions. In the ultraviolet, for example, these lines have been used to study outflows in objects ranging from stars to quasars. Prior to the launch of the Chandra X-ray Observatory (e.g., Weisskopf et al. 2000), however, X-ray spectrometers generally lacked the resolution and sensitivity needed to study resonance line absorption by the outflows from most cosmic objects. Resonance absorption line studies in the X-ray band offer a number of attractive features: (1) the continuum driving the ionization is often directly visible, (2) the atomic physics of the relevant ions is often relatively simple, and (3) unlike ultraviolet photons, X-rays are not easily destroyed in dusty environments.

The mysterious X-ray binary Cir X-1 is thought to contain a neutron star that can radiate at up to super-Eddington levels near the periastron passage (‘zero phase’) of its 16.5-day orbital cycle (e.g., Tennant, Fabian & Shafer 1986; Inoue 1989), and its accretion disk is probably viewed in a relatively edge-on manner (e.g., Brandt et al. 1996). The fact that it can radiate with high luminosity relative to its Eddington luminosity ($`L/L_{\mathrm{Edd}}`$) makes it a natural system in which to expect observable outflows, due to the larger amount of photon pressure available per unit gravitational mass. Indeed, the broad component (up to 2000 km s-1 FWHM) of the strong, asymmetric H$`\alpha `$ line from Cir X-1 has been interpreted as arising in an outflow from regions near the neutron star (e.g., Whelan et al. 1977; Johnston, Fender & Wu 1999); the broad component of the H$`\alpha `$ line is blueshifted relative to the narrow component. Furthermore, the system appears to emit radio jets and is embedded in a synchrotron nebula (e.g., Stewart et al. 1993; Fender et al. 1998). Finally, we note a convergence of opinion that the system is probably a low-mass X-ray binary based on photometry and variability of the optical counterpart (e.g., Stewart et al. 1991; Glass 1994), its orbital period and apparent space velocity (Brandt & Podsiadlowski 1996; Tauris et al. 1999), its correlated X-ray spectral and timing properties (e.g., Shirey, Bradt & Levine 1999), and the observed type 1 X-ray bursts (Tennant et al. 1986).

## 2 Observations and Data Analysis

Cir X-1 was observed using the High Energy Transmission Grating Spectrometer (HETGS; C.R. Canizares et al., in preparation) on Chandra starting at 22:09:50 UT on 2000 February 29 (see Figure 1).111For additional information on the HETGS and ACIS see the Chandra Proposers’ Observatory Guide at http://asc.harvard.edu/udocs/docs. The total integration time was 32 ks, and the observation was continuous. The observation start time corresponds to phase 0.99 using the ephemeris discussed by Glass (1994).

The HETGS carries two types of transmission gratings, the Medium Energy Grating (MEG) and the High Energy Grating (HEG). It allows high-resolution spectroscopy from $``$ 1.2–31 Å ($``$ 0.4–10 keV) with a peak spectral resolution at $`12`$ Å of $`\lambda /\mathrm{\Delta }\lambda 1000`$ for HEG 1st order. The resolution for higher orders improves by a factor of $`n`$ for the $`n`$th order, but the spectral bandpass and efficiency are reduced for these orders. In this analysis we mostly utilize the MEG and HEG 1st order as well as the MEG 3rd order spectra.

The dispersed spectra were recorded with an array of 6 Charge Coupled Devices (CCDs) which are part of the Advanced CCD Imaging Spectrometer (ACIS; G.P. Garmire et al., in preparation).3 To avoid CCD frame ‘dropouts’ due to telemetry saturation, we blocked the 0th order image with a ‘spatial window’ and shut off the two peripheral CCDs (S0 and S5). The omission of the peripheral CCDs could be tolerated because Cir X-1 is heavily absorbed by the interstellar column density ($`N_\mathrm{H}2\times 10^{22}`$ cm-2; e.g., Predehl & Schmitt 1995) and offers hardly any flux above $`15`$ Å (e.g., Morrison & McCammon 1983). The brightness of Cir X-1 also required additional mitigation efforts for ‘photon pile-up’ effects. We applied a ‘subarray’ during the observation that reduced the CCD frame time to 1.7 s. The remaining pile-up effects caused the 1st order spectra below $`10`$ Å (MEG) and $`6`$ Å (HEG) to be depleted of photons, which then populate the higher order spectra at shorter wavelengths. Below 7 Å the MEG 1st order spectrum is so depleted that we will not use it further in this analysis. The MEG 3rd order spectrum is pile-up free from 2.3–9 Å (1.4–5.4 keV). Most of the lines we present below are not or are only moderately affected by pile-up.

Due to a misalignment of the applied subarray during the observation, only half of the dispersed spectra were fully recorded; we only have negative orders in the MEG and positive orders in the HEG complete. The Chandra X-ray Center (CXC) provided aspect-corrected level 1 event lists via its standard pipeline processing.

The determination of the 0th order image position is crucial for the calibration of the wavelength scales because it defines their zero point. We calculated this position by fitting the dispersed images of the MEG and HEG and by determining the intersections of these fits with the 0th order ‘read-out trace’ on the S3 CCD. The fits of the MEG and HEG were consistent to within 0.3 detector pixels, which translates into a zero point accuracy of 0.004 Å for the MEG $`1`$st order and 0.002 Å for the HEG $`+1`$st order. The overall wavelength calibration is good to $`0.1`$% and mostly depends on uncertainties in the positions of the CCD chip gaps, which currently are $`\pm 0.5`$ pixels. For most of the spectral features analyzed here there is only one gap involved leaving us with an overall worst case uncertainty in the scale of 0.008 Å for MEG $`1`$st order, 0.005 Å for HEG $`+1`$st order, and 0.003 Å for MEG $`3`$rd order.

We then processed the events into final event lists using CXC software, and we used custom software and ftools to produce our final grating spectra. After standard grade selection, we have a total of $`1.35\times 10^6`$ events in MEG $`1`$st order, $`1.52\times 10^6`$ events in HEG $`+1`$st order, and $`3.84\times 10^5`$ events in MEG $`3`$rd order. We corrected our final grating spectra for instrumental effects (e.g., energy dependent photon detection efficiency and gaps between the CCDs) using aspect-corrected exposure maps.

From the MEG $`3`$rd order spectrum we determined a mean 2–8 keV flux of $`1.8\times 10^8`$ erg cm-2 s-1. This was calculated by fitting the ASCA continuum model of Brandt et al. (1996) to the 2.3–6.2 Å (2–5.4 keV) spectrum (where pile-up is unimportant) and then extrapolating the result to the 2–8 keV band. The observed flux is within the range that has recently been seen by RXTE (Shirey et al. 1999) at zero phase. Even for a fairly small distance of $`6`$ kpc (see §3 of Stewart et al. 1993 and Case & Bhattacharya 1998), the 2–8 keV luminosity is $`L_{28}\stackrel{>}{}9.6\times 10^{37}`$ erg s-1 (corrected for interstellar absorption); hence for a $`M=1.4`$ $`M_{}`$ neutron star with $`L_{\mathrm{Edd}}=1.3\times 10^{38}(M/M_{})`$ erg s-1 we find $`L_{28}/L_{\mathrm{Edd}}\stackrel{>}{}0.5`$ (note there is certainly substantial X-ray luminosity below 2 keV as well). Figure 2 shows the resulting HEG $`+1`$st and MEG $`1`$st order spectra. These reveal a rich set of highly significant lines from H-like and/or He-like Ne, Mg, Si, S and Fe. The lines are broad ($`\pm 2000`$ km s-1) and show P Cygni profiles (see Figure 3).

Table 1 gives a list of the lines identified with P Cygni profiles; many of these have two or more independent detections (see Figure 3 for an example). We used the line list presented in Mewe, Gronenschild & van den Oord (1985) to determine the line wavelengths and transition types. Most measured wavelengths agree to within 0.01 Å with the predicted wavelengths. Some wavelengths are measured to be slightly lower, which may be attributed to the fact that some of the profiles are broadened and not entirely symmetric. The ‘peak-to-valley’ distances for the P Cygni profiles are a measure of the broadening of the lines. Table 1 shows that the widths are quite significant (0.02–0.10 Å), well above the resolution of the HETGS. While there is evidence for line variability over time, P Cygni profiles are persistently present.

## 3 Discussion

The observed X-ray P Cygni profiles are most naturally interpreted as arising in a high-velocity outflow from the Cir X-1 system. Several properties of the profiles deserve note. First of all, the absorption troughs extend to low velocities; this suggests that we are seeing the region in which the outflow is accelerated, rather than, for example, just a shell of material that has already attained its terminal velocity. Secondly, the strengths of the emission and absorption components are comparable when integrated over the whole observation. Finally, the best-defined profiles look fairly smooth and continuous down to the resolution limit of the HETGS; at least down to this limit there is no evidence for strong shell-like density enhancements or other inhomogeneities in the outflow. The size of the Cir X-1 system is uncertain but is thought to be $`5\times 10^7`$ km (e.g., Tauris et al. 1999), so the characteristic time for the outflow to cross the system is $`0.5`$ days.

As mentioned in §1, an outflow from the Cir X-1 system has already been invoked to explain the broad component of its H$`\alpha `$ line (Johnston et al. 1999). The widths of the X-ray P Cygni profiles are quite consistent with the FWHM of the broad component of H$`\alpha `$, suggesting that these two phenomena may be related; ultraviolet P Cygni profiles and optical H$`\alpha `$ emission are often seen together from stars with high mass loss rates (see §2.3 of Lamers & Cassinelli 1999). Johnston et al. (1999) argued that the broad H$`\alpha `$ line component arises in a strong, anisotropic outflow driven by the radiation pressure from the neutron star. In their model, this outflow must be optically thick (to H$`\alpha `$ photons) to explain the broad blue wing of the line (the redshifted emission being obscured by the outflow itself).

The observed P Cygni profiles are broadly consistent with such an outflow model, although this outflow must be fairly dense to avoid becoming completely ionized by the large X-ray luminosity of Cir X-1. In particular, the ionization parameter $`\xi =L/nr^2`$ must be $`\stackrel{<}{}1000`$ erg cm s-1 so that the ions listed in Table 1 can survive (see Kallman & McCray 1982).222In the definition of $`\xi `$, $`L`$ is the ionizing luminosity, $`n`$ is the electron number density, and $`r`$ is the distance to the ionizing radiation source. If we take the launching radius (i.e., onset radius) of the outflow to be $`r_{\mathrm{launch}}10^5`$ km, so that its observed terminal velocity, $`v_\mathrm{t}`$, is comparable to the escape velocity, $`v_{\mathrm{esc}}`$, at $`r_{\mathrm{launch}}`$, we require $`n_{\mathrm{launch}}\stackrel{>}{}10^{15}`$ cm-3 (we consider the plausible case of outflows launched from larger radii where $`v_\mathrm{t}>v_{\mathrm{esc}}`$ below). For comparison, the expected density of the accretion disk at this radius is $`2\times 10^{18}`$ cm-3, so the outflow-to-disk density contrast seems plausible (compare with Figure 2 of Raymond 1993).333We have calculated the disk density following §5.6 of Frank, King & Raine (1992) with a neutron star mass of $`M=1.4`$ $`M_{}`$, a mass accretion rate of $`\dot{M}=1.5\times 10^{18}`$ gm s-1, and a viscosity parameter of $`\alpha =0.1`$. To allow the survival of H i (for H$`\alpha `$) in the outflow, $`\xi \stackrel{<}{}1`$ erg cm s-1 and hence $`n_{\mathrm{H}\mathrm{I}}\stackrel{>}{}10^{18}`$ cm-3 is required; such a high density outflow would be surprising, so it seems likely that most of the H$`\alpha `$ emission is formed at a larger distance in the outflow or in the accretion disk (see Mignani, Caraveo & Bignami 1997).

Based on the X-ray spectral variability properties of Cir X-1, Brandt et al. (1996) argued that its accretion disk is flared and viewed in a relatively edge-on manner. X-rays emitted in the inner part of such a system can heat the surface of the disk farther out, producing a corona and thermally driven wind (e.g., Begelman, McKee & Shields 1983). The outflow revealed by our X-ray P Cygni profiles might well be this wind, although for this high $`L/L_{\mathrm{Edd}}`$ system radiation pressure acting on lines (in addition to Compton heating) is likely to be important. Raymond (1993) has examined atomic processes in the context of this scenario, and his calculations show that at distances of $`10^5`$$`10^6`$ km there is a significant (subtending $`3`$–20), moderate temperature ($`5\times 10^6`$ K) region where atomic heating and cooling processes dominate over Compton processes. Significant line emission is expected from this region, and in fact most of the lines we list in Table 1 are those predicted to be strong by Raymond (1993). Thus, our favored interpretation for the observed X-ray P Cygni profiles is that they arise in the intermediate temperature region of the wind from an accretion disk viewed in a relatively edge-on manner. Cir X-1 then becomes an X-ray binary analog of a Broad Absorption Line quasar. At smaller radii, wind and coronal material are likely to be heated to the Compton temperature and thus completely ionized. An appealing physical possibility is that the electron scatterer discussed in §4.1 of Brandt et al. (1996) is just this material. Similar highly ionized gas has been invoked in some models for Broad Absorption Line quasars (e.g., the ‘hitchhiking gas’ of Murray et al. 1995).

Despite the general attractiveness of the above picture, we must note a potential difficulty: for a launching radius with $`v_\mathrm{t}v_{\mathrm{esc}}`$, the quantity $`N_{\mathrm{H},\mathrm{launch}}=_{r_{\mathrm{launch}}}^{\mathrm{}}n𝑑r=n_{\mathrm{launch}}r_{\mathrm{launch}}\stackrel{>}{}10^{25}`$ cm-2 is so large that the implied column density through the wind is optically thick to electron scattering \[here we assume that the wind has a radial extent $`r_{\mathrm{launch}}`$ and that $`n=n_{\mathrm{launch}}(r_{\mathrm{launch}}/r)^2`$\]. Most line photons attempting to traverse the wind would then be Compton scattered out of the line. Shielding of the wind from the full X-ray continuum may alleviate this problem by allowing a significant reduction in $`n_{\mathrm{launch}}`$ (see Begelman & McKee 1983). Alternatively, increasing $`r_{\mathrm{launch}}`$ helps because $`N_{\mathrm{H},\mathrm{launch}}r_{\mathrm{launch}}^1`$ when we require $`\xi \stackrel{<}{}1000`$ erg cm s-1 and thus $`n_{\mathrm{launch}}\stackrel{>}{}L/(r_{\mathrm{launch}}^21000\mathrm{erg}\mathrm{cm}\mathrm{s}^1)`$; in this case significant radiation pressure driving may be needed since $`v_\mathrm{t}>v_{\mathrm{esc}}`$. $`r_{\mathrm{launch}}`$ cannot be $`10^6`$ km as the disk probably has an outer radius $`\stackrel{<}{}3\times 10^6`$ km (see Figure 3 of Tauris et al. 1999). Finally, clumping of the wind might also help because when clumping is present $`\xi =Lf/nr^2`$ where $`f`$ is the volume filling factor. We note that Iaria et al. (2000) have recently found evidence for a large ionized column density ($`10^{24}`$ cm-2) along the line of sight even when Cir X-1 is radiating at high luminosity; this material may be the same that makes the P Cygni lines. We are presently examining these issues in further detail.

While we cannot rigorously rule out the possibility that the X-ray P Cygni lines arise in a wind from the companion star, we consider this unlikely. First of all, the large wind velocity ($`\pm 2000`$ km s-1) implied by the profiles would require a high-mass secondary for the system (for an overview see §2.7 of Lamers & Cassinelli 1999), provided the presence of the radiating compact object does not lead to a substantially faster wind from the companion than would otherwise be expected. However, as mentioned in §1, the bulk of the evidence suggests a low-mass X-ray binary nature. Furthermore, the large X-ray luminosity of Cir X-1 should completely ionize an O star wind out to $`\stackrel{>}{}5\times 10^6`$ km (compare with Boroson et al. 1999), while the observed P Cygni profiles suggest that we are seeing the acceleration region of the outflow.

To our knowledge, these are the first reported X-ray P Cygni profiles from an X-ray binary. Hopefully X-ray P Cygni profiles will be identified and studied in other systems to provide geometrical and physical insight into the flows of material near Galactic compact objects.

We thank all the members of the Chandra team for their enormous efforts, and we thank J. Chiang, A.C. Fabian, S.C. Gallagher, S. Kaspi, R.A. Wade, and an anonymous referee for helpful discussions. We gratefully acknowledge the financial support of CXC grant GO0-1041X (WNB, NSS), the Alfred P. Sloan Foundation (WNB), and Smithsonian Astrophysical Observatory contract SV1-61010 for the CXC (NSS).

TABLE 1
X-ray P Cygni lines from Cir X-1     Ion Isoelectronic seq. Predicted $`\lambda `$ Measured $`\lambda `$ $`\mathrm{\Delta }\lambda `$  Name and transitiona (Å)b (Å)c (Å)d    Fe XXVI H-like (Ly$`\alpha `$) 1.78 1.78 0.02  Fe XXV He-like 1.85 1.85 0.02  S XVI H-like (Ly$`\alpha `$) 4.73 4.72 0.03  S XV He-like 5.04 5.04 0.03  Si XIV H-like (Ly$`\beta `$) 5.22 5.21 0.03  Si XIV H-like (Ly$`\alpha `$) 6.18 6.17 0.03  Mg XII H-like (Ly$`\gamma `$) 6.74 6.73 0.04  Mg XII H-like (Ly$`\alpha `$) 8.42 8.40 0.06  Ne X H-like (Ly$`\delta `$) 9.48 9.48 0.09  Fe XXIV Li-like 10.63 10.63 0.10  Ne X H-like (Ly$`\alpha `$) 12.13 12.11 0.06     aWe list the transition here if it can be written compactly; for the other transitions see Tables 1 and 3 of Mewe et al. (1985).
bFrom Mewe et al. (1985).
cMeasured intersection of the P Cygni profile (between the emission maximum and absorption minimum) with a local continuum fit. The measurement accuracy corresponds to about 1 data bin (i.e., $`\pm 0.0025`$ Å).
dMeasured ‘peak-to-valley’ distances for the P Cygni profiles.


text:
# Detection of weakly interacting light vector bosons by coherent scattering

## Abstract

We propose a novel experimental method to produce and detect weakly interacting light vector bosons using coherent processes in refractive media. A light vector boson would be produced by a laser photon scattered at a plane interface between two media of different dielectric properties and will be converted back to photon in the similar scattering process. The effect depends strongly on the indices of refraction of the two media. If incident and recovered photons were traveling through extremely high vacuums, limits on the existence of several light vector bosons could be substantially improved.

PACS: 04.80; 14.80.-j; 29.90.+r

Keywords: Paraphotons; Leptonic photons; Antigravitons

and thanks: E–mail address: Ante.Ljubicic@irb.hr

The existence of weakly interacting light bosons could be of considerable significance in improving our understanding of many phenomena. For example, such pseudoscalar particles as axions  and majorons , could make contributions to the missing mass of the universe and could also be involved in explanations of possible neutrino masses. Hypothetical spin 1 baryonic  and leptonic  photons are candidates to be the quanta of any "fifth" force field, and photon-paraphoton oscillations  could resolve the discrepancy between the observed cosmic background and pure blackbody radiation . The discovery of gravitons, tensor particles with spin 2, or antigravitons, vector particles with spin 1, would be particularly exciting and would be fundamental in the construction of a quantum theory of gravitation.

For the detection of massless or nearly massless weakly interacting light vector bosons we propose a novel experimental method based on coherent processes in refractive media. The behaviour of visible photons in refractive media is determined by Rayleigh scattering and a Feynman diagram of this process is shown in Fig. 1a. Rayleigh scattering is an elastic process and the total amplitude for the process is a coherent sum of scattering amplitudes on many scattering centers. In addition to the Rayleigh scattering there is also an elastic process in which a photon is absorbed and a weakly interacting light vector boson is emitted. A Feynman diagram for this process is shown in Fig. 1b. Weakly interacting light spin 1 particles are denoted with $`\xi `$. In refractive media the scattering amplitudes for this process also add coherently and the weakly interacting vector bosons behave like photons in the sense that they are emitted in both directions at the interface between two refractive media, the backward emission corresponding to reflection and the forward emission to refraction.

We would like to calculate the intensities of reflected and refracted beams for the processes shown in Fig. 1. If $`a_{fi}`$ is the scattering amplitude for a single scatterer placed at the origin of a coordinate system then the scattering amplitude for a single scatterer at an arbitrary position $`\stackrel{}{r}`$ will be $`a_{fi}(\stackrel{}{r})=a_{fi}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r})`$, where $`\stackrel{}{q}=\stackrel{}{k}_i\stackrel{}{k}_f`$ is a vectorial change of the wave vector. If the scatterers are densely packed so that the average distance $`d`$ between two scatterers satisfy the condition $`qd1`$, the total amplitude could be obtained by integrating contributions from all scatterers,

$$A_{fi}=a_{fi}d^3r𝒩(\stackrel{}{r})\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r})$$
(1)
In relation (1) $`𝒩(\stackrel{}{r})`$ is a realistic distribution of the scatterer number density. For two refractive media with well defined scatterer number densities $`N`$ and $`N^{}`$, $`𝒩(\stackrel{}{r})`$ could be represented with the well-known step function $`\theta `$. If the plane surface between two media is in the $`x,y`$ plane, than the density distribution $`𝒩(\stackrel{}{r})=N\theta (z)+N^{}\theta (z)`$ acquires constant values $`N`$ for $`z<0`$, and $`N^{}`$ for $`z>0`$. By integrating relation (1) and using well known representation of $`\theta `$ function

$$\theta (\alpha )=\frac{1}{2\pi i}_{\mathrm{}}^+\mathrm{}\frac{\mathrm{exp}(i\alpha \tau )}{\tau i\epsilon }𝑑\tau $$
(2)
we obtain the result

$$A_{fi}=a_{fi}i(2\pi )^2\delta (q_x)\delta (q_y)\frac{NN^{}}{q_z}$$
(3)
Wave numbers $`k`$ and $`k^{}`$ must be evaluated for particles in $`N`$ and $`N^{}`$ refractive media, and not for particles in free space. For particles of energy $`\omega `$ in media $`N`$ the wave number is $`k=(\omega /c)n`$, where $`n=1+Na_{fi}(0^o)/k_i^2`$ is index of refraction, $`a_{fi}(0^o)`$ is the forward scattering amplitude and $`k_i`$ is the wave number of particles in free space. Eq. (3) does not depend on the type of particles involved in the scattering processes and will be the basis of all our further investigations.

For photons we can obtain the reflection and Snell refraction laws and also the reflection and refraction coefficients $`P_{\gamma \gamma }^{refl}`$ and $`P_{\gamma \gamma ^{}}^{refr}`$. We can extend this to the situation where photons incident from medium N scatter at the interface between the $`N`$ and $`N^{}`$ media, and both photons and weakly interacting vector bosons are emitted. For the process in which weakly interacting vector bosons are emitted into medium $`N^{}`$, in analogy with the similar process in conventional electrodynamics, we define the vector boson refraction coefficient

$$P_{\gamma \xi ^{}}^{refr}=\frac{k_{\xi ^{}z}A_{\xi ^{}\gamma }^{refr}^2}{k_{\gamma ^{}z}A_{\gamma ^{}\gamma }^{refr}^2+k_{\gamma z}A_{\gamma \gamma }^{refl}^2+k_{\xi ^{}z}A_{\xi ^{}\gamma }^{refr}^2+k_{\xi z}A_{\xi \gamma }^{refl}^2}$$
(4)
In Eq. (4) $`A_{\gamma \gamma }^{refl}`$ and $`A_{\gamma ^{}\gamma }^{refr}`$ are the total amplitudes for the reflection and refraction of the photon beam, and are obtained from Eq. (3) by inserting $`q_z=2k_\gamma \mathrm{cos}\theta _i`$ and $`q_z=k_\gamma \mathrm{cos}\theta _ik_\gamma ^{}\mathrm{cos}\theta _\gamma ^{}`$, respectively. $`\theta _i`$ and $`\theta _\gamma ^{}`$ are the angles of incident and refracted photon beams defined relatively to the normal on the plane surface between two media. $`k_\gamma `$ and $`k_\gamma ^{}`$ are wave numbers of photons in $`N`$ and $`N^{}`$ refractive media, and $`k_{\gamma z}`$ and $`k_{\gamma ^{}z}`$ are their projections on the $`z`$-axis. The total scattering amplitudes $`A_{\xi \gamma }^{refl}`$ and $`A_{\xi ^{}\gamma }^{refr}`$ for the processes in which weakly interacting vector bosons are emitted into $`N`$ and $`N^{}`$ refractive media could be obtained from Eq. (3) by inserting $`q_z=k_\gamma \mathrm{cos}\theta _ik_\xi \mathrm{cos}\theta _\xi `$ and $`q_z=k_\gamma \mathrm{cos}\theta _ik_\xi ^{}\mathrm{cos}\theta _\xi ^{}`$, respectively. $`k_\xi `$ and $`k_\xi ^{}`$ are wave numbers of weakly interacting vector bosons emitted at $`\theta _\xi `$ and $`\theta _\xi ^{}`$ angles into $`N`$ and $`N^{}`$ refractive media, and $`k_{\xi z}`$ and $`k_{\xi ^{}z}`$ are their projections on the $`z`$-axis.

The main results of our investigations will become clearer if we make some approximations that will simplify Eq. (4). First we will assume that the rest masses of the weakly interacting vector bosons are very small compared to their total energies. If we further assume that the angle of incidence $`\theta _i0^o`$ we eliminate the dependence of scattering amplitudes $`a_{fi}`$ on initial polarization and scattering angles. For the ratio of relevant scattering amplitudes we can make a reasonably good approximation $`|a_{\xi \gamma }|^2/|a_{\gamma \gamma }|^2\alpha _\xi /\alpha <<1`$. In that case from Eqs. (3) and (4) we obtain

$$P_{\gamma \xi ^{}}^{refr}\left(\frac{\alpha _\xi }{\alpha }\right)\frac{4n_\gamma n_\xi ^{}}{(n_\gamma n_\xi ^{})^2}\left(\frac{n_\gamma ^{}n_\gamma }{n_\gamma ^{}+n_\gamma }\right)^2$$
(5)
where $`n_\gamma `$ and $`n_\gamma ^{}`$ denote indices of refraction for photons, and $`n_\xi ^{}`$ for new vector bosons, respectively; $`\alpha _\xi `$ and $`\alpha `$ are dimensionless constants of interaction for new vector bosons and photons with electrons, respectively.

From Eq. (5) we see that if photons were incident from the high vacuum region, the number of light vector bosons produced in the scattering process at the surface between two refractive media could be very high, because $`n_\xi ^{}1`$ and therefore $`P_{\gamma \xi ^{}}^{refr}1/(n_\gamma 1)^2`$, and for high vacuum $`n_\gamma 1`$. On the other hand the intensities of new vector bosons emitted in the reflection processes are always by an approximate factor $`\alpha _\xi /\alpha `$ lower than the intensities of the reflected photons. Therefore we would be most interested in refraction processes in which photons are absorbed and new vector bosons are emitted in the forward direction. Sensitivity of our experimental method is limited by the masses of light vector bosons produced in the scattering process, because for non-negligible masses the change of the $`z`$-component of the wave vector is $`q_z=k_{\gamma z}k_{\xi ^{}z}\left[n_\gamma 1+(m_\xi c^2/E_\gamma )^2/2\right]`$ for $`\theta _i0^o`$ scattering.

For the experimental investigations of weakly interacting light vector bosons we propose the "shining light through the wall" type of experiment. The initial photon from the laser is incident from the high vacuum region and scatters at a plane surface of a thin glass plate with an index of refraction $`n_\gamma ^{}`$ for photons. The photon beam is then blocked to eliminate everything except the light vector bosons, which pass through because of their extremely weak interaction with ordinary matter. The light boson then interacts at the surface between another glass plate and the high vacuum region to produce a photon, whose detection is the signal for the production of the light boson. A single photon counting system using a low noise photomultiplier could be very efficient for the detection of photons. The sensitivity of the system can be increased if we use $`L`$ and $`L^{}`$ glass plates in front and behind of the blocking piece. In order to minimize the absorption of photons in the glass, the plates should be tilted at the Brewster angle. The expected number of photons recorded in the single photon counting system will be

$$N_c=\epsilon LP_{\gamma \xi ^{}}^{refr}L^{}P_{\xi ^{}\gamma }^{refr}N_{laser}$$
(6)
where $`\epsilon `$ is the efficiency of the photomultiplier, $`N_{laser}`$ is total number of photons emitted from the laser during the experiment, and $`P_{\xi ^{}\gamma }^{refr}`$ is the coefficient for the refraction process, given by an expression similar to Eq. (5), in which weakly interacting vector bosons produce photons.

To estimate the sensitivity of our experiment we will assume that the pulsed laser emits the $`532`$ nm wavelength photons with the total energy of $`30`$ mJ per $`10`$ nsec width pulse at the repetition rate of $`20`$ Hz. We will further assume $`\epsilon =0.25`$, $`L=L^{}=1000`$, $`n_\xi =n_\xi ^{}=1`$,and $`n_\gamma ^{}=1.9`$. Index of refraction $`n_\gamma `$ will be expressed as the function of residual pressure in the experimental apparatus. Using coincidence techniques we can expect the background to be less than $`10`$ counts during $`10^6`$ seconds of laser operation and then from Eqs. (5) and (6) we easily calculate limits for $`\alpha _\xi /\alpha `$ at different vacuum values. Results of our estimates are shown with the full line in Fig. 2.

Present limits for several weakly interacting vector bosons are also shown in Fig. 2. For even modest vacuums of $`10^8`$ Torr the upper limits on coupling constants for paraphotons  with rest masses of $`m_\mathrm{\Gamma }10^7`$ eV could be lowered by an enormous factor of $`10^{28}`$. For very high vacuums of $`10^{12}`$ Torr, which are still available in earth-bound laboratories, the upper limit on coupling constants for leptonic photons  could be improved by a factor of $`10^4`$.

We have also made an estimate for detecting some possible gravitational effects. If a gravity theory would allow the existence of a vector field then antigravitons , i.e. spin 1 gravitational field quanta, could be detected in our apparatus. The probability that the transition between the two electron states will proceed by emission of gravitational, rather than electromagnetic, radiation is typically of order $`\alpha _{grav}/\alpha Gm^2/e^2`$, where $`G`$ is gravitational constant, and $`m`$ and $`e`$ are electron mass and charge respectively. In that case antigravitons could be detected even at vacuums of $`10^8`$ Torr.

The authors gratefully acknowledge A. Peršin for many helpful discussions concerning experimental realities, along with M. Stipčević and B.A. Logan. Early assistence of B. Lakić is also acknowledged. This work was supported by the Ministry of Science and Technology of Croatia.

Figure captions


text:
# The Outbursts of Dwarf Novae

## 1 Introduction

Dwarf novae are a class of cataclysmic variable which undergo regular but aperiodic phases lasting several days during which the system brightness increases by two to four magnitudes. These are the well known normal outbursts of dwarf novae, and they recur on time-scales of weeks to months. There is a well-known bimodal distribution of orbital periods of cataclysmic variables, with a dearth of systems having orbital periods in the range 2.2 $`P_{\mathrm{orb}}`$ 2.8 hours. This is known as the period gap. Dwarf novae which lie above the period gap and only display normal outbursts are classified as U Geminorum systems, after their template. There is a roughly linear period-mass ratio relation for interacting binaries (Frank, King and Raine, 1995). The SU Ursae Majoris systems, which lie below the period gap and have more extreme mass ratios $`q=\frac{M_2}{M_1}<0.25`$ show longer, slightly brighter superoutbursts lasting for ten days or more. These occur in addition to normal outbursts (typically one superoutburst occurs for every 5 to 15 normal outbursts) and show a superimposed periodic variation in brightness at supermaximum (superhumps). In this paper we introduce a numerical method for studying dwarf nova outbursts and apply it to a system of the U Gem class. In a future paper we will present the application to the SU UMa class (Truss et al, 2000).

Historically, two theories have been put forward to explain the outbursts of dwarf novae. Paczyński et al (1969) suggested that a low mass Roche-lobe-filling secondary is potentially unstable and mass transfer may occur quickly enough that its convective envelope would become radiative. The enhanced radiation susequently stabilises the mass transfer and in this way an outburst cycle is initiated. However, observations of the hot-spot (the bright point at which material enters the disc) do not show an increase in brightness during outburst, as would be predicted by such a mass-ransfer instability model (see for example, Rutten et al (1992)). The currently favoured model, and the model which we apply here, is that of an instability in the accretion disc itself. Integration of the density profile normal to the plane of the disc $`\rho (z)`$ yields a relationship between surface density $`\mathrm{\Sigma }`$ and temperature T (or mass-transfer rate) for an annulus in the disc. This is the well-known S-curve and is shown schematically in Figure 1. It represents the locus of points for which the annulus remains in thermal equilibrium. The physical basis which underpins the S-curve is the onset of hydrogen ionization. Hence the upper and lower branches of the curve correspond to the high and low opacity states of ionized and neutral hydrogen. The opacity function $`\kappa (T)`$ is very steep in the range 6000 - 7000 K; such an abrupt change leads to the instability. We can define the local viscous energy generation rate per unit area by

$$Q^+=q^+(z)𝑑z$$
(1)
where $`q^+`$ is the local viscous heat generation rate per unit area, assumed to be concentrated in the centre of the disc. Heat losses can be assumed to be from black-body dissipation from the disc surface. If the surface temperature is $`T_\mathrm{s}`$, then the local energy dissipation rate is

$$Q^{}=2\sigma T_{\mathrm{s}}^{}{}_{}{}^{4},$$
(2)
where each surface of the disc contributes $`\sigma T_{\mathrm{s}}^{}{}_{}{}^{4}`$. In thermal equilibrium, $`Q^+(T_\mathrm{c})=Q^{}(T_\mathrm{s})`$, where $`T_c`$ is the central disc temperature.

An annulus residing on the lower (quiescent) or upper (outburst) branch of the curve subject to perturbations in local surface density can evolve on a viscous time-scale and remain in thermal equilibrium. No stable solution is available to an annulus in the middle branch, however. Here, if the mass-transfer rate (and hence $`T_\mathrm{c}`$) is increased, $`Q^+`$ increases faster than $`Q^{}`$and $`T_\mathrm{c}`$ rises yet further. Therfore, at the inflection points of the curve, the disc will evolve on a thermal time-scale (that is, very rapidly) to the opposite branch. In this way a limit cycle of quiescence and outburst is established. We interpret the upper and lower branches of the curve as high and low viscosity states. The critical surface densities at the points of inflection have been calculated by Cannizzo et al (1988) :

$$\mathrm{\Sigma }_{\mathrm{max}}=11.4R_{10}^{1.05}M_1^{0.35}\alpha _{\mathrm{cold}}^{0.86}\mathrm{gcm}^2$$
(3)
$$\mathrm{\Sigma }_{\mathrm{min}}=8.25R_{10}^{1.05}M_1^{0.35}\alpha _{\mathrm{hot}}^{0.8}\mathrm{gcm}^2$$
(4)
where $`R_{10}`$ is the radius in units of $`10^{10}`$cm, $`M_1`$ is the mass of the primary in solar masses and $`\alpha _{\mathrm{hot}}`$ and $`\alpha _{\mathrm{cold}}`$ are the Shakura-Sunyaev viscosity parameters in the high and low states. Note that these are very close to linear radial dependences. An excellent review of the dwarf nova disc instability is given by Cannizzo (1993a).

We begin with a full explanation of our model and the numerical techniques involved in applying a smoothed particle hydrodynamics (SPH) approach to the problem. In Section 3 we demonstrate that our model is physically consistent with a real dwarf nova and present a wide range of simulated observables.

## 2 Modelling Dwarf Nova Outbursts

### 2.1 The Smooth Particle Hydrodynamics Technique

In this section we discuss substantial developments made to an existing smooth particle hydrodynamics accretion disc code that enable us to model the complete outburst cycle of a dwarf nova. The principal change is the modification of the dissipation term to make it a function of local disc conditions. Smooth particle hydrodynamics (SPH) is a Lagrangian method for modelling the dynamics of fluids. A continuous medium is modelled by a collection of particles that each move with the local fluid velocity. Fluid properties at any given point are determined by interpolating from the particle positions. The interpolation takes the form of a simple summation over the particles with each term weighted according to distance from the point in question. The weighting function is known as the interpolation kernel. For example, the interpolated value for the density at some point $`𝐫`$ in the fluid

$$\rho (𝐫)=\underset{i}{\overset{n}{}}m_iW(𝐫𝐫_i,h),$$
(5)
where the sum is over all the particles. Here, $`m_i`$ and $`𝐫`$ are the mass and position of particle $`i`$ respectively. $`W`$ is the interpolation kernel which has a characteristic length-scale $`h`$, commonly called the smoothing length. For a general introduction to SPH the reader is referred to Monaghan (1992).

In Murray (1996) an SPH code specifically modified for accretion disc problems was described. The key feature of this code was the use of an artificial viscosity term in the SPH equations to represent the shear viscosity $`\nu `$ known to be present in observed discs. The artificial viscosity term introduces, in the continuum limit, a fixed combination of shear and bulk viscosities to the fluid. The viscous force per unit mass,

$$𝐚_\mathrm{v}=\kappa \zeta cL(^2𝐯+2(𝐯)).$$
(6)
$`\kappa `$ is an analytically determined constant that depends upon the kernel. $`\kappa =1/8`$ and $`1/10`$ for the cubic spline kernel in two and three dimensions respectively. $`c`$ is the sound speed, $`𝐯`$ is the fluid velocity and $`L`$ a viscous length scale. In previous work $`L`$ was taken to be equal to the smoothing length $`h`$. Here however we relax that constraint. $`\zeta `$ is a dimensionless parameter. Note that in most SPH papers this parameter is denoted $`\alpha `$ but we have followed Murray (1996) and renamed it to avoid a confusion of subscripts.

In the interior of Keplerian discs we can neglect the velocity divergence and we see that the artificial viscosity term generates a shear viscosity

$$\nu =\kappa \zeta cL.$$
(7)
We can control the shear viscosity throughout the disc by modifying $`\zeta `$ and $`L`$, and so obtain a functional form very similar to the Shakura-Sunyaev form used in accretion disc theory.

Several tests of this code were presented in Murray (1996). The code has been used to look at tidally unstable discs (Murray 1998, 2000), tilted discs (Murray & Armitage 1998) and counter-rotating discs (Murray, deKool & Li, 1999) around accreting pulsars. Kornet & Różycka (2000), using an Eulerian code (quite distinct algorithmically from the SPH code we use), have reproduced several of the results of Murray (1998).

### 2.2 Outbursts

As mentioned in the previous section, quiescence is associated with a very low value of the Shakura-Sunyaev shear viscosity parameter ($`\alpha `$) but also a low temperature. In outburst, angular momentum transport is much more rapid because the temperature is much higher and $`\alpha `$ is much larger. In previous papers (Murray 1998, Armitage & Murray 1998) preliminary attempts were made to model outbursts by instantaneously increasing the value of the shear viscosity throughout the entire disc. Such an approach only enabled us to model the most basic features of an outburst, and we could not of course follow the propagation of state changes through the disc. The principal modification to the code made for this work was to allow the viscosity to change locally in response to disc conditions. This was easy to do, simply requiring that each particle carry a variable $`\zeta `$ that determined its ’viscosity’. To determine the viscosity of the interaction of any particular pair of particles we use the harmonic mean of the two $`\zeta `$ values. Cleary & Monaghan (1999) found that such a form gave good results for heat conduction between materials with vastly different properties.

All that remained was to determine how changes in viscosity were to be triggered. The simplest approach is to let the shear viscosity of each SPH particle be determined by the local surface density $`\mathrm{\Sigma }`$. If, for some region of the disc in the quiescent state, the surface density is less than some critical value $`\mathrm{\Sigma }_\mathrm{c}`$ that region of the disc is assumed to be stably quiescent, and the particle’s shear viscosity remains at a value appropriate for a cool disc. Should $`\mathrm{\Sigma }`$ increase to be greater than $`\mathrm{\Sigma }_\mathrm{c}`$ then we consider that region of the disc to have been ‘triggered’ into the hot state. We then let the particle’s shear viscosity increase to a value appropriate for a hot disc (the details of the transition are described below). The viscosity continues adjusting to its new value even if $`\mathrm{\Sigma }`$ subsequently drops below $`\mathrm{\Sigma }_\mathrm{c}`$. Conversely a hot region of the disc is stable as long as its surface density remains above a second critical value $`\mathrm{\Sigma }_\mathrm{h}<\mathrm{\Sigma }_\mathrm{c}`$. However if a portion of the disc has $`\mathrm{\Sigma }<\mathrm{\Sigma }_\mathrm{h}`$ then that particle’s shear viscosity parameter will reduce to the quiescent value (See figure 2). $`\mathrm{\Sigma }_\mathrm{h}`$ and $`\mathrm{\Sigma }_\mathrm{c}`$ correspond to $`\mathrm{\Sigma }_{\mathrm{min}}`$ and $`\mathrm{\Sigma }_{\mathrm{max}}`$. We preserve their radial dependence but reduce the magnitude of their gradients to reduce the run-time of the code. A full analysis of the scaling is given in Section 3.2.

In previous work (Murray 1998, Armitage & Murray 1998) $`\alpha `$ was set to change instantaneously between quiescent to outburst values. In this work, we assume that when the instability is triggered, the viscosity changes on the thermal time scale. As we have limited knowledge of the mechanism responsible for accretion discs’ high viscosity, we can only propose an appropriate functional form for the transition. We use the hyperbolic tangent function as it allows us to capture both the initial exponential change in $`\nu `$ and a smooth asymptotic approach to its final value. In terms of the Shakura-Sunyaev viscosity parameter, the equation for the transition from quiescence ($`\alpha =\alpha _{\mathrm{cold}}`$) to outburst ($`\alpha =\alpha _{\mathrm{hot}}`$) is

$$\alpha (t)=\frac{(\alpha _{\mathrm{hot}}+\alpha _{\mathrm{cold}})}{2}+\frac{(\alpha _{\mathrm{hot}}\alpha _{\mathrm{cold}})}{2}\mathrm{tanh}(\frac{t}{t_{\mathrm{th}}}\pi )$$
(8)
and the converse, for the transition from $`\alpha _{\mathrm{hot}}`$ to $`\alpha _{\mathrm{cold}}`$ is

$$\alpha (t)=\frac{(\alpha _{\mathrm{hot}}+\alpha _{\mathrm{cold}})}{2}\frac{(\alpha _{\mathrm{hot}}\alpha _{\mathrm{cold}})}{2}\mathrm{tanh}(\frac{t}{t_{\mathrm{th}}}\pi ).$$
(9)
The functional form of the trigger is shown in Figure 3 Of course, $`\mathrm{tanh}`$ doesnt have compact support, so as soon as a transition has been triggered we add a small artificial offset to $`\alpha `$.

In the following section we describe simulations completed with the new variable viscosity. As well as using a density trigger in which single SPH particles could change state, we implemented a more coarse grained approach. The disc was divided into a set of concentric annuli (typically 100 were used). When the mean surface density in a given annulus met the triggering conditions, then all particles in that annulus changed state (again on the thermal time scale). This intermediate, azimuthally smoothed, approach improved the speed of the code by $``$ 50%, and provides a suitable basis for incorporation of a full thermodynamic treatment of the thermal instability. The calculations described in the next section demonstrate that these two approaches to the viscosity switching are consistent with one another.

As it is written, the triggering routine can simply be replaced by a more thermodynamically sophisticated routine. Preliminary calculations have indeed been made. However, as these simulations are proving more computationally demanding, we leave them to a later paper.

## 3 Results

We show that our fast azimuthally-smoothed approximation gives a good representation of the physical behaviour of the system and move on to present a wide range of results and analysis for a simulation of the dwarf nova SS Cygni. SS Cygni is a very well-observed dwarf nova of the U Gem type. It has an orbital period of 0.275130 days and a mass ratio (defined as mass of the mass-losing secondary to that of the accreting white dwarf primary) of 0.59 $`\pm `$ 0.02, with $`M_1`$ = 1.19 $`\pm `$ 0.02 $`\mathrm{M}_{}`$ (Friend et al, 1990).

All the results presented in this paper are for steady state discs, which we define to be those which return to the same quiescent level between outbursts. We take a constant mass-transfer rate of $`10^9\mathrm{M}_{}\mathrm{yr}^1`$ throughout the simulations, following Cannizzo (1993b).

### 3.1 Comparison of triggering regimes

The simulation was performed with both local triggering and azimuthally smoothed triggering as detailed above. In the case of local triggering, two regimes are contrasted; (i)$`\tau _{\mathrm{trigger}}`$ = 250 s, (ii)$`\tau _{\mathrm{trigger}}`$ = 5000 s. For the azimuthally smoothed case, we choose a long trigger timescale of 7500 s. The thermal time-scale $`t_{\mathrm{th}}`$ and the dynamical (Kepler) time-scale $`t_\varphi `$ are related by

$$t_\varphi \alpha t_{\mathrm{th}},$$
(10)
SS Cygni has an orbital period of almost 400 minutes and a Kepler timescale of around 250 seconds at the circularisation radius. This is the radius at which the gas has the same specific angular momentum as it had on passing through the $`L_1`$ point (i.e. as it is injected into the model). It is the radius at which infalling gas would first orbit the primary before formation of an accretion disc. Hence for realistic values of $`\alpha (<<1.0)`$, $`t_{\mathrm{th}}`$ will be of the order of a few thousand seconds, which we use in the azimuthally smoothed code and case(ii).

Figure 4 contains plots of total viscous dissipation against time for the different regimes. The outbursts in both cases are coherent, that is a large region of the disc is transformed to the hot state. It is clear that the azimuthally smoothed result is consistent with a local trigger time-scale somewhere between these two cases. The peak-to-peak time-scale between outbursts is around 30 days in the smoothed code, 20 days in the shorter timescale code and 40 days in the longer time-scale code. Our fast azimuthally smoothed code is consistent with the physically realistic case of a thermal time-scale greater than the Kepler time-scale.

### 3.2 Light curves

Figure 5 shows light curves calculated from the SS Cygni simulation in the U,B,V,J,H and K bands. These are calculated by assuming that the disc is optically thick and each annulus of gas in the disc radiates as a black body. We then simply integrate the Planck function over different wavebands and sum over the annuli to yield the light curves.

The light curves show regular normal outbursts with a comparatively rapid rise to maximum and a slightly longer fall to quiescence. The amplitude of the normal outbursts in the V band varies around 1.5 magnitudes on our diagram, with a rise-time of $``$ 2.5 days and a decay-time of $``$ 4.5 days. SS Cygni is observed to brighten from 12th magnitude in quiescence to 8th magnitude in outburst. Our result is suppressed by the artificially high value of the Shakura-Sunyaev alpha parameter in the two states, which we use to reduce the run-time of the code. We can make an estimate of the compression factor of the amplitude as follows. The dissipation rate of an accretion disc D(R) is directly proportional to the torque exerted by one annulus on another G(R) :

$$D(R)=\frac{G(R)}{4\pi R}\frac{d\mathrm{\Omega }}{dR}$$
(11)
where $`\mathrm{\Omega }`$ is the angular velocity of the gas in the annulus. The viscous torque is

$$G(R)=2\pi R\nu \mathrm{\Sigma }R^2\frac{d\mathrm{\Omega }}{dR}$$
(12)
where $`\nu `$ is the viscosity, so

$$D(R)=\frac{1}{2}\nu \mathrm{\Sigma }(R\frac{d\mathrm{\Omega }}{dR})^2$$
(13)
The dissipation rate is therefore proportional to $`\nu \mathrm{\Sigma }`$. Equation 7 gives the relationship between $`\nu `$ and the alpha-parameter ($`\zeta `$ in Eq.7). In an isothermal disc, the ratio of total dissipation in outburst to that in quiescence is

$$\frac{D_{\mathrm{outburst}}}{D_{\mathrm{quiescence}}}=\frac{(\alpha \mathrm{\Sigma })_{\mathrm{outburst}}}{(\alpha \mathrm{\Sigma })_{\mathrm{quiescence}}}.$$
(14)
We use a simple 1:10 ratio of $`\alpha `$ irrespective of surface density $`\mathrm{\Sigma }`$. The density triggers are set at

$$\mathrm{\Sigma }_{\mathrm{max}}=16.67\frac{R}{a}\mathrm{gcm}^2$$
(15)
$$\mathrm{\Sigma }_{\mathrm{min}}=6.250\frac{R}{a}\mathrm{gcm}^2$$
(16)
where a is the binary separation. The surface density at the peak of outburst will be close to $`\mathrm{\Sigma }_{\mathrm{min}}`$, while the surface density in quiescence just before an outburst will be close to $`\mathrm{\Sigma }_{\mathrm{max}}`$. This is clearly shown in Fig. 9. Therefore, in our isothermal disc simulation,

$$\frac{D_{\mathrm{outburst}}}{D_{\mathrm{quiescence}}}=\frac{\alpha _{\mathrm{hot}}\mathrm{\Sigma }_{\mathrm{min}}}{\alpha _{\mathrm{cold}}\mathrm{\Sigma }_{\mathrm{max}}}3.7.$$
(17)
Real discs, however, are not isothermal. Rewriting Eq.7 in terms of more familiar quantities for a thin Shakura-Sunyaev disc,

$$\nu =\alpha c_\mathrm{s}H\alpha c_\mathrm{s}^2(\frac{R^3}{GM})^{1/2}$$
(18)
where $`c_\mathrm{s}`$ is the sound speed and H is the scale height. In the limit where we can neglect radiation pressure,

$$c_\mathrm{s}^2=\frac{P}{\rho }T_\mathrm{c}$$
(19)
where $`T_\mathrm{c}`$ is the midplane temperature. The critical surface densities will be given by (3) and (4) and in a real disc at fixed radius,

$$\frac{D_{\mathrm{outburst}}}{D_{\mathrm{quiescence}}}=\frac{8.25\alpha _{\mathrm{hot}}T_{\mathrm{c},\mathrm{hot}}\alpha _{\mathrm{hot}}^{0.8}}{11.4\alpha _{\mathrm{cold}}T_{\mathrm{c},\mathrm{cold}}\alpha _{\mathrm{cold}}^{0.86}}$$
(20)
Using Cannizzo’s (1993b) values of $`T_{\mathrm{c},\mathrm{hot}}`$ 60,000K and $`T_{\mathrm{c},\mathrm{cold}}`$ 4,000K, and setting $`\alpha _{\mathrm{cold}}`$=0.01 and $`\alpha _{\mathrm{hot}}`$=0.1, this ratio is about 13.1. The compression factor introduced in dissipation amplitude by the code is about 3.5, bringing our result closer to the eqivalent observed amplitude in SS Cygni.

In the quiescent phase, there appears to be an additional modulation in luminosity. These aperiodic minioutbursts appear in clusters and are particularly blue, suggesting that only the hot inner portion of the disc is involved in producing the additional luminosity. These variations are indeed observed in the V-band light curve of SS Cygni (see, for example, the AAVSO light curves).

It is possible to correlate the number of particles being accreted onto the white dwarf primary with an expected X-ray luminosity for the system by simple consideration of gravitational energy released. For simulation particles of mass $`M_\mathrm{p}`$, the X-ray luminosity produced is

$$L_\mathrm{X}=ϵ\frac{GM_\mathrm{p}}{R}\frac{dM}{dt}$$
(21)
where $`ϵ`$ is an unknown efficiency factor expected to be $`<<1`$.

An X-ray light curve can be built up in this way and is presented in Figure 6 for our simulation. In quiescence and on the decline from outburst, the X-ray emission follows the optical emission very closely; indeed our X-rays appear to be very sensitive to fluctuations in the V-band in quiescence. However, there is a noticeable and significant delay in the onset of the outbursts. This is shown in Figure 7 for each of the four outbursts. Each plot shows 20 days; referring back to Figure 5 we have A : 275-295 days, B : 305-325 days, C : 335-355 and D : 360-380 days. The optical rise leads the X-ray rise by up to a day in all our outbursts. This can be interpreted as the time required for the gas in the part of the disc which initially goes into outburst to migrate through the inner disc before it is accreted onto the primary. The measurement of the delay is in good agreement with observation; Jones $`\&`$ Watson (1992) found a typical delay of 0.5 - 1.1 days from EXOSAT observations. The rise-time is also seen to be much faster in X-rays; we find a rise-time of $``$12 hours in X-ray and $``$2.5 days in the optical. The climb in V-band emission is a gradual process as more and more gas in the disc is transformed to the high state, whereas the onset of enhanced accretion onto the white dwarf itself is instantaneous, triggered when the additional material from the outbursting disc arrives at its surface.

### 3.3 Disc Analysis

In this section we follow the evolution of an accretion disc through an outburst cycle. The light curves give much information about the physical processes taking place in the disc during quiescence and outburst, but a much clearer representation is obtained by following the dynamical behaviour of the accretion disc itself. Figure 8 is a montage of images of the accretion disc through a normal outburst. We find that even in quiescence (t=334.2), a small fraction of the inner disc remains permanently locked in the hot state. The linear relation between the critical surface density triggers and the disc radius ensures that particles in the inner disc are continually cycling between the low and high states (the triggers are very close together here and also very low - these conditions are easily satisfied by particles close to the primary). The rise to outburst is rapid and is intitiated in the inner part of the disc. A density wave (directly analagous to a heating wave in a full thermodynamic treatment) moves both outwards and inwards, transforming a large fraction of the disc to the hot, high viscosity state. Subsequently the disc falls (less rapidly) back into quiescence by another density wave (analagous to a cooling wave). This pattern is repeated for all the normal outbursts. Approximately the same fraction of the disc is seen to reach the high state in each outburst. The maximum extent of the hot portion is seen to vary to a small degree, but appears ultimately to be limited to a radius near the position of the spiral shock arms. This can be readily observed in the light curve - some outbursts reach a slightly higher peak luminosity than others, although the maximum variation is no more than a few percent.

The outburst is, of course, triggered by the local surface density in the disc and as the gas in the disc is accreted inwards on a viscous timescale there is a corresponding change in that local surface density. This response is shown in Figure 9. The disc clearly drains during the outburst phase until the local surface density reaches the lower trigger level, after which the outburst dies away and the disc is replenished. These density profiles are extremely similar to those found by Stehle (1999) in a full thermodynamic one-dimensional treatment (not SPH).

Figure 10 shows a disc in minioutburst, and confirms our observation from the light curve in section 3.2 that it is the inner disk which is hot in a minioutburst. The outburst is initiated very near the inner edge of the disc and never propagates very far. The reason for this can be seen in Figure 9. Just before a normal outburst (top panel), the surface density profile follows the line of the upper trigger very closely. Any part of the disc that goes into outburst will lead to a coherent normal outburst throughtout most of the disc. After a normal outburst (bottom panel), the only part of the disc near the upper trigger is near the disc centre. All the remainder of the disc has been drained. Hence while the centre of the disc can go into outburst, the majority of the disc cannot. This is what happens in a minioutburst.

Disc spectra, temperature profiles and luminosity profiles clearly show the contrast in the physical properties of the disc in the different states. The spectra, in Figure 11, have been computed for the disc in outburst, quiescence and minioutburst. The blue, hot inner disk behaviour in minioutburst becomes clear in comparison with the quiescent spectrum.

The luminosity and temperature profiles in Figure 12 have been calculated on the rise to a normal outburst, and correspond to the discs in the first three frames of Figure 8. All the profiles show an increased region of viscous dissipation (and temperature) near the outer edge of the disk, where the surface density peaks as matter enters the disc from the accretion stream; this is the hot-spot. This normal outburst is apparently triggered at a radius of approximately 0.15a (near the circularisation radius - 0.11a in SS Cyg) and the heating wave is seen to propagate in both directions. We find that in quiescence the disc is typically approximately isothermal at radii greater than 0.15a at temperatures around 4000 K. In outburst we find temperatures ranging from 6000 K to 12000 K in the inner half of the disc and 4000 K to 6000 K in the outer half. Bobinger et al (1997) have perfomed eclipse mapping of the dwarf nova IP Peg , which has strikingly similar parameters to SS Cyg ($`q=0.58,P=0.158d`$), on the decline from outburst. They found an inner disc at 7000 K to 9000 K and an outer disc with temperatures declining to a quiescent level of 3000 K to 4000 K. Our simulated results are in good agreement with these observations. Figure 12 also shows that the hot-spot is consistently around 1500 K hotter than the surrounding regions, but the increase in luminosity from these regions is much more marked - the viscous dissipation is 30 times higher here than in the part of the disk just inside them. The accretion rate from the secondary is constant in this simulation, so we find a constant hot-spot temperature. It is also interesting to note, however, that the luminosity generated in the hot-spot of the quiescent disc is comparable to that generated by the very hot inner part of the disc.

## 4 Discussion

We have performed the first two-dimensional treatment of dwarf nova outbursts. The method gives light curves in different wavelength bands, disc spectra,density and temperature profiles which all agree well with observations and full thermodynamic analyses.

The advantages of using a two-dimensional code to model dwarf nova outbursts are immediately obvious. Tidal forces are not approximated as in one dimensional codes and therefore much more detail is revealed in the results of these current simulations. One dimensional simulations have always predicted extremely regular outburst behaviour, both in recurrence pattern and outburst profile. There is now much more scope for exploring the behaviour of these systems; the simulated observational characteristics presented here are just a selection of the possible areas which can be explored in future work.

We have applied these methods successfully to a study of the SU Ursae Majoris systems (Truss et al, 2000) in addition to the dwarf novae with less extreme mass ratios considered here. A two-dimensional code is ideal for exploring tidal effects and instabilities in such systems.

Observations of dwarf novae are by no means comprehensive or continuous, but this situation should be remedied in the near future with more telescope time dedicated to them. In consequence, it will be possible to test several of the results in this paper. Although our knowledge of the behaviour of dwarf novae in the V-band is good thanks to the efforts of organisations such as the AAVSO, there is little data available in other bands, so new observations will be able to probe the inner regions of the disc and (via X-ray observations) boundary layer. Our findings relating to a hot inner disc and the occurrence of minioutbursts should have some testable observational consequences.

The next stage of development of any numerical model such as this should be twofold : the incorporation of full thermodynamics into the simulation and the extension to three dimensions. This is not a difficult task with the present scheme, just a lengthy one in computational terms. The increasing availability of parallel machines and improvements in CPU perfomance should make these refinements a workable proposition.

## Acknowledgments

Research in theoretical astrophysics at the University of Leicester is supported by a PPARC rolling grant. Many of the calculations were performed using GRAND, a high performance computing facility funded by PPARC and based at the University of Leicester. The authors are grateful to Dr R. Stehle, Prof. J. Pringle, Prof. A. King and an anonymous referee for helpful discussions or comments. MRT acknowledges a PPARC research studentship and a William Edwards Charitable Trust bursary for postgraduate study.