UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Accreting white dwarfs : a theoretical analysis of nuclear burning Pillay, Samara 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2013_fall_pillay_samara.pdf [ 1.25MB ]
Metadata
JSON: 24-1.0074134.json
JSON-LD: 24-1.0074134-ld.json
RDF/XML (Pretty): 24-1.0074134-rdf.xml
RDF/JSON: 24-1.0074134-rdf.json
Turtle: 24-1.0074134-turtle.txt
N-Triples: 24-1.0074134-rdf-ntriples.txt
Original Record: 24-1.0074134-source.json
Full Text
24-1.0074134-fulltext.txt
Citation
24-1.0074134.ris

Full Text

Accreting White Dwarfs: ATheoretical Analysis of NuclearBurningbySamara PillayB.Sc., The University of the Witwatersrand, 2004B.Sc. (Hons.), The University of the Witwatersrand, 2005M.Sc., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2013c? Samara Pillay 2013AbstractAccreting white dwarfs can exhibit a variety of thermonuclear phenomena,such as shell flashes, classical and recurrent novae, as well as Type Ia su-pernovae. To better understand these processes, we consider the accretionof hydrogen-rich material onto the surface of a white dwarf. Our analy-sis is based on a semi-analytical approach that allows the investigation ofproperties of nuclear burning on accreting white dwarfs. In particular, wedetermine steady-state solutions and evaluate the stability of these solu-tions. As a first step, we follow Paczyn?ski?s one-zone model and confirm hisresults by following his analysis independently. We extend the framework toa sophisticated multi-zone model encompassing a variety of detailed physics.We determine accretion rates that may lead to stable or to unstable burn-ing. Regimes of stable burning may result in mass increase and potentiallyidentify progenitors of Type Ia supernovae. Unstable burning may lead tonova-like outbursts. The identification of both burning regimes is impor-tant, as these thermonuclear events influence the chemical and dynamicalevolution of the Universe.iiPrefaceThe work presented in this thesis was independently conducted by theauthor, S. Pillay, under the guidance of supervisors, Professors JaymieMatthews and Jeremy Heyl. The work contained within Chapter 3 was inde-pendently conducted, but followed the work outlined in Bohdan Paczyn?ski,?A One-Zone Model for Shell Flashes on Accreting Compact Stars?, TheAstrophysical Journal, 264:282-295, 1983 January 1. Figure (3.2b), whichappeared in Paczyn?ski?s publication as Figure 1, has been reproduced withpermission from the American Astronomical Society (AAS).The Runge-Kutta-Fehlberg method used in this thesis, based on secondand third order Runge-Kutta routines, was adapted from an algorithm pro-vided by Professor Anthony Peirce. A code for neutron stars was providedby Professor Jeremy Heyl, my co-supervisor, as a tool with which to com-pare my code for white dwarfs. I modified the inherited code for neutronstars to suit the purposes of the research contained within this thesis. Theresults generated from this modified code appear in the appendix.The remainder of the research contained within this thesis, barring ref-erence material, is unpublished and original.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction and Review . . . . . . . . . . . . . . . . . . . . . 11.1 White Dwarf Characteristics . . . . . . . . . . . . . . . . . . 21.2 Cataclysmic Variables . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Dwarf Novae . . . . . . . . . . . . . . . . . . . . . . . 31.2.2 Classical and Recurrent Novae . . . . . . . . . . . . . 41.2.3 Type Ia Supernovae . . . . . . . . . . . . . . . . . . . 41.2.4 Temperature and Mass in Binary Systems . . . . . . 51.3 Shell Flashes . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3.1 Relevant Literature . . . . . . . . . . . . . . . . . . . 61.4 Motivation and Goals . . . . . . . . . . . . . . . . . . . . . . 101.5 Outline and Description of Model . . . . . . . . . . . . . . . 112 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1 Model Framework . . . . . . . . . . . . . . . . . . . . . . . . 132.2 Mean Molecular Weight . . . . . . . . . . . . . . . . . . . . . 132.3 Energy Generation Rates . . . . . . . . . . . . . . . . . . . . 152.3.1 Nuclear Energy Sources . . . . . . . . . . . . . . . . . 152.3.2 Electron Capture . . . . . . . . . . . . . . . . . . . . 182.3.3 Gravitational Energy Sources . . . . . . . . . . . . . 192.4 Opacity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20ivTable of Contents2.4.1 Radiative Opacities . . . . . . . . . . . . . . . . . . . 202.4.2 Conductive Opacity . . . . . . . . . . . . . . . . . . . 222.4.3 Total Opacity . . . . . . . . . . . . . . . . . . . . . . 222.5 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.6 Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.7 Eddington Limit . . . . . . . . . . . . . . . . . . . . . . . . . 252.8 Lagrangian Derivative . . . . . . . . . . . . . . . . . . . . . . 262.9 Core Temperature, Mass and Surface Gravity . . . . . . . . . 262.10 Thermal Instability . . . . . . . . . . . . . . . . . . . . . . . 263 A One-Zone Model . . . . . . . . . . . . . . . . . . . . . . . . . 303.1 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . 303.1.1 Dimensionless Quantities . . . . . . . . . . . . . . . . 333.2 Steady-State Solutions . . . . . . . . . . . . . . . . . . . . . 343.2.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . 353.2.2 Numerical Scheme . . . . . . . . . . . . . . . . . . . . 353.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . 373.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454 A Multi-Zone Approach . . . . . . . . . . . . . . . . . . . . . . 554.1 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . 554.2 Steady-State Equations . . . . . . . . . . . . . . . . . . . . . 564.2.1 Surface Boundary Conditions . . . . . . . . . . . . . . 574.2.2 Inner Boundary Condition . . . . . . . . . . . . . . . 594.3 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . 614.3.1 Imposing the Inner Boundary Condition . . . . . . . 624.4 Steady-State Solutions . . . . . . . . . . . . . . . . . . . . . 654.4.1 XCNO = Zs/3 . . . . . . . . . . . . . . . . . . . . . . 664.4.2 XCNO = 0.8Zs . . . . . . . . . . . . . . . . . . . . . . 714.5 Other Accretion Rates with XCNO = 0.8Zs . . . . . . . . . . 754.5.1 2? 10?9 M yr?1 ? M? ? 6.3? 10?7 M yr?1 . . . . . 764.5.2 M? < 2? 10?9 M yr?1 . . . . . . . . . . . . . . . . . 794.5.3 Accumulated Masses . . . . . . . . . . . . . . . . . . 804.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81vTable of Contents5 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 1115.1 Equations for Perturbed Variables . . . . . . . . . . . . . . . 1125.2 Equations for Real and Imaginary Parts of Perturbed Vari-ables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1145.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 1185.4 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . 1195.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1226 Conclusions, Future Work and Other Remarks . . . . . . . 1236.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 1276.2 Potential Observational Confirmations of Theory . . . . . . . 128Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130AppendixA Neutron Star . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134viList of Figures2.1 Ratio of CN to pp as a function of temperature, with X = 0.7and Z = 0.03. The red line indicates the point at which theenergy generation rates are equal. . . . . . . . . . . . . . . . . 293.1 Column density as a function of log a?, where a? is the di-mensionless accretion rate. The steady-state solutions areshown for a range of dimensionless heat flux from the core,fb. The hydrogen and helium mass fractions at the surfaceare X = 0.7 and Y = 0.27, respectively. The surface gravityis g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . . . . . 483.2 Comparison of results . . . . . . . . . . . . . . . . . . . . . . 493.3 Steady-state solutions for fb = 0.03, 0.01, 0.003 with X =0.7, Y = 0.27 and g = 108.5 cm s?2. . . . . . . . . . . . . . . . 503.4 Variation of eigenvalues, ?1 and ?2, computed using the time-dependent system of equations, equations (3.21) and (3.22),with dimensionless accretion rate for fb = 0.01, with X =0.7, Y = 0.27 and g = 108.5 cm s?2. In the vicinity of theturning points, the eigenvalues become complex conjugatesof each other. Thus, only the red line indicating ?2 is visiblenear the turning points, as ?1 has the same real part overthis portion of the domain. The real part of both eigenvaluescrosses from a negative to a positive quantity at the localmaximum and vice versa at the local minimum, indicating atransition to instability and stability, respectively. . . . . . . . 513.5 Variation of the thermal and nuclear eigenvalues with dimen-sionless accretion rate for fb = 0.01, with X = 0.7, Y = 0.27and g = 108.5 cm s?2. These eigenvalues are valid in the limit?n/?th  1, where |?n/?th|  1 [28], corresponding to regionsaway from the turning points. . . . . . . . . . . . . . . . . . . 523.6 Gradient of the linear series as a function of log a? for fb =0.01, with X = 0.7, Y = 0.27 and g = 108.5 cm s?2. . . . . . . 52viiList of Figures3.7 Variation of the eigenvalues, including the thermal and nu-clear eigenvalues, with log a? for fb = 0.01, with X = 0.7, Y =0.27 and g = 108.5 cm s?2. This is a close-up of the local min-imum for the curve fb = 0.01. The real parts of ?1 and ?2pass through zero at the local minimum, indicating a transi-tion to stability. The thermal eigenvalue passes through zeroat a slightly different point to the local minimum while thenuclear eigenvalue diverges. This is a consequence of the factthat the thermal and nuclear eigenvalues are not valid in thevicinity of the turning points, and ?1 and ?2 are required togive an accurate description of the behavior at the turningpoints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.8 Variation of temperature with dimensionless accretion ratefor fb = 0.01, with X = 0.7, Y = 0.27 and g = 108.5 cm s?2. . 543.9 Variation of density with dimensionless accretion rate for fb =0.01, with X = 0.7, Y = 0.27 and g = 108.5 cm s?2. . . . . . . 544.1 Variation of column density, ??, with dimensionless fluxfrom the surface, fs, for a dimensionless accretion rate ofa? = 0.0131. The surface hydrogen, helium and CNO massfractions are Xs = 0.7, Ys = 0.27 and XCNO = Zs/3, respec-tively. The surface gravity is g = 108.5 cm s?2. . . . . . . . . . 864.2 Solutions for fs = 0.0016 = 0.12a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The redline indicates the core temperature. There are at least twosteady-state solutions here. . . . . . . . . . . . . . . . . . . . 874.3 Solutions for fs = 0.0064 = 0.49a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The redline indicates the core temperature. There are two steady-state solutions, separated by a larger column density than infigure (4.2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.4 Solutions for fs = 0.013 = a?, a? = 0.0131, with Xs = 0.7, Ys =0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line in-dicates the core temperature. There is only one steady-statesolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 894.5 Solutions for fs = 0.026 = 2a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. Thered line indicates the core temperature. There is only onesteady-state solution. . . . . . . . . . . . . . . . . . . . . . . . 90viiiList of Figures4.6 Variation of column density with fs for a? = 0.013, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The redline divides the figure horizontally into a top region and abottom region. The red text above the red line describes thesolutions in the top region while the red text below the red linedescribes characteristics of the solutions in the bottom region.The dashed black lines split the figure vertically indicatingthree regions corresponding to no solution, two solutions andone solution from left to right. . . . . . . . . . . . . . . . . . 914.7 Comparison of temperature variation with fs at different col-umn densities, with a? = 0.0131, Xs = 0.7, Ys = 0.27, XCNO =Zs/3 and g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . 924.8 Solutions for fs = 0.02 = 1.5a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. Thered line indicates the core temperature. The plateaus in theflux are distinctive of an exhaustion of nuclear fuel. . . . . . . 934.9 Variation of column density with fs for a? = 0.0131, withXs = 0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2, inthe extended model, which includes the proton-proton chains,electron capture, conduction and degeneracy pressure in boththe relativistic and non-relativistic regimes. . . . . . . . . . . 944.10 Variation of column density with fs for a? = 0.0131, withXs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, inthe extended model. . . . . . . . . . . . . . . . . . . . . . . . 954.11 Variation of column density with fs for a? = 0.0131, withXs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, inthe extended model. . . . . . . . . . . . . . . . . . . . . . . . 954.12 Variation of column density with fs for a? = 0.0131, withXs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, inthe extended model. . . . . . . . . . . . . . . . . . . . . . . . 964.13 Solutions for fs = 0.01280 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. . . . . . . . . . . . . . . . 974.14 Solutions for fs = 0.01285 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. There are at least twosteady-state solutions. . . . . . . . . . . . . . . . . . . . . . . 984.15 Variation of hydrogen mass fraction with column density fora? = 0.0131, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2, in the extended model. . . . . . . . . . . . . 99ixList of Figures4.16 Solutions for fs = 0.0130 ? 0.99a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. There are two steady-state solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.17 Solutions for fs = 0.0136 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. There are two steady-state solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . 1014.18 Solutions for fs = 0.0139 ? 1.06a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. There are two steady-state solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.19 Solutions for fs = 0.0141 = 1.08a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The redline indicates the core temperature. There are at least twosteady-state solutions. . . . . . . . . . . . . . . . . . . . . . . 1034.20 Variation of column density with fs for a? = 0.013, withXs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2.The red line divides the figure horizontally into three regions,a bottom, middle and top region. The red text describes fea-tures of the solutions in the respective regions. The dashedblack lines split the figure vertically indicating four regionscorresponding to no solution, two solutions, two solutions andone solution from left to right, indicated by black numerals. . 1044.21 Comparison of temperature with fs at different column den-sities, with a? = 0.0131, Xs = 0.7, Ys = 0.27, XCNO = Zs/3and g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . . . . 1054.22 Steady-state solutions for a? = 1.31, corresponding to M? =6.3? 10?7 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zsand g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . . . . 1064.23 Steady-state solutions for a? = 0.13, corresponding to M? =6.3? 10?8 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zsand g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . . . . 1064.24 Steady-state solutions for fs, corresponding to a? = 0.0041,M? = 2 ? 10?9 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO =0.8Zs and g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . 1074.25 Steady-state solutions for a? = 0.0013, corresponding to M? =6.3?10?10 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zsand g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . . . . . 107xList of Figures4.26 Steady-state solutions for a? = 1.31 ? 10?4, corresponding toM? = 6.3? 10?11 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO =0.8Zs and g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . 1084.27 Steady-state solutions for a? = 1.31 ? 10?5, corresponding toM? = 6.3? 10?12 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO =0.8Zs and g = 108.5 cm s?2. . . . . . . . . . . . . . . . . . . . 1084.28 Variation of hydrogen mass fraction with column density fora? = 1.31, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2, in the extended model. . . . . . . . . . . . . 1094.29 Variation of hydrogen mass fraction with column density fora? = 0.0041, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2, in the extended model. . . . . . . . . . . . . 1094.30 Variation of hydrogen mass fraction with column density fora? = 0.0013, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2, in the extended model. . . . . . . . . . . . . 1104.31 Accumulated mass in accreted layer as a function of columndensity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110A.1 Comparison of steady-state results using two different codesfor the neutron star with Xs = 0.7, Ys = 0.28, XCNO = 0.8Zsand g = 1014.3 cm s?2. . . . . . . . . . . . . . . . . . . . . . . 138A.2 Variation of column density with log fout, with Xs = 0.7, Ys =0.28, XCNO = 0.8Zs and g = 1014.3 cm s?2 . . . . . . . . . . . 139xiAcknowledgementsI would like to thank my supervisors, Professors Jaymie Matthews andJeremy Heyl, for their guidance.xiiDedicationTo my sister, a person of great strength and courage.xiiiChapter 1Introduction and ReviewWhite dwarfs are extremely dense stellar remnants that are held up againstgravitational collapse by electron degeneracy pressure. White dwarfs areformed through the final evolutionary stages of low mass stars, with a massless than eight solar masses [9]. As a result of this mass limit, and the factthat high mass stars are rare, the majority of stars in the Galaxy will termi-nate their lives as white dwarfs [9]. These stars spend most of their lifetimesburning hydrogen on the main sequence [16]. Once the star?s hydrogen fuelis exhausted, and after several evolutionary steps, helium burning is ignitedwithin the core [16]. The products of helium burning, carbon and oxygen,generally do not undergo burning [16]. Consequently, these stars collapseunder gravity until the collapse is halted by electron degeneracy pressure[9], at a radius comparable to that of the Earth?s.As a result of this evolutionary process, the majority of white dwarfs?cores are composed of carbon and oxygen, also referred to as a CO core[9], [16]. Furthermore, the CO core is enveloped by a thin layer of non-degenerate gas, a remnant of the evolutionary stage in which the outerlayers of the star are ejected [9], forming a planetary nebula. Hereafter,evolution of an isolated white dwarf is characterised by cooling, as the starslowly releases the energy attained through nuclear burning, via the non-degenerate envelope [16]. However, if the white dwarf is located in a binarysystem, in which it can accrete or accumulate mass from a companion star,the evolutionary path becomes more interesting [36]. In particular, a va-riety of thermonuclear phenomena may occur, which range from relativelylow energy events to explosions accompanied by mass loss [16]. Howeverviolent or non-violent the eruption, it is clear that these events influencethe chemical and dynamical evolution of the Universe. In this thesis, wefocus on hydrogen-accreting white dwarfs and the conditions that may leadto thermonuclear phenomena.11.1. White Dwarf Characteristics1.1 White Dwarf CharacteristicsWhite dwarfs are extremely dense and hot objects. The average surfacegravity is 108 cm s?2 and the average density is approximately 106 g cm?3[9]. The average mass of an isolated white dwarf is approximately 0.6 M[33], where M is one solar mass or approximately 1.989? 1033 g [16]. Withthese characteristics, white dwarf radii are approximately 0.01 R, whereR is one solar radius or approximately 6.96? 1010 cm [16].As already mentioned, electron degeneracy pressure supports the coreagainst gravitational collapse. Degenerate electrons are efficient conductorsof heat, rendering the core essentially isothermal [9]. Heat is conductedout of the core relatively quickly, but must then diffuse through the non-degenerate envelope. Approximately 99% of the mass of a white dwarf iscontained within the core, while the non-degenerate envelope contains theremainder [9]. These thin outer layers, composed of hydrogen and helium,are extremely opaque to the radiation exiting the core, thereby acting as aninsulator [16]. As a result, isolated white dwarfs typically cool over extremelylong timescales, on the order of billions of years [33]. Furthermore, there isa large temperature difference between the core temperature and that of theeffective surface temperature, as a consequence of the opaque layers whichregulate the heat flow from the core [9]. White dwarf effective temperaturesrange from approximately 4,000 K to 150,000 K indicating a large range insurface luminosities as a result of the Stefan-Boltzmann law, L = 4piR2?T 4effwhere R is the radius of the white dwarf, Teff is the effective temperature,and ? is the Stefan-Boltzmann constant [9].The existence of degenerate electrons in the core is also responsible for avariety of unique characteristics. The degeneracy pressure is entirely depen-dent on the chemical composition and density, and is decoupled from thetemperature [16]. As a result of this equation of state, the electrons withinthe core cannot be packed any closer together and the white dwarf evolvesat essentially constant radius [9]. As a further consequence of degeneracypressure, white dwarfs satisfy an interesting mass-radius relationship, whichindicates that as they increase in mass, they decrease in radius [9]. Finally,the maximum mass of a stable, degenerate white dwarf is approximately1.46 M [21]. This limit is known as the Chandrasekhar limit [21].21.2. Cataclysmic Variables1.2 Cataclysmic VariablesWhite dwarfs that are found in a close binary system with a companion starexhibit a variety of interesting evolutionary properties and events. The com-panion star may be a red giant, a main sequence star or in some rare cases,another white dwarf [16]. As the binary evolves, the donor or companionstar may overflow its Roche lobe and transfer mass to the white dwarf [37],[36]. Consequently, mass from the donor star will accrete onto the surfaceof the white dwarf, adding to the non-degenerate layer. As a result of theaccretion, cataclysmic variables erupt periodically, with a sharp increase inbrightness over a day followed by a decline from peak brightness rangingfrom weeks to years [16]. The periodicity of eruptions varies from days tothousands of years, and depends largely on the physics behind the eruptions[36]. There are several types of eruptions that can occur, namely, dwarf no-vae, classical novae, recurrent novae and shell flashes. In the most extremecase, a Type Ia supernova resulting in the explosion of the white dwarf mayoccur.Typically, this accretion process occurs through an accretion disk, whichefficiently transfers angular momentum outwards in the disk while depositingmatter onto the equatorial plane of a non-magnetic white dwarf [37]. In thisthesis, we will not explicitly deal with the dynamics of the accretion disk.We will briefly discuss dwarf novae as well as classical and recurrent novae,which typically occur in non-magnetic cataclysmic variables [37]. We willalso give a brief overview of Type Ia supernovae and progenitors.1.2.1 Dwarf NovaeDwarf novae are caused by thermal instabilities in the accretion disk and arenot accompanied by mass loss [16]. Typically, white dwarfs exhibiting dwarfnovae phenomena accrete at low time-averaged rates of M? < 10?9 M yr?1[39]. However, once every month to a year, rapid mass transfer occurs fora week [39], resulting in an eruption. The white dwarf brightens sharplyover less than a day, increasing in brightness by a factor of 100, and thenfades over the next few weeks [37]. At the peak of the outburst, the opticalluminosity is approximately 1034 ergs s?1 [16]. Over the duration of theburst, the total energy released is approximately 1038?1039 ergs [16]. Dwarfnovae recur on a timescale of approximately weeks to months [36].31.2. Cataclysmic Variables1.2.2 Classical and Recurrent NovaeClassical novae recur over much longer periods than dwarf novae and are gen-erally accepted to occur as a result of unstable nuclear burning of hydrogen-rich material on the surface of a white dwarf [39], [19]. The orbital periodsof close binaries exhibiting classical novae range from 1.4 to 16 hours [39].These white dwarfs accrete matter from low mass companions at a time-averaged rate of M? ? 10?11?10?9 M yr?1 [39].Classical novae recur on a timescale of thousands to hundreds of thou-sands of years [36], [37]. Over the duration of an eruption, novae mayrelease between 1044?1045 ergs, with approximately 1038 ergs in the opticalband [16]. Classical novae are hydrodynamic events, accompanied by massejection of between 10?5?10?4 M [16]. During a nova, the white dwarfbrightens up to 10 million times over a few days and then fades over a pe-riod of months to years [37]. Classical novae can be further divided into fastor slow novae depending on the duration of decline from peak brightnessand expansion velocity [16]. Classical novae are defined as objects that haveerupted only once in recorded history, but may erupt again in the future[16].Recurrent novae have erupted more than once in recorded history, andtherefore have smaller periods. As a result, recurrent novae brighten muchless in magnitude than classical novae. RS Ophiuchi is a recurrent nova thathas erupted several times in recent history, with the last eruption takingplace in 2006 [17].1.2.3 Type Ia SupernovaeType Ia supernovae are the brightest of all supernovae [32]. Although it iswidely accepted that they occur as a result of the nuclear detonation of aCO white dwarf, uncertainty remains in the progenitor systems that lead tothese events [17]. The single degenerate scenario proposes a binary system inwhich a white dwarf accretes non-degenerate matter from a companion [40].In this case, a white dwarf must accrete enough mass from its companion toreach the Chandrasekhar limit of approximately 1.4 M [40]. At this limit,critical density is reached within the CO core, resulting in runaway andexplosion through carbon ignition [17]. It is postulated that stable burning,or scenarios that lead to mass increase of a CO white dwarf, are progenitorsof Type Ia supernovae [35]. However, the nature of the mass donor remainsuncertain and it is a challenge to simulate the growth of a white dwarf upto the Chandrasekhar limit.41.3. Shell Flashes1.2.4 Temperature and Mass in Binary SystemsWe briefly discussed the effective temperatures and masses of isolated whitedwarfs. However, these general properties are slightly different for whitedwarfs in mass-transferring binaries. As previously mentioned, the averagemass of the vast majority of isolated white dwarfs is 0.6 M. However,white dwarfs range in mass from 0.3 to 1.2 M, possibly as a result of binaryevolution [9]. It is reasonable to assume that white dwarfs in a binary systemare heavier [13], or at least exhibit a far wider range in mass than isolatedwhite dwarfs, which are sharply peaked around 0.6 M [36].Core temperatures are estimated to range from approximately 5 ? 106K to 2? 107 K [9]. Townsley and Bildsten [39] investigated the effect of anaccumulated layer of mass 5?10?5 M on the equilibrium core temperatureof a 0.6 M white dwarf for a range of accretion rates. The authors foundthat the core temperatures range from approximately 106 K to 4?107 K foraccretion rates ranging from approximately 10?11 to 10?9.4 M yr?1. Salariset al. [33] investigated cooling timescales for a variety of white dwarf masses.Their cooling sequences began at core temperatures of 5? 107 K [33].There is observational evidence that the core temperature of an accret-ing white dwarf is higher than that of an isolated white dwarf as a result ofaccretion [39], [36]. In a system exhibiting dwarf novae, the internal whitedwarf luminosity dominates the ultraviolet band between eruptions, as aresult of the drop in accretion rate, thus permitting measurement of the ef-fective temperature, Teff > 10,000 K [36], [39]. This elevated temperature ismost likely the result of compressional heating as the accreted layer becomesthicker. Indeed, as the layer becomes thicker, compressional heating heatsthe core [39].1.3 Shell FlashesShell flashes are essentially periodic increases in luminosity as a result ofthe nuclear burning of hydrogen or helium, known as hydrogen and heliumshell flashes, respectively. They may lead to classical and recurrent novae,or may lead to an increase in brightness not accompanied by mass loss,a non-hydrodynamical shell flash [4]. In general, shell flashes not leadingto mass loss are regimes of stable burning while novae are those for whichunstable burning occurs. Indeed, unstable burning models have producedhydrodynamic thermonuclear runaways resulting in a bolometric light curvethat is comparable, to some extent, to that of actual classical novae [19]. It isaccepted that classical novae eject most of the mass accumulated prior to the51.3. Shell Flasheseruption. Thus, it is the stable and steady-burning regimes that presumablyallow mass increase and may be progenitors of Type Ia supernovae. Due tothe importance of both novae and supernovae, many authors have attemptedto elucidate the parameter regimes and conditions that lead to mass loss ormass accumulation.While a detailed review of the literature pertaining to this problem isoutside the scope of this thesis, we briefly review some of the pertinent workin light of the analysis that will be presented in this thesis. We will compareour results to other authors where possible.1.3.1 Relevant LiteratureIn 1983, Paczyn?ski studied the problem of accretion of hydrogen and helium-rich matter onto the surface of compact stars [28]. In particular, he devel-oped a semi-analytical one-zone model to study the steady-state propertiesof nuclear burning and the stability of these solutions. In the one-zonemodel, a layer of hydrogen-rich material accumulates on the surface of thewhite dwarf. The layer is assumed to be thin and contains very little mass.As such, Paczyn?ski fixed the surface gravity in the layer and employed aplane-parallel geometry.The one-zone model describes the time evolution of the layer with twoordinary differential equations, and the steady-state behaviour with two non-linear equations. After determining the steady-state solutions, Paczyn?skiconducted a linear stability analysis to determine if the system is stable orunstable against perturbations of the underlying equilibrium solutions. Wewill review this work in great detail in chapter 3, however, we state theimportant results here in the case of white dwarfs.Paczyn?ski found that for very high or low accretion rates, the columndensity increased with increasing accretion rate. In these regimes, thesteady-state solutions were found to be stable. In regions where the columndensity decreased with increasing accretion rate, the steady-state solutionswere found to be unstable [28]. The analysis is restricted to accretion ratesM? < 10?7 M yr?1, as the model inherently prohibits surface fluxes largerthan or equal to the Eddington flux. The Eddington flux represents an upperlimit, above which gravitational forces are overcome by outwardly directedpressure and the surface layers of the star may no longer be bound [16].Following Paczyn?ski?s one-zone model, Jose? et al. implemented a two-zone model to study the properties of nuclear burning on accreting whitedwarfs [20]. The top layer consists of material that is accumulated throughaccretion, while the bottom layer is composed of the products of nuclear61.3. Shell Flashesburning. In the case of hydrogen accretion, the top layer is composed ofhydrogen, while the second layer is composed of helium. The authors used aplane-parallel geometry and specified the surface gravity and the flux fromthe core in accordance with Paczyn?ski. The layer is assumed to be thin andthus a plane-parallel geometry is used. The authors used a semi-analyticalapproach coupled with a computational algorithm, a fourth-order Runge-Kutta method, to investigate the evolution of the layer through a series ofshell flashes not leading to mass loss [20]. Neutrino losses and convectionwere neglected in this model. Both effects are generally negligible as a resultof the temperatures and densities reached within the accreted layer [20].The authors considered white dwarf masses in the range 1.0?1.2 M asthey were interested in massive white dwarfs as progenitors for Type Iasupernovae. Consequently, accretion rates in the range 2 ? 10?10 ? M? ?2 ? 10?7 M yr?1 were used, so as to study weak flashes and avoid regionsthat lead to novae or steady burning. The range of accretion rates chosenby the authors were motivated by previous investigations of shell flashes andnovae, which we will discuss shortly. To reiterate, the chosen range of ac-cretion rates do not lead to mass loss events, such as novae, but shell flashesand potentially mass accumulation. The authors went on to determine prop-erties of flashes, such as recurrence periods and peak layer temperatures. Inparticular, the authors found that as the hydrogen accretion rate decreasedfor a given white dwarf mass, the recurrence period as well as the peaktemperatures and violence of the eruptions increased. The authors also con-sidered helium flashes, both by direct and indirect accretion of helium. Inthe latter case, hydrogen burning results in the accumulation of helium.The authors found that recurrence periods were shorter for indirect accre-tion than for direct helium accretion. The recurrence periods were muchshorter for hydrogen flashes than for indirect helium flashes, with the firsthelium flash occurring much later than the first hydrogen flash.In 2007, Bildsten et al. derived the stability criterion for hydrogen-accreting white dwarfs by following the one-zone analysis of Paczyn?ski [35].These authors conducted a linear stability analysis. As in Paczyn?ski?s model,an explicit treatment of the flux produced from helium burning, accumu-lated as a product of hydrogen burning, was neglected. The authors founda very narrow range of accretion rates for which stable burning can occuras a result of radiation pressure stabilization. A one solar mass white dwarfcan stably accrete down to 1.7 ? 107 M yr?1 when the core luminosity isexcluded [35]. The authors also considered the effect of the ?-limited hotCNO cycle, which is of importance in accreting neutron stars as a result ofthe densities and temperatures reached within the layer. The authors found71.3. Shell Flashesthat including the hot CNO cycle with saturation increases the range of ac-cretion rates, leading to stable burning only for very low metallicities and istherefore, unlikely to be of significance in Type Ia progenitor scenarios [35].The limits and ranges of accretion rates leading to stable and unstableburning have been investigated by many other authors. We now present abrief overview of previous studies.Steady and Stable Burning Accretion RatesPaczyn?ski and Z?ytkow [29] investigated, through computations, the proper-ties of a 0.8 M white dwarf accreting hydrogen-rich matter with a solar com-position. The authors found that accretion rates ranging from 1.067? 10?7to 2.7?10?7 M yr?1 lead to stable burning in which the accumulating mat-ter is processed as fast as it is accreted. Higher accretion rates were foundto result in envelope expansion and the formation of a red giant containinga degenerate core [29]. In 1982, Fujimoto and Taam [13], who conducted ananalysis of helium flashes as a result of hydrogen burning in the context ofType Ia supernovae, found accretion rates leading to steady burning whichagreed with the results of Paczyn?ski and Z?ytkow.Nomoto confirmed that for a solar mass white dwarf, a narrow range ofaccretion rates 1?4? 10?7 M yr?1 lead to steady burning [27]. This limiton the accretion rate leading to stable burning was also confirmed by Iben,in which an analytical expression for the steady-burning regime accretionrate was found M?steady ? 1.32 ? 10?7(M/M)3.57M yr?1 [19]. Jose? etal. determined an upper bound on the accretion rate leading to flashes,M? = 4.5? 10?7 M yr?1 [20]. Later estimates confirm that large accretionrates are required to produce conditions in which steady burning can occur,in which matter is burned as fast as it is accreted. Townsley and Bildstendetermined that the lower limit for stable burning for a solar mass whitedwarf is 1.7 ? 10?7 M yr?1 [35]. The aforementioned results agree witheach other to within an order of magnitude. The discrepancies among theestimates for the accretion rate for steady and stable burning may possiblybe a result of differing assumptions and the variety of methods employed.Steady hydrogen burning is presumed to be responsible for supersoft X-ray sources, which emit soft X-rays with a luminosity close to the Eddingtonluminosity [17]. For larger accretion rates than those leading to steady burn-ing, approaching or exceeding the Eddington limit, further accretion may beinhibited as the outwardly directed pressure is larger than the gravitationalforce [17].81.3. Shell FlashesShell Flashes and Novae RatesAccretion rates that satisfy M? < 10?9 M yr?1 are expected to lead to novae[35], [23]. In particular, hydrogen accretion rates lower than 10?10?10?9M yr?1, with the exact value depending on the mass and luminosity of thewhite dwarf, result in novae [23], [20], [17].The upper limit, 10?9 M yr?1 was determined by Fujimoto in 1982 [11].Fujimoto used a semi-analytical approach to determine properties of shellflashes on accreting white dwarfs, and found that for accretion rates lowerthan 10?9 M yr?1 novae occur, as they are a consequence of slow accretiononto white dwarfs [11]. Fujimoto found that the mass of the white dwarfand the mass of the overlying layer uniquely determined properties of shellflashes. Fujimoto characterised his findings in terms of the pressure of theaccreted layer, calculated by multiplying the column density by the surfacegravity. In a companion paper, Fujimoto determined the stability of steady-state solutions through a perturbation analysis [12]. Fujimoto found that foraccretion rates in between two limits, approximately 10?6 and 10?5 M yr?1for a one solar mass white dwarf, stable burning in a steady-state results,which corroborates the rates in the discussion above [12].According to Nomoto, accretion rates below the threshold of steady hy-drogen burning, 1?4 ? 10?7 M yr?1, lead to flashes [27]. Furthermore,Nomoto found that below an accretion rate of 10?8 M yr?1 novae occurwith the violence of flashes increasing with decreasing accretion rate [27].Fujimoto and Taam identified hydrogen-rich accretion rates of M? <10?9 M yr?1 and white dwarfs of mass 0.7 M with classical novae [13].General Characteristics of FlashesAccretion rates smaller than those leading to steady burning but larger thanthe limit leading to novae ? 10?9 M yr?1, result in periodic shell flashesthat do not result in mass loss [17]. As the accretion is lowered below10?9 M yr?1, the flashes become more intense and eventually lead to nova-like bursts or hydrodynamic events accompanied by mass ejection [4]. Toreiterate, the intensity of flashes increases with decreasing accretion rate,until the threshold value is reached, and reducing the accretion rate furtherproduces a nova [17]. Accretion rates close to or larger than the Eddingtonlimit, result in envelope expansion and a red giant configuration [4]. Loweraccretion rates lead to a narrow range of steady burning. The lower limitof accretion rates leading to steady, stable burning decreases as the whitedwarf mass decreases [27], [35]. The ignition masses, or mass required to91.4. Motivation and Goalsaccumulate through accretion to trigger a nova is much smaller than thewhite dwarf mass [39]. The ignition mass decreases with increasing whitedwarf mass for a given accretion rate [27].1.4 Motivation and GoalsThe majority of stars in our Galaxy will terminate their lives as white dwarfs.As such, many authors have investigated the problem of hydrogen accretiononto the surface of a white dwarf. Novae, which are accompanied by massloss, directly eject hydrogen, helium or carbon into the interstellar medium.Stable burning regimes that allow the white dwarf to accumulate mass maybe progenitors of Type Ia supernovae. Both these events influence the chem-ical and dynamical evolution of our Galaxy. In this thesis, we will embark ona path that allows the identification of accretion rates that allow for stableor unstable burning of both hydrogen and helium.It is of particular importance to investigate the stable burning regimesthat may possibly identify progenitors of Type Ia supernovae. These eventsare quite significant for a number of reasons. In particular, the explosionof a CO white dwarf releases heavy elements, most notably iron [40], intothe interstellar medium, through a remnant or nebula. These heavy nucleiinfluence the element abundance of the host galaxy [40]. Furthermore, as aresult of the uniformity in luminosity of these events, they can be used aspowerful observational tools to determine cosmological distances and param-eters [32], [40]. The 2011 Nobel Prize in Physics was awarded to researcherswho, through the observation of these supernovae, discovered the accelerat-ing expansion of the Universe [32]. Consequently, it is an important problemto understand the progenitors of these events, so as to constrain both thedynamical and chemical evolution of the Universe.Whilst it is clear that stable burning may lead to mass increase, thereis an alternate scenario that may lead to mass accumulation as well. RSOphiuchi is a recurrent nova, accreting from a red giant, which has eruptedseveral times in the past [17]. The last eruption took place in 2006, 21 yearsafter the previous eruption [18]. The outbursts are presumed to be relatedto nuclear burning of accreted hydrogen-rich material. However, these novaeseem to have the property that more mass is accreted than ejected duringa nova, resulting in mass increase [17]. Theoretical models have reproducedsalient features of the outbursts indicating an accretion rate in the 10?8 to2?10?7 M yr?1 range [18]. The white dwarf increases in mass by 10?6 Mafter every eruption [17].101.5. Outline and Description of ModelIt is evident that an understanding of the accretion rates resulting invarious types of burning is still an important problem and one which wetackle in this thesis. It is worth noting that there are two distinct method-ologies that can be applied to this problem. Certainly, one can follow theevolution of these scenarios with computationally intensive, one, two andthree-dimensional simulations. These simulations are able to describe thephysics in great detail. For example, the FLASH code was developed tostudy shell flashes on accreting compact stars [10]. A sophisticated stellarevolution code, Modules for Experiments in Stellar Astrophysics or MESA, isable to tackle a wide range of problems [30]. MESA star, a one-dimensionalevolution module, is able to investigate classical novae and Type Ia super-novae [30]. While this approach is necessary to produce a comprehensivedescription of the physics at work, it is time-consuming and examining alarge parameter regime is subsequently difficult.The alternate approach, and the approach which we will employ in thisthesis, is based on a semi-analytical method. The advantage is that thisapproach is not particularly time-consuming but allows for the identifica-tion of parameter regimes which lead to unstable and stable burning. Inparticular, the dependence of burning on various aspects, such as surfacegravity, chemical composition and accretion rate can be explored.1.5 Outline and Description of ModelPaczyn?ski?s pioneering work set the stage for the investigation of stabilitythrough a perturbation analysis of steady-state solutions [28]. Indeed, muchof the literature that deals with such an analysis on accreting compact starsutilises his work as a reference. It is for this reason that we begin ouranalysis by independently producing and confirming the results of Paczyn?skiin detail, as a first step to investigating problems of hydrogen-rich accretionon the surface of white dwarfs. Paczyn?ski?s work also forms the basis ofour subsequent analysis and a number of important assumptions and inputphysics arise in the context of Paczyn?ski?s model.We then extend Paczyn?ski?s one-zone model to a multi-zone model incor-porating hydrogen and helium burning. The model can also be consideredas an extension of Jose? et al. [20]. However, we do not explicitly zone the hy-drogen and helium layers, but allow them to evolve naturally. Our analysiswill closely follow the work of Narayan and Heyl [26], in which the authorsconducted an analysis of the thermonuclear stability of hydrogen-rich ma-terial accreting onto the surface of a neutron star. The authors determined111.5. Outline and Description of Modelsteady-state solutions, and conducted a global linear stability analysis todetermine the stability of the steady-state solutions. The authors also wenton to identify conditions that trigger bursts on neutron stars, and identifiedvarious regimes of bursting behaviour. In an earlier paper, the same authorsconducted a similar linear stability analysis of accreting matter on the sur-face of a neutron star [25]. However, the steady-state solutions were foundthrough an alternate application of a boundary condition at the base of theaccreted layer [25]. We adapt their approach, outlined in [25] and [26], soas to be suitable in the context of white dwarfs.In our analysis, the layer of accumulated material is studied in a quasi-static approximation. We will focus on finding steady-state solutions andthereafter present a framework, a linear perturbation analysis, with whichto determine the stability of these solutions. We aim to identify regimes ofstable and unstable burning by accretion rate and flux from the surface. Wewill keep certain parameters such as the mass of the white dwarf and thecore temperature fixed. The identification of possible stable and unstableburning parameter regimes may serve as input for detailed computationalsimulations. Furthermore, our analysis produces results that can potentiallybe tested by experiment.We begin this thesis with a discussion, contained in chapter 2, of theinput physics that is relevant to the analysis of accreting white dwarfs. Inchapter 3, we verify the results of Paczyn?ski?s one-zone model by followinghis approach in detail. We then present an extension of the one-zone model, amulti-zone model in chapter 4. We determine the steady-state solutions andattempt to physically interpret the results and compare the results to thoseof other authors. This is followed by the presentation of a framework thatcan be used to determine the stability of steady-state solutions numericallyin chapter 5. Finally, we conclude with a discussion of our results, potentialobservational confirmations of the theory and future work.12Chapter 2Preliminaries2.1 Model FrameworkThe framework we will use to model the accreted layer consists of a numberof coupled differential equations. The standard stellar structure equations,namely the hydrostatic equilibrium, radiative transfer and energy equations,form the basis of the model. These equations describe the thermal andstructural conditions that must be satisfied by the layer for our model to bephysically correct [16]. The hydrostatic equilibrium equation enforces thecondition that the gravitational force must be balanced by the outwardlydirected pressure [21]. If this condition were not satisfied, the star andthe accreted layer would not be bound. The radiative transfer equationindicates how energy is transferred in the system, and the sources whichhinder the flow of radiation [16]. The opacity, ?, specifies the processesthat obstruct the path of photons in the layer [5]. The conservation ofenergy equation specifies how energy is produced and lost in the system[16]. We will consider two sources of energy generation, , namely nuclearand gravitational sources. The accreted layer is hydrogen-rich, and thereforehelium and other heavy elements will be produced through nuclear burning.Thus, our model will consist of equations that describe the evolution ofhydrogen and helium. We will also need to specify the equation of state andentropy in our model. Thus, the equations that describe the evolution ofthe accreted layer contain a variety of physics. We will discuss these aspectsin this chapter to allow for an easier reading and understanding of materialcontained within subsequent chapters.2.2 Mean Molecular WeightGases in stellar modelling are mixtures of atoms, ions and electrons. Weconsider each of these components as free particles. The ion mean molecularweight, ?I, in which we include neutral atoms, describes the average mass ofan average ion in the gas [16]. In a mixture, we assign to each species, i, of132.2. Mean Molecular Weightions a mass fraction number, Xi, such that?iXi = 1 [16]. For example, if60% of a gas is composed of hydrogen, then Xi = 0.6. Typically, we denotethe hydrogen mass fraction by X, the helium mass fraction by Y and heavierelements by Z, such that X + Y + Z = 1.Each ion or atom is associated with a charge number, Zi, indicating thenumber of protons, and a mass number in atomic mass units, Ai, giving thetotal mass [16]. As an example, carbon has a charge number of Z = 6 anda mass number of A ? 12. We define the total number density of ions as [5]nI = ?NA?iXiAi, (2.1)where NA = 6.022?1023 per mole is Avogadro?s number, and ? is the densityof the mixture. We can define the ion mean molecular weight as [5]?I =(?iXiAi)?1. (2.2)To define an electron number density and mean molecular weight, weneed to know the ionisation fraction, yi, which tells us how many free elec-trons are associated with a particular species, i. The number density ofelectrons is defined by [16]ne =?NA?e, (2.3)where the free electron mean molecular weight is [16]?e =(?iZiXiyiAi)?1. (2.4)The total mean molecular weight is [16]? =( 1?I+ 1?e)?1. (2.5)In the accreted layer that we will consider, as well as within the interiorof a white dwarf, we expect the atoms or ions to be fully ionised so thatyi = 1. If this were not the case it would be quite a difficult task to calculatethe ionisation fractions. In stellar regions that are completely ionised, andcomposed mainly of hydrogen and helium, we can approximate the electronmean molecular weight as?e =21 +X , (2.6)142.3. Energy Generation Rateswhere we have approximated the nuclear charge of major heavy elements asZ ? A/2 [5].We will make use of the following expression for the total mean molecularweight in regions of complete ionisation [5]? = 21 + 3X + 0.5Y . (2.7)Equation (2.7) approximates the average atomic weight of a heavy elementas AZ = 2Z+2 [5]. This is an adequate approximation in the accreted layerwhere the initial heavy element fraction is very small.2.3 Energy Generation RatesFor a star to produce luminosity, or a power output, sources of energy mustexist. We will consider energy generation from both nuclear sources andgravitational sources.2.3.1 Nuclear Energy SourcesThe accreted material we consider in this thesis is hydrogen-rich. As such,we need to describe the processes that liberate energy through hydrogen fu-sion and other relevant reactions. Stars burn hydrogen into helium throughtwo major reaction cycles, namely the proton-proton chains, and the CNO(Carbon Nitrogen Oxygen) cycle [5]. The importance of each of these reac-tion chains in our model depends largely on the temperature of the accretedmaterial.The proton-proton chains consist of three major cycles that are respon-sible for converting hydrogen to helium in main sequence stars [16]. Theeffective energy generation rate for the proton-proton chains is [16]pp = 2.4? 106?X2T?2/36 exp(?33.8T?1/36)erg g?1 s?1, (2.8)where T6 = T/106 K.The CNO cycle is also responsible for converting hydrogen into heliumbut at higher temperatures than the proton-proton cycle [16]. As a result,this reaction cycle usually occurs in more massive main sequence stars. TheCNO cycle can be split into two cycles. The first cycle is known as the CNcycle, which requires carbon and nitrogen as catalysts. We will refer to these152.3. Energy Generation Ratesheavy element catalysts as CN elements. The energy generation rate of theCN cycle, using equilibrium abundances, is [5], [28]CN = 8? 1027?XXCNT?2/36 exp(?152.31T?1/36)erg g?1 s?1, (2.9)where XCN is the CN element mass fraction. This rate can also be used forthe CNO cycle provided we use XCNO instead of XCN. The estimation ofXCN or XCNO requires careful consideration, as it depends on the numberdensity of 14N [16]. In equilibrium, the number density of 14N is the sumof the number densities of the CN or CNO isotopes before burning. It iscrucial to decide on which cycle, CN or CNO, should be used to represent14N in the energy generation rate [16]. A detailed calculation of XCN isbeyond the scope of this thesis. Paczyn?ski set XCN = Z/3 in his one-zonemodel [28]. However, from stellar models dealing with solar-like stars, onecan conclude that it is also reasonable to use XCN = Z/2, which introducesan error of approximately 25% [16].To discern which one of these energy generation rates, CN or pp, willbe most important for our model, we consider the temperature dependenceof each of these rates. It is common to express the energy generation ratesas power-laws of the form [16] = 0??T ? , (2.10)where 0, ? and ? are constants to be determined over a restricted range oftemperature and density. We are predominantly concerned with the tem-perature exponent ?, which can be deduced as follows? =( ? ln ? lnT6). (2.11)From expression (2.8) for the proton-proton chains we find that?pp = ?23 + 11.27T?1/36 , (2.12)and from equation (2.9) for the CN cycle we find that?CN = ?23 + 50.77T?1/36 . (2.13)The temperatures that are relevant to our problem are generally above2?107 K. Substituting this into the equations above we find that ?pp = 3.49and ?CN = 18.04. Thus, the energy generation rate of the CN cycle is162.3. Energy Generation Ratesmuch more sensitive to changes in temperature in the regimes that we areinterested in. However, at lower temperatures, the proton-proton chainsbecome the dominant energy source. Figure (2.1) depicts the ratio of thetwo rates for a range of temperatures on a logarithmic scale. In this plot,we set X = 0.7 and Z = 0.03. The CN cycle dominates for temperatureslarger than approximately 107.5 K, while the proton-proton chain dominatesthe energy generation at temperatures lower than this threshold.While the CN cycle is sufficient for an analysis of burning on whitedwarfs, it is not sufficient for burning on neutron stars. The density andtemperatures will be much higher in the accreted layer due to the surfacegravity, which is approximately seven orders of magnitude larger for a neu-tron star than for a white dwarf. In this case, the hot CNO cycle withsaturation must be used for the energy generation rate. The energy gener-ation rate can be written as [26]CNO = 4E?HrCNOXCNO14 , (2.14)where E?H is the energy liberated from the burning of one gram of hydrogen,XCNO is the mass fraction of the CNO elements and rCNO is the ?rate ofreactions per CNO nucleus? [26]. The following results for rCNO, ?13 and ?14assumes equilibrium abundances for all the nuclear species, along with thefurther assumption that 14N constitutes the majority of these equilibriumabundances [5]. We state the results for rCNO, ?13 and ?14 [26], [24]:rCNO =( 1?13+ 1862)?1 ( 1?13( 1?13 + ?14 + 278.2)+ 1862( 1?14 + 1038)),(2.15)?13 = (X?)?1(3.35? 107T?2/39 ? exp(?15.202T?1/39 ? 0.8702T 29)?(1 + 0.027T 1/39 + 0.9T2/39 + 0.173T9 + 4.61T4/39 + 2.26T5/39)+3.03 ?105T?2/39 exp(?6.348T?19))?1, (2.16)?14 = 3.1? 1010 (X?)?1 T 2/36 exp(152.313T?1/36), (2.17)where T9 = T/109 K. The details of this reaction rate can be found in [26].If temperatures reach 108 K, helium will undergo burning, through thetriple-? reaction, to produce carbon and oxygen. The energy generationrate is [16]3? = 3.5? 1011Y 3?2T?38 exp(?43.2/T8), (2.18)172.3. Energy Generation Ratesand the temperature exponent ?3? is?3? = ?3 + 43.2T?18 . (2.19)For a temperature of 108 K, ?3? = 40.2. This is much larger than thecorresponding exponent for the CN cycle at this temperature, ?CN = 8.01.This implies that helium burning can be quite explosive when a thermalinstability cannot be quenched.2.3.2 Electron CaptureElectron capture reactions occur when an atomic electron is captured by aproton within the nucleus [22],e? + p ?? n+ ?e. (2.20)For the reaction to be energetically favourable, the parent nucleus must begreater in mass than the combined mass of the daughter nucleus and bindingenergy of the electron, Be, or mathematically, [22]Q = (mparent ?mdaughter)c2 ?Be > 0. (2.21)Consequently, the Q-value for the reaction is positive, and the reaction willbe exothermic. If densities exceed 107 g cm?2 in the accreted layer, electroncapture reactions will occur on hydrogen and energy will be liberated [2].Due to the surface gravity of the white dwarf, the density will rarelyreach this threshold in the accreted layer. However, this scenario is likely inthe accreted layer on a neutron star. We state the electron capture rate [2]Rec =ln 21065I(EF, T ), (2.22)where EF is the Fermi energy and T is the temperature. Here I(EF, T ) isan integral that has been evaluated explicitly by Bildsten and Cumming [2],with the following resultI ? Q2(mec2)5????2(kBT )3 exp(EF?QkBT)EF < Q+ kBT ln(92),13(EF +(3? ln 92)kBT ?Q)3EF > Q+ kBT ln(92),(2.23)where Q = 1293.318 keV is the Q-value for the reaction and me is the elec-tron mass. Note that the integral result is an approximation which has been182.3. Energy Generation Ratescalculated using certain limits to attain the analytic fitting formula above.To calculate the electron capture rate accurately, one would have to calcu-late the complete phase space integral over the electron energy, I(EF, T ).The fitting formula agrees with the rate given by the full integral to withina factor of two [2]. The energy generation rate can be written asec = 2E?HRec. (2.24)The factor of two, and, in fact, the contribution of electron capture to theenergy generation rate, is inconsequential for accreting white dwarfs, as aresult of the densities reached in the accreted layer. The Fermi energy canbe calculated provided that the electron number density is known. Sincetemperatures in the accreted layer are well above 104 K, we can safely assumethat the elements are fully ionised [16]. The Fermi energy can be written as[34]EF = mec2(1 + x2F)1/2, (2.25)where xF = pF/(mec) and pF is the Fermi momentum associated with theFermi energy. To calculate xF we use the following relationship [16]:x3FBNA =?NA?e, (2.26)where B = 9.739 ? 105 g cm?3 for electrons and ?e is given by equation(2.6). Using this relationship, we find thatxF =(3pi2 ?NA(1 +X)2)1/3B?, (2.27)where B? =(3pi2BNA)?1/3 = 3.8617? 10?11 g?1/3 mole1/3 cm. Substitutionof this result into (2.25) allows the calculation of the Fermi energy and hencethe electron capture energy generation rate. At the threshold between therelativistic and non-relativistic regimes, xF ? 1 [16].2.3.3 Gravitational Energy SourcesEnergy release from gravitational sources arise as a result of a departurefrom adiabatic behaviour during contraction or expansion. In our model,gas is compressed as material accretes onto the star. Energy is releasedwhen compression is accompanied by a rise in pressure which is less than192.4. Opacityadiabatic [16]. In this case, the energy generate rate, grav, due to gravitywill be greater than zero. We can express this mathematically asgrav = ?Tdsdt > 0, (2.28)where s is the specific entropy and ddt is the Lagrangian derivative [16].2.4 OpacityOpacity describes the resistance encountered by the flow of radiation. Inother words, the opacity affects how much radiation actually reaches thesurface of the star. Generally, to adequately describe the opacity, one mustinclude all processes that hinder the motion of a photon through the stellarmedium. There are two main sources of opacity, namely, radiative andconductive opacity sources.2.4.1 Radiative OpacitiesWithin a star, radiative opacity is described very well by four processes,namely bound-free, free-free and bound-bound absorption and lastly, scat-tering from free electrons [5]. In each of the first three cases, the opacitycan be approximated by Kramers? opacity, which is a power-law of the form?(?, T ) = ?0?T?3.5, where ? is the opacity depending on density, ?, temper-ature, T , and a constant, ?0 [5].Electron Scattering OpacityPhotons are scattered by free electrons in the stellar gas in a process usu-ally known as Compton scattering. In the non-relativistic regime, which istypical of stars, it is known as Thomson scattering [5]. In areas where thegas is almost completely ionised, the opacity produced by the scattering ofphotons from free electrons is important. The electron scattering opacity is[5]?e = 0.2(1 +X), (2.29)where X is the hydrogen mass fraction.The ionisation temperature of hydrogen is approximately 104 K [16]. Be-low this temperature, there are fewer free electrons left in the mixture. As aresult, the electron scattering opacity decreases rapidly from the value givenby the expression (2.29) for temperatures below the ionisation temperature.202.4. OpacityBound-Free AbsorptionIn this process, a photon is absorbed by a bound electron. The photon hasenough energy to ionise or remove the electron from the atom or ion, andhence the process is also referred to as photo-ionisation [5]. The opacity canbe written in Kramers? form as [16]?bf = 4? 1025Z(1 +X)?T?3.5, (2.30)where Z is the heavy element mass fraction. Bound-free absorption becomesan important source of opacity when temperatures exceed 104 K. It is inthis temperature regime that photons within the gas have enough energy toionise electrons.Free-Free AbsorptionThis is the absorption of a photon by a free electron in the vicinity of anion [5]. If the ion were not present, this process would be prohibited byconservation of energy [16]. The inverse of this process is Bremsstrahlung,in which an electron moving past an ion results in the emission of a photon.Once again, the opacity can be written in Kramers? form as [16]?ff = 4? 1022(X + Y )(1 +X)?T?3.5. (2.31)Bound-Bound AbsorptionAn atom or ion absorbs a photon, which results in a bound electron makinga transition to a higher energy bound state [5]. The opacity produced bythis process can be written in Kramers? form. However, the opacity fromfree-free and bound-free transitions is typically an order of magnitude largerthen the bound-bound opacity, and we omit the result here [16].Radiative OpacityWe use a radiative opacity that is the sum of electron scattering opacityand Kramers? opacity for bound-free, free-free and bound-bound transitions[28]:?rad = ?e(1 + 2? 1026(0.001 + Z)?T?3.5). (2.32)212.5. Pressure2.4.2 Conductive OpacityHeat can also flow through conduction, and thus we can define a conductiveopacity. Degenerate electrons support the white dwarf core against gravity.These degenerate electrons also facilitate energy transport through electronheat conduction across a temperature gradient. We can define the conduc-tive opacity as [16]?cond =4acT 33De?, (2.33)where a is the radiation density constant, c is the speed of light and De isa diffusion coefficient. In strongly degenerate, non-relativistic regions, thediffusion coefficient is [5]De = 2.36? 103?T?eZ?, (2.34)where ? is a slowly varying function, depending on electron velocity, thatdescribes Debye screening effects. In our approximations, we set ? = 1,thereby neglecting screening effects.2.4.3 Total OpacityThe total energy flux is the sum of the energy flux produced by radiationand conduction. In our analysis, we neglect convective effects. The totalopacity can be calculated as follows [5]1?tot= 1?rad+ 1?cond. (2.35)If ?rad  ?cond, then ?tot ? ?rad. If If ?cond  ?rad, then ?tot ? ?cond. Innon-degenerate stellar material, radiative opacity dominates. However, indegenerate material, conductive opacity usually dominates the total opacity.2.5 PressureThe pressure we use to describe the accreted material, and the inner regionsof the white dwarf, is the sum of the ideal gas law, radiation pressure, anddegeneracy pressure:P = Pgas + Prad + Pdeg. (2.36)We use the total mean molecular weight, ?, given by equation (2.7) in theideal gas lawPgas =kBNA?T? . (2.37)222.5. PressureThe radiation pressure is given byPrad =a3T4, (2.38)where a is the radiation density constant.The degeneracy pressure for a zero temperature electron gas can be cal-culated in both the relativistic and non-relativistic regimes through the fol-lowing integral [34],Pdeg = 8(pim4ec53h3)? xF0x41 + x21/2= A?f(xF), (2.39)where A? = 24pi2 pim4ec53h3 = 1.42180? 1025 dyne cm?2 and f(x) is given by [34]f(x) = 18pi2(x(1 + x2)1/2(23x2 ? 1)+ ln(x+ (1 + x2)1/2)). (2.40)Here xF = pF/(mec) and pF is the Fermi momentum associated with theFermi energy EF. The quantity xF can be calculated, provided the numberdensity of the electrons is known, through equation (2.27), which assumescomplete ionisation. Note that ln(x+ (1 + x2)1/2)= sinh?1 x.In the non-relativistic case, xF  1, the electron degeneracy pressure is[16],Pdeg = 1.004? 1013( ??e)5/3, (2.41)where ?e, the electron mean molecular weight, is given by equation (2.6).Upon substitution of ?e = 21+X , we findPdeg = 3.16? 1012 (1 +X)5/3 ?5/3. (2.42)The equation of state is decoupled from the temperature. For convenience,we define a coefficient of density, K, asK = 3.16? 1012 (1 +X)5/3 . (2.43)Paczyn?ski uses a slight variation of the parameter K [28]:K = 3.12? 1012 (1 +X)5/3 . (2.44)The expression for pressure that we use in our analysis isP = kBNA?T? +a3T4 +K?5/3. (2.45)232.6. EntropyWe note that we do not include magnetic pressure, as a consequenceof the fact that our model does not extend to magnetic white dwarfs. Inaddition, in our model, the pressure given by equation (2.45) would be largerin magnitude than magnetic pressure, if we were to include weak magneticwhite dwarfs.2.6 EntropyThe specific internal energy is given by [28]u = 3kBNAT2? +3K?2/32 +aT 4? erg g?1, (2.46)which is the sum of the specific energy of an ideal gas, a non-relativisticdegenerate electron gas and radiation pressure.The fundamental thermodynamic relation for infinitesimal and reversibleprocesses isdQ = Tds = du+ Pd(1?), (2.47)where Q is the heat added to or removed from the system per unit mass,du is the change in internal specific energy, Pd(1?)is the work done bythe system on its surrounding if the specific volume changes, and s is thespecific entropy. All the quantities in equation (2.47) are specific quantities,or quantities per unit mass.Since u is a function of density and temperature, we can write the fun-damental thermodynamic relation asTds = ?u?T dT +?u??d??P?2d?, (2.48)where P is given by equation (2.45).We calculate the partial derivatives of u from equation (2.46)?u?T =3kBNA2? +4aT 3? , (2.49)?u?? = K??1/3 ? aT4?2 . (2.50)Substituting these into equation (2.48) and simplifying, we find thatTds =(3kBNA2? +4aT 3?)dT +(?4aT43?2 ?kBNAT??)d?. (2.51)242.7. Eddington LimitThe specific entropy is a function of density and temperature, s(?, T ). Fromdifferential calculus, we can express ds asds = ?s?T dT +?s??d?. (2.52)Using this, we find that?s?T =3kBNA2?T +4aT 2? , (2.53)?s?? = ?4aT 33?2 ?kBNA?? . (2.54)Integrating these equations we finds = 3kBNA2? ln(T ) +4aT 33? + f(?), (2.55)s = 4aT33? ?kBNA? ln(?) + g(T ). (2.56)Combining these two expressions, we find that the specific entropy, s, iss = 4aT33? +kBNA? ln(T 3/2?), (2.57)to within an additive constant.2.7 Eddington LimitThe Eddington luminosity is [35]LEdd =4picGm? , (2.58)where c is the speed of light, G is the gravitational constant, m is the massof the star. The Eddington limit provides a threshold luminosity such thatif it is exceeded, radiative forces exceed the gravitational forces that keepthe star bound. Surpassing this limit is connected with mass loss [34]. Thislimit is approached when the luminosity, and therefore temperature, is veryhigh. Consequently, the opacity used in equation (2.58) is that of electronscattering opacity, ? = ?e.252.8. Lagrangian Derivative2.8 Lagrangian DerivativeThe Lagrangian or material derivative, which follows a particular parcel ofgas, is [6]ddt =??t + v ? ?, (2.59)where v is the fluid velocity and ? is the gradient operator. In the Euleriandescription, the time derivative is equal to the partial time derivative in theLagrangian description,(ddt)Eulerian=(??t)Lagrangian[6].2.9 Core Temperature, Mass and Surface GravityTo conduct the multi-zone analysis presented in chapter 4, we will needvalues for the mass and core temperature of a white dwarf found in a binarysystem. In light of the evidence presented in the introductory chapter, wechoose a core temperature of 4 ? 107 K. Whilst slightly hotter than thetypical central temperature of an isolated white dwarf, it seems to be areasonable choice, as we will be considering white dwarfs in mass-transferringbinaries. In our calculations we use a 0.929 M white dwarf, which was usedin Paczyn?ski?s analysis [28]. The higher mass is substantiated by the factthat the white dwarf is in a binary system, and fits within the range ofwhite dwarf masses chosen by the authors in the literature review. Withthis mass, and a radius of approximately 0.01 R, which is typical of a whitedwarf, the surface gravity is approximately 108.5 cm s?2 [28]. This value isslightly higher than the average surface gravity of an isolated white dwarf108 cm s?2, as a result of the larger mass.2.10 Thermal InstabilityThere are two types of stability that we can investigate. Firstly, dynamicalstability probes properties of a star when hydrostatic equilibrium is per-turbed. Thermal, or secular stability, investigates the effects of perturbingthermal equilibrium. The focus of this thesis is thermal stability. In partic-ular, we will investigate the effects of perturbing thermal equilibrium in theaccreted layer on accreting white dwarfs. There will be stable regimes, wherethe perturbation is quenched, and unstable regimes, where the perturbationgrows exponentially with time.If we exclude gravitational energy sources, the energy source within astar is solely attributable to nuclear energy generation. For a star in ther-262.10. Thermal Instabilitymal equilibrium, the rate of energy production must be balanced by therate of energy loss by radiative transfer, excluding conductive effects. Inother words, the rate of change of energy within the star should be zero, ormathematically [21]dUdt = Lnuc ? Lrad = 0, (2.60)where U is the total energy. This equation implies that the rate of energyproduced, or luminosity, within the star, Lnuc , is balanced by the powerleaving the star, Lrad, at radius, r. Thus, the total energy remains constant.If U? > 0, the total energy increases and thermal equilibrium is perturbed.The relative importance of the various components of the pressure, namely,the ideal gas, radiation and degeneracy pressure, will determine if the systemis able to restore thermal equilibrium. In a system dominated by the idealgas law, thermal equilibrium will be restored. The virial theorem for a sys-tem with the ideal gas law as the equation of state, states that 2E + ? = 0,where E is the internal energy and ? is the gravitational potential energy[21]. The total energy satisfies U = ?E = ?/2. For a bound star, the totalenergy is less than zero, U < 0, as the gravitational potential energy is lessthan zero, ? < 0. If thermal equilibrium is perturbed, the total energy willincrease, which means that it becomes less negative, or smaller in absolutevalue. The gravitational potential energy also becomes less negative, whichsignals an expansion. By the virial theorem, if U and ? increase, the internalenergy must decrease. Thus, an increase in total energy results in expansionand a decrease in temperature, restoring thermal equilibrium [21]. In thecase of a non-relativistic degenerate electron gas, the pressure is decoupledfrom the temperature in a constant density star. If U? > 0, the star willstill expand, but the expansion is accompanied by an increase in tempera-ture [21]. Since the nuclear energy generation rates are proportional to ?T ?in power-law form, an increase in temperature results in further increasednuclear energy generation and a further increase in temperature and so on.Since the temperature in the energy generation rates are usually raised tovery high powers, an increase in temperature initiates a thermonuclear run-away event. If the temperature becomes sufficiently high for the gas tobehave as an ideal gas, thermal equilibrium can be restored.Thermonuclear runaways occur on the surface of accreting white dwarfsthrough ignition of the accreted hydrogen-rich material. This nuclear ig-nition gives rise to hydrodynamic flashes, more commonly known as novaein white dwarfs. To reiterate, we will investigate the stability of novae orflashes in this thesis. In particular, we intend to identify regimes of stableand unstable nuclear burning. Stable regimes do not lead to thermonuclear272.10. Thermal Instabilityrunaways. In stable regimes, the perturbations can be quenched and ther-mal equilibrium is restored. In unstable regimes, nuclear reactions lead tothermonuclear runaways and novae.We have now discussed all the relevant physics that will be required todescribe the nuclear burning of accreted material on the surface a whitedwarf. In the next chapter we investigate the problem using a one-zonemodel and eventually progress to a more sophisticated multi-zone model.282.10. Thermal Instability4 5 6 7 8 9?250?200?150?100?50050log T (K)log(? CN/? pp)Figure 2.1: Ratio of CN to pp as a function of temperature, with X =0.7 and Z = 0.03. The red line indicates the point at which the energygeneration rates are equal.29Chapter 3A One-Zone ModelPaczyn?ski studied the problem of shell flashes on accreting compact stars.In 1983, he published a paper outlining a one-zone model that described thebasic properties of accreting white dwarfs and neutron stars [28]. We followhis analysis independently and confirm his results, as a starting point forour investigation of the properties of accreting white dwarfs.3.1 Model EquationsAs previously discussed, the timescale for evolution, or cooling, of a whitedwarf core is on the order of billions of years [1], [33]. The timescale forthe evolution of the accreted layer of hydrogen-rich material will be muchsmaller than the timescale for cooling, on the order of thousands to hundredsof thousands of years [36], [37]. As such, we can treat the white dwarf coreas if it were in a steady-state compared with the accreted layer. FollowingPaczyn?ski?s approach, we treat the accreted layer as a thin zone of hydrogen-rich material overlying the inert white dwarf core. As a result, we assumethat the height of the one-zone layer is much smaller than the radius of thewhite dwarf, and that the layer contains very little mass in comparison tothe mass of the white dwarf. Consequently, we can fix the surface gravityin the layer. Accordingly, we use a plane parallel geometry to describe thelayer. This model does not include magnetic affects, and assumes a negligiblemagnetic field associated with the white dwarf, such that magnetic pressurecan be ignored.We specify the independent variable as ?, the column density of thelayer in mass per unit area,???z = ?, (3.1)where z is the height of the layer and ? is the density.The standard stellar structure equations, which include the structuraland thermal conditions required to describe the layer, along with an equationfor the evolution of the hydrogen-mass fraction, constitute the set of model303.1. Model Equationsequations, which are [16], [28]?P?? = ?g, (3.2)?Prad?? = ??Fc , (3.3)?F?? = H ? Tdsdt , (3.4)dXdt = ?HE?H, (3.5)where P is the pressure, g is the gravitational acceleration, Prad is the radia-tion pressure, ? is the opacity, F is the flux, c is the speed of light, H is theenergy generation rate of hydrogen burning, s is the specific entropy, andT is the temperature. The time derivatives here are Lagrangian derivativesand E?H is the energy released by burning one gram of hydrogen, approxi-mately 6 ? 1018 ergs cm?2 s?1. In these equations, we integrate from thebottom of the layer to the surface in ?.The first equation describes the condition of hydrostatic equilibrium,indicating that the outwardly directed pressure is balanced by the inwardlydirected gravitational force. The second equation, the radiative transferequation, indicates that energy is transferred through diffusion. This meansthat the flux is driven by gradients in the energy density of the radiation field[16], and accordingly, the opacity is the radiative opacity ? = ?rad given byequation (2.32). The equation for the flux, F , is the energy balance equation.There are two sources of energy generation indicated in the equation forthe flux, namely the energy generated per gram per second by hydrogenburning, H, and gravitational compression, expressed as grav = ?T dsdt > 0.Finally, the last equation is essentially a conservation of mass equation forthe hydrogen mass fraction.To solve this set of equations, we must specify the appropriate boundaryconditions. At the surface of the layer, we specify the following conditions[28]:? = ?s, P ? 0, Prad ? 0, X = Xs, Y = Ys, (3.6)where Xs and Ys specify the composition of the accreted material. In otherwords, Xs and Ys specify the hydrogen and helium composition of the donorstar. These boundary conditions are consistent with zero surface boundaryconditions. At the bottom of the layer, we have the following conditions onthe flux and hydrogen mass fraction [28]:? = ?b, F = Fb, X = 0, (3.7)313.1. Model Equationswhere Fb is the flux from the bottom of the layer, which is the flux emanatingfrom the white dwarf core. The X = 0 boundary condition implies that allthe hydrogen must be burned upon reaching the bottom of the layer. Thecolumn density of the accreted layer is ?? = ?s ? ?b [28].The accretion rate, ??, in g cm?2 s?1, can be related to the surface columndensity through [28]d?sdt = ??. (3.8)The system of equations, equations (3.2) to (3.5), must be integratedover the column density of the layer. However, in Paczyn?ski?s analysis,several simplifications were made to circumvent the numerical solution ofthe system of differential equations. Specifically, Paczyn?ski carried out theintegration analytically by replacing the integrands by their values at theeither the top or bottom of the layer, and arrived at the following set ofequations [28]:P = g??, (3.9)Prad =?c F??, (3.10)F = Fb +(H ? Tdsdt)??, (3.11)Xd(??)dt ?X?? = ?HE?H??, (3.12)where P, Prad, T, ?, s, ? and H refer to the values of the associated variablesat the bottom of the layer, while X and F refer to quantities at the topof the layer. The Lagrangian derivatives have been replaced here by totaltime derivatives. This approximation neglects the heat and matter that iscarried to the bottom of the layer. Despite the simplifications, this modelproduces well-known characteristics of accreting compact stars found bydetailed stellar models [28].By some rearrangement we can derive two differential equations thatencapsulate the details contained within the equations above. Firstly, if wesubstitute ?? = Pg into equation (3.11), we find that [28]T dsdt = H ? (F ? Fb)gP . (3.13)Similarly, if we substitute ?? = Pg into the hydrogen mass fraction equation,equation (3.12), divide by X and multiply by g, we find that [28]dPdt = g???HPXE?H. (3.14)323.1. Model EquationsTo solve this set of differential equations numerically, the accretion rate, ??,the flux from the white dwarf core, Fb, the surface gravity, g, the chemi-cal composition at the surface, Xs, Ys, the energy generation rate, H, theopacity law, ?rad, and the equation of state must be specified. The energygeneration rate will be that of the CN cycle as specified by equation (2.9).We will use the radiative opacity given by equation (2.32), and finally theequation of state will be the sum of the ideal gas, radiation and degeneracypressures, which are specified in equation (2.45). The corresponding specificentropy is given by equation (2.57). The entropy, opacity, equation of stateand energy generation rate will be functions of density and temperature aswell as the hydrogen and helium mass fractions.As is customary, we define 1 ? ? as the ratio of radiation pressure tototal pressure. By using equations (3.9) and (3.10) we find that1? ? = PradP =??e?eFcg . (3.15)If ? ? 0 then radiation pressure dominates. Note that in this stellar model,the flux can approach the Eddington limit but cannot exceed it since 0 <? < 1.3.1.1 Dimensionless QuantitiesBefore proceeding, we introduce some dimensionless quantities that makethe forthcoming analysis clearer. We introduce the Eddington fluxFEdd =cg?e, (3.16)where ?e is the electron scattering opacity, given by equation (2.29). Thisequation is found by dividing the Eddington luminosity, equation (2.58),by the surface area of the star. We define a critical accretion rate [28], inaccordance with Paczyn?ski,??c =FEddE?HX. (3.17)This is essentially the nuclear Eddington accretion rate for total hydrogenburning [35]. The dimensionless accretion rate in units of the nuclear Ed-dington accretion rate isa? = ????c. (3.18)333.2. Steady-State SolutionsWe also define a dimensionless heat flux from the white dwarf core in unitsof the Eddington flux, fb,fb =FbFEdd. (3.19)We can find an equation for the flux at the surface by combining equa-tions (3.15) and (3.16),F = FEdd?e? (1? ?) . (3.20)By using this equation for the flux and the dimensionless quantities, we canwrite equations (3.13) and (3.14) as [28]T dsdt = H +cg2?eP(fb ??e? (1? ?)), (3.21)E?HXPdPdt = ?H +cg2?ePa?. (3.22)3.2 Steady-State SolutionsWe proceed by solving for the steady-state solutions, as Paczyn?ski did. Weset the time derivatives in equations (3.21) and (3.22) to zero, and obtainthe following system of non-linear equations [28]a? = HP?ecg2 , (3.23)a?+ fb =?e? (1? ?) . (3.24)These equations are functions of density and temperature, once the chemicalcomposition, accretion rate, surface gravity and heat flux from the core havebeen specified. We have two unknowns, ? and T , and two equations, and wecan therefore solve this system numerically. Note that the right-hand sideof equation (3.24) is the dimensionless flux from the surface F/FEdd, whichcan be deduced from equation (3.20). The flux from the surface cannotexceed the Eddington flux in our model and therefore a? + fb < 1 so thatF/FEdd < 1.In our analysis, we will vary the accretion rate, a?, while keeping theother parameters constant, namely, X,Y, g and fb. In fact, we will have adistinct set of solutions for each set of parameters, also known as a linearseries of stellar models.343.2. Steady-State Solutions3.2.1 ParametersWe specify the following set of parameters, [28]g = 108.5 g cm?2, (3.25)M = 0.929 M, (3.26)R = 0.009 R, (3.27)X = 0.7, (3.28)Y = 0.27, (3.29)XCN = Z/3. (3.30)Following Paczyn?ski?s analysis, fb will range from 0 to 0.3 and a? will varyover a wide range of values. The chemical composition is that of a solarcomposition star [16]. The Eddington flux and critical accretion rate overthe whole star for this set of parameters is 2.79 ? 1019 ergs cm?2 s?1 and3.26? 1019 g s?1, respectively.To reiterate, we use the following expressions for the opacity, energygeneration rate of the CN cycle and the equation of state [28]?rad(T, ?,X, Y ) = ?e(1 + 2? 1026(0.001 + Z)?T?3.5),CN(T, ?,X, Y ) = 8? 1027?XXCNT?2/36 exp(?152.31T?1/36),P (T, ?,X, Y ) = kBNA?T? +a3T4 +K?5/3,P (T, ?,X, Y ) = g??,given by equations (2.32), (2.9), (2.45) and (3.9), where ?e is given by equa-tion (2.29) and K is given by equation (2.44).3.2.2 Numerical SchemeThe numerical scheme we employ to solve the system of non-linear equationsutilises both Newton?s method and the bisection method [3]. The numericalscheme, which uses a nested loop structure, is outlined below.Loop 11. Choose ??, which fixes P through equation (3.9)353.2. Steady-State SolutionsLoop 22. Choose an initial guess for T3. Choose an initial guess for ? by approximating the pressure asthe sum of the ideal gas pressure and radiation pressure, ? =?TkBNA(P ? Prad)4. Use a one-dimensional Newton?s method to solve for ? from thecomplete equation of state, equation (2.45)5. Once ? is found, calculate ?rad,  and a? from equations (2.32),(2.9) and (3.23)6. Substitute all necessary quantities into equation (3.24) to definea residual, R = ?e? (1? ?)? a?? fb7. If R 6= 0, choose a new value of T , Tnew, and repeat the processfrom step 38. If R(T ) < 0 and R(Tnew) > 0, or vice versa, a zero-crossing hasbeen found9. Use the bisection method to find the value of T for which R = 0,identifying a steady-state solutionEnd Loop 21. Return to step 1, choose a new value of ?? so as to iterate throughall relevant values of ??End Loop 13.2.3 ResultsWe present our results for the steady-state solutions and compare them toPaczyn?ski?s results. Figure (3.1) depicts the change in column density withdimensionless accretion rate. Note that log denotes the logarithm to baseten, while ln refers to the natural logarithm.The accretion rates, a?, range from 10?9 to 1, which translates into M? =5.1608?10?16 M yr?1 to 5.1608?10?7 M yr?1. This is the accretion rateover the whole star, such that ?? = M?/(4piR2). These curves represent aone-parameter family of curves, since each value of fb produces a distinctsolution. For fb = 0.0001, 0.001, 0.01 and 0.1, the curves have two turning363.3. Linear Stability Analysispoints indicated by a relative maximum and minimum within a sufficientlydefined neighbourhood. Between these turning points, the column densitydecreases with increasing a?. The peaks or relative maxima tend to occur at alower value of a? as fb decreases, while the troughs, relative minima, occur atapproximately the same value of a?. The curve for fb = 0.3 is monotonicallyincreasing while the curve for fb = 0 has a trough but no peak. We compareour results to Paczyn?ski?s [28] in figure (3.2).Figures (3.2a) and (3.2b) are plotted against the same axes and for thesame values of the heat flux from the core, fb. The results appear to bequalitatively similar, except for the curve for fb = 0.01. The curve forfb = 0.01 in Paczyn?ski?s figure is offset in log a? from the other curves. Inour results, figure (3.2a), the curve for fb = 0.01 lies in between the curvesfor fb = 0.001 and fb = 0.1. The general trend seems to be that the localmaxima occur at a smaller value of a? as fb decreases. There is no indicationthat the solution for fb = 0.01 should deviate from this.To investigate this hypothesis further, we plotted two additional curvesfor fb = 0.03 and fb = 0.003. The curve for fb = 0.03 should sit below thatfor fb = 0.01 and above fb = 0.1. The curve for fb = 0.003 should sit belowthat for fb = 0.001 and above fb = 0.01. The results are depicted in figure(3.3).As we suspected, the curve for fb = 0.1 lies in between the curves forfb = 0.03 and fb = 0.003. We conclude that the curve for fb = 0.01 shouldnot be offset to the right in the log a?-axis as it is in Paczyn?ski?s figure. Weassume that this was the result of a plotting error.To facilitate a better understanding of the behaviour of the linear se-ries, we conduct a linear stability analysis, which was also conducted byPaczyn?ski [28].3.3 Linear Stability AnalysisIn this section, we perform a linear stability analysis that allows us to in-vestigate the stability of the linear series. The perturbations in density andtemperature, ?T and ??, can be expressed as [6]T = T0 + ?T ?? T0(1 + ?TT0), (3.31)? = ?0 + ?? ?? ?0(1 + ???0). (3.32)373.3. Linear Stability AnalysisIn the linear theory,??? ?TT0???  1 and??? ???0???  1. Higher powers of these per-turbations and products of perturbations are assumed to be negligible in alinear stability analysis. We assume that the perturbations have an expo-nential time dependence associated with an eigenvalue ??TT0= T? (?) exp(?t), (3.33)??T0= ??(?) exp(?t), (3.34)where ?, T? and ?? can, in general, be complex. The real part of the eigenval-ues are indicative of the stability of the steady-state solutions. In particular,a real part that is greater than zero indicates that the perturbation growsexponentially with time. Thus, we have an instability. If the real part ofthe eigenvalue is less than zero, then the perturbation will be quenched anddecay with time, indicating a stable solution.Note that ?TT0 and??T0can be written as ? lnT and ? ln ?, respectively. Wewill make use of this equivalent notation in our analysis. The perturbationsin the other variables, ? ln s, ? ln?, ? ln  and ? lnP can be expressed in termsof ? lnT and ? ln ?.As we discussed in chapter 2, the opacity and energy generation rates canbe written in power-law form. We can extend this to the pressure as well,and we have the following power-law forms of the opacity, energy generationrate and pressure [16]? = ?0???T ?T , (3.35) = 0???T ?T , (3.36)P = P0???T?T , (3.37)where ?0, 0 and P0 are constants. The exponents can be calculated asfollows?? =? ln?? ln ? , ?T =? ln?? lnT , (3.38)?? =? ln ? ln ?, ?T =? ln ? lnT , (3.39)?? =? lnP? ln ? , ?T =? lnP? lnT . (3.40)We will write the results of the linear stability analysis in terms of theseexponents. In differential form, the power-laws ared ln? = ??d ln ?+ ?Td lnT, (3.41)383.3. Linear Stability Analysisd ln  = ??d ln ?+ ?Td lnT, (3.42)d lnP = ??d ln ?+ ?Td lnT. (3.43)These relationships provide the infinitesimal changes in the respective vari-ables, and can therefore be used to express the perturbations ? ln?, ? ln  and? lnP . In accordance with Paczyn?ski?s approach we define a dimensionlessentropy S such that [28]ds = P?T dS, ST =?S? lnT , S? =?S? ln ?. (3.44)Using this, we can write dS asdS = STd lnT + S?d ln ?. (3.45)We also define [28]? = fba? , (3.46)and the characteristic thermal and nuclear timescales [28]?th =P?, (3.47)?n =E?HXH. (3.48)We return to the time-dependent system, equations (3.21) and (3.22),and linearise to allow for calculation of the eigenvalues, ?. To illustratethe linearisation process, we perturb the equations isobarically, dPdt = 0 and? lnP = 0, and calculate the thermal eigenvalue, ?th. This is the usualprocedure for conducting a thermal stability analysis.We perturb equation (3.21) isobarically and linearise to find, in terms ofthe dimensionless variable dS,P?d?Sdt = ? ln ?cg2?P (1? ?) (?? ln k + 4? lnT ) . (3.49)Note that in this notation, any variable that is not preceded by ? is anunperturbed quantity satisfying the steady-state equations (3.24) and (3.23).These quantities should be subscripted with a 0. However, we omit thesesubscripts to simplify the notation. Simplifying, using equations (3.24) and(3.23), we findd?Sdt =?P(? ln ? (a?+ fb)a? (4? lnT ? ? ln k)), (3.50)393.3. Linear Stability AnalysisWe substitute for the differentials, ?th and ? to findd?Sdt = ??1th (??d ln ?+ ?Td lnT ? (1 + ?) ((4? ?T )? lnT ? ??? ln ?)) .(3.51)After simplification and using the fact that d ln ? = ??T /??d lnT fromd lnP = 0, we findd?Sdt = ??1th ??1? (???T ? ???T ? (1 + ?) (4?? ? ?T?? + ?T??)) ? lnT.(3.52)The left-hand side of this equation, using the differential for S and lnP , is??1? (ST?? ? ?TS?)d ln ?Tdt . (3.53)Combining the last two equations, we find a differential equation in ? lnTd? lnTdt = ??1tha1 ? (1 + ?)a2a3? lnT, (3.54)wherea1 = ?T?? ? ?T ??, (3.55)a2 = 4?? ? ?T?? + ?T??, (3.56)a3 = ST?? ? ?TS?. (3.57)The solution of this differential equation yields ? lnT = T? exp(?tht), whichis precisely the form of the perturbation, where?th = ??1tha1 ? (1 + ?)a2a3, (3.58)and T? does not depend on t. The eigenvalue can also be deduced by substi-tuting for ? lnT = T? exp(?tht) in the differential equation (3.54) and solvingfor ?th. This result agrees with that of Paczyn?ski?s.To find the nuclear eigenvalue, we perturb equations (3.21) and (3.22)adiabatically, so that dsdt = 0, dSdt = 0 and ?S = 0. By perturbing equation(3.21) adiabatically, we find0 = ?+ ? lnP ? (1 + ?)(?? ln?? ? lnP + 4? lnT ). (3.59)403.3. Linear Stability AnalysisFrom this, we can find an expression for ? ln ? as a function of ? lnT ,? ln ? = ?(?T + ?T ? (1 + ?)(4? ?T ? ?T ))(?? + ?? + (1 + ?)(?? + ??))? lnT. (3.60)Perturbing equation (3.22), we findE?HXd? lnPdt = ?(? ln + ? lnP ). (3.61)Rewriting ? lnP as ??d ln ? + ?Td lnT using (3.43) and ? ln  using (3.42),and using the expression for ? ln ? in terms of ? lnT , we eventually find thefollowing differential equation in ? lnT for the nuclear eigenvalued? lnTdt =HE?HX(1 + ?)(a1 + a2 + a4)a1 ? (1 + ?)a2? lnT, (3.62)wherea4 = ?T?? ? ?T ?? + 4??, (3.63)We also define a5 asa5 = ST ?? ? ?TS?. (3.64)Solving for the eigenvalue by substituting for ? lnT = T? exp(?nt), wefind?n = ??1n(1 + ?)(a1 + a2 + a4)a1 ? (1 + ?)a2, (3.65)where ?n is defined in equation (3.48). This result also agrees with the resultfound by Paczyn?ski.Now, we use both equations (3.21) and (3.22) to derive a single differen-tial equation for the eigenvalues, ?1 and ?2. Perturbing equation (3.21), wefindd?Sdt = ??1th (? ln + ? lnP ? (1 + ?)(?? ln?? ? lnP + 4? lnT )) . (3.66)Substituting for ? ln , ? lnP and ? ln? using equations (3.42), (3.43) and(3.41), we find an expression for ? ln ? in terms of ? lnT and d?Sdt? ln ? = ?thd?Sdt?? + ?? + (1 + ?)(?? + ??)+(1 + ?)(4? ?T ? ?T )? (?T + ?T ))?? + ?? + (1 + ?)(?? + ??? lnT. (3.67)413.3. Linear Stability AnalysisNow, we can find an expression for d? lnPdt = ddt(?T ? lnT +??? ln ?) in termsof ? lnT , by substituting for ? ln ? from the expression above. Combiningthis with equation (3.66), we find thatd?Sdt =??1th?? + ?? + (1 + ?)(?? + ??)(?n(a1 ? (1 + ?)a2)? lnT??th?n??d2?Sdt2 ? (1 + ?)(a1 + a2 + a4)? lnT+(1 + ?)(?? + ??)?thd?Sdt). (3.68)The left-hand side of equation (3.68) upon rearranging becomes?th (?? + ?? + (1 + ?)(?? + ??))d?Sdt , (3.69)in which the term (1 + ?)(?? + ??)d?Sdt cancels with the term on the right-hand side of equation (3.68). Using equations (3.43), (3.42) and (3.45), andwith a little manipulation we find??dS = S?? lnP + a3? lnT, (3.70)and??dS = S?? ln + a5? lnT. (3.71)Thus,??d?Sdt = S?d? ln dt + a5d? lnTdt , (3.72)??d?Sdt = S?d? lnPdt + a3d? lnTdt , (3.73)??d2?Sdt2 = ?S???1n(d? ln dt +d? lnPdt)+ a3d2? lnTdt2 . (3.74)Substituting these expressions into equation (3.68), we see that the expres-sion S?(d? ln dt + d? lnPdt)on the left and right-hand sides cancel. Eventually,we are left with?n?tha3d2? lnTdt2 + (?th(a3 + a5) + ?n(a2(1 + ?)? a1))d?lnTdt + . . .+(1 + ?)(a1 + a2 + a4)? lnT = 0. (3.75)423.4. ResultsSubstituting for ? lnT = T? exp(?t), we obtain a quadratic equation for ?which is identical to Paczyn?ski?s?n?tha3?2 +(?th(a3 + a5) + ?n(a2(1 + ?)? a1))?+(1+?)(a1 +a2 +a4) = 0.(3.76)We solve for the eigenvalues using the quadratic formula,?1,2 = ?b???b?2 ? 4a?c?2a? , (3.77)wherea? = ?n?tha3, (3.78)b? = (?th(a3 + a5) + ?n(a2(1 + ?)? a1)) , (3.79)c? = (1 + ?)(a1 + a2 + a4). (3.80)The results we obtained for the thermal eigenvalue, nuclear eigenvalue and?1 and ?2 corroborate Paczyn?ski?s findings.3.4 ResultsWe present the results of the eigenvalue calculations for a dimensionless heatflux from the core of fb = 0.01. This is the curve that has been shifted tothe right in the log a? axis in Paczyn?ski?s figure, figure (3.2b). The resultsconform to the general trend of the eigenvalue results presented by Paczyn?skiin 1983 [28] for fb = 0.001. This gives further credence to our result thatthe curve for fb = 0.01 should not be shifted.Figure (3.4) depicts the dimensionless eigenvalues, ?th?1 and ?th?2, ver-sus the logarithm of the dimensionless accretion rate a?. The eigenvalues arepurely real for most of the domain except near the turning points, locatedat log a? ? ?2.1 and log a? ? ?0.24. In the vicinity of the turning points,the eigenvalues become complex conjugates of each other [28]. This is thereason that only the red line indicating ?2 is visible near the turning points.The quantity ?1 has the same real part over this portion of the domain. Wealso note that the real part of the eigenvalues cross from Re(?1, ?2) < 0 toRe(?1, ?2) > 0 at log a? ? ?2.1, which corresponds to the local maximum ofthe linear series for fb = 0.01. The real part of the eigenvalues cross fromRe(?1, ?2) > 0 to Re(?1, ?2) < 0 at log a? ? ?0.24, which corresponds tothe local minimum of the linear series for fb = 0.01. From this observation,we can conclude that the local maximum indicates the onset of thermalinstability, while the local minimum indicates the onset of thermal stability.433.4. ResultsThe thermal and nuclear eigenvalues, equations (3.58) and (3.65), arevalid in the limit that ?n/?th  1 and |?n/?th|  1 [28]. These limitscorresponds to regimes away from the turning points, as can be seen infigure (3.5). In the vicinity of the local extrema, the thermal and nucleareigenvalues become comparable in magnitude, as illustrated in figure (3.5),and ?1 and ?2 are required to adequately describe the behaviour of thesteady-state solutions.We can calculate the gradient of the linear series d log ??d log a? from equations(3.24) and (3.23). We find thatd lnPd ln a? =d ln ??d ln a? =(1 + ?)a2 ? a1(1 + ?)(a1 + a2 + a4). (3.81)The numerator, (1+?)a2?a1, corresponds to the negative of the denomina-tor of the nuclear eigenvalue, excluding ?n. The numerator also correspondsto the negative of the numerator of the thermal eigenvalue. At the turn-ing points, the gradient of the linear series passes through zero and we candeduce that (1 + ?)a2 ? a1 = 0 at the turning points. Thus, the thermaleigenvalue will also be zero, while the nuclear eigenvalue will diverge at theturning point. Figure (3.6) depicts the gradient of the linear series versusthe dimensionless accretion rate. The behaviour of ?th and ?n is depictedin figure (3.5).Our result for the gradient of the linear series, equation (3.81), agreeswith Paczyn?ski?s, up to a negative sign. The sign in our result is correct,which can be seen from figure (3.6). The gradient is positive before wereach the local maximum and is negative afterwards, indicating that thefunction is increasing before the turning point and decreasing afterwardswhich is consistent with a local maximum. A similar argument for the localminimum indicates that the sign of the gradient in our result is correct.Figure (3.7) depicts the variation of the eigenvalues, ?th(?1, ?2, ?th, ?n),with the accretion rate, a?, in the vicinity of the local minimum. The realparts of ?1 and ?2 pass through zero at the turning point, from a positivereal part to a negative real part. This indicates a transition to stability. Thethermal eigenvalue also passes through zero, but at a slightly different pointto that of ?1 and ?2. This is consequence of the fact that the limit used tocalculate the thermal eigenvalue is not valid in the vicinity of the turningpoint, although it is a good approximation. The nuclear eigenvalue divergesnear the turning point, as expected.We have shown that the steady-state solution is stable when the columndensity increases with the accretion rate and unstable between the turningpoints when the column density decreases with increasing accretion rate.443.5. DiscussionWe can analyze this behaviour in terms of the relative magnitudes of fband a?. The dimensionless accretion rate actually indicates the amount offlux produced by nuclear burning. Specifically, a? = ??E?HX/FEdd, which isthe dimensionless flux produced through the nuclear burning of hydrogen.When the heat flux from the core, fb, is much larger than the flux producedfrom nuclear burning, fb  a?, the model is stable. In this case, very lit-tle heat is produced from nuclear burning, and the accretion rate is verylow. This corresponds to the branch of the linear series to the left of thelocal maximum. The branch between the turning points in the linear se-ries correspond to regimes where a? and fb are similar in magnitude. Theenergy generation rate is very sensitive to temperature. For a temperatureof 2 ? 107 K, the temperature in the CN cycle is raised to an exponentof 18.04. Heat is lost in the system through radiation pressure, which isproportional to T 4. Thus, any small increase in temperature results in netheating causing the temperature to rise again, and so forth. In other words,we have a thermonuclear runaway event and instability. The density startsto decrease at the local maximum, as can be seen from figure (3.9). Theenergy generation rate is also proportional to ?. Eventually, the densitydecreases enough to restore stability to the system despite the increasingtemperature. Thus, in regions to the right of the local minimum, for whichfb  a?, stability is restored. The temperature at the bottom of the layerincreases with increasing accretion rate, which is depicted in figure (3.8).3.5 DiscussionWe have investigated the nuclear burning of hydrogen on accreting whitedwarfs using Paczyn?ski?s one-zone model [28]. Through a number of sim-plifications implemented by Paczyn?ski, the most significant being the ap-proximation of the integration of the model equations over the layer, wehave arrived at a set of two differential equations, (3.21) and (3.22). Weset the time derivatives in these equations to zero to allow investigation ofthe steady-state behaviour of the model. This resulted in two non-linearequations in T and ?, namely, equations (3.24) and (3.23).To solve this system of non-linear equations, we employed an iterativescheme utilising both Newton?s method and the bisection method. We solvedthis system for a range of values of the dimensionless heat flux from the coreover a range of accretion rates for X = 0.7, Y = 0.27 and g = 108.5 cm s?2.Each value of the dimensionless heat flux from the core, fb, produces adistinct solution. Thus, we have a linear series of stellar models for each453.5. Discussionfb. We found that curves with fb = 0.1, 0.01, 0.001, 0.0001 have two turningpoints, in agreement with Paczyn?ski. The local maxima of these curves tendto occur at a lower value of a? as fb decreases. The curve for fb = 0.3 ismonotonically increasing, and represents a limiting value, above which noturning points exist [28].Our results agree with those of Paczyn?ski?s except for the model withfb = 0.01. In figure (3.2), we compare our results to that of Paczyn?ski?s.The curve for fb = 0.01 is shifted towards increasing log a? in the log a?-axisin Paczyn?ski?s figure 1 which appeared in [28]. The curve for fb = 0.01 liesbetween the curve for fb = 0.001 and fb = 0.1 in our results. In figure(3.3), we plot two additional curves, fb = 0.03 and fb = 0.003, along withfb = 0.01. The curve for fb = 0.01 lies between the curves for fb = 0.03 andfb = 0.003. This follows the general trend of the curves, and we concludethat there was a plotting error in Paczyn?ski?s graph.To better understand the behaviour of these curves, we conducted a lin-ear stability analysis, which was also conducted by Paczyn?ski. We assumedan exponential time dependence in the perturbations associated with aneigenvalue, ?. We solved for the eigenvalues, ?1 and ?2, of the system bylinearizing equations (3.21) and (3.22). These eigenvalues are plotted in fig-ure (3.4) for fb = 0.01, which agrees with the behaviour of the eigenvaluesfor fb = 0.001, plotted by Paczyn?ski in [28]. This corroborates our findingthat the model for fb = 0.01 should follow the trend of the other steady-statemodels. We see that near the local extrema the eigenvalues become com-plex conjugates of each other. The real parts of the eigenvalues pass throughzero from a negative value to a positive value near the local maximum. Thisindicates that the peak signals the onset of instability where the perturba-tion grows exponentially with time. The real parts of the eigenvalues passthrough zero again near the local minimum, but from a positive real part toa negative real part. This indicates the onset of stability, characterised byan exponential decay of the perturbation.We also solved for a thermal and nuclear eigenvalue by perturbing equa-tions (3.21) and (3.22) isobarically and adiabatically respectively. This ap-proximation to the full calculation of the eigenvalues is not valid near theturning points where the thermal and nuclear eigenvalues are similar in mag-nitude. The approximation holds when |?n/?th|  1 [28]. We also foundthe gradient of the linear series and noted that when the gradient passesthrough zero, corresponding to a turning point, that the thermal eigenvaluealso passes through zero while the nuclear eigenvalue diverges.In figure (3.7), we plotted ?1, ?2, ?th and ?n near the local minimum ofthe stellar model fb = 0.01. It is clear from this figure that the real parts of463.5. Discussion?1 and ?2 become negative as they pass through zero. The point at which?1, ?2 = 0 corresponds to the exact point of the local minima, as expected.The thermal eigenvalue also passes through zero, but at a slightly differentvalue of log a? than that of ?1 and ?2. The nuclear eigenvalue diverges nearthis point. The point of divergence of the nuclear eigenvalue, and the pointat which the thermal eigenvalue passes through zero is slightly differentto where ?1, ?2 = 0. This is a consequence of the approximation used tocalculate the thermal and nuclear eigenvalues, which is not valid near theturning points.We observe that when the column density increases with increasing ac-cretion rate, the models for fb < 0.3 are stable. Between the turning points,these models are unstable. The column density always increases with in-creasing accretion rate for the curve fb = 0.3. We can also conclude from ourstability analysis that this model is always stable. In this respect, fb = 0.3represents a limiting value of fb, above which stellar models are stable.We could explore various other aspects of this model, including the effectsof changing the chemical composition and surface gravity. However, theseaspects were discussed in Paczyn?ski?s paper [28]. We have used his modelas a base with which to begin our investigation. We now move on to a moresophisticated model that more accurately describes the aspects of nuclearburning of interest to us.473.5. Discussion?8 ?6 ?4 ?2 08.899.29.49.69.8log a?log4?( gcm?2)  f b =0.30.10.010.0010.00010Figure 3.1: Column density as a function of log a?, where a? is the dimen-sionless accretion rate. The steady-state solutions are shown for a range ofdimensionless heat flux from the core, fb. The hydrogen and helium massfractions at the surface are X = 0.7 and Y = 0.27, respectively. The surfacegravity is g = 108.5 cm s?2.483.5. Discussion?6 ?4 ?2 08.899.29.49.69.8log a?log4?( gcm?2)  f b =0.30.10.010.0010.00010(a) Steady-state solutions for a range of fb and a?with X = 0.7, Y = 0.27 and g = 108.5 cm s?2. Thisfigure is figure (3.1) replotted so as to compare withPaczyn?ski?s figure [28].(b) Figure 1, Paczyn?ski, B. (1983) A One-Zone Modelfor Shell Flashes on Accreting Compact Stars. TheAstrophysical Journal, 264:282-295. Reproduced bypermission of the AAS.Figure 3.2: Comparison of results493.5. Discussion?6 ?5 ?4 ?3 ?2 ?1 08.88.999.19.29.39.49.59.6log a?log??( gcm?2)  f b = 0.03f b = 0.003f b = 0.01Figure 3.3: Steady-state solutions for fb = 0.03, 0.01, 0.003 with X =0.7, Y = 0.27 and g = 108.5 cm s?2.503.5. Discussion?2.2 ?2 ?1.8 ?1.6 ?1.4 ?1.2 ?1 ?0.8 ?0.6 ?0.4 ?0.2 0?0.4?0.200.20.40.60.811.21.4log a?? th? 1,2  ?1?2Figure 3.4: Variation of eigenvalues, ?1 and ?2, computed using the time-dependent system of equations, equations (3.21) and (3.22), with dimen-sionless accretion rate for fb = 0.01, with X = 0.7, Y = 0.27 andg = 108.5 cm s?2. In the vicinity of the turning points, the eigenvalues be-come complex conjugates of each other. Thus, only the red line indicating?2 is visible near the turning points, as ?1 has the same real part over thisportion of the domain. The real part of both eigenvalues crosses from a neg-ative to a positive quantity at the local maximum and vice versa at the localminimum, indicating a transition to instability and stability, respectively.513.5. Discussion?2.5 ?2 ?1.5 ?1 ?0.5?1.5?1?0.500.511.5log a?? th? th,n  ?th?ndata3data4Figure 3.5: Variation of the thermal and nuclear eigenvalues with di-mensionless accretion rate for fb = 0.01, with X = 0.7, Y = 0.27 andg = 108.5 cm s?2. These eigenvalues are valid in the limit ?n/?th  1, where|?n/?th|  1 [28], corresponding to regions away from the turning points.?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0?0.4?0.3?0.2?0.100.10.20.30.40.5log a?dlog??dloga?Figure 3.6: Gradient of the linear series as a function of log a? for fb = 0.01,with X = 0.7, Y = 0.27 and g = 108.5 cm s?2.523.5. Discussion?0.45 ?0.4 ?0.35 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1?0.3?0.2?0.100.10.20.3log a?? th? 1,2,th,n  ?1?2?th?nFigure 3.7: Variation of the eigenvalues, including the thermal and nucleareigenvalues, with log a? for fb = 0.01, with X = 0.7, Y = 0.27 and g =108.5 cm s?2. This is a close-up of the local minimum for the curve fb =0.01. The real parts of ?1 and ?2 pass through zero at the local minimum,indicating a transition to stability. The thermal eigenvalue passes throughzero at a slightly different point to the local minimum while the nucleareigenvalue diverges. This is a consequence of the fact that the thermal andnuclear eigenvalues are not valid in the vicinity of the turning points, and?1 and ?2 are required to give an accurate description of the behavior at theturning points.533.5. Discussion?5 ?4 ?3 ?2 ?1 07.57.67.77.87.988.1log a?logT(K)Figure 3.8: Variation of temperature with dimensionless accretion rate forfb = 0.01, with X = 0.7, Y = 0.27 and g = 108.5 cm s?2.?5 ?4 ?3 ?2 ?1 0010203040506070log a?log?(gcm?3)Figure 3.9: Variation of density with dimensionless accretion rate for fb =0.01, with X = 0.7, Y = 0.27 and g = 108.5 cm s?2.54Chapter 4A Multi-Zone ApproachIn the previous chapter, we made use of Paczyn?ski?s one-zone model. Inthis chapter, we present a more sophisticated multi-zone approach. Onceagain, we assume that hydrogen-rich material accretes onto the surface of aspherical CO white dwarf. As with Paczyn?ski?s model, we write our equa-tions in terms of the column density, ?, or mass per unit area. The columndensity is measured from the top of the layer in this model, accounting forthe difference in signs in the set of model equations. We assume that thewhite dwarf accretes material from a companion star at a rate ??, which isthe mass accreted per unit area per second. The accretion rate over thewhole star isM? = 4piR2??, (4.1)where R is the radius of the star. We neglect the gravitational redshiftas it is small for a white dwarf. Once again, we assume that the accretedlayer is geometrically thin in comparison to the radius of the white dwarf,and thus we make use of plane-parallel geometry with a fixed gravitationalacceleration, g. We exclude any magnetic effects.4.1 Model EquationsOnce again, we use the standard stellar structure equations along with equa-tions describing the evolution of the hydrogen and helium mass fractions toobtain the following system of differential equations:?P?? = g, (4.2)?T?? =3?F16?T 3 , (4.3)?F?? = ?Tdsdt ? (H + He) , (4.4)dXdt = ?HE?H, (4.5)554.2. Steady-State EquationsdYdt =HE?H? HeE?He, (4.6)where ddt = ??t + ?? ??? is the Lagrangian derivative [15], which follows aparticular element of accreted material. As before, P is pressure, T is tem-perature, ? is the density, F is flux, X is the hydrogen mass fraction, Y isthe helium mass fraction, g is the surface gravity, ? is the opacity, s is thespecific entropy, H and He are the energy generation rates for hydrogenand helium, respectively. E?H = 6.4 ? 1018 and E?He = 5.8 ? 1017 is theenergy released by hydrogen and helium burning, respectively, in ergs pergram. The quantity ? is the Stefan-Boltzmann constant and a = 4?/c isthe radiation density constant.Equation (4.2) is the equation of hydrostatic equilibrium excluding anylocal acceleration effects. Equation (4.3), which appears here in a slightlydifferent form to that used in Paczyn?ski?s work, is the radiative transferequation [16]. We can derive this form of the equation from that used inPaczyn?ski?s model. Recall that Prad = a3T 4, and therefore?Prad?? = 4aT33?T?? .Equating this expression to ?Fc and rearranging yields equation (4.3). Wecan include the effects of conduction in this equation by replacing ? by ?tot,which contains both conductive and radiative opacity sources. Equation(4.4) is the heat balance or conservation of energy equation. In particular,(H + He) is the energy generated per gram per second by nuclear burningand grav = ?T dsdt > 0 is the energy generated through compression. Equa-tions (4.5) and (4.6) describe the evolution of the hydrogen and helium massfractions due to nuclear burning.4.2 Steady-State EquationsTo determine the steady-state solutions, we must set the time derivatives ??tto zero, which results in the following equations for the flux, hydrogen andhelium mass fractions:dFd? = ?T ??dsd? ? (H + He) , (4.7)dXd? = ?1??HE?H, (4.8)dYd? =1??(HE?H? HeE?He). (4.9)The equation of state is a function of temperature, density and the chem-564.2. Steady-State Equationsical composition, P = P (T, ?,X, Y ). Using the chain rule, we find thatdPd? =?P?TdTd? +?P??d?d? +?P?XdXd? +?P?YdYd? . (4.10)We can express dsd? in a similar manner.After some rearrangement, we find the following set of differential equa-tions for the density, ?, temperature, T , flux, F , hydrogen mass fraction, Xand helium mass fraction, Y :d?d? =(g ? ?P?TdTd? +?P?XdXd? +?P?YdYd?)(?P??)?1, (4.11)dTd? =3?F16?T 3 , (4.12)dFd? = ?T ??( ?s?TdTd? +?s??d?d? +?s?XdXd? +?s?YdYd?)? (H + He) ,(4.13)dXd? = ?1??HE?H, (4.14)dYd? =1??(HE?H? HeE?He). (4.15)4.2.1 Surface Boundary ConditionsTo proceed with the integration of the set of steady-state equations, wemust specify five surface boundary conditions. We identify the surface ofthe accreted layer with the photosphere, located at an optical depth of 23[16]. The photosphere is the visible surface of the layer, where the majorityof the radiation escapes.We choose a surface flux, Fs. The precise value of Fs is determinedthrough the solution of the steady-state equations, a process which will bedescribed shortly. We note that the flux at the surface in the steady-state isthe sum of the flux from the bottom emanating from the white dwarf core,Fb, the compressional flux, Fc, expressed through the term grav = ?T dsd?in equation (4.13), and the flux from nuclear burning. The total flux fromnuclear burning, Fnuc, isFnuc = ?? (E?HXs + E?He(Xs + Ys)) . (4.16)In Paczyn?ski?s model, the compressional flux was neglected, as illustratedby equation (3.24). In general, the compressional flux is orders of magnitude574.2. Steady-State Equationssmaller than the flux produced by nuclear burning [35]. However, if the layerbecomes thick, the compressional flux will heat the core [39].Once Fs is chosen, we can determine an effective temperature throughthe Stefan-Boltzmann law,Fs + Facc = ?T 4s , (4.17)where Facc is the flux released by material accreting onto the surface ofthe white dwarf. In other words, this is the gravitational potential energyreleased per unit area per second when material accretes onto the surface.Specifically,Facc = ??c2z(1 + z), (4.18)where z is the redshift factor. We can solve for the temperature at thesurface such thatTs = ((Fs + Facc)/?)1/4 . (4.19)To calculate the surface density, ?s, we use the definition that the surfaceis located at an optical depth, ? ,? = 23 = ?rad?s, (4.20)where ?rad is the radiative opacity, given by equation (2.32) at the surface,and ?s is the column density at a depth where ? = 23 . We define ?0 such that?0 = ?e(1 + 2? 1026(0.001 + Z)), so that ?rad = ?0?T?3.5. We integratethe hydrostatic equilibrium equation down to a depth ?s, so that Ps = g?s.At the surface, degeneracy effects are negligible, and the pressure can beapproximated by the ideal gas law. Solving for ?s, in terms of the pressureyields?s =kBNA?sTs?g . (4.21)We substitute this expression for ? into equation (4.20) to find? = ?radkBNA?sTs?g =23 . (4.22)Substituting for ?rad, we find, after some simplification, that?s =(2?gT 2.5s3?0kBNA)1/2. (4.23)We can solve for ?s, which is where we start our integration, using equation(4.21).584.2. Steady-State EquationsFinally, we specify the composition of the accreted material, the surfacegravity and mass of the white dwarf. We assume the donor star is a solartype star, so that Xs = 0.7, Ys = 0.27 and Zs = 1?Xs?Ys = 0.03. We alsoset g = 108.5 cm s?2 and M = 0.929M. Thus, we use the same chemicalcomposition, surface gravity and mass as we did in the previous chapter,and which was used by Paczyn?ski [28].The surface boundary conditions could of course be specified differently.For example, we could use zero boundary conditions, specifying P = 0,T = 0 at the surface, ? = 0 as in [20]. If were to specify just P = 0 at ? = 0,we could then integrate the radiative transfer equation, equation (4.13),to solve for the surface density, as in [26]. In fact, there are many otherways in which the surface boundary conditions can be specified. However,the steady-state solutions are generally insensitive to the surface boundaryconditions, and converge to the same solutions for a specified Fs.4.2.2 Inner Boundary ConditionThe inner boundary condition is defined by the temperature of the whitedwarf core, such that the temperature at a specified depth in the layermust be equal to the core temperature. In Paczyn?ski?s analysis, the flux atthe bottom of the layer was fixed. We evaluated the effects of using thisboundary condition within this model. However, fixing the temperatureis a more natural and physically intuitive boundary condition to enforce,and we proceed in this manner. We set the white dwarf core temperatureto Tcore = 4 ? 107 K, a choice that was substantiated in the introductorychapter.The boundary condition can be satisfied through two methodologies.The first method is confined to the accreted layer, while the second methodrequires integration into the stellar substrate, the white dwarf core. Theresults presented in this chapter use only the method confined to the layer.However, in the interest of completeness, we present the details of bothmethods.Accreted LayerWe integrate the steady-state equations using the five surface boundary con-ditions, Fs, Ts, ?s, Xs, Ys, to some depth, ?bottom. At this depth, we requirethe temperature at the bottom of the layer to be equal to the temperatureof the white dwarf core, T (?bottom) = Tcore. When the inner boundarycondition is satisfied, we have found a steady-state solution. We allow the594.2. Steady-State Equationshydrogen and helium mass fractions to evolve as necessary to match thetemperature of the white dwarf core at the bottom of the layer. This con-trasts with Paczyn?ski?s analysis, in which the hydrogen mass fraction wasset to zero at the bottom of the layer.Degenerate electrons in the white dwarf core efficiently conduct heatthus rendering the core essentially isothermal. As such, we assume that thesubstrate beneath the layer is isothermal. This assumption is also consistentwith the observation that the evolutionary timescale of the white dwarf coreis much longer than the evolutionary timescale of the layer.We confine our analysis to the accreted layer when matching the innerboundary condition in this manner. This methodology was used by Narayanand Heyl in [25]. In practice, there are several numerical techniques that canbe used to impose this boundary condition, which we will discuss shortly.Note that the column density depth is approximately ?bottom, since ?? =?bottom ? ?s ? ?bottom.Stellar SubstrateOnce again, we integrate the steady-state equations using the five surfaceboundary conditions, Fs, Ts, ?s, Xs, Ys, to some depth, ?bottom. Thereafter,we integrate into the stellar substrate from this depth. We add a differentialequation for the thermal diffusion timescale, ?diff , to the set of five coupleddifferential equations [26],??diff?? =9?kBNA (?? ?bottom)128?T 3? , (4.24)with ?diff = 0 at ?bottom. The form of the equation above is attained byregarding equation (4.13) for the flux as a diffusion problem [26].We introduce the accretion timescale, defined as the time required toproduce a layer of column density ?bottom, at an accretion rate of ??, in thesteady state,tacc =?bottom??. (4.25)We integrate the set of six differential equations from a depth ?bottom untilthe diffusion timescale is twice the accretion timescale at a depth ?diff , [26]?diff = 2tacc. (4.26)At ?diff , we require the temperature to be equal to the core temperature,T (?diff) = Tcore. If this condition is not satisfied, we choose a different604.3. Numerical Schemevalue for the flux at the surface, Fs, and repeat the integration to a depth?diff . Once the inner boundary condition is satisfied, we have converged ona steady-state solution. The substrate is assumed to be isothermal below?diff .A CO white dwarf core contains no hydrogen or helium, and conductionby degenerate electrons is the dominant mechanism of energy transport. Amethod that integrates into the core must take these physical aspects intoaccount. Thus, this methodology is more involved than the method confinedto the accreted layer.4.3 Numerical SchemeWe have a set of five coupled ordinary differential equations to integratedown to a specified depth. To do this, we make use of an adaptive Runge-Kutta-Fehlberg method [3], [8]. This method adapts the step size at eachstep so that the local truncation error stays within a specified bound.We use Runge-Kutta-2 and Runge-Kutta-3, which are second and thirdorder explicit ordinary differential equation solvers, to perform the step sizeadjustment. In particular, we use the difference in the accuracy between thetwo methods, 4, at every step to attain an estimate of the local truncationerror. If the error is unacceptable, in terms of a specified tolerance, then thestep size is recalculated to ensure that the error, 4, is within the specifiederror bound. The procedure is repeated until we find a step size that yieldsan error within the specified error bounds, or the step size falls below theminimum specified step size. If the latter situation occurs, then the routineterminates. We use Runge-Kutta-3 to advance the solution at each step.The Runge-Kutta-Fehlberg [8] routine outlined here was adapted from analgorithm provided by Professor Anthony Peirce [31]. There are higher-orderRunge-Kutta-Fehlberg methods that use, for example, Runge-Kutta-4 andRunge-Kutta-5 to estimate the local truncation error [3], [8]. However, thephysics in our model is accurate to at most one-decimal place and the lowerorder methods are sufficient for the accuracy that we require.Generally, a numerical ordinary differential equation solver will only bestable if the step size is below a certain value. To calculate the upper limit onthe step size, one would have to conduct a stability analysis of the numericalmethod for the system of equations at hand. Since our routine adjusts thestep size to remain within a specified tolerance, the step size is generallybelow the upper bound for stability.The adaptive method outlined here is more efficient than a standard614.3. Numerical SchemeRunge-Kutta-3 method for this problem. A Runge-Kutta-3 method woulduse a uniform step size which could potentially be very small in order toavoid instabilities in the numerical routine. The adaptive method decreasesthe step size when the gradients in the solution change rapidly and increasesthe step size if the specified tolerance can be attained with a larger step size.Our numerical routine rejects any non-physical solutions, such as neg-ative or imaginary parts in the temperature, density, hydrogen or heliummass fractions. If the routine encounters such solutions, it goes back a stepand recalculates the solution with a smaller step size. This methodology issufficient to avoid instabilities in the numerical method, and non-physicalsolutions.The solution of these types of differential equations, equations (4.11)to (4.15), following the evolution of an accreted layer, whether for a whitedwarf or neutron star, are notoriously stiff. The temperature, density, fluxand mass fractions can change on very different length scales. Some ofthe variables may change rapidly with ? while others change slowly with?. From the solutions, this is evidently true of the mass fractions, whichchange slowly over most of the domain, while the density and temperaturechange more rapidly. Numerical methods require increasingly smaller stepsizes to meet integration tolerances when dealing with solutions that changeon different length scales. Due to limits on computational time, however, theimplementation of a minimum step size is necessary. As such, the routineterminates, as we have mentioned, once the solver tries to reduce the stepsize beyond the minimum step size to meet the specified tolerance. Oursolutions are, therefore, limited in certain regimes where the stiffness ofthe problem precludes further investigation. Other numerical approachesto the problem, including those that are computationally intensive can beemployed to follow the evolution of the layer in detail. A treatment of theproblem in this manner is outside the scope of this thesis. As previouslymentioned, the FLASH and MESA codes can be used to investigate burningon the surface of white dwarfs in detail [10], [30].4.3.1 Imposing the Inner Boundary ConditionThere are several numerical methodologies we could employ to apply theinner boundary condition in the accreted layer. We outline two of thesemethods. In our analysis, we use values of the flux at the surface, Fs, inunits of the Eddington flux, FEdd, such that fs = Fs/FEdd.The first method iterates through values of fs for each value of ?, specif-ically ?bottom, that we select as the maximum layer depth in steady-state.624.3. Numerical SchemeMethod 1Loop 11. Set up a grid in column density, ?bottom, and iterate through a rangeof column densities in fixed incrementsLoop 22. Set up a grid in fs and iterate through fs in fixed increments3. Integrate the differential equations for each value of fs, and de-termine if the temperature at the bottom of the layer, specifiedby the value of ?bottom in loop 1, matches the core temperature4. If the temperature condition is satisfied, save the solution, oth-erwise move on to the next value of fsEnd Loop 25. Return to step 1, choose a new value of ?bottom so as to iterate throughall relevant values of ?bottom.End Loop 1This routine integrates the equations for each fs multiple times, performingone iteration for each value of ?bottom on the grid. We could also implement abisection method from step 2, in which the method steps in fs for a particular?bottom, until a zero crossing in the temperature is found. A bisectionmethod could then find the solution fs that matches the temperature at adepth ?bottom. However, we would have to step along the whole grid in fsto avoid missing multiple solutions that might occur at the same value of?bottom.The second method uses a search routine to find steady-state solutionsafter integrating the equations for a particular value of fs to a specifiedmaximum depth.Method 21. Choose a maximum column density, for example ?max = 1010 g cm?2Loop 12. Set up a grid in fs, and iterate through each value of fs634.3. Numerical Scheme3. For each value of fs, integrate the equations down to ?max4. Use a search routine to search the temperature solution, T (?) where? ? [?s,?max], for values that match the inner boundary condition5. If the search routine finds a matching temperature value, save thevalue of ? at which it occurs, thereby identifying ?bottomEnd Loop 1Thus, we integrate the equations for each fs once, and determine if and atwhat depth the inner boundary condition is satisfied. After conducting nu-merical experiments, we determined that this is the most efficient method asit identifies all solutions for each fs without having to integrate the equationsto multiple depths.To utilise this method, Method 2, we must specify a temperature toler-ance that matches the temperature to a specified number of decimal places.Depending on the tolerance, there may be multiple solutions at a particularvalue of fs which lie very close to each other in column density. These so-lutions are approximations to a single solution point, and occur as a resultof the tolerance, the maximum step size in the Runge-Kutta routine, andthe manner in which the inner boundary condition is imposed through thesearch routine. We could increase the tolerance and decrease the maximumstep size in the Runge-Kutta method to avoid these additional solutions.This would increase the computational time. To circumvent this, we imple-ment a simple routine that excludes solutions that are too close to each otherin column density at the same value of fs. In other words, we essentiallytake the first instance of a solution in each branch at each fs. There maybe regions, however, in which there are still multiple solutions. These donot detract from the behaviour of the solutions, if they are indeed approx-imations to a single data point. This is reasonable as the input physics isat most correct to one decimal place. Furthermore, the logarithmic columndensities we will consider will only be significant to one decimal place. In-deed, this will be evident in the graphs of log ?? versus log fs. To reiterate,the general behaviour of the solutions is unaltered with this method, andplotting at least the first instance of the solution in each branch is sufficientfor our purposes of determining the unstable and stable regimes and othershell flash properties. A comparison of the steady-state results producedby this method and a more accurate bisection method is presented for aneutron star in the appendix.644.4. Steady-State SolutionsAs with all numerical methods, a more efficient method is usually lessaccurate than a method that is less efficient. There will always be a compro-mise between speed and accuracy. We have found that the method describedabove is sufficient for our purposes, satisfying both our accuracy and speedrequirements.4.4 Steady-State SolutionsThe steady-state model and method of solution we have outlined thus fardiffers from Paczyn?ski?s model in several respects. In particular, we do notapproximate the Lagrangian derivatives or the subsequent integration overthe column density of the accreted layer, ??. Furthermore, we include he-lium burning through the triple-? reaction, equation (2.18), without screen-ing effects. We neglect screening effects as the densities generally do notreach high enough values, at the temperatures encountered in the layer, forscreening effects to be significant [14]. We present our results for the steady-state solutions in the log fs? log ?? plane, where fs is the dimensionless fluxfrom the surface, for a fixed accretion rate ??, in contrast to the one-zonemodel where the steady-state results were presented in the log a? ? log ??plane for a fixed flux from the core, fb. This presentation enables us tocompare our results to similar models for compact stars. In other words,we would like to compare the solutions for different accretion rates. In ourplots, ?? = ?bottom ? ?s ? ?bottom, such that log ?? ? log ?bottom.In all our analysis, we use equation (2.7) for the total mean molecularweight, equation (2.29) for the electron scattering opacity, equation (2.32)for the radiative opacity, equation (2.45) for the equation of state, equation(2.57) for the entropy, equation (2.9) for the energy generation rate of theCNO cycle, and equation (2.18) for the energy generation rate for heliumburning.In order to draw comparisons with Paczyn?ski?s work, we will begin ouranalysis using the same parameters and expressions for the non-relativisticdegeneracy pressure, using equation (2.43), opacity and energy generationrate for hydrogen burning as we did in the previous chapter for the one-zonemodel. We will refer to this model as the basic model.We will also conduct a more sophisticated analysis based on a more gen-eral approach, which we will refer to as the extended model. In particular, weinclude energy generation through electron capture and the proton-protonchains, conduction and a degeneracy pressure encompassing both the non-relativistic and relativistic realms. Specifically, we use equation (2.33) for the654.4. Steady-State Solutionsconductive opacity and equation (2.35) for the total opacity. For the totalenergy generation rate for hydrogen burning, we use the sum of the energygeneration rates of the proton-proton chains, the CNO cycle and electroncapture, equations (2.8), (2.9) and (2.24) respectively, H = pp +CNO +ec.While we do not expect the electron capture energy generation rate to beof significance in this problem, as the density of the layer rarely exceeds 107g cm?3, we include it for completeness. We use an equation of state whichis the sum of the ideal gas, radiation and degeneracy pressures, equation(2.45), in which the last term for the degeneracy pressure is replaced bydegeneracy pressure calculated through equation (2.39). Entropy is givenby equation (2.57).Another aspect to consider is the fraction of CNO elements, XCNO, tobe used in CNO, equation (2.9). In Paczyn?ski?s model, we used XCN = Z/3,which remained constant, since CNO element production was excluded inthe one-zone model. In this model, helium undergoes burning which resultsin the production of carbon and oxygen. We assume that all the heavyelements produced are CNO elements such that XCNO = nZs + (Z ? Zs),where Zs is the fraction of heavy elements at the surface and n is the fractionof Zs assumed to be CNO elements. We will investigate both n = 13 , as inPaczyn?ski?s model, and n = 0.8 as used by Narayan and Heyl [26].To reiterate, we use g = 108.5 cm s?2, M = 0.929 M, R = 0.009 R,Xs = 0.7 and Ys = 0.27 in our analysis.4.4.1 XCNO = Zs/3We set the accretion rate over the whole star to M? = 6.3 ? 10?9 M yr?1,a value which was determined through numerical experiments to illustratekey aspects of the solutions. This translates into ?? = M?/(4piR2) = 0.0813 gcm?2 s?1. We consider a flux from the surface fs = Fs/FEdd in units of theEddington flux, FEdd = cg/ke, ranging from 10?3 to 2a? = 2??E?HX/FEdd.In flux units, Fs ranges from 10?3FEdd to 2??E?HXs. The Eddington flux isgiven by equation (3.16) which has a value of 2.79 ? 1019 ergs cm?2 s?1 forthe surface gravity we have chosen.We investigate the solutions in several regions, comparing the flux pro-duced from the nuclear burning of hydrogenFH = ??E?HXs = a?FEdd, (4.27)with that at the surface. Specifically, we investigate fs  a?, fs < a?, fs = a?and fs > a?.664.4. Steady-State SolutionsFrom the total flux produced by nuclear burning, equation (4.16), wecan define a second parameter, a?nuc, such thata?nuc =????nuc, (4.28)where ??nuc is a critical accretion rate defined as??nuc =FEddE?HXs + E?He(Xs + Ys). (4.29)Basic Model, ? ? 1010 g cm?2We conduct our analysis for the same set of parameters as Paczyn?ski, withthe same equations of state, entropy, opacity and energy generation rate forhydrogen burning, with the addition of helium burning through the triple-?reaction. We search for column densities ranging from 108.7 to 1010 g cm?2,which fits within the range of column densities found in Paczyn?ski?s model.For the value of ?? chosen, the dimensionless accretion rate, a?, is 0.0131 anda?nuc = 0.0147.We search a grid of dimensionless heat flux from the surface in the rangeindicated above, 10?3 to 2a?, in increments of 3? 10?4. We use a tolerancethat matches the core temperature to two decimal places. Since the grid sizein fs is quite small and temperature tolerance quite low, there will be caseswhere we have two solutions at the same value of fs which are very close toeach other. We exclude solutions at the same value of fs if the differencein the logarithmic column densities is less than 0.015, a value determinedthrough numerical experiments. The behaviour of the steady-state solutionsis depicted in figure (4.1).Firstly, we note that the solution has a completely different behaviourto that found in the previous chapter. Besides the fact that we are solving asystem of differential equations instead of a system of non-linear equations,the shape of the solution is attributable to the manner in which we haveapplied the boundary conditions, specifically the inner boundary condition.We can identify three distinct branches. There is a bottom branch, wherethe column density decreases with increasing fs, a top branch in whichthe column density increases with increasing fs and a near vertical branchsignifying a change in the slope of the solutions.The lower branch has the property that the column density decreases asthe flux from the surface increases. Furthermore, the solutions on the lowerbranch are those for which nuclear burning has essentially not begun. For674.4. Steady-State Solutionscolumn densities less than 109 g cm?2, the mass fractions have not changedto two decimal places. At approximately 109.49 g cm?2, approximately 0.046of the hydrogen has burned. It is only above 109.5 g cm?2 that a significantfraction of the hydrogen starts to burn. Given that the accuracy of ouranalysis is one decimal place, it is reasonable to conclude that no hydrogenburning, to any noticeable degree, occurs in this branch. Thus, the fluxemanating from the surface is essentially that emanating from the whitedwarf core in this branch. The layer is not particularly thick in this branch,indicating that the compressional flux is negligible. The degenerate coreis an excellent conductor of heat, while the outer non-degenerate envelope,the accreted layer in this case, is opaque to the radiation exiting the core.It therefore follows that the lower the column density, the less opaque thelayer is to the exiting radiation, the higher the flux from the surface andcorresponding effective surface temperature.The top branch has the property that the column density increases as thesurface flux increases. Furthermore, the solutions are those for which nuclearburning has begun to a significant degree, indicating that the temperaturein the layer reaches or exceeds 107.6 K, the approximate temperature atwhich hydrogen burning begins. The higher the surface flux, the larger theamount of hydrogen burning. Consequently, solutions on the top branch arecharacterised by partial hydrogen burning to complete hydrogen burning inthe layer. However, temperatures in the layer do not become high enough toignite helium burning, which occurs at approximately 108.1 K. Furthermore,the solutions on the top branch are accompanied by a flux reversal, indicat-ing that the temperature reaches a peak exceeding the core temperature atsome depth ?, and decreases thereafter to match the core temperature atthe bottom of the layer. This is consistent with the fact that temperaturesexceed 107.6 K in the layer, igniting hydrogen burning, and must eventuallydecrease to match the boundary condition. The reversal in the flux directionresults in the heating of the white dwarf core.The deeper the layer, the higher the peak temperature reached in thelayer. A higher peak temperature indicates that more hydrogen is burned,thereby producing a larger flux from nuclear burning. As a result of thisincreased nuclear burning with increasing column density, we find that theflux emanating from the surface is larger for larger column densities. Thetop branch of solutions exist only for fs < a?. Beyond this, at fs ? a? = 0.013,only the lower branch of solutions exist. It is interesting to note that f = a?is the total flux that can be produced from hydrogen burning.The vertical branch of solutions occur for one value of the flux, fs =0.0016 < a?. There are no solutions for fs < 0.0016 or fs  a?. The temper-684.4. Steady-State Solutionsature in the layer does not get high enough to match the inner boundarycondition for fs  a?. For example, for fs = 10?3, the peak temperaturein the layer is approximately 3.82 ? 107 K. We can explain the existenceand non-existence of certain solutions, and elucidate the behaviour of thesolutions explained above by considering X,Y, ?, T and f = F/FEdd as afunction of ? for certain values of fs.The solutions for fs = 0.0016 = 0.12a? are depicted in figure (4.2). Thered line indicates the core temperature. We see that the temperature in thelayer matches the core temperature for at least two values of ?, which arevery close to each other. These solutions indicate the turning point, andare likely approximations to a single data point, if we were to increase thetolerance. These are the points at which the hydrogen mass fraction dropsto 0.6 to one decimal place.The solutions for fs = 0.0064 = 0.49a? are depicted in figure (4.3). Thetemperature increases up to a maximum and then decreases sharply, as wasthe case for fs = 0.0016. We have two solutions, one where hydrogen burninghas not begun, X = 0.69 ? 0.7, and another where it has begun, X = 0.1,corresponding to a steady-state solution on the lower branch and top branch,respectively. We note that as fs increases, the gap in column densitiesbetween the two solutions on the lower and top branch increase. This isa consequence of the fact that the flux produced from hydrogen burningincreases with fs, and thus the peak temperature in the layer increases.A larger column density is required for the temperature to reach a higherpeak temperature and to decrease from a higher peak temperature to matchTcore, thus explaining the larger gap between solutions at the same fs as fsincreases.At fs = a?, the temperature at the end of the domain flattens out, butdoes not decrease. Thus we have returned to a region in which there isonly one solution for the range of column densities we have chosen. Asfs increases beyond a?, the temperature increases with increasing log ??,obeying a power-law-like relationship. The solutions at fs = 0.013 andfs = 0.026, at a? and 2a?, are depicted in figures (4.4) and (4.5), respectively.We can conclude from figures (4.2) to (4.5) that for 0.0016 < fs < a?there are two solutions, one that appears on the lower branch and one thatappears on the upper branch. All solutions for which 0.0016 < fs < a? areaccompanied by a temperature inversion for solutions on the upper branch.For fs < 0.0016, there are no solutions and for fs ? a?, there is only onesolution appearing on the lower branch.We summarise the behaviour of the solutions in figure (4.6), in whichwe indicate how many solutions exist for values of fs satisfying fs < 1.2a?,694.4. Steady-State Solutions1.2a? < fs < a?, and fs ? a?. The three solution regions are indicated by blackdashed lines. In regions where there are two solutions, one solution occurswhere hydrogen burning has essentially not begun, and another where partor all of the hydrogen has been burned. To reiterate, no helium burningoccurs at the column densities considered.As a further confirmation of the steady-state results, we investigate thechange in temperature at a specified depth for a range of fs. The results aredepicted in figure (4.7). Figures (4.7a) and (4.7b) depict the temperaturevariation with fs at a depth of 109.2 g cm?2 and 109.7 g cm?2 respectively.Figure (4.7a) indicates only one solution at a depth of 109.2 g cm?2, confirm-ing the steady-state findings. The temperature increases with increasing fsat this depth. This verifies our explanation that on the lower branch, whereminimal or no hydrogen burning occurs, the flux emanating from the sur-face is that emanating from the white dwarf core. Figure (4.7b) depictsthe variation at a depth 109.7 g cm?2, a column density that lies on theupper branch of the steady-state solutions in which hydrogen burning hasbegun. The temperature changes more interestingly at this depth, decreas-ing until fs ? 0.005, and then increasing thereafter. Once again, there isonly one solution matching the core temperature, which corroborates thesteady-state solution found at fs ? 0.0085 for ? = 109.7 g cm?2. At thisdepth and for the range of fs chosen, there will be a temperature inversionin the layer until fs = a?. Thus, the inverted portion of the temperaturesolution increases with increasing fs until fs ? 0.005, and thereafter theinverted portion decreases until at fs = a? it has disappeared. Beyond a? notemperature inversion occurs.Helium burning can be ignited in this model. As an example, we setfs = 0.02 = 1.5a?, and we integrate down to a depth of ? = 1013 g cm?2.The nuclear fuel consisting of all the hydrogen and helium is exhausted inthis system at the specified depth. The depth at which all the hydrogen isburned corresponds to the point where the flux first plateaus. The point atwhich the second plateau occurs corresponds to the depth at which all thehelium is completely burned. This two step plateau behaviour is distinctiveof the exhaustion or quenching of a fuel source. As all the nuclear fuelhas been burned, Z = 1, or rather XCNO = 1, at a depth just above 1012g cm?2. The almost vertical drops in the flux correspond to the amountof flux produced through nuclear burning. In particular, the first verticalsection corresponds to hydrogen burning while the second vertical sectioncorresponds to helium burning. At fs = 0.02, there is only one steady-statesolution that lies on the lower branch.704.4. Steady-State SolutionsExtended Model, ? ? 1010 g cm?2As a comparison, we conduct our analysis using a more general model inwhich we include electron capture and the proton-proton chains in the energygeneration rate for hydrogen burning. We also include conductive opacityand an electron degeneracy pressure which is valid in both the relativisticand non-relativistic regions. We find that the steady-state solutions, de-picted in figure (4.9), do not change, confirming our earlier assumption thatthe densities reached at the column densities we have considered are not highenough for electron capture, conduction and relativistic degeneracy effectsto be of significance on a white dwarf.4.4.2 XCNO = 0.8ZsWe proceed to investigate steady-state solutions for which XCNO = 0.8Zswithin the more general framework described above. To reiterate, we includeelectron capture and the proton-proton chains in the energy generation ratefor hydrogen burning as well as conductive opacity and an electron degen-eracy pressure which is valid in both the relativistic and non-relativisticregions. As before, we use equation (2.57) for the entropy. Furthermore, weexpand the range of the column densities to a maximum of 1016 g cm?2.Search MethodologyWe conduct the search in stages. Firstly, we search for solutions with acolumn density less than 109.9 g cm?2 on the same grid that we used inthe previous section. Thereafter, we search for solutions with a columndensity larger than 109.9 g cm?2, between the same ranges of fs with a gridsize of 5 ? 10?5. The reasoning behind this is that the tolerance for thetemperature matching and the methodology for the exclusion of solutionswhich are too close to each other in column density is different in each ofthese realms. The adaptive Runge-Kutta method uses larger steps whenintegrating down to a larger column density. Therefore, it is possible tostep over certain solution values if the temperature matching condition istoo stringent. Thus, we reduce the temperature matching condition to onedecimal place, and double the exclusion gap for two logarithmic columndensities at the same value of fs to 0.03 for column densities larger than109.9 g cm?2. To reiterate, for the most part this is equivalent to taking thefirst instance of a solution satisfying the temperature matching condition ineach solution branch, and is sufficient for our accuracy and computationalspeed requirements.714.4. Steady-State SolutionsBehaviour of SolutionsFigure (4.10) shows the solutions in approximately the same range as thesolutions discussed in the previous section. The behaviour is very similar tothose of last section, except that we have more solution points satisfying themodel and numerical method criteria. This is evident in the vertical portionof the solution. Figure (4.11) shows the steady-state solutions above a col-umn density of 109.9 g cm?2. Figure (4.12) depicts the entire solution space,with three branches, a lower, middle and upper branch. Figure (4.11) dis-plays unexpected characteristics, where the column density jumps from ap-proximately 1010.7 g cm?2 at fs = 0.01275 to 1014.06 g cm?2 at fs = 0.01285.There appears to a value of fs for which there is only one solution on thelower branch in this vicinity at approximately fs = 0.0128, which is approx-imately 0.98a?. The solutions, X,Y, T, ? and fs as a function of ?, for thisvalue of fs are depicted in figure (4.13). The solutions for fs = 0.01285 aredepicted in figure (4.14).There seems to be a second value of fs, fs = 0.0144 ? 0.98a?nuc, forwhich there is only one solution on the lower branch. We note that if wewere to exclude the compressional flux in this model, then these values of fsfor which there is only one solution on the lower branch would be situatedat a? and a?nuc. In other words, these values of fs correspond to the total fluxproduced from hydrogen burning and the total flux produced from hydro-gen and helium burning respectively. The addition of the compressional fluxshifts these solutions to slightly smaller values in the flux. Recall that theflux at the surface is the sum of the flux from the core, the flux from nuclearburning and the compressional flux. We are unable to integrate deeper intothe layer for fs = 0.98a? or fs = 0.98a?nuc due to the stiffness of the prob-lem. It is unclear whether these values of fs will eventually lead to a secondsolution at a much higher column density than we have considered. Alter-natively, it is possible that approximations in the input physics prohibitsus from seeing these solutions. On theoretical grounds, we do not expectsolutions above the lower branch for these values of fs if the flux from othersources is non-zero. If the flux from the core and the compressional fluxwere to add to approximately zero, then we would see solutions at highercolumn densities for these values of fs. We leave the investigation of thesepossible solution points to future work.We have explained the behaviour of the lower and middle branch in theprevious section for XCNO = Zs/3. We note that increasing the surfaceCNO element mass fraction has some effect on the solutions in the lowerand middle branch. In particular, the amount of hydrogen burning in the724.4. Steady-State Solutionslower branch increases to a small extent. At a column density of approx-imately 108.8 g cm?2, X = 0.6951, and at 109.3 g cm?2, X = 0.6562. Theturning point occurs at fs = 0.0025, which is larger than fs = 0.0016 in thecase of XCNO = Zs/3. The turning point sits slightly lower in column den-sity, but remains in the range 109.4?109.5 g cm?2. Thus the bottom branchis characterised by minimal hydrogen burning, and solutions for which notemperature inversion occurs in the layer. The middle branch is charac-terised by partial to complete hydrogen burning, indicating a temperatureinversion in the layer, and a flux reversal. This is consistent with the be-haviour found for XCNO = Zs/3. Figure (4.15) depicts the evolution of thehydrogen mass fraction with column density.Third BranchWe are left to consider the behaviour of the solutions in the highest branch,consisting of solutions that decrease and then increase with increasing fs,separated by a turning point at fs ? 0.0136. This point is situated roughlyhalf-way between fs = a? and fs = a?nuc, at approximately 0.98 (a?+a?nuc)2 .The solutions in the third branch all exhibit helium burning. This it themost salient feature. Furthermore, the amount of helium burning increaseswith increasing fs. Indeed, solutions on the highest branch range from theignition of helium burning with very little helium burned to essentially com-plete helium burning, from left to right. This implies that peak temperaturesexceed 108.1 K in this branch, and thus the steady-state solutions are at-tained through a temperature inversion in the layer. This contrasts withthe solutions in the middle branch, for which there was no helium burning,only hydrogen burning. In the highest branch, no hydrogen burning occursas all the hydrogen is burned before reaching these column densities.The column density decreases with increasing fs in the branch to theleft of the turning point. The point at the which the temperature matchesthe core temperature occurs at a shallower depth as one moves towards theturning point in fs. The amount of helium burning increases, but occurs ata shallower depth as we step in increasing fs in this left top branch until wereach the turning point. These features are evident from figures (4.16) and(4.17) for fs = 0.0130 and fs = 0.0136 respectively.In Paczyn?ski?s one-zone model, an unstable region was characterised bya decreasing column density with increasing accretion rate. Narayan andHeyl conducted a similar analysis investigating the thermonuclear stabilityof hydrogen and helium burning on the surface of accreting neutron stars.The authors found that where the column density decreased with increas-734.4. Steady-State Solutionsing surface flux, the steady-state solutions were unstable. Furthermore, inregions where the column density increased with increasing surface flux, thesteady-state solutions were found to be stable. We extrapolate from theresults of Paczyn?ski and Narayan and Heyl, and postulate that in regionswhere the column density decreases with increasing fs, the steady-state so-lutions are unstable, and in regions where the column density increases withincreasing fs, the steady-state solutions are stable.Thus, while we have not conducted a stability analysis to identify stableor unstable regions, it is reasonable to assume that this branch is unstable,leading to helium flashes. In this branch, we have a decreasing columndensity with increasing fs. Recall that for temperatures exceeding 108 K,the temperature is raised to a power of approximately 40.2 in the energygeneration rate for triple-? reaction. The energy generation rate is alsoproportional to the cube of the helium mass fraction. As helium burningignites, energy is released, which in turn raises the temperature, resulting inmore helium burning and so on. In other words, we have a thermonuclearrunaway. Furthermore, as helium burning raises the temperature in thisbranch, more helium burns at a lower column density. The column densitydoes not play a significant role in raising the temperature near the bottom ofthe layer, as the triple-? reaction is responsible for this. Most of the heliumis burned in this branch, before reaching the turning point, fs = 0.0136,which is evident from figure (4.17). Thus, the energy generation rate of thetriple-? reaction decreases to some extent as the helium mass fraction hasdropped from approximately 0.8 to approximately 0.3 in the left slope of thetop branch.As a result, beyond the turning point, the column density must increaseto assist in raising the temperature in the layer sufficiently to burn theremainder of the helium. The column density therefore increases with in-creasing fs in the branch to the right of the turning point, indicating stableburning. After the turning point, helium burning increases at approximatelythe same depth until almost all the helium has been burned. The peaktemperature increases at this depth with increasing fs, as a result of the in-creasing column density. The higher the column density, the higher the peaktemperature, resulting in more helium burning. These features are evidentwhen comparing figures (4.17) and (4.18) for fs = 0.0136 and fs = 0.0139 re-spectively. For fs = 0.0141, helium is almost completely burned, indicatedby figure (4.19). Complete helium burning occurs once the flux exceedsfs = 0.0141. For fs > 0.0141, the temperature as a function of ? no longerexhibits a temperature inversion, indicating that there will be only solutionon the lower branch for fs > 0.0141.744.5. Other Accretion Rates with XCNO = 0.8ZsSummary of Steady-State Solution BranchesIn figure (4.20) we illustrate the salient features of the steady-state solu-tions. In particular, three distinct horizontal sections can be identified. Thebottom branch, where the column density decreases with increasing fs, ispresumably unstable. Minimal hydrogen burning occurs in this branch, andthe flux emanating from the surface is that from the white dwarf core. Themiddle branch exhibits stable solutions, in which significant hydrogen burn-ing is ignited and eventually completely depleted. No helium burning occursin this branch. This branch is composed of both hydrogen and helium, withthe fraction of helium depending on the amount of hydrogen that has beenburned. It is possible that solutions in this region lead to periodic hydrogenshell flashes that are not accompanied by mass loss. The third branch hasan unstable and a stable portion. Helium burning is ignited in this branch,and increases with increasing fs until almost all the helium is burned. Thisbranch is therefore a mixture of helium and CNO elements, with the fractionof CNO elements depending on the amount of helium burned. There arelikely stable and unstable helium flashes in this solution region. The amountof helium burning increases with increasing fs, which is reasonable as theflux emanating from the surface should increase as the flux produced fromhelium burning increases.In figure (4.21), we compare the variation of temperature with fs at twodepths located in the third branch. In figure (4.21a), at a depth of ? = 1012.1g cm?2, the temperature attains a minimum above 4? 107 K, indicating nosolution at this depth. This corresponds to a depth just below the turningpoint in the third branch where no solutions exist. At ? = 1012.15 g cm?2,depicted in figure (4.21b), the temperature dips below the core temperature,satisfying the inner boundary condition before eventually increasing withincreasing fs. In fact, there are two solutions at this depth, correspondingto a data points to the left and right of the turning point in the steady-statesolutions.4.5 Other Accretion Rates with XCNO = 0.8ZsWe present the steady-state solutions for a range of a?. In particular, we de-pict a? = 1.31, 0.13, 0.0041, 0.0013, 1.31?10?4 and 1.31?10?5 correspondingto accretion rates over the whole star of M? = 6.3 ? 10?7, 6.3 ? 10?8, 2 ?10?9, 6.3?10?10, 6.3?10?11 and 6.3?10?12 M yr?1 in figures (4.22), (4.23),(4.24), (4.25), (4.26) and (4.27), respectively.These accretion rates span the range of accretion rates that were dis-754.5. Other Accretion Rates with XCNO = 0.8Zscussed in the introductory chapter leading to non-hydrodynamical shellflashes as well as classical novae. Recall that for accretion rates larger thanthose leading to stable burning, M? ? 1.7 ? 10?7 M yr?1 for a solar masswhite dwarf [35], a red giant configuration may result [29]. Accretion rateshigher than the Eddington limit may preclude further accretion [17]. Inconstructing our model, we assumed that the layer is thin, and thereforethat the gravitational acceleration is fixed in the layer. Thus, we cannotinvestigate accretion rates much larger than 10?7 M yr?1, which is the Ed-dington limit. Our model is not valid if the layer is comparable in radius tothe white dwarf, or if the gravitational acceleration changes appreciably inthe layer, as in a red giant configuration. Thus, our analysis is restricted toaccretion rates M? ? 6.3? 10?7 M yr?1.Accretion rates smaller than those leading to steady burning result innon-hydrodynamical flashes, while accretion rates below 10?10?10?9 Myr?1 lead to classical novae or regimes of unstable burning [23], [20], [17]. Asdiscussed in the introductory chapter, as the accretion rate is decreased fromthe steady burning limit, the violence of the non-hydrodynamical flashesincreases until a nova occurs. Furthermore, the smaller the accretion rate,the longer the recurrence period as it takes more time at a lower accretionrate to accumulate the mass required to trigger a burst [20].4.5.1 2? 10?9 M yr?1 ? M? ? 6.3? 10?7 M yr?1We find that the steady-state solutions exhibit the same general behaviourfor accretion rates in the range 2?10?9 M yr?1 ? M? ? 6.3?10?7 M yr?1.In particular, there is a lower branch, where the column density decreaseswith increasing fs. From our analysis of the steady-state solutions of M? =6.3?10?9 M yr?1, we can conclude minimal or no hydrogen burning occursin these branches for 6.3 ? 10?9 M yr?1 ? M? ? 6.3 ? 10?7 M yr?1. Theluminosity emanating from the surface is that emanating from the whitedwarf core. This branch begins at higher column densities as we decreasethe accretion rate, thus flattening this branch before the turning point. Inaddition, the minimal hydrogen burning in the lower branches decreases withincreasing accretion rate. No hydrogen burning occurs below the turningpoint to two decimal places for M? = 6.3 ? 10?8 M yr?1. For M? = 6.3 ?10?7 M yr?1, no hydrogen burning occurs until the turning point. Thisbehaviour is depicted in figure (4.28) for M? = 6.3 ? 10?7 M yr?1. Thescenario is different for M? = 2? 10?9 M yr?1. For this accretion rate, thelower branch exhibits significant hydrogen burning, where X has dropped to0.65, above approximately 109.1 g cm?2, which is depicted in figure (4.29).764.5. Other Accretion Rates with XCNO = 0.8ZsThis possibly results in unstable hydrogen flashes. In fact, the increase inhydrogen burning in the unstable branch is indicative of a behaviour changethat occurs below accretion rates of 10?9 M yr?1, which we will discussshortly.The middle branch, in which hydrogen burning occurs, is located at ap-proximately the same column density, ?? = 109.5 g cm?2, for accretion ratesM? = 6.3 ? 10?9, 6.3 ? 10?8, 6.3 ? 10?7 M yr?1. The amount of hydrogenburning increases as fs increases. However, as we increase the accretion rate,this branch becomes convex instead of a monotonically increasing function offs. The slight convexity that appears at higher accretion rates results in thecolumn density decreasing with fs for a narrow range of fs. While a stabilityanalysis is required to test the stability of these solution points, it is reason-able to postulate that these points are unstable. The steady-state solutionsimmediately to the left and right of the local minimum are likely inaccessibleto the system. Indeed, the system will relax to the next equilibrium datapoint at the same column density, at approximately ?? = 109.5 g cm?2,which appears on the increasing slope of the hydrogen burning branch. Thesystem will then move up this stable slope.The height or column density of the hydrogen burning branch is approx-imately constant, with the branch beginning at approximately log ?? ? 9.5and ending at log ?? ? 11. The characteristic accretion timescale, calcu-lated as ?/?? through equation (4.25), provides an estimate to the recurrenceperiod of flashes. The characteristic accretion timescale to ignite hydrogenburning decreases with increasing accretion rate as tacc ? ???1. In partic-ular, the characteristic accretion timescale to ignite hydrogen burning at acolumn density of ?? ? 109.5 g cm?2 is approximately 1.2?10, 1.2?102 and1.2 ? 103 years for accretion rates 6.3 ? 10?7 M yr?1, 6.3 ? 10?8 M yr?1and 6.3 ? 10?9 M yr?1 respectively. For a 1.2 M white dwarf, Jose? etal. found that the recurrence periods for stable hydrogen flashes were 9, 50and 832 years for accretion rates 2 ? 10?7 M yr?1, 5 ? 10?8 M yr?1 and5? 10?9 M yr?1 respectively [20]. Whilst our analysis has been conductedin the steady-state, the characteristic accretion timescales are comparablein magnitude to the recurrence timescales found by Jose? et al. The hy-drogen burning branch is entirely stable for M? = 6.3 ? 10?9 M yr?1, andmostly stable for accretion rates larger than this limit. Thus, we expectnon-hydrodynamical flashes, or flashes resulting in mass increase.The top branch, in which helium burning occurs, appears at a lower col-umn density as the accretion rate is increased. Furthermore, the height ofthe branch in column density decreases with increasing accretion rate. Thelocal minimum of this branch is located at approximately 1012.1 g cm?2 for774.5. Other Accretion Rates with XCNO = 0.8ZsM? = 6.3?10?9 M yr?1 and 1011.4 g cm?2 for M? = 6.3?10?7 M yr?1. Thehelium burning branch begins at 1011.7 g cm?2 for M? = 6.3? 10?7 M yr?1,after a characteristic accretion time scale of 2.0 ? 103 years, and at ap-proximately 1014.1 g cm?2 for M? = 6.3? 10?9 M yr?1 after a characteristicaccretion timescale of 4.9?107 years. For M? = 6.3?10?8 M yr?1, the char-acteristic accretion timescale for helium burning is approximately 4.9? 104years. Thus, the characteristic accretion timescale for helium burning ismuch shorter with a larger accretion rate than with a smaller accretion rate.It is reasonable to assume that a larger accretion rate results in an earlieronset of a flash, and a smaller periodicity of flashes. For a 1.2 M, Jose? etal. found recurrence timescales of 1,195 years for M? = 10?7 M yr?1 and4,290 years for M? = 5? 10?8 M yr?1 [20].The branch of the helium burning solutions which decrease in columndensity with increasing fs are presumably unstable and lead to heliumflashes. The right slope, where the column density increases with fs, pre-sumably result in stable burning and non-hydrodynamical helium flashes.It is possible that only some of the solutions are accessible to the system, asit may relax to an equilibrium solution at the same column density.The accretion rate M? = 6.3?10?7 M yr?1 is very close to the Eddingtonlimit. In fact, for fs at the end of the hydrogen burning branch and thehelium burning branch, where fs > 0.94, the addition of the accretion fluxand the surface flux exceeds the Eddington limit, fs + facc > 1. It is,therefore, possible that accretion is actually switched off at these values offs, since the outwardly directed pressure will prohibit further deposition ofmatter onto the white dwarf surface. As a result, some of the solutions inthe hydrogen branch, and all of the helium branch may not be accessible asfurther accretion may not occur at these values of fs.Besides the fact that more hydrogen burns in the lower branch for M? =2 ? 10?9 M yr?1, this accretion rate has another interesting property. Atfs = 0.00405 ? a?, there are solutions at approximately 1015 g cm?2. Thisis difficult to discern from figure (4.24) since this value of fs is very closeto the fs in the helium burning branch. The flux at the bottom of thelayer at fs = 0.00405 is approximately zero. All the hydrogen has beenburned, but no helium burning occurs. As we mentioned previously, forthere to be a steady-state solution above the lower branch at fs = a?, whichsignifies complete hydrogen burning, the sum of the flux from the core andthe compressional flux would have to be zero. It seems that we have realisedthis situation for this value of fs at this accretion rate. To better understandthis behaviour, we would have to conduct a stability analysis, or at leastinvestigate further. However, we leave this to future work, and focus our784.5. Other Accretion Rates with XCNO = 0.8Zsattention on lower accretion rates.4.5.2 M? < 2? 10?9 M yr?1The behaviour of the steady-state solutions changes quite dramatically foraccretion rates less than M? < 2? 10?9 M yr?1. The solutions decrease incolumn density with increasing fs indicating that these are likely unstablesolutions. The amount of hydrogen burning in these branches increaseswith decreasing accretion rate. For M? = 6.3 ? 10?10 M yr?1, completehydrogen burning occurs at smaller fs where the column density starts toclimb steeply at approximately log ?? = 9.7 until approximately log ?? =10.3. However, partial hydrogen burning, X < 0.34, occurs all along thebranch with the amount of hydrogen burning increasing with decreasingfs. As the column density increases, the temperature near the base of thelayer will increase, resulting in more hydrogen burning. However the largercolumn density is more opaque to the exiting radiation resulting in a lower fs.The evolution of the hydrogen mass fraction with column density is depictedin figure (4.30). For M? = 6.3?10?11 M yr?1 and M? = 6.3?10?12 M yr?1,all the data points exhibit complete hydrogen burning. However, none ofthe solutions exhibit helium burning or a flux reversal.The column density at which these branches begin occurs at higher col-umn densities as the accretion rate is lowered. Or rather, as we decreaseaccretion rate beyond 2?10?9 M yr?1, the unstable branch starts at highercolumn density. This is reasonable, as a higher column density is requiredto ignite hydrogen burning at lower accretion rates. Furthermore, the datapoints climb to higher column densities as the accretion rate is decreased.The height of the branch, in column density, is larger with decreasing ac-cretion rate. It is possible, that at either lower accretion rates or highercolumn densities, helium burning may be ignited in these branches. Weleave to future work an investigation of steady-state solutions with heliumburning as well as a detailed exploration of the solutions, as presented forM? = 6.3? 10?9 M yr?1.These accretion rates agree with the general literature indicating un-stable hydrogen burning resulting in classical novae or nova-like outburstsfor accretion rates less than 10?9 M yr?1 [20]. The behaviour in the lowerbranch for M? = 2 ? 10?9 M yr?1 is indicative of a regime change fromsteady-state solutions with stable hydrogen branches to no stable hydrogenbranches. As we decrease the accretion rate, the characteristic accretiontimescale for hydrogen burning to occur increases as the column densityincreases. For M? = 6.3 ? 10?11 M yr?1, the branch begins at a column794.5. Other Accretion Rates with XCNO = 0.8Zsdensity of approximately 1010 g cm?2, indicating a characteristic accretiontimescale of 3.9? 105 years. For M? = 6.3? 10?12 M yr?1, the branch be-gins after a characteristic accretion timescale of 1.2? 107 years. The latterestimate is much larger than the upper limit to the recurrence timescales ofclassical novae of 105 years [37]. However, the general behaviour is consis-tent with the literature indicating that as the accretion rate is lowered therecurrence periods increases [20].Since the amount of hydrogen burning in the branch and the height ofthe branch in column density increases with decreasing accretion rate, it isreasonable to conclude that the flashes become more violent with decreasingaccretion rate. This is consistent with the literature indicating that theviolence of flashes increases with decreasing accretion rate [4].4.5.3 Accumulated MassesThe mass in the accreted layer can be calculated by multiplying the col-umn density by the surface area of the white dwarf. Figure (4.31) depictsthe mass, in solar masses, for the range of column densities encountered inour analysis. At a column density of 109.5 g cm?2, the mass in the layer isapproximately 7.8?10?6 M. This is the approximate mass that must be ac-cumulated for hydrogen shell flashes to occur for accretion rates in the range6.3 ? 10?9 M yr?1 ? M? ? 6.3 ? 10?7 M yr?1. According to Nomoto, fora solar mass white dwarf, accretion rates in the range 10?9?10?8 M yr?1lead to flashes if 3 ? 10?5 M is accreted [27]. Thus our estimate is com-parable to that of Nomoto?s. The mass of hydrogen-rich material requiredfor a flash increases with decreasing accretion rate [27]. Nomoto foundthat for an accretion rate of 10?10 M yr?1, 10?4 M must be accreted,and for 10?11 M yr?1, 3 ? 10?4 M must be accreted [27]. In our analy-sis, as we decrease the accretion rate beyond 10?9 M yr?1, the height ofthe hydrogen burning branch increases, and thus the accumulated mass in-creases. For M? = 6.3? 10?11 M yr?1, the branch begins at approximately1010 g cm?2, which amounts to an accumulated mass of 2.5? 10?5 M. ForM? = 6.3?10?12 M yr?1, the branch begins at approximately 1010.5 g cm?2,which amounts to an accumulated mass of 7.8?10?5 M. Thus, we can inferthat the violence of eruptions increases with decreasing accretion rate as aresult of the higher accumulated masses.For accretion rates 6.3?10?9 M yr?1 ? M? ? 6.3?10?7 M yr?1, heliumburning begins once a substantial amount of mass has been accumulated.In particular, the mass in the accreted layer must reach approximately 1.2?10?3 M for M? = 6.3 ? 10?7 M yr?1, and approximately 3.1 ? 10?1 M804.6. Discussionfor M? = 6.3 ? 10?9 M yr?1. In the latter case, the white dwarf mass isapproaching the Chandrasekhar limit.4.6 DiscussionIn this chapter, we extended the one-zone model presented in the first chap-ter to include both hydrogen and helium burning. In contrast to the one-zone model, we did not approximate the Lagrangian derivatives in the time-dependent system of equations, nor did we approximate the integration overthe height of the layer. The steady-state equations were supplemented withthe appropriate boundary conditions at the surface. We chose a particularvalue of fs at the surface and added the accretion flux to the value of fs todetermine the total flux at the surface. We then determined the temperatureat the surface through the Stefan-Boltzmann law, and solved for the densityand column density at an optical depth of 2/3, coinciding with the photo-sphere. We solved the steady-state equations, equations (4.11) to (4.15),using an adaptive Runge-Kutta method. We integrated the steady-stateequations down to a specified depth, ?max, and identified solutions througha search algorithm where the temperature matched that of the core. Thechoice of the core temperature, Tcore = 4? 107 K, was substantiated in theintroductory chapter. To reiterate, it is reasonable to conclude that the coretemperature of a white dwarf in a mass-transferring binary will be hotterthan that of an isolated white dwarf.The use of the adaptive Runge-Kutta method and the search algorithmwere quite efficient, certainly more efficient than a bisection method. As withall numerical methods, the more efficient the routine, the less accurate theroutine is likely to be. The physical inputs in our model are at most correctto one decimal place, and therefore the numerical methods we employedwere suitable with regards to our accuracy and speed requirements.The steady-state solutions were presented in the log fs? log ?? plane,where fs is the dimensionless flux from the surface, in units of the Ed-dington flux, for a fixed accretion rate. We conducted our analysis for?? = 0.0813 g s?1 cm?2, or M? = 6.3 ? 10?9 g cm?2. We initially conductedour investigation for the same parameter inputs, including the same CNOelement mass fraction, XCNO = Zs/3, and the same input physics as inPaczyn?ski?s one-zone model. The results, within a range of column densitiescoinciding with the solutions in Paczyn?ski?s one-zone model, do not resemblethe steady-state solutions found in the one-zone model. These steady-statesolutions are depicted in figure (4.1) and (4.6). There is a lower branch, in814.6. Discussionwhich minimal hydrogen burning occurs, a turning point, and a top branchin which hydrogen burning occurs. The lower branch is a decreasing func-tion of fs in column density, while the top branch is an increasing functionof fs.We conducted the analysis in the extended model, including electroncapture, the proton-proton chains, conduction and both relativistic and non-relativistic degeneracy pressure. The shape of the solutions did not change,and are depicted in figure (4.9). This is consistent with the understandingthat the effects of electron capture and conduction will only be significantfor densities higher than approximately 106?107 g cm?3. We then conductedthe analysis for a surface CNO element fraction of XCNO = 0.8Zs. This CNOfraction was used by Narayan and Heyl in their investigation of neutronstars accreting from a solar composition star [26]. The steady-state resultschanged slightly when the surface CNO fraction was increased. In particular,the vertical portion of the branch was extended to a small extent, with theturning point occurring at a larger value of fs. The solutions are depictedin figure (4.10). In all our calculations, the elements produced by heliumburning were considered to be CNO elements such that XCNO = 0.8Zs +(Z ? Zs).We then searched for solutions at higher column densities in the ex-tended model with XCNO = 0.8Zs. The solutions are depicted in figure(4.12). There are three distinct branches, the lower branch in which min-imal burning occurs, the middle branch in which hydrogen burning occursand the upper branch in which helium burning occurs. The lower branch isa decreasing function of fs in column density. It follows that, since minimalburning occurs, the higher the column density, the more opaque the layeris to the radiation exiting the core, and hence the smaller the surface flux.Indeed, the temperatures and densities reached at these column densitiesare not high enough to ignite significant hydrogen burning in this branch.The amount of hydrogen burning in the middle branch increases withincreasing fs until all the hydrogen has been burned. The black dashedline in figure (4.20) appears at fs ? a?, or the flux produced from totalhydrogen burning. There is only one solution on the lower branch at thisvalue of fs. There is similar black dashed line at fs ? a?nuc, or the fluxproduced from total hydrogen and helium burning. The middle branch isan increasing function of fs, indicating that the solutions are likely to bestable here. If these solutions are indeed stable, which will only be confirmedafter conducting a stability analysis, there will be periodic flashes, but nomass loss in this regime. Thus, this branch may likely lead to mass increase.No helium burning occurs in this branch, and thus we can conclude that824.6. Discussiontemperatures exceed 107.6 K in the layer, but do not reach 108.1 K, thetemperature at which helium begins burning. The temperature of the core is107.6 K, and thus all the solutions in this branch are a result of a temperatureinversion in the layer. The temperature rises to a peak, exceeding the coretemperature, and then must sharply decline near the bottom of the layer tomatch the core temperature.The amount of helium burning in the third branch increases with in-creasing fs until essentially all the helium is burned. The left slope, beforethe local minimum in this branch, is a decreasing function of fs while theright slope is an increasing function of fs. It is reasonable to assume that theleft slope is unstable whilst the right slope is stable. In the left slope, wherehelium begins to burn, Y ? 0.8, as all the hydrogen has been burned in thelower branch resulting in an accumulation of helium. The energy generationrate of the triple-? reaction is proportional to Y 3, and the temperature israised to a power of approximately 40.2. Thus, as helium begins to burn inthis branch, energy is released initiating a thermonuclear runaway. Thus,more helium can burn at lower column densities, as the rise in the tempera-ture is fueled by the flux produced by helium burning. At the turning point,the helium fraction has dropped to Y ? 0.2, the energy released throughburning has decreased and therefore, the column density must climb againto raise the temperature sufficiently to burn the remainder of the heliumin the right slope of this branch. The right slope is an increasing functionof fs and solutions are likely stable here, resulting in non-hydrodynamical,periodic, helium flashes potentially leading to mass increase. Note that notall the points in the valley on the left and right slope may be accessible tothe system.We then investigated the behaviour of the steady-state solutions for arange of accretion rates. The accretion rates were M? = 6.3 ? 10?7, 6.3 ?10?8, 2?10?9, 6.3?10?10, 6.3?10?11 and 6.3?10?12 M yr?1, which corre-spond to dimensionless accretion rates a? = 1.31, 0.13, 0.0041, 0.0013, 1.31 ?10?4 and 1.31 ? 10?5, respectively. The findings are depicted in figures(4.22), (4.23), (4.24), (4.25), (4.26) and (4.27), respectively. There aretwo distinct regimes. Accretion rates satisfying 2 ? 10?9 M yr?1 ? M? ?6.3 ? 10?7 M yr?1 exhibited similar behaviour while those for M? < 2 ?10?9 M yr?1 exhibited similar behaviour. We leave to future work thedetailed exploration of solutions for accretion rates M? = 6.3 ? 10?7, 6.3 ?10?8, 2?10?9, 6.3?10?10, 6.3?10?11 and 6.3?10?12 M yr?1, as conductedfor M? = 6.3? 10?9 M yr?1.Steady-state solutions for accretion rates in the range 2?10?9 M yr?1 <M? ? 6.3?10?7 M yr?1 contained a lower branch where minimal or no burn-834.6. Discussioning occurred, a middle branch where hydrogen burning occurred and a topbranch where helium burning occurred. For M? = 2 ? 10?9 M yr?1, theonly difference to the aforementioned behaviour is that significant hydrogenburning occurs in the lower branch. We believe that this aberration is in-dicative of the change to unstable hydrogen burning regimes that occurs forlower accretion rates.We find that the column density at which hydrogen burning begins,log ?? = 9.5, is approximately constant across the range of accretion rates2 ? 10?9 M yr?1 < M? ? 6.3 ? 10?7 M yr?1. Thus, for larger accretionrates, the characteristic accretion timescale to reach this column density isshorter, since tacc is proportional to ???1. The lower branch begins at highercolumn densities as the accretion rate is reduced, thereby flattening thelower branch of solutions. The hydrogen burning branch is approximatelyconstant in height across the range of accretion rates, with the amountof hydrogen burning increasing with fs. For M? = 6.3 ? 10?9 M yr?1,the branch increases in column density with increasing fs. However, foraccretion rates larger than 6.3?10?9 M yr?1, the hydrogen burning branchbecomes slightly convex with a small portion which is a decreasing functionof fs, before reaching a local minimum and increasing in column densitywith fs. Once the system reaches log ?? = 9.5, it will likely move towardsthe next steady-state solution with the same column density, which is onthe right stable slope of the hydrogen burning branch where most of thehydrogen has been burned. For the most part, the solutions in the hydrogenbranch are likely stable. This agrees with the literature indicating non-hydrodynamical flashing behaviour for accretion rates smaller than thoseleading to steady burning, approximately 1?4 ? 10?7 M yr?1 for a solarmass white dwarf [27], but larger than those leading to novae, approximately10?10?10?9 M yr?1 [20], [17], [23].The helium burning branch begins at a lower column density as theaccretion rate is increased, and becomes shorter in height in column densitywith increasing accretion rate. Thus, the characteristic accretion timescalerequired for helium burning to begin decreases rapidly as the accretion rateis increased, faster than ???1. This is consistent with the literature indicatinga larger time to flash, and a larger recurrence period for smaller accretionrates [20].For accretion rates smaller than 2 ? 10?9 M yr?1, the steady-state so-lutions lie on one branch, which decreases in column density with increasingfs and are thus likely unstable. Furthermore, the amount of hydrogen burn-ing increases in the branch with decreasing accretion rate, with the entirebranch exhibiting total hydrogen burning for M? = 6.3? 10?11 M yr?1 and844.6. DiscussionM? = 6.3? 10?12 M yr?1. For M? = 6.3? 10?10 M yr?1, partial hydrogenburning, X < 0.34, occurs in the lower portion of this branch, but total hy-drogen burning occurs once the solutions start to climb in column density atsmaller fs. This is reasonable as a larger column density is required to raisethe temperature in the lower layers sufficiently for total hydrogen burningto occur. Furthermore, the decrease in fs is a result of the higher opacity ata higher column density. There is neither helium burning in these branchesnor a flux reversal. The steady-state solutions for M? < 2?10?9 M yr?1 re-sult in unstable behaviour with the violence of the eruption likely increasingwith decreasing accretion rate. This agrees with the literature indicatingthat novae occur for accretion rates lower than 10?10?10?9 M yr?1 [20],[17], [23].Our model is not valid for accretion rates much larger than the Ed-dington limit, 10?7 M yr?1. For accretion rates in this vicinity or larger,envelope expansion resulting in a red giant configuration with a degener-ate core may result [29], [4]. Our model assumes a thin layer with a fixedgravitational acceleration and therefore breaks down if the gravitational ac-celeration changes appreciably in the layer. For accretion rates approachingor larger than the Eddington limit, further accretion may be inhibited asthe outwardly directed pressure inhibits the deposition of material on thesurface of the white dwarf [17].We believe that our results agree quite well with previous authors. How-ever, we have postulated the stability of solutions based on Paczyn?ski?s anal-ysis [28] as well as that of Narayan and Heyl [26]. To confirm our resultsand insights, the next step is to conduct a global linear stability analysis todefinitively determine the stability of solutions. We present this frameworkin the next chapter.854.6. Discussion0 0.005 0.01 0.015 0.02 0.0258.899.29.49.69.810log??( gcm?2)fsa? = 0 .013XCNO = Z s /3Figure 4.1: Variation of column density, ??, with dimensionless flux fromthe surface, fs, for a dimensionless accretion rate of a? = 0.0131. The surfacehydrogen, helium and CNO mass fractions are Xs = 0.7, Ys = 0.27 andXCNO = Zs/3, respectively. The surface gravity is g = 108.5 cm s?2.864.6. Discussion0 2 4 6 8 1000.20.40.60.81X,Y   X Y 0 2 4 6 8 105678logT(K)0 2 4 6 8 10?4?2024log?( gcm?3)0 2 4 6 8 10?1012 x 10?3f sfs = 0 .0016a? = 0 .013 XCNO = Z s /3log ?(g cm?2 )Figure 4.2: Solutions for fs = 0.0016 = 0.12a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are at least two steady-state solutions here.874.6. Discussion0 2 4 6 8 1000.20.40.60.81X,Y   X Y 0 2 4 6 8 105678logT(K)0 2 4 6 8 10?4?2024log?( gcm?3)0 2 4 6 8 10?4?202468 x 10?3f slog ?(g cm?2 )fs = 0 .0064a? = 0 .013 XCNO = Z s /3Figure 4.3: Solutions for fs = 0.0064 = 0.49a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are two steady-state solutions, separated by alarger column density than in figure (4.2).884.6. Discussion1 2 3 4 5 6 7 8 9 1000.20.40.60.81X,Y   X Y 1 2 3 4 5 6 7 8 9 105678logT(K)1 2 3 4 5 6 7 8 9 10?4?2024log?( gcm?3)1 2 3 4 5 6 7 8 9 10051015 x 10?3f slog?(g cm?2 )a? = 0 .013 XCNO = Z s /3fs = 0 .013Figure 4.4: Solutions for fs = 0.013 = a?, a? = 0.0131, with Xs = 0.7, Ys =0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line indicates the coretemperature. There is only one steady-state solution.894.6. Discussion1 2 3 4 5 6 7 8 9 1000.20.40.60.81X,Y   X Y 1 2 3 4 5 6 7 8 9 105678logT(K)1 2 3 4 5 6 7 8 9 10?4?2024log?( gcm?3)1 2 3 4 5 6 7 8 9 100.010.0150.020.0250.03f slog ?(g cm?2 )a? = 0 .013 XCNO = Z s /3fs = 0 .026Figure 4.5: Solutions for fs = 0.026 = 2a?, a? = 0.0131, with Xs = 0.7, Ys =0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line indicates the coretemperature. There is only one steady-state solution.904.6. Discussion0 0.005 0.01 0.015 0.02 0.0258.899.29.49.69.810log??( gcm?2)fs2 10No Heburning Partial to completeH burning,accompanied by TinversionMinimalH burningNo HeburningNumber of solutionsa? = 0 .013XCNO = Z s /3Figure 4.6: Variation of column density with fs for a? = 0.013, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line dividesthe figure horizontally into a top region and a bottom region. The red textabove the red line describes the solutions in the top region while the red textbelow the red line describes characteristics of the solutions in the bottomregion. The dashed black lines split the figure vertically indicating threeregions corresponding to no solution, two solutions and one solution fromleft to right.914.6. Discussion0 0.005 0.01 0.015 0.02 0.02533.544.555.5 x 107f sT(K)a? = 0 .013 XCNO = Zs /3log ?? = 9 .2(a) Variation of T (K) with fs at ?? = 109.2 g cm?2.0 0.005 0.01 0.015 0.02 0.025234567 x 107f sT(K)a? = 0 .013 XCNO = Zs /3log ?? = 9 .7(b) Variation of T (K) with fs at ?? = 109.7 g cm?2.Figure 4.7: Comparison of temperature variation with fs at different columndensities, with a? = 0.0131, Xs = 0.7, Ys = 0.27, XCNO = Zs/3 and g =108.5 cm s?2.924.6. Discussion2 4 6 8 10 1200.20.40.60.81X,Y   X Y 2 4 6 8 10 1256789logT(K)2 4 6 8 10 12?4?20246log?( gcm?3)2 4 6 8 10 120.0050.010.0150.020.025f slog ?(g cm?2 )a? = 0 .013 XCNO = Z s /3fs = 0 .02Figure 4.8: Solutions for fs = 0.02 = 1.5a?, a? = 0.0131, with Xs = 0.7, Ys =0.27, XCNO = Zs/3 and g = 108.5 cm s?2. The red line indicates the coretemperature. The plateaus in the flux are distinctive of an exhaustion ofnuclear fuel.934.6. Discussion0 0.005 0.01 0.015 0.02 0.0258.899.29.49.69.810log??(gcm?2 )fsa? = 0 .013XCNO = Z s /3Extended modelFigure 4.9: Variation of column density with fs for a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = Zs/3 and g = 108.5 cm s?2, in the extended model,which includes the proton-proton chains, electron capture, conduction anddegeneracy pressure in both the relativistic and non-relativistic regimes.944.6. Discussion0 0.005 0.01 0.015 0.02 0.0258.88.999.19.29.39.49.59.69.79.8log??(gcm?2 )fsa? = 0 .013XCNO = 0 .8 Z sFigure 4.10: Variation of column density with fs for a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, in the extended model.0.0115 0.012 0.0125 0.013 0.0135 0.014 0.01459.51010.51111.51212.51313.51414.5log??(gcm?2 )fsa? = 0 .013XCNO = 0 .8 Z sFigure 4.11: Variation of column density with fs for a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, in the extended model.954.6. Discussion0 0.005 0.01 0.015 0.02 0.02599.51010.51111.51212.51313.51414.5log??(gcm?2 )fsa? = 0 .013XCNO = 0 .8 Z sFigure 4.12: Variation of column density with fs for a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, in the extended model.964.6. Discussion2 4 6 8 10 12 14 1600.20.40.60.81X,Y   X Y 2 4 6 8 10 12 14 1656789logT(K)2 4 6 8 10 12 14 16?505log?( gcm?3)2 4 6 8 10 12 14 1600.0050.010.0150.02f slog ?(g cm?2 )a? = 0 .013 XCNO = 0 .8 Z sfs = 0 .01280Figure 4.13: Solutions for fs = 0.01280 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature.974.6. Discussion2 4 6 8 10 12 1400.20.40.60.81X,Y   X  Y2 4 6 8 10 12 1456789logT(K)2 4 6 8 10 12 14?505log?( gcm?3)2 4 6 8 10 12 1400.0050.010.0150.02f slog?(g cm?2 )a? = 0 .013 XCNO = 0 .8 Z sfs = 0 .01285Figure 4.14: Solutions for fs = 0.01285 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are at least two steady-state solutions.984.6. Discussion8.5 8.75 9 9.25 9.5 9.75 1000.10.20.30.40.50.60.70.8log ?? (g cm?2 )Xa? = 0 .013XCNO = 0 .8 Z sFigure 4.15: Variation of hydrogen mass fraction with column density fora? = 0.0131, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2,in the extended model.994.6. Discussion2 4 6 8 10 1200.20.40.60.81X,Y   X  Y2 4 6 8 10 1256789logT(K)2 4 6 8 10 12?505log?( gcm?3)2 4 6 8 10 1200.0050.010.015f slog ?(g cm?2 )a? = 0 .13 XCNO = 0 .8 Z sfs = 0 .0130Figure 4.16: Solutions for fs = 0.0130 ? 0.99a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are two steady-state solutions.1004.6. Discussion2 4 6 8 10 1200.20.40.60.81X,Y   X  Y2 4 6 8 10 1256789logT(K)2 4 6 8 10 12?505log?( gcm?3)2 4 6 8 10 12051015 x 10?3f slog ?(g cm?2 )a? = 0 .013 XCNO = 0 .8 Z sfs = 0 .0136Figure 4.17: Solutions for fs = 0.0136 ? 0.98a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are two steady-state solutions.1014.6. Discussion2 4 6 8 10 1200.20.40.60.81X,Y   X Y 2 4 6 8 10 1256789logT(K)2 4 6 8 10 12?505log?( gcm?3)2 4 6 8 10 1200.0050.010.015f slog ?(g cm?2 )a? = 0 .13 XCNO = 0 .8 Z sfs = 0 .0139Figure 4.18: Solutions for fs = 0.0139 ? 1.06a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are two steady-state solutions.1024.6. Discussion2 4 6 8 10 12 1400.20.40.60.81X,Y   X Y 2 4 6 8 10 12 1456789logT(K)2 4 6 8 10 12 14?505log?( gcm?3)2 4 6 8 10 12 1400.0050.010.0150.02f slog ?(g cm?2 )a? = 0 .013 XCNO = 0 .8 Z sfs = 0 .0141Figure 4.19: Solutions for fs = 0.0141 = 1.08a?, a? = 0.0131, with Xs =0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red line indicatesthe core temperature. There are at least two steady-state solutions.1034.6. Discussion0 0.005 0.01 0.015 0.02 0.02599.51010.51111.51212.51313.51414.5log??( gcm?2)fsAll HburnedNo HeburningPartial to completeHe burning asf s  increasesT inversionT inversion0 2 2 1No TinversionPartial to completeH burning asf s  increasesa? = 0 .013XCNO = 0 .8 Z sMin H burning No He burningFigure 4.20: Variation of column density with fs for a? = 0.013, withXs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2. The red linedivides the figure horizontally into three regions, a bottom, middle and topregion. The red text describes features of the solutions in the respectiveregions. The dashed black lines split the figure vertically indicating fourregions corresponding to no solution, two solutions, two solutions and onesolution from left to right, indicated by black numerals.1044.6. Discussion0.0125 0.013 0.0135 0.014 0.0145 0.0150.511.522.5 x 108f sT(K)a? = 0 .013 XCNO = 0 .8 Zslog ?? = 12 .1(a) Variation of temperature with fs at ?? =1012.1 g cm?2.0.0125 0.013 0.0135 0.014 0.0145 0.01500.511.522.5 x 108f sT(K)a? = 0 .013 XCNO = 0 .8 Zslog ?? = 12 .15(b) Variation of temperature with fs at ?? =1012.15 g cm?2.Figure 4.21: Comparison of temperature with fs at different column den-sities, with a? = 0.0131, Xs = 0.7, Ys = 0.27, XCNO = Zs/3 and g =108.5 cm s?2.1054.6. Discussion0 0.25 0.5 0.75 1 1.25 1.577.588.599.51010.51111.512log??( gcm?2)f sa? = 1 .31XCNO = 0 .8 ZsFigure 4.22: Steady-state solutions for a? = 1.31, corresponding to M? =6.3 ? 10?7 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g =108.5 cm s?2.0 0.05 0.1 0.15 0.2 0.2588.599.51010.51111.51212.5log??( gcm?2)f sa? = 0 .13XCNO = 0 .8 Z sFigure 4.23: Steady-state solutions for a? = 0.13, corresponding to M? =6.3 ? 10?8 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g =108.5 cm s?2.1064.6. Discussion2 6 10 14x 10 ?38910111213141516log??( gcm?2)f sa? = 0 .0041XCNO = 0 .8 ZsFigure 4.24: Steady-state solutions for fs, corresponding to a? = 0.0041,M? = 2 ? 10?9 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2.1 1.5 2 2.5 3 3.5 4x 10 ?39.29.49.69.81010.210.410.6log??( gcm?2)f sa? = 0 .0013XCNO = 0 .8 ZsFigure 4.25: Steady-state solutions for a? = 0.0013, corresponding to M? =6.3 ? 10?10 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g =108.5 cm s?2.1074.6. Discussion1.2 1.4 1.6 1.8 2 2.2 2.4 2.6x 10 ?41010.51111.51212.5log??( gcm?2)f sa? = 1 .31 ? 10 ?4XCNO = 0 .8 ZsFigure 4.26: Steady-state solutions for a? = 1.31 ? 10?4, corresponding toM? = 6.3 ? 10?11 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2.1 1.5 2 2.5 3 3.5 4x 10 ?510111213141516log??( gcm?2)f sa? = 1 .31 ? 10 ?5XCNO = 0 .8 ZsFigure 4.27: Steady-state solutions for a? = 1.31 ? 10?5, corresponding toM? = 6.3 ? 10?12 M yr?1, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs andg = 108.5 cm s?2.1084.6. Discussion7 7.5 8 8.5 9 9.5 1000.10.20.30.40.50.60.70.8log ?? (g cm?2 )Xa? = 1 .31XCNO = 0 .8 ZsFigure 4.28: Variation of hydrogen mass fraction with column density fora? = 1.31, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2, inthe extended model.8.8 9 9.2 9.4 9.6 9.8 10 10.2 10.400.10.20.30.40.50.60.7log ?? (g cm?2 )Xa? = 0 .0041XCNO = 0 .8 ZsFigure 4.29: Variation of hydrogen mass fraction with column density fora? = 0.0041, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2,in the extended model.1094.6. Discussion9.2 9.4 9.6 9.8 10 10.2 10.4 10.600.050.10.150.20.250.30.35log?? (g cm?2 )Xa? = 0 .0013XCNO = 0 .8 ZsFigure 4.30: Variation of hydrogen mass fraction with column density fora? = 0.0013, with Xs = 0.7, Ys = 0.27, XCNO = 0.8Zs and g = 108.5 cm s?2,in the extended model.7 8 9 10 11 12 13 14 15 16?8?7?6?5?4?3?2?1012log ?? (g cm?2 )logM ?Figure 4.31: Accumulated mass in accreted layer as a function of columndensity.110Chapter 5Stability AnalysisWe now present a framework with which to evaluate the stability of thesteady-state solutions found in the previous chapter. The analysis will notbe as simple as it was in the case of Paczyn?ski?s model. In fact, we willhave to solve for the eigenvalues numerically, as a complete analytic treat-ment is not possible due to the complexity of the problem. To reiterate, weconduct a linear stability analysis in order to find the eigenvalues associatedwith a particular steady-state solution found in the previous chapter. Thiswill allow us to determine the instability or stability of that solution point,enabling us to evaluate and potentially confirm our insights regarding thestability of the steady-state solutions. The identification of stable regimeswill allow us to identify parameter regimes that may lead to mass increaseand possible progenitors of Type Ia supernovae. The identification of unsta-ble regimes possibly lead to novae. The stability analysis is, therefore, keyto our understanding and applications of the steady-state solutions.To begin our analysis, we must return to the original, time-dependentset of equations, which we state again here,?P?? = g, (5.1)?T?? =3?F16?T 3 , (5.2)?F?? = ?Tdsdt ? (H + He) , (5.3)dXdt = ?HE?H, (5.4)dYdt =HE?H? HeE?He. (5.5)The perturbation expansion for the temperature can be written in thefollowing form,T = T0 + ?T ?? T0(?)(1 + ?TT0), (5.6)where??? ?TT0 (t, T0)??? =??T?1(?) exp(?t)??  1, and ? is the eigenvalue, which, in1115.1. Equations for Perturbed Variablesgeneral is complex. Alternatively, we can write this perturbation expansionin the formT = T0 + ?T ?? T0(?) + T1(?) exp(?t), (5.7)where T1(?) = T0T?1(?). Similarly, we can write the perturbations for theother variables as? = ?0(?) + ?1(?) exp(?t), (5.8)X = X0(?) +X1(?) exp(?t), (5.9)Y = Y0(?) + Y1(?) exp(?t), (5.10)F = F0(?) + F1(?) exp(?t). (5.11)Note that T, T1, ?, ?1, X,X1, Y, Y1, F and F1 are, in general, complex. T0, ?0,X0, Y0 and F0 are the solutions to the unperturbed problem, or the steady-state solutions. The perturbed quantities are much smaller than the steadystate solutions, |T1|  |T0|, |?1|  |?0|, |X1|  |X0|, |Y1|  |Y0| and|F1|  |F0|.We substitute these expansions into the set of time-dependent differentialequations, equations (5.1) to (5.5), to find a system of equations for theperturbed variables.5.1 Equations for Perturbed VariablesAny function of the perturbed variables, given by equations (5.7) to (5.11),can be expanded in a Taylor expansion, to first order, asf(T, ?,X, Y ) = f0 + . . .+ exp(?t)(T1( ?f?T)T0+ ?1(?f??)?0+X1( ?f?X)X0+ Y1( ?f?Y)Y0),(5.12)where f0 = f(T0, ?0, X0, Y0). For simplicity of notation, we will define ?fsuch that?f =(T1( ?f?T)T0+ ?1(?f??)?0+X1( ?f?X)X0+ Y1( ?f?Y)Y0). (5.13)We substitute the perturbed variables into the radiative transfer equa-tion, equation (5.3). Upon substitution, the left hand side of the equationbecomes?T?? =?T0?? + exp(?t)?T1?? . (5.14)1125.1. Equations for Perturbed VariablesThe right hand side, to first order in exp(?t), is3?F16?T 3 =316?T 30(?0F0 + exp(?t)(?0F1 + ??F0 ? 3T1T0?0F0)), (5.15)where ?0 = ?(T0, ?0, X0, Y0) and?? =(T1( ???T)T0+ ?1(????)?0+X0( ???X)X0+ Y0( ???Y)Y0). (5.16)Matching terms of order exp(?t), we find two equations,dT0d? =3?0F016?T 30, (5.17)?T1?? =316?T 30(?0F1 + ??F0 ? 3T1T0?0F0). (5.18)The first equation is the equation for the steady-state solution, which wehave already dealt with in the previous chapter. The equation for the per-turbed temperature, T1, to first order in exp(?t), is given by the secondequation. Combining these equations, we find that the equation for thetotal perturbed temperature, given by equation (5.7), to first order, is?T?? =3?F16?T 3 , (5.19)which takes the form of equation (5.2). This is a more convenient form towork with, in which we solve for the total perturbed quantity, instead ofjust T1.The equation for the evolution of the hydrogen mass fraction involves atime derivative. Recall that the Lagrangian derivative satisfies ddt = ??t +?? ??? . Explicitly, equation (5.4) is?X?t + ???X?? = ?HE?H. (5.20)The steady-state solution, X0, does not have any time dependence, thus?X0?t = 0 and?X0?? can really be written as a total derivative, as we did inthe steady-state case, in equation (4.14). Upon substitution of X, given byequation (5.9), we find, to first order in exp(?t)? exp(?t)X1 + ??(?X0?? + exp(?t)?X1??)= ?H,0 + ?HE?H, (5.21)1135.2. Equations for Real and Imaginary Parts of Perturbed Variableswhere H,0 = H(T0, ?0, X0, Y0) and ?H is defined by equation (5.13). Wecould split this equation into two equations, one for the steady-state solutionand one for the perturbed variable, X1. However, it is more convenient towork with X, such that, to first order?(X ?X0) + ???X?? = ?HE?H, (5.22)where we have replaced exp(?t)X1 by (X ?X0) since X ?X0 = exp(?t)X1through equation (5.9).Similarly, for the helium evolution equation to first order, we find?(Y ? Y0) + ???Y?? =HE?H? HeE?He, (5.23)where Y is given by equation (5.10).The equations for the pressure and flux take the same form as that givenin equations (4.11) and (4.13). The derivatives, however, are now functionsof the perturbed variables. In the steady-state case, the derivatives are re-stricted to evaluation at unperturbed data points. We will explicitly considerthe real and imaginary parts of ?, T, ?,X, Y, and F in the next section, andexpress the pressure and flux equations in terms of these variables instead.5.2 Equations for Real and Imaginary Parts ofPerturbed VariablesBy splitting the equations for the perturbed variables into real and imagi-nary parts, we will derive a framework that allows for the solutions of theeigenvalues. We begin by writing all quantities in terms of the their real andimaginary parts.The eigenvalues, ? are complex? = ?R + i?I. (5.24)Thus, the quantities, T, ?,X, Y and F are complex and we can rewrite theperturbed variables in terms of real and imaginary parts. The perturbedtemperature, in terms of real and imaginary parts, isT = TR + iT?I, (5.25)TR = T0 + Re(exp(?t))T1, (5.26)T?I = Im(exp(?t))T1, (5.27)1145.2. Equations for Real and Imaginary Parts of Perturbed Variableswhere TR is the real part of T and T?I is the imaginary part of T . We makean adjustment to the imaginary part, T?I, so that we can write it in the sameform as TR. We define TI asTI = T0 + T?I. (5.28)The first order Taylor expansion for a function of T isf(T ) = f(T0) + exp(?t)T1( ?f?T)0. (5.29)The real part of f(T ) can be expressed as a Taylor expansion to firstorder in the variable TR,Re(f(T )) = f(TR) = f(T0) + Re(exp(?t)T1)( ?f?T)0. (5.30)The imaginary part can be expressed similarly. However, we have to addf(T0) to both sides such thatIm(f(T )) = (TI ? T0)( ?f?T)0, (5.31)Im(f(T )) + f(T0) = f(TI), (5.32)where TI ? T0 = T?I = Im(exp(?t)T1).We have considered a function of one variable here for illustration. Thisframework can be extended to several variables such that, to first order, wehaveRe(f(T, ?,X, Y )) = f(TR, ?R, XR, YR), (5.33)Im(f(T, ?,X, Y )) + f(T0, ?0, X0, Y0) = f(TI, ?I, XI, YI), (5.34)where ?R, ?I, XR, XI, YR and YI are natural extensions of TR and TI. Westate these belowTR = T0 + Re(exp(?t)T1), (5.35)TI = T0 + Im(exp(?t)T1), (5.36)?R = ?0 + Re(exp(?t)?1), (5.37)?I = ?0 + Im(exp(?t)?1), (5.38)XR = X0 + Re(exp(?t)X1), (5.39)XI = X0 + Im(exp(?t)X1), (5.40)YR = Y0 + Re(exp(?t)Y1), (5.41)YI = Y0 + Im(exp(?t)Y1). (5.42)1155.2. Equations for Real and Imaginary Parts of Perturbed VariablesThe radiative transfer equation, equation (5.2), can be expressed in termsof real and imaginary parts as?T?? =?TR?? + i?TI?? = fR + ifI, (5.43)where fR = f(TR, ?R, XR, YR), fI = f(TI, ?I, XI, YI) and f = 3?F16?T 3 . Equat-ing real and imaginary parts, we find?TR?? =( 3?F16?T 3)R, (5.44)?TI?? =( 3?F16?T 3)I. (5.45)Equations (5.22) and (5.23) can be written as?R(XR ?X0)? ?I(XI ?X0) + ???XR?? = ?H(TR, ?R, XR, YR)E?H,(5.46)?R(XI ?X0) + ?I(XR ?X0) + ???XI?? = ?H(TI, ?I, XI, YI)E?H,(5.47)and?R(YR ? Y0)? ?I(YI ? Y0) + ???YR?? =H,RE?H? He,RE?He, (5.48)?R(YI ? Y0) + ?I(YR ? Y0) + ???YI?? =H,IE?H? He,IE?He, (5.49)respectively, where we have simplified the notation such thatH,R = H(TR, ?R, XR, YR), He,R = He(TR, ?R, XR, YR), (5.50)H,I = H(TI, ?I, XI, YI) He,I = He,I(TI, ?I, XI, YI). (5.51)The hydrostatic equilibrium equation can be written as??R?? =(g ?(?P?T)R?TR?? +( ?P?X)R?XR?? +(?P?Y)R?YR??)(?P??)?1R,(5.52)??I?? =(g ?(?P?T)I?TI?? +( ?P?X)I?XI?? +(?P?Y)I?YI??)(?P??)?1I.(5.53)1165.2. Equations for Real and Imaginary Parts of Perturbed VariablesThe energy equation can be written as?FR?? = ?T0(( ?s?X)0(?R(XR ?X0)? ?I(XI ?X0)))?T0(( ?s?Y)0(?R(YR ? Y0)? ?I(YI ? Y0)))?T0(( ?s?T)0(?R(TR ? T0)? ?I(TI ? T0)))?T0((?s??)0(?R(?R ? ?0)? ?I(?I ? ?0)))?T0??(( ?s?X)R?XR?? +( ?s?Y)R?YR?? +( ?s?T)R?TR??+(?s??)R??R??)?(TR ? T0)??(( ?s?X)0?X0?? +( ?s?Y)0?Y0?? +( ?s?T)0?T0??+(?s??)0??0??)? (H(TR, ?R, XR, YR) + He(TR, ?R, XR, YR)) .and?FI?? = ?T0(( ?s?X)0(?I(XR ?X0) + ?R(XI ?X0)))?T0(( ?s?Y)0(?I(YR ? Y0) + ?R(YI ? Y0)))?T0(( ?s?T)0(?I(TR ? T0) + ?R(TI ? T0)))?T0((?s??)0(?I(?R ? ?0) + ?R(?I ? ?0)))?T0??(( ?s?X)I?XI?? +( ?s?Y)I?YI?? +( ?s?T)I?TI??+(?s??)I??I??)?(TI ? T0)??(( ?s?X)0?X0?? +( ?s?Y)0?Y0?? +( ?s?T)0?T0??+(?s??)0??0??)? (H(TI, ?I, XI, YI) + He(TI, ?I, XI, YI)) .Equations (5.44), (5.45), (5.46), (5.47), (5.48), (5.49), (5.52), (5.53),(5.54) and (5.54) are the first-order equations for the perturbed equations1175.3. Boundary Conditionsthat we must solve, along with the steady-state equations, equations (4.11)to (4.15).5.3 Boundary ConditionsWe must solve the system of fifteen equations, five for the steady-state so-lutions, five for the real part of the perturbations and five for the imaginarypart of the perturbations. We have already specified the five boundary con-ditions for the steady-state system. Thus, we need to specify ten boundaryconditions for the real and imaginary set of equations.We begin by specifying the flux at the surface, Fs,p = Fs,0 + Fs,1, whereFs,1  Fs,0 is some small real perturbation and Fs,p is the surface flux forthe perturbed problem. We refer to the surface flux for the unperturbedproblem as Fs,0. We then solve for the temperature, Ts and density, ?s asoutlined in the previous chapter. Since, Fs,1 is real, Ts and ?s will be real.Note that the subscript s, 0 refers to the surface boundary conditions of thesteady-state problem. Through our definitions for the perturbed variables,we can define TR,s, TI,s, ?R,s, ?I,s, FR,s and FI,s asFR,s = Re(Fs,0 + Fs,1) = Fs,p, FI,s = Fs,0, (5.54)TR,s = ((Fs,p + Facc) /?)1/4 , TI,s = Ts,0, (5.55)?R,s =(2?gT 2.5R,s3?0kBNA)1/2, ?I,s = ?s,0, (5.56)Xs,R = Xs,I = Xs,0 = 0.7, (5.57)Ys,R = Ys,I = Ys,0 = 0.27, (5.58)Zs,R = Zs,I = Zs,0 = 0.03, (5.59)X + Y + Z = 1, (5.60)?R,s = ?I,s = ?s,0, (5.61)where we have made use of the same methodology, as outlined in the pre-vious chapter, to calculate Ts and ?s. Facc is the energy released from thedeposition of accreted material on the surface of the white dwarf, whichis specified by equation (4.18). At the surface, we have specified the samecomposition of the accreting material as in the steady-state case, Xs,0 = 0.7,Ys,0 = 0.27. Note that X+Y +Z = 1, and therefore the perturbed quantitiesmust satisfy X1 + Y1 + Z1 = 0.We will also have to specify the eigenvalue, and we do this by choosingeigenvalues that are purely real initially. We search for the eigenvalues1185.4. Numerical Schemethrough an iterative scheme, which will be outlined in the following section.The inner boundary condition is still the core temperature of the whitedwarf. The temperature of the perturbed quantity, T = T0+exp(?t)T1, mustequal the core temperature at the point in ? where we do the matching. Inother words, the perturbation must vanish at the matching point to satisfythe inner boundary condition.5.4 Numerical SchemeSince the steady-state solutions appear in the equations for the perturbedvariables, we integrate the total system of fifteen differential equations to-gether. In the previous chapter, we determined steady-state solutions at var-ious depths, ?bottom. At this depth, the steady-state temperature matchedthat of the white dwarf core, T (?bottom) = Tcore, and an equilibrium solutionwas found. We can use this information to solve the perturbed problem. Inparticular, we have an idea of which values of the surface flux, Fs,0 will yieldsolutions, and the depth, ?bottom, at which those solutions appear.Our aim is to find eigenvalues for which the perturbation, T1, vanishesat the matching point so that T = T0 = Tcore. To this end, we define afunction of ? asf(?) = T ? T0 = exp(?t)T1. (5.62)To match the boundary condition, f(?) = 0, implying that T1 = 0, asexp(?t) cannot be zero. Using the definition that T = TR + iT?I, we find thatT?T0 = (TR?T0)+iT?I. Recall that we defined a new variable, TI = T0 + T?I,in equation (5.28). In fact, we have used the variables, TI, ?I, XI, YI and FIin our system of equations. Substituting for T?I, we findf(?) = (TR ? T0) + i(TI ? T0) = 0, (5.63)to satisfy the inner boundary condition. Recall that the steady-state solu-tion, T0, already satisfies the inner boundary condition, and thus the pertur-bations must vanish at the bottom of the layer. We are essentially searchingfor zeros of f(?). To do this, we use Mu?ller?s method, a root-finding tech-nique solving equations of the form f(?) = 0, published in 1956 [3]. Mu?ller?smethod requires three function values to choose the next value of ?. Thus,we have to integrate the equations for three initial guesses of ?. There-after, Mu?ller?s method finds the zeros of f(?). Once a zero is found, wehave found an eigenvalue of the perturbed problem. Since we are interestedin identifying regions of stability, we will only be interested in eigenvaluesthat have a positive real part that allows for the instability to grow on a1195.4. Numerical Schemetimescale shorter than the accretion timescale. The accretion timescale istacc = ?bottom/??, as given by equation (4.25). The perturbations of interestto us will have a real part specified byRe(?) ? 3tacc, (5.64)where the factor of three ensures that the instabilities grow faster than theaccretion timescale. The factor of three was based on Narayan and Heyl,[26], in which the authors considered the growth of instabilities on accretingneutron stars. If the real part of the eigenvalue is negative, then we havestability.Equation (5.63) potentially has an infinite number of solutions, giventhat we are dealing with continuous functions, which are not polynomials.The equations for the real and imaginary parts of each of the variablesexhibit opposite signs for the coefficient of ?I. As a result of this symmetry,the eigenvalue solutions will be complex conjugate pairs. Thus, to findmultiple eigenvalues associated with a particular Fs,0, we need to factor outthe complex conjugate pairs from f(?). In particular, if ?1 is a solution,then ??1 will also be a solution, and to find another solution, we solve foreigenvalues satisfyingf(?)(?? ?1)(?? ??1)= 0. (5.65)The numerical scheme is outlined below. We would have to perform thealgorithm described below for each value of Fs,0 of interest.Loop 11. Choose three values of ? real, first iteration as initial guesses, thereafterthrough Mu?ller?s method1205.4. Numerical SchemeLoop 22. Choose the surface flux, Fs,p, where Fs,0 is a solution of thesteady-state problem3. Calculate the surface boundary conditions for the real, imagi-nary and steady-state systems4. Integrate down to a depth ?bottom5. At ?bottom, evaluate if f(?) = (TR ? T0) + i(TI ? T0) = 0, savethe solution if satisfiedEnd Loop 26. Choose a new value of ? through Mu?ller?s method, until f(?) = 0End Loop 1We briefly outline Mu?ller?s method below.Mu?ller?s MethodMu?ller?s method, an adaptation of the secant method, uses three initialguesses, (?1, f(?1)), (?2, f(?2)) and (?3, f(?3)) to determine the next ap-proximation, ?4, by considering the intersection of the ?-axis with theparabola through (?1, f(?1)), (?2, f(?2)) and (?3, f(?3)) [3]. It has theadded advantage that it will take ? into the complex plane if necessary.This is essential for the problem at hand as the eigenvalues can be complex.Mu?ller?s method uses the three initial guesses to calculate three coefficients,a, b and c [3]a = (?2 ? ?3)(f(?1)? f(?3))? (?1 ? ?3)(f(?2)? f(?3))(?1 ? ?3)(?2 ? ?3)(?1 ? ?2), (5.66)b = (?1 ? ?3)2(f(?2)? f(?3))? (?2 ? ?3)2(f(?1)? f(?3))(?1 ? ?3)(?2 ? ?3)(?1 ? ?2), (5.67)c = f(?3), (5.68)to determine the parabola P = a(???3)2 + b(???3) + c. The next approx-imation to the zero of f(?), ?4 is calculated from [3]?4 = ?3 ?2cb+ sgn(b)(b2 ? 4ac)1/2 . (5.69)1215.5. DiscussionThis formula is derived by calculating the zero of the parabola, P , by usingthe quadratic formula. The algorithm chooses the sign of the square rootto be the same as the sign of the coefficient b, so that it chooses the ap-proximation to the zero closest to ?2. The procedure is then repeated with(?2, f(?2)), (?3, f(?3)) and (?4, f(?4)) as the initial guesses to determinethe next approximation to the zero. We terminate the routine when f(?)is sufficiently close to zero or when two successive approximations of ? aresufficiently close to each other.5.5 DiscussionWe have outlined the methodology with which to test the steady-state so-lutions for stability. Stability is an important aspect to understanding thesteady-state solutions and identifying stable and unstable regimes. Whilstwe have postulated the stability of the steady-state solutions in the previouschapter, it is necessary to conduct the linear stability analysis to evaluateour insights. As previously mentioned, stable regimes could elucidate pro-genitors of Type Ia supernovae, while unstable regimes may illuminate theconditions that lead to novae. We leave the determination of the stabilityof the steady-state solutions to future work.122Chapter 6Conclusions, Future Workand Other RemarksWe have investigated the problem of nuclear burning of hydrogen-rich ma-terial leading to shell flashes on the surface of accreting white dwarfs. Shellflashes are essentially periodic increases in luminosity as a result of nuclearburning of accreted material. These events may be hydrodynamic or non-hydrodynamic. Non-hydrodynamic shell flashes, caused by stable burning,do not incur any mass loss and lead to mass increase. As a result, it ispossible that white dwarfs undergoing stable burning may reach the Chan-drasekhar limit, resulting in carbon detonation and a Type Ia supernova.On the other hand, unstable burning leading to mass loss is the likely causeof classical and recurrent novae. Regardless of the violence of the flash,these thermonuclear phenomena are significant, as they influence the chem-ical and dynamical evolution of the Universe to some extent. Our aim wasto elucidate parameter regimes, essentially the accretion rate, that lead tostable and unstable burning or flashes.As a first step, we followed Paczyn?ski?s one-zone model [28], developedin 1983, and confirmed his results by following his analysis independently.After several approximations, Paczyn?ski obtained a system of two ordi-nary differential equations describing the time evolution of the layer anda system of two non-linear equations describing the steady-state behaviour.Paczyn?ski?s model required the choice of a flux emanating from the bottomof the layer or core, fb = Fb/FEdd. Paczyn?ski used a 0.929 M white dwarfwith log g = 8.5, accreting from a solar type star with X = 0.7, Z = 0.03.Furthermore, Paczyn?ski stipulated that all the hydrogen must be burnedupon reaching the bottom of the layer. We solved the non-linear system ofequations for the steady-state solutions using a numerical routine that em-ployed both a bisection method and Newton?s method. Our results and thoseof Paczyn?ski?s are depicted in figure (3.2). Our results matched Paczyn?ski?squite well. There are three branches that can be identified for fb < 0.3.There are two branches which increase in column density with increasingaccretion rate, separated by a branch which decreases in column density123Chapter 6. Conclusions, Future Work and Other Remarkswith increasing accretion rate. Therefore, a local maximum and a localminimum separates the branches. The local maximum tends to occur at alower value of the accretion rate as fb is decreased. However, Paczyn?ski?scurve for fb = 0.01 seems to be shifted to the right in increasing log a?. Weconclude that this was a plotting error. For accretion rates larger than orequal to fb = 0.3, the curves are monotonically increasing. To understandthese solutions better, Paczyn?ski conducted a linear stability analysis usingthe system of time-dependent equations. We followed his analysis indepen-dently and produced results that match those of Paczyn?ski?s. In particular,we found that portions of the solutions which were increasing functions ofthe accretion rate resulted in stable solutions, where the real parts of theeigenvalues were negative. Branches in which the column density decreasedwith increasing accretion rate resulted in unstable solutions. The turningpoints coincided with the points at which the real parts of the eigenvaluespassed through zero. In particular, at the local maxima, the real part ofthe eigenvalues passed from a negative quantity to a positive quantity, indi-cating a transition to instability. At the local minima, the real part of theeigenvalues passed from a positive quantity to a negative quantity indicat-ing the onset of stability. Whilst the one-zone model illustrates the basicbehaviour of stable and unstable solutions, a more sophisticated model isrequired to produce realistic results.In chapter 4, we extended the one-zone model to a more sophisticatedmulti-zone model including hydrogen and helium burning. Furthermore, wedid not simplify the integration of the model equations over the layer or theLagrangian derivatives, as in the one-zone model. We included diffusion andconduction in the radiative transfer equation, as well as energy generationdue to the proton-proton chains and electron capture in addition to theCNO cycle for hydrogen burning. We used the triple-? reaction for heliumburning. We adopted the same surface gravity, surface chemical compositionand mass as that used in the one-zone model. In addition, we investigatedthe surface CNO element mass fraction forXCNO = Zs/3 andXCNO = 0.8Zs.To find the steady-state solutions, we set the time derivative in the La-grangian derivative to zero, resulting in equations (4.11) to (4.15). Wesupplemented these equations with five boundary conditions at the surfaceof the layer, coinciding with the photosphere. To solve for the depth of thelayer, we employed an inner boundary condition enforcing that the tem-perature at the bottom of the layer matched that of the white dwarf core,Tcore = 4? 107 K. To solve this system of equations, we employed an adap-tive Runge-Kutta-Fehlberg method [8], using a third order Runge-Kuttamethod to advance the solutions at each step [31]. We also made use of a124Chapter 6. Conclusions, Future Work and Other Remarkssearch algorithm to identify the bottom of the layer by matching the coretemperature to at least one decimal place. The numerical scheme was suf-ficient for our speed and accuracy requirements. In contrast to Paczyn?ski?sanalysis, the multi-zone model allows the hydrogen and helium mass frac-tions to evolve as necessary to match the core temperature. Furthermore,we produced results in the fs? log ?? plane, where fs is the dimensionlesssurface flux in units of the Eddington flux.We found the steady-state solutions for the following accretion rates,6.3 ? 10?12, 6.3 ? 10?11, 6.3 ? 10?10, 2 ? 10?9, 6.3 ? 10?9, 6.3 ? 10?8 and6.3 ? 10?7 M yr?1. The range of accretion rates chosen span those thathave been investigated by other authors, and are expected to lead to bothstable and unstable burning for white dwarfs accreting hydrogen-rich matter.Our model is not valid for accretion rates much larger than 10?7 M yr?1,which is the Eddington limit. For accretion rates equal to this limit or larger,a red giant configuration results [4], [29]. In addition, further accretion mayactually be inhibited at these accretion rates [17].We found three distinct solution branches for accretion rates in the range2? 10?9 M yr?1 < M? ? 6.3? 10?7 M yr?1. In particular, there is a lowerbranch in which minimal or no burning occurs as the temperatures anddensities reached are not high enough at column densities less than approx-imately 109.5 g cm?2 to ignite significant burning. In fact, the flux ema-nating from the surface is that emanating from the white dwarf core. ForM? = 2?10?9 M yr?1, the lower branch exhibits significant hydrogen burn-ing indicating the onset of unstable hydrogen burning that occurs for loweraccretion rates. The second or middle branch for accretion rates in the range2 ? 10?9 M yr?1 ? M? ? 6.3 ? 10?7 M yr?1 is characterised by hydrogenburning, with the amount of hydrogen burning increasing with increasingsurface flux. As a result, temperatures in the layer exceed the core tem-perature and must decrease sharply near the bottom of the layer to matchthe core temperature. This branch increases in column density with surfaceflux for accretion rates 2?10?9 M yr?1 and 6.3?10?9 M yr?1. Whilst wehave not conducted a stability analysis as yet, we postulate that this branchis stable for these accretion rates, leading to periodic flashes without massloss. For accretion rates 6.3 ? 10?8 M yr?1 and 6.3 ? 10?7 M yr?1, thehydrogen branch becomes slightly convex. However, for the most part, thisbranch also increases in column density with increasing surface flux. It ispossible that for these accretion rates, the hydrogen branch leads to bothunstable and stable burning. However, it is more likely that the systemupon reaching this branch will relax towards the right stable slope, where?? ? 109.5 g cm?2.125Chapter 6. Conclusions, Future Work and Other RemarksIt is noteworthy that the hydrogen burning branch begins at approx-imately the same column density, 109.5 g cm?2, for accretion rates in therange 2 ? 10?9 M yr?1 < M? ? 6.3 ? 10?7 M yr?1. This implies that thecharacteristic accretion timescale will decrease with increasing accretion rateas ???1. The fact that this branch is generally stable for the range of ac-cretion rates discussed agrees quite well with the literature. Other authorsindicate that accretion rates smaller than those leading to steady burning,but larger than those leading to unstable burning 10?10?10?9 M yr?1 [20],[17], [23], result in stable hydrogen burning.We also identified a third branch, characterised by helium burning, forthe accretion rates 2 ? 10?9 M yr?1 ? M? ? 6.3 ? 10?7 M yr?1. Heliumburning increases with increasing surface flux in this branch, which has atrough-like shape. The left slope is a decreasing function of the surfaceflux, while the right branch increases in column density with surface flux.Most of the helium is burned upon reaching the turning point as a resultof a thermonuclear runaway in the left branch. It is plausible that the leftbranch leads to unstable helium flashes resulting in mass loss, whilst the rightbranch leads to stable helium flashes. Temperatures in this branch exceed108 K, and thus a temperature inversion occurs in the layer in order to matchthe inner boundary condition. The helium branch begins at lower columndensity with increasing accretion rate. Thus, the characteristic accretiontimescale required to ignite helium burning decreases faster than ???1. Itis reasonable to conclude that the time to initiate helium burning and therecurrence time increases with decreasing accretion rate. This agrees withthe results found by Jose? et al. for helium flashes [20].For accretion rates 6.3?10?12 M yr?1 ? M? ? 6.3?10?10 M yr?1, theshape of the solutions changes quite dramatically. There is only one branchcharacterised by hydrogen burning. The amount of hydrogen burning in thebranch and the total height of the branch in column density increases withdecreasing accretion rate. Furthermore, the branch decreases in columndensity with increasing surface flux. It is reasonable to postulate that thesebranches lead to unstable hydrogen burning or nova-like outbursts, and thatthe violence of these flashes increases with decreasing accretion rate. Thisagrees well with other authors, [20], [23], [17] and [4]. The branch beginsat higher column densities as the accretion rate is decreased. Thus thecharacteristic accretion timescale increases as we lower the accretion rate.This agrees with Jose? et al., where the authors found that decreasing theaccretion rate increased the initial time to flash and recurrent timescale forhydrogen burning [20]. We leave a detailed investigation of the behaviourof these solutions and those potentially including helium burning to future1266.1. Future Workwork.Our steady-state results agree quite well with previous literature. Infact, we can reasonably conclude that for accretion rates in the range 2 ?10?9 M yr?1 < M? ? 6.3 ? 10?7 M yr?1, stable hydrogen burning occurs,as well as unstable and stable helium flashes. It is plausible that the re-gions of stable hydrogen and helium burning result in mass increase, andpresent a set of parameters that could lead to progenitors of Type Ia su-pernovae. Accretion rates smaller than this lower limit, or in the range6.3? 10?12 M yr?1 ? M? ? 6.3? 10?10 M yr?1, result in unstable hydro-gen burning, most likely leading to nova-like outbursts. However, some ofour explanations were based on an inference of stability in the branches.To confirm our intuition, we must conduct a global linear stability analysis.The outline of such an analysis was presented in chapter 5. Our main goalfollowing this work is to test our steady-state solutions for stability.As we have already noted, the evolution of the accreted layer can betreated in great detail using intensive computer simulations. However, thisapproach is time-consuming, making the exploration of a large parameterset difficult. The approach employed in this thesis was based on a semi-analytical method, in which we determined the steady-state solutions andprovided a framework with which to determine stability. The numericalroutines employed in our analysis were not time-consuming, thereby allowingthe exploration of a large parameter set. Indeed, we were able to attain anunderstanding of the dependence of the steady-state solutions on a range ofaccretion rates, surface fluxes and surface CNO element mass fractions.6.1 Future WorkAs previously mentioned, we would like to conduct a global linear stabilityanalysis to confirm our intuitive beliefs of the stability of the steady-statesolutions. Once this analysis is complete, we would like to investigate thedependence of the burning regimes on the core temperature, the white dwarfmass and the chemical composition of the donor star. In particular, wewould like to use a range of masses and core temperatures Tcore < 4 ? 107K, to evaluate if and how the steady-state solutions change. We could alsoalter the input physics, to include, for example the hot CNO cycle withsaturation to evaluate if the range of accretion rates leading to unstable andstable burning changes.Narayan and Heyl investigated the problem of hydrogen-rich accretionon the surface of neutron stars [26]. In their analysis, they used an inner1276.2. Potential Observational Confirmations of Theoryboundary condition based on a diffusion timescale as outlined in chapter4. We would like to extend our model to utilise this boundary conditioninstead, and compare the results with those found in this thesis. We haveconducted an analysis utilising both methodologies for the inner boundarycondition for a neutron star. The reader is directed to the appendix for adiscussion. From this analysis, it is evident that the boundary conditionalters the shape of the solutions.We would like to implement a more accurate routine to identify thesteady-state solutions when the inner boundary condition is applied at thebottom of the layer. For example, we could employ a bisection method toevaluate if the steady-state solutions change to any noticeable degree. Inthe appendix, steady-state solutions are illustrated for a neutron star, whichwere found by using both a bisection method and our search routine. Theresults produced by both numerical routines agree incredibly well. However,we would like to investigate the effects in the case of a white dwarf.6.2 Potential Observational Confirmations ofTheoryLastly, and most importantly, we would like to compare our theoreticalfindings to observational data where possible. There are obvious difficultieswith obtaining observational data. Classical novae are expected to recur ona timescale of thousands to hundreds of thousands of years [36], [37]. Thus, ifa classical nova was recorded previously, it is unlikely that a second eruptionwill be observed in the near future. However, there are some aspects of ouranalysis that could conceivably be compared to observational data.The total fluence produced by nuclear burning is likely to be an observ-able quantity during a flash. In future work, we would like to calculate thetotal fluence, EH +EHe, by integrating the energy generated per gram overthe layer,?E?HX(?)+E?He(X(?)+Y (?)) d?. An effective upper limit to theflash duration can be calculated by dividing the fluence by the Eddingtonluminosity [26]. A time-dependent integration of the model equations willelucidate how much of the fuel is burned and emitted during a flash, whichis a calculation separate to the current semi-analytical analysis.Using near-infrared spectra, Das et al. estimated the mass ejected fromthe 2006 eruption of the recurrent nova RS Ophiuchi to be approximately3?10?6 M [7]. It is possible to determine the accreted mass in our analysis.Assuming that all of this mass is ejected during an unstable flash, we couldpossibly compare certain unstable flash regimes to the mass ejection estimate1286.2. Potential Observational Confirmations of Theoryof ? 10?6 M, based on theoretical estimates of the accretion rate for RSOphiuchi found by Hernanz and Jose? in [18].According to Townsley and Bildsten, approximately 30 non-magneticwhite dwarf systems exhibiting dwarf novae were observed during quies-cence, allowing for spectral detection of the photosphere and measurementsof the effective surface temperatures [38]. The authors used these observa-tions along with theoretical work to determine local time-averaged accretionrates [38]. In our analysis, we solve for the effective surface temperaturesusing the Stefan-Boltzmann law, for a fixed accretion rate. Based on theaccretion rates found by Townsley and Bildsten for a range of surface tem-peratures, we could conceivably compare our effective surface temperaturesat specified accretion rates to their findings. Despite the fact that theiranalysis was based on observations of dwarf novae systems, a comparisonwith our analysis may lead us to useful insights.There are a number of modifications that can be made to the inputphysics in our model. However, it is only through a comparison with obser-vational data that we can evaluate the suitability of these inputs. Whilst anumber of further investigations can be performed within our model frame-work, we believe that our analysis produces the well-known accretion regimesdiscovered by other authors reasonably well. We aim to further validate ourmodel with relevant observational data in the future.129Bibliography[1] P. Bergeron, M. T. Ruiz, and S. K. Leggett. The chemical evolutionof cool white dwarfs and the age of the local galactic disk. The Astro-physical Journal Supplement Series, 108:339?387, 1997.[2] L. Bildsten and A. Cumming. Hydrogen electron capture in accret-ing neutron stars and the resulting g-mode oscillation spectrum. TheAstrophysical Journal, 506:842?862, 1998.[3] Richard L. Burden and J. Douglas Faires. Numerical Analysis.Brooks/Cole-Thomson Learning, seventh edition, 2001.[4] S. Cassisi, I. Iben, Jr., and A. Tornambe?. Hydrogen-accreting carbon-oxygen white dwarfs. The Astrophysical Journal, 496:376?385, 1998.[5] D. D. Clayton. Principles of Stellar Evolution and Nucleosynthesis.The University of Chicago Press, 1983.[6] J. P. Cox. Theory of Stellar Pulsation. Princeton University Press,1980.[7] R. Das, D. P. K. Banerjee, and N. M. Ashok. A near-infrared shock wavein the 2006 outburst of recurrent nova RS Ophiuchi. The AstrophysicalJournal, 653:L141?L144, 2006.[8] E. Fehlberg. Klassische Runge-Kutta-formeln vierter und niedrigererordnung mit schrittweiten-kontrolle und ihre anwendung auf wa?rmelei-tungsprobleme. Computing, 6:61?71, 1970.[9] G. Fontaine, P. Brassard, and P. Bergeron. The potential of whitedwarf cosmochronology. Publications of the Astronomical Society ofthe Pacific, 113:409?435, 2001.[10] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb,P. MacNeice, R. Rosner, J. W. Truran, and H. Tufo. Flash: An adaptivemesh hydrodynamics code for modeling astrophysical thermonuclear130Bibliographyflashes. The Astrophysical Journal Supplement Series, 131:273?334,2000.[11] M. Y. Fujimoto. A theory of hydrogen shell flashes on accreting whitedwarfs. I. Their progress and the expansion of the envelope. The As-trophysical Journal, 257:752?766, 1982.[12] M. Y. Fujimoto. A theory of hydrogen shell flashes on accreting whitedwarfs. II. The stable shell burning and the recurrence period of shellflashes. The Astrophysical Journal, 257:767?779, 1982.[13] M. Y. Fujimoto and R. E. Taam. On the secular evolution of accret-ing white dwarfs and Type I supernovae. The Astrophysical Journal,260:249?253, 1982.[14] I. Fushiki and D. Q. Lamb. S -matrix calculation of the triple-alphareaction. The Astrophysical Journal, 317:368?388, 1987.[15] I. Fushiki and D. Q. Lamb. New insights from a global view of X-raybursts. The Astrophysical Journal, 323:L55?L60, 1987.[16] C. J. Hansen and S.D. Kawaler. Stellar Interiors: Physical Principles,Structure, and Evolution. Springer-Verlag New York, Inc., 1994.[17] M. Hernanz and J. Jose?. Accreting white dwarfs. Proceedings of Science(Supernova), 014, 2007.[18] M. Hernanz and J. Jose?. The recurrent nova RS Oph: A possiblescenario for type Ia supernovae. New Astronomy Reviews, 52:386?389,2008.[19] I. Iben, Jr. Hot accreting white dwarfs in the quasi-static approxima-tion. The Astrophysical Journal, 259:244?266, 1982.[20] J. Jose?, M. Hernanz, and J. Isern. Hydrogen and helium shell flashes onmassive accreting white dwarfs. Astronomy and Astrophysics, 269:291?300, 1993.[21] R. Kippenhahn and A. Weigert. Stellar Structure and Evolution.Springer-Verlag Berlin Heidelberg, 1990.[22] K. S. Krane. Introductory Nuclear Physics. John Wiley & Sons, Inc.,1988.131Bibliography[23] J. MacDonald. CNO abundances and the strengths of nova outburstsand hydrogen flashes on accreting white dwarfs. The AstrophyscialJournal, 267:732?746, 1983.[24] G. J. Mathews and F. S. Dietrich. The 13N(p, ?)14O thermonuclear re-action rate and the hot CNO cycle. The Astrophysical Journal, 287:969?976, 1984.[25] R. Narayan and J. S. Heyl. On the lack of type I X-ray bursts in blackhole X-ray binaries: Evidence for the event horizon? The AstrophysicalJournal, 574:L139?L142, 2002.[26] R. Narayan and J. S. Heyl. Thermonuclear stability of material ac-creting onto a neutron star. The Astrophysical Journal, 599:419?449,2003.[27] K. Nomoto. Accreting white dwarf models for Type I supernovae. I.Presupernova evolution and triggering mechanisms. The AstrophysicalJournal, 253:798?810, 1982.[28] B. Paczyn?ski. A one-zone model for shell flashes on accreting compactstars. The Astrophysical Journal, 264:282?295, 1983.[29] B. Paczyn?ski and A. N. Z?ytkow. Hydrogen shell flashes in a white dwarfwith mass accretion. The Astrophysical Journal, 222:604?611, 1978.[30] B. Paxton, L. Bildsten, A. Dotter, F. Herwig, P. Lesaffre, andF. Timmes. Modules for experiments in stellar astrophysics (MESA).The Astrophysical Journal Supplement Series, 192:3, 2011.[31] Anthony Peirce. Private communication, 2012.[32] S. Perlmutter and A. Reiss. Cosmological parameters from supernovae:Two groups? results agree. American Institute of Physics ConferenceProceedings, 478:129?142, 1999.[33] M. Salaris, I. Dom??nguez, E. Garc??a-Berro, M. Hernanz, J. Isern, andR. Mochkovitch. The cooling of CO white dwarfs: Influence of theinternal chemical distribution. The Astrophysical Journal, 486:413?419,1997.[34] S. L. Shapiro and S. A. Teukolsky. Black Holes, White Dwarfs, andNeutron Stars: The Physics of Compact Objects. John Wiley & Sons,Inc., 1983.132[35] K. J. Shen and L. Bildsten. Thermally stable nuclear burning on ac-creting white dwarfs. The Astrophysical Journal, 660:1444?1450, 2007.[36] E. M. Sion. White dwarfs in cataclysmic variables. Publications of theAstronomical Society of the Pacific, 111:532?555, 1999.[37] R. C. Smith. Cataclysmic variables. Contemporary Physics, 47(6):363?386, 2006.[38] D. M. Townsley and L. Bildsten. Measuring white dwarf accretion ratesvia their effective temperatures. The Astrophysical Journal, 596:L227?L230, 2003.[39] D. M. Townsley and L. Bildsten. Theoretical modeling of the thermalstate of accreting white dwarfs undergoing classical nova cycles. TheAstrophysical Journal, 600:390?403, 2004.[40] S.-C. Yoon and N. Langer. The first binary star evolution model produc-ing a Chandrasekhar mass white dwarf. Astronomy and Astrophysics,412:L53?L56, 2003.133Appendix ANeutron StarWe use the code developed for a white dwarf to conduct a short analysisof nuclear burning on accreting neutron stars. We compare the results pro-duced by our code to the results produced by a code that was developedfor a neutron star. The code was obtained from my co-supervisor, ProfessorJeremy Heyl. This code uses a bisection method in the surface flux and astandard Runge-Kutta 2 method that does not adjust the step size accord-ing to the local truncation error. Instead, the step size is adjusted througha very simple routine, making the numerical method susceptible to instabil-ities for some values of the input parameters. This code applies the innerboundary condition at the bottom of the accreted layer, which agrees withthe methodology used in our code.The necessary adjustments were made to the code received in order toproduce the same output as our code for comparison. We also made someadjustments to our code for the extended model. In particular, we replacedthe energy generation rate of the CNO cycle with that of the hot CNOcycle with saturation, given by equation (2.14). Furthermore, we changedthe sign of grav in our model to be consistent with the code received. Thesign change does not affect the solutions to any significant degree, with thegeneral features of the solutions remaining unchanged. This implies that theenergy generated through nuclear burning is far larger than that producedthrough compression or expansion in the model. This is consistent with theobservation in white dwarfs that the flux produced by compression of theaccreted material is orders of magnitude smaller than the flux produced bynuclear burning [35]. We also alter our expression for the entropy and surfaceboundary conditions slightly to be consistent with the code received. Lastly,the code received approximates the energy generation rate of hydrogen by =  exp(?10?4/X), (A.1)in which the X can be replaced by Y for the energy generation rate of heliumburning. If X or Y reaches 10?7, the energy generation rates are set to zero.In our white dwarf code, we set the energy generation rates to zero if theywere lower than 10?16 or X < 10?10. We adjust our methodology so as134Appendix A. Neutron Starto better compare the results produced by the codes. Note that this doesnot change the behaviour of the solutions for the white dwarf, but sinceX reaches zero at a shallower depth in the neutron star due to the largersurface gravity, this is a more efficient method of keeping the code stableand preventing X or Y from becoming arbitrarily small.We conduct the analysis for the following set of parameters,Xs = 0.7, Ys = 0.28, (A.2)XCNO = 0.8Zs, XCNO = 0.8Zs + (Z ? Zs), (A.3)Mns = 1.4 M, R = 10.4 km, (A.4)lacc =M?c2z(1 + z)LEdd, fout =FsFnuc, (A.5)log g = 14.3, Tcore = 1? 108 K. (A.6)The chemical composition, mass and radius are the same as those used inNarayan and Heyl [26]. The steady-state solutions for log lacc = ?1.4 aredepicted in figure (A.1). We plot the solutions against fout where fout isthe surface flux, Fs, divided by the total flux produced by nuclear burning,Fnuc, given by equation (4.16). Figure (A.1a) was produced using the codewe received and figure (A.1b) was produced using our code. It is clear thatthe solutions produced by both codes are consistent with each other. It isnoteworthy that the features of the solutions produced by the respectivecodes are dependent on certain inputs. For example, the solutions producedby the bisection method depend on the values of ? chosen. If we did not usea refined grid in log ? between log 9.8 and log 10, the features of the peakwould be obscured or missed entirely. Furthermore, the bisection routine infout terminates for a particular value of ? once a solution has been found,and thus this routine will never produce two solutions of fs for one value of?. In our code, the features are dependent on the grid size in fout and thetolerance used for the temperature matching. We match the temperatureto three decimal places and use a grid size of 8? 10?4 in fout. Despite thedifferences in the code, they produce nearly identical solutions.The solutions have some behaviour in common with those that werefound for a white dwarf. Specifically, there is a branch which is a decreasingfunction of fout, there is a turning point followed by a branch in which thecolumn densities increase with increasing fout. It is interesting that thesesolutions do not resemble the solutions in Narayan and Heyl [26], which werefound by matching the inner boundary condition with a diffusion timescale,for log lacc = ?1.25 or log lacc = 1.5. The salient features in our solutionappear at column densities between 105 and 1013 g cm?2 and the range of135Appendix A. Neutron Starlog fout for which there are steady-state solutions run from approximately?1.4 to ?0.1. The prominent features of the solutions in [26] range from 107to 108.5 g cm?2 in column density and ?2.3 to 0 in log fout. Both solutionsexhibit a local maximum near log fout ? ?0.4.It is evident that the shape of the solutions depend significantly on themanner in which the inner boundary condition is applied. To investigate thishypothesis further, we developed a code that applies the boundary conditionat ?diff , where the diffusion timescale, determined by the integration ofequation (4.24), is twice the accretion timescale. We discussed the procedureof implementing this boundary condition in chapter 4. In particular, theroutine uses an iterative method to find solutions. Specifically, for each valueof ?bottom, the code iterates through each value of the flux, fout, integratesto ?bottom, then into the outer crust until ?diff = 2tacc. The routine thendetermines if the temperature condition is satisfied at a depth ?diff . Thenumerical scheme matches the diffusion timescale to one decimal place andthe temperature condition to at least one decimal place.We use the same parameters and equations as those used in our ex-tended model, and replace the energy generation rate of the CNO cyclewith equation (2.14) for the energy generation rate of the hot CNO cyclewith saturation. To integrate into the stellar substrate, the neutron starouter crust in this case, we set X = Y = 0 in the equation for the meanmolecular weight (2.7). This is an adequate approximation to the meanmolecular weight of the neutron star outer crust, which is composed of 56Fe[26]. In actuality, however, the outer crust will also consist of the productsof Type I bursts. The solutions are depicted in figure (A.2). The solutionswere produced using two different grids, one between log fout = ?2.5 to 1.2and a less refined grid from log fout = ?1.2 to 0.We note that matching the boundary condition in this manner producesthe same qualitative behaviour as that for log lacc = ?1.25 in [26], whichis different to the behaviour found by matching the boundary conditionin the accreted layer. In particular, there are two peaks. According to[26], these correspond to the hydrogen and helium burning peaks. Thetemperature at the bottom of the layer increases with increasing fout. Atapproximately 107.7 K, hydrogen burning begins [26]. The hydrogen peakappears shortly after this temperature is reached in fout. Thus, the solutionsare accompanied by a flux reversal just past the hydrogen burning peak, sothat the temperature drops in the substrate to match the inner boundarycondition.We conclude that the behaviour of the steady-state solutions is directlydependent on the implementation of the inner boundary condition. We136Appendix A. Neutron Starfurther conclude that the code used to produce the white dwarf steady-stateresults has been independently corroborated through the comparison of theresults for the neutron star, figures (A.1a) and (A.1b), produced by the codereceived from J. Heyl and our code respectively. In future work, we intendto determine solutions for the white dwarf using the diffusion timescale toimplement the inner boundary condition.137Appendix A. Neutron Star?1.4 ?1.2 ?1 ?0.8 ?0.6 ?0.4 ?0.256789101112131415log??( gcm?2)log foutlog l acc = ?1 .4(a) Variation of column density with log fout, results pro-duced with code provided by J. Heyl.?1.4 ?1.2 ?1 ?0.8 ?0.6 ?0.4 ?0.256789101112131415log??( gcm?2)log foutlog l acc = ?1 .4(b) Variation of column density with log fout, results pro-duced with code developed for white dwarf.Figure A.1: Comparison of steady-state results using two different codesfor the neutron star with Xs = 0.7, Ys = 0.28, XCNO = 0.8Zs and g =1014.3 cm s?2.138Appendix A. Neutron Star?2.5 ?2 ?1.5 ?1 ?0.5 05.566.577.588.5log??( gcm?2)log fout  log l acc = ?1 .4Figure A.2: Variation of column density with log fout, with Xs = 0.7, Ys =0.28, XCNO = 0.8Zs and g = 1014.3 cm s?2139

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0074134/manifest

Comment

Related Items