1. law, govt and politics
2. armed forces
3. navy

How to find the mass donor in
accreting double white dwarf binaries?
Master thesis by
Remco Groenendijk
Supervisor
Dr. Paul Groot
18 August 2004
Department of Astrophysics
University of Nijmegen
Toernooiveld 1
6525 ED Nijmegen
The Netherlands
Archived under code AF 10
Typeset using LATEX
Email: remcogro@sci.kun.nl
1
1
Summary
In this thesis a theoretical model for interacting double white dwarf binaries is
constructed. These are systems containing two white dwarfs orbiting each other.
These stars are so close to each other that mass is being transferred between
them. They are also called AM CVn stars, named after the first system of this
kind to be detected.
These systems emit a very specific spectrum. The observations of periodically changing radial velocity signatures in a helium dominated accretion disk
tell us that they are binaries and that the donor star must have mainly helium
in its outer layers. Because a helium star would outshine the other components,
it is concluded that the donor is a helium white dwarf or a semi-degenerate
remnant of a giant star. In the spectrum of these systems no radiation has yet
been identified as coming from the donor. The goal of this thesis is to make
a model which can be used to predict under what circumstances and how we
can observe the donor in these kind of systems. This gives observers a handle
on how to find these objects. Once found the masses of both white dwarfs can
be measured. The anticipated mass distribution of the two components in AM
CVn stars puts strong constraints on evolutionary models of binaries.
I have constructed a model of an AM CVn star, taken a few free parameters
and have tried to find out for which values of these parameters the donor star
can be easily observed. I found that in most cases the donor star can not be
seen. It might become visible in the first ∼0.1 or ∼1 Myr of mass transfer by
irradiation coming from an impact spot of the accreted matter on the primary.
I have calculated for different white dwarf masses how long this phase lasts.
Under the most favourable conditions the donor might also be more luminous
than the accretor after ∼ 108 yr. This could happen when the donor is a helium
white dwarf and the accretor is a carbon white dwarf, because helium white
dwarfs cool more slowly than carbon white dwarfs do. But for this to happen
the heating of the accretor, as a result of mass transfer, needs to be negligible.
2
Contents
1 Summary
2
2 Introduction
2.1 The evolution of single stars . .
2.2 The evolution of binaries . . . .
2.3 The evolution of AM CVn stars
2.4 The stability of AM CVn stars
2.5 Why we look at AM CVn stars
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4
4
7
8
10
11
3 Method and approach
15
3.1 Irradiation of the donor . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Time evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 White dwarf cooling . . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Results
23
4.1 The effect of different white dwarf masses . . . . . . . . . . . . . 24
4.2 The effect of the cooling time . . . . . . . . . . . . . . . . . . . . 28
4.3 The effect of the accretor’s spin . . . . . . . . . . . . . . . . . . . 28
5 Discussion and conclusions
32
6 Future work
34
7 Appendix
38
3
2
Introduction
In this paper I have constructed a theoretical model for double white dwarfs
with mass transfer. In this section first a short introduction is given on the
evolution of these kinds of systems, starting with evolution of single stars and
continuing with binary stars. Then the evolution of double white dwarfs with
mass transfer is discussed. Finally I explain why it is important to study these
systems in detail and how they help us to better understand binary evolution
and accretion physics.
2.1
The evolution of single stars
Stars are born out of clouds of gas. These clouds contract under their own
gravity. This happens for a spherical cloud when its mass is bigger than the
so-called Jeans Mass MJ given in units of solar masses (M¯ ) by (see for instance
Prialnik 2000 [27])
3
T2
(1)
MJ ≈ 105 √ M¯ ,
n
where T is the temperature of the cloud in Kelvin and n the number of gas
particles per m3 .
Figure 1: M16, a star forming region.
As a cloud contracts its density increases and smaller parts of the cloud
will also exceed their Jeans Mass and start contracting themselves. At first,
the gravitational energy that is released in the collapse of the cloud is radiated
away. But as the small clouds get more and more dense, they get more opaque
and the gravitational energy is converted more and more into heat. This heat
causes an outward pressure given approximately by the ideal gas law:
Pgas =
<
ρT
µ
where
4
1
1
1
+ ,
≡
µ
µI
µe
(2)
where Pgas is the gas pressure, < = 8.3 × 103 J kg−1 K−1 is the ideal gas
−1
constant, ρ the density, µ−1
I the mean atomic mass and µe the average number
of free electrons per nucleon. For solar composition µ=0.61. This pressure slows
down the contraction. But without an internal heat source, the cloud has to
keep contracting to create a pressure that balances the force of gravity. The
small cloud is then called a protostar. In the end at the center of these clouds,
depending on their mass, the density and temperature will eventually rise to
about 105 kg m−3 and 107 K. At this point hydrogen fusion starts and a new
star is born.
The fusion causes the temperature at the center of the star to remain high.
Therefore the pressure, as given by equation 2, remains high too. This stops the
contraction for a long time. The star is then said to be on the main sequence
and it remains there for about 90% of its active life. Consequently, most stars
including our Sun are on the main sequence. Hydrogen is turned into helium
and slowly the mass percentage of helium in the core increases. After hydrogen
is exhausted in the core, there is no more heat supply in the inert He-core. As
a result of that it starts to contract. This causes its temperature to increase
and hydrogen starts to burn in a shell surrounding the core. While the core is
contracting the outer layers of the star expand to giant dimensions; the star is
then called a red giant. In this process the core heats up and finally helium can
ignite returning the star to a semi-stable phase of He-core burning or the He
main sequence.
Figure 2: The Crab nebula, the remnant of a supernova that was observed
almost 1000 years ago.
After helium is exhausted, the temperature of the core can rise again if the
mantle of the star is massive enough. In practice this happens for stars that
were more massive than 8 M¯ at the start of their main sequence. Heavier and
heavier elements can fuse until finally iron is formed, after which a supernova
(SN) explosion forms a neutron star or a black hole. The amount of energy
released in a SN is huge. During a SN a single star can produce more light than
5
the entire galaxy it is in (some 1011 stars!). SN are classified observationally in
SN Type I where no hydrogen lines are observed and in Type II where there is
hydrogen present.
Most stars (the ones with mass M < 8M¯ ) do not finish this evolutionary
path. Once a carbon- oxygen-core is formed, the temperature and pressure in
the core never increases to a point where carbon can fuse and that is why the
evolution stops. Now the core could start contracting further to increase the
temperature and the pressure again to a point where carbon fusion would start.
However this doesn’t happen. The reason for this has to with the Pauli exclusion
principle.
The Pauli exclusion principle states that two fermions can not be in the
same state. This happens for instance in atoms; the outer electrons have to
have a higher quantum number for their energy, because all the lower energy
states are already filled with other electrons. The same thing can happen in a
compact star. The electrons are packed so closely together that some of them
have to have a higher energy and subsequently a higher corresponding velocity.
This causes a pressure called degeneracy pressure. This pressure can prevent a
further collapse of a stars. If a burnt-out star is supported by this pressure, it
slowly cools and fades away from our view. Such a star is called a white dwarf
(WD). Now it is easy to see why in general the radius of a WD is smaller for
a bigger mass (see equation 6). When a WD is more massive, it requires more
pressure to prevent it from collapsing. This means the electrons must have a
higher average velocity. This happens if they are packed together even more
tightly. So the WD must have a smaller radius.
The degeneracy pressure depends on the speed of the electrons. And of
course this speed can not be larger than the speed of light c. This causes a
maximum to this pressure. This maximum pressure corresponds to a maximum
WD mass. This mass is called the Chandrashekar mass ' 1.4 Solar masses
(= 1.4 M¯ ).
Figure 3: Mass distribution of DA white dwarfs, based on multichannel spectroscopy for 70 stars done by Weidemann and Koester (1984) [38]
During their lives, stars blow away a large fraction of their mass in a wind.
6
That is why it is estimated that stars with a mass up to ±8M¯ (Prialnik [27]
page 165 ) can evolve to white dwarfs. From observations it is known that the
mass distribution of white dwarfs shows a clear peak, with two tails (see figure
3). The peak is at a mass of ∼ 0.6M¯ a tail is found between 0.2M¯ and
0.4M¯ and to higher masses (> 0.6M¯ ). The first tail belongs to white dwarfs
with a helium core. Their progenitors were not heavy enough to ever ignite
helium. The peak belongs to carbon oxygen white dwarfs, whose progenitors
did ignite helium. In the final phases of their lives these stars undergo thermal
pulses caused by the alternate burning of hydrogen and helium in thin shells.
This causes a large stellar wind that blows away most of their mass and lets
them evolve to ∼ 0.6M¯ carbon white dwarfes. Finally there are some very
rare oxygen-neon-magnesium white dwarfes most of which have masses near the
Chandrashekar mass.
2.2
The evolution of binaries
In the past century our understanding of the evolution of single stars has increased dramatically. There are still some loose ends, but most evolutionary
phases, especially the main sequence, are well understood. This creates a new
opportunity for astrophysicists. Knowing that more than 50 % of all the observed stars are in a binary system (Duquennoy and Mayor 1991[3]), it is an
obvious step to go and study their evolution in more detail.
L1
M
2
~R 2
M1
a
~R 1
Figure 4: The Roche lobe; R1 and R2 are defined as the radii of spheres having
a volume equal to the volume of the corresponding Roche lobe.
The simplest binary systems are wide binary systems. These stars orbit
one another, but their separation is so large that they do not interact during
their evolution. This allows them to evolve just like two single stars. The more
interesting ones are close binary systems. In these systems the stars are so
close to each other that mass is transferred between at least once during their
evolution. A good way of understanding what happens then, is by looking at
the so-called Roche lobe of such a system (see figure 4). The heavier star has a
mass M1 and the lighter a mass M2 .
7
Looking at the system in a co-rotating frame there is a point where all the
forces exactly balance and where the potential goes through a saddle point.
This point is called the inner Langragian point (L1 ). The Roche lobe is an
equipotential line through this point. If either of the stars fills its Roche lobe,
than mass will be transferred from that star to the other. So both the heavier
as well as the lighter star can, in principle, accrete. During binary evolution
Roche lobe overflow can occur several times.
A star can fill its Roche lobe when either it expands or when its Roche lobe
shrinks. Expansion of a star can for instance occur when the star is in its red
giant phase. Shrinking of the Roche lobe occurs when the orbital separation (a)
between the stars decreases. This happens when orbital angular momentum is
lost by, for instance, magnetic braking or gravitational wave emission. There are
many evolutionary paths for binaries. Not all of these paths will be discussed
here, but just the ones leading to the formation of AM CVn stars.
The two stars in a close binary usually cannot be resolved, but the systems
still can be recognised as binaries. This is because the stars in a binary revolve
around each other. The light of the stars is then periodically Doppler shifted to
the blue (or red) if a star is coming towards (or going away from) the observer.
Another way to find semi detached binaries is looking for an accretion disk.
These disks can be found in the spectrum, because they are often hot (producing
UV or X-rays) and the profile of their emission lines has a very specific shape.
Observations tell us that the radius of the disk is about half the radius of the
Roche lobe of the primary. Very often a hot spot can be seen on the disk; this
is the spot where the stream of mass from the donor hits the disk.
2.3
The evolution of AM CVn stars
White dwarf binaries with mass transfer are also called AM CVn stars after
the first system of this kind to be detected by Humason and Zwicky (1947)
[8]. Smak (1967) [32] discovered that AM CVn shows photometric variability
with a period of about ∼18 minutes and suggested that it is a binary. Soon
after that Paczy´
nski (1967) [22] realised that this could be a semi-detached pair
of degenerate white dwarfs in which mass transfer is driven by loss of angular
momentum due to gravitational wave radiation.
Right now there are 13 of these systems known with periods varying from 5.3
to 65.1 minutes (see table 1). From the spectra of AM CVn stars, the absence of
hydrogen features strongly suggests that donors in these systems are extremely
helium rich, because for helium rich accretion discs Balmer emission lines are
stronger than He lines, even if the He/H ratio is ∼100.
Several formation scenarios have been proposed for AM CVn stars. The first
way to form an AM CVn star is discussed by Faulkner et al. (1972) [4]. They
suggested two possibilities for the helium rich donor in AM CVn: (i) a helium
star with a mass of 0.4-0.5M¯ and (ii) a helium white dwarf. The first possibility
is excluded because the helium star would dominate the spectrum and cause
the accretor to have large radial velocity variations, none of which are observed.
Thus they concluded that AM CVn stars are interacting double white dwarfs.
Compact detached double WD binaries are taken as the direct progenitors of
AM CVn stars. These detached binaries are expected to form after a common
envelope phase. In this phase the orbital distance of the system shrinks from
typically 100 R¯ to 1 R¯ (see Nelemans et al. (2001a) [16], hereafter NPVY).
8
Table 1: list of all the known and suspected AM CVn systems.
Object
Porb (s)
Reference
RX J0806
321.25 ∗ Israel et al.(2002) [10];
Ramsay, Hakala and Cropper (2002) [28]
V407 Vul
569.38 ∗ Cropper et al. (1998) [2]
ES Cet
620.26
Warner and Woudt (2002) [36]
AM CVn
1028.7
Solheim et al. (1998) [33];
Skillman et al. (1999) [31]
HP Lib
1102.7
O’Donoghue et al. (1994) [21];
Patterson et al. (2002) [25]
CR Boo
1471.3
Wood et al. (1987) [37];
Patterson et al. (1997) [23]
KL Dra
1500 ∗
Wood et al. (2002) [39]
V803 Cen
1612.0
Patterson et al. (2000) [24]
CP Eri
1701.2
Abbott et al. (1992) [1]
‘2003aw’
2040 ∗
Woudt and Warner (2003) [41]
SDSS J1240-01 2241
Roelofs private communication
GP Com
2974
Nather, Robbinson and Stover (1981) [15];
Marsh, Horne and Rosen (1991) [11]
CE-315
3906
Ruiz et al. (2001) [29];
Woudt and Warner (2001) [40]
Notes: ∗ It is not sure if these periods are orbital periods. They might
be the super hump periods.
The timescale in which the stars spiral in is estimated to be .1 year (Iben (1993)
[9] and references therein). After this phase the two WDs cool and spiral in as
a result of gravitational wave emission. The AM CVn star is formed when the
lighter WD fills its Roch lobe and mass transfer starts. By then the secondary
is expected to be cool enough (see section 3.3) that it can be approximated with
a zero temperature WD.
A second way to form a helium transferring binary in the right period range
was first outlined by Savonije et al. (1986) [30]. They describe a system with
a neutron star as an accretor, but the idea remains the same. Consider a lowmass, non-degenerate helium burning star, with a white dwarf companion. If
the components are close enough, loss of angular momentum via gravitational
wave radiation may result in Roche lobe overflow before helium exhaustion in
the stellar core. When the mass of the helium star decreases below ∼0.2M ¯ ,
core helium burning stops and the star becomes semi-degenerate at a period of
∼10 minutes. This causes the mass-radius relation to invert, which means (see
section 3.2 and equation 24) that their orbital period will start to increase as a
result of mass loss. Such a system will also have the observed characteristics of
an AM CVn star.
A third way to form an AM CVn star is to form it out of a cataclysmic
variable (CV). This possibility was dismissed by NPVY, because its probability
would be very low. Recently Podsiadlowski et al. (2002) [26] have argued that
the probability might not be so low after all. A CV is a binary consisting of a
WD and a main sequence star (or a red giant). This formation route requires
that a normal H-rich star (∼1M¯ ) starts to transfer mass near the end of or just
9
after hydrogen core burning. Such systems become ultra compact binaries with
a period as short as ∼10 minutes where the secondary is initially non-degenerate
and hydrogen-rich, but increasingly becomes degenerate and helium-rich during
its evolution.
2.4
The stability of AM CVn stars
Most of the AM CVn stars are expected to have a helium accretion disk, but at
periods of a few minutes the size of the stars becomes comparable to the orbital
distance. Marsh and Steeghs (2002) [12] have argued that V407 Vul is a system
where the mass stream of the donor hits the accretor directly. To be able to
get a small period, the donor must fill its Roche lobe very late. So one needs
a small radius for the donor. The donors of the first formation scenario (the
detached double WD route) are colder and thus smaller. Therefore the double
white dwarf route is expected to produce the highest number of systems in a
direct impact phase. It is instructive to see the relative sizes of AM CVn stars
for two typical cases (See figure 5). The WDs are assumed to be cold, so that
the zero temperature mass-radius relation (equation 6) can be used.
Sun
R= 7.0 x 10
10
cm
a = 5.8 x 10 9 cm
P= 3 x 102
WD1
s
WD 2
WD1
Earth
a = 3.3 x 10 10 cm
P= 4 x 103 s
WD 2
a = 3.8 x 10 10cm
P= 2.6 x 10 6 s = 1 month
Moon
Figure 5: The relative sizes of AM CVn systems with masses: M1 = 0.6M¯ and
M2 = 0.1M¯ . Their orbital distances a and periods P are plotted in the figure.
Once mass transfer starts, it can be either stable or unstable. The first type
of instability occurs when the accretion is super-Eddington. When the mass
stream hits the accretor, it will produce radiation. The radiation pressure can
prevent more mass from falling in. This amount of luminosity is called the
Eddington luminosity. When that happens a common envelope will form in
which the stars spiral in and finally merge. The Eddington accretion rate is
given by (Marsh 2004 [13])
8πGmp cM1
,
M˙ Edd =
σT (φL1 − φa )
10
(3)
where G is the gravitational constant, mp the mass of a proton, c the speed
of light, M1 the mass of the accretor, σT the Thomson cross-section of the
electron and (φL1 − φa ) the potential drop between the inner Langragian point
and the accretor. In figure 7 the systems that are stable against the Eddington
instability are the ones below the dot-dashed line.
Mass loss can cause the orbital distance to increase or decrease. This depends on the speed at which the angular momentum can be transferred back
to the orbit. When the orbital distance decreases this leads to a second kind of
instability. Marsh et al. 2004 [13] calculate the criteria for stability and they
parameterize the speed at which the angular momentum is fed back to the orbit
by the synchronisation timescale τs . τs is the timescale on which the accretor
synchronises its spin to the orbit of the system. Figure 6 and 7 summarize the
results of Marsh et al. 2004 [13]. The conclusion is that the number of systems,
that remain stable after the onset of mass transfer is highly dependent on τ s ,
which is poorly known. This means that there are still large uncertainties in
modelling the expected number of AM CVn stars.
Figure 6: This figure is from Marsh et al.(2004) [13]. The criteria for stability
for different values of M1 and M2 are plotted. In the upper left part mass
transfer is always unstable. In the lower right part always stable. In between
the stability of the mass transfer depends on the value of τs . Lines that indicate
the stable region for different values of τs are plotted as dotted lines with the
value of τs in years next to them.
2.5
Why we look at AM CVn stars
As mentioned in the last section astrophysicists are trying to develop evolutionary models for all kinds of binaries. Currently, uncertainties in population
models are estimated to be as high as an order of magnitude. Processes such
as common envelope evolution and the onset of mass transfer are poorly understood. Unfortunately these phases last a very short time, so there is almost no
chance of finding a system in such a phase. The only thing that can be done
to improve our understanding of these processes is observing systems that are
11
Figure 7: This figure is from Marsh et al. (2004) [13]. Figure 6 is overplotted
with the masses of double white dwarfs at birth from the binary population
models of NPVY that survive in the case of strong coupling (τs = 0). The surviving models do not reach the limit of guaranteed stability in this case (upper
dashed curve) because super-Eddington accretion becomes the more stringent
constraint (as marked by the dot-dashed line).
going to go into, and have left these phases. AM CVn stars are clearly recognisable end products of binary evolution. So they can systematically be searched
for. Finding the total number of AM CVn systems up to a certain magnitude
puts serious constraints on evolutionary models. If the mass ratio for AM CVn
stars could be determined as well, the constraints would be even stronger.
Up until now no direct emission coming from the donor has ever been measured. If a spectral line could be found in which the donor dominates the
spectrum, the velocity of the donor could be found by looking at the Doppler
shift in that line. With that velocity the mass estimates for both the components of the system can be improved. That is why the goal of my thesis has
been to model how we can find a signal from the donor in the spectrum of an
AM CVn star.
Some AM CVn stars might be progenitors of SN explosions. If the mass
transfer in an AM CVn star pushes the accretor over the Chandrasekhar mass,
then the accretor will collapse. Because there is no more hydrogen in the system,
the explosion will be detected as a SN Type Ia. If it is known how many AM
CVn stars should explode in a SN Type Ia explosion, we know how many of
these events still have to be explained.1
Another thing that should be noted is that AM CVn stars are strong sources
of gravitational waves. The energy loss by means of gravitational radiation is
1 Supernovae Type Ib and Ic also exist, it is expected that these types of explosions occur
when giant stars collapse. These stars have such a strong wind, that they have no more
hydrogen envelopes left at the moment they collapse. That is why they are still detected as a
SN type I.
12
Figure 8: The emission of gravitational waves.
given by the formula (see for instance Paczy´
nski (1967) [22])
32 G4 (M1 M2 )2 (M1 + M2 )
dEgr
=−
dt
5 c5
a5
(4)
for the two masses M1 , M2 moving on a circular orbit with a separation a. c
is the speed of light and G the gravitational constant. This means that heavy
objects with a short orbital distance are the strongest source of gravitational
waves.
Although compact binaries with neutron stars or black holes are even stronger
sources, it turns out that these kind of binaries are far more rare. Consequently
they are on average much further away from us. The detection of gravitational
waves has recently become a hot topic, because of the Laser Interferometer
Space Antenna (LISA). LISA will be the first gravitational wave detector in
space. Its launch is planned for 2012 and its total lifetime is estimated to be
five years. Nelemans et al. (2001c) [17] and Nelemans (2004) [19] find that the
majority of the systems that can be detected by LISA will be double WD binaries: They expect LISA to find ∼11.000 resolved detached double white dwarfs
as well as ∼11.000 resolved AM CVn systems.
AM CVn stars also provide a unique opportunity to study accretion physics.
Accretion disks show a wide variety of (partly) unexplained phenomena such as
dwarf novae outbursts and jets (see figure 9). That is why they are a hot topic
in astrophysics. Many models for these disks are being made. AM CVn stars
provide an excellent opportunity to test these models. Their orbital period is
convenient for observations. And most of all because almost all accretion disks
are dominated by hydrogen, AM CVn stars provide a unique opportunity to
test disk models on helium disks.
Another topic that can be studied with AM CVn stars are super humps.
The small orbital distance in AM CVn stars allows the secondary to perturb
the disk. The disk will become eccentric and will start to precess. This motion
13
Figure 9: A simulation of a disk with jets around a black hole made by Rob
Hynes.
results in super humps in the lightcurve. Simulations of Hirose and Osaki (1990,
1993) [6] and [7] have shown that there is a relation between the orbital (P orb )
and the super hump period (Pprec ) given by:
1+q
Pprec
=A
,
Porb
q
(5)
where A ' 3.73 for systems of q < 0.1. Where q is the mass ratio of the system
(M2 /M1 ). This is why some observed periods might be the super hump period
(see footnote of table 1).
14
3
Method and approach
The goal of my research has been to find a method in which the donor of an
AM CVn star can be identified observationally, because in all the AM CVn
systems there has never been a signal found coming from the donor. So I want
to know what kind of spectrum the donor emits and how that spectrum can be
distinguished from the other components of the system. The donor in AM CVn
systems is rather cool, large and near the accretor. That is why irradiation by
for instance an impact spot might make it visible. I have written a program
that simulates such a system to see when this happens (The complete code is
added in the appendix). I will explain what it does here; first in an overview
and afterwards in more detail.
In subsection 3.1 I show how with two white dwarf masses the Roche lobe and
orbital period of a binary can be calculated. With this information I calculate
where the mass stream from the donor hits the primary and how much energy is
released in that process. Then I calculate how much of the energy gets irradiated
on to the donor.
In subsection 3.2 I show how those parameters evolve with time. I calculate
the mass transfer rate, and with that the values of M1 and M2 as a function of
time.
In order to calculate the evolution of AM CVn systems we first need the
masses and temperatures of the WDs at the birth of the AM CVn system.
There are many possibilities for these numbers (see figure 7 and 13). Therefore
I have tracked the evolution of AM CVn systems for a number of different initial
parameters. Apart from those different cases I distinguish between two extreme
cases:
1. The radiation that is falls on both the primary and the secondary is completely absorbed and causes them to heat up. So I assume that the atmosphere is perfectly conducting all the heat to the core.
2. The radiation that is falls on both the primary and the secondary is processed completely locally and re-emitted there. No heat is conducted to
the cores, so the WDs do not heat up.
Using these cases I finally show in subsection 3.3 how I model the cooling of
both the accretor and the donor.
3.1
Irradiation of the donor
First I model an AM CVn system in the direct impact phase. I assume that
the WDs of the AM CVn system have formed out of a pair of detached double
WDs (the first formation channel; see section 2.3). This means that we can use
Eggleton’s zero temperature mass-radius relation (equation 24 of Marsh et al.
(2004) [13]), to calculate a radius R for a given mass M :
R
R¯
=
0.0114
"
"µ
M
MCh
× 1 + 3.5
µ
¶− 23
M
Mp
15
−
¶− 32
M
MCh
¶ 32 # 21
µ
¶−1 #− 23
µ
+
M
Mp
,
(6)
Figure 10: The ballistic stream as plotted by my program for M1 = 0.5 M¯ and
M2 = 0.25 M¯ , 60,000 years after the onset of mass transfer. The irradiated
surface on the donor has been filled.
where MCh = 1.44M¯ and Mp = 0.00057M¯ . This relation is valid for all white
dwarf masses in the range 0 < M < MCh . We find the orbital separation a
using an equation for the Roche lobe quoted by Warner (1995) ([35] equation
2.5c):
1
2
R2 (0.6q 3 + log(1 + q 3 )
.
(7)
a=
2
0.49q 3
With this a plot of the Roche lobe is made. The orbital period P is given by
Kepler’s third law:
4π 2 a3
P2 =
.
(8)
G(M1 + M2 )
To find out where the stream of mass hits the accretor the trajectory of this
stream is calculated. Besides the two gravitational forces F~g of both stars, there
is also a coriolis force F~cor and a centrifugal force F~cen working on a particle of
the stream of mass m. This is because the Roche lobe is taken in a co-rotating
frame. The forces are given by:
GM2 mˆ
r2
GM1 mˆ
r1
+
,
F~g =
2
2
r1
r2
F~cor = −2m~
ω × ~v ,
F~cen = −~
ω × (~
ω × ~r), (9)
where ~r1 and ~r2 are vectors from the centers of the accretor and the donor to
the particle, ~r is the vector from the center of mass to the particle, v its velocity
r
(~v = d~
~ is its angular velocity (~v = ω
~ × ~r).
dt ) and ω
The mass-stream is then plotted and the place of the impact spot on the
primary is calculated. We want to see if the donor gets irradiated by the impact
spot. To do this we assume that all the energy is converted into radiation at
that spot and that that radiation then spreads out in all directions. This means
for instance that we assume that half of the radiation falls on the accretor.
16
We now calculate how many of the 360◦ are covered by the donor. This
gives a percentage of the radiation that falls on the donor. We approximate the
Figure 11: The projection of the irradiated surface onto a sphere.
donor by a sphere that is partly being irradiated. We can use a y-parameter to
indicate upto how far this sphere is being irradiated. Setting y to y = 0 at the
center of the donor (see figure 11) we get
y=
tan(δ − α/2)
,
tan(α/2)
(10)
where α is the angle in which the donor is seen in the sky from the impact spot
and δ is the angle in which the irradiated surface of the donor is seen in the sky
from the impact spot (figure 11). So y runs from y = 1 for complete irradiation
to y = −1 if there is no irradiation. Now we can calculate the percentage Q of
the cross section of the donor that is being irradiated. This is given by:
Z
2 yp
Q=
1 − y 02 dy 0 .
(11)
π −1
In this way we approximate the projection with a flat circle. This is a good
approximation, because the distance d between the center of the donor and the
impact spot is more than three times the radius of the donor R2 for common
mass ratios (see figure 11). Now we can find out what percentage of the radiation
emitted at the impact spot is being absorbed by the donor. Taking the radiation
to be emitted in X-rays this percentage is given by:
η=
Pab
(1 − A)QπR22
=
,
Pem
4πd2
(12)
where Pab is the absorbed power on the donor and Pem the emitted power of
the radiation from the impact spot, d is the distance between the two stars and
A is the X-ray albedo. Albedo is the fraction of light that is reflected by a body
or surface. It is taken to be constant and equal to 0.4 which roughly agrees with
Milgrom and Salpeter (1974) [14].
The emitted power is calculated using the difference in potential φ between
the start (at L1) and the end of the ballistic stream (at the impact spot on the
17
accretor), because the motion of the particles in the stream is well approximated
by a free-fall. Setting the potential energy U equal to the kinetic energy K, we
get for the velocity of the particles in the stream at the impact point in the
frame of the Rochelobe:
p
vabs = φ(L1 ) − φ(impact spot).
(13)
If the accretor is spinning at an angular velocity ω1 , then the relative velocity
of the particles with respect to the accretor, vrel , will be less than their absolute
velocity. Taking the direction of the normal on the primary at the impact spot
ˆ we can decompose the ~vabs and ~vrel in
as rˆ and the perpendicular direction φ,
these directions:
~vabs
=
~vrel
=
ˆ
vabs r rˆ + vabs φ φ,
ˆ
vrel r rˆ + vrel φ φ.
(14)
(15)
The components of ~vrel are calculated by:
vrel φ
=
vrel r
=
vabs φ − ω1 R1 ,
vabs r
(16)
(17)
and the emitted power follows from:
2
,
Pab = M˙ vrel
(18)
with M˙ as the mass transfer rate. As a result of mass transfer the accretor will
be spun up, but the tidal coupling to the orbit will cause it to spin down. The
strength of this coupling is not well known (Marsh et al. (2004) [13]). That is
why ω1 is taken as a free parameter that can vary between 0 and its break-up
spin rate ω1 max , given by:
s
GM1
.
(19)
ω1 max =
R13
3.2
Time evolution
I have shown how to calculate the amount of radiation from the impact spot
that reaches the donor. Now I will show how that radiation changes with time.
The evolution of a binary is driven by the loss of angular momentum from
the orbit. The orbital angular momentum is given by:
r
Ga
Jorb =
M1 M2 ,
(20)
M
with M = M 1 + M 2. This angular momentum slowly decreases. So taking the
mass to be conserved (M˙ 1 = −M˙ 2 ) we get:
1 a˙
d log(Jorb )
J˙orb
M˙ 2
=
= (1 − q)
+
,
dt
J
M2
2a
(21)
with q = M2 /M1 . The size of the Roche lobe can be approximated by (Paczy´
nski
(1967) [22]):
¶1
µ
a˙
1 M˙ 2
R˙ 2
M2 3
= +
.
(22)
⇒
R2 ' 0.46a
M
R2
a 3 M2
18
Filling in
a˙
a
for equation 22 in equation 21 gives:
à ! ·
¸−1
ζ(M2 ) 5
J˙
M˙ 2
×
,
=
+ −q
M2
J
2
6
(23)
where ζ(M2 ) is given by:
ζ(M2 ) =
R˙ 2 M2
d log(R2 )
.
=
˙
R 2 M2
d log(M2 )
(24)
Using equation 6, it follows that:
ζ2 (M )
=
1
−
3
"µ
M
Mch
¶− 32
−
µ
M
Mch
¶ 23 #−1 µ
M
Mch
¶− 32
"
µ
¶ #−1 " µ
¶#
¶− 23 µ
¶− 23 µ
2
7 M
M
Mp
Mp
×
+ 1 + 3.5
(25)
.
+
+
3
Mp
M
3 Mp
M
Assuming that the angular momentum is only lost from the orbit by the emission
of gravitational waves we get (Landau and Lifshitz (1971) [20]):
à ! à !
J˙
J˙
32 G3 M1 M2 M
.
(26)
=
=−
J
J
5 c5
a4
orb
With these numbers equation 23 gives the mass transfer rate for each value of
M1 and M2 and the evolution an AM CVn system can be calculated. 2 Knowing
the mass transfer rate as a function of the masses allows for the calculation of
the masses as a function of time. The mass transfer rate does not depend on
direct impact or disk accretion in our treatment where the angular momentum
of the primary is being neglected. The evolution of the period as a function of
time is given in figure 12.
In the next few sections I will often plot the orbital period on the X-axis.
There is a one-on-one correspondence between the orbital period and the age of
an AM CVn star for given masses. This is plotted in figure 12.
3.3
White dwarf cooling
To find most probable values for the masses and temperatures WDs at the
onset of mass transfer, models of NPVY and Tutukov (1996) [34] are used.
These models expect these values to vary of a large range see figures 7, 13 and
14.
Therefore we take 3 different combinations of masses and 3 different ages of
the WDs at the birth of the AM CVn star. We use figure 7 to get reasonable
masses. First we assume, that the value of τs ' 1000 yr. We then take two
extreme cases of AM CVn stars, that would still survive the onset of mass
2 A more thorough treatment is given in Marsh et al. (2004) [13], where the stability of
the mass transfer is considered and the spin of the accretor is also taken into account. The
angular momentum stored in this spin is however small. That is why it is neglected here. It
is however very important for the stability of the mass transfer at the moment of Roche lobe
overflow.
19
Figure 12: The evolution of the orbital period for three different sets of masses
at the start of mass transfer: M1 = 0.2 M¯ and M2 = 0.1 M¯ (solid), M1 = 0.5
M¯ and M2 = 0.25 M¯ (dashed), M1 = 1.0 M¯ and M2 = 0.2 M¯ (dotted).
t = 0 yr when the mass transfer starts.
transfer: M1 = 0.2 M¯ , M2 = 0.1 M¯ for low masses and M1 = 1.0 M¯ ,
M2 = 0.2 M¯ for high masses. If we assume that τs ' 0.1 year then most
systems will have masses in the range: M1 = 0.5 M¯ , M2 = 0.25 M¯ . So from
now on I will only track the evolution of these three systems.
To get the temperatures of the WDs at the birth of the AM CVn system
we use figure 13 from Nelemans (private communication) and figure 14. I take
systems that are 1 Myr, 0.1 Gyr and 1 Gyr old at the start of mass transfer.
In order to get a temperature of the WDs at these different ages, I use an
approximation of the cooling curves of Hansen 1999 [5] made by Nelemans et
al. (2004) [19] as shown in figure 15. These curves are a rough translation of
fig. 10 of Hansen (1999) [5].
The heat content U of a WD is stored in the nondegenerate ions.
U∝
TM
,
Amp
(27)
where T is the temperature and M the mass of the WD, A is the baryon number
and mp the proton mass. This means that helium cores (baryon number A=4)
contain more heat than do carbon or oxygen cores (A=12 and 16), thus helium
core WDs are brighter at the same age.
The luminosity of a star is related to its temperature by:
L = 4πR2 σT 4
⇒
1
T = αL 4 ,
(28)
where σ = 5.67 × 10−8 Jm−2 s−1 K−4 is Stefan-Boltzmann‘s constant and α =
1
(4πR2 σ) 4 . If a WD is not heated by accretion, then:
κα − 3 ˙
L 4 L,
−L = U˙ = κT˙ =
4
20
(29)
0.4
1e+10
1e+09
tcool (yr)
m2 (Msun)
0.3
0.2
1e+08
1e+07
0.1
1e+06
0
0.5
1
M1 (Msun)
1.5
2
0
1e+10
1e+10
1e+09
1e+09
tcool (yr)
tcool (yr)
0
1e+08
1e+07
1e+06
1e+06
0.1 0.2 0.3 0.4 0.5 0.6
q = m2/M1
0.2
m2 (Msun)
0.3
0.4
0.5
1
M1 (Msun)
1.5
2
1e+08
1e+07
0
0.1
0
Figure 13: This figure, from the models from NPVY, shows the time the WD
have been cooling, tcool , at the moment that mass transfer starts.
Figure 14: This figure is taken from Tutukov and Yungelson (1996) [34]. The
birthrate of AM CVn stars as a function of the age of the secondaries is shown.
21
where κ is the proportionality factor in equation 27 (U = κT ). Equation 29 is
solved by:
µ
¶− 34
3
L=
,
(30)
t+C
ακ
where C is an integration constant. This solution is fitted to the cooling curves
of Hansen 1999 [5] by Nelemans (2004) [19] as:
log(L) = log(Lmax ) − 1.33log(t),
(31)
where L is given in units of L¯ and t in years. log(Lmax ) is fitted by3 :
½
9 − (0.9 − MWD )
for MWD > 0.5M¯ .
log(Lmax ) =
9.4 − 1.33(0.45 − MWD ) for MWD < 0.5M¯ .
Because carbon WDs stars their lifes warmer the helium WDs, the luminosity
of carbon WDs is cut off at log(L/L¯ )= 2 and for helium WDs at log(L/L¯ )=
0.5. From this fit it follows that:
µ ¶− 34
3
.
(32)
Lmax =
ακ
This is how the WDs cool if there is accretional heating (case 2 in section 3).
If the WDs heat up (case 1) then equation 29 gives:
κα − 3 ˙
L 4L
P − L = U˙ =
4
⇒
3
4 − 34
(P (t) − L) L 4 ,
L˙ = Lmax
3
(33)
where P (t) is the power of the accretional heat in units of L¯ that heats up the
core of the WD for a given time t. This is then integrated numerically to find
L as a function of age t for the accretor. Taking for the accretor:
½
1/2Pem (t) in the direct impact phase.
P (t) =
1/2Pbl (t) in the disk phase.
Where Pbl is the power released in the boundary layer. Pbl is calculated by:
µ
¶
Pbl = 1/2Pacc = 1/2M˙ φ(L1 ) − φ(acc) ,
(34)
where Pacc is total accretion power. We assume that half of it is released in the
disk and the other half is released in the boundary layer. φ(acc) is the potential
taken on the surface of the accretor.
For the donor the accretional heating is taken to be Pab , which is taken as the
irradiation from the impact spot in de direct impact phase and the irradiation
from the hot spot on the disk in the disk phase. The disk radius is taken to be
half the radius of the Roche lobe of the accretor.
3 no cristalization effects are taken into account, which would release a latent heat. This is
because that happens only after roughly 5 Gyr
22
Figure 15: This figure is taken from Nelemans et al. (2004) [19]. It presents the
cooling curves of white dwarfs.
4
Results
In the introduction I have explained what AM CVn stars are and why it is
important to look at them. I have made a model describing the evolution of
an AM CVn star. In the theoretical background I explained the formulae I
have used and the assumptions I have made. In this section the results will be
discussed.
In my program I have taken a few free parameters:
• The accretion can cause the WDs to warm up or not.
Case 1 : The WDs heat up as a result of accretion.
Case 2 : The WDs do not heat up as a result of accretion.
• The masses of the WDs at the onset of mass transfer. As explained in
section 3.3, I have looked at the evolution of these three systems:
1. M1 = 0.2 M¯ and M2 = 0.1 M¯ .
2. M1 = 0.5 M¯ and M2 = 0.25 M¯ .
3. M1 = 1.0 M¯ and M2 = 0.2 M¯ .
• The time tcool the WDs have been cooling at the moment the AM CVn
star is formed (see figure 14 and 13).
1. tcool = 106 yr.
2. tcool = 108 yr.
3. tcool = 109 yr.
• The spin of the accretor ω1 is taken to vary between 0 and ω1 max
1. ω1 = 0 s−1 .
2. ω1 = 1/2 ω1 max .
3. ω1 = ω1 max .
23
4.1
The effect of different white dwarf masses
The goal of my thesis is to find the secondary in an AM CVn star. So I want
to know in which systems they are relatively bright. Making plots for the three
systems with and without heating up of the WDs (case 1 and case 2) we get
6 plots: figure 16, 17, 18, 19, 20 and 21 (taking tcool = 108 yr and ω1 = 0
s−1 ). The luminosity of the accretor LWD1 , the luminosity of the donor LWD2 ,
the accretion power released at the impact spot Pem and the accretion power
received by the secondary Pab are plotted against time in the direct impact
phase. In the disk phase Pem is replaced by the power released at the hot spot
and Pab is then the power of the hot spot that falls on the donor. In the disk
phase the power released in the boundary layer Pbl is also plotted. The power
released in the accretion disk can be found in the same figures too. This is
because the power released in the disk (including the power released at the hot
spot on the disk) is taken equal to Pbl .
Looking at these plots some conclusions can be drawn for all plots: The
irradiation of the donor in the direct impact phase lasts only a short time,
because the impact spot moves behind the accretor as seen from the donor.
When that happens, Pab = 0 L¯ so log10 (Pab /L¯ ) = −∞. Using this, it is easy
to see in the plots when the disk phase starts, by looking at the sharp increase
in Pab (see figure 16). A lighter systems stays in the direct impact phase much
longer (figure 16 and 17), the third system (M1 = 1.0 M¯ and M2 = 0.2 M¯ ,
figure 20 and 21) even starts with a disk. It is also clear that even in the most
favourable scenario the heating up of the donor is negligible, but the heating up
of the accretor can be very important.
Looking at the first system (M1 = 0.2 M¯ and M2 = 0.1 M¯ , figure 16 and
17), the donor is always going to be less luminous than the accretor, because
both the donor and the accretor are taken to be made up out of helium. This
means that the heavier star (the accretor) will cool slowest. I still see two
possibilities how one might find the donor in such a system:
Until 106 years the irradiated part of the donor might be seen when the
impact spot is blocked from our view by the accretor. This part of the donor
could heat up locally and produce a different spectrum and different spectral
lines as the accretor.
After 108 years the donor is not a lot less luminous than the other components and it might be seen. But it is still going to be hard to distinguish the
spectrum of donor in the total spectrum or the system, because the spectrum
of the accretor will resemble that of the donor and it will be more luminous.
Looking at the second system (M1 = 0.5 M¯ and M2 = 0.25 M¯ , figure 18
and 19), we see that without heating up, the donor is brighter than the accretor
at all ages. This is because helium WDs cool slower than carbon WDs (The
WDs are carbon WDs if their mass is bigger than 0.5 M¯ ). After 108 years the
luminosity of the donor will also be comparable to that of the accretor even if
there is a moderate heating of the accretor. The chance of seeing a hot irradiated
surface (like in the first system) seems very small, because of the short lifespan
of the direct impact phase.
Looking at the third system (M1 = 1.0 M¯ and M2 = 0.2 M¯ ), figure 20 and
21), we see that without heating up the donor is not brighter than the accretor.
Although the donor is made of helium and thus cools more slowly, the large
mass of the accretor has the same effect on its cooling. Even if the accretor is
24
Figure 16: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.2 M¯ and M2 = 0.1 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
Figure 17: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 2, M1 = 0.2 M¯ and M2 = 0.1 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
25
Figure 18: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
Figure 19: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 2, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
26
Figure 20: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 1.0 M¯ and M2 = 0.2 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
Figure 21: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 2, M1 = 1.0 M¯ and M2 = 0.2 M¯ , tcool = 108 yr and ω1 = 0 s−1 .
27
heated up only a little, it is still going to be very hard to see a sign of the donor
in the spectrum of this system.
4.2
The effect of the cooling time
If the WDs are older when mass transfer starts they will be cooler. In order
to see this effect I have made a plot of the second system (M1 = 0.5 M¯ and
M2 = 0.25 M¯ ; the most promising system to find the donor). For that system
I have made plots for a cooling time of 106 years and 109 years. With and
without letting the stars heat up (case 1 and case 2) with ω1 = 0. See figure
22, 23, 24 and 25. The plots for tcool = 108 years are already plotted in figure
18 and 19
If the AM CVn star is born after 106 years of cooling of the WDs, then the
WDs are still luminous at the moment mass transfer starts. Because carbon
WDs start their lifes warmer the accretor will be the most luminous star. In
the case when the accretion only moderately heats up the WDs, the donor will
be more luminous after ∼ 107 years. This could happen even sooner if the
accretor is much older than the donor at the onset of mass transfer. If that
is also the case, then the donor will always be the most luminous star of the
system. These seem to be the most favourable circumstances to find the donor.
Note that the cooling curves of the donor in figure 22 and 23 are not accurate
before t = 107 years, because the intrinsic luminosity of helium WDs is cut off
at log10 (L/L¯ ) = −0.5 (see section 3.3).
If both, the donor and the accretor, are 109 years old at the onset of mass
transfer, then the best chances of finding the donor are in a very old system
t & 108 yr and only if the accretor does not heat up significantly.
4.3
The effect of the accretor’s spin
If the accretor is spun up, then the accretion flow will hit the accretor at a
relatively lower velocity producing less energy. In this we assume that the
energy released depends n the velocity difference between the stream and the
WD surface. Two systems are plotted, both for the second system, ((M1 = 0.5
M¯ and M2 = 0.25 M¯ ) with tcool = 108 years only with heating up (case 1)
for ω1 = 1/2ω1 max (figure 26) and ω1 = ω1 max (figure 27). The figure with
ω1 = 0 s−1 is already plotted (figure 18).
In this model a higher ω1 only has an effect as long as the system is in the
direct impact phase. It decreases the radiation produced at the impact spot
most, when the mass stream hits the donor at an angle of almost 90◦ with
respect to the normal to the surface of the accretor. So this effect is most
important when the impact spot is on the back side of the accretor as seen from
the donor. This means that it will hardly influence the amount of irradiation
on the donor. It also hardly seems to influence the heating of the accretor.
Therefore it can be concluded that the effect of spin of the accretor in finding
the donor is negligible.
28
Figure 22: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 106 yr and ω1 = 0 s−1 .
Figure 23: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 2, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 106 yr and ω1 = 0 s−1 .
29
Figure 24: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 109 yr and ω1 = 0 s−1 .
Figure 25: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 2, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 109 yr and ω1 = 0 s−1 .
30
Figure 26: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 108 yr and ω1 = 0.5ω1 max .
Figure 27: The power that the different components of an AM CVn system emit
as a function of time. LWD1 (solid), LWD2 (dashed), Pem (dash-dot-dot-dot),
Pab (dotted), and Pbl (dot-dash-dot-dash).
For: case 1, M1 = 0.5 M¯ and M2 = 0.25 M¯ , tcool = 108 yr and ω1 = ω1 max .
31
5
Discussion and conclusions
In order to find out when the donor can be seen in AM CVn stars I have made
a model for the evolution of this kind of system, to determine how the spectrum
of such a system changes during its evolution. I have looked at the effects of
different parameters: heating up of the WDs as a result of accretion, different
masses of the WDs, different ages of the WDs at the start of mass transfer and
at the effect of the spin of the accretor on the luminosity of the system. This
gives a large range of parameters in which an AM CVn star can be found. I
have made a model to predict out how and when the donor star can be best
found. There are two cases in which this happens.
In very young systems (t . 105 or 106 yr), where the system is in a direct
impact phase (see figure 6 for the possible mass ratios for this phase) and the
impact spot is still facing the donor. The X-rays produced at the impact spot
will irradiate the donor with ionizing radiation. This could change the spectrum
coming from this part of the donor and produce different spectral lines. If there
is a spectral line in which the donor dominates the spectrum, then the speed
of the donor can be found by looking at the red- and blue-shift of that line. It
is clear from figure 6 that the number of systems in a direct impact phase is
highly sensitive to the value of τs . A small value for τs would make it easier to
find the donor in this way. It is important to note that a large part of the time
that a system is in a phase of direct impact, the impact spot is not facing the
donor (see for instance figure 19). This will significantly decrease the chances
of finding the donor in this way. In this same time span the donor might also
be more luminous than the accretor if it is a very young star and the accretor
is an old one. The accretor can be cold if it has not had enough time to warm
up or if is not warming up effectively.
Another promising scenario to find the donor is in systems that contain a
carbon WD as an accretor and an helium WD as a donor. If the accretor is
also not much heavier than the donor (not more then 3 to 4 times) then the
donor will have a larger heat content than the accretor and thus cool more
slowly (see equation 27). Therefore a system with a large value for q would be
favourable for the detection of the donor in this way, but the existence of these
kind of systems depends strongly on the value of τs . A small value of τs would
also increase the chances of finding the donor in this way. In most cases the
spectrum of the donor will not be observable, because due to the mass transfer
the donor will have a similar composition in its surface as the accretor.
There are two aspects that have not been modelled, that I expect to have a
major influence on the luminosity of an AM CVn star. In old AM CVn systems
both of these effects would decrease the luminosity of the donor much further.
The first aspect is that for the cooling of both the accretor and the donor I have
taken their heat content to be constant. Since mass is lost the heat content will
in fact change, especially for the donor, whose mass relatively changes the most.
The second aspect is that I have not taken the effect of compressional heating
into account. As the masses of the WDs change, so does their radius. Since the
donor grows when mass is lost, it will cool. Likewise the accretor will heat up.
Some other aspects that might also play a roll, but have not been considered
in the model are: i) Other possible formation scenarios. If the donor is formed
out of a helium burning star, then the donor will still be warm when its core
burning stops. At the same time the accretor might also be warm as a result of
32
past accretion. Note that the mass radius relation for the donor will be different.
The system will start as an AM CVn star at a period of ∼ 10 minutes, which
is for most combinations of WD masses already in the disk phase. ii) The tidal
heating of the donor, deformation by the gravity of the accretor, might heat
it up. I expect that this effect is not very strong, because these forces will be
zero if the donor is synchronised to the orbit. iii) The rapid spin of the accretor
might smear out the impact spot. If the smearing is small it will move the
impact spot further out of the view of the secondary. In the most extreme case
the spin of the accretor could smear the impact spot over a narrow band on its
equator, moving it back into sight again for the donor. The spin of the accretor
would also decrease the light produced in the disk phase and thus lower the
luminosity released in the boundary layer.
33
6
Future work
The next step would be to add compressional heating and changing heat contents
to this model. The compressional energy can be estimated by using the fact that
for a star of mass M and radius R, with a constant density, the total internal
energy U is given by (see Prialnik (2000) [27] section 2.4):
U =−
3 GM 2
.
5 R
(35)
Together with the mass-radius relation (equation 6), this gives the total internal
energy of a WD for a given mass. Also a better estimate for the efficiency at
which the accretion heats up the core of the accretor is required. We must
calculate how far the accretion stream penetrates the atmosphere of the accretor.
This can be estimated by looking at where the atmospheric pressure is equal
to the ram pressure of the stream. Also the thickness of the atmosphere and
the efficiency at which heat is conducted through the atmosphere need to be
calculated. All this is required to check if the donor could in any scenario
outshine the accretor after ∼ 108 yr.
In order to find out how the donor can be seen in the direct impact phase,
it is important to find out what kind of emission lines are produced in a helium
atmosphere which gets irradiated by X-rays. The absorbtion for helium of Xrays should be calculated. This will depend on the angle at wich the irradiation
hits the atmosphere. To calculate this, the spatial plot of my program (see
figure 10) can be used. In order to estimate what kind of X-rays are produced
at the impact spot, observations of polars can be used, or a theoretical model
of the impact spot can be made.
34
References
[1] Abbott, Timothy M. C.; Robinson, Edward L.; Hill, Gary J.; Haswell, Carole A., 1992, ”Two unusual cataclysmic variables at high Galactic latitude
- CP Eridani and AL Comae”, ApJ, 399, 680
[2] Cropper M.; Harrop-Allin M. K.; Mason K. O., Mittaz J. P. D.; Potter
S. B.; Tamsay G., 1998, ”RX J1914.4+2456 - The first double-degenerate
polar?”, MNRAS, 293, L57
[3] Duquennoy, A.; Mayor, M. 1991, ”Multiplicity among solar-type stars in
the solar neighbourhood. II - Distribution of the orbital elements in an
unbiased sample”, A & A, 248, 485
[4] Faulkner, John; Flannery, Brian P.; Warner, Brian, 1972, ”UltrashortPeriod Binaries. II. HZ 29 (=AM CVn): a Double-White Semidetached
Postcataclysmic Nova?” ApJ, 175, 79
[5] Hansen, Brad M. S., 1999, ”Cooling Models for Old White Dwarfs”, ApJ,
520, 680
[6] Hirose, Masahito; Osaki, Yoji, 1990, ”Hydrodynamic simulations of accretion disks in cataclysmic variables - Superhump phenomenon in SU UMa
stars”, PASJ, 42, 135
[7] Hirose, Masahito; Osaki, Yoji, 1993, ”Superhump periods in SU Ursae
Majoris stars: Eigenfrequency of the eccentric mode of an accretion disk”,
PASJ, 45, 595
[8] Humason, M., Zwicky, F., 1947, ”A Search for Faint Blue Stars.”, ApJ,
105, 85
[9] Iben, Icko, Jr.; Livio, Mario, 1993, ”Common envelopes in binary star
evolution”, PASP, 105, 1373
[10] Israel G. L. et al., 2002, ”RX J0806.3+1527: A double degenerate binary
with the shortest known orbital period (321s)”, A&A, 386, L13
[11] Marsh, T. R.; Horne, Keith; Rosen, Simon, 1991, ”Evidence for CNO processed material in the accretion disk of GP Comae”, ApJ, 366, 535
[12] Marsh, T. R.; Steeghs, D., 2002, ”V407 Vul: a direct impact accretor”,
MNRAS, 331, L7-L11
[13] Marsh, T. R.; Nelemans, G.; Steeghs, D., 2004, ”Mass transfer between
double white dwarfs”, MNRAS, 350, 113-128
[14] Milgrom, M.; Salpeter, E. E., 1974 ”Models for X-ray illuminated atmospheres”, ApJ, 196, 583-588.
[15] Nather, R. E.; Robinson, E. L.; Stover, R. J., 1981, ”The twin-degenerate
interacting binary G61-29”, ApJ, 244, 269
[16] Nelemans, G.; Yungelson, L. R.; Portegies Zwart, S. F.; Verbunt, F., 2001,
”Population synthesis for double white dwarfs I. Close detached systems”,
A&A, 365, 491-507
35
[17] Nelemans, G.; Yungelson, L. R.; Portegies Zwart, S. F., 2001, ”The gravitational wave signal from the Galactic disk population of binaries containing
two compact objects”, A&A, 375, 890-898
[18] Nelemans, G.; Yungelson, L. R.; Portegies Zwart, S.F. 2003, ”Short-period
AM CVn systems as optical, X-ray and gravitational-wave sources”, MNRAS, 349, 181-192
[19] Nelemans, G.; Yungelson, L. R.; Portegies Zwart, S. F., 2004, ”Shortperiod AM CVn systems as optical, X-ray and gravitational-wave sources”,
MNRAS, 349, 181-192
[20] Landau, Lev Davidovich; Lifshitz, E. M., 1971, ”The classical theory of
fields” (Book)
[21] O’Donoghue, D.; Kilkenny, D.; Chen, A.; Stobie, R. S.; Koen, C.; Warner,
B.; Lawson, W. A., 1994, ”EC:15330-1403 and the Am-Canum Stars”,
MNRAS, 271, 910
[22] Paczy´
nski, B. (1967), ”Gravitational Waves and the Evolution of Close
Binaries” Acta Astron., 17, 287
[23] Patterson, J. et al., 1997, ”Superhumps in Cataclysmic Binaries. XII. CR
Bootis, a Helium Dwarf Nova”, PASP, 109, 1100 (P1997)
[24] Patterson, Joseph; Walker, Stan; Kemp, Jonathan; O’Donoghue, Darragh;
Bos, Marc; Stubbings, Rod, 2000, ”V803 Centauri: A Helium-rich Dwarf
Nova”, PASP, 112, 625 (P2000)
[25] Patterson J.; Fried R. E.; Rea R.; Kemp J.; Espaillat C. et al.; 2002, ”Superhumps in Cataclysmic Binaries. XXI. HP Librae (=EC 15330-1403)”,
PASP, 114, 65
[26] Podsiadlowski, PH; Han, Z; Rappaport S., 2002, ”Cataclysmic variables
with evolved secondaries and the progenitors of AM CVn stars”, MNRAS,
340, 1214-1228
[27] Prialnik, D. 2000, ”Stellar Structure and Evolution” (Book), Cambridge
University Press
[28] Ramsay G.; Hakala P.; Cropper M., 2002, ”RX J0806+15: the shortest
period binary?”, MNRAS, 332, L7
[29] Ruiz, Maria Teresa; Rojo, Patricio M.; Garay, Guido; Maza, Jose, 2001,
”CE 315: A New Interacting Double-Degenerate Binary Star”, ApJ, 552,
679
[30] Savonije, G. J.; de Kool, M.; van den Heuvel, E. P. J., 1986, ”The minimum
orbital period for ultra-compact binaries with helium burning secondaries”,
A&A, 155, 51-57
[31] Skillman D. R. et al., 1999, ”Superhumps in Cataclysmic Binaries. XVII.
AM Canum Venaticorum”’, PASP, 111, 1281
[32] Smak, J. (1967), ”Light Variability of HZ 29”, A&A, Vol 17, No 3
36
[33] Solheim J. -E et al. 1998, ”Whole Earth Telescope observations of AM
Canum Venaticorum - discoseismology at last”, A&A, 332, 939
[34] Tutukov, A.; Yungelson, L, 1999, ”Double-degenerate semidetached binaries with helium secondaries: cataclysmic variables, supersoft X-ray
sources, supernovae and accretion-induced collapses”, MNRAS 280, 10351045
[35] Warner B., 1995, ”Cataclysmic variable stars” (Book), Cambridge Astrophysics Series
[36] Warner B.; Woudt P. A., 2002,”KUV 01584-0939: A Helium-transferring
Cataclysmic Variable with an Orbital Period of 10 Minutes”, PASP, 114,
129
[37] Wood, M. A.; Winget, D. E.; Nather, R. E.; Hessman, Frederic V.; Liebert,
James; Kurtz, Donald W.; Wesemael, F.; Wegner, G., 1987, ”The exotic
helium variable PG 1346 + 082”, ApJ, 313, 757
[38] Weidemann V.; Koester D., 1984, ”Mass distribution of DA white dwarfs
and atmospheric parameters of ZZ Ceti stars”, A&A, 132, 195-202
[39] Wood, Matt A.; Casey, Matthew J.; Garnavich, Peter M.; Haag, Brooke,
2002, ”Superhumps in the helium dwarf nova KL Draconis”, MNRAS, 328,
159
[40] Woudt, P. A.; Warner, B., 2001 ”High-speed photometry of faint cataclysmic variables - I. V359 Cen, XZ Eri, HY Lup, V351 Pup, V630 Sgr,
YY Tel, CQ Vel and CE-315” MNRAS, 328, 159
[41] Woudt, A.; Warner, B. 2003, ”The new AM CVn star in Hydra”, MNRAS,
345, 1266-1270
37
7
Appendix
In this section the c-program is included with which the calculations have been
done:
/********************************************************************************************************/
/* A program that calculates and plots graphicly the impactpoint of a balistic stream of mass, */
/* that is transferred in an AM CVn star. Further it calaculates the amount of radiation of this */
/* impact point, that is absorbed on the secondary. (the irradiated surface is also shown graphicly). */
/* M1 and M2 are taken as free variables. Coordinate system: positive x to the right, positive y */
/* up, positive z out of paper, zero in CM, M2 on negative x-axis, M1 on positive x-axis. */
/* */
/* This can be done sevaral times for different values of M2 which are characterize by values of m. */
/* The evolution of an AM CVn system is calculated for given masses of M1 and M2. The luminosity */
/* of the different components is calculated and plotted against the time. This can be done with and */
/* without heating up of the WDs. */
/********************************************************************************************************/
#include "cpgplot.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi
3.1415926536 // Rarely a program uses a more accurate value of pi.
()
#define e
2.7182818
// this one doesn’t need a discription, right? ()
#define c
2.998e8
// velocity of light through vacuum (m/s)
#define G
6.673e-11
// Gravitational constant (m^3 kg^-1 s^-2)
#define Msol 1.989e30
// Solar mass (kg)
#define Rsol 6.9598e8
// Solar radius (m)
#define Lsol 3.8515e26
// Solar luminosity (J s^-1)
#define Mch
1.44
// Chandrasekhar mass (Msol)
#define year 3.156e7
// Number of seconds per year (s)
#define Kb
1.38e-23
// Boltsmanns constant (J/K)
#define sigma 5.67e-8
// Stefan-Boltsmanns constant (J m^-2 s-1 K-4)
#define m_H
1.66e-27
// Atomic mass unit (kg)
#define Mp
0.00057
// a constant used in Eggletons mass-radius relation. (Msol)
/********************************************************************************/
/* There are a few free parameters: (they are set in main() ) */
/* */
/* -M1, M2 at the birth of the AM CVn system (Mch>M1>M2>0) */
/* -w1 can be chanced (from 0 to w1_max) by variing factorw1 between 0 and 1 */
/* -factor, can be chanced to increase accuracy (=1.01 accurate; =1.5fast */
/* -autoscale can be chanced to automatically scale plotspatial. */
/********************************************************************************/
char deviceSpatial[] = "10/xserve"; // "10/xserve" or "/ps" or "?"; the device in which the spatial plot is made. ()
char deviceXvsY[]
= "11/xserve"; // "11/xserve" or "/ps" or "?"; the device in which X[m] is plotted vs Y[m]. ()
float tbegin; // the time WD2 can cool before masstransfer starts (year)
float t_ev; // the evolution time; a running parameter (year)
float tdisk; // the time at which the diskphase starts (year)
float t3[4001]; // runs with t_ev; is used to get a better resolution (year)
float timestep; // how much time passes in one step=4e5 (year)
float q;
// mass ratio ()
float M1, M2;
// primary mass, secondary mass (Msol)
float M1initial, M2initial;
// primary mass, secondary mass (Msol)
float M2dot; // The masstransferrate; the time derivative of M2 (<0) (Msol year^-1)
float R1; // radius of primary (m)
float R2; // radius of secondary (m)
float a; // binary seraration (m)
float Porb; // orbital period (s)
float w; // omega angularvelocity of the orbit (s^-1)
float R_disk; // radius of the accretiondisk. (only present in the second part of the program (m)
float w1; // omega angularvelocity of WD1 around it’s own axis (is set in getDynamics) (s^-1)
float w1_max; // maximum anglarvelocity of WD1 around it’s own axis. (s^-1)
float factorw1; // w1=factorw1*w1_max
float L1[2];
// inner Lagrange point distance vector to centre of mass (m)
float KL1[2];
// inner Lagrange point velocity vector to centre of mass (m)
float a1[2], a2[2];
// primary, secondary distance vectors to centre of mass
(m)
float K1[2], K2[2];
// binary velocity vectors (m/s)
float xval[1000]; // the x values of the balistic stream [steps!!!!] (m)
float yval[1000]; // the y values of the balistic stream (m)
float roche1x[2][1001]; // rochelobepoints of the primary
[x or y][1+rocheBins] ()
float roche2x[2][1001]; // rockelobepoints of the secondary [x or y][1+rocheBins] ()
float impact[2]; // impact position[x,y]; either on WD1 or the hotspot on the disk (m,m)
float phi; // place of impact on WD1; define by angle with resect to horizontal line (degrees)
float impactangle; // impactangle of the balisticstream and the normal of the primary (degrees)
float irradiatedsurface[2][1001]; // point which mark the surface being radiated [x/y][rochebins]
(m,m)
float partofcircle; // Taking the projection of WD2 as a circle, this is the part of the circle...
// ...that gets irradiated. 1 => completely; -1 => not at all. ()
float openingangle; // How much of the 180 degrees irradiation hits WD2
(degrees)
float Pout[1001]; // Outgoing Power emitted at the impactspot [calculated at getdynamics] (J/s)
float Pin[1001]; // Incomming Power received at WD2 [calculated at getdynamics] (J/s)
float P_bl[1001]; // The energy releades at the boundrylayer [calculated at getdynamics] (J/s)
float X1[1001],X2[1001],X3[1001],X4[1001],X5[1001];// X is plotted vs Y for different values of m.... (?)
float Y1[1001],Y2[1001],Y3[1001],Y4[1001],Y5[1001];// ...this way the memory is only reserved once. (?)
float factor; // timestep*=factor; So it sets the speed; factor=1.3 fast =1.01 accurate. () ()
int
heating; // when the WDs heat up then heating=1; if not heating =0. ()
int
autoscale=1; // scaling of plotspatial; =1=> zoom to fit screen; =0=>manual scaling. ()
int
n1,n5; // count number of point to plot Y vs X. should be less the 1001. ()
int
rocheBins = 100; // roche lobe phase bins. The number of point of a rochelobe = 100 ()
int
plotSpatial; // plot IDs ()
int
plotXvsY; // plot IDs ()
int
noimpact=1; // is set to 0, when there is impact. (In getBallistics) ()
int
irradiation; // Does WD2 see the impactspot; 2=completely, 1=partially, 0=not at all ()
38
int
long
long
p; // counts rochebins of WD2, which are irradiated ()
steps = 1000;
// iterative steps in several calculations =1000 ()
steps2 = 1000;
// iterative steps in ballistics calculation (should be changed there)=step1 ()
/************************************************************************/
/* This function uses Eggleton’s zero temperature mass-radius relation */
/* (formula 24 of Mass Transfer between Double White Dwarfs to get R */
/* for all values of M. It’s called at getDynamics() */
/************************************************************************/
float getRadius(float M)
{
float R; /*radius (m)*/
R=0.0114*Rsol*sqrt(pow((M/Mch),-0.6667)-pow((M/Mch),0.6667))*pow((1+3.5*pow((M/Mp),-0.6667)+(Mp/M)),-0.6667);
if(M>Mch)
printf("Mass > Mch");
return(R);
}
/************************************************************************************************/
/*This function calculates the basic values of an AMCVn system. It’s only imput are M1 and M2 */
/************************************************************************************************/
void getDynamics(void)
{
R1=getRadius(M1);
R2=getRadius(M2);
q=M2/M1;
a=R2*(0.6*pow(q,0.667)+log(1+pow(q,0.333)))/(0.49*pow(q,0.667)); /*2.5c Warner*/
Porb=pow(((4*pi*pi*a*a*a)/(G*Msol*(M1+M2))),0.5); /*Keplers law 2.1a Warner*/
w = 2*pi/Porb;
w1_max=sqrt(G*M1*Msol/(R1*R1*R1));
w1=factorw1*w1_max;
a1[0] = a*M2/(M1+M2);
a1[1] = 0.;
a2[0] = -a*M1/(M1+M2);
a2[1] = 0.;
K1[0] = 0.;
K1[1] = -2*pi*a1[0]/Porb;
K2[0] = 0.;
K2[1] = -2*pi*a2[0]/Porb;
L1[0] = -a / (1.0015 + pow(q,0.4056)) + a1[0]; /* Warner 2.4c accurate to 1% in 0.04<q<1 */
L1[1] = 0.;
KL1[0] = 0.;
KL1[1] = -2*pi*L1[0]/Porb;
// R_disk=pow(((-a2[0]*(-a2[0]-R2)-a1[0]*a1[0])*w),2)/(G*Msol*(M1+M2)); /*calculates the radius of a possible accretiondisk*/
R_disk=0.5*(a1[0]-L1[0]);
// printf(" a1(0)=%4.2le a2(0)=%4.2le Porb=%4.2le a=%4.2le \n R1=%4.2le R2=%4.2le M1=%g M2=%g w=%4.2le w1=%4.2le w1_max=%4.2le\n",
// a1[0], a2[0], Porb, a, R1, R2, M1, M2, w, w1, w1_max);
}
/************************************************************************/
/* This function calculates Zeta_2 for a given value of mass (M). */
/* Zeta_2 = d log(R2)/ d log(M2) */
/* It is called for in getAngularmomentum(). */
/************************************************************************/
float getZeta_2 (float M)
{
float Zetafirst; /*dummy; the first part of the expression of Zeta_2*/
float Zetasecond; /*dummy; the second part of the expression of Zeta_2*/
float Zeta_2; /*= d log(R2)/ d log(M2) */
Zetafirst = -0.333*(1/(pow((M/Mch),-0.667)-pow((M/Mch),0.667)))*pow((M/Mch),-0.667);
Zetasecond = (0.667/(1+3.5*pow((M/Mp),-0.667)+(Mp/M)))*(2.333*pow((M/Mp),-0.667)+(Mp/M));
Zeta_2=Zetafirst+Zetasecond;
// Y[m]=Zetafirst;
// printf("zetafirst=%4.2le zetasecond=%4.2le \n",Zetafirst, Zetasecond);
return(Zeta_2);
}
/************************************************************************/
/* This function calculates the masstransferrate in an AMCVn system. */
/************************************************************************/
void getAngularmomentum()
{
double J_orb; /*The angular momentum stored in the orbit of the binary*/ // (kg m^2 s^-1)
double Jdot_gr; /*The change of J_orb per second due to gravitational waves*/ // (kg m^2 s^-2)
float adot; /*The timedirivative of a*/ // (
m
s^-1)
float Zeta_2; /*Dummy variable; See Marsh et ai 2004 for a definition*/
float Zeta_rL; /*Dummy variable; See Marsh et ai 2004 for a definition=0.333*/
float k=0.2; /*The moment of inertia of the WD = 0.2 in equation 6 of Nelemans 2004 ai*/ //(?)
J_orb=sqrt(G*a/((M1+M2)*Msol))*M1*Msol*M2*Msol;
Jdot_gr=-6.4*((G*G*G)/(pow(c,5)))*(M1*Msol*M2*Msol*(M1+M2)*Msol/(pow(a,4)))*J_orb;
Zeta_2 = getZeta_2(M2);
Zeta_rL=0.333; /*0.333*(1+q)*2*log(1+pow(q,0.333)-pow(q,0.333)/(1+pow(q,0.333)))/(0.6*pow(q,0.667)+log(1+pow(q,0.333))); in Marsh 2004 (almost equal to 0.3)*/
M2dot=year*M2*Jdot_gr/(J_orb*(1+0.5*Zeta_2-0.5*Zeta_rL-q)); /*be careful with the units of year and Msol*/
// Y[m]=log10(-M2dot);
// printf("Zeta_2=%4.2lf Zeta_rL=%4.2lf J_orb=%4.2le Jdot_gr=%4.2le M2dot=%4.2le \n", Zeta_2, Zeta_rL, J_orb, Jdot_gr, M2dot);
}
/********************************************************************************************************/
/* This function calculates the potential for given x and y coordinates. It’s called at getRocheLobes() */
/********************************************************************************************************/
float getRochePotential(float x, float y)
{
return (-G*M1*Msol / sqrt(pow(x-a1[0],2)+pow(y-a1[1],2))
-G*M2*Msol / sqrt(pow(x-a2[0],2)+pow(y-a2[1],2))
-0.5*pow(w,2)*(x*x + y*y));
}
/****************************************************************************************/
/* This function calculates the coordinates of the Roche Lobe. Method: move outward */
39
/* radially from the primary and the secondary as long as the potential is less than */
/* the potential in L1. Then repeat at a different angle, theta. */
/****************************************************************************************/
void getRocheLobes(void)
{
float potentialL1; /*the potential in L1
*/
float potential; /*the potential in a point
*/
float theta; /*going round M1 and M2 by changing theta
*/
float r,x,y; /*distance to the center of M1 or M2; x,y:Coordinates
*/
int n; /*a running variable to sweet around M1 and M2. It runs from 0 to rocheBins */
potentialL1 = getRochePotential(L1[0],L1[1]);
theta = 0.;
for (n=1; n<rocheBins; n++)
{
r = 0.9*R2;
theta = 2*pi*n/rocheBins;
do
{
r += fabs(a1[0]-L1[0])/steps; /*fabs=absolute value*/
x = a1[0] - r*cos(theta); /* starting at L1 viewed from primary, going counter-clockwise */
y = a1[1] - r*sin(theta);
potential = getRochePotential(x,y);
//
printf("x-a0=%4.2le y=%4.2le r=%4.2le potential=%4.2le potentialL1=%4.2le \n "
//
, x-a1[0], y, r, potential, potentialL1);
}
while (potential <= potentialL1);
roche1x[0][n] = x;
roche1x[1][n] = y;
//
printf("y1=%4.2le n=%d \n",roche1x[1][n], n);
}
roche1x[0][0] = L1[0];
roche1x[1][0] = L1[1];
roche1x[0][rocheBins] = L1[0]; /*coordinates of the Roche Lobe, they .... */
roche1x[1][rocheBins] = L1[1]; /*...are plotted in plotrochelobe(). */
theta = 0.; /*for the secondary the same...*/
for (n=1; n<=rocheBins; n++)
{
r = 0.9*R2; /*the rochelobe of the primary can’t be smaller than the secondary */
theta = 2*pi*n/rocheBins;
do
{
r += fabs(a2[0]-L1[0])/rocheBins;
x = a2[0] + r*cos(theta); /* starting at L1 viewed from the secondary, */
y = a2[1] + r*sin(theta); /* going counter-clockwise
*/
potential = getRochePotential(x,y);
}
while (potential < potentialL1);
roche2x[0][n] = x;
roche2x[1][n] = y;
roche2x[0][rocheBins] = L1[0];
roche2x[1][rocheBins] = L1[1];
//
printf("y2=%4.2le n=%d \n",roche2x[1][n], n);
}
roche2x[0][0] = L1[0];
roche2x[1][0] = L1[1];
roche2x[0][rocheBins] = L1[0];
roche2x[1][rocheBins] = L1[1];
}
/**************************************************************************/
/* This function calculates the coordinates of the balistic stream. This */
/* is a free fall trajectory from L1 until either WD1 or the disk is hit */
/* by the stream. This is determined by impactradius.
*/
/**************************************************************************/
void getBallistics(float impactradius) /*the ballistic stream runs until it closer to WD1 than impactradius*/
{
float t = 0.; /*a time variable (s)*/
float dt; /*timestep*/
float mx[2], mv[2], ma[2]; /*co rotating -position, -velocity and -acceleration in SI*/
float distance; /*Distance to M1 or M2 */
int n=0; /*counts the points of the trajectory. There are step2 point, with a maximum of step1*/
mx[0] = L1[0]+0.01*fabs(L1[0]); /* start test mass just outside of L1 point to start falling*/
mx[1] = 0;
mv[0] = 0;
mv[1] = 0;
noimpact=1;
for(n=0;n<steps;n++)
{
ma[0] = G*M1*Msol*(a1[0]-mx[0]) * pow(pow(a1[0]-mx[0],2) + pow(mx[1],2),-1.5) /* gravity from M1
*/
+ G*M2*Msol*(a2[0]-mx[0]) * pow(pow(a2[0]-mx[0],2) + pow(mx[1],2),-1.5) /* gravity from M2
*/
+ mx[0] * pow(w,2) /* Centrifugal force */
+ 2*w*mv[1];
/* Coriolis force
*/
ma[1] = G*M1*Msol*(-mx[1]) * pow(pow(a1[0]-mx[0],2) + pow(mx[1],2),-1.5)
+ G*M2*Msol*(-mx[1]) * pow(pow(a2[0]-mx[0],2) + pow(mx[1],2),-1.5)
+ mx[1] * pow(w,2)
+ -2*w*mv[0];
dt = 2*sqrt(pow(a/steps,2)) / sqrt((pow(mv[0],2) + pow(mv[1],2))); /* time it takes with previous speed to cross 2*a/steps */
if (dt > 2*Porb/steps) /* (dt>2*Porb/steps)*/
{ /* Dont change these 2’s. The maximum length of the ballistic stream... */
dt = 2*Porb/steps; /* is set in such a way, that is will hit WD1 in the first passage... */
} /* else not at all. */
mv[0] += ma[0] * dt;
mv[1] += ma[1] * dt;
mx[0] += mv[0] * dt;
mx[1] += mv[1] * dt;
t += dt;
xval[n] = mx[0];
40
//
yval[n] = mx[1];
distance=sqrt(pow((mx[0]-a1[0]),2)+pow((mx[1]-a1[1]),2));
printf("xval =%4.2le and yval=%4.2le and distance=%4.2le R1=%4.2le\n", xval[n], yval[n],distance, R1);
if(distance<impactradius)
{
steps2=n;
impact[0]=(mx[0]-a1[0]);
impact[1]=(mx[1]-a1[1]);
noimpact=0;
break;
}
}
phi=180/pi*atan(impact[1]/impact[0]); /*place of impact on WD1; define by angle phi with resect to horizontal line*/
if(phi<0) /*So -pi/2<atan<=pi/2 !!!*/
{
phi+=180;
}
}
/****************************************************************************************/
/* Function calculates the coordinates of the donor that are being irradiated. It also */
/* calculates the % of radiation emitted at the impactspot, that is being absorbed by */
/* the donor. Use Plotconstructionlines() and appended figure to how the angles work. */
/* Note: *counterclockwise is positive for all angles! */
/*
*angle=0 means horizontal. */
/* */
/* Fist it is checked if the impactspot "sees" the secondary. This done with alpha1, */
/* alpha2 and beta. Then the openingangle is found between gammamin and gammamax. */
/* If there is irradiation, than it is calculated upto what height (y or placeoncircle) */
/* WD2 is being irradiated. */
/****************************************************************************************/
void getIrradiatedsurface()
{
float directionnormal[2]; /* the x- and y-components of a vector between a1[] and impact[]
*/
float alpha1; //(degrees) /* the angle between the horizontal and the lineA1 between impact and the upperpart of WD1 */
float alpha2; //(degrees) /* the angle between the horizontal and the lineA2 between impact and the lower part of WD1 */
float alpha; //(degrees) /* the complete openingangle for WD2: alpha=alpha2-alpha1
*/
float beta;
//(degrees) /* the angle between the horizontal and the lineB. lineB goes through impact and is....
*/
/* ...perpendicular with the normal. It’s the horizon for a person on WD1 at impact[]
*/
float gamma; //(degrees) /* the angles between the horizontal and the lines between impact and the Roch lobe of WD2 */
float gammamax;//(degrees) /* angle between the horizontal and the line C1 between impact and the lower part of the rochelobe of WD2 */
float gammamin;//(degrees) /* angle between the horizontal and the line C2 between impact and the upper part of the rochelobe of WD2 */
float placeoncircle; //() /* runs from 1 to -1; the irradiation hits a unit-circle at this y-coordinate
*/
float dummy; // () /* this dummy is used to solve an integral
*/
int n; /* counts the rocheBins
*/
int min, max; /* impact[][min]/impact[][max] is the minimum/maximum irradiated point on WD2.
*/
directionnormal[0]=impact[0]/(sqrt(impact[0]*impact[0]+impact[1]*impact[1]));
directionnormal[1]=impact[1]/(sqrt(impact[0]*impact[0]+impact[1]*impact[1]));
alpha1=(180/pi)*atan(((a2[1]-R2*directionnormal[1])-(impact[1]+a1[1]) )/( (a2[0]-R2*directionnormal[0])-(impact[0]+a1[0]) ) );
alpha2=(180/pi)*atan(((a2[1]+R2*directionnormal[1])-(impact[1]+a1[1]) )/( (a2[0]+R2*directionnormal[0])-(impact[0]+a1[0]) ) );
beta= (180/pi)*atan(directionnormal[0]/(-1*directionnormal[1]));
/*first see what kind of irradiation there is */
if(beta<alpha1)
{
irradiation=2; /*full irradiation from impactspot*/
}
else
{
if(beta<alpha2)
{
irradiation=1; /*partial irradiation from impactspot*/
}
else
{
irradiation=0; /*no irradiation from impactspot */
}
}
// printf("alpha1=%4.2f alpha2=%4.2f beta=%4.2f irradiation=%d \n",alpha1, alpha2, beta, irradiation);
/*then find the points of roche2x[][] which are irradiated and at Plotirradiatedsurface() a polygon is drawn between them*/
if((irradiation==1)||(irradiation==2))
{
p=0;
gammamax=-90;
gammamin=90;
min=0;
max=0;
for(n=0;n<rocheBins;n++)
{
gamma=(180/pi)*atan(((impact[1]+a1[1])-roche2x[1][n])/((impact[0]+a1[0])-roche2x[0][n]));
//
printf("n=%d gamma=%4.2f beta=%4.2f switch=%d \n",n, gamma, beta, (roche2x[0][n]>a2[0]));
if((gamma>beta)&&(roche2x[0][n]>a2[0])) /*gamma>beta => WD1 is not blocking the view AND roche2x[0][n]>a2[0] => WD2 is not blocking the view*/
{
if(gammamax<gamma)
{
gammamax=gamma;
max=n;
}
if(gammamin>gamma)
{
gammamin=gamma;
min=n;
}
irradiatedsurface[0][p]=roche2x[0][n];
irradiatedsurface[1][p]=roche2x[1][n];
/*if you want to see which points are found; you let the lines below be drawn, remember to move this function below Plotrochelobes!
*/
//
cpgsls(1);
41
//
//
//
//
//
cpgmove((impact[0]+a1[0]), (impact[1]+a1[1])); /*draws upper line ...
*/
cpgdraw(roche2x[0][min], roche2x[1][min]); /*... to mark impact area: line C1 */
cpgmove((impact[0]+a1[0]), (impact[1]+a1[1])); /*draws lower line ...
*/
cpgdraw(roche2x[0][max], roche2x[1][max]); /*... to mark impact area: line C2 */
printf("n=%d irradiatedsurface[0][%d]=%4.2le irradiatedsurface[0][%d]=%4.2le \n", n, p,irradiatedsurface[0][p],p, irradiatedsurface[1][p]);
p++;
}
}
}
else
{
/*no irradiation from impactspot */
//
printf("no irradiation \n");
gammamin=0;
gammamax=0;
}
openingangle=gammamax-gammamin;
if(openingangle==-180)
{
openingangle=0; /*to correct the small difference between the first and second approach; R2 is not exactly R_{Rochelobe}*/
} /*because the radius of the rochelobe of R2; it’s not exactly circular.*/
// printf("openingangle=%4.2f gammamin=%4.2f gammamax=%4.2f \n", openingangle, gammamin, gammamax);
alpha=gammamax-alpha1;
if(openingangle>0)
{
placeoncircle= (tan((pi/180)*(openingangle-0.5*alpha)))/tan((pi/180)*0.5*alpha); /*runs from 1 to -1*/
if(placeoncircle>1)
{
placeoncircle=1; /*to correct for small mistakes due to rounding (again because R2 is nog R_{Rochelobe}), they’d blow up later.*/
}
//
printf("openingangle=%4.2f gammamax=%4.2f alpha1=%4.1f beta=%4.1f alpha=%4.2f placeoncircle=%4.2f \n \n",
//
openingangle, gammamax, alpha1, beta, alpha, placeoncircle);
}
else
{
placeoncircle=-1;
}
partofcircle=0;
for(dummy=-0.999; dummy<=placeoncircle; dummy+=0.003) /*this could perhaps be done analytically, but this works fine*/
{
partofcircle+=(2/pi)*sqrt(1-dummy*dummy)*0.003; /*we’re integrating a function that’s half a circle, so we do x2 and divide...*/
}
/*... by the area of an entire unit-circle:pi*/
// printf("placeoncircle=%4.2f partofcircle=%4.2f openingangle=%4.2f \n \n",placeoncircle, partofcircle, openingangle);
}
/****************************************************************************************/
/* This function calculates the impactangle between the balistic stream and the normal */
/* of the surface of WD1. Note: that if the impact is on the upper half of WD1 things */
/* will go wrong! */
/* It also calculates some variables connected the impact, such as what percentage of
*/
/* the radiation emitted at the impactspot falls onto WD2. (=Pin/Pout) */
/****************************************************************************************/
void getImpactdynamics()
{
float psi; //(degrees)/* angle of impact not relative to the primary */
float v_abs, v_abs_r, v_abs_phi; //(m/s)/* absolute velocities in CM-frame*/
float v_rel, v_rel_r, v_rel_phi; //(m/s)/* relative velocities to rotation WD1*/
float Mdot; //(Msol/year)/*should go to getDynamics; around 1e-8 Msol*/
float distance; //(m)
/*distance between the center of WD2 and the impactspot*/
float irradiatedpart; //() /*what part of the 4pi degrees of radiation fall on WD2*/
float Tout; //(K) /*T_max of the blackbody radiation assumed to be coming from the impactpoint*/
float nu_max, lambda_max; //(Hz), (m)/*The maximum of a blackbody spectrum for resp. the frequency and the wavelength*/
float albedo=0.4; //(), /*The X-ray albedo ~0.4*/
Mdot=-M2dot;
psi=180/pi*atan((yval[steps2]-yval[steps2-1])/(xval[steps2]-xval[steps2-1])); /*be careful with pi/2<atan<pi/2 */
impactangle=phi-psi;
// printf("psi=%4.2f phi[m]=%4.2f, impactangle=%4.2f \n", psi, phi[m], impactangle);
v_abs
= sqrt(2*(getRochePotential(L1[0],L1[1])-getRochePotential((impact[0]+a1[0]),(impact[1])+a1[1]))); /*U=K; energy conservation */
v_abs_r
= v_abs*cos(impactangle*pi/180);
v_abs_phi
= v_abs*sin(impactangle*pi/180);
v_rel_r
= v_abs_r;
v_rel_phi
= v_abs_phi-w1*R1;
v_rel
= sqrt(v_rel_r*v_rel_r+v_rel_phi*v_rel_phi);
// printf("m=%d v_abs=%4.2le v_abs_r=%4.2le v_abs_phi=%4.2le v_rel=%4.2le v_rel_r=%4.2le v_rel_phi=%4.2le, L1=%4.2le,impact=%4.2le \n",
// m, v_abs, v_abs_r, v_abs_phi, v_rel, v_rel_r, v_rel_phi, getRochePotential(L1[0],L1[1]), getRochePotential(impact[0],impact[1]));
// printf("w1*R1=%4.2le impact[0]=%4.2le, impact[1]=%4.2le \n", w1*R1, impact[0], impact[1]);
Pout[n1]=0.5*Mdot*(Msol/year)*v_rel*v_rel;
Tout=(2*m_H)*v_rel*v_rel/(3*Kb); /*1/2 m v^2 = 3/2 K T; Helium mass=2*m_H*/
lambda_max=0.002898/Tout;
nu_max=c/lambda_max;
distance=sqrt(pow((impact[0]-a2[0]),2)-pow((impact[1]-a2[1]),2));
irradiatedpart=partofcircle*R2*R2/(4*distance*distance); /*we take the projection of WD2 to be a flat circle, it should be round, but*/
Pin[n1]=albedo*Pout[n1]*irradiatedpart;
/*since distance>3*R2 this is a minor effect*/
P_bl[n1]= 0.5*Mdot*(Msol/year)*(getRochePotential(L1[0],L1[1])-getRochePotential((a1[0]),(a1[1]+R1))); /*mind the 0.5* !!!! */
// printf("R2=%4.2le, distance=%4.2le a=%4.2le irradiatedpart=%4.2f \n",R2, distance, a, irradiatedpart);
// printf("Pout=%4.2le Pin=%4.2le Tout=%4.2le lambda_max=%4.2le nu_max=%4.2le \n", Pout, Pin, Tout, lambda_max, nu_max);
}
/************************************************************************/
/* This function is used to see in what range variables are plotted. */
/* It finds minium- and maximum value of a float data set. */
/************************************************************************/
void getminmax(float *a, int na, float *amin, float *amax) /*data points, number of points, returned min, returned maximum) */
{
int i;
*amax = a[0];
*amin = *amax;
for(i=0;i<na;i++)
{
42
if (*amin>a[i]) *amin = a[i];
if (*amax<a[i]) *amax = a[i];
}
return;
}
/********************************************************************************/
/* This function plots the axis for the spatial plot. This can be done
*/
/* by manual imput of the range (in this function) then autoscale = 0,
*/
/* or it can automatically scale the plot to fit in the screen. (autoscale=1)
*/
/* It also puts to stars at the centre of WD1 and WD2. */
/********************************************************************************/
void Plotspatialaxis()
{
float factor=1.1;
float xmin1, xmax1; /*minimum xvalues for the rochelobe around M1*/
float xmin2, xmax2; /*minimum xvalues for the rochelobe around M2*/
float ymin1, ymax1; /*minimum yvalues for the rochelobe around M1*/
float ymin2, ymax2; /*minimum yvalues for the rochelobe around M2*/
float xmin, xmax, ymin, ymax; /*min and max values for the axis*/
cpgsls(1); /*select linestyle; solid */
if(autoscale==1)
{
getminmax(roche1x[0],rocheBins,&xmin1, &xmax1);
getminmax(roche1x[1],rocheBins,&ymin1, &ymax1);
getminmax(roche2x[0],rocheBins,&xmin2, &xmax2);
getminmax(roche2x[1],rocheBins,&ymin2, &ymax2);
if(xmin1<xmin2)
xmin=xmin1;
else
xmin=xmin2;
if(xmax1>xmax2)
xmax=xmax1;
else
xmax=xmax2;
if(ymin1<ymin2)
ymin=ymin1;
else
ymin=ymin2;
if(ymax1>ymax2)
ymax=ymax1;
else
ymax=ymax2;
cpgsls(1);
cpgenv(factor*xmin, factor*xmax, factor*ymin, factor*ymax, 1, 0); /* sets the axis 4th:1=same scale */
}
else
{ /* =(a1[0]-10.7e7, a1[0]+5e7, -5e7, 5e7, 1, 0)
*/
//
cpgenv(a1[0]-10.7e7, a1[0]+5e7, -5e7, 5e7, 1, 0); /* with this scale the direct impact phase fits just in the screen
cpgenv(a1[0]-3.2e8 ,a1[0]+2.2e8,-2.2e8 ,2.2e8 , 1, 0); /*(a1[0]-3.2e8 ,a1[0]+2.2e8,-2.2e8 ,2.2e8 , 1, 0) */
} /* with this scale the diskfase upto t=1e9 yr fits just in the screen */
cpgpt1(a1[0],a1[1],3); /*Plots CM of M1 */
cpgpt1(a2[0],a2[1],3); /*Plots CM of M2 */
cpglab("x (m)", "y (m)", "Balistic stream");
}
*/
/****************************************************************/
/* This function draws the Roche Lobe by drawing lines between */
/* the point of the Roche Lobe */
/****************************************************************/
void Plotrochelobes() /*Plots the axis and the Roche lobes */
{
Plotspatialaxis();
cpgsls(2); /*set line style dashed */
cpgline((rocheBins+1), roche1x[0], roche1x[1]); /*plots lines around M1 */
cpgsls(1); /*set line style solid */
cpgline((rocheBins+1), roche2x[0], roche2x[1]); /*plots lines around M2 */
// printf("ymin=%4.2le ymin1=%4.2le ymin2=%4.2le rocheBins=%d \nand roche1x[0][0]=%4.2le roche1x[0][rocheBins/2]=%4.2le \n",
// ymin, ymin1, ymin2,rocheBins, roche1x[0][0], roche1x[0][rocheBins/2]);
// printf("xmin=%4.2le xmax=%4.2le ymin=%4.2le ymax=%4.2le \n", xmin, xmax, ymin, ymax);
}
/***********************************/
/* This function plots the primary */
/***********************************/
void PlotM1()
{
cpgsfs(2); /*set fillstyle*/
cpgcirc(a1[0],a1[1],R1); /*draw circle
}
*/
/*********************************************/
/* This function plots the ballistic stream */
/*********************************************/
void Plotballistics()
{
cpgsls(1);
cpgline(steps2, xval, yval);
}
/*************************************************************************/
/* This function plots some construction lines (A1, A2 and B), that are */
/* mentioned in getIrradiatedsurface(). The use of this function is */
/* optional and insightful. */
/*************************************************************************/
void Plotconstructionlines() /*if this is turned on lines A1, A2, B, and a horizontal line are drawn */
{ /*it helps to understand getIrradiation. */
float directionnormal[2]; /*the x- and y-components of a vector between a1[] and impact[] */
directionnormal[0]=impact[0]/(sqrt(impact[0]*impact[0]+impact[1]*impact[1]));
43
directionnormal[1]=impact[1]/(sqrt(impact[0]*impact[0]+impact[1]*impact[1]));
cpgsls(4);
/*set linestyle*/
cpgmove(a1[0],a1[1]);
/*moves pen position to CM of M1; is used together with the next line...*/
cpgdraw((3*impact[0]+a1[0]), (3*impact[1]+a1[1]));
/*...draws the normal on the surface of WD1
*/
cpgmove(impact[0]+a1[0]+5*R1*directionnormal[1], impact[1]+a1[1]-5*R1*directionnormal[0]); /*draws line B through... */
cpgdraw(impact[0]+a1[0]-5*R1*directionnormal[1], impact[1]+a1[1]+5*R1*directionnormal[0]); /*...impact and perpendicular to the normal */
cpgmove(a2[0]-R2*directionnormal[0],a2[1]-R2*directionnormal[1]);
/*draws line through a2[] parallel to the normal... */
cpgdraw(a2[0]+R2*directionnormal[0],a2[1]+R2*directionnormal[1]);
/*...to normal, with a length of 2*R2.
*/
cpgmove((impact[0]+a1[0]-4*R1), (impact[1]+a1[1]));
/*draws horizontal line through... */
cpgdraw((impact[0]+a1[0]+4*R1), (impact[1]+a1[1]));
/*impact[] */
cpgsls(2);
cpgmove((impact[0]+a1[0]), (impact[1]+a1[1])); /*draws lineA1 between impact and... */
cpgdraw(a2[0]-R2*directionnormal[0],a2[1]-R2*directionnormal[1]); /*lower part of WD1
*/
cpgmove((impact[0]+a1[0]), (impact[1]+a1[1]));
/*draws lineA2 between impact and... */
cpgdraw(a2[0]+R2*directionnormal[0],a2[1]+R2*directionnormal[1]); /*upperpart of WD1
*/
}
/************************************************/
/* This function plots the irradiated surface */
/************************************************/
void Plotirradiatedsurface()
{
cpgsfs(4); /*set fill style; to color the irradiated surface.*/
if((irradiation==1)||(irradiation==2))
{
cpgpoly(p, irradiatedsurface[0],irradiatedsurface[1]);
/*draws a polygon */
}
}
/************************************************************************/
/* This function make the plot between X1[m]...X5[m] and Y1[m]....Y5[m] */
/* m changes for different values of M1 and M2. */
/************************************************************************/
void PlotYvsX() /*Makes a plot of the X1 vs Y1, X2 vs Y2 and X3 vs Y3, for different m, i, j */
{ /*X, Y are set in getXvsY() and the labels of the axis are set in main() */
float xmin1, xmax1, ymin1, ymax1; /*minimum and maximum of X1 and Y1*/
float xmin2, xmax2, ymin2, ymax2; /*minimum and maximum of X2 and Y2*/
float xmin3, xmax3, ymin3, ymax3; /*minimum and maximum of X3 and Y3*/
float xmin4, xmax4, ymin4, ymax4; /*minimum and maximum of X4 and Y4*/
float ymin5, ymax5; /*minimum and maximum of Y5*/
float xmin, xmax, ymin, ymax; /*min and max values for the axis*/
float deltax, deltay; /*range of X and of Y */
cpgslct(plotXvsY); /*selects the figure to plot*/
getminmax(X1, n1, &xmin1, &xmax1); /*gat the minimum and maximum values of 1,2 and 3*/
getminmax(Y1, n1, &ymin1, &ymax1);
getminmax(X2, n1, &xmin2, &xmax2);
getminmax(Y2, n1, &ymin2, &ymax2);
getminmax(X3, n1, &xmin3, &xmax3);
getminmax(Y3, n1, &ymin3, &ymax3);
getminmax(X4, n1, &xmin4, &xmax4);
getminmax(Y4, n1, &ymin4, &ymax4);
getminmax(Y5, n5, &ymin5, &ymax5);
xmin=xmin1; /*compare the minimum and maximum values of 1,2 and 3 and find the extremes of them*/
xmax=xmax1;
ymin=ymin1;
ymax=ymax1;
if(xmin2<xmin)
xmin=xmin2;
if(xmax2>xmax)
xmax=xmax2;
if(ymin2<ymin)
ymin=ymin2;
if(ymax2>ymax)
ymax=ymax2;
if(xmin3<xmin)
xmin=xmin3;
if(xmax3>xmax)
xmax=xmax3;
if(ymin3<ymin)
ymin=ymin3;
if(ymax3>ymax)
ymax=ymax3;
if(xmin4<xmin)
xmin=xmin4;
if(xmax4>xmax)
xmax=xmax4;
if(ymin4<ymin)
ymin=ymin4;
if(ymax4>ymax)
ymax=ymax4;
if(ymin5<ymin)
ymin=ymin5;
if(ymax5>ymax)
ymax=ymax5;
deltax=xmax-xmin;
deltay=ymax-ymin;
// printf("xmin1=%4.2f and xmax1=%4.2f and ymin1=%4.2f and ymax1=%4.2f \n", xmin1, xmax1, ymin1, ymax1);
// printf("xmin2=%4.2f and xmax2=%4.2f and ymin2=%4.2f and ymax2=%4.2f \n", xmin2, xmax2, ymin2, ymax2);
// printf("xmin3=%4.2f and xmax3=%4.2f and ymin3=%4.2f and ymax3=%4.2f \n", xmin3, xmax3, ymin3, ymax3);
// printf("xmin=%4.2f and xmax=%4.2f and ymin=%4.2f and ymax=%4.2f \n", xmin, xmax, ymin, ymax);
cpgenv(xmin, xmax, ymin+0.1*deltay, ymax, 0, 0); /*sets the axis 4th:0=different scale; 0.9=>1.1 if xmin<0 */
if(heating==0)
{
cpgsls(1); /*1=solid*/ //WD1 without heating up
cpgline(n1,X1,Y1);
cpgsls(2); /*2=dashed*/ //WD2 without heating up
cpgline(n1,X2,Y2);
}
cpgsls(4); /*4=dotted*/ //Pin
44
cpgline(n1,X3,Y3);
cpgsls(5); /*5=dash-dot-dot-dot*/ //Pout
cpgline(n1,X4,Y4);
cpgsls(3); /*3=dot-dash-dot-dash*/ //P_bl
cpgline(n5,X5,Y5);
}
/****************************************************************/
/* For a given a WD mass (M) and cooling time (t) this */
/* function calculates the luminosity of the star. */
/* It is also used to get right boundary conditions to solve */
/* the differential equation in PlotnewLwd() */
/****************************************************************/
float get Luminosity(float M, float t)
{
float Lmax; /*in Lsol*/
float logL;
if(M<0.5)
{
Lmax=9.4-1.33*(0.45-M);
logL=Lmax-1.33*log10(t);
//
if(logL>-0.5)
//
{logL=-0.5;}
}
if(M>=0.5)
{
Lmax=9-(0.9-M);
logL=Lmax-1.33*log10(t);
//
if(logL>2)
//
{logL=2;}
}
// printf("log10L=%4.2f ",logL);
return(pow(10,logL));
}
/************************************************************************/
/* This function gets the X1[m]...X5[m] coordinates that are plotted */
/* against Y1[m]...Y5[m]. */
/************************************************************************/
void getXvsY()
{
float Xvalue, Yvalue;
Xvalue=log10(t_ev); /*Select X-axis for the second plot. it is not time PlotnewLwd() should not be plotted*/
t3[n1]=t_ev;
X1[n1]=Xvalue;
Y1[n1]=log10(getLuminosity(M1initial,t_ev+tbegin)); /*cooling of WD1 without heating up*/
X2[n1]=Xvalue;
Y2[n1]=log10(getLuminosity(M2initial,t_ev+tbegin)); /*cooling of WD2 without heating up*/
X3[n1]=Xvalue;
if(log10(Pin[n1]/Lsol)>-7)
{Y3[n1]=log10(Pin[n1]/Lsol);}
else
{Y3[n1]=-7;}
X4[n1]=Xvalue;
if(log10(Pout[n1]/Lsol)>-7)
{Y4[n1]=log10(Pout[n1]/Lsol);}
else
{Y4[n1]=-7;}
printf("n1=%d X1=%4.2f Y1=%4.2f Y2=%4.2f Y3=%4.2f Y4=%4.2f \n", n1, X1[n1], Y1[n1], Y2[n1], Y3[n1], Y4[n1]);
if(noimpact==1)
{
X5[n5]=Xvalue;
Y5[n5]=log10(P_bl[n1]/Lsol);
//
printf("n5=%d X5=%4.2f Y5=%4.2f Y2=%4.2f\n", n5, X5[n5],Y5[n5], Y2[n1]);
n5++;
}
n1++;
}
/************************************************************************/
/* This function chops the Pin-graf up and it roughly ‘integrates‘ */
/* between the calculated point, to get a better resolution for 4000 */
/* in stead of ~40 steps.
*/
/************************************************************************/
float getPout(float t2)
{
int n=0;
float P; /*in J/s*/
float stop=0;
while(stop==0)
{
if(t2>t_ev)
{
printf("error in Pout()");
break;
}
if(t2>t3[n])
{n++;}
else
{stop=1;}
}
// printf("n=%d",n);
if(t2<tdisk) /*in the direct impact phase Pout is the power released in the impact spot.*/
{return(Pout[n]/Lsol);}
else /*in the disk phase Pout is the power released in the boundary layer. */
{return(P_bl[n]/Lsol);}
}
45
/************************************************************************/
/* This function chops the Pin-graph up and it roughly ‘integrates‘ */
/* between the calculated point, to get a better resolution for 4000 */
/* in stead of ~40 steps.
*/
/************************************************************************/
float getPin(float t2)
{
int n=0;
float P; /*in J/s*/
float stop=0;
while(stop==0)
{
if(t2>t_ev)
{
printf("error in Pin()");
break;
}
if(t2>t3[n])
{n++;}
else
{stop=1;}
}
return(Pin[n]/Lsol);
}
/****************************************************************/
/* This function calculates Lmax for a given M */
/****************************************************************/
float getLmax(float M)
{
float logLmax;
if(M<0.5) /*first get Lmax for different values of M*/
{logLmax=9.4-1.33*(0.45-M);}
if(M>=0.5)
{logLmax=9-(0.9-M);}
return(pow(10, logLmax));
}
/****************************************************************/
/* This function calculates and plots the new WD luminosities */
/* taking into account the heating of the WDs. */
/* getPin() and getPout() are used, because this function needs */
/* more steps to solve the integral than the original program */
/* uses. */
/****************************************************************/
void PlotnewLwd(float M, int linestyle)
{
float Lmax; /*in Lsol*/
float logLmax;
float L; /*in Lsol*/
float t2; /*in years*/
float t2step;
float Ldot; /*an infinidecimal step in L*/
float x[4001],y[4001], z[4001];
int n;
t2=1e4;
Lmax=getLmax(M);
t2step=t2/100;
L=getLuminosity(M,tbegin);
logLmax=log10(Lmax);
// printf("logLmax=%4.2le, Lmax=%4.2le, L=%4.2le \n", logLmax, Lmax, L);
n=0;
while(t2<t_ev/factor)
{
if(linestyle==1)
{
Ldot=1.333*pow(Lmax,-0.75)*(0.5*getPout(t2)-L)*pow(L,0.75);
}
if(linestyle==2)
{
Ldot=1.333*pow(Lmax,-0.75)*(getPin(t2)-L)*pow(L,0.75);
}
L+=Ldot*t2step;
t2+=t2step;
x[n]=log10(t2);
y[n]=log10(L); /*calculated luminosity with heating*/
z[n]=log10(getLuminosity(M,t2)); /*the old luminosity without heating*/
t2step*=1.005;
//
printf("n=%d x=%4.2g y=%4.2g, z[n]=%4.2g Ldot=%4.2le, -log10Ldot=%4.2le
// n, x[n], y[n], z[n], Ldot, log10(-Ldot), getPin(t2));
n++;
}
// printf("logLmax=%4.2le, Lmax=%4.2le, L=%4.2le \n", logLmax, Lmax, L);
cpgsls(linestyle);
cpgline(n,x,y);
// cpgline(n,x,z);
}
Pin=%4.2g\n",
int main()
{
/*checking if PGplot is working*/
if ((plotXvsY = cpgopen(deviceXvsY)) <= 0)
/* Open Y vs X plot*/
{
return EXIT_FAILURE;
}
if ((plotSpatial = cpgopen(deviceSpatial)) <= 0) /* Open Spatial Coordinates plot*/
{
46
return EXIT_FAILURE;
}
printf("test \n");
/*initializing some parameters*/
/*these variables can be changed freely:*/
tbegin=1e8; /*cooling time*/
M1=0.5;
M2=0.25;
factor=1.5; /*1.3=fast. 1.01 accurate*/
heating=1; /*do the WDs heat up? 1=> yes; 0=>no.*/
factorw1=0.5;
/*these variables should stay the same:*/
n1=0; /* n1..n4 change for different values of M1 and M2 */
M1initial=M1;
M2initial=M2;
timestep=7e3;
t_ev=timestep;
/*calculation the evolution of of the direct impact phase.*/
while(1)
{
noimpact=1;
M2+=M2dot*timestep; /* mass is conserved; Note: M2dot<0
*/
M1-=M2dot*timestep;
getDynamics();
getAngularmomentum();
getRocheLobes();
getBallistics(R1); /* noimpact can be set to 0 here if there is impact
*/
if(noimpact==1)
{
tdisk=t_ev;
printf("no more impact and tdisk=%4.2le \n", tdisk);
break; /* if there is no more impact, the diskphase ends.
*/
}
getIrradiatedsurface(); /*it calculates both the openingangle and the irradiatedsurface*/
getImpactdynamics();
getXvsY();
t_ev+=timestep;
timestep=timestep*factor;
//{} /*Close the loop here to make only one plotSpatial
*/
cpgslct(plotSpatial); /*Selects the figure to plot*/
cpgask(0); /*0=> new spatial graphs are automatically drawn. 1=>press enter for a new graph */
Plotrochelobes();
PlotM1();
Plotballistics();
Plotirradiatedsurface();
//
Plotconstructionlines(); /*This function is optional
*/
} /*Close the loop here to make a plotSpatial for all m */
/*calculation of the evolution of the disk impact phase */
while(t_ev<1e10)
{
getDynamics();
M2+=M2dot*timestep;
M1-=M2dot*timestep;
getAngularmomentum();
getRocheLobes();
getBallistics(R_disk); /*Note: there should always be an impact on the disk.
*/
getIrradiatedsurface(); /*it calculates both the openingangle and the irradiatedsurface
*/
getImpactdynamics();
noimpact=1; /*just to values for P_bl are taken... */
getXvsY(); /*...in getimpact */
t_ev+=timestep;
timestep=timestep*factor;
//{} /*Close the loop here to make only one plotSpatial
*/
cpgslct(plotSpatial); /*Selects the figure to plot*/
cpgask(0); /*0=> new spatial graphs are automatically drawn. 1=>press enter for a new graph */
Plotrochelobes();
PlotM1();
Plotballistics();
Plotirradiatedsurface();
cpgsfs(2); /*set fill style*/
cpgcirc(a1[0],a1[1],R_disk); /*Plotdisk, plots the disk*/
} /*Close the loop here to make a plotSpatial for all m (or i or j) */
cpgclos(); /*Closes plotSpatial*/
PlotYvsX();
if(heating==1)
{
PlotnewLwd(M1initial,1); /*gets and plots a cooling curve...*/ /* (mass of WD, accretor=1; donor=2)*/
PlotnewLwd(M2initial,2); /*... with accretional heating
*/
}
cpglab("log\\u10\\d(t) (years) ", "log\\u10\\d(L) (Solar luminosity)","");
cpgclos(); /*Closes plotXvsY*/
printf("heating=%d, M1=%4.2f M2=%4.2f log10(tbegin)=%4.2f, factorw1=%4.2f \n", heating, M1initial, M2initial,
log10(tbegin), factorw1);
return EXIT_SUCCESS;
}
47