Stellar Nursery — how it works

1 · What this is, and what it is not

It is a zero-age, one-dimensional, hydrostatic snapshot. You give it a total mass and a uniform composition. It integrates the structure equations outward from the centre, finds the central density and temperature that make the model hold together with the mass you asked for and radiate exactly what its nuclear burning supplies, and then names what came back.

It is not any of the following:

• Not a stellar evolution code. There is no time coordinate. Nothing ages, contracts, expands, mixes, or loses mass. MESA solves a sequence of models linked by time; this solves exactly one model, and it is always the zero-age one.
• Not a menu of star types. There is no table of outcomes keyed on mass. The words at the top of the results panel are chosen after the integration, from the numbers the integration produced. Set Z = 0.99 and you do not get a hard-coded "black hole" row — you get whatever a ball of metal at that mass turns out to be.
• Not an atmosphere model. The integration stops where the pressure runs out. There is no photosphere, no opacity table, no model atmosphere; T_eff is inferred from L and R, not integrated to.
• Not a relativistic code. It solves Newtonian hydrostatic equilibrium. It can tell you that a configuration is unstable to collapse, and where the horizon would be, but it cannot follow anything through collapse or model what comes out the far side.
• Not accurate. Factor-of-two thinking throughout. The Sun lands within a few percent because one opacity constant was set to make it, and everything else follows from the physics with no further tuning — which means everything else is as right or as wrong as the physics is.

2 · Inputs and normalisation

Total mass M, on a log slider from 10^-3 to 10⁸ M_☉. Type a value into the box instead if you want a particular one; the unit beside it is yours to pick, solar masses or Jupiter masses, and the line underneath always gives the same mass written the other way. The upper end runs far past any real star because that is where the collapse checks in §5 start firing — everything above a few times 10⁴ M_☉ is there to be watched failing, not to be believed as an object.

Mass fractions X, Y, Z — hydrogen, helium, everything else. Slide them or type them; the three always sum to exactly 1, and setting one rescales the other two in proportion, so you never have to make room first. Composition is uniform: one mix, everywhere, all the way through. No onion shells, no burning-modified core, no time slider, because there is no time.

Composition feeds the physics through three mean molecular weights:

The first is the form the brief asked for. The ion weight normally follows as $1/\mu -1/{\mu }_e=X+Y/4$, which quietly assumes metal ions contribute nothing; that is fine at Z = 0.02 and absurd at Z = 0.99, where it would leave the gas with no ion pressure at all, so a Z/16 term is added for metals at <A> ≈ 16.

The presets are just composition shortcuts: solar (0.70 / 0.28 / 0.02), primordial (0.75 / 0.25 / 0), metal-poor, metal-rich, helium star (Y = 0.98), C/O (all metals), and almost-pure metals (Z = 0.99). None of them changes the solver's behaviour — they only move the sliders.

3 · The equations it actually uses

Structure

Mass continuity and hydrostatic equilibrium, integrated outward in radius:

Rather than carry P and invert the equation of state at every step, the code carries ρ directly and uses the chain rule, which removes the root-find from the inner loop and is why the whole solve runs in tens of milliseconds:

Equation of state

Ideal gas, radiation, and electron degeneracy:

The degenerate branch is the full Chandrasekhar expression, so it runs continuously from non-relativistic through relativistic and gives sensible white-dwarf radii without a special case. The electron pressure is the smooth maximum (a⁴ + b⁴)^1/4 of the ideal and degenerate branches — a cartoon that overshoots by at most 2^1/4 ≈ 19% right at the crossover, in exchange for a P(ρ) that stays smooth and monotonic, which is what the integrator needs.

Thermal stratification

The transport equation is not solved. Instead the temperature is tied to the density by an effective polytrope, which is Eddington's standard model:

This is the single biggest simplification in the program, and section 8 says what it costs. The adiabatic gradient reported in the diagnostics is Chandrasekhar's for a gas + radiation mix, ${\nabla }_{\text{ad}}=\frac{2(4-3\beta )}{32-24\beta -3{\beta }^2}$, which runs from 0.4 for pure gas to 0.25 for pure radiation.

Opacity

Electron scattering plus Kramers, with the textbook coefficients converted to SI:

The bound-free guillotine factor g/t is uncertain at the factor-of-ten level in the literature. It is fixed here at 0.033, and that is the one number in this program tuned to an answer rather than derived: it is the value that puts a 1 M_☉ solar-composition model on the Sun. Everything else — the mass scaling, the composition scaling, the whole main sequence — comes out of the equations after that one constant is set.

Nuclear burning

No network. Four one-line rates, each a power law times a Gamow exponential, evaluated on the local ρ, T and mix. They exist to answer two questions — is the centre hot enough to burn, and on which channel — and are used for nothing else.

The reported channel is whichever of the four is largest at the centre. The carbon rate is the crudest thing in the file: the exponential has roughly the right shape and the coefficient was chosen to make ignition land near T ≈ 6–8 × 10⁸ K. Read it as a flag, not as a rate.

4 · The shooting method, and why failing is an answer

Two central unknowns, ρ_c and T_c, so two conditions are needed. The integration starts just off the singular centre — at the radius enclosing 10^-9 of the target mass, where the series corrections are still O(r²) and irrelevant — and runs outward with RK4 and an adaptive step of a few percent of the local scale length. The surface is where the density runs out; the code stops at ρ = 10^-8 ρ_c and calls that R.

Condition 1 — mass

m(R) = M. At fixed T_c, the mass runs down as ρ_c rises while thermal support dominates and back up once degeneracy takes over, so M(ρ_c) has a minimum and there can be two roots. Which side you want is a physical question, not a numerical one: the gas-supported root is the star, the degenerate root is the brown dwarf or the white dwarf. The solver states which branch it is on and refuses gas-branch answers that come back strongly degenerate.

Condition 2 — thermal equilibrium

Nuclear burning has to supply exactly what the interior radiates away. Integrating the radiative diffusion equation from the centre out to a deep reference point gives the luminosity the structure can actually carry:

The solver finds the T_c where L_nuc = L_diff. This is what pins a star onto the zero-age main sequence, and it is well behaved to solve because the two sides depend on temperature so differently: burning rates go as T⁴ to T¹⁸, transport roughly as T⁴. The residual is steep and monotonic, which also means the answer is insensitive to the factor-of-two uncertainty in the opacity — a factor of four error in L_diff moves T_c by about 20%.

Both conditions are solved by bracketed bisection with false-position acceleration, nested one inside the other. Bisection cannot run away, which matters when the inputs are allowed to be silly.

The cold branch

When nothing burns there is no thermal-equilibrium condition to impose, and the second condition has to come from somewhere else. What sets a non-burning object's central temperature is its contraction history: it heats as it shrinks and stops when degeneracy takes over the support. The solver takes the smaller of the virial estimate kT_c ≈ ½ μ m_H GM/R and the hottest T_c at which a hydrostatic solution of that mass exists at all. The second of those is the interesting one — it is the peak of contraction, found by walking T_c up until the mass minimum climbs past M and there is no solution left.

Failing to converge is physics

If the shooting method cannot close — no ρ_c gives the mass, no T_c balances the energy, the integration never reaches a surface — the program says NO HYDROSTATIC SOLUTION and stops. That is a result. It means the equations, with those inputs, do not describe a static self-gravitating object, and the honest thing to print is that they do not. It is not a black hole; nothing here computes collapse, and the word will never appear on a results panel.

In practice this branch is hard to reach inside the slider's range. Across a sweep of the whole input domain, essentially every combination produced some hydrostatic solution. That is itself worth knowing: a static ball of gas is easy to build, and what varies is whether it is a star.

5 · From diagnostics to an English verdict

Everything below reads the solved model. None of it can produce an outcome the integrator did not produce. The tests are applied in this order, and the first one that matches wins:

1. $2GM/R{c}^2\ge 1$ → BLACK HOLE — INSIDE ITS OWN HORIZON. The structure came back smaller than its own Schwarzschild radius, so it is not a structure.
2. ${\Gamma }_1<\frac{4}{3}+1.12(2GM/R{c}^2)$ → COLLAPSES TO A BLACK HOLE — GR INSTABILITY. See below.
3. L > L_Edd → NO EQUILIBRIUM — RADIATION DRIVEN. Its own luminosity exceeds 4πcGM/κ_es, so radiation pressure beats gravity in the envelope and the static solution is not one anything could hold.
4. Burning pp or CNO, X ≥ 0.1 → main-sequence analogue, with an MK letter and an H-burning lifetime.
5. Burning pp or CNO, X < 0.1 → hydrogen burning but hydrogen poor. Real burning, real lifetime, but no MK class: it is not on anybody's main sequence.
6. Burning 3α → helium burner. No MK class, no H-burning lifetime.
7. Burning carbon → carbon burner, flagged as the crudest rate in the program.
8. M < 0.0124 M_☉ (~13 M_J) → planetary-mass object. Below the deuterium limit, nothing burns at all.
9. Not burning, X ≥ 0.1 → brown dwarf. It has hydrogen and it never lit it; degeneracy stopped the contraction first.
10. Not burning, degenerate, M < M_Ch(μ_e) → white-dwarf-like. "Degenerate" means degenerate electrons carry more than 40% of the central pressure — a statement about pressure, not about η alone.
11. Not burning, X + Y < 0.05 → metal ball. No fuel, no spectral class.
12. Anything else that holds itself up but does not burn → cold hydrostatic object.

When it cannot hold together at all

Newtonian hydrostatics will support any mass you like if you make it hot enough — there is no Newtonian limit — so a solver that stops at dP/dr = −Gmρ/r² will never tell you anything has collapsed. Three checks are applied on top of the solution to catch what the integrator cannot see by itself.

Γ₁ runs from 5/3 for a gas-supported star down to 4/3 for a radiation-supported one, and 4/3 is the knife edge: squeeze a sphere with Γ₁ = 4/3 and it neither resists nor rebounds. General relativity raises that threshold by an amount proportional to the compactness (Chandrasekhar's coefficient, 1.12 for an n = 3 polytrope), so a star that is radiation-dominated and compact has nothing left holding it against radial collapse. This is the direct-collapse route to a black hole, and the object is still burning happily right up to the moment it goes — what the results panel shows is the last static model before the instability, not the remnant. Following it further needs hydrodynamics and general relativity, neither of which is here.

Composition decides where this happens, and it is worth seeing why. Radiation pressure takes over at lower mass when the mean molecular weight is high, so a ball of metal turns unstable far sooner than a ball of hydrogen. Running this solver up its own mass slider:

Those numbers are output, not input — nothing in the program knows about supermassive stars. Note also that the instability always arrives long before the horizon does. A C/O ball turns unstable at ~6 × 10³ M_☉, while 2GM/Rc² does not reach 1 until around 10⁹ M_☉ — five orders of magnitude later. The horizon test is therefore a guard that the reachable range never trips, and that ordering is the right way round: real supermassive stars collapse because Γ₁ fails them, not because they were ever compact.

Neutron stars: no, and here is why

This program cannot produce one, for two separate reasons, and it is worth keeping them apart.

The formation route is dynamical. Neutron stars are made in core collapse: a degenerate core grows past M_Ch, collapses on a free-fall timescale, electron capture drives the matter to neutrons, and the infall bounces off nuclear-density material. Every step of that is time evolution. This is a static snapshot with no time coordinate, so it cannot represent a process — only a configuration.

The equation of state has no neutrons. That second point is the sharper one, because a neutron star is a hydrostatic object, and a static solver with the right physics could perfectly well find one. This one cannot: it has ideal ions, ideal and degenerate electrons, and radiation. No neutron degeneracy pressure, no nuclear forces, no weak interactions. Its densest supported object is electron-degenerate and therefore capped at M_Ch.

So instead of inventing one, the solver reports the boundary. Above ρ_c ≈ 10¹⁰ kg/m³ electrons begin being captured onto nuclei and the mix stops being the one you typed in; the results panel flags that and says the description has stopped applying. And when a fuel-free object is handed a mass above M_Ch, the shooting method finds nothing, and the verdict says degeneracy cannot support it and names core collapse as the regime rather than guessing at the remnant.

In practice that neutronisation flag never fires either, and the reason is instructive. The densest thing this solver ever supports is a near-Chandrasekhar white dwarf at about 4 × 10⁹ kg/m³ — a factor of a few short. Anything denser would have to be a heavier degenerate object, and by then the carbon has lit. A genuinely inert composition, iron say, would let the white-dwarf branch run right up to M_Ch and past neutronisation; but the burning rates here treat every metal as carbon, so there is no inert mix to hand it.

The metal-ball branch turns out to be almost impossible to reach. Every fuel-free ball the solver could build in the slider's range was either degenerate enough to be called white-dwarf-like (below about 0.89 M_☉ at Z = 1) or contracted hot enough to light carbon (above it). The transition is sharp and it is the physics, not a threshold anyone typed in.

6 · Where the class, the burn temperature and the lifetime come from

MK class

Only for main-sequence analogues. T_eff is binned, and that is the whole method — no colours, no spectra, no lines:

Everything gets luminosity class V, because everything here is zero-age by construction. There are no giants in this program and there is no way to make one: a giant is a thing that happens to a star over time, and there is no time.

Burn temperature

T_c is reported as the burn temperature, because it is: it is the central temperature of the converged model, and the channel named beside it is whichever of the four rates is largest there. Both are read off the solution.

Hydrogen-burning lifetime

Fuel over burn rate, calibrated so the Sun comes out around 10 Gyr. It assumes L and the composition stay exactly as they are in this snapshot, which is false for every star that has ever existed — luminosity climbs as the core composition changes, so real main-sequence lifetimes are shorter than this. It is a factor-of-a-few estimate, not a death clock. If the model is not burning hydrogen, no lifetime is shown, because the formula would not mean anything.

The H-burning minimum mass

This is the one place a number from outside the integration decides which branch gets reported. The solver will happily find a thermal-equilibrium solution for a 0.05 M_☉ ball — such a state does satisfy the equations — but it is not a state contraction can reach, because degeneracy halts the object long before it gets there. Deciding that from first principles needs a brown-dwarf equation of state this program does not have, so the textbook value is used and scaled with the hydrogen content:

Less hydrogen, more mass needed, with a small nudge upward for a very metal-poor mix. Below M_min the cold structure is reported instead, and the object is called a brown dwarf. The 13 M_J deuterium limit is used the same way — as a label boundary, not as physics.

7 · Worked examples

Numbers below were read straight off the solver, not from a reference table.

8 · Lies we tell on purpose

In rough order of how much damage each one does:

• Transport is not solved. The temperature is tied to the density by a fixed polytrope rather than by radiative diffusion and convection. This is the big one, and it is why:
  • Low-mass stars come out badly wrong. Below about 0.4 M_☉ real stars are convective all the way down and their radii are set by the adiabat and the photosphere. Under an n = 3 stratification a 0.2 M_☉ model comes out roughly three times too large and far too faint, so its T_eff lands near 800 K instead of 3200 K. The results panel raises a note whenever you are in that regime. Do not trust the M dwarfs.
  • Massive stars come out too cool and too big — around 30 M_☉ by roughly 40% in radius and a factor of a few in luminosity, which costs about one spectral class. Their cores are convective and this model treats them as radiative.
  A Schwarzschild check runs on the finished model and reports what fraction of the mass is convectively unstable, so the model at least tells you when its own assumption is failing.
• No atmosphere. The integration stops where the pressure runs out; there is no photosphere and no τ = 2/3 surface. T_eff comes from $L=4\pi R^2\sigma T_{\text{eff}}^4$ after the fact. The outermost points of the T(r) plot are therefore meaningless — that curve is the polytrope running out of gas, not a real temperature structure, and the model's T_eff is nowhere on it.
• Kramers opacity. A power law standing in for an opacity table, with one uncertain constant tuned to the Sun. It has no H⁻, no molecular opacity, no conduction, and no bump anywhere it should have one. It diverges at low temperature where real Rosseland means saturate, which is a large part of why the low-mass end is so bad.
• Uniform composition. Real stars build a composition gradient the moment they start burning, and μ in the core is what drives most of main-sequence evolution. Here the mix is the same everywhere, forever, which is exactly true at zero age and false ever after.
• Crude ε. Four rates, no network, no screening corrections, no branching. Metals are treated as carbon for the burning and as CNO catalysts at Z/2 for the cycle. The carbon rate is a shape with a plausible ignition temperature and should be read as an on/off flag.
• Smooth-max electron pressure. Ideal and degenerate electron pressures are blended with (a⁴ + b⁴)^1/4 instead of a proper Fermi–Dirac integral. Correct in both limits, about 19% high right at the crossover — which is exactly where brown dwarfs live.
• ∇_ad ignores ionisation. The gas + radiation formula is used everywhere. Real partial-ionisation zones push ∇_ad down to 0.1 or so, which is a large part of why simple models cannot reach a sensible photosphere — and part of why this one does not try.
• Gravity is Newtonian. The collapse tests in §5 are relativistic corrections bolted onto a Newtonian solution, which is fine for deciding whether a configuration is stable and useless for describing one that is not. Once 2GM/Rc² is more than about 0.1 the reported radius and temperature should be read as order-of-magnitude at best, and the results panel says so.
• No neutrons. No neutron degeneracy pressure, no nuclear forces, no electron capture. The equation of state stops describing matter somewhere around 10¹⁰ kg/m³ and the program flags that rather than extrapolating quietly.
• No rotation, no magnetic fields, no mass loss, no binarity, no time. Not approximated: absent.

9 · Symbols