Mathematical Solar Model: Core, Radiative Zone, and Convective Zone

[Orion1]

I am interested in building a three shell mathematical solar model based upon the core, the radiative zone, convective zone.

I have been researching to determine a way to build my first mathematical solar model. The best mathematical model for the core density function that I have located is from reference 1, model 4.

Note that the models in (4.1.13) and (4.1.1) are really valid only in the interior core of the Sun. The convective zone has an entirely different behavior for the density distribution. But for computing the total mass, pressure, temperature and luminosity, we have integrated out over the entire length of the solar radius $R_{\odot}$. This is not appropriate. The integration should have been done only in the interior core of the Sun. Hence a more appropriate model is of the following form:
$$\rho(r) = \rho_0 \left( 1 - a \left( \frac{r}{R_{\odot}} \right)^{\delta} \right)^{\gamma} \; \; \; \delta > 0 \;, \gamma > 0 \;, a > 0 \; \; \; \text{(Model 4)}$$

Since $1 - ax^{\delta} > 0 \;, x = \frac{r}{R_{\odot}}$, we have $0 \leq x \leq \frac{1}{a^{\frac{1}{\delta}}}$. Hence for the total integral, the range should have been $0 \leq r \leq \frac{R_{\odot}}{a^{\frac{1}{\delta}}}$.

This density function contains three dimensionless parameters $a, \delta, \gamma$ and the reference does not describe any other equation solutions for model 4.

However, the impression I receive from reference 1, is that this density function is valid only in the interior core of the Sun.

I am uncertain that the other models that I have researched have addressed the issue of the different behavior for the density distribution functions.

What are the correct density functions for the core and the radiative zone and the convective zone of the Sun?

And should the dimensionless parameter $a$ be written as $\alpha$ in the formal equation?

According to reference 5, radiative transport of energy is described by the radiative temperature_gradient equation:
$$\frac{dT}{dr} = - \frac{3 \kappa \rho(r) L(r)}{64 \pi r^2 \sigma T(r)^3}$$

According to reference 4, I should model the convection zone as a polytrope:
$$P(r) = K \rho(r)^{\gamma}$$

Ideal gas law:
$$P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_u}$$

Polytropic convective zone pressure:
$$P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_u} = K \rho(r)^{\gamma}$$

Polytropic ideal gas convective zone density:
$$\boxed{\rho(r) = \left(\frac{K \bar{\mu} m_u}{k_B T(r)} \right)^{\frac{1}{1 - \gamma }}}$$

Is this equation correct?

Reference 7 defines the polytropic_pressure as:
$$P(r) = \frac{N_A k_B T(r)}{\bar{\mu}}$$

Why is this pressure_equation missing a volume dimension?

Solar core range:
$$r_1 = (0 \to .25) R_{\odot} = .25 R_{\odot}$$

Radiative zone range:
$$dr_2 = (.25 \to 0.71) R_{\odot} = 0.46 R_{\odot}$$

Convective zone range:
$$dr_3 = (0.71 \to 1) R_{\odot} = 0.29 R_{\odot}$$

Reference:
http://neutrino.aquaphoenix.com/ReactionDiffusion/SERC4chap4.pdf"
http://en.wikipedia.org/wiki/Sun"
http://www.nasa.gov/worldbook/sun_worldbook.html"
http://astro.berkeley.edu/~eliot/Astro252/hw2.pdf"
http://en.wikipedia.org/wiki/Standard_Solar_Model#Equations_of_state"
http://en.wikipedia.org/wiki/Ideal_gas_law#Alternative_Forms"
http://en.wikipedia.org/wiki/Polytrope"

 

[Orion1]

SSM temperature function...


I am attempting to derive the temperature function for the Solar Standard Model (SSM), however the published internet papers that I have already researched simply state that it is an 'iterative process' and then provide a list of differential equations.

According to my research, this particular temperature gradient equation is valid in the solar core and radiative zone only.

Reference 1 defines the temperature gradient equation as:
$$\frac{dT}{dr} = - \frac{3 \kappa(r) \rho(r) L(r)}{64 \pi \sigma r^2 T(r)^3}$$

Reference 2 defines the temperature gradient equation as:
$$\frac{dT}{dr} = - \frac{3 \overline{\chi} \rho(r) (L(r) - L_c)}{16 \pi a c r^2 T(r)^3}$$

Which equation is correct?

Is $\overline{\chi}$ defined as the 'average opacity'?

Is $$L_c$$ defined as the core luminosity from mixing length theory?

What is the key definition for [itex]a c[/tex] in the second equation definition?

Is this temperature gradient equation separable?

$$\int_0^{r_1 + dr_2} T(r)^3 dT = - \frac{3}{64 \pi \sigma} \int_0^{r_1 + dr_2} \frac{\kappa(r) \rho(r)(L(r) - L_c)}{r^2} dr$$

Is this equation and integration correct with respect to to the Solar Standard Model (SSM)?

Any collaboration in building this model would be appreciated.

Reference:
http://en.wikipedia.org/wiki/Standard_Solar_Model#Energy_transport_in_the_Sun"
http://astro.elte.hu/~kris/napfiz/ssm.pdf"

 

[Helios]

Dear Orion,
Here's a good method for solving Emden's equation here
http://www.vikdhillon.staff.shef.ac.uk/teaching/phy213/phy213_le.html

Beyond this you could then go for a composite polytrope with a radiative core ( polytropic index = 3 ) and a convective envelope ( polytropic index = 3/2 ). The core is taken to be an "E-function" ( a sol'n to Emden's eq. with normal boundry conditions ). The envelope is taken to be an "M-function" ( a sol'n to Emden's eq. with special boundry conditions ). The tricky part is to solve the EQUATIONS OF FIT, so the E-function merges with the M-function. I haven't totally figured this all out myself though.

The inner luminosity core ( place of fusion ) is a whole other story. I've heard of a triple composite polytrope where the core of luminosity was given a very high polytropic index. I don't know why though. I've never found an estimate of its radius with any explanation.

 

[Orion1]

Lane-Emden differential equation...


If I understand the calculus correctly, the first step in modeling an index n = 3 polytropic solar Eddington Standard Model (ESM) requires solving the Lane-Emden differential equation for an index n = 3 and $\theta(\xi)$.

Lane-Emden differential equation:
$$\frac{1}{\xi^2} \frac{d}{d\xi} \left({\xi^2 \frac{d\theta}{d\xi}}\right) + \theta(\xi)^n = 0$$

However, analytical solutions are known for only for index n = 0,1, and 5.

[tex]\theta(\xi) = 1 - \frac {\xi^2}{6} \; \; \; n = 0[/itex]
[tex]\theta(\xi) = \frac{\sin\xi}{\xi} \; \; \; n = 1[/itex]
[tex]\theta(\xi) = \left(1 + \frac{\xi^2}{3}\right)^{-\frac{1}{2}} \; \; \; n = 5[/itex]

Unfortunately, the reference 6 definition does not define the density function $\rho(\xi)$:
$$\theta(\xi) = \frac{\rho(\xi)}{\rho_c}$$

How was this graph numerically plotted with all the indexes from 0 to 5?

I also notices another difference in formulation regarding the alternate density form of the Ideal gas law, all the models I have studied so far define the Ideal gas law as:
$$P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_H}$$

However, reference 3 defines the Ideal gas law as:
$$P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_u}$$

The difference is:

Hydrogen mass:
$$m_H = 1.67372 \cdot 10^{-27} \; \text{kg}$$

Atomic mass constant:
$$m_u = 1.660 538 782(83) \cdot 10^{-27} \; \text{kg}$$

If using the atomic mass constant definition is correct, then the Hydrogen mass definition is incorrect.

Which equation definition is correct?

Reference:
http://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation"
http://mathworld.wolfram.com/Lane-EmdenDifferentialEquation.html"
http://en.wikipedia.org/wiki/Ideal_gas_law#Alternative_Forms"
http://en.wikipedia.org/wiki/Atomic_mass_constant"
http://en.wikipedia.org/wiki/Atomic_mass_unit"
http://www.vikdhillon.staff.shef.ac.uk/teaching/phy213/phy213_polytropes.html"

 

[Orion1]

Numerically solving the Lane-Emden differential equation...


Numerically solving the Lane-Emden differential equation...

Lane-Emden differential equation:
$$\frac{1}{\xi^2} \frac{d}{d\xi} \left({\xi^2 \frac{d\theta}{d\xi}}\right) + \theta(\xi)^n = 0$$

However, analytical solutions are known for only for polytropic index n = 0,1, and 5.
$$\theta(\xi) = 1 - \frac {\xi^2}{6} \; \; \; n = 0$$
$$\theta(\xi) = \frac{\sin\xi}{\xi} \; \; \; n = 1$$
$$\theta(\xi) = \left(1 + \frac{\xi^2}{3}\right)^{-\frac{1}{2}} \; \; \; n = 5$$

Solutions for all other values must be solved numerically.

All solutions are subject to boundary conditions:
$$\frac{d \theta}{d \xi} = 0 \; \; \; \theta(\xi) = 1 \; \text{at} \; \xi = 0$$

Express the Lane-Emden differential equation in the form:
$$\frac{d^2 \theta}{d \xi^2} = - \frac{2}{\xi} \frac{d \theta}{d \xi} - \theta(\xi)^n$$

The numerical integration technique involves stepping outwards in radius from the center of the star and evaluate density at each radius $\theta(\xi)$. At each radius, $\theta(\xi)_i$ plus the change in density over the step $\Delta \xi$:

$$\theta(\xi)_{i + 1} = \theta(\xi)_i + \Delta \xi \left( \frac{d \theta}{d \xi} \right)_i$$

Now $\frac{d \theta}{d \xi}$ is unknown, however we can write:
$$\left( \frac{d \theta}{d \xi} \right)_{i + 1} = \left( \frac{d \theta}{d \xi} \right)_i + \Delta \xi \left( \frac{d^2 \theta}{d \xi^2} \right)_i$$

Then we can replace the second derivative term in the above by the rearranged form of the Lane-Emden differential equation:
$$\left( \frac{d \theta}{d \xi} \right)_{i + 1} = \left( \frac{d \theta}{d \xi} \right)_i - \left( \frac{2}{\xi} \frac{d \theta}{d \xi} + \theta(\xi)^n \right)_i \Delta \xi$$

Now we can adopt a value for $n$ and integrate numerically. We have the boundary conditions at the center:
$$\frac{d \theta}{d \xi} = 0 \; \; \; \theta(\xi) = 1 \; \text{at} \; \xi = 0$$

Starting at the center, we determine:
$$\left( \frac{d \theta}{d \xi} \right)_{i + 1}$$

Which can also be used to determine:
$$\theta(\xi)_{i + 1}$$

The radius is then incremented by adding $\Delta \xi$ to $\xi$ and the process is repeated until the surface of the star is reached at which point $\theta(\xi)$ becomes negative.

Can anyone here numerically demonstrate the first two steps of this particular numerical integration process?

Reference:
http://cc.oulu.fi/~jpoutane/teaching/STELL09/l07.pdf"

 

[Orion1]

Lane-Emden series expansion...


Express the Lane-Emden differential equation in the form:
$$\frac{d^2 \theta}{d \xi^2} = - \frac{2}{\xi} \frac{d \theta}{d \xi} - \theta(\xi)^n$$

We have the boundary conditions at the center:
$$\frac{d \theta}{d \xi} = 0 \; \; \; \theta(\xi) = 1 \; \text{at} \; \xi = 0$$

We calculate starting values using a series expansion for small [itex]\xi[/tex]:
$$\theta(\xi) = 1 - \frac{\xi^2}{6} + \frac{n \xi^4}{120} \; \text{. . .}$$

$$\frac{d \theta}{d \xi} = -\frac{\xi}{3} + \frac{n\xi^3}{30} \; \text{. . .}$$

Can anyone here demonstrate how this series expansion was generated?

Reference:
http://www.astro.uu.nl/~pols/education/stev/practicum/practicum2.pdf"

 

[Orion1]

expansion series equation...

For other values of the polytropic index $n$, it is possible to develop a series solution which is useful for starting many numerical methods for the solution. The first few terms in the solution are:

$$\theta_n(\xi) = 1 - \frac{\xi^2}{6} + \frac{n \xi^4}{120} - \left( \frac{8n^2 - 5n}{15120} \right) \xi^6 + \; \text{. . .}$$
However, none of the reference provide an explanation as to how this equation solution was derived.

Unfortunately, when I plot this solution for all polytropic indexes from 0 to 5, the graph does not match the correct functions, listed as first attachment.

Power series solutions about $\xi = 0$ are given, to facilitate such computations, according to Sugimoto:

$$\theta_n(\xi) = 1 - \frac{\xi^2}{6} + \frac{n \xi^4}{120} - \frac{n(8n - 5)}{42 \cdot 360} \xi^6 + \frac{n(122n^2 - 183n + 70)}{42 \cdot 360 \cdot 216} \xi^8 \; \text{. . .}$$

Unfortunately, when I plot this solution for all polytropic indexes from 0 to 5, the graph does not match the correct functions, listed as second attachment.

The correct functions are listed in third attachment.

Reference:
http://ads.harvard.edu/books/1989fsa..book/AbookC02.pdf"
http://www.mpa-garching.mpg.de/~weiss/Cox_Vol_II_CD/ch23.pdf"

 

[Orion1]

Numerically solving the Lane-Emden differential equation...



Lane-Emden differential equation:
$$\frac{1}{\xi^2} \frac{d}{d\xi} \left({\xi^2 \frac{d\theta}{d\xi}}\right) + \theta(\xi)^n = 0$$

How do I numerically integrate the Lane-Emden differential equation with polytropic index n = 3?

Reference:
http://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation"

 

[Orion1]

I am attempting to complete the Equation of State for hydrostatic equilibrium for the sun:
$$\frac{d P(r)}{dr} = G \frac{M(r) \rho(r)}{r^2}$$

And completing the temperature gradient equation:
$$\frac{dT}{dr} = - \frac{3 \kappa(r) \rho(r) L(r)}{64 \pi \sigma r^2 T(r)^3}$$

I am specifically searching for an equation solution for $$T(r)$$ that obeys the Solar Standard Model (SSM).

Is this temperature gradient equation separable?
$$\int_0^{r_1 + dr_2} T(r)^3 dT = - \frac{3}{64 \pi \sigma} \int_0^{r_1 + dr_2} \frac{\kappa(r) \rho(r) L(r)}{r^2} dr$$

Science Advisor requested.

 

[sameerpaisari]

hey guys..i need any input i can get on "Solar Tsunami"..ty..:)

 

[Orion1]

Science Advisor requested.

 

[Orion1]

Equation of State for hydrostatic equilibrium for a star:
$$\frac{dP(r)}{dr} = - G \frac{m(r) \rho(r)}{r^2}$$

Mass continuity equation:
$$\frac{dm}{dr} = 4 \pi r^2 \rho(r)$$

Mass continuity function:
$$m(r) = 4 \pi \int_0^{r} r^2 \rho(r) dr$$

Integration via substitution:
$$\frac{dP(r)}{dr} = -G \frac{m(r) \rho(r)}{r^2} = - \frac{4 \pi G \rho(r)}{r^2} \int_0^{r} r^2 \rho(r) dr$$

Equation of State for hydrostatic equilibrium:
$$\boxed{\frac{dP(r)}{dr} = - \frac{4 \pi G \rho(r)}{r^2} \int_0^{r} r^2 \rho(r) dr}$$

Is this equation correct?

Reference:
http://en.wikipedia.org/wiki/Stellar_structure#Equations_of_stellar_structure"

 

[Chronos]

The wiki EOS for hydrostatic equilibrium, IMO, does not properly account for relativistic effects.

 

[Orion1]

The wiki EOS for hydrostatic equilibrium, IMO, does not properly account for relativistic effects.

According to my understanding, the mainstream Standard Solar Model (SSM) uses the Newtonian EOS for hydrostatic equilibrium that is listed on Wikipedia, which would fail for the highest stellar masses, where General Relativity effects become mathematically significant.

Tolman–Oppenheimer–Volkoff equation:
$$\frac{dP(r)}{dr}=-\frac{G}{r^2}\left[\rho(r)+\frac{P(r)}{c^2}\right]\left[M(r)+4\pi r^3 \frac{P(r)}{c^2}\right]\left[1-\frac{2GM(r)}{c^2r}\right]^{-1}$$

The Tolman–Oppenheimer–Volkoff equation under General Relativity should be used for the Equation of State for hydrostatic equilibrium for stellar modeling?

Reference:
http://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation"

 

[cvvvarma]

what will be the values of the parameters,ie, a ,delta ,gamma which are mentioned in the density function at the very beginning of this forum?

 

[twofish-quant]

The Tolman–Oppenheimer–Volkoff equation under General Relativity should be used for the Equation of State for hydrostatic equilibrium for stellar modeling?

Not for any "reasonable" stellar model. Even high mass stellar routines usually don't use GR. It turns out that GR doesn't make a difference.

Also what you are trying to do by including a convective region is going to be very difficult to do non-numerically. The problem is that the temperature gradients determine where the convective regions are which then impact the temperature gradients. What you do is to dump everything into a computer and then let it figure out a consistent solution.

You might start with a preexisting code

http://mesa.sourceforge.net/
http://astro.ensc-rennes.fr/index.php?pw=ycesam

 

[twofish-quant]

Express the Lane-Emden differential equation in the form:
$$\frac{d^2 \theta}{d \xi^2} = - \frac{2}{\xi} \frac{d \theta}{d \xi} - \theta(\xi)^n$$

We have the boundary conditions at the center:
$$\frac{d \theta}{d \xi} = 0 \; \; \; \theta(\xi) = 1 \; \text{at} \; \xi = 0$$

We calculate starting values using a series expansion for small [itex]\xi[/tex]:
$$\theta(\xi) = 1 - \frac{\xi^2}{6} + \frac{n \xi^4}{120} \; \text{. . .}$$

$$\frac{d \theta}{d \xi} = -\frac{\xi}{3} + \frac{n\xi^3}{30} \; \text{. . .}$$

Can anyone here demonstrate how this series expansion was generated?

Reference:
http://www.astro.uu.nl/~pols/education/stev/practicum/practicum2.pdf"

Try...

http://articles.adsabs.harvard.edu//full/1980Ap&SS..73..227M/0000230.000.html

 

[twofish-quant]

Which equation definition is correct?

Neither or both. The difference is much smaller than the other errors that you have in the model so it doesn't matter which one you choose.

If you are doing high precision work it may matter, but in that situation the person that hands you the values of mu that you are working with will tell you what units those things are in. Astronomers will usually quote mu in terms of the mass of hydrogen, whereas I think people in other fields will quote in terms of atomic mass units.

 

[twofish-quant]

[The Tolman–Oppenheimer–Volkoff equation under General Relativity should be used for the Equation of State for hydrostatic equilibrium for stellar modeling?

A bad idea unless you have something in your doctoral dissertation that requires it.

The problem is that if you treat hydrostatic equilibrium with GR, you'll have to treat everything with GR. If your pressure equations are GR, and your radiation equations are Newtonian you end up with a big, big mess that's worse than if everything was Newtonian.

Intuitively you can think about this by imagining a black hole. When a black hole forms, the stuff inside of the event horizon stops radiating, and are treating pressure with GR and radiation Newtonian, you'll get a big mess.

One thing that you should realize that the TOV doesn't work if you have a radiating body. If the thing you are trying to study is radiating, then the T01 component of the stress-energy tensor is non-zero, and you have to rederive everything.

There is a "standard" way of bolting on GR to stellar codes.

http://adsabs.harvard.edu/abs/1979ApJ...232..558V

What you do is to run an Newtonian code and then set up a "GR time correction" factor and a "GR volume correction" factor that corrects for everything. If you do this for any sort of stellar code, you'll find that both correction factors are so close to one that it doesn't matter.

 

[twofish-quant]

Also there are two "philosophies" of stellar modelling. There is the "cartoon" approach which were leave out everything accept for the important parts so that you have something that you can understand. The thing about doing things via cartoon is that it helps a lot in understanding how things work, but you aren't going to get anything that is close to something matches the sun.

You can also try a photo realistic approach in which you dump everything into a computer. You get realistic output, but there is a loss of insight.

Something that would be interesting would be try to create the models in python or mathematica.

 

[Orion1]

I note the advice from the Science Advisors and provide my current 'cartoon' model 'college trial' for core pressure, for the benefit of insight...

Experimentally determined parameters:
Total solar radius:
$$R_{\odot} = 6.955 \cdot 10^8 \; \text{m}$$
Total solar core density:
$$\rho_c = 1.622 \cdot 10^5 \; \frac{\text{kg}}{\text{m}^3}$$

Stellar density function:
$$\rho(r) = \rho_c \left( 1 - \alpha \left( \frac{r}{R_{\odot}} \right)^{\delta} \right)^{\gamma} \; \; \; \delta > 0 \;, \gamma > 0 \;, \alpha > 0$$

Integration via substitution (ref. 2):
$$\frac{dP(r)}{dr} = - \frac{4 \pi G \rho(r)}{r^2} \int_0^{r} r^2 \rho(r) dr$$

Equation of State for hydrostatic equilibrium:
$$\frac{dP(r)}{dr} = - \frac{4 \pi G \rho_c^2 r}{3} \left( 1 - \alpha \left( \frac{r}{R_{\odot}} \right)^{\delta} \right)^{\gamma} \cdot _2F_1 \left( - \gamma, \frac{3}{\delta}, \frac{3 + \delta}{\delta}, \alpha \left( \frac{r}{R_{\odot}} \right)^{\delta} \right) \; \; \; \; \; \; \alpha = 1 \; \; \; \delta = 2 \; \; \; \gamma = 1$$

Total stellar core pressure integration:
$$P_c = \int_0^{R_{\odot}} \left( \frac{dP(r)}{dr} \right) dr$$

Orion1 solar model core pressure:
$$\boxed{P_c = -7.115 \cdot 10^{17} \; \frac{\text{N}}{\text{m}^2}}$$

Standard Solar Model core pressure (SSM):
$$P_c = -2.477 \cdot 10^{16} \; \frac{\text{N}}{\text{m}^2}$$

Known problems with 'cartoon' modeling approach:
1. The model is static and does not evolve.
2. The model only uses experimentally determined parameters.
3. There is a functional discontinuity in the SSM, ZAMS, and Lane-Emden (Eddington) models on the radiative-convective shell boundary.
4. Newtonian laws are incompatible with General Relativity.


Reference:
http://en.wikipedia.org/wiki/Sun"
https://www.physicsforums.com/showpost.php?p=3226531&postcount=12"
http://en.wikipedia.org/wiki/Stellar_structure#Equations_of_stellar_structure"
http://home.earthlink.net/~dckennedy/pubs/Sun991110ApJ.pdf"
http://reference.wolfram.com/mathematica/ref/Hypergeometric2F1.html"
http://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html"

 

[Orion1]

Would adding a positive differential radiation pressure that opposes negative differential gravitational pressure to the Newtonian Equation of State for hydrostatic equilibrium produce a more accurate stellar core pressure reading?

$$\frac{dP(r)}{dr} = \frac{dP_g(r)}{dr} + \frac{dP_r(r)}{dr}$$

Stellar Equation of State for hydrostatic equilibrium:
$$\frac{dP_g(r)}{dr} = - G \frac{m(r) \rho(r)}{r^2}$$

What would the equation be for stellar interior differential radiation pressure?

Reference:
http://en.wikipedia.org/wiki/Radiation_pressure#In_stellar_interiors"

 

[Orion1]

You can also try a photo realistic approach in which you dump everything into a computer. You get realistic output, but there is a loss of insight.

Something that would be interesting would be try to create the models in python or mathematica.

My Equation of State for hydrostatic equilibrium equation seems sophisticated enough to use a photo realistic approach to modulate the constants α, δ, γ using a realistic output computer data dump into Mathematica of the mass and density functions for the Standard Solar Model (SSM) at present time.

Would the constants α, δ, γ become stellar evolution functions of time using such an approach?

Is Mesa star capable of doing this?

Reference:
Stellar structure - Equations of stellar structure - Wikipedia
Mesa star - stellar evolution code

 

[Vanadium 50]

This thread has turned into one person talking to himself.