Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

A solar model

  1. Feb 26, 2010 #1

    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.
    This density function contains three dimensionless parameters [itex]a, \delta, \gamma[/itex] 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 [itex]a[/itex] be written as [itex]\alpha[/itex] in the formal equation?

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

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

    Ideal gas law:
    [tex]P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_u}[/tex]

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

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

    Is this equation correct?

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

    Why is this pressure_equation missing a volume dimension?

    Solar core range:
    [tex]r_1 = (0 \to .25) R_{\odot} = .25 R_{\odot}[/tex]

    Radiative zone range:
    [tex]dr_2 = (.25 \to 0.71) R_{\odot} = 0.46 R_{\odot}[/tex]

    Convective zone range:
    [tex]dr_3 = (0.71 \to 1) R_{\odot} = 0.29 R_{\odot}[/tex]

    http://neutrino.aquaphoenix.com/ReactionDiffusion/SERC4chap4.pdf" [Broken]
    http://en.wikipedia.org/wiki/Sun" [Broken]
    http://www.nasa.gov/worldbook/sun_worldbook.html" [Broken]
    http://astro.berkeley.edu/~eliot/Astro252/hw2.pdf" [Broken]
    http://en.wikipedia.org/wiki/Polytrope" [Broken]
    Last edited by a moderator: May 4, 2017
  2. jcsd
  3. Mar 6, 2010 #2
    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:
    [tex]\frac{dT}{dr} = - \frac{3 \kappa(r) \rho(r) L(r)}{64 \pi \sigma r^2 T(r)^3}[/tex]

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

    Which equation is correct?

    Is [itex]\overline{\chi}[/itex] defined as the 'average opacity'?

    Is [tex]L_c[/tex] 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?

    [tex]\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[/tex]

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

    Any collaboration in building this model would be appreciated.

    http://astro.elte.hu/~kris/napfiz/ssm.pdf" [Broken]
    Last edited by a moderator: May 4, 2017
  4. Mar 7, 2010 #3
    Dear Orion,
    Here's a good method for solving Emden's equation here

    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.
    Last edited by a moderator: Apr 24, 2017
  5. Mar 7, 2010 #4
    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 [itex]\theta(\xi)[/itex].

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

    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 [itex]\rho(\xi)[/itex]:
    [tex]\theta(\xi) = \frac{\rho(\xi)}{\rho_c}[/tex]

    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:
    [tex]P(r) = \frac{k_B \rho(r) T(r)}{\bar{\mu} m_H} [/tex]

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

    The difference is:

    Hydrogen mass:
    [tex]m_H = 1.67372 \cdot 10^{-27} \; \text{kg}[/tex]

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

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

    Which equation definition is correct?

    http://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation" [Broken]
    http://mathworld.wolfram.com/Lane-EmdenDifferentialEquation.html" [Broken]
    http://en.wikipedia.org/wiki/Atomic_mass_constant" [Broken]
    http://en.wikipedia.org/wiki/Atomic_mass_unit" [Broken]
    http://www.vikdhillon.staff.shef.ac.uk/teaching/phy213/phy213_polytropes.html" [Broken]
    Last edited by a moderator: May 4, 2017
  6. Mar 10, 2010 #5
    Numerically solving the Lane-Emden differential equation...

    Numerically solving the Lane-Emden differential equation...

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

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

    Solutions for all other values must be solved numerically.

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

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

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

    [tex]\theta(\xi)_{i + 1} = \theta(\xi)_i + \Delta \xi \left( \frac{d \theta}{d \xi} \right)_i[/tex]

    Now [itex]\frac{d \theta}{d \xi}[/itex] is unknown, however we can write:
    [tex]\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[/tex]

    Then we can replace the second derivative term in the above by the rearranged form of the Lane-Emden differential equation:
    [tex]\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 [/tex]

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

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

    Which can also be used to determine:
    [tex]\theta(\xi)_{i + 1}[/tex]

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

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

    http://cc.oulu.fi/~jpoutane/teaching/STELL09/l07.pdf" [Broken]
    Last edited by a moderator: May 4, 2017
  7. Mar 12, 2010 #6
    Lane-Emden series expansion...

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

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

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

    [tex] \frac{d \theta}{d \xi} = -\frac{\xi}{3} + \frac{n\xi^3}{30} \; \text{. . .}[/tex]

    Can anyone here demonstrate how this series expansion was generated?

    Last edited by a moderator: Apr 24, 2017
  8. Mar 13, 2010 #7
    expansion series equation...

    [tex]\theta_n(\xi) = 1 - \frac{\xi^2}{6} + \frac{n \xi^4}{120} - \left( \frac{8n^2 - 5n}{15120} \right) \xi^6 + \; \text{. . .}[/tex]
    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.

    [tex]\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{. . .}[/tex]

    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.

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

    Attached Files:

    Last edited by a moderator: May 4, 2017
  9. Jun 9, 2010 #8
    Numerically solving the Lane-Emden differential equation...

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

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

    http://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation" [Broken]
    Last edited by a moderator: May 4, 2017
  10. Aug 2, 2010 #9

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

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

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

    Is this temperature gradient equation separable?
    [tex]\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[/tex]

    Science Advisor requested.
    Last edited: Aug 2, 2010
  11. Aug 8, 2010 #10
    hey guys..i need any input i can get on "Solar Tsunami"..ty..:)
  12. Mar 27, 2011 #11
    Science Advisor requested.
  13. Apr 3, 2011 #12

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

    Mass continuity equation:
    [tex]\frac{dm}{dr} = 4 \pi r^2 \rho(r)[/tex]

    Mass continuity function:
    [tex]m(r) = 4 \pi \int_0^{r} r^2 \rho(r) dr[/tex]

    Integration via substitution:
    [tex]\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[/tex]

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

    Is this equation correct?

    Last edited by a moderator: Apr 25, 2017
  14. Apr 4, 2011 #13


    User Avatar
    Science Advisor
    Gold Member

    The wiki EOS for hydrostatic equilibrium, IMO, does not properly account for relativistic effects.
  15. Apr 4, 2011 #14

    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:
    [tex]\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}[/tex]

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

    http://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation" [Broken]
    Last edited by a moderator: May 5, 2017
  16. Apr 19, 2011 #15
    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?
  17. Apr 19, 2011 #16
    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://astro.ensc-rennes.fr/index.php?pw=ycesam [Broken]
    Last edited by a moderator: May 5, 2017
  18. Apr 19, 2011 #17
    Re: Lane-Emden series expansion...


    Last edited by a moderator: Apr 25, 2017
  19. Apr 19, 2011 #18
    Re: Lane-Emden differential equation...

    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.
  20. Apr 19, 2011 #19
    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.


    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.
  21. Apr 19, 2011 #20
    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook