# Balick et al. 2013

Old site https://clover.pas.rochester.edu/trac/astrobear/blog/crl618Figures

**Paper Section 4, COMPARISON TO MODEL SIMULATIONS**

**4.1. Models**

Here we describe the methodology and implementation of our CRL618 models in detail. We have carried out two Eulerian-grid numerical simulations to follow the formation of nebular lobes via the propagation of a bullet or a jet into a stratified and ridged ambient medium. Fluid dynamics is followed in 3 dimensions using the equations of radiation hydrodynamics. The eﬀects of optically thin cooling have been included using the cooling tables of Dalgarno & McCray (1972). The hydrodynamic equations are solved with the adaptive mesh reﬁnement (AMR) numerical code AstroBEAR2.0 [FOOTNOTE: https://clover.pas.rochester.edu/trac/astrobear/wiki]. In particular, the Euler equations with cooling source terms are solved using a second-order MUSCL Hancock shock capturing scheme and Marquina ﬂux functions (Cunningham et al. 2009). While AstroBEAR is able to solve the equations of magnetohydrodynamics (MHD) and to compute several microphysical processes, such as gas self-gravity and heat conduction, we do not consider these processes in the present study. The computational domain is a rectangle with dimensions 0<*x*<24000 AU, -4000<*y,z*<4000 AU. We use a coarse grid with 200x100x100 cells along with two adaptive refinement levels which increases the grid resolution by a factor of 4; the simulations attain a maximum resolution of 20 AU. Typical simulation ﬂow times are of order 400 yr. We use BlueHive [FOOTNOTE: http://www.circ.rochester.edu/wiki/index.php/BlueHive_Cluster] and Blue Gene/P [FOOTNOTE: http://www.circ.rochester.edu/wiki/index.php/Blue_Gene/P] —IBM's parallel
cluster and supercomputer, respectively— which are maintained by the Center for Integrated Research Computing of the University of Rochester. We ran simulations for about 2days, using 256 processors.

**4.1.1. Initial conditions**

We set the ambient medium with a static velocity (**V**=0) field and a density profile which decreases with distance form the origin (0,0,0) and has a series of periodic spherical ridges spaced by ~333 AU. In detail,

*n_amb(r)= 300 / (r/500AU) ^{2} *
.5 + SIGMA_i exp(-[6d0*{(r/500AU)-5/3-i}]^{2})* particles / cm

^{3},

where *r* is the radial distance form the origin and *i* is an integer number. Our choice for such structure is based on Fig. 1 of NF+07, where it is clear that the AGB “spherical halo” containing the lobes of CRL 618 has a stratified structure with semi-periodic ridges of enhanced density separated by ≈1″ (section 3). Observational studies of AGB winds suggests that they expand isotropically with mass-loss rates and velocities of order 10^{-5} M_{sun}/yr and 20km/s, respectively (see e.g. Hrivnak et al. 1989; Bujarrabal et al. 2001). Our model AGB circumstellar medium is static, however, and we do not expect this difference to affect the results of our model because (i) we explore short, ~200yr, nebular expansion times, (ii) the axial velocity of the bullet/jet is more than an order of magnitude faster (300km/s; below) that the expected AGB wind terminal velocities.

**BULLET MODEL**

We base our bullet model in that of Dennis et al. (2008). At time=0 yr, a spherical bullet is placed in the grid at coordinates (*r_b*, 0, 0), where *r_b*=500AU and it is the bullet's, and the jet's, radius, with an axial velocity of 300km/s and a density profile given by

*n_b (r) = min( 100*n_amb(r) , r_b ^{2}/r^{2} )* particles / cm

^{3}.

**JET MODEL**

For the jet model, we continually inject gas at the grid cells which are located at the bottom face of our computational domain, within *r<r_b*. The injected gas has a constant axial velocity of 300km/s and a constant density of

*n_j = max(n_b(r)) particles / cm ^{3}*.

**4.1.2. Evolution**.

We present the results of our CRL 618 numerical simulations using 2d maps of density and velocity field (Figure 3, columns 1 and 3, respectively), and 1d plots of density (Figure 3, columns 2 and 4, top panel) and velocity (Figure 3, columns 2 and 4, bottom panel). All the images in Figure 3 were done using the computational data located at the middle plane of the computational doman, thus no integration of density or velocity along the line of sight was carried out.

The 2d maps in columns 1 and 3 of Figure 3 show the following structures. The bullet/jet (left/right) are the densest whitest central features. Below the bullet, or to the right of the jet (without loss of generality), we see the lobe that is formed by the bullet/jet, which is separated from the ambient medium by a contact discontinuity. Beyond the lobe we see the ambient medium which is stratified and shows a series of concentric spherical shells given by the initial conditions (above). The vertical short red, or blue, lines at the top of the maps show the positions at which the profiles in Figure 3 (columns 2 and 4) have been taken.

The bullet and the jet propagate at 300km/s, away from the ambient medium's densest region (upwards), forming an elongated lobe. We see that the expansion of the lobes is predominantly axial; the major to minor lobes' axes ratio is always > 2. We find that the lobes formed by the bullet and the jet are quite similar (compare left vs. right panels in columns 1 and 3, Figure 3), except for the material located within a bullet/jet radii from the axis. As the bullet (jet) penetrates the ambient's shells, these form regularly separated vertebrae-looking features along the lobe. The axial velocity of these features decreases with radial distance from the axis.

We see that the bullet shrinks as it propagates though the ambient medium (compare columns 1 and 3, Figure 3). Similarly, the head of the jet seems to adopt a conical up-ward pointing geometry as it propagates. Such effect is due to ablation of the material located at the working surface of the bullet/jet as it interacts with the ambient medium. Such effect is consisten with the models of Dennis et al. (2008).

The density profiles in Figure 3 (columns 2 and 4, top panel) show the following.

the evolution of the bullet and the jet quantitatively, as a function of time and position inside the lobes.

*Bruce's wish**: that Martin/Adam explain why the speeds of the gas in the body of the lobes follow a Hubble-like pattern since Bruce makes the case in section 5 that this is where the H_2 emission is likely to arise.*

# Poloidal Collimation

Matt, Winglee & Bohm, 2003, MNRAS, 345, 660, http://adsabs.harvard.edu/abs/2003MNRAS.345..660M

Spruit, H.~C., Foglizzo, T., & Stehle, R., 1997, mnras, 288, 333 http://adsabs.harvard.edu/abs/1997MNRAS.288..333S

## 19 dec '12

## 18 dec '12

Ditto but for the run of Nov. 21 |

# Meeting Update Dec.18

**Triggered Star Formation** Found a place in my problem module where the clump radius is put in by number not by reference. This caused a discrepancy in terms of clump radius between the problem module and the problem.data file that I used to run the simulations. Upon fixing this bug, I was able to get different distributions that can match Boss' numbers by tuning xi, clump radius, central density to match Boss' numbers (previously I always get one setup because the changes in clump radius were only in problem.data but not reflected in the problem module). In this way, I can get fairly stable setups that has the same density distribution and temperature inside the clump as Boss', although I still have to artificially increase or decrease the ambient density to what they were using. Following is one of such settings in 1/10 of clump crossing time (This setup is confirmed to stay for at least one crossing time). The ambient density is reduced by a factor of 100. The xi in this case is about 6.78 so it's still in the "unstable" regime. Currently working on preliminary runs with shocks.

2D:

http://www.pas.rochester.edu/~shuleli/meeting_1210/bss.gif

Cut:

http://www.pas.rochester.edu/~shuleli/meeting_1210/bsc.gif

**LLE stuff** Ported the implicit solver from a version a year ago to the latest gold version. Made the code to compile successfully and run in fixed grid, still needs work in AMR and testing (considering self-similar heat conduction test as 1-D quantitative test). Future minor change could also include making self gravity and implicit heat conduction to work together (currently you cannot turn on both at the same time). I expect to be able to have a fully functioning version by our next meeting.

# Meeting Update 12/18/2012 - Eddie

- did a few preliminary runs of 1D radiative shocks with magnetic fields. All three runs are with shock velocity 50 km/s and preshock density 100 cm
^{-3. In uG, black is 1, blue is 30, and red is 100. }

- working on a script to automatically do all the runs that Pat wants
- constructed an Al cooling table to use for lab jet sims for Francisco, and compared it to the plot in Post et al., 1977…

Post Al Cooling Rate | My Constructed Cooling Table |
---|---|

# Christine's Visit

==============

Attached is an ApJ letter that I've published on my experiment and the astrophysical motivation for it. What I would like to do with AstroBEAR is create the Cataclysmic Variable (CV) detailed in the first paragraph of the introduction. The CV related to this work is considered non-magnetic because the secondary star in the binary system donates mass to the white dwarf (WD) in the orbital plane, (not along field lines to a polar axis).

The localized area of interest is the collision region at which the accreting stream impacts the formed accretion disk around the WD. I've also attached an ApJ paper that did simulations on this (comparing isothermal and adiabatic EOS's) because they detail parameters of the CV system. You will see in my paper the connection of this astrophysical system to my laser experiment, but to re-iterate the "problem I want to solve" is the dynamics of the shocked stream at the accretion disk edge. I have new data, which I showed to your group in July, that hasn't been published yet, but shows perhaps some stagnation and a diverted stream in an oblique shock scenario. This corresponds most likely to some optical thickness in the shock system, so it would be interesting to do the CV simulation looking more closely at the shocks that form in it and how mass moves around the collision region, as a function of scale height of the disk.

I can offer more details with questions/concerns, but my paper in combination with the Armitage and Livio one (sections 1 through 2.3) offer a good overview of the system.

I am currently on a different experiment in Livermore, but I am going to start exploring the AstroBEAR wiki, etc, this week! Thank you for setting up an account for me.

=====================

APJ letter http://www.pas.rochester.edu/~bliu/Christine/ApJL452485p4.pdf

Armitage and Livio http://www.pas.rochester.edu/~bliu/Christine/36557.pdf

# Meeting Update 12/18/2012 -- Baowei

- Golden Version & Machines

- Youtube Channel for AstroBEAR

- Worked on
- coding and testing Revision 1174:9b0df1e0242b
- Help install for new machines and new users
- Testing Shule's AblativeRT module

# Meeting update

I updated my page with things currently on my plate, as well as pages for them:

https://clover.pas.rochester.edu/trac/astrobear/wiki/u/erica

Namely,

- I updated a BE paper figures page: https://clover.pas.rochester.edu/trac/astrobear/wiki/u/erica/BERunsFigures
- Made a coding page for the exact riemann solver I am studying: https://clover.pas.rochester.edu/trac/astrobear/wiki/u/erica/ExactRiemannSolver
- Made a page where I can publish drafts for my paper: https://clover.pas.rochester.edu/trac/astrobear/wiki/u/u/erica/BEPaper

The paper introduction is nearly done (not published on that page yet as I need to do some polishing of the text and my references to papers/results). I have begun writing a Fortran program for the Riemann solver following the outline I posted on the page above. I am writing it in the fashion of an astrobear module, for good practice. I also thought about plotting the B&P run. B&P normalized their data to unperturbed values, but I feel it is more clear to use perturbed values (which will both normalize the density to the perturbed central density, as well as change the truncation non-dimensional radius to beyond the critical xi). This will most clearly show the effect of perturbation - to increase xi beyond its critical value. Also, to compare the density to the unperturbed central density seems to be less intuitive as comparing to the perturbed central density. As far as normalizing the times in the plots, I derived a sound crossing time for the BE sphere, consisting of variables specific to the BE problem. I have normalized all the times to this time for my plots.

# Meeting Recap 12/17/2012

- Spent some time working on Christina's setup and improving refinement criteria. 'I' shape due to density gradients produced due to bending magnetic field lines and strong shear. In 3D this I is a cylinder. So I modified module to use three refinement objects.
- First maintain refinement to level 3 within a 2 pc wide region around the initial interface (nested inside of coarser patches)
- Next allow for density gradient based refinement to level 3 within a 10 pc region centered on the initial interface (nested inside of coarser patches)
- And finally allow for jeans length based refinement (64 cells per jeans length) all the way to level 4 so that collapsing structures will be better resolved.

- 3D MHD with self-gravity on bluestreak runs at 1475 cell updates/second/core
- On 32 bluestreak cores, at 32
^{3}+ 3 it was spending 76% of the time advancing, 23% solving elliptic equations, and < 1% everything else… and was slotted to finish in 8 days (30 Myr). - In principle if everything scales up - the initial resolution of 128
^{3}+4 could require between 256 and 4096 times as many resources depending on the filling fraction of the highest level… - At 256 this would give 16 days on ¼ of bluestreak (4000) which seems doable? This is 1.5 Mhrs of CPU time on bluestreak - or probably .4 Mhrs on kraken

- Things to work on over break…
- MHD convergence - fast waves now converge - slow waves still do not.
- PPM convergence
- FLD solver
- Documentation ???
- Developer's manual ???

- Bash script for running parameter studies…
interporders=( 1 2 3 ) Gmxs=( 100 200 400 ) Waves=( 1 2 3 4 5 6 7 ) mkdir out mkdir testdata for (( i=0;i<${#interporders[@]};i++)); do interporder=${interporders[i]} ../../scripts/subst.s InterpOrder $interporder solver.data for (( j=0;j<${#Waves[@]};j++)); do wave=${Waves[j]} ../../scripts/subst.s WaveFamily $wave problem.data for (( m=0;m<${#Gmxs[@]};m++)); do Gmx=${Gmxs[m]} ../../scripts/subst.s Gmx $Gmx,1,1 global.data mpirun -n 1 ../../astrobear > output.out grep dx output.out | awk '{print $3, $4}' >> scaling_${interporder}_${j}.curve mv out/data.curve data_${interporder}_${j}_${Gmx}.curve cat output.out done done done

# Meeting Update - Jason and Ivan 1217

Made a plot of accretion rates.

Each model has 2 lines: Total accretion (rate at which the central particle accretes), Tracer accretion (rate at which the central particle accretes disk material).

Two things worth noticing: 1) The central particle is accreting a big deal of material from the ambient and much less from the disk.
2) The 20 Jupiter Mass behaves differently compared to the other two models.

# Martin's update, 18 dec '12

**CRL 618**, new data and plots @ top of https://clover.pas.rochester.edu/trac/astrobear/blog/crl618Figures

# Meeting Update Dec.10

**cloud stability** Using my setup to run the cloud stability for one crossing time:

http://www.pas.rochester.edu/~shuleli/meeting_1210/mysetup.gif

We can see that the cloud breathing, and there is no significant expansion or contraction. There is a movie for the cut-through density here:

http://www.pas.rochester.edu/~shuleli/meeting_1210/mysetupcut.gif

This setup is obtained by feeding Boss 2010's central and ambient density to the BE module in the code. The problem is that the density at the edge of the clump is not the same as Boss' number. We have 1000 density contrast at the edge of the cloud (radius = 0.058pc), they have 100 density contrast. Also, our temperature is 124k inside the cloud while they have 10k.

If we increase the cloud radius to 4 times the original, we can get the edge density drop to Boss' value. But the cloud seems to collapse under 1 crossing time scale:

http://www.pas.rochester.edu/~shuleli/meeting_1210/myset2.gif

In order to match exactly their profile, I tried increasing the density contrast (making the profile steeper). By doing this, the supposed ambient density will be 10 times lower than the original (1e-23 g/cc vs 1e-22 g/cc). I pad the ambient by the Boss' ambient density and match the temperature at the cloud edge. This results in a profile exactly the same as theirs, with the same in-cloud temperature 10K. However, this setup collapses for under 1/10 crossing time (about 1 shock crossing time):

http://www.pas.rochester.edu/~shuleli/meeting_1210/bosssetup.gif

Here is a movie for the cut-through density:

http://www.pas.rochester.edu/~shuleli/meeting_1210/bosscut.gif

If we do not pad the ambient by the original value and just use what is calculated in the code, we get:

http://www.pas.rochester.edu/~shuleli/meeting_1210/boss2.gif

At this point, I think using our own setup should be just fine since it produces a fairly stable cloud for 1 crossing time scale. However, since the density at the edge of the cloud is about 1000 times the ambient and only 1.5 times below the central density, the incoming wind will have a peak density close to the central density of the cloud if we still want to match the wind density with the cloud edge density. I'm not sure currently how important this condition is to the result. I think we need to do some wind-cloud interaction simulations next to verify how to setup the wind.

**Papers** Finished the revised version of the clump paper and the referee response. The only tiny problem is that in his comments 5 and 6, the referee suggested the restructuring of the material, which is not done in the revision (we just argued that our version is the best).

http://www.pas.rochester.edu/~shuleli/meeting_1210/paper_shule_1210.pdf

http://www.pas.rochester.edu/~shuleli/meeting_1210/referee_response_1210.pdf

# Meeting Update 12/10/2012 -- Baowei

- Golden Version & Blue Streak
- Updated Information about Performance of Blue Streak https://clover.pas.rochester.edu/trac/astrobear/blog/johannjc11132012#comment-5
- Considering the speed difference of blue gene Q & P, current AstroBEAR code should be able to run faster than now. Testing individual library.

- Tickets
- New ticket: #270(standard out "walltime remaning" accuracy)

- Worked on
- Ticket #265
- parallel hdf5
- MUSCL-Hancock scheme

# Meeting update

Last week I finished reading Toro ch.'s 1-4 and met with Adam on the first program I will be writing. Basically it will be setting up a so-called Exact Riemann solver, which will take 2 initial data states (left and right) and solve for the pressure of an intermediate state in the "star region". This pressure will be found using an iterative approach from an initial approximation to the actual pressure. Once this is found, the remaining fluid variables can be solved for once the types of waves on either side of the star region have been identified as either shock or rarefaction. This is because each type of non linear wave has a set of conditions that apply for the density and velocity across them (only pressure is necessarily constant across the middle wave - the contact discontinuity). With this in mind, once P is found, the program will follow a flow chart to determine which conditions must be met, and hence which equation or set of equations should be solved to find the remaining fluid variables of the problem.

I also began working on writing the introduction and thinking about redoing the plots for the paper so that the curves show the actual resolution (which increases over time of simulation). I am thinking that for the plots, I can adjust my program to use just a constant line out in the center region of the sphere (where I think the angular distribution of data was symmetric), and an averaging scheme over the outer shells. In this way, the data sampled from Visit should automatically adjust for smaller radii as more refinement is added later on in the simulation.

Here is what I have for the Intro so far:

- Various properties of the “Bonnor-Ebert” (BE) sphere, a hydrostatic sphere in pressure equilibrium with its ambient environment, make it a good candidate for numerical modeling of protostellar collapse. First, as a candidate star forming structure is envisaged as gravitationally bound and unstable, it is easy to imagine a protostar evolving from an initially hydrostatic configuration. Indeed, spherical clumps have been observed in or near hydrostatic equilibrium, such as the Bok Globule, B68 (Myers). Second, the stability criterion against gravitational collapse has been worked out analytically. Third, pushing the sphere out of the stability regime with various physical perturbations illuminate collapse characteristics. Such features of the collapse may help advance single star formation theory as well as provide clues to observational astronomers in identifying potential star forming sites.

- While the collapse of a BE sphere has been studied extensively over the years, the literature reveals studies of the BE sphere in precarious and unphysical conditions. The BE sphere has largely been modeled as residing in low density ambient mediums, so to seemingly isolate the collapse of the sphere from the environment. Examples of previous works - in response to Shu’s claim that all collapse would be inside out - the BE sphere does not collapse in this way…

- Although modeling the collapse in this way does illuminate unique features of the BE collapse, it inadvertently constrains the model. If the BE sphere might be considered as an early basic model of star formation, it should be examined in a more physically plausible setting. That is, the BE sphere should be tied to its parent cloud as it would be in nature, by being modeled with no discontinuous boundary between BE sphere and parent clump. The collapse properties of such a simulation would be more physically accurate to any actual situation involving the collapse of a protostellar BE sphere. Discuss those papers that began to address this - Myers and Hannebelle.

- How this paper is different from previous similar work.

Self Gravity papers to read? https://clover.pas.rochester.edu/trac/astrobear/wiki/ExternalLinks#Self-GravityCoreCollapsePapers https://clover.pas.rochester.edu/trac/astrobear/wiki/SelfGravity

# Jason and Ivan - Meeting Update 1012

- We started analyzing the runs and making plots.

The accretion disks start with the same pressure and different densities, as suggested in paper meeting we had back then, as a consequence they also differ in temperatures.

Let's remember the initial conditions briefly:

Mass | Density | Temperature |

20 | ~20.0 g/cc | 5d3 K |

10 | ~10.0 g/cc | 1d4 K |

50 | ~5.0 g/cc | 5d4 K |

10 | ~1.0 g/cc | 1d5 K |

Although the disks start with different initial conditions, they appear to settle at the similar temperatures. We decided to measure the disk mass using temperature to determine what material belongs to the disk, as of right now this seems to be the best option, suggestions are welcome.

The threshold we decided to use is 2d6K. To visualize this value the plot below can come in handy:

This being said, below are two plots of our 4 runs, the first one is relative mass, the second is absolute mass. Part of the 5 Jupiter Mass line is dotted (predicted mass) since we decided to extend that run a bit to get it to 10 orbits.

- There is a vortex-like feature going on with the 20 Jupiter Mass. It does not compromise the general picture, however, it is something to think about. What could be causing it ? Movie here

- We need to check accretion rates since they seem a little bit too high. I will check them in detail later on this week.

# Martin's update, 11 Dec '12

**Binary disks paper**. Working on Noam's comments.

**CLR 618**.

- New plots: https://clover.pas.rochester.edu/trac/astrobear/blog/crl618Figures
- Writing "numerics" section of the paper.

**PN jet power**. Helping Scott to run PN jets. He's already running in grass and playing with module parameters. Movies/images soon.

**AGN**.

- New data: http://www.pas.rochester.edu/~martinhe/2012/agnTrunc.pdf
- Writing Milan's proceedings and running sims.

# Meeting Update 12/03/2012 - Baowei

- Golden Version & Blue Streak
- AstroBEAR runs a little slower on Blue Gene/P than on blue streak: https://clover.pas.rochester.edu/trac/astrobear/blog/johannjc11132012#comment-4 ;
- Tried to compile AstroBEAR with the fftw3 of IBMmath but haven't succeeded yet.

- Tickets:

- Users:
- IO's chombo files: http://www.pas.rochester.edu/~bliu/Yat_Tien_UCLA/CND_chombo/
- Christine will visit us during the week of Jan 14th

# Meeting Update 12/03/2012 - Eddie

- made progress on understanding cylindrical source terms and algorithm in the Skinner & Ostriker paper on Athena in cylindrical coordinates. ehansen11272012
- my token is activated now…how are we splitting up our resources?
- currently rerunning some mach stem simulations for longer time to see if they reach a steady state with mach reflection

### This Week

- finish rerunning mach stems
- add aluminum cooling table that Francisco gave me and run his jet set up…will probably need kraken for this
- start running 1D radiative shocks with magnetic fields using parameters that Pat sent

# Meeting Update Dec.3

Did the VISA thing last week and passed. The internet connection here is fairly good so I can work as usual. I'll be back after about 4 weeks.

**Triggered Star Formation: Cloud Stability** Did the Boss setup stability test with pressure matching on the boundary. I changed xi so that the central density and the ambient density can match their number exactly (6.2e-19 and 3.6e-22). One thing I found is that using this setup, our code gives a different isothermal temperature (124K) comparing to theirs (10K). Also, at the clump radius (0.058 pc), the density jump seems to be 1000 times rather than 100 times. The following movie is for a run to about 1.3e+5 years, which is the supposed shock crossing time scale (1e+5 yrs). The sound crossing time is 1e+6 yrs, so we need to wait a bit longer to see the long term evolution.

http://www.pas.rochester.edu/~shuleli/collapse/bss1.gif

**Resistive Paper** Been trying to write up the introduction section. I'm currently looking into previous literature so that I can correctly place the paper's goal.

**Clump Paper** Will be sending Adam and Eric the final revision this week.

HEDLA paper was accepted last week.

Registered the Cancun meeting in February. I used the previous clump paper's abstraction.

# Meeting Update - Jason and Ivan 1202

Runs on Kraken are complete. We are only running the 5 Jupiter mass for a little bit longer to see it disrupt.

I begun analyzing evolution of disk mass over time and sent a few plots to Jason.

We are currently deciding what criteria to use to determine what material is part of the disk.

# Meeting update

Not too much to post on the blog this week as I have mostly been learning from Toro this past week, and spent some time working on BE paper stuff. I am in chapter 4 now of Toro, and it seems like I will be writing up the exact Riemann solver code this week.