Mass Accretion in Cataclysmic Variable Stars

Compute and visualize the impact point of stellar mass transfer streams onto accretion disks using event detection

Cataclysmic variable stars often involve binary stars where mass is being accreted onto a degenerate star (often a white dwarf star) by a primary donor star. The mechanics of these systems can be modeled using the restricted three body problem. The shape of the donor star is constrained and defined by the shape of the gravitational potential of the system. Event detection can be used to determine when and where the mass transfer stream impacts the accretion disk.

First, initialize the system with masses in solar masses, mass ratio, orbit period and Coriolis term:

In[1]:=

M1=1.37;M2=1.12;q=1.12/1.37;μpotential=M2/(M1+M2);a=1(*normalizedforsimplicity*);G=1(*normalizedforsimplicity*);Porb=Sqrt[4*Pi^2*(a^3/(G*(1+q)))];Ωorb=2*(Pi/Porb);

Define the effective gravitational potential of the system:

In[2]:=

Φ[x_,y_,z_]:=-((G*M1)/Sqrt[(x+q/(1+q))^2+y^2+z^2])-(G*M2)/Sqrt[(x-a+q/(1+q))^2+y^2+z^2]-(1/2)*Ωorb^2*((x-μpotential+(q/(1+q))*a)^2+y^2);

Solve for critical Lagrange points, L1, L2 and L3 and the potential at that location:

In[3]:=

deriv=D[Φ[x,0,0],x];posAtL1=FindRoot[deriv0,{x,.3}];ΦAtL1=Φ[x/.posAtL1,0,0];posAtL2=FindRoot[deriv0,{x,1.5}];ΦAtL2=Φ[x/.posAtL2,0,0];posAtL3=FindRoot[deriv0,{x,-1}];ΦAtL3=Φ[x/.posAtL3,0,0];

The position of the white dwarf is found at {-μ,0,0}:

In[4]:=

μ=mu/.Quiet[Solve[qmu/(1-mu),mu],Solve::ratnz]〚1〛

Out[4]=

0.449799

The surface of the red giant is generated using the following:

In[5]:=

grL1=ContourPlot3D[Φ[x,y,z],{x,x/.posAtL1,x/.posAtL2},{y,-0.7,0.7},{z,-0.6,0.6},Contours{ΦAtL1},ContourStyle{RGBColor[1,.3,0]},MeshNone,Lighting{AmbientLight[GrayLevel[.7]],PointLight[White,{-μ,0,0}]}]

Out[5]=

The accretion disk is represented using a textured 3D disk centered on the white dwarf at {-μ,0,0}:

In[6]:=

texture=ColorReplace

,BlackRGBColor[0,0,0,0];

In[7]:=

pdisk[mu_,rad_]:=ParametricPlot3D[{rCos[t]-mu,rSin[t],0},{r,.02,rad},{t,0,2π},PlotPoints75,PlotRange{Automatic,Automatic,{-.1,.1}},ImageSize600,PlotStyleTexture[texture],TextureCoordinateFunction({#1,#2}&),MeshNone,Lighting{DirectionalLight[GrayLevel[.35],{0,0,1}],DirectionalLight[GrayLevel[.35],{0,0,-1}]}]

Define the initial velocity of the mass transfer stream in the x and y directions:

In[8]:=

Vx0=-0.01;Vy0=0.00333300000000002;

Solve the differential equations of motion and use event detection to determine when the mass transfer stream impacts the accretion disk:

In[9]:=

event[x_,y_,t_,drad_]:=WhenEvent[N[Norm[{x[t]+μ,y[t],0}]]drad,Sow[t];"StopIntegration"]

In[10]:=

eqs[x_,y_,t_,drad_]:=

′′

[t]2

′

[t]+x[t]-

(1-μ)(x[t]+μ)

3/2

(

(x[t]+μ)

y[t]

)

μ(x[t]-1+μ)

3/2

(

(x[t]-1+μ)

y[t]

)

′′

[t]-2

′

[t]+y[t]-

(1-μ)y[t]

3/2

(

(x[t]+μ)

y[t]

)

μy[t]

3/2

(

(x[t]-1+μ)

y[t]

)

,y[0]0,x[0]Evaluate[x/.posAtL1],

′

[0]Vy0,

′

[0]Vx0,event[x,y,t,drad];

In[11]:=

solution[drad_]:=Reap@NDSolve[eqs[x,y,t,drad],{x,y},{t,0,13}]

Extract the time of the impact from the solution:

In[12]:=

diskimpacttime[drad_]:=Last[solution[drad]]〚1,1〛;

The point where the mass transfer stream impacts the accretion disk is the location of a bright spot as the infalling gas changes velocity and is shocked:

In[13]:=

bspot[drad_]:=Flatten[{Map[#[diskimpacttime[drad]]&,({x,y}/.solution[drad]〚1〛)〚1〛],0}];

Here is a graphics representation of the mass transfer stream:

In[14]:=

transferstream[drad_]:=ParametricPlot3D[Evaluate[{x[t],y[t],0}/.solution[drad]〚1〛],{t,0,diskimpacttime[drad]},PlotStyleDirective[RGBColor[1,.3,0],Thickness[.003]]];

Define the orbital inclination of the system:

In[15]:=

inclination=67;

We can visualize the system by specifying the orbital phase. The white dwarf is tiny, on the order of.01 solar radii, so would be difficult to see if rendered to scale, so its size is magnified here:

In[16]:=

mypr={{(x/.posAtL3),x/.posAtL2},{-.5,.5},{-0.6,0.6}}

Out[16]=

{{-1.27938,1.3102},{-0.5,0.5},{-0.6,0.6}}

Define the radius of the accretion disk in scaled units of the semimajor axis:

In[17]:=

diskrad=.06;

Setup the 3D scene:

In[18]:=

visualize[orbtime_,drad_]:=Graphics3D[{{Glow[White],Sphere[bspot[drad],.01]},{White,Lighting{AmbientLight[GrayLevel[.1]],DirectionalLight[White,ImageScaled[{0,0,1}]]},Sphere[{-μ,0,0},.005]},grL1〚1〛,Lighting{DirectionalLight[GrayLevel[.35],{0,0,1}],DirectionalLight[GrayLevel[.35],{0,0,-1}]},pdisk[μ,drad]〚1〛,transferstream[drad]〚1〛},BackgroundBlack,BoxedFalse,AxesFalse,ImageSize{600,400},ViewPoint{4Cos[-2πorbtime]Sin[inclinationDegree],4Sin[-2πorbtime]Sin[inclinationDegree],4Cos[inclinationDegree]},PlotRangemypr,SphericalRegionTrue,ViewAnglePi/25,Lighting{AmbientLight[GrayLevel[.3]]}]

The mass transfer stream leaves the L1 point and falls towards the compact primary star where it impacts the accretion disk (indicated with a white sphere) and merges with it:

In[19]:=

visualize[.25,diskrad]

Out[19]=

Compute the solution again for a larger accretion disk radius:

In[20]:=

diskrad=.4;

In[21]:=

visualize[.25,diskrad]

Out[21]=

Publisher Information

Contributed by: Jeff Bryant

Mass Accretion in Cataclysmic Variable Stars

Related Symbols

Publisher Information