Mon. Not. R. Astron. Soc. 000, 1–18 (2005) Printed 18 November 2021 (MN LATEX style file v2.2)
Comprehensive simulations of superhumps
Amanda J. Smith1,⋆, Carole A. Haswell1, James R. Murray2,
Michael R. Truss3 and Stephen B. Foulkes1
1Department of Physics & Astronomy, The Open University, Walton Hall, Milton Keynes, MK7 6AA, UK
2Department of Astrophysics & and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia
3Department of Physics, Durham University, South Road, Durham, DH1 3LE, UK
ABSTRACT
We use 3D SPH calculations with higher resolution, as well as with more realistic
viscosity and sound-speed prescriptions than previous work to examine the eccentric
instability which underlies the superhump phenomenon in semi-detached binaries. We
illustrate the importance of the two-armed spiral mode in the generation of super-
humps. Differential motions in the fluid disc cause converging flows which lead to
strong spiral shocks once each superhump cycle. The dissipation associated with these
shocks powers the superhump. We compare 2D and 3D results, and conclude that 3D
simulations are necessary to faithfully simulate the disc dynamics. We ran our sim-
ulations for unprecedented durations, so that an eccentric equilibrium is established
except at high mass ratios where the growth rate of the instability is very low.
We collate the observed data on superhumps. Our improved simulations give a
closer match to the observed relationship between superhump period excess and bi-
nary mass ratio than previous numerical work. The observed black hole X-ray transient
superhumpers appear to have systematically lower disc precession rates than the cat-
aclysmic variables. This could be due to higher disc temperatures and thicknesses. No
high-resolution 3D disc with mass ratio q > 0.24 developed superhumps, in agreement
with analytical expectations.
The modulation in total viscous dissipation on the superhump period is over-
whelmingly from the region of the disc within the 3 : 1 resonance radius. The pre-
cession rates of our high resolution 3D discs match the single particle dynamical pre-
cession rate at 0.87R3:1. As the eccentric instability develops, the viscous torques are
enhanced, and the disc consequently adjusts to a new equilibrium state, as suggested
in the thermal-tidal instability model. We quantify this enhancement in the viscosity,
which is ∼ 10 per cent for q = 0.08. The disc motions can be described as superpo-
sitions of the S(k,l) modes, and the disc executes complex standing wave dynamics
which repeat in the inertial frame on the disc precession period. We characterise the
eccentricity distributions in our accretion discs, and show that the entire body of the
disc partakes in the eccentricity.
Key words: accretion, accretion discs — hydrodynamics — instabilities — methods:
numerical — binaries: close — novae, cataclysmic variables
1 INTRODUCTION
Cataclysmic variables (CVs) are semi-detached binaries with
a Roche-lobe filling low-mass donor star, mass M2, transfer-
ring matter onto a white dwarf (WD) primary, mass M1,
via an accretion disc. The SUUMa-type dwarf novae (DNe)
are short-period cataclysmic variables which display two
distinct modes of outburst. The normal outbursts are at-
tributed to a thermal–viscous limit–cycle between low and
⋆ E-mail: amanda.smith@open.ac.uk
high viscosity states (Osaki 1974; Hōshi 1979; Meyer &
Meyer-Hofmeister 1981), and consequently between low and
high mass transfer states (see Lasota 2001, for a review).
The larger amplitude and longer-lasting superoutbursts are
characterised by a periodic photometric modulation known
as superhumps. Superhumps are attributed to an eccentric
apsidally-precessing accretion disc (Vogt 1982). The super-
hump period, Psh, is a few percent longer than the orbital
period, Porb: the orientation of the mass donor star relative
to the progradely precessing eccentric disc repeats on Psh.
Whitehurst (1988) and Lubow (1991a) explained that a disc
c© 2005 RAS
http://arxiv.org/abs/0704.1519v1
2 A. J. Smith et al.
which encounters a 3:1 eccentric inner Lindblad resonance
with the tidal potential of the secondary star may become
eccentric and precess. A mass ratio q = M2/M1 . 1/4
is required for the tidal truncation radius, Rtides, to lie
outside the 3:1 resonance radius, R3:1 (Paczynski 1977).
The thermal–tidal instability (TTI) model (Osaki 1989) at-
tributes the increased brightness and duration of superout-
bursts over normal outbursts to an enhanced viscous torque
acting once the disc becomes eccentric.
Superhumps arise in several guises, summarised by Pat-
terson et al. (2002a). ‘Common’ or ‘normal’ superhumps en-
sue after the onset of superoutburst in SUUMa systems.
They are powered by the periodically varying tidal interac-
tion which modulates dissipation in the disc (e.g. Foulkes
et al. 2004). ‘Late’ superhumps sometimes follow common
superhumps and are roughly anti-phased to them. The mod-
ulation in energy dissipation at accretion stream’s impact on
the non-axisymmetric disc powers late superhumps (Rolfe,
Haswell & Patterson 2001). ‘Negative’ superhumps, with Psh
slightly less than Porb, are sometimes observed simultane-
ously with the more usual ‘positive’ superhumps, and may
be related to retrograde precession of a warped accretion disc
(Patterson et al. 1993b; Murray & Armitage 1998; Foulkes,
Haswell & Murray 2006). ‘Permanent’ superhumps are seen
in high mass transfer rate systems: the nova-like variables,
old novae and some AMCVn systems. Some low mass X-ray
binaries (LMXBs, the neutron star and black hole analogues
of CVs) also show superhumps (Bailyn 1992; Haswell 1996;
O’Donoghue & Charles 1996). In LMXBs optical emission
arises overwhelmingly from the reprocessing of X-rays, and
superhumps arise from a modulation in reprocessing caused
by the changing solid angle subtended by the tidally flexing
eccentric disc (Haswell et al. 2001). Recently superhumps
were reported in the microquasar GRS1915+105, which has
an orbital period exceeding 30 days (Neil, Bailyn & Cobb
2006).
In all the above cases, the fractional superhump excess,
ǫ = (Psh − Porb)/Porb varies with q, with |ǫ| increasing with
higher values of q. The exact relationship has proved difficult
to determine. We have performed unprecedentedly compre-
hensive numerical simulations of apsidally precessing accre-
tion discs, and we present them in the context of previous
numerical work, the observational data and salient analyt-
ical theory. Section 2 describes our simulations and exam-
ines the growth rates of the eccentric instability; the en-
hanced viscous torques; the superhump light curves which
result from the eccentric instability; and compares 2D and
3D simulation results. In section 3 we use two methods to
quantify the eccentricity distributions in our simulated discs.
Section 4 focuses on the ǫ − q relationship. In section 5 we
discuss our findings. Section 6 gives a summary list of our
principle conclusions.
2 SPH SIMULATIONS OF A PRECESSING
ACCRETION DISC
SPH is a Lagrangian method which models fluid flow as a set
of moving particles. A detailed review is given by Monaghan
(1992). SPH simulations by Murray (1998) provide consid-
erable support for the TTI model, showing that the energy
released from a disc that has become tidally unstable is suffi-
cient to account for the excess luminosity of a superoutburst.
Foulkes et al. (2004) carried out 2D SPH simulations of a
binary system with mass ratio 0.1. They show an eccentric,
non-axisymmetric precessing disc of changing density, which
is continuously flexing and relaxing on the superhump pe-
riod. Very clear too in the surface density maps are tightly
wrapped spiral density waves which extend from the outer-
most regions to small radii. They produce shear and dissi-
pation in the outer disc, and propagate angular momentum
outwards, allowing disc gas to move inward.
2.1 Simulation details
For the calculations presented here, an SPH code is used
which has been designed specifically for accretion disc prob-
lems. Detailed description of the code can be found in Mur-
ray (1996, 1998). The code is normalised in units of a, the
binary separation, for distance; Mt = M1+M2 for mass; and
Porb/(2π) for time. This code has since been updated to in-
clude adaptive spatial resolution, allowing the SPH smooth-
ing length λ to vary in both space and time (Murray, de Kool
& Li 1999). Here, λ is set to a maximum value of 0.005 a. The
code has also been extended to three dimensions (Murray &
Armitage 1998).
Simulations were run for a range of mass ratios, de-
tailed in Table 1. The simulations were built up from zero
mass and a single particle was injected at the L1 point each
timestep ∆t, so simulating the mass-transfer stream from
the secondary. We ran simulations at different mass resolu-
tions, where ∆t is between 0.01Ω−1orb and 0.0025 Ω
orb. This
latter value is at resolution higher than previous calculations
(Murray 1996, 1998; Foulkes et al. 2004), and we find our
results stable to a further reduction in ∆t. The number of
particles in the accretion disc in these high-resolution simu-
lations is of order 100 000, and the average number of ‘neigh-
bours’, i.e. the average number of particles found within 2λ
of each other that are used in the SPH update equations, is
between 20 and 30.
The stream boundary conditions at L1 are a function
of q, calculated by Lubow & Shu (1975) from perturbation
analysis of L1. We follow their calculations to determine
the direction prograde of the binary axis at which particles
are to be injected, θinj. The initial speed of the particles,
vinj = 0.1 aΩorb, is determined from cs at L1, where the
z-velocity is an arbitrary fraction of this. The actual value
is not critical as the gas will rapidly accelerate to become
supersonic, and the stream is not well-resolved in our sim-
ulations. For some simulations the third dimension is sup-
pressed. For the inner boundary condition a hole of radius
R1 centred on the primary was used, with R1 set to values
between 0.03 a and 0.1 a, and particles entering this were
removed from the simulation.
The calculations presented in Murray (1998) and Mur-
ray, Warner & Wickramasinghe (2000) were of cool isother-
mal discs. Here we model a more realistic steady-state disc
where cs is a function of disc radius, r, and is given by
cs = C r
−3/8. C is a constant and in general we have set
it equal to 0.05 aΩorb. This means that cs at the resonance
radius in each of our simulated discs will be≃ 0.067 aΩorb as
detailed in Table 2. The SUUMa systems ZCha and OYCar
in outburst have brightness temperatures, TBR, in the outer
disc of ∼ 6000 – 7000K (Horne & Cook 1985; Bruch, Beele
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 3
Table 1. Summary of simulations and results. The first column denotes the run number. Columns 2 to 6 describe the simulation
parameters, recording respectively whether the simulation is conducted in 2D or 3D, the mass ratio, the injection time step, the constant
describing the sound speed (cs = C r
−3/8) and the radius of the primary (the central hole). The time at which the simulation was
terminated is given in column 7. The remaining columns record outcomes of the simulations: the mean superhump excess as measured
from the simulated lightcurve, the final total number of particles in the simulation, the final average number of neighbours, the strengths
of the eccentric and 2-armed spiral modes averaged over the final 10 superhump periods in each simulation, and, as a measure of when
the disc initially encounters the resonance, the time at which S(1,0) = 0.01. ζ = 1.0 in all cases.
Run 2D/ q ∆t C R1 tend mean ǫ np nne S(1,0) S(2,2) tecc
3D (Ω
) (aΩorb) (a) (Ω
1 3D 0.3333 0.0025 0.050 0.03 3798.0575 — 123995 28.28 0.001 0.104 —
2 3D 0.2422 0.0025 0.050 0.03 8616.6650 — 141474 28.50 0.001 0.099 —
3 3D 0.2346 0.0025 0.050 0.03 15496.7600 1 0.062 143274 28.48 0.007 0.099 —
4 3D 0.2270 0.0025 0.050 0.03 14340.9425 1 0.058 143088 28.67 0.065 0.0948 6089.7625
5 3D 0.2195 0.0025 0.050 0.03 15487.8575 0.056 143730 28.18 0.092 0.092 4542.1675
6 3D 0.2121 0.0025 0.050 0.03 8738.4900 0.054 144674 28.62 0.107 0.090 3188.7950
7 3D 0.1765 0.0025 0.050 0.03 4590.9525 0.043 150628 27.41 0.151 0.081 1271.1125
8 3D 0.1429 0.0025 0.050 0.03 4173.9300 0.034 160459 27.63 0.161 0.078 777.2475
9 3D 0.1111 0.0025 0.050 0.03 3134.1350 0.024 171971 26.92 0.172 0.075 770.7275
10 3D 0.0811 0.0025 0.050 0.03 4323.1750 0.017 175989 24.76 0.227 0.062 477.1300
11 3D 0.0526 0.0025 0.050 0.03 3535.2575 0.015 2 227047 26.76 0.001 0.080 426.7250
12 3D 0.0256 0.0025 0.050 0.03 3491.8375 — 257057 24.07 0.001 0.070 —
13 3D 0.2400 0.0025 0.050 0.05 2000.0000 0.05 118947 21.87 0.001 —
14 3D 0.1900 0.0025 0.050 0.05 1204.2430 0.04 128135 21.56 0.003 —
15 3D 0.1500 0.0025 0.050 0.05 2000.0000 0.04 141277 22.09 0.074 1028.375
16 3D 0.1000 0.0025 0.050 0.05 2135.4880 0.03 145893 20.1 0.186 622.850
17 3D 0.0300 0.0025 0.050 0.10 942.0200 — 127310 27.06 0.001 —
18 3D 0.0700 0.0025 0.045 0.05 1014.2300 0.01 431782 98.03 0.008 —
19 2D 0.5385 0.0100 0.050 0.03 3081.36 — 18148 14.95 0.008 0.088 —
20 2D 0.4815 0.0100 0.050 0.03 21798.61 0.093 17973 14.41 0.23 0.057 1858.91
21 2D 0.4286 0.0100 0.050 0.03 18638.44 0.079 18293 14.35 0.30 0.043 974.19
22 2D 0.3333 0.0100 0.050 0.03 18285.39 0.055 18493 13.85 0.43 0.020 764.33
23 2D 0.2500 0.0100 0.050 0.03 18660.42 0.037 18822 14.89 0.53 0.019 591.20
24 2D 0.1765 0.0100 0.050 0.03 17355.41 0.026 20438 13.27 0.56 0.021 400.43
25 2D 0.1111 0.0100 0.050 0.03 14731.33 0.016 24766 13.97 0.54 0.020 308.67
26 2D 0.0526 0.0100 0.050 0.03 8872.88 — 41753 16.00 0.001 0.081 —
27 2D 0.0256 0.0100 0.050 0.03 1608.00 — 46093 15.47 0.001 0.077 —
1 equilibrium not yet reached, simulation continues to run 2 temporary appearance of superhumps
& Baptista 1996). For the eight SUUMa eclipsing systems
in Table 3 we obtain 0.031 aΩorb . cs . 0.036 aΩorb at
the disc midplane at R3:1, assuming a mass transfer rate
of Ṁ = 10−9 M⊙ yr
−1 and a fully ionised cosmic mixture
of gases. Hence our simulations have a sound speed within
a factor of two of that prevailing in reality, a much better
match than previously possible.
We have also improved on the viscosity used in previous
calculations. Our code includes an artificial viscosity term
which generates a shear viscosity in the disc,
ν(r) = κ ζ cs H, (1)
where ζ is the dimensionless artificial viscosity parameter.
ζ = 1 here, and H is the disc scaleheight. κ may be found
analytically, and for a standard cubic spline kernel in three
dimensions, κ = 1/10, whilst in two dimensions, κ = 1/8.
The bulk viscosity is fixed to be twice the shear viscosity.
Using the Shakura–Sunyaev viscosity parametrisation (ν =
α cs H) (Shakura & Sunyaev 1973), our simulated 3D and
2D discs have α(R3:1) = 0.1 and 0.125 respectively. This
matches estimates of α ∼ 0.1−0.2 derived from observation
of systems in the high viscosity state (Smak 1999).
2.2 Simulation Results
Table 1 summarises our simulations. In some cases these
have run for > 2000 orbits without reaching mass-transfer
equilibrium; one has run for over a year of elapsed time.
The rate of energy dissipation (i.e. viscously-generated lu-
minosity with units M a2 Ω3orb) from different regions in the
disc is recorded each timestep and used to produce simula-
tion lightcurves. Psh was determined from timings of maxima
in dissipation in a smoothed lightcurve for the disc region
r > 0.3 a. Column 8 of Table 1 gives the mean value of ǫ
so obtained, over the time in which the system has reached
equilibrium where applicable. Figure 1 shows the disc evo-
lution for the 3D simulations 1 to 12. Here the disc mass is
taken to be the number of SPH particles in the disc. The disc
eccentricity is estimated from the eccentric mode strength,
c© 2005 RAS, MNRAS 000, 1–18
4 A. J. Smith et al.
Table 2. Characterisation of the simulated accretions discs at the 3:1 resonance radius. The 3:1 resonance radius is recorded in the
fourth column, followed by the sound speed, the scale height, the characteristic value of the ratio of the disc semi-thickness to the radius
as given in Goodchild & Ogilvie (2006), the shear viscosity, the bulk viscosity, the Shakura–Sunyaev parameter for the shear viscosity
and that for the bulk viscosity.
Run 2D/ q R3:1 cs(R3:1) H(R3:1) h(R3:1) νsh(R3:1) νbk(R3:1) αsh(R3:1) αbk(R3:1)
3D (a) (aΩorb) (a) (a) (a
2Ωorb) (a
2Ωorb)
1 3D 0.3333 0.437 0.068 0.023 0.039 1.55× 10−4 3.10× 10−4 0.100 0.200
2 3D 0.2422 0.447 0.068 0.023 0.038 1.52× 10−4 3.05× 10−4 0.100 0.200
3 3D 0.2346 0.448 0.068 0.023 0.037 1.52× 10−4 3.04× 10−4 0.100 0.200
4 3D 0.2270 0.449 0.068 0.023 0.037 1.52× 10−4 3.04× 10−4 0.100 0.200
5 3D 0.2195 0.450 0.067 0.022 0.037 1.52× 10−4 3.03× 10−4 0.100 0.200
6 3D 0.2121 0.451 0.067 0.022 0.037 1.51× 10−4 3.03× 10−4 0.100 0.200
7 3D 0.1765 0.455 0.067 0.022 0.037 1.50× 10−4 3.01× 10−4 0.100 0.200
8 3D 0.1429 0.460 0.067 0.022 0.036 1.49× 10−4 2.98× 10−4 0.100 0.200
9 3D 0.1111 0.464 0.067 0.022 0.036 1.48× 10−4 2.96× 10−4 0.100 0.200
10 3D 0.0811 0.468 0.066 0.022 0.035 1.47× 10−4 2.94× 10−4 0.100 0.200
11 3D 0.0526 0.473 0.066 0.022 0.035 1.46× 10−4 2.92× 10−4 0.100 0.200
12 3D 0.0256 0.477 0.066 0.022 0.034 1.45× 10−4 2.91× 10−4 0.100 0.200
13 3D 0.2400 0.447 0.068 0.023 0.041 1.52× 10−4 3.05× 10−4 0.100 0.200
14 3D 0.1900 0.454 0.067 0.022 0.040 1.51× 10−4 3.02× 10−4 0.100 0.200
15 3D 0.1500 0.459 0.067 0.022 0.039 1.49× 10−4 2.99× 10−4 0.100 0.200
16 3D 0.1000 0.466 0.067 0.022 0.038 1.48× 10−4 2.96× 10−4 0.100 0.200
17 3D 0.0300 0.476 0.066 0.022 0.037 1.45× 10−4 2.91× 10−4 0.100 0.200
18 3D 0.0700 0.470 0.060 0.020 0.034 1.19× 10−4 2.38× 10−4 0.100 0.200
19 2D 0.5385 0.416 0.069 0.023 0.042 2.01× 10−4 4.02× 10−4 0.125 0.250
20 2D 0.4815 0.422 0.069 0.023 0.041 1.99× 10−4 3.98× 10−4 0.125 0.250
21 2D 0.4286 0.427 0.069 0.023 0.040 1.97× 10−4 3.94× 10−4 0.125 0.250
22 2D 0.3333 0.437 0.068 0.023 0.039 1.94× 10−4 3.88× 10−4 0.125 0.250
23 2D 0.2500 0.446 0.068 0.023 0.038 1.91× 10−4 3.82× 10−4 0.125 0.250
24 2D 0.1765 0.455 0.067 0.022 0.037 1.88× 10−4 3.76× 10−4 0.125 0.250
25 2D 0.1111 0.464 0.067 0.022 0.036 1.85× 10−4 3.70× 10−4 0.125 0.250
26 2D 0.0526 0.473 0.066 0.022 0.035 1.83× 10−4 3.66× 10−4 0.125 0.250
27 2D 0.0256 0.477 0.066 0.022 0.034 1.82× 10−4 3.63× 10−4 0.125 0.250
S(1,0), where the strength of the (k θ−lΩorb t) mode, S(k,l), is
obtained by Fourier decomposing the simulated disc density
distributions in azimuth and time (Lubow 1991b; Murray
1996).
2.2.1 The growth of the eccentricity
Figure 1 shows that initially, the eccentricity grows expo-
nentially. Lubow (1991a) found an exponential eccentricity
growth rate which is proportional to the square of the mass
ratio; in contrast we see a very low growth rate for high
mass ratios (q & 0.2) (Figure 2). Our results can be prof-
itably compared with those of Goodchild & Ogilvie (2006)
who formulated a single equation to describe the resonant
excitation, propagation and viscous damping of the eccen-
tricity in a 2D accretion disc. They showed that the res-
onance may have the effect of locally suppressing the ec-
centricity which consequently leads to extremely low eccen-
tricity growth rates. We find for q = 0.2195, for example,
a steady state is not reached until ∼ 1700 orbital periods,
and for q = 0.2346 a steady state is still not achieved after
∼ 2500 orbital periods.
Figure 2. Eccentricity growth rate (growth rate of the strength
of the (1, 0) eccentric mode) as a function of mass ratio for sim-
ulations presented in this work, and in previous works (Murray
1998). Simulation parameters are as indicated in the legend.
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 5
Table 3. Parameters of the eight eclipsing SUUMa systems: mass ratio, primary mass, total system mass (M1 + M2), observed disc
precession rate (from the measured ǫ value), dynamical precession rate as calculated from Equation 6 at the location of the 3:1 resonance,
inferred pressure contribution to precession (ωobs − ωdyn) and pitch angle of the spiral wave (see Section 4 for details). Errors on ωobs
and ip correspond to the range of superhump periods observed. In the last but one column the sound speed, in SPH units, is calculated
for the midplane of each disc at the 3:1 resonance radius, assuming a mass transfer rate of 10−9 M⊙ yr
−1 and a fully ionised ‘cosmic’
mixture. References for M1 and Mt are provided in the final column. Details and references for all other observational data can be found
in Table 5.
System q M1 Mt ωobs ωdyn(R3:1) ωpr ip cs(R3:1) Ref
(M⊙) (M⊙) (rad d
−1) (rad d−1) (rad d−1) (◦) aΩorb
OY Car 0.102 ± 0.003 0.685 ± 0.011 0.755± 0.011 2.189
+0.447
−0.181 3.649 −1.460 12.69
+2.44
−0.70 0.0334
XZ Eri 0.1098 ± 0.0017 0.767 ± 0.018 0.851± 0.018 2.697
+0.328
−0.0035 4.015 −1.319 13.02
+1.91
−0.16 0.0321
IY UMa 0.125 ± 0.008 0.79± 0.04 0.89± 0.04 2.208
+0.329
−0.120 3.714 −1.506 11.17
+1.42
−0.41 0.0322
Z Cha 0.1495 ± 0.0035 0.84± 0.09 0.965± 0.091 3.161+0.634
−0 4.282 −1.121 14.31
+6.85
−0 0.0360
HT Cas 0.15 ± 0.03 0.61± 0.04 0.70± 0.04 2.725
−0 4.343 −1.618 11.62
−0 0.0347
DV UMa 0.1506 ± 0.0009 1.041 ± 0.024 1.198± 0.024 2.349
+0.080
−0.286 3.738 −1.389 10.54
+0.31
−0.92 0.0314
OU Vir 0.175 ± 0.025 0.90± 0.19 1.06± 0.19 2.732+0.033
−0 4.987 −2.255 8.76
+0.06
−0 0.0305
V2051 Oph 0.19 ± 0.03 0.78± 0.06 0.93± 0.07 2.748
+0.815
−0.319 6.201 −3.453 7.83
+1.11
−0.33 0.0312
1 Wood et al. (1989) 2 Feline et al. (2004c) 3 Steeghs et al. (2003) 4 Wade & Horne (1988) 5 Horne, Wood & Stiening (1991)
6 Feline et al. (2004a,b) 7 Baptista et al. (1998)
2.2.2 The enhanced viscous torques of an eccentric disc
For all calculations in which the disc (eventually) becomes
sufficiently eccentric we see similar behaviour but on dif-
ferent timescales. As Figure 1 shows, initially the disc mass
builds, exponentially approaching a steady-state value. Then
as eccentricity increases, the disc mass tends to a new lower
steady state value. This reveals the non-eccentric accretion
disc approaching an equilibrium between the tidal removal
of angular momentum and the angular momentum added
by material from the L1 point. As the disc becomes eccen-
tric it readjusts to a new equilibrium in which tidal removal
of angular momentum is more efficient. This is exactly the
premise of the TTI model.
Table 4 quantitatively examines this. As a measure of
the increase in efficiency of tidal removal of angular momen-
tum in the eccentric disc, we took the decrease in disc mass
between maximum and the value it finally reached in equi-
librium (column 4 of Table 4). For a steady state accretion
disc the total mass is
Mtot =
Z Rout
3ν(R)
dR (2)
(Frank, King & Raine 1985). The mass transfer rate through
L1 remains constant, and we are comparing the equilibria
with and without an eccentric precessing disc. In both cases,
the mass accretion at all disc radii must equal the mass
transfer rate at L1. Thus, the change in disc mass can then
be related to a change in viscosity by
Mb −Ma = ∆M = A
f (R) dR
νb(R)
f (R) dR
νa(R)
, (3)
where A is a constant, and the subscripts b and a respec-
tively refer to before and after the disc became fully eccen-
tric. Defining
f(R)dR
νb(R)
f (R) dR
, (4)
where ν̄b is some unknown weighted mean value of νb(R),
and similarly for 1
, then we have
f (R) dR
. (5)
This quantity is given in column 5 of Table 4. The radial de-
pendence of H is cs(r/µ)
1/2r, where µ = 1/(q+1). Together
with the sound speed prescription that we use, then this
allows us, via Equation 1, to give explicitly the fractional
change in viscous torque necessary to bring about the de-
crease in disc mass we see (column 6 of Table 4). We assume
Rout = Rtides. Here, we use the formulation Rtides ≃ 0.9R1
(Frank, King & Raine 1985), where R1 is the effective Roche
lobe radius (Eggleton 1983), though we note the limitation
of this approximation (e.g. Murray, Warner & Wickramas-
inghe 2000; Truss 2007). For runs 6, 7, 9 and 10 we measured
the average outer disc radius, averaging the position of the
outermost particle as a function of azimuth over a super-
hump period. We found in each case the difference between
this value and Rtides to be < 0.006 a.
2.2.3 The superhump
The development of the superhump for two selected simula-
tions is shown in Figure 3. For q = 0.0526, superhumps were
only temporarily present at a time when the eccentricity was
highest. For q = 0.2121 we see that the superhump profile
c© 2005 RAS, MNRAS 000, 1–18
6 A. J. Smith et al.
Figure 1. Disc mass (red) and disc eccentricity (green) as a function of time for the 3D simulations 1 to 12. Every 5000th and 2000th
timestep is plotted respectively. Also plotted is the superhump period (blue) in the form of O-C. This is calculated relative to the mean
superhump period from each run, derived from the point at which the superhump signal becomes well-formed. The dissipation (smoothed)
from the disc for r > 0.3 a is shown in grey. The mass ratio is shown in the upper left of each panel.
evolves with time as the disc is reaching eccentric equilib-
rium. In the case of other more extreme mass ratios (smaller
values of q) there are higher frequency periodic components
present in the early development of the superhumps (Smith,
PhD thesis in prep.). Our discs accumulate from zero mass
with no switch between viscosity states, so the development
of their superhump will differ from that of discs observed in
superoutburst.
Variations in the superhump period are displayed in the
form of O-C (‘observed’ minus ‘calculated’) in Figure 1. In-
tervals of non-zero 2nd derivative in the value of O-C indi-
cate period changes. Generally Psh decreases as the system
approaches eccentric steady-state. For q = 0.0526 a steady-
state eccentric disc is not achieved, and the period evolution
differs.
Each simulated lightcurve was folded on the mean su-
perhump period once equilibrium was reached (Figure 4).
These are asymmetric with, for most, a steep rise and slower
decline. A secondary hump structure is seen, the profile dif-
ferent for different mass ratios. In general, these superhump
profiles resemble those observed in CVs cf. figure 8 of Patter-
son et al. (1995), figure 4 of Imada et al. (2006) and figure 6
of Maehara, Hachisu & Nakajima (2006).
In Figure 5 we show the density distributions of the
accretion disc for run 9 at selected superhump phases, and
compare these with the density profile of a disc which does
not show superhumps. The disc is not maximally distorted
at the time of maximum energy production (φsh = 0) as
might be expected, but at φsh = 0.71. For q = 0.1111, this
coincides with a small secondary peak seen in Figure 4. The
superhumping disc is most similar to the non-superhumping
disc at superhump maximum, the most visible difference be-
ing in the strongly enhanced spiral density waves in the su-
perhumping disc. Furthermore, the appearance of the spiral
density wave changes dramatically over Psh. This illustrates
the crucial importance of the spiral density waves to the
superhump phenomenon, as suggested by Lubow (1991a)
and Osaki (2003) who considered analytic theory and obser-
vations respectively. Comparing superhump maximum with
superhump minimum, then we see more a more open spiral
structure at superhump maximum with outer reaches show-
ing enhanced density. After superhump maximum, differen-
tial motion in the superhumping disc causes the spiral to
become eccentric and more tightly wrapped with lower den-
sity contrast (see the bottom two panels in Figure 5). Each
radius in the disc has its own characteristic Keplerian angu-
lar velocity, and eccentric mode precession rate. Differential
precession in the eccentric fluid disc cannot persist, because
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 7
Table 4. Enhanced viscous torque resulting from the accretion
disc encounter with the eccentric resonance. The fourth column
records the fractional decrease in disc mass before and after the
disc becomes fully eccentric. The fifth column gives a measure of
the change in viscous torque indicated by this mass decrease as
described in the text and given by Equation 5. The final column
assumes the radial dependence of viscosity is given by Equation 1,
and records, in that case, the fractional change in viscous torque
necessary to bring about the decrease in disc mass.
Run 2D/ q ∆M/M 1
νa/νb
3D 1/(a2Ωorb)
5 3D 0.2195 0.0223 175 1.019
6 3D 0.2121 0.0270 211 1.023
7 3D 0.1765 0.0464 360 1.040
8 3D 0.1429 0.0360 273 1.031
9 3D 0.1111 0.0465 353 1.040
10 3D 0.0811 0.1063 800 1.098
20 2D 0.4815 0.0735 423 1.060
21 2D 0.4286 0.1145 665 1.097
22 2D 0.3333 0.205 1200 1.190
23 2D 0.2500 0.267 1538 1.258
24 2D 0.1765 0.257 1392 1.231
25 2D 0.1111 0.165 823 1.127
Figure 3. Evolution of the simulated lightcurve for runs 6 and
11. Each panel covers 5 orbits. The lightcurve, which has been
smoothed somewhat, is given by the rate of energy dissipation
for radii r > 0.3 a. The mass ratio is shown in the uppermost
panel of each.
this would cause widespread orbit crossings. This converg-
ing fluid motion causes the strong spiral shock shown in the
upper right panel of Figure 5, and it is the dissipation associ-
ated with this shock which powers the observed superhump.
Figure 5 illustrates the mechanism by which spiral waves
power superhumps.
Figure 4. Equilibrium lightcurves (energy dissipation rate for
radii r > 0.3 a) folded on the derived mean superhump period for
runs 5 to 10. The superhump cycle is repeated for clarity.
Figure 5. Density profiles comparing the non-superhumping disc
of run 2 (upper left panel) with the superhumping disc of run 9
at 3 selected superhump phases, φsh. The coordinates are centred
on the binary system centre of mass. The solid line is the primary
Roche lobe and the dashed line is the 3:1 resonance radius. The
upper left panel also shows our definition of azimuth used in Sec-
tion 3.1. The secondary star is at an azimuth of π radians with
respect to the primary, the position of which is marked with a
cross.
c© 2005 RAS, MNRAS 000, 1–18
8 A. J. Smith et al.
Figure 6. Comparison of results for 2D and 3D simulations with q = 0.1111, runs 24 and 9 respectively. The following are shown above
and below respectively for each: the disc evolution, details of the simulated lightcurves, folded superhump lightcurves and density profiles.
Here symbols are as in Figure 5, and the colour scale used is the same for both the 2D and 3D discs. The simulated lightcurves cover 5
orbital cycles once the disc mass and eccentricity had stabilised. The raw lightcurve is shown in the top panel and increasingly smoothed
lightcurves below.
2.2.4 2D versus 3D
In Figure 6 we compare the results of q = 0.1111 calcu-
lations in 2D and 3D. In the 2D simulations the accretion
disc achieves a far higher eccentricity more quickly. This
in turn affects the superhump profile. We see from Table 1
that these 2D accretion discs (runs 19 to 27) can become
eccentric at much larger values of mass ratio than in the 3D
simulations, and well beyond the bounds that theory sug-
gests, indicating that three dimensions and the high mass
resolution (∆ t = 0.0025) is necessary to accurately simu-
late the processes involved. In this case, (3D, ∆ t = 0.0025)
no disc develops superhumps where q & 0.24. We also see
that for the same value of q, the 2D discs precess at smaller
rates than the 3D discs; compare for example run 24 with
run 7.
3 ECCENTRICITY DISTRIBUTION IN
ACCRETION DISCS
An eccentric accretion disc underlies the superhump phe-
nomenon, but how eccentric does the disc become? How does
the eccentricity vary throughout the disc, and how does it
change with superhump phase?
We plotted an estimate of the disc eccentricity in our
simulations in Figure 1. Calculating the strength of the
(1,0) mode takes into consideration the disc as a whole.
For those 3D discs which reached equilibrium and showed
superhumps, final values of eccentricity e between ∼ 0.09
and ∼ 0.23 were seen, with more extreme mass ratios har-
bouring more eccentric discs. This agrees with modelling
of line profiles in AMCVn assuming a constant eccentric-
ity throughout the disc: Patterson, Halpern & Shambrook
(1993c) found e = 0.1 − 0.2. In their analytic study Good-
child & Ogilvie (2006) examined the spatial distribution of
eccentricity, but explicitly avoided examination of the be-
haviour on the orbital timescale. Their eccentricity distri-
bution was locally suppressed by the presence of dynamical
resonances. Next we examine the eccentricity distributions
of our simulated discs, using two complementary methods
to characterise their spatial and temporal variation.
3.1 Eccentricity distribution from particle
trajectories
An instantaneous snapshot of the radial eccentricity distri-
bution was found by projecting the elliptical orbit of each
particle in the disc using its position and velocity. Particles
in the mass transfer stream were discounted. Each parti-
cle was assigned a radius given by an average along this
elliptical path, weighted by the time spent at each radius.
This method calculates the trajectories particles would have
if they orbited an isolated primary star. This simplification
may introduce certain artefacts into the calculated eccentric-
ity distributions, but does give an adequate approximation.
Figure 7 shows the results for each of the 3D simula-
tions 1 to 12, at superhump maximum where applicable. In
this and subsequent similar figures, each point represents
the mean eccentricity of particles, obtained in the manner
described above, in each of a number of azimuthal and ra-
dial bins of size 0.1π rad and 0.1 a respectively. Overplot-
ted on these figures are the 3:1, 4:1 and 5:1 resonance radii
from right to left respectively. For simulations 11 and 12 the
2:1 resonance is also shown. Each particle is colour coded
according to the azimuthal position of the particles it rep-
resents, as shown in the key, where θ is the angle defined
in Figure 5. Purple and blue particles are approaching the
mass donor star, while green and yellow particles are re-
ceding from it. The eccentricity distributions at different
azimuths are clearly distinct. Caution is required, however,
in interpreting the radii assigned in this section. The parti-
cles at r > R3:1 are concentrated at azimuths close to the
x-axis, where in fact the disc edge is relatively close to the
compact object: these particles have larger radii assigned to
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 9
Figure 7. Instantaneous eccentricity distribution plotted as a function of radius for 3D runs 1 to 12, once equilibrium is reached (see
text for details). Each point represents the average of all disc particles in each of 20 equal azimuthal sections and binned in radius. These
are colour coded by their azimuthal position in the disc, as defined in Figure 5, according to the scale on the right. The time is shown
in the upper right of each panel, together with the superhump phase where applicable.
them than those found further out at other azimuths be-
cause they have relatively circular projected orbits.
For all the discs with q > 0.0811 the eccentricity in the
outer disc has a wide spread around e ∼ 0.2. This agrees
reasonably with the eccentricities found in the outer discs of
OYCar: Hessman et al. (1992) found e = 0.38±0.10; and IY
UMa: Patterson et al. (2000a) estimate e = 0.29± 0.06. For
0.2121 6 q 6 0.3333 the eccentricity distributions are low
for the disc within the 5:1 resonance radius, then increase
at larger radii to e up to ∼ 0.3, this apparent eccentric-
ity likely due to tidal distortions in the outer disc. For this
range of q, if every particle were individually plotted, the ec-
centricity distribution outside R5:1 is “double horned”, with
a concentration of points near the upper and lower edges
of the distribution and a dearth of points midway between
the envelope. These two “horns” correspond to the halves of
the disc which are approaching and receding from the mass
donor star.
Figure 7 demonstrates that the body of the disc within
the 5:1 resonance appears eccentric for 0.081 6 q 6 0.2195.
For these discs the entire mass appears to partake in the ec-
centric instability. This roughly flat region of the eccentricity
distribution has an eccentricity dependent on q, increasing
to ∼ 0.12 as we move to smaller q.
Figure 8 shows how the eccentricity distribution devel-
ops in run 9 as the disc becomes eccentric. Before super-
humps set in (top left panel) the eccentricity distribution
resembles those seen in the less extreme mass ratio systems
i.e. the top row of panels in Figure 7. The remaining 5 pan-
els of Figure 8 span the interval where the disc eccentricity
is growing and the disc is emptying as it approaches the
new mass transfer equilibrium caused by the enhanced vis-
Figure 8. Development of accretion disc radial eccentricity dis-
tribution for run 9. Azimuthal colour-coding is as in Figure 7.
cous torque (c.f. Figure 1). In each case superhump max-
imum is plotted. Within R5:1 the eccentricity of the body
of the disc grows, while the range of eccentricities at any
given radius remains more or less unchanged. Outside R5:1
the eccentricity distribution becomes more scattered, as the
range of eccentricities at any given radius increases. The last
panel of Figure 8 shows the distribution when the strength
of the (1,0) mode approaches maximum, but mass equilib-
rium is not yet reached. It resembles the q = 0.1111 panel
c© 2005 RAS, MNRAS 000, 1–18
10 A. J. Smith et al.
Figure 9. As Figure 7, showing the changing radial eccentricity
distribution with superhump phase, for run 9.
in Figure 7, which occurs almost 230 orbits later, after mass
equilibrium is established.
The eccentricity distribution is shown as a function of
superhump phase in Figure 9. At some phases and radii
there is little spread in e, while at other phases there is a
large spread in e at the same radius. This behaviour may
be related to the cusps at the resonance points discussed by
Goodchild & Ogilvie (2006), but their treatment explicitly
excluded the behaviour of the disc on timescales comparable
to or shorter than the orbital period. Numerical simulation
facilitates examination of the flexing of the disc during a
single orbit, which is important as it is this flexing which
gives rise to the modulation in viscously-generated or repro-
cessed light which constitutes an observed superhump. The
flexing of the disc over the superhump period appears most
dramatic at radii between R5:1 and R3:1, in accordance with
Pearson’s (2006) point that the dynamic precession rate at
the 4:1 resonance actually agrees far better with the obser-
vations.
Figure 9 neatly illustrates the changes which occur as
the binary’s gravitational potential moves relative to the
disc. Even as far in as r=0.15a we see the more eccentric
particles belonging to the approaching side of the disc at
superhump maximum and the receding part of the disc at
superhump minimum.
3.2 Eccentricity from the mass distribution
How good a representation of the disc eccentricity are the
distributions we calculated in section 3.1? If we were to fol-
low a single particle as it orbits in the accretion disc it would
not move on the elliptical orbit we extrapolated from its in-
stantaneous velocity and position. We now use an alterna-
tive way of examining the eccentricity of the accretion disc,
looking at the mass distribution.
The disc in run 9 was split into 100 azimuthal sections.
For each, we recorded the number of particles contained
within each radial step outwards. In Figure 10, we show
the contour maps which result. Every step of 24 particles
was recorded and colour-coded, so each colour represents a
‘contour’ of enclosed mass. To each contour we fitted an el-
lipse which has one focus at the WD. The fitted parameters
are the semi-major axis, asemi, the eccentricity, e, and the
Figure 11. As Figure 10 for the non-superhumping disc of run
angle the semi-major axis makes with the positive x-axis in
an anti-clockwise direction, α0, measured in radians. This is
analogous to our definition of θ in Section 3.1. For the out-
ermost complete contour, that is the outermost contour for
which each azimuthal section is represented, marked in dark
blue, these fitted parameters are displayed in the bottom
right of each panel in Figure 10 and apply to the overplot-
ted magenta ellipse.
Figure 10 allows us to look at the disc mass distribution
in a more quantitative way than simply looking at the den-
sity distribution. The eye is drawn to the pronounced non-
axisymmetry outside the 3:1 resonance, though this mass
constitutes less than 9 per cent of the mass in the disc, and
is the lowest density region. This mass contributes only 4
per cent of the total dissipation, and this contribution to
the modulation on Psh is not in phase with the overall su-
perhump. The region of the disc within R3:1 is overwhelm-
ingly responsible for generating the dissipation-powered su-
perhump in the simulations.
The four panels in Figure 10 are equally spaced in su-
perhump phase. It is noticeable that the outermost disc does
not simply precess as the binary frame moves. Instead the
flexure of the disc combines with the orbital motion of the
binary frame to leave the outer edge of the disc almost fixed
in the binary frame between φsh = 0.02 and φsh = 0.27;
similarly the outer edge of the disc remains almost the same
in the two lower panels at φsh = 0.52 and φsh = 0.77. Fig-
ure 11 shows that a non-superhumping disc is extended at
similar azimuths. These extended disc edges are analogous
to the raised tides in Earth’s oceans. As Figure 11 shows,
this effect produces an elliptical shape centred on the pri-
mary. This illustrates the distinction that needs to be made
between different contributions to the eccentricity distribu-
tions found here and in Section 3.1. In the outer disc tidal
distortions are important, and are present in discs at all
mass ratios. This is to be distinguished from the (1,0) ec-
centric mode eccentricity found in the superhumping discs
which has the primary at one focus. Truss (2007) finds that
for discs which have not yet become tidally unstable, the
exact azimuth of the extended ‘wing’ depends on the mass
ratio, viscosity parameter and sound speed of the gas (the
major axis of the outer streamline moves clockwise with de-
creasing q, or increasing α and cs). The interior regions of
the disc (e.g. around R4:1) more closely approximate simple
relative motion between a slowly apsidally precessing disc
orientation and a rapidly moving orbital frame.
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 11
Figure 10. Pictorial representation of the mass distribution in the disc for run 9 at four superhump phases. The disc is split into
100 equal azimuthal sections (dashed lines). Each colour graduation represents ‘contours’ of particle numbers. Moving outwards, each
contour represents an increase of 24 particles in each azimuthal section (90 contours in all). Further details are described in the text.
The outermost dark red points mark the positions of the outermost particles in each azimuthal section. Also plotted are three circles at
the 3:1, 4:1 and 5:1 resonance radii moving inward respectively (black, dark and light grey) which also act as a guide to the eye to show
the non-circularity of the contours.
Figure 12 shows the fitted ellipse parameters for the
four phases shown in Figure 10. The average radius of all
points making up the contour was used. The general trends
in Figs. 12 and 10 are seen in analogous plots for q = 0.1765,
including the peaks and troughs in the radial distribution
of the eccentricity parameter. Within R3:1, the orientation
of the semi-major axis, α0, changes systematically over the
superhump period and in the opposite sense to the motion
of the gas in the disc i.e. precession of the slowly-moving disc
as viewed from the rapidly rotating binary reference frame.
Figure 12 is complementary to Figure 9. Figure 9 shows
the eccentricity distributions deduced from instantaneous
velocities, while Figure 12 shows the eccentricity distribution
deduced from the instantaneous mass distribution. Since the
disc is continuously flexing in a complex way, these are not
the same. The disc motions can be described as superpo-
c© 2005 RAS, MNRAS 000, 1–18
12 A. J. Smith et al.
Figure 12. Fitted ellipse parameters for each ‘contour’ in Fig-
ure 10, as a function of radius. See Figure 10 and text for details.
Results for four superhump phases are presented corresponding
to the four panels in Figure 10. Vertical lines indicate the 3:1, 4:1
and 5:1 resonance radii from right to left respectively. Red points
indicate contours which are not complete.
sitions of the S(k,l) modes, and the resonance radii act as
nodes and antinodes in the complex standing wave dynamics
the disc executes over a full precession period. To summarise
contributions to these disc motions we compare, in table 1,
the strengths of the (1, 0) and (2, 2) modes.
For comparison, we show the equivalent of Figures 10
and 12 for q = 0.3333 in Figures 11 and 13, a disc which does
not show superhumps. Here, as we would expect, the mass
distribution remains approximately constant over time, and
the eccentricity is much lower.
4 PERIOD EXCESS VERSUS MASS RATIO:
DRAWING TOGETHER OBSERVATION,
THEORY AND SIMULATION
If a reliable relationship between ǫ and q can be deduced,
this would be immensely useful. ǫ can be easily measured
using relatively modest equipment, while the mass ratio q
is more fundamental and less easily determined. Patterson
Figure 13. As Figure 12 and corresponding to the panels in
Figure 11.
et al. (2005) fitted observations of eclipsing systems with
ǫ = 0.18 q + 0.29 q2. Figure 14 collates observation, theory
and simulations of positive superhumps. All observed sys-
tems with q determined by some means independent of ǫ
are plotted. Errors in q are formal errors given by the au-
thors and do not necessarily reflect the uncertainty in the
method. The eclipsing systems (red circles) should therefore
be given more weight. We note, however, the scatter of the
eclipsing systems seems typical of the scatter of the other
points. Error bars for ǫ denote the range of values observed
rather than errors in individual values, except in cases where
only one measurement has been made. Observational data
is tabulated in Table 5.
Our high-resolution 3D simulations, which are repre-
sented by blue squares in Figure 14, provide a far better
match with observed systems than previous studies by Mur-
ray (1998, 2000). As in Murray (2000) we see that simple
dynamical precession as given by
ωdyn = p(r)
1 + q
Ωorb, (6)
poorly represents observed systems. This is true even if the
location of the resonance for a gaseous disc, rather than
for isolated particles, is used; the location of the resonance
changes only by . 1 per cent. In a real gaseous disc the
retrograde effect that pressure forces have on disc preces-
sion rates must be taken into account (Lubow 1992; Murray
2000). In a gaseous disc, the excited eccentricity propagates
through the disc as a wave and is wrapped into a spiral by
the differential precession of the gas. Murray (2000) assumed
that the hydrodynamical precession is given by
ω = ωdyn + ωpr, (7)
where ωpr is the pressure contribution to the precession, and
showed that, under the assumption that the eccentricity is
tightly wound (that is if it is wound up on a length-scale
much smaller that the disc radius) then this pressure con-
tribution at the 3:1 resonance radius is
ωpr ≃ −
Ωorb a
tan ip
, (8)
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 13
Table 5. Observed superhump systems with independently determined mass ratio. All periods are given in days. Columns 5 and 6 list
the minimum and maximum observed superhump periods respectively. Column 7 gives the superhump excess where the errors indicate
the range of values as calculated from min Psh and max Psh. In column 8, the object type, (E) denotes a WD eclipsing system, while (e)
denotes a system which shows eclipse of the accretion stream/disc impact region only.
System q Porb Psh min Psh max Psh ǫ Type Ref
OY Car 0.102(3) 0.0631209180(2) 0.06454(2) 0.064245 0.06466 0.0225
+0.0019
−0.0047
SU UMa(E) 1,2,3,4,5
XZ Eri 0.1098(17) 0.061159491(5) 0.062808(17) 0.062603 0.06283 0.0270
+0.0003
−0.0034
SU UMa(E) 6,6,7,7,8
IY UMa 0.125(8) 0.07390897(5) 0.07588(1) 0.07558 0.07599 0.0267
+0.0015
−0.0041
SU UMa(E) 9,9,10,11,11
Z Cha 0.1495(35) 0.074499 0.07740 0.0768 0.0389
−0.0080
SU UMa(E) 12,13,14,14
HT Cas 0.15(3) 0.07364720309(7) 0.076077 0.0330 SU UMa(E) 15,16,17
DV UMa 0.1506(9) 0.0858526521(14) 0.08870(8) 0.0886 0.08906 0.0332
+0.0042
−0.0012
SU UMa(E) 6,6,18,18,18
OU Vir 0.175(25) 0.072706113(5) 0.07508(9) 0.07505 0.0327
−0.0005
SU UMa(E) 19,19,8,20
V2051 Oph 0.19(3) 0.0624278634(3) 0.06418(16) 0.06365 0.06439 0.0281
+0.0033
−0.0085
SU UMa(E) 21,22,23,23,23
WZ Sge 0.060(7) 0.0566878460(3) 0.05726(1) 0.05716 0.05738 0.0101
+0.0021
−0.0018
WZ Sge(e) 24,25,26,5,27
VY Aqr 0.11(2) 0.06309(4) 0.06437(9) 0.0642 0.06489 0.0203
+0.0082
−0.0027
SU UMa 28,29,5,5,5
CU Vel 0.115(5) 0.0785(2) 0.08085(3) 0.0799 0.0299
−0.0121
SU UMa 30,30,20,31
SW UMa 0.14(4) 0.056815(1) 0.058182(7) 0.05790 0.05833 0.0241
+0.0026
−0.0050
SU UMa 32,33,34,35,36
HS 2219+1824 0.19(1) 0.0599 0.06184 0.0324 SU UMa 37,37,37
EK TrA 0.20(3) 0.06288(5) 0.06492(10) 0.0648 0.0324
−0.0019
SU UMa 38,38,38,39
EI Psc 0.21(2) 0.044572(2) 0.04654 0.04579 0.0442
−0.0169
SU UMa 40,41,42,42
VW Hyi 0.21
+0.03
−0.02
0.074271038(14) 0.07714(5) 0.07621 0.07824 0.0386
+0.0148
−0.0125
SU UMa 43,44,44,45,45
YZ Cnc 0.22 0.0868(2) 0.09204 0.0905 0.0604
−0.0178
SU UMa 46,46,47,47
WX Hyi 0.23
+0.07
−0.04
0.0748134(2) 0.07737 0.0783 0.0342
+0.0124
SU UMa 43,48,49,50
T Leo 0.71(15) 0.0588190(5) 0.06022(2) 0.06021 0.06025 0.0238
+0.0005
−0.0002
SU UMa IP? 32,51,52,53,52
U Gem 0.357(7) 0.1769061898(30) 0.20 0.197 0.203 0.131
+0.017
−0.017
U Gem(e) 54,55,56,56,56
V603 Aql 0.22(3) 0.13809(12) 0.14640(6) 0.144854 0.14686 0.0602
+0.0033
−0.0112
Fast nova 57,58, ,59, 60
UU Aqr 0.30(7) 0.163580429(5) 0.17510(18) 0.0704 NL(E) 61,61,8
DW UMa 0.39(12) 0.136606527(3) 0.1454(1) 0.1461 0.0644
+0.0051
NL(E) 62,63,64,63
MV Lyr 0.43
+0.19
−0.13
0.132335 0.1377(4) 0.1487 0.0405
+0.0832
NL 65,66,67,68
AM CVn 0.18(1) 0.011906623(3) 0.012167 0.012164 0.012169 0.0218
+0.0002
−0.0002
AM CVn 69,70,70,70,70
KV UMa 0.037(7) 0.1699339(2) 0.170529(6) 0.17049 0.17073 0.0035
+0.0012
−0.0002
BHXRT 71,72,73,74,8
(XTE J1118+480)
QZ Vul 0.042(12) 0.3440915(9) 0.3469(1) 0.3474 0.0082
+0.0014
BHXRT 75,75,76,77
(GS 2000+2)
V1487 Aqr 0.058(33) 30.8(2) 31.4 31.2 31.6 0.0195
+0.0065
−0.0065
BHXRT 78,79,79,79,79
(GRS 1915+105)
V518 Per 0.111
+0.027
−0.033
0.2121600(2) 0.2157(10) 0.0167 BHXRT 80,80,81
(GRO J0422+32)
GU Mus 0.13(2) 0.432602(1) 0.4376(10) 0.0116 BHXRT 82,83,77
(N Mus 1991)
1 Wood et al. (1989) 2 Pratt et al. (1999) 3 Bruch et al. (1996) 4 Schoembs (1986) 5 Patterson et al. (1993a) 6 Feline et al.
(2004c) 7 Uemura et al. (2004) 8 Patterson et al. (2005) 9 Steeghs et al. (2003) 10 Uemura et al. (2000) 11 Patterson et al.
(2000a) 12 Wood et al. (1986) 13 Baptista et al. (2002) 14 Warner & O’Donoghue (1988) 15 Horne et al. (1991) 16 Feline et al.
(2005) 17 Zhang, Robinson & Nather (1986) 18 Patterson et al. (2000b) 19 Feline et al. (2004a) 20 Kato et al. (2003) 21 Baptista
et al. (1998) 22 Baptista et al. (2003) 23 Patterson et al. (2003) 24 Skidmore et al. (2002) 25 Patterson et al. (1998) 26 Ishioka
et al. (2002) 27 Patterson et al. (2002a) 28 Augusteijn (1994) 29 Thorstensen & Taylor (1997) 30 Mennickent & Diaz (1996) 31
Vogt (1981) 32 Shafter (1983) 33 Howell & Szkody (1988) 34 Nogami et al. (1998) 35 Semeniuk et al. (1997) 36 Robinson et al.
(1987) 37 Rodŕıguez-Gil et al. (2005) 38 Mennickent & Arenas (1998) 39 Vogt & Semeniuk (1980) 40 Thorstensen et al. (2002)
41 Uemura et al. (2002b) 42 Skillman et al. (2002) 43 Smith, Haswell & Hynes (2006) 44 van Amerongen et al. (1987) 45 Vogt
(1983) 46 Shafter & Hessman (1988) 47 Patterson (1979) 48 Schoembs & Vogt (1981) 49 Bailey (1979) 50 Walker, Marino &
Freeth (1976) 51 Shafter & Szkody (1984) 52 Kato (1997) 53 Lemm et al. (1993) 54 Naylor, Allan & Long (2005) 55 Smak (1993)
56 Smak & Waagen (2004) 57 Arenas et al. (2000) 58 Patterson et al. (1993b) 59 Haefner & Metz (1985) 60 Patterson et al. (1997)
61 Baptista, Steiner & Cieslinski (1994) 62 Araujo-Betancor et al. (2003) 63 Stanishev et al. (2004) 64 Patterson et al. (2002b)
65 Hoard et al. (2004) 66 Ritter & Kolb (2003) 67 Skillman, Patterson & Thorstensen (1995) 68 Pavlenko & Shugarov (1999) 69
Roelofs et al. (2006) 70 Skillman et al. (1999) 71 Orosz (2001) 72 Torres et al. (2004) 73 Uemura et al. (2002a) 74 Zurita et al.
(2002) 75 Harlaftis, Horne & Filippenko (1996) 76 Charles et al. (1991) 77 O’Donoghue & Charles (1996) 78 Harlaftis & Greiner
(2004) 79 Neil et al. (2006) 80 Webb et al. (2000) 81 Kato, Mineshige & Hirata (1995) 82 Orosz et al. (1996) 83 Casares et al.
(1997)
c© 2005 RAS, MNRAS 000, 1–18
14 A. J. Smith et al.
Figure 14. Superhump period excess plotted as a function of binary mass ratio for both observed systems and for SPH simulation. Also
plotted are the dynamical and hydrodynamical theoretical predictions. Simulations presented in this paper are displayed as filled squares,
as detailed in the top right legend. The second legend down refers to previous works by Murray. The third legend refers to theoretical
predictions, and the legend for observational data is in the bottom right. The inset shows data over a large range of q, whilst the main
panel focuses on the data at low q where most of the points are clustered, and is plotted on a logarithmic x-axis.
where ip is the pitch angle of the spiral wave, and cs is the
sound speed.
For each eclipsing SUUMa system we have calculated
an inferred pressure contribution to the precession rate, ωpr
(column 7 in Table 3) in the manner of Murray (2000). A
weighted mean gives ωpr = −1.36 rad d−1. This is combined
with the mean mass of the WD in systems below the period
gap (Smith & Dhillon 1998) to calculate predicted hydro-
dynamical precession rates assuming that the precessional
pressure contribution is similar for all systems.This is plot-
ted as a solid curve on Figure 14. We note that Pearson
(2006) considers the inclusion of pressure effects in an al-
ternative way. Whilst a reasonable fit to the observations is
achieved, the fact that there is a distribution in the eclipsing
systems above and below this curve, demonstrates that the
situation is not so simple. Possible contributory factors are
distributions in primary mass or in disc temperature. We
have plotted two further curves on Figure 14, one represent-
ing the hydrodynamical prediction for larger WD masses,
and the other encompassing both larger M1 and higher ac-
cretion disc temperature. These are the mean value of M1
for the eight SUUMa eclipsing systems and the value of ωpr
found by Murray (2000) respectively. We tabulate M1, the
total system mass, Mt, and the mid-plane sound speed at
the resonance radius for each of the eclipsing SUUMa sys-
tems (Table 3). These numbers do not seem to provide any
clue as to the true cause of the spread in observed systems.
In calculating cs, we assumed the same mass transfer rate
throughout (10−9 M⊙ yr
−1). Undoubtedly this, and conse-
quently the temperature, varies from system to system.
Figure 14 also shows the best fit that Goodchild &
Ogilvie (2006) found to their analytic curve. They found
that retrograde pressure terms are in fact only ∼ 1 per cent
of the dynamical term, far smaller than the values we infer
above. They suggest that the offset of observation from dy-
namical precession rate be due instead to averaging over the
disc, the eccentricity being distributed throughout the disc
rather than being sharply peaked at the resonance itself.
This is in accordance with our findings that dissipation-
powered superhump overwhelmingly originates in the disc
regions within R3:1. Goodchild & Ogilvie (2006) found that
their eccentricity distributions peaked close to 0.37 a, a value
put forward by Patterson (2001) to match observation. We
have plotted the dynamical curve as evaluated at 0.37 a (dot-
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 15
dashed light-blue line). We have also plotted two further
lines which show dynamical curves evaluated at some frac-
tion of the resonance radius, one to give the best weighted
fit to the eclipsing SUUMa systems (red) and the second
to give the best fit to our 3D simulations that have reached
equilibrium, runs 5 to 10, (blue). We find values of 0.83R3:1
and 0.87R3:1 respectively.
5 DISCUSSION
The precession rates for the simulated discs provide a much
closer match with observation than has been achieved previ-
ously. There are a number of reasons for this. The improved
mass resolution changes the stream–disc impact, leading to
a different angular momentum distribution in the disc and a
different precession rate. The radial sound speed distribution
is correct for a steady-state disc (as opposed to an isother-
mal one as in the case of Murray (1998)), which means
that the variation of density with radius is more realistic,
which in turn will determine whether the disc precesses and
at what rate. These updated simulations have hotter discs.
This means the retrograde effect of pressure on precession
will be greater (Lubow 1991a). The viscosity is more in line
with what is inferred for an α-disc in the high state. The
simulations of Murray (1998) were too viscous, which al-
lowed the disc to grow and penetrate the resonance too eas-
ily. The propagation of the eccentricity inward through the
disc would also have been inhibited. If the precession rate
is dictated by a weighted average of the eccentricity as sug-
gested by Goodchild & Ogilvie (2006), then the inhibition of
inward eccentricity propagation would lead to high preces-
sion rates. Finally, the extension to 3D changes the character
of the resonance. We can see in Figure 6 that the accretion
disc look very different in 2D and 3D. One reason for this is
that the character of the stream-disc interaction differs: for
2D mass injection all the particles follow one another exactly
and so the stream tends to punch through the outer disc and
deposit its angular momentum at r < R3:1. Conversely in
the 3D case the particles are more easily captured by the
outer disc and so they effectively reduce the specific angu-
lar momentum of the outer disc. The growth and decline of
the eccentric mode is strongly influenced by the interaction
with the mass transfer stream. The picture is further com-
plicated by the results of Kunze, Speith & Hessman (2001)
who found substantial stream overflow in their simulations.
An important extension of this work would be to system-
atically isolate the importance of each of the effects which
contribute to the disc precession rate.
It appears that the offset of both observation and simu-
lation from dynamical expectations of precession rates is due
to the averaging over the radii interior to R3:1 which partici-
pate in the disc precession, and the inclusion of a retrograde
effect of pressure forces is insufficient (Goodchild & Ogilvie
2006). In Section 3 and Figure 5 we show that the eccen-
tric instability is indeed manifest at radii as small as 0.15a,
so the majority of the disc area contributes to the mean
precession rate. We note that Goodchild & Ogilvie (2006)
required a very low semi-thickness parameter (h = 0.003) to
fit observations, whereas in our simulations this is an order
of magnitude higher (Table 2) and in better agreement with
observational constraints on, and theoretical expectations
of, disc thickness.
Empirically, systems which show superhumps generally
have q < 0.24. There are two classes of exceptions: nova-likes
and UGem, and magnetic systems. For the former, Osaki
(2005) pointed out that Paczynski’s derivation of Rtides as-
sumes the disc to be cold. A sufficiently hot accretion disc
may expand beyond this owing to a weakening of shocks
by strong pressure effects at the last non-intersecting orbit.
Osaki (2005) argues the persistently high-state nova-likes,
and the unusually long 1985 outburst of UGem in which su-
perhumps were reported, may satisfy this temperature cri-
terion. The magnetic systems, TVCol and possibly TLeo,
may have discs which are pushed out by magnetic forces
(Retter et al. 2003).
The high resolution 3D simulations presented here re-
produce this distribution well; the upper limit for which su-
perhumps are observed is also q ≃ 0.24, albeit at this mass
ratio it is a very slow process indeed. This also compares well
with the theoretical upper limit of q . 0.25. For q = 0.2422 it
appears that the eccentricity starts to grow more than once,
but fails each time and the eccentricity falls away again,
presumably because the encounter with the resonance was
only marginal. The final eccentricity is highest for q ∼ 0.1,
and declines gradually as q becomes less extreme. This is
expected as the disc has more room to grow beyond the 3:1
resonance before being tidally truncated for extreme mass
ratios. For the two most extreme mass ratio simulations,
with q = 0.0526 and q = 0.0256 the final eccentricity is neg-
ligible. In the q = 0.0256 case there is no eccentricity growth.
The q = 0.0526 case is interesting: as the disc is growing, the
eccentricity initially increases and superhumps are appar-
ent, albeit weak and ill-formed (Figure 3). This eccentricity
is, however, damped away again. What is the source of the
damping? One possibility could be the action of the 2:1 reso-
nance which may be excited in ultra-low mass ratio systems
and which can act to damp eccentricity (Lubow 1991b). Per-
haps this results in a close competition between the 2:1 and
3:1 resonance where the 2:1 resonance is marginally excited
in the case of q = 0.0526, and the 2:1 resonance only comes
out on top once the disc has grown and this resonance is
sufficiently populated. In reality, the most extreme mass ra-
tio system in which superhumps are observed is the black
hole X-ray transient KVUMa (XTE J1118 +480), which has
q = 0.037± 0.007. In 2D, the range for which the simulated
discs become eccentric extends to much higher mass ratios
(q 6 0.4815) (Table 1, Figure 14). The confinement to 2D
means that the character of the resonance is different and
the strength of the resonance is increased.
We find that the growth rate of the eccentricity (the
growth rate of S(1,0)) is highly dependent on q, with high
mass ratio systems taking a very long time indeed to be-
come eccentric (Table 1). We see in Figure 2 that our new
simulations have slower growth rates than previous 2D sim-
ulations by Murray (1998). Observationally, the rise times
of superoutbursts is of order a day or so, and superhumps
are generally detected within a day or so. The outburst on-
set is, however, probably governed by the thermal-viscous
instability and hence its timescale is independent of the disc
eccentricity. Generally there are no observations suitable for
assessing whether or not the disc was eccentric before the
superoutburst began, with intensive observations beginning
c© 2005 RAS, MNRAS 000, 1–18
16 A. J. Smith et al.
after the rise to outburst, so timescales on which SUUMa
discs typically develop their eccentricity is ill-constrained.
We note further that our simulations do not include the
thermal-viscous instability, so their evolution differs from
that of SUUMa discs. Lubow (1991a) found an analytical
expression for the growth rate of the eccentricity, finding it
to be proportional to q2. This is applied to an ideal narrow
fluid ring. Clearly the eccentricity growth rate of our simu-
lated discs is not proportional to q2. Osaki (2005) uses the
dependence of eccentricity growth on q2 to propose, as a re-
finement to the TTI model, an explanation for type A/type
B superoutbursts, namely those which show a precursor and
those which do not. He suggests that it depends on whether
the eccentricity growth rate is large enough to excite a signif-
icantly eccentric disc within the duration of a normal out-
burst, and that this is why most SUUMa systems having
relatively low mass ratio show only Type B superoutbursts.
However, if it is the case that high q means lower growth
rates as we see, then this argument fails, and probably other
factors contribute.
Goodchild & Ogilvie (2006) also find extremely low
growth rates in their analytical work. They use the growth
term found by Lubow to describe the rate at which eccentric-
ity is created at the resonance, but consider further how this
eccentricity propagates through the disc. They explain their
low growth rates as due to the eccentricity being strongly
suppressed at the resonance itself. They find that the growth
rates depend most strongly on mass ratio and on bulk viscos-
ity, with further weaker dependence on disc semi-thickness.
In column 7 of Table 2 we list the parameter, h, referred
to by Goodchild & Ogilvie (2006) as the characteristic disc
semi-thickness, and in column 11 the bulk viscosity at R3:1 in
our simulations. For the 3D simulations, these are ∼ 0.036
and 0.2 respectively. Comparing our Figure 2 with Good-
child & Ogilvie (2006)’s results, we see that the trend in our
3D points matches quite well with the h = 0.01 line (the
highest value of h given) in their figure 10. In particular we
see in both cases low growth rates at high q, a maximum
at q ∼ 0.08 and a steep drop in growth rate at mass ratios
below this. Our growth rates, though, are about a factor of
10 higher. Looking to their figure 9, this could be explained
by our high semi-thickness parameter. Our bulk viscosity,
though, is also rather high. It would be very interesting to
make a study of empirical growth rates and their depen-
dence on q and on other known parameters to compare with
these findings. This would require systematic monitoring of
dwarf novae to catch the onset of outburst.
There are, however, distinct differences between the
work of Goodchild & Ogilvie (2006) and the simulations
that we present. Tidal modes, which would presumably act
to truncate the eccentric mode, are not included in their
work. The outer boundary conditions differ. We see too in
Section 3 that the eccentricity distribution as a function of
radius appears quite different from their findings. The situa-
tion in the simulations is further complicated by the presence
and importance of the tidal 2-armed spiral structure which
is not included in Goodchild & Ogilvie (2006)’s work.
In our simulations we see the superhump period de-
creasing (Figure 1) as is often observed over the course of
a superoutburst (Patterson et al. 1993a). As the period
changes, the disc eccentricity is increasing, and, in most
cases, the disc mass has begun to decrease in response to
the enhanced tidal torques. The superhump period decrease
in the simulations can then be explained by the eccentric
wave propagating inward, and additionally by radial shrink-
ing of the disc. We are unsure how to explain the exception,
q = 0.0526, which shows an opposite behaviour for the ini-
tial part of the simulation. Perhaps it is related to the ideas
of Uemura et al. (2005), where they suggest an explanation
for +ve Ṗsh observed in a few cases. They suggest Ṗsh is
related to the amount of matter beyond R3:1, so allowing
for an outward propagation of the eccentric wave. For low
q, the distance between R3:1 and Rtides is greater.
The observational data points to a many-valued ǫ(q)
relation (Figure 14). In particular 3 of the BHXRTs show
systematically lower precession rates than those of CVs.
QZVul (GS 2000+2) might be an exception to this. However
the superhump period measurement is uncertain as it has
been sparsely observed. The microquasar V1487 Aqr (GRS
1915+105) is a much longer period system and may not
be comparable to the other BHXRTs. The accretion discs
in BHXRTs are irradiated by the central X-ray source. We
would expect these discs to be both hotter and also thicker
due to the bloating effect of irradiation. Both of these would
act to reduce the precession rate, due to the retrograde ef-
fect of pressure and due to the dependence on semi-thickness
found by Goodchild & Ogilvie (2006). It would be very inter-
esting if ǫ and q for further BHXRTs could be determined.
The only ultra-compact helium binary included, AMCVn,
also lies below the main cluster of points in Figure 14. As
Roelofs et al. (2006) noted, the helium accretion disc in
AMCVn could be thicker than its hydrogen-rich counter-
parts.
6 SUMMARY
The main findings in this work can be summarised as follows:
• We present improved accretion disc simulations for a
range of q. The main improvements are both numerical and
physical: a higher mass resolution, extension to 3D, more
realistic disc temperature and viscosity, and a radial depen-
dence of sound speed appropriate to a steady-state accretion
disc. We ran the simulations until equilibrium was reached.
For 0.08 < q < 0.24 the 3D discs reach an eccentric equilib-
rium and show a superhump signal in their energy dissipa-
tion rate (which we refer to as a simulated lightcurve).
• The ǫ(q) dependence for the SPH simulations presented
in this work shows a greatly improved match with observa-
tion than previous simulations.
• No high resolution 3D disc with q > 0.24 developed
superhumps. This agrees with theoretical expectations and
matches the majority of observations.
• The region of the disc within R3:1 is overwhelmingly re-
sponsible for generating the dissipation-powered superhump
in the simulations.
• If the difference between observed precession rates and
dynamical precession rates calculated at the 3:1 resonance
radius is due to averaging over the disc as Goodchild &
Ogilvie (2006) suggest, then we find that the best fit char-
acteristic radius of the eccentricity distribution at which the
dynamical precession rate is evaluated to be 0.87R3:1 and
0.83R3:1 for the 3D simulated discs and the observed eclips-
ing systems respectively. The differences between these two
c© 2005 RAS, MNRAS 000, 1–18
Comprehensive simulations of superhumps 17
best-fit radii may be partly due to the differing surface den-
sity distributions in the two cases.
• Our simulations show the effect of the increased effi-
ciency of tidal return of angular momentum to the binary
for an accretion disc which has become eccentric. The disc
mass approaches a new lower steady-state value as the disc
becomes eccentric. This is exactly as asserted by the TTI
model. With the assumption of a radial dependence of vis-
cosity, we deduce an effective ∼ 4 per cent increase viscous
torque between a disc which is circular and one that is ec-
centric. The increase depends on q.
• As the eccentricity grows and the disc mass falls, the
superhump period decreases.
• The dependence of eccentricity growth rates on q that
we see in the simulations presented here is comparable to
the work of Goodchild & Ogilvie (2006). Particularly, we
find that for high mass ratios the growth rates are very low
indeed, in contrast to the result of Lubow (1991a). This
needs to be reconciled with observation.
• We show that superhumping discs have noticeable ec-
centricity even in their inner regions (r ∼ 0.15a). Conversely,
non-superhumping discs are seen to be eccentric only in
their outer regions. In this case however, this ‘eccentricity’ is
steady-state and has origin in tidal distortions, being there-
fore different from that which dominates the main body of
the superhumping discs. We characterise the eccentricity dis-
tributions using two different methods.
• The disc motions can be described as superpositions of
the S(k,l) modes, and the resonance radii act as nodes and
antinodes in the complex standing wave dynamics the disc
executes over a full precession period. We characterise the
disc motions on Psh, the key timescale for the powering of
the observed superhumps.
• The 4:1 and 5:1 resonances may play roles in the dy-
namics of eccentric discs for q < 0.24. This may explain
why the observed precession rates are closer to the dynamic
precession rate at the 4:1 resonance than they are to the
dynamic precession rate at the 3:1 resonance.
• The observational data shows a multi-valued ǫ(q) rela-
tion. In particular, the BHXRTs show systematically lower
precession rates than those of the CVs, which may be ex-
pected when the higher temperature and thickness of their
irradiated discs is considered.
7 ACKNOWLEDGEMENTS
Helpful comments from the referee were much appreciated.
We acknowledge the use of the supercomputing facilities
at the Centre for Astrophysics and Supercomputing, Swin-
burne University of Technology. AJS was supported by a
PPARC studentship.
REFERENCES
Araujo-Betancor S., et al., 2003, ApJ , 583, 437
Arenas J., Catalán M. S., Augusteijn T., Retter A., 2000, MN-
RAS , 311, 135
Augusteijn T., 1994, A&A , 292, 481
Bailey J., 1979, MNRAS , 188, 681
Bailyn C. D., 1992, ApJ , 391, 298
Baptista R., Steiner J. E., Cieslinski D., 1994, ApJ , 433, 332
Baptista R., Catalan M. S., Horne K., Zilli D., 1998, MNRAS ,
300, 233
Baptista R., Jablonski F., Oliveira E., Vrielmann S., Woudt P. A.,
Catalán M. S., 2002, MNRAS , 335, L75
Baptista R., Borges B. W., Bond H. E., Jablonski F., Steiner
J. E., Grauer A. D., 2003, MNRAS , 345, 889
Bruch A., Beele D., Baptista R., 1996, A&A , 306, 151
Casares J., Martin E. L., Charles P. A., Molaro P., Rebolo R.,
1997, New Astronomy, 1, 299
Charles P. A., Kidger M. R., Pavlenko E. P., Prokof’eva V. V.,
Callanan P. J., 1991, MNRAS , 249, 567
Eggleton P. P., 1983, ApJ , 268, 368
Feline W. J., Dhillon V. S., Marsh T. R., Stevenson M. J., Watson
C. A., Brinkworth C. S., 2004a, MNRAS , 347, 1173
Feline W. J., Dhillon V. S., Marsh T. R., Stevenson M. J., Watson
C. A., Brinkworth C. S., 2004b, MNRAS , 354, 1279
Feline W. J., Dhillon V. S., Marsh T. R., Brinkworth C. S., 2004c,
MNRAS , 355, 1
Feline W. J., Dhillon V. S., Marsh T. R., Watson C. A., Littlefair
S. P., 2005, MNRAS , 364, 1158
Foulkes S. B., Haswell C. A., Murray J. R., Rolfe D. J., 2004,
MNRAS , 349, 1179
Foulkes S. B., Haswell C. A., Murray J. R., 2006, MNRAS , 366,
Frank J., King A. R., Raine D. J., 1985, Accretion power in as-
trophysics. Cambridge and New York, Cambridge University
Press, 1985, 283 p.
Goodchild S., Ogilvie G., 2006, MNRAS , 368, 1123
Haefner R., Metz K., 1985, A&A , 145, 311
Harlaftis E. T., Greiner J., 2004, A&A , 414, L13
Harlaftis E. T., Horne K., Filippenko A. V., 1996, PASP , 108,
Haswell C. A., 1996, in IAU Symp. 165: Compact Stars in Bina-
ries, van Paradijs J., van den Heuvel E. P. J., Kuulkers E.,
eds., pp. 351–371
Haswell C. A., King A. R., Murray J. R., Charles P. A., 2001,
MNRAS , 321, 475
Hessman F. V., Mantel K.-H., Barwig H., Schoembs R., 1992,
A&A , 263, 147
Hoard D. W., Linnell A. P., Szkody P., Fried R. E., Sion E. M.,
Hubeny I., Wolfe M. A., 2004, ApJ , 604, 346
Horne K., Cook M. C., 1985, MNRAS , 214, 307
Horne K., Wood J. H., Stiening R. F., 1991, ApJ , 378, 271
Hōshi R., 1979, Progress of Theoretical Physics, 61, 1307
Howell S., Szkody P., 1988, PASP , 100, 224
Imada A., et al., 2006, PASJ , 58, 143
Ishioka R., et al., 2002, A&A , 381, L41
Kato T., 1997, PASJ , 49, 583
Kato T., Mineshige S., Hirata R., 1995, PASJ , 47, 31
Kato T., Nogami D., Moilanen M., Yamaoka H., 2003, PASJ , 55,
Kunze S., Speith R., Hessman F. V., 2001, MNRAS , 322, 499
Lasota J.-P., 2001, New Astronomy Review, 45, 449
Lemm K., Patterson J., Thomas G., Skillman D. R., 1993, PASP ,
105, 1120
Lubow S. H., 1991a, ApJ , 381, 259
Lubow S. H., 1991b, ApJ , 381, 268
Lubow S. H., 1992, ApJ , 398, 525
Lubow S. H., Shu F. H., 1975, ApJ , 198, 383
Maehara H., Hachisu I., Nakajima K., 2006, preprint
(astro-ph/0611519 v1)
Mennickent R. E., Arenas J., 1998, PASJ , 50, 333
Mennickent R. E., Diaz M., 1996, A&A , 309, 147
Meyer F., Meyer-Hofmeister E., 1981, A&A , 104, L10
Monaghan J. J., 1992, ARA&A , 30, 543
Murray J. R., 1996, MNRAS , 279, 402
Murray J. R., 1998, MNRAS , 297, 323
Murray J. R., 2000, MNRAS , 314, L1
c© 2005 RAS, MNRAS 000, 1–18
http://arxiv.org/abs/astro-ph/0611519
18 A. J. Smith et al.
Murray J. R., Armitage P. J., 1998, MNRAS , 300, 561
Murray J. R., de Kool M., Li J., 1999, ApJ , 515, 738
Murray J. R., Warner B., Wickramasinghe D. T., 2000, MNRAS ,
315, 707
Naylor T., Allan A., Long K. S., 2005, MNRAS , 361, 1091
Neil E. T., Bailyn C. D., Cobb B. E., 2006, preprint
(astro-ph/0610480)
Nogami D., Baba H., Kato T., Novak R., 1998, PASJ , 50, 297
O’Donoghue D., Charles P. A., 1996, MNRAS , 282, 191
Orosz J. A., 2001, The Astronomer’s Telegram, 67, 1
Orosz J. A., Bailyn C. D., McClintock J. E., Remillard R. A.,
1996, ApJ , 468, 380
Osaki Y., 1974, PASJ , 26, 429
Osaki Y., 1989, PASJ , 41, 1005
Osaki Y., 2003, PASJ , 55, 841
Osaki Y., 2005, Proceedings of the Japan Academy, Ser. B: Phys-
ical and Biological Sciences, Vol. 81, p. 291-305., 81, 291
Paczynski B., 1977, ApJ , 216, 822
Patterson J., 1979, AJ , 84, 804
Patterson J., 2001, PASP , 113, 736
Patterson J., Bond H. E., Grauer A. D., Shafter A. W., Mattei
J. A., 1993a, PASP , 105, 69
Patterson J., Thomas G., Skillman D. R., Diaz M., 1993b, ApJS ,
86, 235
Patterson J., Halpern J., Shambrook A., 1993c, ApJ , 419, 803
Patterson J., Jablonski F., Koen C., O’Donoghue D., Skillman
D. R., 1995, PASP , 107, 1183
Patterson J., et al., 1997, PASP , 109, 1100
Patterson J., Richman H., Kemp J., Mukai K., 1998, PASP , 110,
Patterson J., Kemp J., Jensen L., Vanmunster T., Skillman D. R.,
Martin B., Fried R., Thorstensen J. R., 2000a, PASP , 112,
Patterson J., et al., 2000b, PASP , 112, 1584
Patterson J., et al., 2002a, PASP , 114, 721
Patterson J., et al., 2002b, PASP , 114, 1364
Patterson J., et al., 2003, PASP , 115, 1308
Patterson J., et al., 2005, PASP , 117, 1204
Pavlenko E. P., Shugarov S. Y., 1999, A&A , 343, 909
Pearson K. J., 2006, MNRAS , 371, 235
Pratt G. W., Hassall B. J. M., Naylor T., Wood J. H., 1999,
MNRAS , 307, 413
Retter A., Hellier C., Augusteijn T., Naylor T., Bedding T. R.,
Bembrick C., McCormick J., Velthuis F., 2003, MNRAS , 340,
Ritter H., Kolb U., 2003, A&A , 404, 301
Robinson E. L., Shafter A. W., Hill J. A., Wood M. A., Mattei
J. A., 1987, ApJ , 313, 772
Rodŕıguez-Gil P., Gänsicke B. T., Hagen H.-J., Marsh T. R., Har-
laftis E. T., Kitsionas S., Engels D., 2005, A&A , 431, 269
Roelofs G. H. A., Groot P. J., Nelemans G., Marsh T. R., Steeghs
D., 2006, MNRAS , 371, 1231
Rolfe D. J., Haswell C. A., Patterson J., 2001, MNRAS , 324, 529
Schoembs R., 1986, A&A , 158, 233
Schoembs R., Vogt N., 1981, A&A , 97, 185
Semeniuk I., Olech A., Kwast T., Nalezyty M., 1997, Acta Astro-
nomica, 47, 201
Shafter A. W., 1983, PhD thesis, UCLA
Shafter A. W., Hessman F. V., 1988, AJ , 95, 178
Shafter A. W., Szkody P., 1984, ApJ , 276, 305
Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337
Skidmore W., Wynn G. A., Leach R., Jameson R. F., 2002, MN-
RAS , 336, 1223
Skillman D. R., Patterson J., Thorstensen J. R., 1995, PASP ,
107, 545
Skillman D. R., Patterson J., Kemp J., Harvey D. A., Fried R. E.,
Retter A., Lipkin Y., Vanmunster T., 1999, PASP , 111, 1281
Skillman D. R., et al., 2002, PASP , 114, 630
Smak J., 1993, Acta Astronomica, 43, 121
Smak J., 1999, Acta Astronomica, 49, 391
Smak J., Waagen E. O., 2004, Acta Astronomica, 54, 433
Smith A. J., Haswell C. A., Hynes R. I., 2006, MNRAS , 369, 1537
Smith D. A., Dhillon V. S., 1998, MNRAS , 301, 767
Stanishev V., Kraicheva Z., Boffin H. M. J., Genkov V., Papadaki
C., Carpano S., 2004, A&A , 416, 1057
Steeghs D., Perryman M. A. C., Reynolds A., de Bruijne J. H. J.,
Marsh T., Dhillon V. S., Peacock A., 2003, MNRAS , 339, 810
Thorstensen J. R., Taylor C. J., 1997, PASP , 109, 1359
Thorstensen J. R., Fenton W. H., Patterson J. O., Kemp J., Krajci
T., Baraffe I., 2002, ApJ , 567, L49
Torres M. A. P., Callanan P. J., Garcia M. R., Zhao P., Laycock
S., Kong A. K. H., 2004, ApJ , 612, 1026
Truss M., 2007, MNRAS , 376, 89
Uemura M., et al., 2000, PASJ , 52, L9
Uemura M., et al., 2002a, PASJ , 54, 285
Uemura M., et al., 2002b, PASJ , 54, L15
Uemura M., et al., 2004, PASJ , 56, 141
Uemura M., et al., 2005, A&A , 432, 261
van Amerongen S., Damen E., Groot M., Kraakman H., van
Paradijs J., 1987, MNRAS , 225, 93
Vogt N., 1981, PhD thesis, Bochum University
Vogt N., 1982, ApJ , 252, 653
Vogt N., 1983, A&A , 118, 95
Vogt N., Semeniuk I., 1980, A&A , 89, 223
Wade R. A., Horne K., 1988, ApJ , 324, 411
Walker W. S. G., Marino B. F., Freeth G., 1976, Informational
Bulletin on Variable Stars, 1185, 1
Warner B., O’Donoghue D., 1988, MNRAS , 233, 705
Webb N. A., Naylor T., Ioannou Z., Charles P. A., Shahbaz T.,
2000, MNRAS , 317, 528
Whitehurst R., 1988, MNRAS , 232, 35
Wood J., Horne K., Berriman G., Wade R., O’Donoghue D.,
Warner B., 1986, MNRAS , 219, 629
Wood J. H., Horne K., Berriman G., Wade R. A., 1989, ApJ ,
341, 974
Zhang E.-H., Robinson E. L., Nather R. E., 1986, ApJ , 305, 740
Zurita C., et al., 2002, MNRAS , 333, 791
c© 2005 RAS, MNRAS 000, 1–18
http://arxiv.org/abs/astro-ph/0610480
	Introduction
	SPH simulations of a Precessing Accretion Disc
	Simulation details
	Simulation Results
	Eccentricity distribution in accretion discs
	Eccentricity distribution from particle trajectories
	Eccentricity from the mass distribution
	Period Excess versus Mass Ratio: Drawing together observation, theory and simulation
	Discussion
	Summary
	Acknowledgements
	REFERENCES

  (Abridged) We use 3D SPH calculations with higher resolution, as well as with
more realistic viscosity and sound-speed prescriptions than previous work to
examine the eccentric instability which underlies the superhump phenomenon in
semi-detached binaries. We illustrate the importance of the two-armed spiral
mode in the generation of superhumps. Differential motions in the fluid disc
cause converging flows which lead to strong spiral shocks once each superhump
cycle. The dissipation associated with these shocks powers the superhump. We
compare 2D and 3D results, and conclude that 3D simulations are necessary to
faithfully simulate the disc dynamics. We ran our simulations for unprecedented
durations, so that an eccentric equilibrium is established except at high mass
ratios where the growth rate of the instability is very low. Our improved
simulations give a closer match to the observed relationship between superhump
period excess and binary mass ratio than previous numerical work. The observed
black hole X-ray transient superhumpers appear to have systematically lower
disc precession rates than the cataclysmic variables. This could be due to
higher disc temperatures and thicknesses. The modulation in total viscous
dissipation on the superhump period is overwhelmingly from the region of the
disc within the 3:1 resonance radius. As the eccentric instability develops,
the viscous torques are enhanced, and the disc consequently adjusts to a new
equilibrium state, as suggested in the thermal-tidal instability model. We
quantify this enhancement in the viscosity, which is ~10 per cent for q=0.08.
We characterise the eccentricity distributions in our accretion discs, and show
that the entire body of the disc partakes in the eccentricity.

Introduction
	SPH simulations of a Precessing Accretion Disc
	Simulation details
	Simulation Results
	Eccentricity distribution in accretion discs
	Eccentricity distribution from particle trajectories
	Eccentricity from the mass distribution
	Period Excess versus Mass Ratio: Drawing together observation, theory and simulation
	Discussion
	Summary
	Acknowledgements
	REFERENCES
