Developed a solver for determining which branch of the local Bondi flow kernel edge corresponds to, and integrates the equations inward, and pastes the soln on the grid.

See the different possible soln classes in attached PDF.

To test the solver have a problem module that specifies an outer boundary condition of the flow (constant values at r≥1), which is internal to the grid.

For this boundary, the code figures out which branch to integrate the equations on and pastes the steady state soln within r<1.

I let the sim run for 20 frames using this set up (holding the boundary values of the problem fixed for r>1), using the Krumholz accretion module (see attached rho.gif for a movie, and the pdf for the initial conditions).

I need to finalize some things with the new sink accretion algorithm to utilize this framework. Things are mostly built with that so shouldnt be much longer. Need to come up with test problems.

​http://www.icehap.chiba-u.jp/IAU_B1/resources.html

New Work

1. Drag force/secondary mass paper:
   1. Movies showing contours of component of force per unit volume exerted on secondary…
      1. along \phi direction
      2. along direction of \bm{v}_2-\bm{v}_1 (=drag if net force is negative and =thrust if net force is positive)

Movies

Movies are in the frame of reference co-orbiting with the particles, with the secondary at the center.

Face-on density with contours showing magnitude of force/volume exerted by gas on particle 2 along direction of velocity of particle 2 relative to particle 1 \bm{v}_2-\bm{v}_1:
​Density with drag force/volume contours (Run 143: No subgrid accretion, M2 = 1 Msun)

Face-on density with contours showing magnitude of force/volume exerted by gas on particle 2 along \phi direction:
​Density with drag force/volume contours (Run 143: No subgrid accretion, M2 = 1 Msun)

Future work to improve the above movies

1. Add vector to show direction of force exerted on particle 2 by the gas \bm{f}_\mathrm{2,gas}.
2. Add vector to show \phi-component of force exerted on particle 2 by the gas f_{\mathrm{2,gas,}\phi}.
3. Add vector to show direction of velocity of particle 2 with respect to particle 1 \bm{v}_2-\bm{v}_1.
4. Add vector to show (\bm{v}_2-\bm{v}_1)-component of force exerted on particle 2 by the gas \bm{f}_\mathrm{2,gas}\cdot(\bm{v}_2-\bm{v}_1)/|\bm{v}_2-\bm{v}_1|.

(adding these is somewhat non-trivial in VisIt so haven't quite gotten to it yet.)

Qual

• What should I be reviewing?
• What should I focus on in the brief (10 journal-style pages)?
  • Should I try to work in all of the papers, tell the story so far in limited detail? Focus on just one?

​Radiation pressure paper

• Sent an email around with information related to our previous discussion. I've begun incorporating some of the information where appropriate. Body of email copied below for reference.

• Could use some brainstorming on alternate ways to make the same point as these measurements:
  • Comparing the average speed of a particle being significantly affected by radiation pressure to the escape speed of the system
  • F_{L-\alpha}, F_{th}, F_{ram} comparison (could also compare pressures, if I can figure out a good length scale to convert volumetric radiation force to radiation pressure)

Even using the high-memory nodes on the visual partition, VisIt either crashes (at end of trying to load Chombo) or hangs (while processing) when trying to make the necessary plot. Possible to load all of the VTKs created when saving 3D plot in python, maybe?

• Cherenkov et al 2017 ionization procedure:
  • This is the entirety of what they say about their ionization procedure: “In the grey approximation, the last term in equation (6), the rate of photoionization, can be written in the form , where h is the Planck constant, σUV = 6.3 × 10−18 cm2 is the photoionization cross-section at the Lyman limit frequency ν0, FUV = F0 (apl2/l2) is the intensity of ionizing radiation at distance l from the star centre, F0 = 884 erg cm−2 s−1 is the radiation intensity at the planet orbital distance, apl, (Schneiter et al. 2016) and τ is the optical depth.”
  • Sections 2.3 and 2.4 talk about the radiation transfer model they use, but only in the context of the radiation pressure (Lyman-alpha line). My assumption would be that they intend for this to extend to the ionizing radiation, but they never make this explicit.
  • None of the previous studies that could reasonably be considered to be by their group use self-consistent ionization
• Launch conditions from Bourrier et al.
  • Bourrier et al 2013 launch metaparticles of neutral hydrogen only, and stop tracking ionized hydrogen—this may have a significant effect on the dynamics, since much of what is being pushed on (indirectly) is ionized
  • Their highest density is at launch radius, and is ~107 /cm3. Around the same location (Mach surface/Hill sphere at ~2.8 Rp), we have a total density of ~6x106 /cm3, with a neutral density of ~8.5x104 /cm3 (for no-flux simulation). So their optical depth to Lyman-alpha radiation is significantly higher.
  • They say “the number of atoms in a meta-particle is calculated to keep a low value of the corresponding optical depth dτ in front of a fraction of the occulted stellar surface.” Unclear to me what exactly this means.
  • Their metaparticles are launched randomly in every direction with thermal speed of 11000 K (much higher than our temperature at this point).
• HD 209458 vs. other stars' Lyman-alpha fluxes
  • From Wood 2005, Lyman-alpha fluxes for other G0 V stars normalized to HD209458b orbital distance, in erg/s/cm2: 18344, 18955, 14374, 23836; HD209458 = 6409
  • So it’s on the low side, by a factor of ~3
• Optical depth in wings (near planet)
  • At 3 Rp in +y, 50 km/s ± 2km/s bin, optical depth of 1.66 (in high flux case)
• What happens as we transition between ionizing flux regimes?
  • From Sec 5.3 of John's paper: While in E-limited, neutral fraction constant, density higher, so total neutral density (and therefore optical depth) higher. Crossing into R-limited, neutral fraction reduces as 1/F1/2, and density increases at the inverse rate, so optical depth becomes constant. At some point, wind may start turning around to create the swirly arms – at this point, optical depth should ~double from previous.
• Local ionization equilibrium
  • Essentially maintained at 0.03 events/sec throughout wind, outside of small region - see attached VisIt plots (​high side, ​high top, ​med side, ​med top, which plot the absolute value of RecombinationRate-IonizationRate (so they're a measure of the distance from ionization equilibrium)
• Effect of averaging over bubbling in high-flux case
  • Very small effect, mostly deeper center transit in the phase where the wind is expanding out of the Roche lobe. See ​attenuation and ​observation comparisons.
• Plots of velocity, neutral density
  • ​Final attached plot is neutral mass density with velocity contours plotted over. Contours are: blue, 1x10^{5 cm/s; red, 5x10}5 cm/s; cyan, 1x10^{6 cm/s; yellow, 5x10}6 cm/s; magenta (none present), 1x10^{7 cm/s}

New Work

1. Jets paper: wrote problem module for phase 1 (no subgrid accretion), compiled on bluehive, but run error. At this point it would be wise to merge my version of astrobear with the version that implements jet feedback.
2. Drag force/secondary mass paper:
   1. Made movies for runs 149 and 151 (M2 = 0.5 Msun and M2 = 0.25 Msun)
   2. Wrote script to compute drag force
   3. Ran script to compute drag force on reduced resolution run 143 and produced various plots (see below)
   4. Started draft of ​paper

Drag force plots for Run 143 (de-resolved version to enable faster analysis using VisIt)

1) ​Force exerted by gas on particles, calculated numerically
Description: Force \bm{f}_i\cdot\bm{v} and components of \bm{f}_i for both particles, with inter-particle separation curve plotted in grey for reference. Velocity components for each particle are also plotted for reference (with an arbitrary linear scale). Note that positive r-component means away from the other particle, while positive \phi-component means in the sense of the orbit.
Comments: Force on secondary in frame of primary \bm{f}_2\cdot\bm{v}/v is shown by black line. Phi-component f_{2,\phi} is shown by orange line. Note that at early times, \bm{f}_2\cdot\bm{v} is dominated by f_{2,r}v_r since f_{2,r} is so large. At later times the orange and black lines coincide, so \bm{f}_2\cdot\bm{v} is dominated by f_{2,\phi} v_\phi.

We see that the \phi-component of the force (orange) is positive during plunge-in (so of opposite sign to the predicted drag force). The secondary is being accelerated around in its orbit by the posterior side of the envelope, which lags the primary particle in its orbit (paper 2). Subsequently, the force is mostly a drag force (-ve \phi component) but the \phi component actually oscillates from positive to negative .

2) ​Force exerted by gas on particles, calculated semi-analytically
Description: Force is now calculated using formula, but with velocity inputted from simulation. Black line is drag force on secondary in the frame of the primary. Density and sound speed are equal to the values in the initial RGB profile, at radius equal to the current inter-particle separation a(t). Density, relative speed and inter-particle separation are plotted for reference (with arbitrary linear scales; density and speed increase with time while separation decreases).

Comments: The drag force on the secondary in the frame of the primary is predicted to be higher in magnitude when a is smaller and density and speed are larger. The low density is predicted to cause the drag force to be negligible at early times.

3) ​Force exerted by gas on particles, numerical normalized by semi-analytical
Description: Numerical solution normalized by semi-analytical solution for a few choices of formulae.
Comments: At late times, the magnitude of the drag force is of order a few per cent of the predicted value.

4) ​Bondi radius
Description: Various definitions of Bondi radius plotted for both particles. Also plotted for reference are the inter-particle separation (grey) and pressure scale height of the original RGB profile at r=a(t).
Comments: The most relevant Bondi radius is probably the one represented by the thick solid red curve, corresponding to the Bondi radius around the secondary in the frame of the primary, including the sound speed in the denominator. However, for all definitions, a uniform medium is assumed. That is, the density and pressure gradients are neglected. Also, the initial sound speed of the original profile is assumed. Even more problematic, the envelope is assumed to be stationary in the frame in which the Bondi radius is computed (lab frame or frame of primary point particle). In any case, looking at the thick red line, we see that R_\mathrm{2,BH}(t) is comparable to a(t) and comparable to the pressure scale height H. Thus, applying the Bondi-Hoyle-Lyttleton formalism to this case is highly questionable. However, it remains to beseen whether applying this formalism would be more justified as M_2 is reduced (in the limit M_2\rightarrow0, R_\mathrm{BH}\rightarrow0 since v\sim (\mathrm{G}M_1/a)^{1/2}, so independent of M_2, in that limit. This can be tested with Runs 149 and 151 which have M_2 equal to ½ and ¼ of the value in the fiducial run (Run 143) plotted here.

Updated Movie Library for Runs 132, 143, 149, 151

Density (zoomed in) in lab frame

​Face-on density (Run 143: No subgrid accretion, M2 = 1 Msun)
​Face-on density (Run 149: No subgrid accretion, M2 = 0.5 Msun)
​Face-on density (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Movies corresponding to figures in Paper 1

Movies are in the reference frame corotating about the secondary with the instantaneous orbital angular speed of the particles, and with the secondary at the center.

Figure 1/Figure 2—-face-on density:
​Face-on density (Run 143: No subgrid accretion, M2 = 1 Msun)

​Face-on density (Run 132: Subgrid accretion, M2 = 1 Msun)

​Face-on density (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Face-on density (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Figure 4 top panel—-edge-on density:
​Edge-on density (Run 143: No subgrid accretion, M2 = 1 Msun)

​Edge-on density (Run 132: Subgrid accretion, M2 = 1 Msun)

​Edge-on density (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Edge-on density (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Figure 4 bottom panel—-edge-on density, zoomed in:
​Edge-on density (zoomed in) (Run 143: No subgrid accretion, M2 = 1 Msun)

​Edge-on density (zoomed in) (Run 132: Subgrid accretion, M2 = 1 Msun)

​Edge-on density (zoomed in) (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Edge-on density (zoomed in) (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Figure 6—-flow around companion:
​Tangential velocity with velocity vectors (Run 143: No subgrid accretion, M2 = 1 Msun)

​Tangential velocity with velocity vectors (Run 132: Subgrid accretion, M2 = 1 Msun)

​Tangential velocity with velocity vectors (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Tangential velocity with velocity vectors (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Movies corresponding to figures in paper 2

Movies are in the lab (~system CM) reference frame with the CM of the particles located at the center of the frame.

Figure 3 top row—-face-on normalized gas binding energy (red means unbound, blue means bound, yellow is density contours, vectors are velocity):
​Face-on normalized energy (Run 143: No subgrid accretion, M2 = 1 Msun)

​Face-on normalized energy (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Face-on normalized energy (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Figure 3 second from top row—-face-on normalized gas kinetic energy (magenta means thermal energy dominates, green means bulk KE dominates, yellow is density contours, vectors are velocity):
​Face-on normalized kinetic energy (Run 143: No subgrid accretion, M2 = 1 Msun)

​Face-on normalized kinetic energy (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Face-on normalized kinetic energy (Run 151: No subgrid accretion, M2 = 0.25 Msun)

Extra Movies

​Temperature (Run 143: No subgrid accretion, M2 = 1 Msun)

​Temperature (Run 132: Subgrid accretion, M2 = 1 Msun)

​Temperature (Run 149: No subgrid accretion, M2 = 0.5 Msun)

​Temperature (Run 151: No subgrid accretion, M2 = 0.25 Msun)

​Sound speed (Run 143: No subgrid accretion, M2 = 1 Msun)

​Sound speed (Run 132: Subgrid accretion, M2 = 1 Msun)

​Mach, lab frame (Run 143: No subgrid accretion, M2 = 1 Msun)

​Mach, lab frame (Run 132: Subgrid accretion, M2 = 1 Msun)

​Mach, frame corotating about secondary (Run 143: No subgrid accretion, M2 = 1 Msun)

​Mach, frame corotating about secondary (Run 132: Subgrid accretion, M2 = 1 Msun)

Side-by-side comparison of Model A/143 (left) and Model B/132 (right) from Paper 1

​Face-on density
​Edge-on density
​Edge-on density, zoomed
​Tangential velocity
​Temperature
​Sound speed
​Mach in lab frame
​Mach in frame corotating about secondary

Questions for further analysis

• What is the relative velocity between the particle and gas in its vicinity?
  • Could make movies of relative gas velocity (had this a while back)
• What is the relative contribution to the gas-particle force as a function of distance from the particle?
  • Could plot net gas-particle force including only gas in spheres of different radii around the particle
• What does the 2D map of gas-particle force look like for each particle?
  • Could make movies with vectors for the particle-gas force at each point

Next steps

• Energy analysis for runs 149 and 151, including 5 figures:
  • Energy components vs time plots (Fig 1a and 1b of Paper 2) for runs 149 and 151 (Yisheng)
  • Mass unbinding vs time plots (Fig 2 of Paper 2) for runs 149 and 151 (Yisheng)
  • Particle CM — Envelope CM relative motion (Fig 6a and 6b of Paper 2) for runs 149 and 151 (Yisheng)
• Stability analysis comparing runs 143 (without damping run) and 132 (with damping run) initial conditions, to go into appendix of Paper 3 (Yisheng)
• Drag force computation for Run 143 (full resolution), Run 149 and Run 151 (may need post-processing)
• Jet test runs (need to merge versions of astrobear)

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

AMR line transfer

Making good progress. I believe I have the implementation almost completely defined, just remains to write some of the organizational code. The algorithm I have isn't as efficient as I'd like, but I haven't yet figured out a way around the limitations. Once I have a working prototype (end of the week?) it may be worth reviewing to see if there are any speedups or reduced memory usage possible.

In some searching, I also found what appears to be a ​decent spherical algorithm for ray tracing. It's an approximation method combining the long-characteristic and short-characteristic methods that was implemented for FLASH. I definitely think it's something we should consider for AstroBEAR (I haven't tried yet, of course, but I believe it's possible to implement it in our code).

Radiation pressure paper

Have looked at Luca's comments. The biggest is ​a paper by Cherenkov et al., who have done something similar but included the Doppler shift in the absorption coefficient. I can't tell how they calculate the optical depth for ionizing radiation. They also include the stellar wind and start with a much higher-temperature (~2x) and higher-density (~100x) planet. Their resolution is significantly smaller, though it changes in a way that's unclear (I believe they just use unevenly-spaced points) from the outer portion to the inner portion of the grid.

There are also a couple of comments on clarity that I could use another set of eyes on.

Finally, the measurements that I'm attempting to make keep crashing VisIt. I'm going to try yet another method this week.

​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

Project Statuses

• Radiation Pressure
  • ​Paper sent around to everyone (astrobear Overleaf account still on free plan)
  • Luca's looked at it and made some small edits and some comments — haven't looked at any of them yet
  • Two measurements still to make (I've commented at the location)
    • Radiation pressure vs ram pressure vs thermal pressure at edge of wind (need to think about this one some more)
    • Average speed of material absorbing significant amounts of flux (currently using between 100% and 1% of maximum radiation pressure as boundary of "significant")
• AMR line transfer
  • Sketch of code written, significant modifications required
  • Main focus for next couple of weeks
• Charge exchange
  • Will run just stellar wind, no ionizing flux, to see effects of just mixing (can run this week)
• HD209458b
  • Need to discuss at some point whether we want to do a paper on just the planetary wind
  • High ionizing flux - effect on mass loss rate/wind structure

Charge Exchange