Summer school

Got to play a little with Gandalf SPH code and simulated the collapse of a singular isothermal sphere. Learned that python can process chombo hdf files, among lots of other cool astrophysics things.

Amira

Today I followed up with Christoph Federrath about the Amira vs. Disperse paper. I largely fleshed out the paper structure and figures and am excited to begin this project.

The goal of the paper will be to address the following 4 questions:

Do these different programs:

1. Find the same # of filaments
2. Find the same filaments
3. Find "realistic" filaments
4. Find "actual" filaments

For 3., Christoph's paper (link: ​http://www.mso.anu.edu.au/~chfeder/pubs/filaments/Federrath_filaments.pdf) talks about the differences in filament count between adding different physics to simulations. He finds that as you add more and more physics (turbulence, gravity, outflows), more and more filaments are formed in the simulations (no real surprise here). However, of greater interest is the fact that the average filament width is r=0.1 pc. He finds this across the board (irrespective of Disperse's sensitivity threshold and simulation resolution), so long as there is turbulence in the simulation. Since this matches observation (see Herschel findings), I called this "realistic" filaments. I think it would of interest if Amira also finds this result.

For 4., I described how we were throwing the idea around of creating a highly idealized simulation of 'definite' filaments and then feeding this through both Disperse and Amira. Do they both track all of the actual filaments we created? (I would imagine they would). Next we add some noise to the filaments progressively and watch if and when Disperse and Amira diverge in their findings.

I told him if he could send me the 2D CDM of one of his simulations with turbulence I could run it through Amira to get this started. I think we should think about how Disperse's 'persistence threshold' might be compared to parameters Amira uses and how to try to normalize the two programs the best we can..

Colliding Flows

Going to work on discussion and turbulent Jeans mass estimate these next couple of days

I tested the new pulsed jet module on Rice's machine DAVinCI, and it seems to work well. This pulsed jet is 3D, MHD, with non-equilibrium cooling.

movie​

movie​

The code ran pretty well, producing chombos at a steady pace with no restarts. Here are some numbers to show how well this run did on DAVinCI:

The run does slow down as it moves along, so it's hard to say how long the full run would take at this resolution. The code estimates 11.2 more hrs after frame 30, but the actual time is probably closer to 24 hrs. If I assume that the time between chombos grows linearly, then I can estimate the total runtime to be about 24 hours.

Thus, the total run at 32 cells/Rjet would take 24 hours on 576 cores, costing 4.70. This is the maximum number of cores that I can use on DAVinCI, and I can only submit for a maximum of 8 hours at a time. So even at this low resolution a run is going to take more than a day if you include wait time and restarting. If we double the resolution to 64 cells/Rjet, this run could take 2 weeks or more instead of 1 day, and the cost would be more than 50.

UPDATE

Two new models: hydro (no MHD), and an MHD with beta = 1 (same as above) but with precession. Precession amplitude is 3 degrees and the precession frequency is 2.5 times the pulsation frequency.

movie​

movie​

• No co-rotation

With Temp profile old

​movie

Mass loss Rate: 2.13E+09 gram per sec

With Temp profile		​movie	Mass loss Rate: 3.04E+09 gram per sec
No Temp profile		​movie	Mass loss Rate: 3.82E+09 gram per sec

• Very Low Ambient Density

• Wind on all boundaries

​20-zones-per-radii

Movie showing streamlines of the velocity field (magenta), magnetic field lines, the mach surface (black) and log density (pseudocolor) as the simulation relaxes after injecting a dipole field.

movie​

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.

New 2-D Mach stem runs with lower M

I tested gamma = 1.4 with 4 different separations: 4.5, 5, 6, 6.5 rclump. I chose these separations because they are right on the edges of the different regimes (see figure below).

Velocity was lowered from 50 km/s to 15 km/s resulting in Mach number being decreased from 5.18 to 1.55. I had to lower the maximum of the density scale from 35 to 16 since these are weaker shocks (number density in computational units of 10,000 cm-3).

movie​

movie​

movie​

movie​

The table below summarizes the Mach stem behavior for the relevant models after 75 years. To characterize this time in terms of clump-crossing times, it would be about 40 for the original M = 5.18 runs and about 12 for the new M = 1.55 runs. (TMS = Transient Mach Stem, SB = Single Bow, SMS = Stable Mach Stem, RR = Regular Reflection)

As you can see, with a lower M, we now get a single bow in every model. I think this is because there is less pre-shock ram pressure to keep the enhanced post-shock thermal pressure at bay. In regards to Mach stems, I see two possibilities: 1) Achieving a SMS is impossible, and the behavior will transition to RR at larger separations, or 2) the SMS region has just shifted to larger separations. Either way, I want to do some more runs to see if I can find the SMS regime for this lower M.

Other Stuff

• Debugging the error in the projections. The error I'm getting now is not where I expected it to be which is confusing. It is in data_declarations.f90 where the Info object is initialized. Maybe I need to run with more cores to get by this part of the code, and trigger the real error?

• Still haven't received any confirmation about my guest account on the Rice machines.

• I'm in the process of moving some data off of BlueHive to make room for some 3D pulsed jet runs.

• PAPER EDITS!!! There is a lot to be done, so I'm mainly busy with this now.

1. Checked out the WireTurbulence data
2. Ran code
3. Used Visit to analyze result
4. Output the movie

And here is schematic showing the star/planet surfaces, their sonic surfaces, and the location of the bow shock, hill radius, and coriolis radius. Note the bow shock will likely be closer to the planet due to the planet wind being sub-parker and that the coriolis radius will probably be larger for the same reason.

And here is the same setup, but with a stellar density 10x higher

• Level 0 resolution is 1600x160x160
• Simulations started at level 2
• Backed off to level 1 and restarted MHD from beginning
• Hydro run continued from frame 62 but at level 1.

Level	Hydro	MHD
2	0-84	0-39
1	62-200	0-42

Hydro run at frame 200

​

movie​

—-

MHD run at frame 38

movie​

• re-applied for a guest account on Rice machines, waiting to hear back from Pat

• I'm working on some low M, 2-D Mach stem simulations for Pat. This is for the referee's report; to see how low M affects Mach stem formation in simulations.

• Working on re-generating some of the paper figures for my 3-D clump/bow shock paper. Mainly want to increase resolution of images, but I'm having some issues. See images below…

The central stellar cluster's wind extracts angular momentum from the circumnuclear disk's inner rim, therefore the inner rim migrates

inwards. In the following I derive the timescale for the inwards migration of the CND's inner rim. Please see our paper draft for more

details.

We assume that angular momentum is extracted from a ring of radius s and radial extension \Delta l, therefore it has the mass M_{\text{D}} = 2 \pi \Sigma s \Delta l. This ring has a total angular momentum of J=j M_{\text{D}} with specific angular momentum j.

After a time dt the wind adds the mass f \dot{M} dt to the ring, where \dot{M} is the wind's outflow rate and f is the

fraction of the wind that interacts with the disk. Thus the ring has the mass M_{D} + f \dot{M} dt and the specific angular momentum

j' = \frac{J}{M_{D} + f \dot{M} dt} = \frac{j}{1+\frac{f \dot{M}}{M_{\text{D}}} dt}.

Furthermore the ring moves to s-ds, and because j = \sqrt{GMs} and j' = \sqrt{GM(s-ds)} we get

dt = \frac{M_{\text{D}}}{f \dot{M}} \left( \sqrt{\frac{s}{s-ds}} - 1\right) = \frac{M_{\text{D}}}{f \dot{M}} \frac{ds}{2 s}

where I used a taylor expansion in the last step.

Using M_{\text{D}} = 2 \pi \Sigma s \Delta l, f=\frac{h}{2s} and integrating this equation gives the time the inner rim needs to

migrate inwards:

\tau = \int dt = \int \frac{M_{\text{D}}}{f \dot{M}} \frac{ds}{2 s} = \frac{2 \pi \Sigma \Delta l}{h \dot{M}} \int s ds = \frac{\pi \Sigma \Delta l}{h \dot{M}} \left( s_0^2 - s_\text{i}^2\right),

where s_0 is the disk's initial inner rim and s_\text{i} the location of the rim after inwards migration.

In our simulations s_\text{i} = 0.5 \text{pc}, for s_0 = 1 \text{pc} we get \tau = 22 \cdot 10^4 \text{yrs}, for s_0 = 2 \text{pc} we get \tau = 112 \cdot 10^4 \text{yrs}. The actual time the inner rim needs to move inwards (as observed in the simulations) is 20 \cdot 10^4 \text{yrs} and 65 \cdot 10^4 \text{yrs}, respectively. For the s_0 = 1 \text{pc} simulation this is pretty close, for the s_0 = 2 \text{pc} simulation the deviation is less than a factor of two, which I find an acceptable deviation considering that the above mentioned derivation is only an order of magnitude estimate.

Thoughts on parker wind

The solution for the parker wind is found by solving the equations below

Note however, that the actual radius of the planet (and lambda) does not really matter…

All that a higher lambda implies, is a smaller planet radius (compared to the sonic radius) - and sampling more and more of a subsonic atmosphere before having a hard boundary.

Also - per our discussion the other day, it seems that

\frac{R_\Omega}{a} = \frac{v_{esc}}{a\Omega} = \sqrt{\frac{2GM_p}{r_p a^2}\frac{a^3}{G\left(M_p+M_s \right) }} = \sqrt{\frac{2aq}{r_p \left(q+1 \right)}}

So this effect is more pronounced when the orbital sepration is small and when the planet radius and mass ratio is close to - ie Hot Jupiters

I believe the planet we want is HD 209458 b

​http://exoplanet.eu/catalog/hd_209458_b/

I'm back from getting married and having a honeymoon !

Next week I travel to Heidelberg for a week long summer school on ISM dynamics (bringing schedule to meeting, looks very good). I just booked flights and am preparing this trip.

I am also working on journal club happening in 2 days (got a fun paper to review) and reading up on literature from the last month until I get the draft back to pick up again

• 1. low-res result
  

\gamma=5/3; \omega_{rot}=0.5 with analytic cooling; Still running…

​8 zones per radii

• 2. zCooling is not working with the current version of cold. Haven't figured out how to fix it yet…

• 3. Higher Resolution runs

qTolerance              = .10,.30,.30,1d30,1d30,1d30,1d30,1d30,1d30