New fiducial run with accretion onto the companion, run021 (J4)

click on each image to play a .gif file

Jet tracer density rendered in pseudocolor with ramp transparency

The labels 0 and 1 mark the center of primary and secondary stars respectively. (Please ignore the time slider in the movies.)

Left: Day 0-14, in lab frame looking down towards the orbital plane. Middle: Day 15-40, in lab frame looking down towards the orbital plane (notice the scale is smaller compare to the left one) Right: Day 0-40, in lab frame side view looking at the x-z plane (this view includes the simulation's outer boundary as the gray frame around it.).

Pseudocolor jet tracer density, compare two cases with accretion (J4 and J5)

Viewed in slices with perspective from the primary to the secondary. The green circle marks softening radius around the companion particle.

Left: J4, the new fiducial model with Jet Mdot = 2e-3 Msun/yr. Middle: J5, 10 times the Jet Mdot, other setups are the same as J4. Both simulations have the companion accreting from their environment.

Right: J1, Jet without accretion onto the companion.

Jets reappear at late time

Compare the time evolution in jet density (top row) and total gas density (bottom row). The view is from core of the primary looking at the companion, zoomed around the particle (softening radius marked by the green circle for reference).

Left: J4, P2 accretes from the environment, Jet Mdot = 2e-3 Msun/yr Middle: J5, P2 accretes from the environment, Jet Mdot = 20e-3 Msun/yr Right: J1, no accretion, P2 loses mass to supply the jets.

A repository of plots for all the runs

​Click_this_link

Common Envelope Physics and Outcomes (CEPO 2021)

​https://phsites.technion.ac.il/common-envelope-physics-and-outcomes/

A virtual meeting, August 30 - September 3, 2021 (13:00-19:00 GMT, 9:00-15:00 EST)

I (Amy) submitted abstract on the jets in CE simulations.

Asymmetrical Post-main-sequence Nebulae 8 e2021 (APN8)

The Shaping of Stellar Outflows

​https://sites.google.com/view/apn8-e2021

A virtual meeting, October 4-8, 2021

Luke put in abstract on accretion and jets during CE evolution, also including the recent AGB simulations, collimation of post-AGB outflows and planets in CE.

​http://www.pas.rochester.edu/~yzou5/CE/Blogpost_20210712/

see attached files

Current jet setup

• Problem: jets are launched inside the particle softening radius, and most of the material is bound relative to particle 2.

• Increasing jet radial speed by a factor of 2 (the green lines below):

Energy relative to P2, with smoothed gravity

Energy relative to P2, with Newtonian gravity

Notes about the Federrath+14 jet/outflow model:

​http://www.pas.rochester.edu/~yzou5/CE/notes_Federrath+14_jet_model.pdf

Jet parameters in our model:

to-do's

1. understand the feedback module. what are the initial profiles of density and velocity inside the spherical cones launching the jets?

2. compute how much of the jet material is unbound, using the spline potential, velocity profile and density profile.

3. plot the jet tracer density

​http://www.pas.rochester.edu/~yzou5/CE/Feedback_Tracer/

Attached to this blog post is my notes on ​Krumholz+2007

Bound vs. unbound energy with density contours:

Frames in each row are chosen as:

• 1st periastron
• 1st apastron
• 2nd periastron
• ~ 30 days (periastron for runs 6, 10, 11; and apastron for runs 13, 14)
• ~ 35 days (periastron)

Run 06, normal jet

Run 10, super jet

Run 11, no jet

Run 13, no jet, small P2

Run 14, normal jet, small P2

Jet

Determining mass-accretion and jet mass-loss rates in post-AGB binary systems

​Bollen+2020

In this work, Bollen et al. modeled two observed post-AGB binary systems with jet associated with the companion (MS). They determined the jet mass-loss rate in the order of 1e-7 to 1e-6 Msun/yr, and the mass accretion rate onto the companion disk is between 1e-6 to 1e-4 Msun/yr. They found that re-accretion from the circumbinary disk should be the major source supplying the accretion disk and the jet.

Ranch

Have access to 2TB of storage on ranch, but the 5 production runs of the jet simulations take a total of 27.4TB. How much quota do we have? Should I backup every other frame instead of all?

Radiation Transfer Test

I tried to run the RadTransferTest module. Got run time error with the photon-ionizing line transfer, which writes 'dt is NaN'. I ran the module with line transfer turned off (iLineTransfter=0).

Conferences coming up

1. APN8-Asymmetrical Post-Main-Sequence Nebulae 8, the Shaping of the Circumstellar Matter

​http://apn8.iaa.csic.es/

October 4-8, 2021, Granada, Spain.

2. IAU Symposium 366-The Origin of Outflows in Evolved Stars

​https://iaus366.be

November 1-5, 2021, Leuven, Belgium

This one is postponed from 2020. I will check if they still accept abstracts.

Jet simulations

1. How does the morphology of the jet (outside the launch region) evolve?

• Density slice in orbital plane ​movies-face-on

TBD: parallel to the orbital plane; Add contours to show percent of jet material.

• Comparison of fiducial and superset runs to see the difference directly

​movies-edge-on

• Slice through the secondary that is perpendicular to the orbital plane (view from P1), superjet run ​superjet-movie (adding velocity vector somehow messed up the particle contours)

2. How sensitive are the properties to the parameters: jet Mdot, secondary mass?

• Compare superjet and normal jet orbital evolution and explain the differences using the other analysis.

• How important are real jets (how realistic is the superjet?).

• Effect of smaller mass companion.

3. How does the presence of the jet affect the orbital evolution, and why?

• Two competing effects: jet blow off effect at the beginning and increasing drag later?

• If necessary, plot force density and drag-force vs time as in the force paper (Luke's paper 3).

4. How does the jet affect the envelope unbinding?

• In-progress: time evolution of bound and unbound ENVELOPE mass

(but also look at jet material to see whether it stays unbound or becomes bound).

• TBD: spatial map of unbound vs bound mass.

Plan for this week besides analyzing the jet runs:

• Do a few runs with the radiation module.

• Gather references for the jet paper.

Forming a qual committee

• I've got all five people confirmed to be on the committee. They are Adam, Eric, Segev, Miki, and Rip. I will check in with them for scheduling (also need to work with Laura and Jean). Probably a week in February/March.

CEJet

• all runs:

• Super jet runs:

• Normal jet runs:

• Small P2 runs:

• Plot and analyze the data this week.

• What happened with the small companion runs?

​Plots for the new jet runs

CEJet

• Now the tracers work correctly tracking the primary envelope, ambient, and jet material.
• Do the three remaining runs this week.
• Plot orbital and energy info when bluehive comes back.

Rad hydro - reading

• ​Krumholz+2007
• ​Ensman 1994

What do we have for the CEJet project

Runs:

Code:

• I'm testing the tracer on ambient. otherwise the astrobear_1008 version is good to go for the remaining runs.

code version	note
_1014_C	particle buffer on both P1 and P2
_1212	use Shape to define refinement region
_1008	updated to fix tracer on jet, good for stp_run_006 & _009
_0112	refine on density criteria

Analysis:

• Orbits, average separation plot. ——> working on the script
• Mass, bind vs. unbind, find a best way to present this. ——→ working on script to read the totals data file
  • maybe something like Fig. 8 in the AGB paper, where more mass gets elevated but still not unbind
• Energy, use data from the totals, do post-process if necessary
  • maybe analyze drag force as well, but orbit plot also shows the effect of drag force
• Present figures similar to the AGB paper
  • Fig. 6, 7, 8
• Update the separation plot
  • use CE run 152 for high mass fiducial, which doesn't change the softening radius

From thermal bomb meeting, 2020-7-16

Emily's new PN project:

• Inject for a brief period but strong energy into the outflow
• over free-fall / orbital time scale
• how much energy does it take to destroy the torus?
• calculate reasonable parameters for the CE binary used in my PN simulation (Emily)
• launch the wind with new parameters in the current PN module for ~6 days (free-fall timescale) (Amy)

○ pressure support for the central region after shutting down the outflow?

GEE reading notes

Soker (2020) "Shaping PNe with jets and the GEE"

General properties of jets in shaping PNe

1. Vjet ~ V_esc, MS V_esc ~ 600 km/s —> V_jet ~ (0.5~2) V_esc = 300 ~ 1200 km/s

2. Jets deposit energy + momentum by heating/expelling ambient —> less mass to accrete

negative feedback cycle —> quench/re-start —> many ejection episodes

3. "Pressure release valve" —> high accretion rates

positive feedback cycle

4. might have a wide open angle —> theta~90 degrees

GEE: companion grazes the giant envelope (~photosphere)

• BHL accretion of surrounding gas
• RLOF from inner dense envelope

ILOT (intermediate-luminosity optical transients)

(​https://phsites.technion.ac.il/soker/ilot-club/)

t_photon_diffusion .le. t_gas_expansion —→ large amount of E_th —→ Radiation lost

GEE forming SNe IIb

• SNe IIb = core collapse, almost no H
• jets remove H in the giant envelop during GEE
• jets are on the MS companion, the giant will end in core collapse SN

Five phases of Jets' launching —- companion launches the jets

1. Jets in wind accretion zone
   • Hillel+2020
   • P2 has jet before enter P1 envelope
   • BHL accretion from slow wind
   • Jets being dragged/bent behind, can't penetrate the accreting wind
2. Jets from pre-CEE RLOF accretion
   • equatorial accretion flow, RLOF mass gain much greater than that from BHL
   • jets penetrate (likely) the wind, form bipolar lobes/point-symmetric
   • can also occur earlier in wind-RLOF
     • Harpaz+1992
     • Mohamed+Podsiadlowski 2007

3. Jets during the GEE (too delicate??!)
   • Eta Carinae?
   • P2 accretes near periastron
   • BHL+RLOF
   • Shiber+2017

4. Jets during the CEE
   • BHL from inside the envelope
     • Shiber+2019
     • jets don't break out the envelope
   • what is jets' role in helping mass removal?

5. Jets during the post-CEE
   • P2 accretes from circumbinary disk
     • Soker2019, Kashi+Soker2011, Chen+Podsiodlowski2017
   • P2 launches jets w/ variable directions + intensities

​EAS 2020 Interactive Schedule S8 is the common envelope section.

Luke's talk on Thursday, July 2, at 11 am Eastern Time ​abstract

Amy's talk on Friday, July 3, at 4:15 am Eastern Time ​abstract

Relevant papers

​Krumholz+2007

​Leenaarts 2020 review paper

Radiation transfer and flux limited diffusion in AstroBEAR

​AstroBEAR wiki page reproduced

I've checked format up to section 2.2 only. Later sections may contain typo and mistakes.

My note for today's meeting

​My note on Leenaarts section 2

​My note on Marshak wave (following Castor chapter 12.1)

CEJet

​Notes

Run time on Stampede2

Frame rate for the jet modal on 120 cores:

Time for one modal with 120 framesis ~170 hours in wall time = 20 400 SUs.

I have 3 runs planned for this week, which will use about 81 200 SUs.

PN talk for EAS2020

Presentation time is 15 min, including a 4-5 min Q&A and 1 min for transition. I will prepare 10 slides max for the talk, and maybe <2 extra slides for questions.

Outline:

1. Introduction to the problem. Show ~2 examples of bipolar PPN with collimated jets.

2. Simulation setup (1): initial condition from the SPH model, choice of physical parameters.

3. Simulation setup (2): technical details in the AMR simulations: relax the ejecta on AMR grid, fallback during the quiescent period, regions of refinement, and application of the cooling curve.

4. Results (1) - model A, the central pillar of jets and inner shock

5. GISW model and the elliptical shock

6. Results (2) - model B and model C

7. Results (3) - model D, compare to OH 231.8+04.2

8. Asymmetries in line plots

9. Momentum budget

10. Conclusions

Meeting on rad hydro

Thursday (6/18) at 11 am. Zoom coordinates to be emailed.

Revision

Revised paper is ready to submit again. ​The last version

Jet project - plan for this week

Runs:

1. Rerun the super jet model, fix tracer, add tracer for ambient
2. Rerun the no-jet fiducial, M2 = 1 Msun, use my jet modual, this is like CE run 152 (or 143, but don't change softening radius), use the initial orbit conditions for the jet runs, add tracer for ambient
3. Half-solar mass companian, fiducial, no jet run
4. (Maybe) Restart, add jet from middle of the fiducials

Analysis:

• Orbits, average separation plot. see the eccentricity formula form Luke's blog post
• Mass, bind vs. unbind, find a best way to present this.
  • maybe something like Fig. 8 in the AGB paper, where more mass gets elevated but still not unbind
• Energy, use data from the totals, do post-process if necessary
  • maybe analyze drag force as well, but orbit plot also shows the effect of drag force
• Present figures similar to the AGB paper
  • Fig. 6, 7, 8
• Update the separation plot
  • use CE run 152 for high mass fiducial, which doesn't change the softening radius

Calculations:

• how much energy do we pump into the simulation (jet)?
• Does the jet material fall back? (tracer plot will also help)

Extending the CEPN simulations

​Link to a summary document

Starting the Rad Hydro

​my plan to start

List of new runs

1. Run like model A, but has rsoft=5, rwind=4 (production runs have rsoft = 5, rwind=6). -stopped at frame 32, too slow. mintracer=1d-7
   • Another run with rsoft = rwind = 5 to capture the central mass. Also, set mintracer=.1 for quiescent period.
2. Run like model A, but without the quiescent phase. then compare the momentum analysis for both. — queue on bluehive.
3. Run like model A, but with larger refinement on level 7, restart on frame 149. — waiting for reservation on bluehive.

Major issues

Issue 1.

• In section 2, paragraph 3. explain the SPH-AMR mapping. Whether it is conservative or not. —basically averaging and interpolating by some standard sph practice. Add inputs from Thomas.

Issue 2.

• Modify section 2, discuss the choice of momentum parameters. explain that the momentum was chosen to be in line with Bujaharrabal.

• Discuss how sensitive are the results to the choice of injection radius.

— New run 1. but had issue with central high density cells.

— also compare the outflow starting time, New run 2.

Issue 3.

• max refinement level was never reduced during the simulation.
• boxy structure: is always leading ahead of the central cylinder of refinement.
• New run 3, doubling refinement radius on level 7.

Issue 4.

• Clearer description about 'tracer' in Astrobear.

The referee was confused about tracer as added material, but it's just a label, or a 'dye'.

Issue 5.

More discussion on the appendix in the main text. The horizontally launched jets enhance the shaping argument.

Issue 6.

• Add input from Bruce.
• Produce density projection figures.

Minor issues

• In table 2, added a line to indicate model names (A,B,C,D) to the corresponding run parameters.
• Fixed typo throughout the paper.

-TBD-

• Replace the text labels in figures with larger font.
• Redo figure 6. Experiment with other color-tables and vector properties. Just very hard to interpret. Try use arrows of variable length on a uniform mesh of points.

​http://www.pas.rochester.edu/~yzou5/CE/Blogpost_20200302/

• Backed up data and code versions from Stampede2 to Bluehive.
• Extended model 6 (half-mass companion) to frame 146.

See detailed note here ​http://www.pas.rochester.edu/~yzou5/CE/Blogpost_20200210/Jet_update.pdf

Link to movies shown last time ​http://www.pas.rochester.edu/~yzou5/CE/Movies_CEJet_20200131/

Jet runs and status:

model A, run002, at 46 frames

Starting with P2 at the edge of the RG envelope.

Gets to the green dot below in the binary separation plot. The next job is in queue on stampede.

Evolution movies: ​rho-edge-on ​(with-v-vec) ​rho-face-on ​(with-v-vec) ​Temperature ​z-speed

​Run002-model-A-Movies-and-plots

model C, run 003, at 11 frames

Starting with P2 at 1.5 times RG radius.

I might pause the run and fix the following problems first:

• RG gets unstable around the edge of its envelope.
• Jets' temperature has a sharp change aligned with the boundary between level 3 to level 4 refinement area. Seems like pressure protection, but the actual temperature is 1000 times higher (~1e-6 K) than the minTemp set in physics.data (1e-10 K).
• May be a bug in the particle buffer code I haven't found out.

Evolution movies: ​rho-edge-on ​(with-v-vec) ​rho-face-on ​(with-v-vec) ​Temperature ​z-speed

​Run002-model-B-Movies-and-plots

Two production runs on Stampede

Run 002, model A

Jet-mass-loss-rate = 2 Msun/year, higher ambient profile (rho0 = 6.69e-9, p0 = 1e5), a0 = 49 Rsun.

Stopped at frame 28, waiting for the next queue.

Run 003, model C

Jet-mass-loss-rate = 2 Msun/year, lower ambient profile (rho0 = 1e-10, p0 = 2e3), a0 = 73.5 Rsun.

Currently running, has got to frame 5.

• Each run is on 16 nodes (768 cores). Recent queue time on Stampede is about 3.5 days, so I submit jobs in a 3-day separation.

• lSkipProfile = False for production quality. This takes time at every restart, (dt~13 hours between the first two chombos files). Later on each frames take about 2 - 1.5 hours, and is getting faster throughout the run time (so I didn't use the restarting at each 5 frames scheme).

• The last run on Stampede was killed with the following message:

any ideas?

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 334217 RUNNING AT c498-002
=   EXIT CODE: 7
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
   Intel(R) MPI Library troubleshooting guide:
      https://software.intel.com/node/561764
===================================================================================
TACC:  MPI job exited with code: 7
TACC:  Shutdown complete. Exiting.
Starting replacement routine
/tmp/slurmd/job4506883/slurm_script: line 69: .true.: syntax error: operand expected (error token is ".true.")

Summary of the four models planned

• A: 2 Msun/yr jet, higher ambient, starting at 49 Rsun.
• B: 2 Msun/yr jet, lower ambient, starting at 49 Rsun.
• C: 2 Msun/yr jet, lower ambient, starting at 73.5 Rsun.
• D: no jet, lower ambient, starting at 73.5 Rsun.

​density-vz-and-Temperature-comparison

CE Jet - Stampede run 1

Setups

• 73 frames done.
• Mass-loss-rate = 2 Msun/year
• Jet_radius = 2.25 Rsun (16 cells)
• working on restart on Stampede

Separation compares to fiducial run

Movies

​Density-edgeon

​Density-faceon

​Compare-Jet-fiducial-edgeon

​Compare-Jet-fiducial-faceon

Things to do

questions:

• how does the jet evolve?

○ quenched? ○ when, and why ○ look at tracer, zoom in of the jet, force balance, ram pressure (balance with gravity and thermal pressure)

• does the ambient matter?

○ compare run with lower ambient (another run)

• does the launch time matter?

○ start from farther out (another run) ○ and a fiducial (no jet) run to compare

• what happens at later times?

runs

1. continue this run, 2 Msun/yr jet
2. a lower ambient run
3. larger separation, with low ambient
4. no jet, larger separation, low ambient

• AGB?

PN paper

Paper is coming along. ​paper I spent some time fixing the calculations for section 4.3. ​notes Need to finish the conclusion.

Movies

​production runs

​tests_without_IW

CEJet

Current run on Stampede

In queue, expected start time is tomorrow morning at 6:30.

Same as the 1st run, but with all the total quantities defined as the current CE settings. Use refinement radius as CE model B in paper II, 4d12 for frames 0-45.

Parameters of the Jet

• mass loss

• launching radius

• speeds

• others

—-

PN

Flows during initial wind period

The total amount of gas flow into the wind radius is 0.19~M_{\odot} . This amount of gas has been removed from the simulation. The inflow through the outer boundaries is almost negligible (1e-7 solar mass).

This means almost half the initial gas (0.4 solar mass) has been removed. There is a stream of inflow through the funnel with a slow speed, but since it’s at the region of highest density, it takes a lot mass in.

​velocity plots

​density plots

​python notebook

​test without IW

The SPH CEE simulation contains

• Primary of mass 0.88 M_sun, where
  • the core mass is 0.392 M_sun, and
  • the envelope is 0.488 M_sun
• Secondary (point particle) of mass 0.6 M_sun
• so the amount of gas for us to start the PN simulation is about the envelope mass of the SPH primary

Frame 0 has

• total mass of the gas about 16 solar mass — M_tot =16.136 M_sun
• if set a density threshold at 1e-8 g/cc, then
  • the high density region sum up to — M_high=15.712 M_sun, while
  • the low density region sum up to — M_low = 0.424 M_sun

NOTE

1. The density threshold is chosen as such, because since frame 1, the maximum density in the simulation is on the order of 1e-8 g/cc.
2. The total mass of the low density gas is about right to the total amount of gas we start with.
3. The high density region is produced from mapping the SPH data onto the AMR grids. The SPH data has local maximum density on the order of 0.1 g/cc, while frame 0 has local maximum of 0.05 g/cc.
4. The mass of the high density gas is confined at a point close to the simulation center with a radius of 7.6 R_sun (about the size of our finest grid). These 4 extremely high density grids are within the wind initialization radius, so once the simulation starts they are removed, and won't affect the simulation any further.
5. I interprete this 16 M_sun region as the final location of the binary core from the SPH simulation. Maybe the core is more confined in the SPH simulation but spreads out in AMR grids and produce this artifect of a huge amount of mass.

Between Frame 0 and Frame 1, there are few sources and sinks of mass:

1. Sink - Lowering the floor density for density protection, the maximum effect has the change of -1.4e-3 M_sun
2. Source-Sink - carving out the central region and replace with a 1 solar mass point particle and fill the rest of the region with initial wind density.

From Frame 1 to Frame 149 (before fast wind starts):

1. Source - inflow from the outer boundary of the simulation. During the entire initial wind period (3000 days), the total amount of inflow gas is 1.11E-7 M_sun, almost negligible
2. Sink - gas flows into the wind initilization radius will be removed from the simulation, and the total change is -0.19 M_sun

CEJet

From last week: 3 test runs planned

4 levels of AMR with 64 base grid resolution. Each run is planed for 10 days in simulation time, and takes about 2 days to complete on bluehive.

Run 45 - low ambient, small separation

Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use lower ambient density and pressure profile: rho_amb = 1e-10 g/cc, p_amb = 1961 dyn/cm2 .

​Movies and frames

​density-edgeon-throughparticles

​density-faceon-box

Run 46 - high ambient, small separation

Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use ambient profile as the CE fiducial run (run 143): rho amb = 6.67e-9 g/cc, p_amb = 1e5 dyn/cm2

This and the first run will compare the effect of ambient material over the jet. Also a comparison with the CE fiducial run.

​Movies and frames

​density-edgeon-throughparticles

​density-faceon-box

​Temp-edgeon-throughparticles

Run 47 - low ambient, large separation

Start with the secondary farther out: a0 = 98 R_sun (meant to do 72.5 R_sun, but I messed up), and use lower ambient profile.

​Movies and franes

​density-edgeon-throughparticles_frame16

​density-faceon-box

Issue: unstable primary, 2 tests without jet

Secondary starts right against the primary envelope: a0 = 49 R_sun. Use lower ambient density and pressure profile: rho_amb = 1e-10 g/cc, p_amb = 1961 dyn/cm2 .

test 1 - low ambient, small separation

no mass in the secondary (M2=0),

but the two particles are given the same initial velocity as the orbiting case.

​Movies and frames

See the primary blows up in the opposite direction of its motion: ​density-faceon

test 2 - low ambient, small separation

with secondary (M2=0.978 Msun),

this is like Run 45 without the jet.

​Movies and frames

Same effect observed in the orbiting case: ​density-faceon

a zoomed-in view: ​density-faceon-zoom

CEJet

Working on branch ticket_442.

Questions to answer for jet project:

• How does the jet affect the morphology of the envelope?
• How does the jet affect the ejection of the envelope?
• How does the jet evolve, does it get quenched?
• What is the dependence of these questions on when the jet gets turned on?
• What is the dependence on accretion rate?
• Opening angle?

movies from test run 44

Starting with secondary right against the primary envelope. No secondary mass. Jet mass loss rate at 2 M_sun/year, and radial velocity 430 km/s (Keplerian).

​density-edge-on ​temperature-edge-on ​density-face-on

3 test runs planned for this week:

4 levels of AMR with 64 base grid resolution. Each run is planed for 10 days in simulation time, and takes about 2 days to complete on bluehive.

1. [Run 45] Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use lower ambient density and pressure profile: rho_amb = 1e-10 g/cc, p_amb = 1961 dyn/cm^{2 .}

2. [Run 46] Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use ambient profile as the CE fiducial run (run 143): rho amb = 6.67e-9 g/cc, p_amb = 1e5 dyn/cm²

This and the first run will compare the effect of ambient material over the jet. Also a comparison with the CE fiducial run.

3. [Run 47] Start with the secondary farther out: a0 = 72.5 R_sun, and use lower ambient profile.

There is also a CE run without jet to be compared.

Some parameters to explore

• Initial separation (a0) (ref. ​Shiber+2019)

49 R_sun, right against the primary envelope. Jet might be quenched after the secondary enters the envelope.

72.5 R_sun, start farther out, longer time for the jet to evolve (although farther means lower mass loss rate if we use accretion rate to determine the mass loss of the jet).

109 R_sun, Roche distance of our system.

• Jet mass loss rate

From CE paper 1, the secondary can have extremely high accretion rates of ~ 0.2 - 2 M_sun/year. This rate corresponds to 100 - 1000 times the Eddington rate for a MS star, or 1e4 - 1e5 times for a WD.

Currently we are using the upper limit (2 M_sun/year). For a more physical value, we plan to compare between 1 vs 100 times Eddington rate of a WD (100 times WD rate would be ~ Eddington rate for a MS).

Check out the Xsede proposal for more suggested tests.

• Jet radial velocity

​Federrath+2014 suggests using Keplerian velocity at the surface of the star.

In the current test runs, this is computed using initial mass of the secondary (0.978 M_sun), and 1 solar radius. So,

jet_vrad = sprt (G*M2/R_sun) = 430.75 km/s

• Jet radius, distance from the secondary particle within which the jet is initialized. 64 cells (70.4 R_sun) for now. Federrath+2014 also tested this effect, and suggests using 16 or 32 cells.

• Refinement radius (compare to the orbit). See CE paper 1 fig. 3, and paper 2 appendix fig. C1.

• Refinement shapes and qTolerance.

Jet in x-direction

This is mostly to address a comment from Eric, if we put in an equatorial outflow, will it still be collimated?

Instead a "disk" outflow, I put a jet in the +- x-directions, with 15-degree open angle, otherwise the same as model A.

​http://www.pas.rochester.edu/~yzou5/PN_xjet_figures.pdf

On momentum (L/c)

In our simulation, the "momentum rate" (L/c) is set by the outflow's density and velocity

L_{mech}/c = 4\pi r_w^2 \rho v_w^3 / c

for high-momentum flux (models A and C), (L/c)_high = 3.835e+34 (g cm s^{-1 year}-1), and

for low-momentum flux (models B and D), (L/c)_high = 1.917e+33 (g cm s^{-1 year}-1)

At the last frame presented in PN paper, the total momentum ejected is then:

P_{sim} = t_{dur}\cdot (L/c)

Observationally, the upper limit on t comes from the size and jet speed: t<R/v, and maximum momentum is estimated as:

P_{measure}< (R/v_w)\cdot(L/c)

Here I use the extend of the top lobe for R, and wind speed v=300 km/s, the momentum calculated in those two ways and their ratio are:

Initial Mass Flow

This seems problematic. I will check my Visit code and come back with another update. Ignore the following for now.

Outer boundary flows inward, floor density raise remove some mass, but both are small.

total change: -0.2 Msun

outer BC: 1.1e-7 Msun

floor raise: -1.4e-3 Msun

Left: including frame 0; Right: excluding frame 0.

PN paper

draft updated on 05/05: ​figures ​paper draft

CE Jet

​note

Last time we got the jet turned on, and added tracer to it.

Should look at the initialization of temperature, and measure the boundary flows.

Eric suggested starting jet from RLOF.

Others

• Start working on the Toro book.
• Radiation transport in CE.
• MESA, and plan for summer travel.

Links to a working draft of the paper:

​figures

​paper draft

​references

on overleaf:

​https://www.overleaf.com/read/tcqgypntptny

PN paper

• Prepared list of references and graphs to show. Writing sections for simulation setup and results. Share on overleaf?
• Forwarded last blog post to Orsola, waiting for her comments.
• High momentum runs are getting slow on Bluehive: 4 frames / 5 days using 120 cores. Maybe end the simulation at different times?
• Currently we have 220 frames (4400 days) for high momentum runs, and 300 frames (6000 days) for low momentum (need another five days for the non-cooling case to get there).
• Movies and graphs from the production runs are now password protected on my PAS server, what I learnt from Yisheng's instruction, same setup as he did.

CEJet

• Made scripts to calculate jet mass loss and boundary outflow. It seems the jet is not turned on, but keep in mind that the recent test runs have only a few frames each, so actually can't make any conclusion from this calculation.
• to-do: add tracer onto the jet.

A quick glance at the four models — density evolution

panels are not time synchronized

• Top left: Model 1, non-cooling, high-momentum
• Top right: Model 2, non-cooling, low-momentum
• Bottom left: Model 3, DM-cooling, high-momentum
• Bottom right: Model 4, DM-cooling, low-momentum (notice the scaling difference)

Numerical setup: ​link to Overleaf write-up

Some take-away points

1. Outflow is highly supersonic, and creates layers of shock.

​Example: Mach number in model 1

2. Outflow is highly collimated. Most of the outflow materials are redirected into the top and bottom lobes.

​Example: Outflow tracer and mesh in model 3

3. Top-bottom asymmetry is most dramatically shown in model 4, see the bottom right panel above (or click to show in full-screen).

Movies and plots from the four models:

​everything

Table of index

Naming rule decoded:

• "ADB" = "adiabatic", no extra cooling in the code
• "DMC" = "DM-cooling", cooling curve applied to T > 30000 K
• "_High"/"_Low" = "high-momentum outflow" or "low-momentum outflow"
• For each model, I plotted the following variables:
  1. Mach number with contours for selected Mach number
  2. Mesh and outflow tracer density
  3. Temperature and contour at 10000 K
  4. Density, labeled as "rho"
  5. Density at the orbital disk, labeled as "zrho"
• Numerals indicate certain zoom-in options:
  0. the entire simulation domain
  1. core area, focused on the CEE ejecta
  2. center, slightly zoomed out from the core
  3. lobes in the early stage
  4. tall, extended lobes in later stage

Current work on PN and CEJet

1. PN runs: restarting non-cooling this week, cooling low-momentum finished at 6000 days.
2. Analyse run time as a function of jets extend (TBD)
3. Write up a first draft paper (TBD)
4. (With Yisheng) build modulated code doing all the visit plots and movies. Can be use for both PN and CE simulations (in progress)
5. Analyse boundary mass loss rate in the jet runs (TBD)

Discuss the CEJet paper, Shilber et al. 2019. (see attachment)

PNe

• I have a new reservation on BlueHive started this morning at 8. Currently running the two cooling models: High-momentum at frame 189/300, and low momentum at frame 231/300.
• Need to extend the time for those cooling case to evolve longer. The morphology hasn't grow far from the core area, and in low-momentum case, the jet hasn't propagate through the funnel yet.
• The non-cooling, high-momentum ran into a resolution problem at frame 212. Redoing those frames now.
• Non-cooling, low-momentum case, currently at frame 259/300. I'm thinking of redo some later frames to keep the high resolution area cover the entire PN lobes.

CEJet

• Did a test run with no secondary mass (run008-amy).
• I made script for measuring total mass and mass-loss rate. in run008, the rate mass being added is not a constant.
• Schedule a paper discussion next week.

others

• the MESA summer school

​non-cooling

​cooling plots

PNe Simulation

• Production runs: initial wind runs to 168/300 frames. Estimated wall time remaining is 10 days, but I've doubled the cores (120) and should take less than that.
• Need to reserve nodes on Bluehive for the rest (high-/low-momentum, adiabatic and cooling for each).
• Low-momentum, refinement issue solved by reducing mintracer to 0.0001 (0.1 originally)
• cooling, issue with tracer density when triggering temperature protection. Should be a quick fix. ​plots
• Inner boundary conditions, look for zero pressure gradient at the inner radius where outflows start.

Moving on

• Begin working on Toro's book
• Pick up recombination

ProtectOutfowTemp subroutine added on Jan 31

SUBROUTINE ProtectOutflowTemp(pos,t,dt,q,Ghost)
  USE GlobalDeclarations
  REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
  REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
  REAL(KIND=qPREC), INTENT(IN) :: t, dt
  LOGICAL, INTENT(IN) :: Ghost
  if (q(particle%outflowobj%itracer) > 1e-4*particle%outflowobj%density) then
     q(iE)=max(q(1)*FloorTemp/TempScale, q(iE))
  end if
END SUBROUTINE ProtectOutflowTemp

A lot movies to show, with links to individual frames. Three runs in comparison:

1. ADB-run003, no cooling, fast wind starts at 3000 days (frame 150).
2. DMC-run001, restart from ADB-run003 frame 149, use DMCool for T > 10000 K. A problem is that outflow temperature was set to floor temperature (10 K).
3. DMC-run002, currently running this version as Jonathan helped me modify source_control.f90, but the temperature problem still presents.

​Links to all movies and frames, breaks down below:

Must see

Comparison between non-cooling (ADB003, on the left) and cooling (DMC001, on the right), showing side by side

Movies on individual variables

Temperature Line Out

left: along x-axis, right: along z-axis (need to fix x-scale)

• Frame 149

• non-cooling
  • frame 150

• Frame 151

• cooling (DMC002)
  • frame 150

• Frame 151

​Link to current simulation movies

From last time:

​Run003-no_cooling-rho ​Run003-no_cooling-Temp

​Gamma=1-rho ​Gamma=1-Temp

Status

• Three runs over the break, each runs for 5 days. See ​figures below.
• Problem with tracer not being placed onto last positions upon restart.
• Next step, add cooling function and turn on self-gravity. Also as a respond to Garcia-Segura's paper.
• Some found reference for journal club talk: ​https://www.overleaf.com/read/yrwnmdcfjmhd

4-level, bigbox

tracer problem upon restart

7-level AMR

gamma = 1

PNe in big box

Use the updated problem module, 4 level AMR, and a cubic box with 128 kR_sun radius (64000 in computational unit, base grid 1000 R_sun). I did two runs over the past week. The first one has minTemp=10K, and mindensity=1d-10 for all. In the second one, I raised those two value to doubled for frame 0, and then lower them back for the rest of the simulation.

• Links to movies and graphs:

​Nov29 ​Dec2

• "Disk" — lower density along equatorial plane
• "cold bubble" in temperature plots

PNe

I'm running the updated problem module from Jonathan. Many thanks.

Initial wind phase, 7-level AMR. Plots of density (g/cc), temperature (K) and Mach number. ​frames

From last meeting with Adam, Luke and Zhuo, we came up some options to explore:

• gamma = 1 (still use the ideal gas EoS?)
• widen box — make it a cube, run on lower resolution
• apply DM cooling where the shocked gas temperature is greater than 10000 K
• wide/narrow angle jets
• precessing jet — need Nozzle module

Things to look into:

• opacity in the CE ejecta vs. the out-flowing wind.
• as temperature goes up and density goes done, what happens to opacity?
• metal cooling and other cooling curves in the code

Planetary Nebula Simulation

1. Protection issue
   • MinDensity = 1d-100 in density protection
   • Regenerated frame 0, running the initial wind phase
   • Now temperature drops below the MinTemp(=1d-10) in pressure protection
   • See ​plots

Mach number (did something wrong?)

2. Cooling (Gamma = 1)
   • Need to learn terms in a(5,128) and b(5,192) matrices in problem modules.
   • Read in quantities from SPH output and mapped to Info%q(i,j,k,1-5). See quoted code below.
   • a,b matrices: 1-density, 2-v_x, 3-v_y, 4-v_z, 5-not sure
   • Info%q: 1-density, 2-p_x, 3-p_y, 4-p_z, 5-energy density
   • Missing gamma?

select case(info%level)
    case(0)
        do k=1,128
           do j=1,128
              read(195,*),a(1,:)
              read(196,*),a(2,:)
              read(197,*),a(3,:)
              read(198,*),a(4,:)
              read(199,*),a(5,:)
              if (j.ge.57.and.j.le.72) then
                  do i=1,16
                      Info%q(i,j-56,k,1)=a(1,56+i)/rscale
                      Info%q(i,j-56,k,2)=a(1,56+i)*a(2,56+i)/rscale/velscale
                      Info%q(i,j-56,k,3)=a(1,56+i)*a(3,56+i)/rscale/velscale
                      Info%q(i,j-56,k,4)=a(1,56+i)*a(4,56+i)/rscale/velscale
                      Info%q(i,j-56,k,5)=a(1,56+i)*(a(5,56+i)+0.5*(a(2,56+i)&
                      *a(2,56+i)+a(3,56+i)*a(3,56+i)+a(4,56+i)*a(4,56+i)))/pscale
                  end do
              end if
           end do
        end do

$\gamma_{7_1} = \frac{5}{3},\,\, \gamma_{7_2} = 1\\ e_{T1}=\gamma_{7_1} P_1\\ P_1 = \frac{e_{T1}}{\gamma_{7_1}}=P_2=P\\ e_{T2} = e_{T1} \frac{\gamma_{7_2}}{\gamma_{7_1}} = \gamma_{7_2}P\\ E = e_{k1}+e_{T2}$

3. AMR
   • Trying on different qTolerance

Planetary Nebula Simulation

Current status and things to be done

1. Tuning AMR
   • Refinements triggered on density gradient.
   • Restarted from high-momentum run02 frame with jets out.
   • Ran 5-level AMR (figure below). However, the outflow was turned off.
   • Current run starts from the same frame with 7-level AMR, and outflow turned on. The job is running on 120 nodes with 20-gb memory per node. It has been running for 4 hours (as I'm writing this blog) to generate 1 frame. Need better tuning on qTolerence, or reduce AMR level.
   • Next, start with no jets and then turn on the outflow. See if AMR capture the jets as we want.

2. Resolution on cavity walls
   • Issues on pressure and density protection. (tbd)
   • Turn on outflow after initial-wind fills the entire simulation domain.

3. Gamma = 1 (Cooling)
   • Set gamma_7 = 1.001 in the code. (where exactly?)
   • Generate Frame0: Read in SPH results (pressure, kinetic and thermal energy, gamma = 5/3), keep pressure fixed, calculate new thermal energy with gamma = 1.
     $\gamma_{7_1} = \frac{5}{3},\,\, \gamma_{7_2} = 1\\ e_{T1}=\gamma_{7_1} P_1\\ P_1 = \frac{e_{T1}}{\gamma_{7_1}}=P_2=P\\ e_{T2} = e_{T1} \frac{\gamma_{7_2}}{\gamma_{7_1}} = \gamma_{7_2}P\\ E = e_{k1}+e_{T2}$
   • Criteria, say, apply cooling if v > 100 km/s ?

level 5 mesh on 5-level AMR

temperature

Taken at the beginning of the simulation and just before/after turning on the fast wind.

High momentum run 01	High momentum run 02	Initial Wind
InterpOpts = 0	InterpOpts = 1	InterpOpts = 1 , and no fast wind

​Movies and ​Graphics

7-level AMR run returns the following error

bhc0074.bluehive.circ.private:UCM:2a62:11a41700: 520618407 us(520618407 us!!!): dapl async_event: DEV ERR 10
bhc0074.bluehive.circ.private:UCM:2a5a:3db57700: 521043336 us(419038 us): dapl async_event: DEV ERR 9

density and pressure protections in physics.data

!================================================================================
! Density Protection Options
!================================================================================
lRestartOnDensityProtections = .false.  ! Do density protections trigger restarts?  [.false.]
lTrackDensityProtections = .false.      ! Track density protections [.false.]
iDensityProtect  = 0                    ! 0-Set to MinDensity, 1-Set to minimum nearby density, [2]-Average nearby densities
iMomentumProtect = 1                    ! 0-Conserve momentum, 1-Set to zero, [2]-Average nearby velocities
MinDensity       = 1d-15                ! Minimum computational density before protection is triggered [1d-10]

!================================================================================
! Pressure Protection Options
!================================================================================
lRestartOnPressureProtections = .false. ! Do pressure protections trigger restarts? [.false.]
lTrackPressureProtections = .false.     ! Track pressure protections [.false.]
iPressureProtect = 0                    ! 0-Set to MinTemp, 1-Set to minimim nearby pressure, [2]-Set to average nearby pressure, 3-Set to minimum nearby temperature, 4-Set to average nearby temperature
MinTemp          = 1d-10                ! [1d-10] minimum allowed temperature for the system in Kelvin before protection is triggered

!================================================================================

PNe

Temperature

• last time, made mistake with Temperature, should be scaled by 9975.67K according to scale.data. Once it's fixed, the low temperature is 1d-11 to 1d-10 K, which is closer to the set MinTemp. I will fix and update other temperature movies
• ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/High_momentum_run02_T_inK.gif

Initial wind

• density ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Initial_wind_D10120_rho.gif
• temperature ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Initial_wind_D10120_T.gif

Recombination

Luke helped me to check the equation of states, the gamma factor is correct. we also draw BCs from Luke's CE simulation as initial state for the Rec equation solver. still working on it. this is the old note: ​https://www.overleaf.com/read/xsrfvtjfcyzs

Directory for all the movies

• ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/

• Table 2: links to individual movie

• Note: Interpolations in Run-02 are set to 1, and 0 in Run-01

PNe

• Will update links to output once they are available

High momentum wind

• Run-02: ran to 1465 days post quiescent phase

Rarefaction issue

• Run-02: a snapshot from high-momentum Run-02 with interpolations = 1.
• Notice the rarefaction on the edge of a square "box" with side = 8 (unit = 1000 R_sun). It's also seen in the previous simulation by Zhuo (see low-momentum case in the pdf)

• Run-01: Here is a comparison between snapshots from the conference proceeding and those reproduced using Run-01 data. ​http://www.pas.rochester.edu/~yzou5/PNe_compare-snapshot.pdf

Low momentum wind

• Reran to 2600 days post quiescent phase

Recombination

• Adding recombination into the isothermal solver (not physically correct, but just a tester)
• Need to understand the roles of gamma in EoS and energy equation. I'm reading MacLeod2017 paper ​https://arxiv.org/abs/1704.02372 and Hansen's text book on the topic.
• ​https://www.overleaf.com/read/xsrfvtjfcyzs

PNe

Rerunning high-momentum wind with interpolations = 1 (Run-02)

Quiescent phase (6000 days from start) is done, continue to see if fast wind phase produces the same rarefactions between cell boundary. Zhuo also suggested removing point mass at center.

Here is a comparison between snapshots from the conference proceeding and those reproduced using Run-01 data.
​http://www.pas.rochester.edu/~yzou5/PNe_compare-snapshot.pdf

• a snapshot from high-momentum Run-02 with interpolations = 1. Notice the rarefaction on the edge of a square "box" with side = 8 (unit = 1000 R_sun). It's also seen in the previous simulation by Zhuo (see low-momentum case in the pdf)

Recombination

​https://www.overleaf.com/read/xsrfvtjfcyzs

To-do list

• recombination equations (priority)
• sort out high momentum wind result
• continue running low momentum wind

Low Momentum Wind

Speed = 300 km/s, density 5e-13 g/cm³

Run to 7400 days (1400 days after quiescent phase) Density movie:
​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/Low_momentum_Day7400_rho.gif

Next, adjust refinement region and level to continue.

High Momentum Wind

Run to 8000 days (2000 days after quiescent phase). Will update a link to the results as soon as they are done.

Maybe redo the runs for better consistency in frame rate and refinement?

self-gravity

To do list:

• sharp density gradient
• try to make the Poisson solver work properly
• ambient density

ratio between Laplacian of potential and density

wind

• started low momentum wind
• high momentum wind extended to 7500 days

Movies

Outflow

370 days after quiescent phase ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/High_momentum_Day6370_rho.gif

Quiescent phase

6000 days ​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/High_momentum_D6000_Qphase.gif

Refinement region

for frames 380-392 (t = 6370 days). it takes ~6 hours to generate 1 frame. dimensions are in solar radius.

level 1-4

(x, y, z) = +-(5000, 5000, 10000)

level 5-7

(x, y, z) = +-(1000, 1000, 6000)

PNe simulation

x-z density profile at 6010.9, 6150, 6196.43, 6250 days.
​http://www.pas.rochester.edu/~yzou5/PNe_simulation_movies/Old_movies/High_momentum_Day6250_rho_y-slice.gif

• ran simulation to 6250 days (250 days after quiescent phase), no self-gravity
• next: adjust refinement region

Recombination

• still working on the equations

isothermal flow and Bondi problem:

​https://www.wolframcloud.com/objects/yxamyz/isothermal-fixed_x.nb

Progress and questions

1. Set up flow equations in Mathematica
2. Try to solve simplified cases:
   • Isothermal, with fixed x=0 or x=1
   • Adiabatic cooling
   Luke suggested there may be analytical solutions for these situations.
3. How to calculate the recombination coefficients?

​equations ​test03 ​isothermal

May be relevant to our work

PDFs in the attachment

• ​Goldfinch1971

• ​Goldfinch&Pack1971

• [Chou&Talbot1967]

• [Gol'dfarb&Kostygova&Luk'yanov1969]

other references are listed in the attached note from Luke