Parker Wind

Differential Equation

The Parker solution to the solar wind assumes a spherically symmetric radial outflow (with velocity u(r)). Then by conservation of momentum,

\rho u \frac{d u}{d r} = -\frac{d p}{d r} -\rho \frac{G M}{r^2}$

Continuity (conservation of particle #) gives

\frac{1}{r^2} \frac{d (r^2 \rho u)}{d r} = 0$

If we assume an isothermal corona, the ideal gas law is

p = \frac{2 \rho T}{m}$

Continuity shows that r^2 \rho u = C, where C is a constant. If we substitute for p and rho into the momentum equation, we get a differential equation in u and r:

\frac{1}{u} \frac{d u}{d r}(u^2 - \frac{2 T}{m_p}) = \frac{4 T}{m_p r} - \frac{G M}{r^2}$

Define r_c = \frac{G M m_p}{4 T} so that u_c^2 = \frac{2 T}{m_p} (i.e., sonic radius); substitute and rearrange to find

\frac{d u^2}{u_c^2} - \frac{d u^2}{u^2} = \frac{4 d r}{r} - \frac{4 r_c d r}{r^2}$

Integrate both sides to find the Parker wind equation,

\psi - \log(\psi) = 4 \log(\xi) + \frac{4}{\xi} + C$, where $\psi = (\frac{u}{u_c})^2$ and $\xi = \frac{r}{r_c}$

Plotting

Solving the Parker wind equation:

\psi - \log(\psi) = 4 \log(\xi) + \frac{4}{\xi} - 3$

There are 4 types of solutions - combinations of sub/supersonic at the surface and v → 0/v_f as xi → infinity. The only physical solution is subsonic at the surface and has a nonzero velocity infinitely far away, and passes through xi = 1, psi = 1.

First attempted to solve equation in Mathematica - solution appears correct up to the sonic radius, but the incorrect solution (v → 0) appears for r > r_s. Next attempted to plot approximations for r << r_s and r >> r_s in Matlab, but as would be expected, the approximations fail near r = r_s. Numerically solving the equation in Matlab gives the opposite case of Mathematica - correct solution for r > r_s, but not for r < r_s.

Looking at Jonathan's code for calculating the Parker wind to find the correct approximation, try to implement it in a different way to test understanding. Works for r > r_s, but not r < r_s.

Pertinent line is: yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i),y,(xx(i) > 1) * xx(i) + (xx(i) ⇐ 1)*exp(3-4*log(xx(i))-4/xx(i)))

Interpreted as:

if xx(i) ≥ 1

yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i));

else

yy(i)=vpasolve(y - log(y) == exp(3-4*log(xx(i))-4/xx(i)));

end

Correct & incorrect plots:

Numerical Integration of DE

Write the differential equation for \frac{d u}{d r} as a parameterized system:

\frac{d u}{d s} = -2uc^2(r-\frac{GM}{2c^2})$ $\frac{d r}{d s} = -r^2(u^2-c^2)$

Calculate the Jacobian and evaluate at the critical point to linearize the system, then move a small distance along the stable eigenvector (i.e. tangent to the stable manifold, the transonic solution of interest). Integrate from these points out.

Jacobian = \begin{bmatrix} 0 & 2c^3 \\ \frac{GM}{2c^3} & 0 \end{bmatrix}