Planetary magnetic fields are dipole fields of the form

\frac{r_0^3}{r^5} \left ( 3\mathbf{r}\left ( \mathbf{B_0} \cdot \mathbf{r} \right) - \mathbf{B_0} r^2 \right )

or

B_{\theta}=-B_0 \cos{\theta}\frac{r_0^3}{r^3}

and

B_r=2 B_0 \sin{\theta}\frac{r_0^3}{r^3}

which (if we assume that the dipole is oriented along z) gives a field of

2 B_0 \hat{z}

at the pole and a field of

-B_0 \hat{z}

at the equator.

This also corresponds to a vector potential

\mathbf{A}=r_0^3 \frac{\mathbf{B_0} \times \mathbf{r}}{r^3}

Now if we assume that

\mathbf{B_0} = B_0\hat{z}

Then our vector potential is

\mathbf{A}=r_0^3 \frac{B_0 \sin \theta}{r^2} \hat{\phi}

Also, the dipole field used in Matsakos was a constant field inside of the planetary core (r<r_p/2)

They quote a constant field of 16 B_0 \hat{z} at the center would require a potential of the form

\mathbf{A}=\frac{16 B_0 r \sin{\theta}}{2} \hat{\phi}

which agrees with the potential outside at r=r_0/2

So we don't need any smoothing… just a piecewise definition for the potential that switches at r=r_0/2

Now, if we want to truncate the field lines so they do not pass out the inflowing boundaries, we need to have the derivative of the potential go to zero. This is tricky since the curl of the potential needs to go to zero in both the theta and r directions. However, the following potential has the desired properties

A_\phi^*=r_0^3\frac{B_0 \sin{\theta}}{r^2} \left (1 -\frac{r}{r_s} \right ) ^2

For r > r_s the potential is just 0.

This gives

B_{\theta}^*=\frac{1}{r}\frac{\partial}{\partial r}\left ( r A_\phi \right ) = -B_0 \sin{\theta}\frac{r_0^3}{r^3} \left ( 1 - \frac{r^2}{r_s^2} \right )

and

B_r^*=\frac{1}{r \sin{\theta}}\frac{\partial}{\partial \theta}\left ( \sin{\theta} A_\phi \right ) = 2B_0 \cos{\theta}\frac{r_0^3}{r^3} \left ( 1 - \frac{r}{r_s} \right )^2

for r < r_s and \mathbf{B}=\mathbf{0} for r>r_s

	Smooth	UnSmoothed
log
linear

Here is an outline of the ideas I am thinking for the ring analysis. The main idea is that the shocked material of the flows acquires some new density and velocity which translates into a ram pressure outward into the ambient medium. We can assume that when this material meets the ambient of equal pressure, the material will come to rest and that is where the ring will form. An issue is — we need a radial dependance of the pressure in the ambient (or from the splashing), so that we can find a radius where these 2 pressures balance.

Using just the shock conditions means we consider the splashed material having constant rho*v². So what are the physical mechanisms that will provide a radial dependance on the pressure?

Here were some of my ideas (by the way, ignoring MHD for now):

-The splashed material slows down as it is being ejected from the collision region due to gravity — perhaps refer to Parker wind solutions (or something similar to Joe's analysis…, parker wind might not be the right term here for what I am thinking, but haven't looked into this yet… )

-The ambient medium is infalling toward the gravitational potential minima, i.e. the center of the box. This infall has some radial behavior —- and if I can back out the radial density and velocity profiles, will have a ram pressure to try and balance to the shock ram pressure

-Then more complicated stuff, like compression of the outgoing splashed material as it rams into the ambient gas, leading to cooling/loss of KE/and a concomitant decrease in velocity.

I focused on the 2nd case as it was the one I was most familiar. Although, it might not be entirely the right regime to do this analysis, depending on whether the ambient is fully jeans unstable. I thought it was a nice place to start, to jump into the equations, and the process of building the model.

Now, I am nearly done for the Hydro, nornmal shock situation, assuming no reflection. As you can see on this page, ways I can modify this model for the head on case:

• making the shock reflected (expect small deviations from the values)
• making the collapse not spherical but rather cylindrical (also don't expect much deviation)
• adding magnetic pressure to the mix… probably here would get the greatest - change in the values… i.e. my model might fail the most because it doesn't include magnetic pressures at this point.

Here is a page describing the shock conditions I am using:

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

Here is a page describing the uniform collapse model I have been developing:

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

Uniform collapse involved:

• understanding the solutions for uniform collapse again
• couldn't find this easily on the web, so dug into the solutions myself
• coded it up in mathematica to visualise the solutions
• got the model nearly where it needs to be to begin using numbers from my sims to check it out

As there was no great description of this that I could find, either in text books or the web, I thought it would be good to write up what is in my notebook on the wiki. Need to try and get this visible in the google universe, as a search for uniform collapse sphere, or homologous collapse sphere, isn't bringing it up on google yet. Will check quickly if the visibility can be increased.

As far as papers, I have been reading astroph everyday, reading the highlights in astro on popular science sites, and have been reading the star formation newsletter.

I applied to Astrohack week but didn't get in the program. I have also been working through the equations and there importance in the Fed. paper on outflows.

XSEDE Proposal Renew

1. highest priority. Due July 15th.
2. ​Current documents

OutflowWind module

1. 3D co-rotating frame: fixed some bugs and latest results— blog:bliu06222015
2. Park wind solution for the star: merged with Jonathan's setup and running a test. Will post results soon

pnStudy

1. Merged with current development branch for Eddie's cooling stuff. Compiled. Running test
2. will do 3D runs with the new merge version on bluehive and install the new code on Spain machines

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

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

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

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

A recent movie​

• Astronum and my talk went very well. I got a lot of good compliments, and talked with a lot of people.
  • Amira vs. Disperse paper with Federrath?

• Finished more 2-D Mach stem runs for Pat and the paper: ehansen06222015

The complete set of models result in the following figure:

• Work continues on generating all the emission maps for my 3-D clump runs. Here is the updated table: MachStems
  • I have a reservation on BlueStreak right now to finish up run E. It should be done by next Monday. Then, I just need to finish up the emission maps for runs D and E.
  • For next week, I will start generating movies of all the emission maps that I do have.

• My goal for July is to finish the 3-D clump paper and start some new 3-D pulsed jet simulations.

Here is a matlab routine to solve for the parker wind given lambda x is in units of planet radius, and y is velocity in units of sound speed.

Also notified Titos about typo in arxiv preprint that made initial attempts to solve the equation nonsensical.

Also note that the original version of this post assumed that

\psi=v/c_s

when in fact

\psi=v^2/c_s^2

So here is the correct version of the matlab routine

lambda=11.5;
N=81;
decs=4;
xx=logspace(0,decs,N);
yy=zeros(N,1);
syms y positive;
for i=1:N
  yy(i)=vpasolve(y - log(y) == -3 - 4*log(lambda/2)+4*log(xx(i))+2*lambda/xx(i),y,(xx(i) > lambda/2) * 100+.001);
end
fid=fopen(['parker_wind_',num2str(lambda),'.dat'],'w')
fprintf(fid, '%i\n',N);
fprintf(fid,'radius, density, velocity\n')
for i=1:N
    fprintf(fid,'%e %e %e\n',xx(i),1d0/xx(i)^2/sqrt(yy(i)/yy(1)), sqrt(yy(i)));
end
fclose(fid);

1. bug and fix

The code of co-rotating frame for results in blog:bliu06222015 and the results before hard-coded the position of particle and pointgravity object as (0,0,0). Since the origin point in the co-rotating frame should be the mass center and the planet position is at (200,0,0) (currently,will make more sophisticated). So this is a bug:

IF(.NOT. lRestart) THEN
      CALL CreateParticle(PlanetParticle)
      PlanetParticle%q(1)=planet_mass
      PlanetParticle%xloc=0
      PlanetParticle%radius=radius
     CALL CreatePointGravityObject(PlanetParticle%PointGravityObj)
      PlanetParticle%lFixed=.true.
      PointGravityObj=>PlanetParticle%PointGravityObj
      PointGravityObj%mass=planet_mass
      PointGravityObj%x0=PlanetParticle%xloc

The fix is set the PlanetParticle position to be Outflow%position —

      PlanetParticle%xloc=Outflow%position

2. Low-Res Results after the fix

After the bug-fix, the low-res results with omega=0 looks promising. The tight sonic surface problem in blog:bliu06222015 seems gone. But it's still not quite right for the gravity/mach number plot inside the planet —trying with higher res.

;​low res density move 5-zone-per-radii	;​low res temperature move 5-zone-per-radii
	​log mach plot

plot of log|vx|
shows the unmatched gravity plot for 2D and 3D seems from the low-res of 3D plot.

3. 3D corotating with omega=0.5

​movie for 8 or 16 zones per radii

This is to try to understand/fix the tight sonic surface found in 3D co-rotating frame with omega=0 as shown in the 3rd part of blog:bliu06112015_2 and the different line plots for the 3D co-rotating and 2.5D as shown in the 4th part of the same blog post… In 2.5D and 3D simulations before the problem module hardcoded PLUMERSOFT and soft radius=1 for the pointgravity.

        SPLINESOFT = 1, & !g ~ r/r^2 for r < r_soft and then goes to 0 as r-> 0 
        PLUMMERSOFT = 2 !g ~ (r^2+r_soft^2)^(-3/2) r

1. Code modification: Splinesoft

Since in the 3D case we use r=1 but in 2.5D we use r=2. So the g in 3D is about twice of the g in 2.5D.

To handle point gravity more properly, we change the code to use SPLINESOFT instead. Here's the new update to the OutflowWind module

        outflow_radius=0 ! outflow radius is 0 always
        outflow_thickness=planet_radius 
        pointgravity_soft=1 !splinesoft always
        pointgravity_r_soft=0.5*planet_radius

2. Testing Results with Splinesoft =

	New Results with SplineSoft (low res)	Old Results with PlummerSoft (high res)
2.5D no stellar wind	;​movie 8 zones per radii	;​high res movie
2.5D stellar wind	;​movie 8 zones per radii	;​high res movie
3D Co-rotating	; ​density 5 to 10 zones per radii	; ​density with 40 to 80 zones per radii
3D Co-rot and 2D plots compare	;

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

1. Model Description

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

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

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

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

The temperature of the simulation is fixed at 100 K.

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

2. Simulation

3AUsideview​

3AUtopview​

​4AUsideview

4AUtopview​

6AUsideview​

​6AUtopview

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

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

3. Conclusion

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

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

Wrote this python script that generates points on a sphere in Cartesian coordinates. It then writes these coordinates with some random values with velocity (vx, vy, vz) and density (n) to a tab-delimited ascii file. You can choose the radius of the sphere, and N is number of points on the sphere. Forwarded from Martin, ​a stackoverflow discussion. Also I am using python3.

import random
from math import pi,sin,cos

#Creating Sphere Dataz                                                                                                                                                                                      
def createSphere(r=5, N=100):
    lst = []
    for phi in [(pi*i)/(N-1) for i in range(N)]:
        M = int(sin(phi)*(N-1))+1
        for theta in [(2*pi*i)/M for i in range(M)]:
            x = r * sin(phi) * cos(theta)
            y = r * sin(phi) * sin(theta)
            z = r * cos(phi)
            lst.append((x, y, z))
    return lst

#Opens/creates new ascii file                                                                                                                                                                               
outfile = open("test.dat", "w")

#Writes the data to the file                                                                                                                                                                                
for x,y,z in createSphere():
    rho = random.random()*1000000
    vx = random.random()*10
    vy = random.random()*100
    vz = random.random()
    print("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}".format(x,y,z,vx,vy,vz,rho), file=outfile)

#Closes the file                                                                                                                                             
outfile.close()

So the pseudo-sphere data seems to visualize just fine in Shape:

3D Module

Render Module

Aside: I am starting to document how to use Shape, ​here.

small window with new parameters:

Tiny stellar wind density to check the planetary wind close to the planet surface.

​movie for 40 zones per radii

2. No rotation omega=0

​density omega=0; 20 zones per radii

​temperature omega=0; 20 zones per radii

​density omega=0; lambda=10

​temperature omega=0; lambda=10

3. No rotation omega=0 High Res longer

Set stellar wind with very low density and run with high resolution. To check if it's possible to reproduce the planetary wind in 2.5D.

​density omega=0; 40-80 zones per radii
​Temperature omega=0; 40-80 zones per radii
​pressure omega=0; 40-80 zones per radii

4. Line plots of density, temperature, Mach number and GM/R^2on dayside (from left boundary to critical radius)

Same run data as 3 above.

3D:

gamma=1.01
3dpressure=(gamma-1)*(E-0.5*(px^2+py^2+pz^2)/rho)
vScale=sqrt(vx^2+vy^2+vz^2)
soundSpeed=sqrt(gamma*3dpressure/rho)
mach=vScale/soundSpeed

2D:

gamma=1.01
2dpressure=(gamma-1)*(E-0.5*(px^2+py^2)/rho)
soundSpeed=sqrt(gamma*2dpressure/rho)
mach=sqrt(vx^2+vy^2)/2dsoundSpeed

change the density to be 10^{-6} smaller to make the density flux on the planet surface to be a reasonable value.

​8 zones per radii

[​8 zones per radii tiny omega

In the latest high-res results of co-rotating frame as in blog:bliu06032015, the planetary wind seems not be able to expand. This test to test the effect of tiny corotating speed omega=1d-5 . All other parameters are same as blog:bliu06032015, namely rho_sw=4d-1…

​density of 8 zones per radii;

​temperature of 8 zones per radii

Column Density Map (m2)	PV Diagram (km/s vs. m)
​GIF	​GIF

The simulation I am trying to visualize in Shape is probably Martin's old poloidal collimation project. Above is a visualization of frame 22 (chombo files cannot be found, so I cannot compare in VisIt), so the beginning of the simulation. Toward the end of the simulation, it should supposedly look like this:

VisIt

Shape (probably not frame 212)

Objectives

• Still think the PV Diagram could use improvement. Not sure if I used the appropriate parameters for the velocity axes. For the above visualization, the axis goes from -500 to 500 km/s. With this arrangement I got something to show in the PV diagram part of the render module. However I think I could mess with it more.
• Document how to do this on the wiki!
• Regrid Bruce's simulations so they are in ascii, and thus fed into Shape.
• Visualize Bruce's simulations.
• Show Bruce how to use Shape.
• Pretty up the output?

The setwind subroutine In the current co-rotating frame sets up the stellar wind density proportional to 1/r^{2}. where r is about 200 C.U. for our runs So to make the stellar wind density same as in the 2.5D and 3D we have to set it 40000 times larger in these co-rotating runs. Otherwise the stand out distance and bow-shock size will be very large as we see in our former runs…

By using larger stellar wind density, the result for \gamma =1.01, Mach 5 runs seems promising. — didn't change the density of planet for this run still 1g/cc also used larger planet as in blog:bliu06022015

Movies

​Density movie -- 16 zones per radius

​Temperature movie -- 16 zones per radius

​Density movie -- 32 zones per radius and twice higher stellar wind density

​Temperature movie -- 32 zones per radius and twice higher stellar wind density

Double Rp and Mp to keep the lambda same. New refinement setup. Here's the results of a 2.5D testing run with gamma=1.01, lambda=5.3, Rp=4, Mp=1Mj, mach=5 and max 3 levels of AMR

Mesh		​mesh movie;​stellar wind mesh movie
Planet wind		​zoomed-in movie to check sonic surface;​stellar wind movie