The number of radiative force is,

f_{rad}/f_{g}

The first row is the top views, the second row is the side view.

Central star mass is 0.4 solar mass

R coordinate is from 0.8 au - 4 au

Inner gas temperature is 2E4 K and the hydrogen atoms are all ionized.

Gas is a mixture of hydrogen and helium with X=0.7 Y=0.3

movie.gif​

1. Ionized parker wind with recombination. X=0.7 Y=0.3, therefore \mu=1.29 is fully molecular state and \mu=0.678 is the fully ionized hydrogen and molecular helium state.

m=0.1 m_{sun}

movie.gif​

2. I am writing the code that solves the 1D Eulerian radiation hydrodynamics with flux-limited diffusion radiation transfer combined with ray-tracing radiation transfer and a general equation of state hydrodynamics. It will also have dust formation, frequency dependent radiation transfer, and dust destruction. I have successfully combined the two radiation transfer method in optically thick and optically thin region. One week away from the completion of the code.

A steady state wind.

movie.gif​

​Riemann solvers

There are two simulations distinguished by the momentum in flux: higher momentum flux and lower momentum flux.

\rho (g\cdot cm^{-3}): density of the outflow

v_{r}(km\cdot s^{-1}): velocity of the outflow, spherical

lower momentum flux

higher momentum flux

convergence.pdf​

1drhd.pdf​

I put a 1 solar mass object at the origin of the simulation box. There is no outflow from 0 to 6000 day. I list the physical quantities during the quiet phase below:

Radius of the object that has outflow: r=46.875R_{\odot}

Mass of the object: 1M_{\odot}

Density of the outflow: \rho=10^{-13}g\cdot cm^{-3}

Temperature of the object: T_{1}=30000K

After the quiet phase, there is an outflow. An important constraint is the momentum budget. It is related to the luminosity of the object by: \dot{p}\le\frac{L}{c}. Model ID 6 is the old simulation with only 500 days of quite phase.

I think the implied luminosity is still too high for the outflow. Also, it remains a question that what kind of driving mechanism is this outflow? Model 7 seems to be a nova outburst and model 6 seems to be a supernova….

I put a 1 solar mass object at the origin of the simulation box. There is no outflow from 0 to 500 day. From 500 to 1000 day, there is outflow for the boundary of the object. I list the physical quantities below:

Radius of the object that has outflow: r=46.875R_{\odot}

Mass of the object: 1M_{\odot}

Density of the outflow: \rho=10^{-11}g\cdot cm^{-3}

Temperature of the object: T_{1}=30000K

The outflow is asymmetric, in the polar direction (theta=0 and theta=pi), the outflow velocity is the maximum; in the equatorial plane, the outflow velocity is 0. The outflow velocity is a function of polar angle, there are two kinds of functions, linear and Gaussian - to vary the open angle. The Gaussian is narrower.

Linear: v_{r}(\theta)=v_{max}\times\frac{\|\pi-2\theta\|}{\pi}

Gaussian: v_{r}(\theta)=v_{max}\times\exp(-10(\frac{\|\pi-\|2\theta-\pi\|\|}{\pi})^{2})

Spherical: v_{r}(\theta)=v_{max}

Below shows the linear and Gaussian velocity v.s. angle.

The outflow temperature: 30000K

Outflow duration: 500 day - 1000 day

Mass loss: \dot{M}=\int_0^{\pi}2\pi\rho r^{2}\sin\theta v_{r}(\theta)d\theta. However, much of the outflow will be impeded by the surrounding gas, the actual mass loss will be lower than the calculated ones.

ID	v_{r}(0)(km\cdot s^{-1})	\rho (g\cdot cm^{-3})	\dot{M}(M_{\odot}\cdot yr^{-1})	simulation	approximate expansion velocity
1	300 (linear)	10^{-11}	2.31\times10^{-4}		120(km\cdot s^{-1})
				movies:​density and temperature
2	400 (linear)	10^{-11}	3.07\times10^{-4}
3	400 (gaussian)	10^{-12}	1.00\times10^{-5}
4	400 (gaussian)	10^{-11}	1.00\times10^{-4}		310(km\cdot s^{-1})
				movies:​density and temperature
5	300 (spherical)	10^{-11}	6.35\times10^{-4}		240(km\cdot s^{-1})
				movies:​density and temperature

Simulations that has indefinite outflow duration. (it is always on)

Some questions: where exactly does the jet emerge? What are the open angle and outflow velocity. What are the outflow density and temperature?

1. Previously, we have this version of paper on the orbital period change in binary system with giant stars.

​paper

The referee asked us to show that our angular momentum evolution do not have serious problem. So I sampled the angular momentum through multiple of spherical shells in a single AGB star test in the co-rotating frame. Ideally, the angular momentum of the wind at any shell is the same, or not varying if we average of enough time. However, simulation has errors and sampling is never infinitely long.

​angular momentum error study

We should also replace the five images of alpha in the appendix in the last version of the paper.

​3au

​4au

​6au

​8au

​10au

We are going to use two figures to show the corrected orbital period change.

​exclude all angular momentum which measured at 0.8d in the AM error study

​only exclude the spin angular momentum of the AGB star

I am rewriting the paper.

2. I studied the Riemann problem with realistic EOS.

​A study of EOS Riemann problem

The conclusion is, realistic gas behave somewhat like ideal gas with a smaller gamma, though not exactly the same. We can in the future run some simulation with gamma=1.1 instead of gamma=1.67.

I recommend to remodel the AGB wind and binary star model with this argument.

In this simulation, I placed an outflow sphere at the center. The boundary condition is 30000K, 1e-13 g/cc and 300km/s from 0 to 100 day (mass loss rate is 7.2e-6 solar mass / yr) at the radii of 50 R_sun then the outflow speed changes to 10km/s. This outflow is not doing much at this timescale.

I launched another outburst from 1000 to 1100 day with 30000K, 1e-10 g/cc and 1000km/s. Thus the mass loss rate is 0.024 solar mass / yr.

​topdensity

​sidedensity

It is unphysical to force the initial condition to be the initial condition that Thomas (Orsola's student) sent me since the density in SPH is sampled via kernel function with its own radius but AstroBEAR has another mesh size. In short, the density at the central point is sampled within a small volume in SPH but our finest grid size is huge compared to that volume. If nothing done to the initial condition, the core will explode since it has around 3.5 solar mass. The density scale is in computational unit (1e-12 g/cc)

Explode, topdensity​

Explode, sidedensity​

After carving out the core with a replacement of 50R_sun radius, 10000 Kelvin, 5E-9 g/cc density and 0 outward velocity spherical boundary. The density is in physical unit (g/cc). Sorry that I did not make two sets of movies in the same scale.

​topdensity

​sidedensity

​topdensityzoomin

The flow looks more smooth and outflow is less violent. The flow structure is preserved within 6 years of simulation.

Initial condition: from Orsola's SPH result.

Boundary condition: a spherically symmetric outflow set at r=156R_sun, the outflow density is 5e-11 g/cc, outflow speed ram from 10km/s to 1000km/s. The duration of the 1000km/s outflow is from time=1day to time=40day. The mass loss rate is 1e-1 solar mass per year and the total mass loss is about 1e-2 solar mass. I have also added a 2 solar mass object at the center of the outflow.

Simulation: I use gamma=1.667 in this simulation. The maximum temperature during the late simulation is about 5e7 Kelvin which is unphysical. That high temperature happens at low density region where the shock heating is strong. I use 120 cores for 25 hours, that is 3000 CPU hours. It is static mesh refinement.

side density zoom in​

side density​

top density​

​Drive an outflow

​A new class of accurate, mesh-free hydrodynamic simulation methods

​Orbital evolution in binary systems with giant stars

Binary torque

​Angular momentum loss and mass-exchange instability in binary stars

​An Eccentric Circumbinary Accretion Disk and the Detection of Binary Massive Black Holes

Eccentricity

​Orbital evolution of mass-transferring eccentric binary systems. I. Phase-dependent evolution

​Orbital evolution of mass-transferring eccentric binary systems. II. Secular Evolution

Soker

​Forming equatorial rings around dying stars

Sweden

​Mass transfer in white dwarf-neutron star binaries

I have setup the components that will be used in this simulation. The figures show the mesh and the boundary of the outflow.

In this simulation, we are going to use static-mesh-refinement. Adaptive-mesh-refinement can be applied when needed. This is because we know where need to be resolved.

The solid sphere in the center of the figures is the outflow boundary. Its radius is roughly 160R_sun.

The next question will be:

1. What outflow do we need? (velocity, gravity, temperature, density)

2. Which machine should we use to carry out this simulation? It is a huge simulation. I estimated that to generate 1 year of simulation, we need roughly 1000 CPU hours.

​Three-dimensional hydrodynamical models of wind and outburst-related accretion in symbiotic systems

​Playlist

3AU, 0 degree and 90 degree

4AU, 0 degree and 90 degree

6AU, 0 degree and 45 degree

8AU, 0 degree and 45 degree

10 AU model will see separating binary stars and others will have orbital decay. 3 AU model will shrink 10% of the separation in 4000 yrs.

Thomas, Professor De Marco's student, communicated with me before. Jonathan and I took a look at their data this morning and find that the effective resolution is 32000. The size of the simulation box is (128000 R_{sun})^{3} and the finest cell is (2.9 R_{sun})^{3}. It need at least 6 levels of AMR.

This job is very computationally intense.

​3au

yproj.gif​

zproj.gif​

singld3dside.gif​

singld3dtop.gif​

burstsidezoomout.gif​

bursttop.gif​

Test an isolated AGB star burst. Below is a plot of its mass loss rate (solid line) measured at 10 AU radii and the initial burst velocity (dash line).

​Staff et al. (2016) calculated the 3-D mass loss rate (ejection) in common envelop evolution.

topdensity.gif​

4aubursttop.gif​

4auburstside.gif​

I can write a conference proceedings (poster).

I got the 1st student award during the meeting.

in attachments.

​https://astrobear.pas.rochester.edu/trac/wiki/u/zchen/summary

This post analyze the 4 AU binary simulation with great pulsation in the co-rotating frame. The pulsation is from t=15yrs to t=25yrs which last for 10 years. The wind speed during the pulsation is 20km/s and the escape speed is 40km/s at that radii.

This image shows the density looking from top. There is a spiral structure that extends outward indefinitely. This structure may be a consequence of the great pulse or just the normal AGB pulsation. But a direct reason for this structure is that the secondary can not accrete all the gas that has been pulled to it. Therefore the secondary will periodically release some gas and absorb the rest.

I do not think the spiral structure is starting from the L2 point. Because the L2 point is not located on the high density arm that extends outward. As I can see, the gas that leaking from the L2 point is gravitationally bound.

This image shows another information. First I will define a quantity, angular momentum per unit mass:

j_{b}=\frac{m_{agb}v_{agb}d_{agb}+m_{com}v_{com}d_{com}}{m_{agb}+m_{com}}

where subscript b stands for binary. d is the distance to the center of mass and v is the velocity with respect to the center of mass in the lab frame. This value for such system is 2.425\times10^{19}g\cdot cm^{2}/s.

Then we calculate the actual angular momentum per unit mass in the simulation and name it j_{g}. This image shows the value of s1=j_{g}/j_{b}. The greater the s1 the higher the angular momentum per unit mass. We can see that this value is above 4 in most area in this image.

Next I calculated the mass of the gas that leave this system, it is \dot{M}_{runaway}=2.16652\times10^{-7}M_{\odot}/yr. The accretion rate of the secondary is \dot{M}_{sec}=1.037\times10^{-7}M_{\odot}/yr. Therefore the massloss rate of the AGB star is M_{agb}=3.2\times10^{-7}M_{\odot}/yr. In the single AGB star case, the massloss rate is around 1.6\times10^{-7}M_{\odot}/yr. So the existence of the secondary double the massloss rate of the AGB star.

I also calculated the angular momentum loss rate, that is, the angular momentum in the z direction that is carrying away by the runaway gas per unit time. It is \dot{J}_{z}=1.52845\times10^{39}g\cdot cm^{2}/s^{2}. Hence the angular momentum per unit mass in the z direction in the runaway gas is:

j_{g}=\frac{\dot{J}_{z}}{\dot{M}_{runaway}}

compared to j_{b}, it is:

\frac{j_{g}}{j_{b}}=4.613

Now, we can see that the runaway gas is carrying more angular momentum than the average angular momentum in the binary system. More specifically, it is 4.613 times of the average angular momentum per unit mass.

If we write j_{b} in terms of m_{agb},m_{sec} and d, it is:

j_{b}=\frac{m_{agb}m_{sec}sqrt(Gd)}{(m_{agb}+m_{sec})^{3/2}}(1)

On the other hand, by the conservation of angular momentum:

j_{b}(m_{agb}+m_{sec})+j_{g}m_{g}=j_{b0}(m_{agb0}+m_{sec0})

where 0 stands for the initial values.

If we let this angular momentum and accretion run for 5\times10^{5}yrs under such condition, then j_{b}=0.727 j_{b0} at that time. Use equation (1) and substitute the correct numbers we can get that the separation will be 2 AU.

Researchers have put much attention to the angular momentum carried away by the gas leaking from the L2 point. The angular momentum loss in this simulation is different from that and it is more effective. However, I can not tell if this situation is going to last because it has only been 9 years after the pulsation. However, considering that the massloss rate for the single AGB star is low. I would argue that if the massloss rate of the AGB star is very high, such expulsion may last because the secondary can not accrete that much of gas. I will run the 4 AU simulation without the great pulsation and check whether if will also be there.

I think when there is no formation of circumbinary disk, there is strong orbital decay.

topdensity.gif​

topdensityzoomin.gif​

sidedensity.gif​

I have done the 6 AU binary simulation. The circumbinary disk forms after I enlarge the AMR region. We can not see the formation of circumbinary disk and the accretion disk in low resolution simulation.

topdensity.gif​

topdensityzoomin.gif​

sidedensity.gif​

Single star in corotating frame. This simulation only take 5 hours on 120 cores. It is a 48AU^3 box and the base grid is 60^3. I use 3 level of AMR. So the finest grid is 0.05AU^3.

singletopdensity.gif​

singlesidedensity.gif​

The orbital period of 6 AU simulation is 12 years, so almost 3 orbits has been done.

6autopdensity.gif​

6ausidedensity.gif​

The density contour of 3 AU separation simulation shows that something very similar to common envelop forms. The mass transfer mode between the two stars is between RLOF and WRLOF. I believe it is closer to the RLOF. The mass is flowing out through the L2 point. The orbital period of 3 AU simulation is 4.96 years. So more than 3 orbits has been done. The result has some similarity to this paper ​Sawada, Matsuda & Hachisu(1986). In that paper, they call the bridge that connects the two stars the elephant trunk.

3autopdensity.gif​

3ausidedensity.gif​

I tested two sets of data.

model 1 has extremely large mass loss rate - 1.5\times10^{-4}M_{\odot}/year and model 2 has 3\times10^{-5}M_{\odot}/year. In previous model, we can just achieve 1.5\times10^{-7}M_{\odot}/year because I may set the abseff a little higher (around 0.45?) and the resolution is poor. Also, we have changed the CFL number. Current CFL number is more appropriate.

Below is the mass loss rate in 3D single AGB star simulation.

We may use model 2 in our simulation. If the mass loss rate is still thought to be too high. I can lower the abseff to 0.45 and surface density to 8\times10^{-9}g/cc.

I can control the initial condition to avoid the strong shock now.

I can also control the variables outside the boundary cylinder to avoid in falling material and shocks.

Below are the simulation from model 1:

topdensityhml.gif​

sidedensityhml.gif​

Below are the simulation from model 2:

topdensitylml.gif​

sidedensity.gif​

Below is model 3:

I found that the resolution of radiation field is very important.

I have run binary simulation with low radiation field resolution and high resolution. The low resolution one has Bondi-Hoyle accretion but the high resolution one has Wind Roche-lobe accretion. The difference comes from the radiation force on the gas. The lower resolution one has higher radiation force on gas thus drive a faster wind. The difference in radiation force is tiny and the terminal velocity should not be too much different (which I will go back and check the single 3D AGB star wind).

Below is the low resolution run. It is faster and more stable than the high resolution run. The image show the side view of the density, velocity and Mach number contour. The result is pretty much in line with previous research on BH accretion ​Shima et al. (1985)​Ruffert (1996). More specifically, it resemble the result of low gamma and high Mach BH accretion. This is very reasonable because the cooling in our model is in fact equivalent to decreasing the gamma.

topdensity.gif​

topdensityzoomin.gif​

sidedensityzoomin.gif​

The temperature of outer region (near boundary) increases unphysically. This can be the inappropriate boundary condition in the co-rotating frame.

toptemp.gif​

sidetemp.gif​

Below is a snapshot of high resolution result. This run also has problem. But it shows that there is a bridge which resemble Roche-lobe accretion. However, the position of this bridge is not at the L1 point (not even close, probably 0.8~1.0 AU away). This may be a sign of wind Roche-lobe overflow.

Actually, when go from Roche-lobe to BH wind accretion. There should be a separation that the bridge appear and disappear as the AGB star pulsate. That state can distinguish the RLOF and BH accretion and be a sign of WROLF.

1. I am half way to calculate a new EOS. It will consider the balance between: electron, hydrogen atom, hydrogen ion, hydrogen molecule, helium atom, HeII and HeIII.

diagnose.gif​

2. I updated the model 2 and model 4 in ​summary

It turns out that model 4 will have unphysical result.

1. Running model 3 now. Model 2 has finished. I will post model 2 and model 3 result on Friday.

2. I am thinking about how to analyze the results. Also, three dimensional single AGB star is not so symmetric (see the wind velocity profile ​summary). Implying that there should be something not right in the code.

I will run model 3, model 4 and model 5.

q is the mass ratio of the two stars, defined by M_{com}/M_{AGB}

z is the tidal force coefficient, defined by (\frac{M_{com}}{(d-r_{AGB})^2}-\frac{M_{com}}{(d+r_{AGB})^2})/\frac{M_{AGB}}{d^2} where d is the binary separation.

​linktosummary

1. Jonathan helped in debugging. The code can restart on bluehive without producing error.

2. I ran a 4 AU simulation this week. This simulation has a pulsating AGB star which described before. The difference is: I injected an outflow during 16.83 year and 26.83 year. So the total outflow time window is 10 years. Then the AGB star is set back to normal pulsating state. The result is very exciting. Below are three density plot and contour plot from the simulation. The burst of the outflow has ended and the AGB star is back to pulsating state. The speed of the outflow is 20km/s (remember the pulsating amplitude is only 5km/s), and the escape velocity (with respect to the AGB star only) where the outflow is ejected is about 40km/s. However, we can not treat the gas as ballistic matter. The pressure and temperature will push much of the gas out.

A large amount of the gas is gravitationally bound. I mistakenly overwrite the result of later time. A new simulation is running. I saw the late time evolution, a smooth circumbinary disk will form around the binary. (click the image to enlarge)

A disk forms around the secondary. (click the image to enlarge)

3. Baowei helped me to set up the code on Bluestreak. I can use 8192 cores on bluestreak. The 3 AU, q=0.1 simulation is schedules to run on bluestreak.

1. Plan for runs

1. I was improving the code's efficiency and debugging. I restrict the highest level only to the AGB star and the secondary. The radii of dust formation zone is less than 2.4 AU. The finest cells consist of two parts, one is the cells that are within the 2.4 AU sphere of the AGB star and the other part is the cells that are within the cylinder which has 2.4 AU height and 3 AU radius. AMR is set to capture the spiral shock, too.

2. The current simulation is: 96AU*96AU*48AU physical size with 5 level of AMR. The base resolution is 60*60*30. The coarsest cell is 1.6 AU and the finest cell is 0.05 AU. The AGB star has 18 cells across its radius and 21 cells across its photosphere. The simulation will take 8 days on 120 cores, correspond to 12 orbits, or 72 years.

3. The 4 AU separation run is still going on.

1. About the simulation zone. The physical size of the current simulation has 48au^3. However, the hydrodynamic behavior suggest that the gas will go out further and will take a long time to fall back. I estimate that a minimum physical size would be 96au^3 for 3au and 4au binary separation. The simulation zone may still need to be larger for 6au and 10au separation. The physical size of the finest cell should be 0.05au or less. Therefore, a 4 level AMR with 120^3 base grid can handle the simulation. The time should take (estimated) more than 15 days on 120 cores (3 level AMR is 5 days). If I can run with 2*120 cores on bluehive. I think the simulation can be done by the end of August.

2. The pulsation period is 1 yr. This is the highest macroscopic frequency in this simulation (compared to orbital period and/or any other). To resolve the mass loss rate, I need at least 8 outputs per pulsation period. Therefore if the simulation is 300 years. The output will have 2400. Each output will take 200 MB at least. So the total disk space needed for one simulation is at least 500 GB. Currently I have 1200 GB on bluehive.

3. Run a 4au simulation with high resolution. 96au*96au*48au

1. High mass loss rate AGB star need a large atmosphere. Its photosphere will be around 1.5au. ​Liljegren et al. (2016) use 1.5au AGB and achieved 3\times10^{-6}M_{\odot}/yr. However, 1.5au will be too large for us. The tidal force will be 50% of the gravitational force from the AGB star. The highest mass loss rate I can do is 10^{-7}M_{\odot}/yr. Below is the mass loss rate in 3-D simulation. I sum the mass flux through a 10 AU spherical shell. The luminosity of this AGB star is 2343L_{\odot} and mass is 1M_{\odot}. So the mass loss rate and size is between RGB and super AGB.

2. A detailed description of this AGB star model is:

AGB star model

There are two most important aspects in binary star simulation. One is the radiative transfer and the other is the AGB star model. The first one determine the large scale motion and morphology while the second one determine the robustness of inner boundary.

A sketch of the AGB star's structure is shown in the picture.

The hard pulsating sphere has velocity variation. Its amplitude is 5km/s and its period is 1 year.

The radius of the photosphere is also fixed. In reality, it must change but it is difficult to calculate. Within the photosphere, reduced gravity (80% of the original gravity) is being used. The reduced gravity has its physical stand - higher opacity inside the star but it is also because the resolution is poor. The simulation is resolving 2~3 order of magnitude density drop within 3 cells.

Between the photosphere and the dust formation zone, it is the gas opacity (2\times10^{-4}cm^{2}/g). The dust formation radii is calculated by

R_d=\sqrt{\frac{\alpha L \exp{(-\tau)}}{4 \pi \sigma T_{sub}^{4}}}

Where \tau=\int_{rphotosphere}^{rmax} \rho\kappa_{gas}dr is first calculated by assuming that there is only gas in the simulation. When the R_d is determined, I update the opacity where there is dust. The dust opacity is 8cm^{2}/g.

A questionable number is \alpha=0.3.

3. A 3 AU separation simulation now will take about 5 days on 120 cores. I tested the resolution on single AGB star and it has a difference. The size of the simulation zone is 48au^3. The previous resolution is 120^3 base grid with 2 level of AMR. The current AMR is 3 while the base grid resolution is the same. The finest physical resolution is 0.05au thus the AGB star has 18 cells across its radii. The photosphere has 21 cells.

4. I calculate the strength of tidal force by:

(\frac{M_{com}}{(d-r_{AGB})^2}-\frac{M_{com}}{(d+r_{AGB})^2})/\frac{M_{AGB}}{r_{AGB}^2}

where d is the binary separation. Substitute d=3au, r=1.5au, M_{AGB}=M_{\odot} and M_{com}=0.5M_{\odot} you will get the value is close to 0.5 which is 50%.

5. Currently I have some 3au and 4au simulations.

1. 3au simulation:

I use T_{eff}=2900K or 2343L_{\odot} at the first 23 years. The AGB star is the same as in the single star case. After the simulation has come to a pattern state - the spiral shock escapes from the system. Below shows the pattern. Left figure shows the density plot and the right figure shows the temperature plot.

After 23 years, I make the AGB star dimmer by setting T_{eff}=2600K or 1513L_{\odot}. The spiral shock is not escaping in this situation but the simulation box is too small so we can not see its fall back.

The movie is here:

densitytop3au.gif​

2. 4au simulation:

Interestingly, the 4au separation binary star will have non-escaping spiral shock when I put the normal AGB star (T_{eff}=2900K) in the simulation. Also, the simulation box is too small.

The movie is here:

densitytop4au.gif​

It seems that 3au simulation has Roche lobe overflow while 4au simulation has wind Roche lobe overflow?

1. I discovered a mistake in the AGB star model (in the code). It is relevant to the optical depth. I will fix it this week.

2. Due to the mistake, the mass loss rate was not correct. I think the mass loss rate is crucial in the formation of circumbinary disk. 10^{-7}M_{\odot}/yr is somewhat small. If there could be a 10^{-5}M_{\odot}/yr AGB wind, the formation of the circumbinary disk will be much more easier. ​Liljegren et al. (2016) did a detailed study on AGB wind recently. I will try their model.

3. I will work with Jonathan to add the self-consistent mass loss calculation in AstroBEAR.

4. Run 4 low res simulations (3au, 4au, 6au, 10au). May not finish all of them since it takes 24hrs to finish a 3au simulation.

There are two radius in this simulation with respect to the AGB star: dust sublimation radius, L1 point radius.

The dust sublimation radius is calculated by:

r_{sub}=\sqrt{\frac{\alpha L}{4\pi \sigma T_{sub}^4}}

where \alpha represent the photon absorption efficiency. L is the AGB star's luminosity. \sigma is the Stefan-Boltzmann constant. T_{sub} is the dust sublimation temperature.

A major question is, what is the value of \alpha? I choose \alpha=0.35 because some dust only absorb optical and ultraviolet wavelength photon e.g. MgSiO3. However, iron based silicates can absorb NIR photon efficiently. On the other hand, the self-extinction can redden the SED a lot. At last, a cool atmosphere should be assumed in the first place, e.g. 3000 K or so.

Usually the dust formation radii is around 2 R_{AGB} \sim 3 R_{AGB} (I saw it somewhere). If I choose the photosphere of the AGB star to be around 0.9AU then the dust formation radii can be around 1.8AU \sim 2.7AU

Let me list my chose of parameters:

\alpha=0.35,T_{sub}=1500,L=2342L_{\odot}

This will give r_{sub}=1.99\ AU

If we choose M_{AGB}=1M_{\odot},M_{com}=0.5M_{\odot} then the radii of L1 points are:

1. d_{sep}=3AU,r_{L1}=1.71AU

2. d_{sep}=4AU,r_{L1}=2.28AU

3. d_{sep}=6AU,r_{L1}=3.42AU

4. d_{sep}=10AU,r_{L1}=5.70AU

Therefore, case a is the only one that may be different from the others. But it is not necessary that case one must have different mass transfer behavior. Because the separation is still larger than the dust sublimation radii, in the current model, the radiative force on the fluid may blow away the disk or any structure around the secondary. If one want to see the disk around the secondary, the radiative ray tracing algorithm should be 3-D. Besides, the dust sublimation around the secondary should also be considered.

Difference between iron deficit and iron rich amorphous dust can be found here.

​Jena dust database

The AGB star is put into the simulation. In high resolution simulation, a circumbinary disk can form with 3AU separation

topdensity.gif​

​low luminosity AGB star

​high luminosity AGB star

​top density

​side density

There is a 1 solar mass AGB and a 0.5 solar mass secondary with 3 AU separation in the simulation. In this simulation, we include physics such as cooling, central luminosity, dust formation zone (LTE condition) and pulsating AGB star. The luminosity of the AGB star is 3500 solar luminosity.

The contour lines still have the meaning as before.

In this simulation, a gravitationally bound circumbinary disk appear and disappear from time to time. The pulsating AGB star extends its atmosphere periodically.

I think the most important parameters include:

1. luminosity/mass

2. binary mass ratio

3. separation

top.gif​

side.gif​

1. Wrote "The creation of AGB fallback shells" and published on MNRAS.

2. Went to WorkPLANS workshop in Leiden, it was a very enjoyable trip. I met some astronomers.

3. Wrote "Three dimensional hydrodynamic simulations of L2 Puppis" and submitted it to MNRAS. The paper is in its second version.

4. Pass the Ph.D qualify exam.

5. Developed the pulsating AGB model, added cooling (400K - 2000K) to AstroBEAR.

​Youtube link

2 AMR level. 0.1 AU the finest resolution.

The radial velocity at the boundary varies as v=v_{piston}\cos{2\pi t/P} where P= 1\ yr.

Dust form at fixed radii, the radius is calculated by R=\sqrt{\frac{L \xi}{4 \pi \sigma T_{df}^4}} where \xi=0.5 is the dust absorption coefficient and T_{df}=1500\ K is the dust formation temperature.

AstroBEARAGBproperties1.pdf​

AstroBEAR_single_AGB001.gif​

AstroBEAR_single_AGB002.gif​

AstroBEAR_single_AGB003.gif​

0411.gif​

constant.gif​

nonconstant.gif​

nonconstantdivvfixed.gif​

oneDstar.gif​

I intend to put this new wind model into our code. This movie shows only a radiative driven pulsating AGB wind. It has cooling.

I am going to apply it to the binary paper.

Simple model and slightly sophisticated model.

I am learning radiation transfer in some depth.

​http://arxiv.org/abs/1506.05121

This is a paper that I am also going to present in the coming journal club.

A very good literature is:

​http://deepblue.lib.umich.edu/bitstream/handle/2027.42/60735/wollaber_1.pdf?sequence=1

I think radiation transfer is a very promising field because.

1. It is under developed. There are many algorithms for optical thick and LTE (flux limited diffusion, variable Eddington tensor or M1 closure), there are also many algorithms for optical thin and non-LTE (analytic approximation or non-analytic approximation). But there is no efficient algorithm for both. Monte Carlo is an very interesting algorithm to learn and as the computer get more powerful, it will eventually become a global algorithm.

2. It is in many problems. To name a few, star formation, stellar wind, ISM, binary evolution and galaxy evolution. If we solve the radiation transfer, we can apply it to those problem easily.

3. In application, ICF and MCF both concern radiation transfer, high temperature engines will also consider radiation.

Model Description:

\chi=\frac{\dot{M}_{AGB}v_{AGB}}{\dot{M}_j v_j}

where \dot{M}_{AGB}=10^{-6}M_{\odot}/yr, \dot{M}_j, \chi and v_j are varibles.

The jet is bipolar jet and has open angle of only 10^{\circ} so the solid angle is only \omega=0.015\times 4 \pi. I calculate v_{AGB} with given \dot{M}_j, \chi and v_j.

The outflow density of the jet is \rho_j=\dot{M}_j/(\omega R_{com}^2 v_j).

The base grid is 80\times80\times160 and I added 3 levels of AMR. AGB star has radius of 4 AU and 1.4 solar mass (as in the paper by Garcia-Arredondo and Frank 2004) and the companion has 0.6 solar mass and radius of 2 AU.

This is the case of strong jet, which is,

\chi=\frac{\dot{M}_{AGB}v_{AGB}}{\dot{M_j}v_j}=0.625

where v_j=400km/s, \dot{M}_j=10^{-7} M_{\odot}/yr so v_{AGB}=25km/s which is just the escape velocity of the AGB star.

m29.gif​

Below is the weak jet case.

\chi=\frac{\dot{M}_{AGB}v_{AGB}}{\dot{M_j}v_j}=5

where v_j=200km/s, \dot{M}_j=2.5\times10^{-8} M_{\odot}/yr so v_{AGB}=25km/s unchanged.

low.gif​

Figures and movies of L2 Puppis. There are lots more in the attachment.

topview.gif​

sideview.gif​

This model is still the dense shell + constant stellar wind model.

AGB mass: 2 solar mass with alpha=0.134

Secondary mass: 0.5 solar mass

Mass loss rate in the burst is: 9.3E-6 solar mass / year last for 11.6 years

Mass loss rate of stellar wind is: 1.65E-7 solar mass / year

The speed of burst is between the escape velocity of the AGB star and the binary system. The speed of stellar wind is above the escape velocity of the binary system.

Simulation is in isothermal state at 400 K.

0.65 micron and 0.55 micron respectively. They looks identical…

SED looks more reasonable. The red line is the SED I got in current model. Note the flat spectrum energy signature from 1 micron to 4 micron. This meets the observation in Kervella et al (2015). Also, the 10 micron bump is there.

Dust composition is 0.8 % of olivine and 0.2 % of pyroxene.

​topview

​sideview

A point under debate is whether L2 pup is a binary system.

\omega=\sqrt{\frac{G M}{r^3}}

If the orbit is 2 AU, and period is half year (which is approximately the case for L2 pup), M should be 32 M_{\odot}. So it is not likely that the variation in apparent magnitude is associated with the orbital period of the binary. So the variation in brightness is not a evidence for the existence of binary system.

The pro argument is, the specific angular momentum of keplerian orbit:

l \sim \sqrt{G M r}

and for a radius of 0.5 AU AGB, this l is 10 times larger at 50 AU radii. AGB star is very slow rotator, the difference of l is huge between the surface of the AGB and the circumstellar disk. So there should be some mechanism to transfer angular momentum to the dust and gas. Binary system is a possible answer.

One may say there can be some poloidal magnetic field that extend to several AU of the AGB star, I do not know if this magnetic field is possible in a post-AGB star. Also, whether this poloidal field is effective enough to give the weakly polarized material angular momentum.

The main points in this paper may include:

1. Binary picture is suitable in L2 pup by the argument of specific angular momentum.

2. Hydrodynamically, stable circumstellar disk is to be found around the binary system.

3. Radiation images matches in V and N bands.

4. A discussion of SED. (maybe include a prediction of long wavelength result?) I am not sure what to claim in this part because there is some difference in my results and the observed results.

I use Mie model to calculate the dust optical properties as follow.

Then I put a mixture of 0.1 um olivine and 0.3 um pyroxene in the hydro simulation. The spatial distribution of the dust is a superposition of the following three distribution. Let me start from the output of hydro simulation:

1. Create a 3 AU radius sphere that has 0 dust since this region is probably too hot for dust. You can imagine we are in a spherical coordinate.

2. Make the distribution along z-axis as a gaussian and the spread increases as radius increase. Here we are in cylindrical coordinate. The functional form is C\rho_{gas}Exp(-z^{2}/s). The purpose of this layer is to create a thin disk of dust.

3. Add a dust distribution that is proportional to the gas density. This dust density at equatorial plane is low compared to the dust density in step 2 but higher in polar or other regions.

In this simulation, the total dust mass is only 1.3\times10^{-8}M_{\odot}, it is 20 times less than the dust mass in Kervella's paper.

You can compare the result with Kervella's result in ​http://www.eso.org/public/archives/releases/sciencepapers/eso1523/eso1523a.pdf at page 9.

I put a 3500 K (effective) AGB star at the center of the hydro simulation result. Its radius is 123 solar radii, equivalent to 2000 solar luminosity. This type of AGB star fits in HR diagram with ZAMS evolution. The SED is from 0.5 um to 100 um.

In Kervella's paper, they take a reddened version of AGB spectrum as its initial SED, which I did not do. If I go through that procedure, I think the SED should looks more similar.

Below are four pictures at 0.65 um, 1.24 um, 2.17 um and 4.05 um respectively, they are in the same logarithmic intensity scale.

Using Mie Model.

I am still using the default silicate data in RADMC3D.

I find a program that can calculate dust absorption and scattering coefficients for irregular shape dust (ofc regular works) from n-k data.

​http://www.ddscat.org/

q=\frac{E_{kin}}{\|\Phi\|}

q Measures the ratio of kinetic energy to gravitational potential. In stellar evolution theory, I should take enthalpy instead of just kinetic energy but since I use isothermal temperature, the enthalpy is not conserved so I did not take it.

By Virial Theorem, if h=enthalpy

h=-0.5\Phi

then the material is stable on its orbit, if it is larger than, it will go to higher orbit, if it is less than, it will go to lower orbit.

h=-\Phi

is the condition that material will go to infinity (unbounded)

So in my case, q=1 absolutely means unbounded, but gas is not bounded even q is close to 1 because I did not take temperature into account.

The figure shows that the gas in polar direction is not gravitational bounded while it is in equatorial plane. (q<0.5)

At this moment, the circumbinary disk has not reached steady state yet.

The first image shows a streamline plot of constant wind. There is no polar outflow, however, there is a circumbinary disk. Fluid parcels start in polar region just fall onto the disk. The disk looks stable as the red streamline means the fluid parcel position in later time.

The second image shows a streamline plot of dense shell wind. The disk is not stable as there is no limit circle around it. There is no red lines in the streamline indicating that the fluid parcel just run away.

Below shows a polytropic wind with gamma=1.35 . The wind start at 0.5 AU, 1500 K with 20km/s. It drops to about 300 K at 10 AU. This simulation is done with my 1D spherical Riemann Solver.

Heat can be transported by conduction, convection and radiation. High order Riemann solver can handle convection very well. Conduction is not important in the exterior region of star but important in interior (very deep, can be estimated). Radiation is the most important heat transfer mechanism in low density region, for Local Thermal Equilibrium (LTE) situation, Monte Carlo (MC) method is the best both for multi-dimension and 1D. In non-LTE situation, people sometimes use particular cooling function, for example, hydrogen cooling or any other molecular cooling, there are also MC model for non-LTE but I need to research more.

In general, MC method is universal in simulating radiative transfer.

L2pup is the a binary system that is located only 64 parsec away from earth. It has an AGB star and secondary that maybe a RG. Assume a 3000 K AGB star (blackbody) with 1140 solar luminosity (SL), we can calculate its radius to be 126 solar radius (SR) by standard procedures. The Spectrum Energy Distribution (SED) at 64 parsec away is shown in the picture below. It peaks at 1 um with 5.867\times10^{-9}W\cdot m^{-2} /\mu m intensity.

Red line is result from Radmc3d with the same condition of pure blackbody radiation in vacuum. Blue line is SED of blackbody radiation. Simulation is "Redder".

Compare to the results from other simulations by eye, the difference in magnitude is not very significant. Thus much of the light reaches earth.

Where can I get the data in Kervella's 2015 paper?

I find the accreting secondary has accreted too much gas material. In other papers talking about L2pup, they usually guestimate the mass of the dust is on order of 1E-7 solar mass with 1% dust gas ratio. In my simulation, the total mass of ejected gas is 1E-5 solar mass but the secondary accreted 50% of the gas as time goes on and another about 40% escape so the total mass in the disk is only about 1E-6 solar mass which make the simulation not realistic.

1. Read ​Winter 2002 A&A

2. Make SED of the radiative transfer simulation.

3. Nail down a proper luminosity for the AGB star. Mass loss rate is not clear now, we can treat it as a free parameter but just need to be reasonable.

4. Try more species of dust.

I used radmc3d to produce a radiative transfer simulation based on the result of my hydro simulation. The radiative transfer code use Monte Carlo method. Below are the wavelength picture correspondingly from observation and my simulation.

This is a synthetic radiative transfer simulation. It corresponds to the dense shell model​dense shell hydro simulation. I also analyzed the density and velocity properties of constant wind and dense shell model here. ​What defines a disk or a torus?

The viewing angle is the same as L2 pup, that is 83 degree. In this radiative transfer simulation, I assumed a 1 AU radius, 3000 K black body in the center of the simulation zone. I also assumed the dust/(dust+gas) mass ratio is 0.05. This is a bit high dust gas ratio but since the actual binary separation is 2 AU (in the simulation it is 4 AU), the optical depth may go higher if the separation shrink to 2 AU.

I do not know the particular species of the dust (probably silicates). I adopted the default dust species setting in radmc3d, I make the scattering to be isotropic even though the dust scattering is likely to be forward.

Further research can be done on:

1. Redo 2 AU separation simulation.

2. Figure out the dust species and their scattering, absorption properties.

3. If there is enough computational resources, I want to do a larger simulation zone. Currently, for 2 AU separation, I can only go 40AU*40AU*40AU with 100*100*100 box plus 3 level AMR and it takes 120CPU*7days.

I picked out the x-axis data (density and velocity in y-direction, which is a measurable data in Doppler shift) then I averaged over time.

The secondary locates at 2.67 AU away from origin and the primary locates at 1.33 AU away from origin. So the density spike within 2.7 AU is natural. The velocity in y direction is very spiky. I do not know why.

I do not think this is a disk or torus, people can be more clear when I compare this data to the dense shell ejection model.

Constant wind result, interpolated, averaged.

Dense shell result, interpolated, averaged.

Isothermal Temperature: 400 K

Primary mass: 1 SM with alpha=0.134 (behave like a 0.134 SM star to gas)

Secondary mass: 0.5 SM

Escape velocity to the primary star: 14.7 km/s

Separation: 4 AU

AGB wind: 8\times10^{-8}M_{\odot}/yr 24 km/s

constwind_sideview.gif​

constwind_topview.gif​

I have written a 1D hydro code using Fortran and python together. Python is used to call the excutable and fortran provide the numerically compiled exe.

This is a 1D parker wind.

My goal is to write a 1D radiative hydro code with dust formation and cooling.

I did a constant wind binary simulation last week. The escape velocity for AGB star is approximately 15 km/s. The wind is 24 km/s.

We see gas is not escaping from the system due to the deepened gravitational well. Gas material is more inclined to stay in equatorial plane and it is infalling. It does not have polar outward flow anymore but the density in polar direction is low.

constwind_sideview.gif​

I used 400 K isothermal gas temperature. In ​The dust disk and companion of the nearby AGB star L2 Puppis, they infer that the inner rim of the circumbinary disk is about 1000 K - 1200 K at 6 AU away from the AGB star. So 400 K is approximately the temperature at 24 AU - 40 AU. Since my simulation extend to 80 AU, I think 400 K is somewhat reasonable.

The gas is assumed 100% hydrogen atom. The wind is all spherical this time. The slow dense shell bear initial velocity of 16 km/s and mass loss rate of 5E-5 SM/yr while the AGB wind has initial velocity of 32 km/s and mass loss rate of 1E-7 SM/yr. The dense shell ejection last for 11.6 yrs and the AGB wind last for the rest of time so the total mass loss in this simulation is approximately 6.5E-4 SM. In the L2 Pup paper, the author estimated the dust weigh 1E-7 SM, if the gas dust mass ratio is 100, it implies that the mass of the disk is about 1E-5 SM. Usually, most of ejected material will leave the system and only a fraction stay in the circumbinary disk, I hope this fraction is 10%.

There is no dust and radiation transfer in my simulation so this initial velocity compensate some absence of the accelerate mechanism.

The AGB star is 1 SM and has radiative pressure on gas (I know the physics is much more complicate, e.g. radiative pressure on dust and dust and gas momentum coupling). The secondary is a compact accreting object which is different from the paper. I adopt the accreting object to make a hole in the center of the circumbinary disk, as I can keep removing gas from the disk. In reality, high temperature near the star can expel gas and create this hole. So temperature gradient is needed if we want to see a self consistent hole in the center.

Isothermal Temperature: 400 K

Primary mass: 1 SM with alpha=0.134 (behave like a 0.134 SM star to gas)

Secondary mass: 0.5 SM

Escape velocity to the primary star: 14.7 km/s

Separation: 4 AU

AGB wind: 10^{-7}M_{\odot}/yr 32 km/s

Dense shell: 5\times10^{-5}M_{\odot}/yr 16 km/s

Wind ejection radius: 1 AU

Dense shell ejection duration: 11.6 yrs

​topview

​sideview

compared to:

I updated some movies ​https://astrobear.pas.rochester.edu/trac/blog/zchen04132015, a new one that has longer z-axis is still computing. That one proves slow AGB wind is not the key to the formation of collimated outflow but it has similar shape to the butterfly PN in Balick's paper (Balick, 1987 ApJ). On the other hand, in this paper​The Necklace: equatorial and polar outflows from the binary central star of the new planetary nebula IPHASX J194359.5+170901, they measured that the kinematic age of polar outflow is older than the equatorial structure, implying that the collimated outflow is more likely to be the polar jets from a compact secondary. Another point to mention is that to have a circumbinary disk, a binary system just require a slow wind from AGB star, either spherical or equatorial will do the work. It does not need to carry angular momentum from the rotation of the AGB star.

I read some papers mentioned by the reviewer, I found that this paper is very very useful ​2D models for dust-driven AGB star winds. The author is an expert in AGB dust formation. In this paper, he used equation of motion and energy equation to solve the fluid motion. Dust evolution is computed by moment equations (see ​Dust formation in stellar winds. IV - Heteromolecular carbon grain formation and growth and their book ​Physics and chemistry of circumstellar dust shells of chapter 14). Woitke also included radiative transfer, he assumed Local Thermal Equilibrium (LTE) to simplify his calculation (​High precision Monte Carlo radiative transfer in dusty media). I do not understand the detailed technique and Monte Carlo simulation in the paper. In his discussion, he compared his results with 1D models and simulated the case without pulsations. The nature of pulsation is nuclear reaction lying deep in the AGB star and result in a sudden jump in luminosity. The physics in this paper is much more complicated. The temperature gradient is an issue in our paper, I think this part should be modified to a more realistic model.

I am now very interested in the dust evolution. It is actually very important to AGB evolution, binary system evolution and PN evolution. Dust formation is deeply linked to the chemical composition of the star (i.e. high metalicity, C type star or O type star). Dusts absorb and scatter 30% of phontons (see ​Interstellar Dust) so it provide the much momentum in many kinds of AGB winds. its evolution is sometimes coupled with the spectral near the AGB star, especially for the optical thick case (but rare). I think a good description in dust formation in our code will let us have more physics.

As for the binary simulation, I have an idea that we can just run a simulation to mimic ​The Necklace: equatorial and polar outflows from the binary central star of the new planetary nebula IPHASX J194359.5+170901. The modification should just be add a polar outflow with high speed from the secondary before the pulsation from the AGB star. A rough estimate is that we need 112 cores * 24 hr * 30 on bluehive (the longer z-axis one need 21 days with slow stellar wind speed.)

A recent movie​

In this specific post, I discuss the impact of separation on the formation of fall back disk around binary stars.

1. Model Description

In this simulation, there are two stars. One is assumed to be an AGB star, it has two states - spherical wind state and equatorial outflow state and the other is a companion.

The equatorial outflow state eject high density but low speed wind (20 km/s) into the space. It has open angle. Open angle means, compared to earth, the latitude. For example, a pi/6 open angle means the outflow is being ejected from pi/6 North to pi/6 South, that is, pi/3 in total angle.

In this simulation, all outflow has open angle of pi/6. The equatorial outflow last 11.2 years, the orbital period for three simulations are 4 years, 6.2 years and 11.6 years respectively, so the outflow can cover at least a full orbital period. At all other time, the AGB star is in spherical wind state. Its spherical wind speed is above the escaping velocity (40 km/s).

The initial specific angular momentum of the outflow is negligible, let's say it is approximately 0. This is in line with the property of AGB and post-AGB stars.

The temperature of the simulation is fixed at 100 K.

The separation are 3 AU, 4 AU and 6 AU (See (5,3,8) data set in ​3AU4AU6AU) respectively in the following simulations, each has topview and sideview.

2. Simulation

3AUsideview​

3AUtopview​

​4AUsideview

4AUtopview​

6AUsideview​

​6AUtopview

The equatorial plane has higher density than the polar direction. Although the AGB star has spherical outflow after the equatorial outflow, the "heavy" equatorial gas is barely supported by the successive wind while the gas in polar direction has been pushed out. However, a fraction of the polar outflow is fallen back to the equatorial plane at large radii then start to orbit the binary stars (see Muhammad Akashi and Noam Soker 2008 ​http://www.sciencedirect.com/science/article/pii/S1384107607000863). This simulation validate their assumption in 3D. We can clearly see the formation of disk in equatorial plane.

We can also see the bipolar shape of PPNe. However, the simulation zone is not big enough to see the full bipolar outflow picture.

3. Conclusion

The conclusion is that the closer the binary stars is, the more likely the fall back disk forms or the smaller the circumbianry disk will be.

Please see my further progress of this topic in this post.

I did it with MATLAB, HLLC Rieman Solver.

rho.gif​

Is it nice?

I have calculated the density variation in equatorial plane of the giant in a binary system. The parameters are:

Temperature: 3000 K

Giant mass: 1.01 SM

Companion mass: 0.5 SM

Separation: 4 AU

This plot shows the density variation of 0.2 AU, 0.4 AU and 1 AU in equatorial plane of the giant compared to the highest density at that radius.

Absolute density profile at 1 AU.

I plotted the isodensity contour of equatorial plane. It is not a circle indeed. The 3d contour has a hole because of the singularity on z-axis.

If we use AstroBEAR to exam the theory, I suggest:

Simulation zone: 10 AU * 10 AU * 5 AU

Giant boundary at 0.2 AU

Base resolution at 100*100*50

AMR level at 5 or 6

the smallest cell will be 1/32 AU = 0.03125 AU or 1/64 AU = 0.015625 AU.

Stellar rotation speed is important in this simulation. Research in the rotation period and gyrochronology has been down in M35 (S., Meibom, R., D., Mathieu and K., G., Stassun, 2009). The picture in the paper shows the period distribution of K and M type of stars.

The sun's rotation period is 24 days and the orbital period at sun's surface is 0.12 days. Therefore the picture above shows that:

\[\frac{\Omega}{\omega_{orb}}=\frac{P_{orb}}{P}\in(0.001,1)\]

and Skumanich's law is the relation between angular speed and the age of the star:

\[\Omega\sim t^{-0.5}\]

The older the star, the slower it rotate. K and M type of stars may be low-mass orange dwarfs or K type giants, but the giants can not have such short rotation period, it must be main sequence low-mass orange stars. That's the tool the paper use to determine the age of the star but it also provide us some data.

Basically, the period of AGB or giants cannot be measured at this time because the Doppler effect would be too small at its photon sphere and they indeed rotation very slow. Thus its rotation speed is negligible compared to the break up speed. Here we assume there is a third star that is so close to the primary and it spin up the envelop of the primary's.

So we need some other mechanism to give the envelop angular momentum. In N., Soker and A., Harpaz, 2000, they give a formula of the ratio of the final angular frequency and break up angular frequency for horizontal branch stars:

\[\frac{\Omega}{\omega_{Kep}}=0.1(\frac{M_p}{0.01M_{env}})(\frac{a_0}{R})^{{0.5}(\frac{\alpha}{0.1})}{-1}\]

where I_{env}=\alpha M_{env} R^2, \alpha is the coefficient that determined by the envelop's density profile.

M_p is the mass of the planet (Jupiter like or BD), M_{env} is the mass of the envelop's. a_0 is the initial separation of the planet and the star. R is the radius of the star.

An example can be:

M_{env}=0.6M_{\odot}, a_0=25R=2AU, \alpha=0.1 and M_p=0.006M_{\odot}

then the envelop will reach half of the break up speed.

I just realize that

\[\frac{\Omega}{\omega_{orb}}\]

should incorporate the effect of \alpha=0.134 for gas. (Or should I? Because observation can not determine the \alpha)

(1, 2, 9) can investigate the effect of open angle.

(2, 3, 7) can investigate the effect of angular momentum of the equatorial outflow.

(3, 4, 10) can investigate the effect of temperature.

(3, 5, 8) can investigate the effect of separation.

(5, 6) can investigate the effect of mass ratio.

Binary separation = 4 AU

Binary orbital period = 6 years

Simulation temp = 100 K

Primary:

Mass = 1 SM

Quiet phase outflow = 1E-7 SM/yr

Equatorial outflow open angle = 20 degree

Equatorial outflow = 2E-5 SM/yr

It is calculated in this way:

\[\dot{M}=\rho S v\]

\[S=2 \pi r^{2 \times 2\int_{7\pi/18}}{\pi/2} \sin\theta d\theta\]

\[v=\sqrt{v_r^{2+(\Omega r)}2}\]

Equatorial outflow density is increased by 1000 times compared to quite phase, varies with time as the picture below shows, it last for 10 years. The total mass ejected from the equatorial outflow is:

\[M=\delta t \times \dot{M}=2\times10^{{-4}M_{\odot}\]}

Equatorial outflow has angular velocity of:

\[\Omega=1.36294\times10^{-7 s}{-1}\]

This is equivalent to a period of 1.46 years. Since the equatorial outflow is launched at 1 AU, this means \Omega<\omega_{orb}.

Secondary:

Mass = 0.5 SM

The simulation:

This is a contour plot of density. The more transparent the color, the low the density. Therefore, red, which is on the top of the color bar, is the most opaque color and represent the highest density, the next highest the grey and so on.

We can see there is a primitive disk like structure around two stars. Spiral shocks compress the disk, push the disk outwards and give it angular momentum.

contour.mov​

Since the angular frequency of the equatorial outflow is below the orbital angular frequency and we still see a circumbinary disk, it means that, it the outflow is induced by the in fall on a planet.

Future works:

1. Include mass, momentum and angular momentum conservation of the sink particles so that we can study the orbital evolution.

2. Add cooling. NK cooling preferred (D., Neufeld and M., Kaufman (1993)).

Assume a third star with mass between 0.001 SM to 0.05 SM (hot Jupiter to Brown Dwarf) around the primary fall into the primary during common envelop phase evolution from an orbit of 1 AU to 0. So it will induce an equatorial outflow J. Nordhaus and E. Blackman 2006. The angular momentum, the potential energy and the kinetic energy of the planet will partially be transferred to the outflow. A detailed calculation will not be presented here.

Let's assume the equatorial outflow is 10^{-3}\sim 10^{-4}M_{\odot}/yr, lasting for 1 - 10 yrs, so the total mass of the outflow is 10^{-2}\sim 10^{-4}M_{\odot}.

Below are some parameters:

a : Binary separation, can vary from 4 AU to 20 AU. A close binary will transfer more angular momentum to emitted gas but will also challenge the small planet around the primary. Because they are basically three body system. (Boldly assume the small planet is 0.8 - 1 AU from the primary such that the secondary will impose so much effect.)

M_p : Primary mass, let's fix it to be 1 SM and low metallicity. Metallicity is important in AGB evolution since it will provide dust thus drive the wind away from the primary.

q : Ratio of secondary's mass to primary's mass. It can be 0.2 to 1, varies in correspond to the binary separation. The greater the ratio the less the orbital period thus the more angular momentum will be transferred to the gas emitted from the primary. Note here when the gas is emitted from the primary, it receives some angular momentum from the primary because the primary has orbital motion.

\theta : Open angle of the equatorial outflow. Let it vary from Pi/9 \sim Pi/6. The more mass the planet (or Brown Dwarf), the wider the open angle. This open angle is closely associated with turbulence driven by the in falling planet rather than viscosity. companion driven dynamo happens in this process.

\Omega : Spinning angular velocity of the primary. The envelop of the primary is spun up by the in fall of the third planet and may induce a dynamo. This \Omega may also come from its own spinning (though very small, data missing here). A value to be noticed is the orbital angular velocity \omega at the radii which the outflow will be emitted. \Omega>\omega.

Given that:

\[\omega=\sqrt{G M/r^3}\]

\[l_z=r\times v=r^{2 \omega = \sqrt{G M r}\]}

This means, any object with angular momentum larger than 0 can find its orbit, mathematically, if the radial velocity is 0.

T : Actual temperature can be around 2000 K, like the temperature in the envelop. In the test simulations, we can choose 200 K simply because it reduce the pressure and the sound speed.

v_w : Radial velocity of the outflow. The best scenario is we can depict its angular dependence. Here we can use step function due to low angular resolution in test runs. This speed should below escaping velocity. Let's say, ¾ of the escaping velocity.

r_0 : The radius which the outflow is emitted. It is set to be 1 AU.

A recent simulation:

Primary star has a constant alpha of 0.1339 throughout space. Simulation temp is 100 K.

density.gif​

After an equatorial outflow from the primary star.

Bounding box is 192 AU * 192 AU * 48 AU

Primary mass = 1 SM

Seconary = 0.5 SM

alpha of primary = 0.1339 throughout

Isothermal temp = 1500 K

post wind has pushed the dense shell outwards

density.gif​

Lz.gif​

contour.gif​

In order to create a 1-D spherical solver, I write a python version of Riemann solver. The long term goal is to add cooling, thermal conduction and radiation to examine theories and provide an accurate benchmark for AstroBEAR.

Simulations are run with 500 points.

Below are results from Roe solver: (each picture consist of rho, v, P)

Below are results from HLLC solver:

I am developing TVD-MUSCL scheme now, Toro's book is wrong on this Chapter. I have succeeded in writing TVD-MUSCL 2 years ago.

Binary seperation = 5 AU

Secondary mass= 0.15 SM

Secondary accrete gas around it

Isothermal temp = 1000 K

Primary's parameters:

Primary mass = 1 SM

Base density is fixed at 1.18E-5 g/cc throughout the primary (even in equatorial outflow)

Alpha of the secondary is 1. Alpha of the primary varies as below:

The velocity at the base of the wind does vary near the equatorial plane. In order to mimic the equatorial outflow induced by the in fall of planet (more massive than Jupiter), its open angle is (1/9)*Pi. Its velocity surges from 0.04 km/s to 2.4 km/s. Then equatorial outflow last only for 0.6 yr and start at frame= 150, i.e. time= 16.5 yr

Density plot + vector velocity of mid-plane:

density.gif​

Z-Angular momentum + density contour of mid-plane:

The whole angular momentum is negative, i.e. anti-clockwise. Which means, angular momentum has been transferred from stars to gas.

Lz.gif​

This simulation is basically based on Parker's model. The disk is not obvious and the pulsation is also not clear. The pulsation is subsonic and smeared out by the high density gas around the primary.

Which is supposed to be seen if the simulation is correct? Disk or spiral shock?

I took luminosity into consideration.

The primary's luminosity is large so the photons on the dust drive the dust and gas away from the primary. The secondary is assumed to be dim - although it actually maybe not given that it is a 0.5M_{\odot} WD. In all, the secondary star has full gravity on gas while the primary has no gravity on gas as soon as the dust form at 4 AU.

I also assumed that, there is a magic that the dust evolves such that the photon-dust interaction changes. The scattering and absorption of photon diminish at 12 AU away from the primary. This is only a test in order to see what will happen.

The primary is pulsating periodically with a period of 36 yrs and pulsation last for 1.1 yr. The pulsation is similar to the last post.

Temperature is isothermal at 2000 K. Separation is 7 AU. At this separation and temperature, the accretion mode is probably WRLOF or RLOF.

Total simulation bos is 320AU\times320AU\times8AU

density.mov​

It looks like a steady state disk as gas is continuously coming out from the primary.

This post is aimed to reproduce G. H. Bowen 1988

\alpha profile is fixed.

The density at the base which is at 1AU radii is also fixed as 7.52E-11 (g/cc).

Initial velocity varies in time. Its velocity in quiet phase is almost zero, the pulsed velocity is much below the escape velocity - 29.7 km/s.

The flux varies from 8\times10^{-9}M_{\odot}/yr to 5\times10^{-6}M_{\odot}/yr

flux.mov​

rho.mov​

Circumbinary Disk

Typical circumbinary disks are mainly made of dusts and accompanied with some gas (Jura & Kahane 1999). The disk is at low temperature, usually 100 K to several hundreds K (Bujarrabal et al. 2005). Our module is purely hydro, and we can not simulate dust formation. However, if the disk is purely gas, it will disperse or be reaccreted onto the stars.

The inner radii of a disk is 2-20 AU, the outer radii of a disk is 100-2000 AU. (I am looking for specific datas)

People has not found a stable circumbinary disk in simulation yet. It is possible that certain physics is not included in all the simulations (M. Akashi & N. Soker 2007).

Tidal Force

Strong tidal force will make the boundary of the primary star invalid. If we set the upper limit ratio of tidal force to gravitational force to be \gamma, the radius of the primary is r, separation is d, primary mass is m_1, secondary mass is m_2, then:

\[\frac{2 r G m_2}{d^{3}∧frac{m_2 G}{d}2}<\gamma\]

That means:

\[\frac{2r}{d}<\gamma\]

We set our r=1 AU, if let \gamma=0.1, then d>20 AU. This is an example of lower limit of binary separation.

Separation, Mass Ratio and Accretion Mode

Different binary separation and mass ratio have impact on shape of the Roche Lobe. If the wind emitted from the primary star is subsonic at L1 point, it will become wind-RLOF, if the wind is supersonic at L1 point, it is more like wind accretion picture, if the envelop its self is filling the Roche Lobe of the primary, it is RLOF.

In simulations, we want to fix the initial mass of the primary and the secondary and to vary the effective mass of the primary to adjust the Roche Lobe while keeping the orbital period the same. In the Temperature section we can see how the sonic radii behaves when we change the effective gravity coefficient - the greater the effective gravity, the farther the sonic radii. Therefore, to get a RLOF, we need a high \alpha, to get a wind accretion mode, we need a low \alpha.

Radiation Pressure

The effective gravity of the primary exerted on the gas is:

\[f_g=\frac{\alpha G m_1}{r^2}\]

Where \alpha \le 1

Prof. Morris pointed out that we can make alpha to be a function of radius, I believe this is true because in reality the dust formation and condensation will change the opacity and momentum coupling.

Another fact is that \alpha will decrease to zero during the helium flash since the luminosity of AGB star is exceeding the Eddington Luminosity but then increase slowly during the decaying phase and reach a certain value at quiet phase, thus making \alpha time dependent, too. If we have enough computational resource to simulation the whole period of a pulsation, I think we can see the envelop pushed out and then fall back.

Feed Back from the Secondary

As more material accrete onto the secondary, more gravitational energy will be released from the star. It is a good reason for us to let the high accretion rate to create a strong photon flux and offset some secondary's gravity (similar to primary's \alpha). This will make set an upper limit of the accretion and slow down the reaccretion onto the secondary.

Temperature

Wind-RLOF require the wind speed to be subsonic at L1 point, this set a lower limit on temperature. According to Parker's solution:

\[r_s=\frac{\alpha G m_1}{c_{s}^{2}=\frac{\alpha G m_1}{k T/m}\]}

Where m is the mean particle mass. At temperature between 800 K to 2000 K atoms start to condensate into molecules and dusts. So we can take m=1.3 m_{proton}.

The temperature in the photonsphere of an AGB star is about 2000 K.

On the other hand, we do not want the dispersion to be too large, we can incorporate some cooling, especially low temperature cooling.

Mass Loss Rate

Typical AGB mass loss rate is between 10^{-7}\sim10^{-5} M_{\odot}/yr. We can make the quiet phase mass loss rate to be 10^{-7} M_{\odot}/yr and He-flash mass loss to be 10^{-5} M_{\odot}/yr. The duration of He-flash is several years.

Wind Speed

Quiet phase wind speed is greater than the combined escape velocity of the binary stars. Pulsation wind speed is less than the primary's escape velocity.

Pulsation

Some reference say the typical duration of pulsation is 200 days to 2000 days (Kyung-Won Suh 2014), which is more friendly for simulation. So, I assume the pulsation duration to be 1000 days at this moment and decays one E fold in 10000 days. The density during pulsation is 500 times higher than the quiet phase, so after 50000 days, the density will drop to quiet phase density.

From my experience, we can choose following parameters:

A typical simulation box for 20 AU separation binary will be 100 AU\times 100 AU \times 20 AU which will be divided by 200 \times 200 \times 40 base cells. Therefore a base cell has physical dimension of 0.5 AU \times 0.5 AU \times 0.5 AU. Using 2 levels of AMR will make the finest grid be {(0.125 AU)}^3, thus the primary star will have a radii of 8 cells. This is actually 8 times larger than my previous simulations.

If the simulation time range for this simulation is 1000 yrs. I guess it could take 200*120 CPU hrs or more.

Computational temperature unit is 2000 K. The gas emitted from primary is 2000 K. Mass of secondary is 0.5 SM, mass of primary is 0.8 SM but has radiation power. Separation is 7 AU. gamma = 1.01, which means very isothermal, but not exact.

Max temperature is about 6700 K, minimum is about 1760 K.

This sim region is 24 AU by 24 AU, if we extend to 100 AU by 100 AU, the minimum temp will drop below 1000 K probably.

temp.gif​

Some references that may be referred in the fall back shell paper.

​http://adsabs.harvard.edu/abs/2000A%26A...353..583O

​http://adsabs.harvard.edu/abs/2002ARA%26A..40..439B

​http://adsabs.harvard.edu/abs/2014apn6.confE..59M

​http://adsabs.harvard.edu/abs/2006MNRAS.370.2004N

​http://www.annualreviews.org/doi/abs/10.1146/annurev.astro.43.072103.150600

​http://adsabs.harvard.edu/abs/1996A%26A...305L..61G

​http://adsabs.harvard.edu/abs/2010A%26A...515A..27O

obliqueview1.gif​

Binary separation = 7 AU

obliqueview2.gif​

Binary separtion = 20 AU

oblique.gif​

Separation: 7AU

Temperature: 2000K

Primary mass = 0.8 SM

Secondary mass = 0.5 SM

This term paper relate to our research. So I would like to present it in group meeting.

Click the picture to see enlarged picture. This is a simple incompressible or steady state flow, mass conserved disk-jet model. They are stream lines.

Initial separation = 4 AU

Temperature = 1500 K

Primary mass = 1.5 solar mass but behave like 0.1379 solar mass star when interact with gas

Secondary mass = 0.5 solar mass

Mass conserved, momentum conserved.

By the way, at this resolution, I did not see the disk tilting. Maybe tilted disk is a resolution issue or illusion?

alpha = 0.12 case

temperature = 2000 K

Binary separation = 7 AU

Temperature = 2000 K

Primary mass = 0.8 AU with alpha = 0.1723 that is, primary star behave like an 0.1379 solar mass star. This is presumably due to radiation pressure on dust grains (so the temperature should be like 1500 K)

I am not sure about the mass loss rate of the primary star.

Secondary accrete gas with initial mass of 0.5 solar mass

First picture shows the mass gain of the secondary and the separation decrease of the two stars

​midplane

This gif file is the mid plane density plot of the whole simulation.

People should not take above result too seriously. Many things, for example, mass conserved primary star, self-consistent boundary condition of the primary star, tidal effect and tidal force, non-isothermal situation, should be taken into consideration.

The following are simulations with mass conserved primary.

Red line is the theoretical angular momentum computed from the instantaneous separation (that is, the blue line) of the two stars and their corresponding mass. It is calculated by:

L=m_1 m_2 \sqrt{\frac{G d}{m_1+m_2}}

Green line is the result of angular momentum from AstroBEAR.

Initial separation is 4 AU, primary mass = 0.8 solar mass with alpha = 0.1723, secondary is 0.5 solar mass. Gas temperature is 1500 K.

The orbit changed to elliptical one since I turn on accretion and mass conservation mechanism suddenly, so it serves as a perturbation. Distance, theoretical angular momentum and computational angular momentum synchronizes afterward. We can see that the mean slope of green line is lower than the red line. That is to say, part of the orbital energy is changed to some other kind of energy. Since the simulation is isothermal, it can only be transformed into gas's kinetic energy or its gravitational potential energy.

I forgot to deduct momentum of the primary in simulation, that is, its momentum is not conserved, that is why angular momentum is increasing.

Fig1 shows the roche-lobe and 5 Lagrangian points. The red ring is the boundary of the primary star. Gif2, gif3 and gif4 are animations. Fig5 is the mass gain of the secondary.

I use thermal-driven wind model, which is parker's solution of isothermal wind. I have checked the wind speed for single star case, it is about 30% higher than theoretical solution at this resolution. So the mass loss rate will not deviate too much.

The mass loss rate of a single primary star (without the secondary) is 2\times10^{-9} solar mass per year. The mass of the primary star is 0.1379 solar mass, secondary is 0.1 solar mass. Isothermal temperature is 2000K, sonic radius is at 7.5 AU, also the wind should reach its escape velocity at about 9 AU according to Parker's solution. The separation of the two stars is 4 AU. The scale in Roche-Lobe picture is 1AU per unit length, that is, red ring has radius of 1 AU.

But when I add the secondary there, the accretion rate is 1\times10^{-8} solar mass per year. I think this is reasonable because the secondary is in the subsonic region and its existence is affecting the structure of the wind, or the stellar structure of the primary, which leads to much higher mass loss rate.​ The reason that the mass loss can behave very differently to the solely thermal-driven case may be the secondary is attracting the gas therefore removing the gas from the boundary much faster than thermal-driven wind. Then the more gas will be pushed from the boundary due to pressure (Ideal gas).

rochelobe.pdf​

45degreeview.gif​

equatorialview.gif​

z0.gif​

massgain.pdf​

3d.gif​

Wind accretion. Boundary set to be outflow only.

Accretion has been turned on.

Feedback from gas to particle is turned off.

tracer.gif​

2.5D breakup.

3D boundary wind object is not good.

z0.gif​

3D Binary simulation can not restart.

z0.gif​

test020.pdf​

escape3_6.pdf​

subplot.pdf​

Below is the numerical result from AstroBEAR.

During time =0.20 and time=0.25, a pulse of wind is launched. Density and velocity vary with time.

3_6.gif​

EVERYTHING IS ABOUT "AS IF"

If we locate the origin at the primary star. The Newton's law is:

F=\frac{G m M}{r^2}=m_{reduce} \omega^2 r

\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho \mathbf{v})=0

\frac{\partial\rho \mathbf{v}}{\partial t}+\nabla\cdot(\rho \mathbf{v} \mathbf{v})=\mathbf{f_1}+\mathbf{f_2}

\mathbf{f_2} is using m instead of m_{reduce}

r=r_1+r_2

r_1 is the distance from primary to COM

r_2 is the distance from secondary to COM

If we locate the origin at COM.

F=\frac{G m M}{r^2}=m \omega^2 r_2=m_{reduce} \omega^2 r

\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho \mathbf{v'})=0

\frac{\partial\rho \mathbf{v'}}{\partial t}+\nabla\cdot(\rho \mathbf{v'} \mathbf{v'})=\mathbf{f_1'}+\mathbf{f_2'}

So the newton's law does not really change if we just use F=\frac{G m M}{r^2}

for the fluid mechanics part, it explicitly has the following relation:

\mathbf{v}=\mathbf{v'}+\mathbf{v}_{com}

and \mathbf{v}_{com} is known - doing circular motion

Therefore, it is the mapping from COM frame to primary frame that is complicated. As long as we always work in the same frame, we are fine.

\frac{d\mathbf{r}}{dt}=\frac{d\mathbf{r'}}{dt}+\frac{d\mathbf{r}_{com}}{dt}

NOTE: All quantities with com subscript are known.

The trick here is to use m_{reduce} in kinematic equation but use m in fluid equations and newton's law. Thing will get more complicated if secondary accrete gas - need to update reduced mass every time step.

Actually you can view every parcel of gas as a body, and it will has its own reduced mass with respect to either star which is approximately its own mass. That justify why we don't need to change the density. If we have to change the density, serious things will happen - EOS changed and that will affect the dynamical picture.

Two simulations have the same following parameters:

Temperature=20K

Density varies as:

Same wind launching radius (1 AU) and resolution (320*320 base grid of one quarter plane with 3 level of AMR)

Wavenumber: 40 waves of perturbation in 2\pi radian

Density of Perturbation: Both are 2.5% of the peak the density.

Two simulations only have tiny difference in velocity:

[[Image(lowvelocity0200.png,width=400)]][[Image(highvelocity0200.png,width=400)]]

lowvelocity.gif​highvelocity.gif​

3D binary star:

Mass loss rate varies from 10^-7 to 10^-5 solar mass per year. An AGB star locates at the origin and the separation of two stars are 7.5 AU.

3Dcolor.gif​

3Dsideview.2.gif​

3Dtopview.gif​

From the side view, we can see BH accretion of the secondary star during the pulsation of the AGB star.

I think the tilting disk has gone. Since it is low temperature simulation and during the evolution we saw a flat plane of high density gas.

The simulation stuck due to "Restart Required"

​2.5D

Escape velocity is 2.64 CU, and 2 CU velocity wind will halt at 7.34597 AU.

1 CU length = 10 AU. 1 CU velocity = 3.9 km/s

Left side: lowest velocity is 1.65 CU

Right side: lowest velocity is 1.89 CU

color.gif​color2.gif​

test.gif​test2.gif​

It is about debuging. The former program let the sink particle move when it interact with gas via gravity(not via accretion), the trajectory of the sink particle is odd.

The result of this is that asymmetric angular momentum of the outflow will be generated. Theoretically, there should not be any angular momentum at all, but I guess it is the boundary effect on xy, xz, yz planes.

jzasy.gif​

Then I fixed the sink particle and do not let it move. Using outflow parameters:

\rho-t

Two different outflow velocity

Jz become symmetric but its value increases.

jzsymmnonesc.gif​jzsymm.gif​

Their line plots are:

test1.gif​test2.gif​

Red line indicate the escape velocity at the boundary.

Green line is the escape velocity at that radius.

test.gif​

color.gif​

The program still has bug in sink particle.

Green line indicates the local escape velocity. Raleigh-Taylor Instability happens at last stage. Temperature is 100K so that wind are all super sonic (low as 8 Mach speed high as 28 Mach speed).

From theoretical calculation, the dense shell should stop at r=7.4AU, in sim, the shell stop around r=6.2AU

test1.gif​

rho.gif​

I lowered the temperature of the wind by a factor of 20, now it is only 100K - freezing cold. So every wind is now super super sonic wind now. The dense shell stopped at r=10AU, from theoretical calculation, it should stop at r=8AU.

rho.gif​

Logrho.gif​

v.gif​

The program stopped there with:

High CFL - restarting step. 0.22934E+03 > 0.10000E+03

forrtl: severe (174): SIGSEGV, segmentation fault occurred

Wind parameter as time dependent variables

If there is no pressure, the second wind should halt at r=4.7AU.

Logrho.gif​

v.gif​

The velocity=0 happens at r=4AU, which seems meet the theoretical result.

rhov.gif​

Very low res, star radius only consist of 8 cells, but the result is good. The zigzags are due to boundary. Red is from Parker's solution.

X axis is the radius, 1 CU = 2.5 AU

The outburst is very obscure. There is no significant shock in linear plot.

rho.gif​

rescale2.gif​

vfield.gif​

rho.gif​

The plot of isosurface of density. The disk is a little tilted.

isotopview.gif​

isosideview.gif​

Start with single star. The mass loss rate vary wish time as (increase to 10 times at time=1.5 and reduce to 1/10 at time=3.0):

and I fix the velocity at certain radius so that mass loss rate is proportional to wind density. The velocity at exact radius is obtained by solving:

(\frac{v^2}{a^2}-1)-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4(\frac{r_c}{r}-1)

numerically. From isothermal condition(2000 K throughout space) we can get the velocity that certain radius.

Under these conditions, the single star simulation become:

(Picture on the left is rho + vector v, picture on the right is speed + vector v)

single.gif​

Pick out the magnitude of velocity on x axis:

lineout.gif​

The relative difference of final velocity to the solution from Parker's stellar wind model is around 5%.

Problem 1: It is not symmetric.

Problem 2: When the star has relatively low mass loss rate, the wind speed profile does not meet the Parker's solution and the stellar structure behaves like a bubble. Either there is some threshold for this boundary condition or something worse is out there.

Problem 3: The speed profile plot of x axis does not match the color plot at small radius.

Phenomenon: It already has some fall back behavior.

Then we put a secondary star near the primary star and let they do circular motion. Their separation is 3.999 AU. The mass of them are 1M_{\bigodot} and 0.5M_{\bigodot} respectively.

zoomout.gif​

zoomin.gif​

Problem: The CFL number blow up.

High CFL - restarting step. 0.43877E+07 > 0.10000E+03

forrtl: severe (174): SIGSEGV, segmentation fault occurred

When there is only the primary star, this does not occur, but when I add the secondary, it does.

The fall back phenomenon is not obvious at this time frame.

This is a zoomed in simulation. There is a disk around the secondary star. I did not make the point particle to actually accrete material. Please pay attention to Time>3 when mass loss rate drop again, it seems there is a fall back process.

Binary.gif​

But I still encounter this (thus the program stuck here at Time=4):

High CFL - restarting step. 0.33520E+03 > 0.10000E+02

forrtl: severe (174): SIGSEGV, segmentation fault occurred

Is my physical resolution too high?

Mass loss rate vary with time as:

while v is constant thus density changes like Mass loss rate with time.

This is a test run. Constant mass loss rate at 10E-4 solar mass per year, the primary star's mass is 1 solar mass and the secondary is 0.5 solar mass. Separation is 3.999 AU, computation box is 32 AU * 32 AU * 16 AU(z dimension). Gas temperature is 2000K. According to Parker's solar wind solution, the sonic radius of the primary star is 2.5 AU. I specify the wind velocity at 0.5 AU which is 0.23v_{sound} thus can get the density at that radius from constant mass loss rate.

circular.gif​

Mass loss rate is 10E-4 solar mass per year, separation is 4 AU but it seems that the orbit is elliptical?

Binary.gif​

Binary simulation encounter segmentation fault. (in the attachment) If I don't let the stars move, CFL number will be OK and will not trigger restart and the program is OK. But when I add velocity to stars, it encounter segmentation fault.

Stellar Winds shows that if we fix mass loss rate while changing wind density and wind velocity, the wind will not excite shock wave. Only when mass loss rate is increased will the successive winds generate a shock. (I had some problem on bluehive2 and can not make the movie, will ask for help, but this is the conclusion.)

I was testing the program on bluehive2.

This is the wind vprofile when it reaches steady state. However, it is very weird that when I use color plot, it shows that the lowest velocity is 0.005691 which is exactly the boundary condition I specified. Besides that, a considerable region is below unity, but in the lineout plot, the lowest velocity is above unity.

The boundary condition would agree with the color plot, but why there is a difference?

The asymptotic behavior of the wind is very close to Parker's solution (

comparing the final velocity, the relative diff is only 5%

), the global quantitative difference has not been determined yet since there maybe a mistake in the lineout plot.

They are not very symmetric…

(\frac{v^2}{a^2}-1)-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4(\frac{r_c}{r}-1)

is solved by Newton-Raphson method.

This is the velocity profile

I partly resolve the firewall issue by using VPN, but the connectivity is still poor.

I read Chapter 4,5 of Astrophysics in a nutshell and get some basic idea and questions about star formation, PN and PPN. It is interesting that the birth and death of stars in a molecular cloud is connected by turbulence, magnetic field and radiation.

I did not manage to see the Dean of Astronomy of Beijing Normal University because I felt uncomfortable in Beijing after staying for 5 days and went back to Wuhan.

Below is the research in stellar wind:

The next goal is to change the intensity of radiation force. Restate the governing equation first. To change the radiation force means to change the sonic radius and perhaps the temperature.

(\frac{v^2}{a^2}-1)-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4(\frac{r_c}{r}-1)

The first way is to solve this algebraic equation by Newton-Raphson method (This is one of the many methods). A given r will correspond to 2 v values as solutions. Then plug the r vs. v, r vs. \rho profile (from the wind base up to sonic point) back to outflow source and then change the wind.

Another way is to solve the hydrostatic equation:

\nabla P=-\frac{GM}{r^2}

M is constant mass because I assume the wind base is at moderate radius s.t. the inner mass dominate. Use isothermal assumption:

P=\frac{2\rho k T}{m_p}

The gas is assumed to be consisted of fully ionized hydrogen. m_p is the proton mass and \rho is the gas density.

The final solution of \rho is:

\rho=\rho_0 Exp(\frac{m_p G M}{2 k T}(\frac{1}{r}-\frac{1}{r_0}))

where \rho_0=\rho(r_0)

Hydrostatic solution is appropriate only when \nabla P \gg \rho u \cdot \nabla u i.e. in much inner region than the sonic radius.

The strategy is to specify the wind base v and \rho profile and let AstroBEAR build self consistent atmosphere and wind via Riemann solver.

Remark: People may think the first way is easier and more computationally economic, but the second way also has certain advantage: it is more self consistent, we can further specify the nonuniform magnetic field at the wind base. If we go deeper into the AGB star I think we will touch the regime of dynamo theory and polytrope.

China has firewall against USA. I can not login into bluehive or alfalfa as smooth as in USA. I have downloaded some astrobear f90 files to read, I am reading Clumps and PlanetaryAtmospheres to learn how to generate density, velocity profiles.

I am also reading a book: Astrophysics in a nutshell. I have read 3 chapters. It is interesting and helpful in understanding the process of radiation and binary evolution (as in chapter 4)

here is its amazon link.

​http://www.amazon.com/Astrophysics-Nutshell-In-Princeton/dp/0691125848/ref=sr_1_1?ie=UTF8&qid=1401716251&sr=8-1&keywords=astrophysics+in+a+nutshell

I am now in Tsinghua University (meeting friends and learning here) and will go to Beijing Normal University next week to meet Zonghong Zhu. He is the dean of Astronomy of Beijing Normal University.

​http://astrowww.bnu.edu.cn/CN/index.php/2010-04-28-10-49-23/122-zzh

I installed visit on my laptop and tuned to be able to use visit on alfalfa remotely.

Plan to modify the outflow module. Since as the radiation force changes, velocity and density profile changes dramatically within the sonic radius and sonic radius also moves. So I probably need to solve this equation and plug the velocity and density profile into outflow module.

(\frac{v^2}{a^2}-1)-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4(\frac{r_c}{r}-1)

wind.swf​

Starting from equation:

(v^2-a^2)\frac{1}{v}\frac{dv}{dr}=\frac{2a^2}{r}-\frac{A}{r^2}

It is:

\frac{dv^2}{2dr}-a^2\frac{dlnv}{dr}=\frac{2a^2}{r}-\frac{A}{r^2}

Integrate with respect to r:

v^2-2a^2lnv=2a^2lnr+\frac{A}{r}+C (1)

Where the sonic point implies:

A=2a^2 r_c

Substitute above relation in to equation (1) we can get:

\frac{v^2}{a^2}-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4\frac{r_c}{r}+C

For wind start subsonic and accelerate to supersonic, apply boundary condition that:

v=a,r=r_c

Then C can be determined to be -3

(\frac{v^2}{a^2}-1)-ln{\frac{v^2}{a^2}}=4ln{\frac{r}{r_c}}+4(\frac{r_c}{r}-1) (2)

Given a certain radius and time dependent A(t), we can find the velocity at the radius by solving (2)

In the limit that r\rightarrow\infty,v^2\gg a^2

v\frac{dv}{dr}=\frac{2a^2}{r}

v^2=4a^2ln\frac{r}{r_c}+C

v will grow indefinitely but very slowly. This implies that if wind is to fall back, right hand side should change sign, or it will just keep growing. So, gravity (and radiation driving force) together should give a non-constantly accelerating wind solution—it should decelerate at some radius beyond the sonic radius.

Circumstances might be: Optically thick AGB? Or binary star located at radius larger than sonic point s.t. deepen the potential well significantly.

Choose the star mass = 1 solar mass. Mass loss rate 1E-5 solar mass per year. Radiation force to be 92.83% of gravitational force and applied to all space. Temperature is 2000 Kelvin. Isothermal wind.

Result from Mathematica is shown as below:

Below is the result from AstroBEAR:

The relative difference is about 16%. The reason we have this difference is that we have to specify the wind speed and density on the edge of the star and AstroBEAR module will linearize them within the star and apply them as ghost cell in Riemann Solver. Thus create self-inconsistent boundary conditions.

Another factor is the resolution, which is obvious.

If we want to get a more realistic result:

Solve the problem to give the boundary condition. However, the problem may become tricky because the solution is singular at sonic point.

1 solar mass, mass loss rate 1e-5 solar mass per year, mean molecular weight 1.3 proton mass, cs=3.1km/s, T=1500K isothermal, dust form around 2.4AU, sonic point at 2.5AU, radiation force off around 25AU(should have turned off gradually and smoothly, here is for simplicity). radiation force strength around 95% of gravity, constant.

Red: Wind velocity

Blue: Escape velocity

Actually wind velocity may be just above the escape velocity. If we make it a binary system, the deepened potential well will trap material and have some phenomenon like fall back disk, accretion disk around secondary.

Simulation requirement: For binary star, 10AU separation use 40 cells to resolve, therefore, finest cell correspond to 0.25AU. 200AU*200AU*100AU simulation box. This can define the roughest cell, like 2AU.

Base cell: 100*100*50 Level: 3

Star usually need extra resolution to see the disk. Then use 4 level AMR.

This post is mainly intended to support my contemplation.

Theoretically, winds can't be isothermal throughout the space. It will cool and condense to other matters. Dusts ramp into ISM and gas and and finally drag force win radiation force (hypothesis) then lost its momentum quickly. Then driven force on gas vanishes.

Blue: Isothermal continuous dust driven wind. Driven force applied to everywhere.

Yellow: Isothermal continuous dust driven wind with driven force turned on at certain radius (0.8 numerically) and extend to infinity.

Green: Non-isothermal continuous dust driven wind with drive force turned on at 0.8 and extend to infinity.

Red: Non-isothermal continuous dust driven wind with drive force turned on at 0.8 and turned off at r=5.

Blue and yellow line coincide at large radius. Here is the enlarged [0,5] radius.

Potential I used: driven force turned on at 0.8 and off at 5.

Temperature profile I used: for non-isothermal case.

Conclusion: Wind material can fall back as shown by red line. Its velocity will converge to 0 (solution blow up at this point) and physically start to fall back. We can guess it will fall back to the outer potential well and stay there. To push one step further, materials will condense there and forming rocks, asteroids and so on?

Red: Radiation driven wind

Blue: Without driving force

Singularity of isothermal, radiation-dust driven wind

\frac{1}{v}\frac{dv}{dr}=\frac{\frac{2a^2}{r}-\frac{GM}{r^2}+\frac{A}{r^2}}{v^2-a^2}

Nature of dust driven wind - momentum coupling

Gas momentum equation

\rho_g(\frac{\partial u_g}{\partial t}+u_g\cdot\nabla u_g)=-\nabla p_g+\rho_g g+f_{drag}

Dust momentum equation

\rho_d(\frac{\partial u_d}{\partial t}+u_d\cdot\nabla u_d)=-\nabla p_d+\rho_g g-f_{drag}+f_{rad}

f_{drag}=\pi b^2\rho_g w_{dr}\sqrt{w_{dr}^2+a_g^2}

f_{drag}: Drag force on gas

w_{dr}=u_d-u_g: Drift velocity

b: Radius of dust grain

a_g: Sound speed of gas

We can argue that as dust grain moving through gas, drag force can "immediately" let dust velocity relax to a proper drift velocity with respect to gas. Therefore it is approximately quasi-static everywhere and drag force should be balanced by radiation force and any other forces. Usually radiation force and drag force can be balanced.

f_{drag}=f_{rad}=\frac{\pi b^2 Q_{rp} L_*}{4\pi r^2 c}

Where Q_{rp} is radiation pressure mean efficiency, it consist of Planck mean absorption efficiency and scattering mean efficiency and is dependent on the properties of light and dust grain. Such as radius and mass of dust grain, constituents of dust grain, wavelength of light and light spectroscopy. Further, formation of dust grain is relevant to the stage of the star.

L_* luminosity of the star

At last, equate drag force and radiation force will give:

w_{dr}=\sqrt{\frac{Q_{rp} L_*}{4\pi r^2 \rho_g c}}=\sqrt{\frac{Q_{rp} L_* u_g}{\dot{M}c}}

I learnt some wind theory and AGB knowledge. I also did some modification to the program, but the program has some weird problems.

​http://astrobear.pas.rochester.edu/trac/astrobear/wiki/u/zchen/2

​Feb 17th

Wind density \propto e^{-t}

Wind speed \propto e^{-t^2}

Gravity \propto Arctan(t)

These are very zoomed in simulations. Actual simulation box is far beyond these movies.

1. Density and velocity vector plot. There is fall back as well as outflow process at the same time.

rhoz0.gif​

2. Velocity magnitude and velocity vector plot. The stationary region (characterized with blue color) is rotating. As the time goes on, stationary region also emerges in outer spiral.

The nature of this rotating stationary region is due to retardation.

vz0.gif​

Our next intention is to consider the radiation effect of the AGB (primary) star in forming the fall back disk.

Roughly speaking, this will modify the hydrodynamical equations:

\[\frac{\partial \rho\mathbf{u}}{\partial t}+\nabla\cdot(\rho\mathbf{u}\mathbf{u})=-\nabla P-\nabla E-\nabla \Psi\]

Where E(\mathbf{r},t) stands for radiation energy density. Radiation driven wind is actually radiation-driven dusty wind. A portion of Radiation energy is transferred to dust's kinetic energy then transferred to gas kinetic or thermal energy. According to Poynting Flux theorem, we know that radiation energy density drop by \frac{1}{r^2}. Neglect the inhomogeneity and non-linearity in energy transfer between radiation and dust, dust and wind, we can combine the radiation effect with gravity.

Assume the total radiation power of the AGB will fall in binary evolution. I use the following prescribed time dependent curve to illustrate this problem.

1. Gravity \propto Arctan(t)
2. Wind speed, wind density \propto e^{-t} or e^{-t^2}

The goal is to explore the parameters to get fall back phenomenon. Unfortunately, after two weeks of trial, the result is still not good.

Test A: gas is gravitationally bounded, no falling back process.

z0bounded.gif​

Test B: gas run away supersonically, no falling back.

The AGB seems to blow out all its envelop at some point, which looks as if I suddenly turned off the wind, but I did not.

z0escape.gif​

Skimmed through:

1.A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids​link

2.Numerical Simulations of Radiatively-Driven Dusty Winds​link

I attended MS graduation test on Friday last week (tested by Professor John Thomas and Professor Chuang Ren). It was good.

I attended NSF workshop today.

Read following 3 papers:

1. "Collimated Outflow Formation Via Binary Stars: Three-Dimensional Simulations of Asymptotic Giant Branch Wind and Disk Wind Interaction" by F. Garcia-Arredondo and Adam Frank ​link

2. "Bipolar Preplanetary Nebular: Hydrodyamics of Dusty Winds in Binary Systems. I. Formation of Accretion Disks" by Nikos Mastrodemos and Mark Morris ​link

3. "Bipolar Preplanetary Nebular: Hydrodyamics of Dusty Winds in Binary Systems. II. Morphology of Circumstellar Envelops" by Nikos Mastrodemos and Mark Morris ​link

They are all very interesting. I am inspired by these three papers.

This is an experiment. A drop of mercury is vibrating in water. It resemble the fractal of accreting sink particle very much.

Actually they have physical resemblance, too. In this experiment, the vibration of mercury and water are coupled. In our simulation, high density part and low density part are also coupled. Waves serve as the driving force in the experiment.

vibratingmercuryinwater.flv​(you need to download it to watch)

It may imply that our simulation is physically right!

1. In fact binary star, but the secondary is 400AU away from the primary and its mass is only 1e-9 solar mass (Pluto), so its impact is negligible.

Primary star use Krumholz accretion criterion. I think rectangular geometry induce the asymmetry and the accreted gas has similar behaviour to the no accretion case.

Why the Krumholz accretion will give more symmetric simulation?

In Lagrangian picture (lab frame):

\rho\frac{D\mathbf{u}}{Dt}=-\nabla(p+\phi)

p=nkT

k is Boltzmann constant, T is isothermal temperature.

therefore:

\rho\frac{D\mathbf{u}}{Dt}=-\nabla(nkT+\phi)

Since gas is not in Keplerian motion - pressure balance a fraction of gravity.

The effective gravity is smaller than the case without isothermal assumption. Gas can stay in higher orbit if there is pressure gradient. In isothermal simulation without accretion, it is really the case.

Krumholz accretion remove gas, therefore remove pressure. Gas will suffer much more "gravity" and fall back - finally ate by the sink particle.

So even if there is some asymmetry, i.e. angular momentum in this simulation, we can not see because gas is eaten. But in no accretion case, the unwanted "angular momentum" is amplified by pressure.

In all, simulation become more insensitive to angular momentum if there is accretion. That is why Krumholz accretion seems to make the fall back symmetric. Actually it is not convincing.

z0.gif​

y0.gif​

2. Plot of isosurface of density around sink particle(half of the MAX density). The data of this movie is from the same simulation of the above one.

iso.gif​

3. Wind density decays exponentially at first. After 1 e folding time, I shut off the wind(at cycle 150, time=50). There is still a simulation that does not shut off the wind running on bluehive. Wind is ejected from the primary at escape velocity from 10 AU.

z0binary.gif​

1. 40AU separation, 1.5 solar mass + 1.0 solar mass. There is no accretion for both stars. I did not assume the primary star has radiative pressure that balance out its gravity in the outflow phase, that is, the primary impose gravity on gas throughout the simulation. outflow (wind) temp=2000K, sound speed is roughly 4km/s, wind velocity=8km/s therefore Mach=2, supersonic flow. Mass loss is supposed to be 10e-5 solar mass/yr but I doubt that because the radius of primary seems larger than 10AU which is supposed to be. But in the code, its initial setting is 10AU indeed. Wind duration is 1200yr, whole simulation time is 2400yr. Simulation box dimension is 1200AU*1200AU*600AU.

Comment: The fall back disk is not so obvious when looking from the side, I think it may be due to the high temperature (2000K) of the outflow.

z=0 plane view

z0.gif​

Side view

y0.gif​

Enlarged view.

enlarge.gif​

2. single star

z0single.gif​

enlargesingle.gif​

1. Separation = 40AU, 1.5 solar mass primary + 1 solar mass secondary. wind initial velocity is 1 c.u. = 4km/s (below escaping velocity). wind density is 2.64e-5 c.u. emitted from 10 AU, mass loss rate 10e-5 solar mass per year. Non-accretion for both stars.

The result is very complicated. It seems the gas is on the transition to fall back or the wave pattern will continue until very far away.

I doubted that the unbinding effect win at this separation, so I did another sim.

z0.gif​

This is the side view of the same sim.

y0.gif​

2. Separation = 400AU, others unchanged.

Instead of falling back directly, there is acceleration region around the primary star which is not due to unbinding process of the secondary. Which is so weird. And the gas decelerated very fast in a small region, but I don't see significant fall back process.

Both star are not accreting gas.

z0large.gif​

1. 40 AU separation, 1.5 solar mass primary, 1 solar mass secondary, wind emitted from 10 AU with velocity of 1km/s, which is far below escape velocity. Mass loss rate is 10E-4 solar mass per year. Something like CE or fall back disk formed and there is a numerical shock. The secondary is actually unbinding the envelope.

CE.gif​

2. 40 AU separation, 1.5 solar mass primary, 1 solar mass secondary, wind emitting from 10 AU with velocity of 5km/s, also below escape velocity. Mass loss is 10E-2 solar mass per year, which is rather fast. A steady disk formed around secondary. We can see the boundary of primary impose some non-physical effect and induce inconsistent behaviour.

disk.gif​

I don't quite understand the meaning of t1,t2 and radiusw in the binary module. Does the radiusw refer to where the wind is emitted in the unit of base cell? That is, if the sim box is 200 AU and I have 256 base grid, then each cell has 0.8 AU and if radiusw=12 then the wind emit at 10AU?

This is a very long time sim. It was terminated around 20 cu. It shows a bug in program because bluehive indicated an unexpected termination which has things to do with pointers.

It evolved to the picture just before the common envelop stage.

infall.gif​

With Baowei's help, I finally made the Binary module run. This is a basic simulation. There is no fall back steady disk but BH accretion. I will adjust the wind velocity and prolong the simulation time to see what happens.

I have also run the basic disk and binary modules. The binary module was not working while the disk is good.

Was reading the objects and modules in astrobear since there is no module simulate the fall back disk I may have to write a new one. The objects and modules are complicated. I think the module that has the most resemblance to the fall back disk is corotating-binary and basic-disk, but neither of them simulate the global picture of binary star.

Basically reading book - the Riemann solver book and astrophysical fluid dynamics book. The 2D programming has some problems.

I thought about the problem of collapse of molecular cloud and hope that Virial theorem may introduce something new by replacing the energy equation in Euler equation. But found nothing significance so far.

I took a rest last week because I was not in good condition. I was still thinking about the 2D problem. I will finish the programming this week.

Read chapter 16 of Toro

Only read chapter 15 of Riemann solver. I apologies about that. Quantum mechanics and E&M theory took a lot of my time. They have too much homework. I am trying to balance the time (or probably exploring new time). As for the 2D problem - the supersonic flow pass through a flat plate - this problem is meant to solve full NS equation and EOS. The difficulty is I don't have a Riemann problem solver to solver NS equations. All previous practice are just the implementation of Riemann solver to Euler equations.

1. Summarized my method which was implemented in Second order Conservative MUSCL-Hancock Method with Total Variation Diminishing principle and HLLC Riemann Solver code. It is in my notebook.
2. Reading this paper​Primitive, Conservative and Adaptive Schemes for Hyperbolic Conservation Laws

wrote a program using MUSCL-Hancock Method with TVD principle (MINBEE slope limiter). Details can be found here ​MHMTVD

1. Updated ​mypage
2. Read chapter 13 and 14 in Toro's book. They make sense to me in general, but the detail, reason, motivation is not clearly understood. Anyway, I know how to implement MUSCL-Hancock Scheme and will write its program to solve Euler equations this week. My interest in these two chapters is about the functional analysis in deriving Total Variation Diminishing (TVD) theorems, in overcoming Gibbs phenomenon (happens in the presence of high gradient of data), in devising specific schemes. Still, I also hope that Fourier analysis and Spectral theory can be implemented in Riemann solver sufficiently.

1. I updated my page. Solution using HLLC scheme is here​https://astrobear.pas.rochester.edu/trac/astrobear/wiki/u/zchen/HLLC
2. I have written Roe scheme Riemann solver with entropy fix. I will post its solution on my page tonight.

1. HLLC using approximate solver to solve P_star and u_star is working and correct. I will update my page later. Next goal is Roe solver.
2. reading books and papers (all previously mentioned).
3. I am trying to model a more general method to initialize the boundary condition for co-rotating binaries (with prescribed elliptic orbit, arbitrary wind and orbital velocity) using interpolation method. I roughly get the idea and I am confident it should work.

1. Read chapter 8,9,10,11. I roughly understand what are their purposes and means. There are some typos in the book and I also did not figure out some abbreviated derivation. I would like to reread and resolve the derivation later (perhaps after a long time) because it does not hinder my understanding now.
2. Updated my page.
3. My father will come to NYC next Tuesday, can I have 3 days off?

1. After days of crazy debugging, I finally get the correct Godunov solver. I will post the correct picture to my page this week. Thank god I can proceed to next chapter!

2. Created my page, added to the group member page. I will add more material to it.

1. Read 1-12 chapters in "Statistical Physics" by Kerson Huang (I did this not just in one week but since the beginning of this summer). It is about ensembles, time-series analysis, stochastic process and thermal dynamics - all of them are very helpful in understanding astrophysics (diffusion, viscosity and heat conduction) and even cosmology (distribution and ensembles).
2. Read chapter 5 and 6 in "Riemann Solver". After reading these, I feel that Fourier analysis can be powerful in Riemann solver, and could be even powerful when combined with AMR.
3. Plan to read some numerical analysis books, a potential one is "Finite difference methods for ordinary and partial differential equations : steady-state and time-dependent problems" by Randall J. LeVeque, UW. It is best know both numerical stuff and physics stuff so that one may tell what is physical and what is unphysical in simulations.

1. Finished programming complete Riemann solver, it is correct.
2. attended HEDP summer school, discovered a lot of fun in fusion technology, high energy density physics, astrophysics and computational science.

1. I will leave for OSU to attend HEDP summer school from July 14th - July 20th
2. I have finished programming Riemann solver, i.e. 1D shock tube. It is correct. But I did not consider the vacuum case.

1. I know what is chapter 1-4 of Riemann Solver talking about now, and I roughly know the flow chart of shock tube problem program.
2. Learnt how to write f90 program, I can use f90 now.
3. Ran several CorotatingBinary on bluehive, the density at stationary point is higher than the companion which is not very reasonable.

1. Read Chapter 3 and 4 of Riemann Solver.
2. Read Martin's paper in depth ​http://arxiv.org/abs/1211.1672v2 "The Formation and Evolution of Wind-capture Disks in Binary system". I built a different model (It could include thermal non-uniformity and I hope to include orbital change in the future) for elliptic orbit, I want to discuss it on Monday's seminar.

There are three references in this model.

\rho_w=\frac{\dot{m}_p}{4\pi\|\mathbf{r'}\|^2\mathbf{v}_w}

Outside the computation box:

\mathbf{a_r}=\mathbf{a_i}-2\mathbf{\Omega}\times\mathbf{v}_w-\mathbf{\Omega}\times(\mathbf{\Omega}\times\mathbf{r'})-\frac{d\mathbf{\Omega}}{dt}\times\mathbf{r'}

\mathbf{r'}-\mathbf{r}=\mathbf{L}

We can get \mathbf{u} field in the reference of co-rotating primary after integration. So that we can get the consistent flux on the boundary of the computation box.

Inside the computation box (include the boundaries): First compute the \mathbf{v} field in the reference of co-rotation companion.

\mathbf{v}=\mathbf{u}+\mathbf{\Omega}\times\mathbf{r}

Special attention should be paid to the boundaries since they may swell or shrink along the orbital motion. Which introduce additional flux.

\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0

\frac{(\partial \rho \mathbf{v})}{\partial t}+\nabla\cdot(\rho\mathbf{v}\mathbf{v})=-\nabla p-\rho\nabla\phi-2\rho\mathbf{\Omega}\times\mathbf{v}-\rho\mathbf{\Omega}\times(\mathbf{\Omega}\times\mathbf{r})-\rho\frac{d\mathbf{\Omega}}{dt}\times\mathbf{r}

\frac{\partial E}{\partial t}+\nabla\cdot[(E+p)\mathbf{v}]=-\rho \dot{Q}_{cool}+\rho\frac{\partial\phi}{\partial t}

E=\rho(\frac{1}{2}\mathbf{v}^2+e)

e=\frac{p}{(\gamma-1)\rho}

T=\frac{p}{\rho R}

To include orbital change:

\frac{1}{2} m(\dot{r}^2+r^2\dot{\theta}^2)+Q(r)=C+\int_{\mathbf{s}(t)}^{\mathbf{s}(t+\triangle t)}\mathbf{f}(\mathbf{s})\cdot \mathbf{ds}

\frac{d(mr^2\dot{\theta})}{dt}=L(\mathbf{s})

m may be changing. \mathbf{s}(t) is the relative location of companion to the primary, and is basically in the inertial reference. It's initial condition is exactly \mathbf{L}.

Therefore, f(\mathbf{s}) and L(\mathbf{s}) is important to incorporate drag force and drag moment. These will give coupled orbital equation.

Question: In low-mass companion binary systems, what is the order in all these phenomenons below :

thermal non-uniformity, orbital change, mass change, frictional force, rotation of the companion and even EM field

Remark: Resolution is important as computational matter, order is the "resolution" in physics and modelling.

1. Read chapter 2 of Riemann Solver.
2. Review the previous papers.
3. ​http://www.netlib.org/lapack/

1. ​http://adsabs.harvard.edu/abs/1944MNRAS.104..273B Bondi and Hoyles' paper on the mechanism of accretion by stars.
2. ​http://adsabs.harvard.edu/doi/10.1093/mnras/stt725 Formation and Evolution in binary by our group
3. ​http://arxiv.org/abs/astro-ph/0406166 A Review of Bondi—Hoyle—Lyttleton Accretion
4. First chapter of Toro's book
5. ran some co-rotatingbinary (can find rho of ambient, cant find rho of the companion.)

1. read ​http://adsabs.harvard.edu/abs/2006MNRAS.370.2004N Common envelopes and their relevance for planetary nebulae. It is written by Jason.
2. successfully run CorotatingBinary.
3. tried to give the governing equations of a specific binary system. I am still working on that. equations include: continuity equation; momentum equation; energy equation; equation of state; Maxwell's equation. gravity (poisson)
4. keep practising f90/95 in linux

I did not forget to learn the theory of AMR, it's on my to do list.

1. Learning Fortran 90/95 in Linux environment.
2. Reading an article "The Evolution of Binary Systems", here is the link: ​http://www-astro.physics.ox.ac.uk/~podsi/binaries.pdf I read some parts of it and gain knowledge in binary stars.
3. Reading "Common Envelops in Binary Systems" ​http://adsabs.harvard.edu/abs/1993PASP..105.1373I I read a few pages.
4. Reviewed R-T instability and K-H instability (still in progress). Because I think K-H instability is important in binary evolution.
5. Learning AstroBEAR, tried to run example.

I am now:

1. reading an introduction to close binary stars, R. W. Hilditch
2. reading principles of astrophysical fluid dynamics, Clarke and Carswell
3. reading Riemann solver, Toro

and papers:)