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: [email protected] 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

© Copyright 2019