============================================================================== DAISYWORLD: A PLANET THAT REGULATES ITS OWN TEMPERATURE Science Journaling Club, Volume 1 Issue 4, Summer 2025 -- replication ============================================================================== This is a simulation. No planet, no daisy and no measurement of any physical object is involved anywhere in this file. Every number below is arithmetic performed by this program. python : 3.12.3 numpy : 2.4.2 master seed : 20250621 ============================================================================== [0] CONFIGURATION ============================================================================== Stefan-Boltzmann sigma = 5.670320e-08 W m^-2 K^-4 solar constant S at L = 1 = 917.0 W m^-2 albedo black/ground/white = 0.25 / 0.50 / 0.75 heat transport q = 2.060e+09 K^4 death rate gamma = 0.30 per unit time growth optimum T_opt = 295.50 K (22.35 C) growth half-width h = 17.50 K -> window 278.00 to 313.00 K seed floor per species = 0.010 RK4 step dt = 0.050, ceiling t = 200 steady-state tolerance = 1e-09 main grid 0.60 to 1.60 in 501 points (dL = 0.0020) scan grid 0.60 to 1.60 in 201 points (dL = 0.0050) MC grid 0.60 to 1.60 in 101 points (dL = 0.0100) Where q = 2.06e9 comes from. Watson and Lovelock write the local temperature in a linearised form, T_i = T_e + q'(A - A_i), with q' = 20 K per unit albedo. Differentiating T^4 gives the quartic version used here, with q = 4 T_opt^3 q'. So: club 4 * 295.50^3 * 20.0 = 2.064251e+09 K^4 published = 2.060000e+09 K^4 difference = 4.251e+06 K^4 (0.206%) The published constant is that arithmetic rounded to three figures. We keep 2.06e9 so that the replication uses the published number. ============================================================================== [1] VALIDATION: AREA CONSERVATION ============================================================================== The three areas are integrated as three independent state variables. Nothing in the integrator forces them to sum to one. If RK4 and the right-hand side are both correct, the sum is a conserved quantity and only floating-point rounding can move it. main sweep: 501 luminosities, 1162915 RK4 steps, 9.0 s luminosity steps that hit the integration ceiling : 124 of 501 worst |a_b + a_w + a_g - 1| at ANY of the 1162915 steps : 8.882e-15 worst at the 501 recorded steady states : 8.438e-15 float64 machine epsilon : 2.220e-16 club value 8.882e-15 required value 0 difference 8.882e-15 that is 40.0 machine epsilons, i.e. rounding and nothing else VERDICT: PASS, area conserved to machine precision Step-size check. The same sweep on a coarser grid at three values of dt. If the answer moves with dt, the integrator is the result. dt regulated L range drift(C) mean T(C) steps 0.200 0.705 to 1.555 10.23454 22.15211 119395 0.100 0.705 to 1.555 10.23454 22.15211 238626 0.050 0.705 to 1.555 10.23454 22.15211 477098 0.025 0.705 to 1.555 10.23454 22.15211 954072 ============================================================================== [2] VALIDATION: BARE-PLANET CONTROL vs ANALYTIC BLACK BODY ============================================================================== Run the same integrator with no daisies seeded at all. The answer must be the Stefan-Boltzmann temperature of a grey sphere of albedo 0.5, which we can write down in closed form: Te = ( S L (1 - A) / sigma )^(1/4) L club Te (K) analytic (K) difference (K) 0.60 263.91906 263.91906 0.000e+00 0.80 283.59949 283.59949 0.000e+00 1.00 299.86995 299.86995 0.000e+00 1.20 313.85442 313.85442 0.000e+00 1.40 326.18573 326.18573 0.000e+00 1.60 337.25853 337.25853 0.000e+00 max |club - analytic| over all 501 luminosities : 0.000e+00 K VERDICT: PASS bare planet warming across the full sweep: L = 0.60 -> -9.23 C ; L = 1.60 -> 64.11 C ; span 73.34 C best straight line: 71.859 C per unit L, worst residual 5.153 C the bare curve is L^(1/4), so 'linear' is an approximation good to 5.15 C over this range. We say near-linear and mean exactly that. ============================================================================== [3] THE MAIN RESULT: DOES IT REGULATE? ============================================================================== regulated interval (daisy cover > 0.01, longest contiguous run) L from 0.704 to 1.556, width 0.852 that is 85.2% of the swept range planetary temperature INSIDE the regulated interval min 19.158 C at L = 1.346 max 29.579 C at L = 1.556 mean 22.147 C residual drift (max - min) = 10.421 C daisy growth optimum for comparison = 22.35 C mean |regulated T - optimum| = 1.903 C standard deviation of regulated T = 2.246 C bare planet OVER THE SAME luminosity interval min 1.529 C, max 61.766 C, span 60.236 C club value vs control: regulated drift 10.421 C against bare span 60.236 C suppression factor = 5.78 (bare span divided by regulated drift) Watson & Lovelock 1983 state this qualitatively: temperature is held close to the optimum for growth over a wide range of luminosity, and their Figure 3 shows the two-species curve flat from roughly L = 0.7 to roughly L = 1.6. Our quantitative version: over L in [0.704, 1.556] the planet stays within 7.23 C of 22.35 C, while an identical bare planet moves through 60.2 C. The drift above is dominated by the two edges, where regulation is in the act of failing. Two narrower statements of the same result: |T - T_opt| < 2.0 C over L = 0.816 to 1.162, width 0.346 drift inside that band = 3.961 C bare planet span over the same L = 26.333 C suppression factor there = 6.65 middle 90% of the regulated interval (L = 0.746 to 1.512): drift 6.488 C against bare span 53.835 C, suppression 8.30 Both edges are jumps, not ramps. One grid step is dL = 0.0020. edge L before L after T before(C) T after(C) jump(C) cover before -> after cold 0.702 0.704 1.335 22.824 +21.489 0.00003 -> 0.69614 hot 1.556 1.558 29.579 61.873 +32.294 0.66492 -> 0.00000 largest single-step change anywhere: +32.294 C between L = 1.556 and 1.558 over that same step the bare planet changes by +0.1076 C The hot edge is the catastrophic collapse of life under gradual forcing described by Ackland, Clark and Lenton; the cold edge is the same discontinuity run the other way, colonisation rather than death. Phases of the sweep: black only : L 0.704 to 0.742 coexistence : L 0.744 to 1.340 white only : L 1.342 to 1.556 peak total cover 0.7000 at L = 1.480 bare ground never falls below 0.3000, because steady state requires a_g = gamma / beta >= gamma = 0.30. FIGURE DATA 1: L, planetary T (C), bare T (C), black, white, bare ground 0.600 -9.2309 -9.2309 0.00000 0.00000 1.00000 0.620 -7.0586 -7.0586 0.00000 0.00000 1.00000 0.640 -4.9382 -4.9382 0.00000 0.00000 1.00000 0.660 -2.8669 -2.8669 0.00000 0.00000 1.00000 0.680 -0.8421 -0.8421 0.00000 0.00000 1.00000 0.700 1.1384 1.1384 0.00000 0.00000 1.00000 0.720 24.2666 3.0770 0.68799 0.00001 0.31200 0.740 25.6445 4.9755 0.66929 0.00515 0.32556 0.760 25.3425 6.8360 0.62851 0.04496 0.32653 0.780 24.9583 8.6601 0.58892 0.08455 0.32653 0.800 24.5957 10.4495 0.55169 0.12178 0.32653 0.820 24.2529 12.2056 0.51661 0.15686 0.32653 0.840 23.9282 13.9299 0.48349 0.18998 0.32653 0.860 23.6203 15.6236 0.45219 0.22128 0.32653 0.880 23.3278 17.2881 0.42255 0.25092 0.32653 0.900 23.0497 18.9245 0.39444 0.27903 0.32653 0.920 22.7849 20.5337 0.36775 0.30572 0.32653 0.940 22.5324 22.1170 0.34237 0.33109 0.32653 0.960 22.2915 23.6752 0.31822 0.35525 0.32653 0.980 22.0614 25.2092 0.29520 0.37827 0.32653 1.000 21.8412 26.7199 0.27323 0.40024 0.32653 1.020 21.6305 28.2082 0.25224 0.42122 0.32653 1.040 21.4286 29.6747 0.23218 0.44129 0.32653 1.060 21.2350 31.1202 0.21298 0.46049 0.32653 1.080 21.0492 32.5454 0.19458 0.47889 0.32653 1.100 20.8706 33.9509 0.17694 0.49653 0.32653 1.120 20.6990 35.3374 0.16001 0.51346 0.32653 1.140 20.5338 36.7055 0.14375 0.52972 0.32653 1.160 20.3748 38.0556 0.12812 0.54535 0.32653 1.180 20.2215 39.3884 0.11309 0.56038 0.32653 1.200 20.0738 40.7044 0.09861 0.57486 0.32653 1.220 19.9313 42.0041 0.08466 0.58881 0.32653 1.240 19.7937 43.2878 0.07122 0.60225 0.32653 1.260 19.6607 44.5561 0.05825 0.61522 0.32653 1.280 19.5322 45.8094 0.04573 0.62774 0.32653 1.300 19.4080 47.0481 0.03364 0.63983 0.32653 1.320 19.2880 48.2726 0.02196 0.65152 0.32653 1.340 19.1760 49.4833 0.01072 0.66280 0.32648 1.360 19.2643 50.6804 0.00265 0.67295 0.32439 1.380 19.7596 51.8645 0.00031 0.68098 0.31871 1.400 20.4417 53.0357 0.00002 0.68739 0.31259 1.420 21.1926 54.1945 0.00000 0.69256 0.30744 1.440 22.0010 55.3411 0.00000 0.69649 0.30351 1.460 22.8749 56.4758 0.00000 0.69906 0.30094 1.480 23.8297 57.5989 0.00000 0.70000 0.30000 1.500 24.8920 58.7106 0.00000 0.69888 0.30112 1.520 26.1160 59.8114 0.00000 0.69478 0.30522 1.540 27.6450 60.9013 0.00000 0.68520 0.31480 1.560 61.9806 61.9806 0.00000 0.00000 1.00000 1.580 63.0496 63.0496 0.00000 0.00000 1.00000 1.600 64.1085 64.1085 0.00000 0.00000 1.00000 1.600 64.1085 64.1085 0.00000 0.00000 1.00000 FIGURE DATA 2: L, planetary albedo A, local T_black (C), local T_white (C) 0.600 0.50000 -2.4898 -16.5322 0.640 0.50000 1.4992 -11.8756 0.680 0.50000 5.3217 -7.4555 0.720 0.32800 25.7819 15.6375 0.760 0.35411 27.3384 17.3618 0.800 0.39252 27.3383 17.3617 0.840 0.42662 27.3383 17.3617 0.880 0.45709 27.3383 17.3617 0.920 0.48449 27.3383 17.3617 0.960 0.50926 27.3383 17.3617 1.000 0.53175 27.3383 17.3617 1.040 0.55228 27.3383 17.3617 1.080 0.57108 27.3383 17.3617 1.120 0.58836 27.3383 17.3617 1.160 0.60431 27.3383 17.3617 1.200 0.61906 27.3383 17.3617 1.240 0.63276 27.3383 17.3617 1.280 0.64550 27.3383 17.3617 1.320 0.65739 27.3385 17.3620 1.360 0.66757 27.5099 17.5516 1.400 0.67184 28.6736 18.8380 1.440 0.67412 30.1522 20.4693 1.480 0.67500 31.8548 22.3439 1.520 0.67369 33.9446 24.6389 1.560 0.50000 65.3501 58.5063 1.600 0.50000 67.4158 60.7009 1.600 0.50000 67.4158 60.7009 ============================================================================== [4] VALIDATION: STEADY-STATE ALGEBRA ============================================================================== At a steady state with a living population, growth must exactly balance death: a_g * beta(T_i) = gamma for every surviving type i. That is an algebraic identity we never imposed on the integrator. black daisies alive at 319 of 501 luminosities club |a_g beta(T_b) - gamma| max = 8.884e-05 required 0 median = 5.242e-09 max away from a live/dead edge = 5.059e-05 (315 points) white daisies alive at 407 of 501 luminosities club |a_g beta(T_w) - gamma| max = 3.951e-04 required 0 median = 1.917e-09 max away from a live/dead edge = 1.971e-04 (403 points) 60 of the 427 live luminosities carry a residual above 1e-6. The worst is at L = 0.746, where white cover is 0.015236 and falling. Our reading is that these are unfinished transients, not errors: a population on its way to zero approaches steady state slowly, and our integration window is t = 200. That is a testable claim. If it is right, the residuals shrink when the window is lengthened and the physical answer does not move. So we lengthened it fourfold. t = 200 t = 800 change worst |a_g beta - gamma| 3.951e-04 1.161e-06 340.38x luminosities above 1e-6 60 1 regulated interval low 0.704 0.704 regulated interval high 1.556 1.556 residual drift (C) 10.4210 10.4815 mean regulated T (C) 22.1472 22.1446 seconds 10.9 The residual falls and the physics does not move. We report the t = 200 numbers throughout, because that is the setting the whole study was run at, and note that the slow points are all dying populations near an edge. VERDICT: PASS at 367 of 427 live luminosities to 1e-6; the rest are transients that shrink by 340.4x when the window is quadrupled. ============================================================================== [5] VALIDATION: THE COEXISTENCE SYMMETRY, AN ANALYTIC PREDICTION ============================================================================== Where both colours coexist, both satisfy a_g beta(T) = gamma with the same a_g and the same gamma, so beta(T_b) = beta(T_w). The growth curve is a parabola about T_opt, and the two local temperatures are not equal, so they must sit symmetrically either side of the optimum: T_black + T_white = 2 T_opt = 591.00 K, exactly. Nobody told the integrator this. It is a consequence of the model and it is the sharpest single check available anywhere in the sweep. coexistence at 299 of 501 luminosities, L = 0.744 to 1.340 club mean (T_b + T_w) = 591.000425687 K analytic 2 T_opt = 591.000000000 K difference of means = 4.257e-04 K worst single luminosity = 3.763e-02 K (at a boundary) worst away from an edge = 1.885e-02 K (295 points) median over the band = 6.847e-07 K same quantity at t = 800 = 7.905e-05 K worst, 6.847e-07 K median VERDICT: PASS; the median point satisfies the identity to 7e-07 K and the worst points are the same edge transients as section [4]. The same algebra goes further, and predicts the numbers rather than a relation between them. Write T_b = T_opt + d and T_w = T_opt - d. Subtracting the two local-temperature equations kills T_e entirely: (T_opt + d)^4 - (T_opt - d)^4 = q (A_w - A_b) 8 T_opt^3 d + 8 T_opt d^3 = q (A_w - A_b) which is one cubic in d with no reference to luminosity, albedo, or the size of either population. Solve it and everything else follows: beta = 1 - (d/h)^2 for both colours, and a_g = gamma / beta. predicted d = 4.988282517 K predicted beta = 0.918749510 club analytic difference T_black (C) 27.338282837 27.338282517 3.209e-07 T_white (C) 17.361717839 17.361717483 3.551e-07 bare ground a_g 0.326530785 0.326530786 -1.732e-09 spread of a_g across the interior of the coexistence band: 1.287e-04 So the bare-ground fraction is pinned at 0.326531 for the whole coexistence band, independent of L. Everything the luminosity does is redistribute the remaining 0.673469 between black and white. VERDICT: PASS sample of the coexistence region: L T_black(K) T_white(K) sum(K) 2T_opt(K) diff(K) 0.744 300.49769 290.52213 591.019823 591.00 +1.982e-02 0.810 300.48828 290.51172 591.000001 591.00 +7.873e-07 0.876 300.48828 290.51172 591.000001 591.00 +5.543e-07 0.942 300.48828 290.51172 591.000001 591.00 +5.393e-07 1.008 300.48828 290.51172 591.000001 591.00 +5.137e-07 1.074 300.48828 290.51172 591.000001 591.00 +5.764e-07 1.140 300.48828 290.51172 591.000001 591.00 +7.278e-07 1.206 300.48828 290.51172 591.000001 591.00 +1.052e-06 1.272 300.48828 290.51172 591.000002 591.00 +1.895e-06 1.338 300.49129 290.51504 591.006330 591.00 +6.330e-03 ============================================================================== [6] THE CONTROL THE ARGUMENT NEEDS: EQUAL ALBEDOS ============================================================================== Set both daisy albedos equal to bare ground. The daisies still grow, compete and die, so the ecology is untouched. What is removed is the only channel by which they can influence temperature. If regulation survives this, the regulation was an artefact of something else. L flat-daisy Te(C) bare control(C) difference(K) daisy cover 0.60 -9.23094 -9.23094 -5.684e-14 0.00000 0.70 1.13841 1.13841 -5.684e-14 0.00000 0.80 10.44949 10.44949 0.000e+00 0.44192 1.00 26.71995 26.71995 -5.684e-14 0.68005 1.10 33.95092 33.95092 5.684e-14 0.46481 1.20 40.70442 40.70442 -1.137e-13 0.00000 1.40 53.03573 53.03573 -1.137e-13 0.00000 1.60 64.10853 64.10853 -1.137e-13 0.00000 club value 2.842e-13 K required value 0 K difference 2.842e-13 K temperature span of the equal-albedo planet : 73.339 C temperature span of the bare planet : 73.339 C daisies are alive at 188 of 501 luminosities, peak cover 0.7000 so the biosphere is present, large, and completely without effect. VERDICT: PASS, regulation vanishes entirely ============================================================================== [7] ROBUSTNESS 1: NARROWING THE GROWTH WINDOW ============================================================================== h is the half-width of the temperature range in which a daisy can grow at all. The published value is 17.5 K. Narrow it and the biota becomes fussier. How fussy can it be and still regulate? h(K) L_low L_high width drift(C) bare span(C) suppression peak cover 25.0 0.640 1.600 0.960 15.798 69.047 4.37 0.7000 22.5 0.660 1.600 0.940 14.807 66.975 4.52 0.7000 20.0 0.685 1.575 0.890 11.629 63.126 5.43 0.7000 17.5 0.705 1.555 0.850 10.235 60.085 5.87 0.7000 15.0 0.730 1.535 0.805 8.623 56.599 6.56 0.7000 12.5 0.755 1.520 0.765 7.544 53.437 7.08 0.7000 10.0 0.775 1.505 0.730 6.456 50.779 7.86 0.7000 7.5 0.800 1.490 0.690 4.613 47.707 10.34 0.7000 5.0 0.960 1.480 0.520 2.102 33.924 16.14 0.7000 4.0 0.970 1.475 0.505 1.670 32.874 19.69 0.6991 3.0 0.980 1.470 0.490 1.928 31.830 16.51 0.6964 2.0 0.990 1.460 0.470 2.254 30.508 13.54 0.6879 1.0 1.005 1.420 0.415 2.613 27.100 10.37 0.6472 FIGURE DATA 4: h, regulated width, drift, suppression 25.00 0.9600 15.7978 4.371 22.50 0.9400 14.8069 4.523 20.00 0.8900 11.6290 5.428 17.50 0.8500 10.2345 5.871 15.00 0.8050 8.6231 6.564 12.50 0.7650 7.5435 7.084 10.00 0.7300 6.4565 7.865 7.50 0.6900 4.6127 10.343 5.00 0.5200 2.1016 16.142 4.00 0.5050 1.6698 19.687 3.00 0.4900 1.9275 16.513 2.00 0.4700 2.2539 13.536 1.00 0.4150 2.6129 10.372 ============================================================================== [8] ROBUSTNESS 2: A DEATH RATE SHOCK ============================================================================== Raise gamma. Steady state needs a_g beta = gamma, and both a_g and beta are at most 1, so gamma = 1 is the hard ceiling at which no population can persist at any temperature whatsoever. gamma L_low L_high width drift(C) suppression peak cover 0.05 0.680 1.600 0.920 9.394 6.91 0.9461 0.10 0.685 1.600 0.915 8.753 7.36 0.8939 0.20 0.695 1.600 0.905 7.596 8.35 0.8000 0.30 0.705 1.555 0.850 10.235 5.87 0.7000 0.40 0.720 1.445 0.725 10.080 5.21 0.6000 0.50 0.730 1.350 0.620 9.652 4.77 0.5000 0.60 0.745 1.270 0.525 11.713 3.39 0.4000 0.70 0.765 1.200 0.435 13.924 2.40 0.3000 0.80 0.785 1.135 0.350 14.263 1.91 0.2000 0.90 0.815 1.085 0.270 17.907 1.18 0.1000 0.95 0.835 0.920 0.085 6.703 1.05 0.0500 0.99 extinct at every luminosity FIGURE DATA 5: gamma, regulated width, drift, suppression 0.05 0.9200 9.3944 6.914 0.10 0.9150 8.7527 7.364 0.20 0.9050 7.5959 8.355 0.30 0.8500 10.2345 5.871 0.40 0.7250 10.0795 5.213 0.50 0.6200 9.6519 4.771 0.60 0.5250 11.7128 3.393 0.70 0.4350 13.9236 2.399 0.80 0.3500 14.2628 1.911 0.90 0.2700 17.9065 1.180 0.95 0.0850 6.7029 1.049 0.99 0.0000 nan nan Now a transient shock instead of a permanent change: hold gamma at the published 0.3, sweep to L = 1.0, then kill a fraction of the standing crop instantly and integrate back to steady state. pre-shock state at L = 1.000: black 0.27323, white 0.40024, bare 0.32653 pre-shock temperature : 21.84125 C kill immediate T(C) recovered T(C) recovery error(C) steps 10.0% 22.3401 21.8412 -5.89e-07 1209 25.0% 23.0836 21.8412 -5.88e-07 1306 50.0% 24.3106 21.8412 -5.90e-07 1399 75.0% 25.5226 21.8412 -5.88e-07 1484 90.0% 26.2427 21.8412 -5.89e-07 1558 99.0% 26.6723 21.8412 -5.90e-07 1697 99.9% 26.7152 21.8412 -5.90e-07 1819 The planet returns to the same steady state after losing 99.9 per cent of its biota, because the attractor is a property of the equations rather than of the standing crop. ============================================================================== [9] MONTE CARLO OVER PARAMETER UNCERTAINTY ============================================================================== Every number so far comes from one point in parameter space, chosen by Watson and Lovelock in 1983. None of it is measured, because there is nothing to measure it against. So the useful question is not what the regulated width is, but how much it depends on choices nobody can check. We draw parameters at random and sweep each drawn planet. trials : 2000 stream : PCG64 spawned from SeedSequence(20250621), child 0 A_black ~ N(0.25, 0.05^2) clipped to [0.02, 0.48] A_white ~ N(0.75, 0.05^2) clipped to [0.52, 0.98] q ~ N(2.060e+09, 3.00e+08^2) clipped to [2e8, 6e9] gamma ~ N(0.30, 0.06^2) clipped to [0.05, 0.95] T_opt ~ N(295.50, 4.00^2) K h ~ N(17.50, 3.50^2) K clipped to [2, 35] vectorised sweep: 2000 planets x 101 luminosities, 151500 RK4 steps, 263.0 s trials with a living population somewhere : 2000 of 2000 (100.00%) trials measurable for drift (>= 2 points) : 2000 regulated luminosity width mean 0.8066 L SD 0.1124 SE 0.0025 n = 2000 median 0.8300 5th pct 0.5900 95th pct 0.9600 residual temperature drift mean 10.0103 C SD 4.1484 SE 0.0928 n = 2000 median 8.9830 5th pct 5.7564 95th pct 17.2258 suppression factor mean 6.4034 x SD 2.2989 SE 0.0514 n = 2000 median 6.2147 5th pct 3.1599 95th pct 9.9049 baseline at the published parameters, same 101-point grid: width 0.8500, drift 9.5085 C, suppression 6.240 baseline width vs Monte Carlo mean: 0.8500 vs 0.8066 +- 0.0025 -> 17.29 SE baseline drift vs Monte Carlo mean: 9.5085 vs 10.0103 +- 0.0928 -> -5.41 SE fraction of drawn planets with suppression > 2 : 0.9930 fraction with suppression > 5 : 0.7070 fraction with suppression > 10 : 0.0445 fraction with residual drift < 5 C : 0.0180 fraction with residual drift < 15 C : 0.9110 fraction that never regulate at all : 0.0000 CONVERGENCE: running mean and standard error as trials accumulate trials mean width SE mean drift(C) SE mean supp SE 10 0.7520 0.0404 9.2182 0.9027 6.265 0.688 25 0.7972 0.0282 8.6517 0.4493 6.805 0.359 50 0.8192 0.0171 9.1078 0.3740 6.701 0.247 100 0.8169 0.0119 9.2805 0.4043 6.773 0.196 200 0.8093 0.0079 9.5026 0.2843 6.658 0.147 400 0.8084 0.0057 9.6719 0.1888 6.543 0.107 700 0.8090 0.0043 9.8710 0.1511 6.462 0.080 1000 0.8071 0.0036 9.9318 0.1237 6.407 0.067 1400 0.8046 0.0030 9.9893 0.1104 6.402 0.063 1700 0.8065 0.0027 9.9814 0.1001 6.418 0.056 2000 0.8066 0.0025 10.0103 0.0928 6.403 0.051 FIGURE DATA 3: convergence, every 25th trial: n, running mean width, SE 25 0.79720 0.02815 50 0.81920 0.01709 75 0.80987 0.01429 100 0.81690 0.01190 125 0.81400 0.01072 150 0.81660 0.00927 175 0.81286 0.00840 200 0.80930 0.00791 225 0.80951 0.00733 250 0.80968 0.00698 275 0.80909 0.00677 300 0.80970 0.00651 325 0.81015 0.00614 350 0.81089 0.00598 375 0.81144 0.00579 400 0.80838 0.00565 425 0.80873 0.00546 450 0.80964 0.00529 475 0.81034 0.00515 500 0.81136 0.00500 525 0.81074 0.00487 550 0.80909 0.00478 575 0.80941 0.00469 600 0.80803 0.00464 625 0.80835 0.00454 650 0.80795 0.00444 675 0.80815 0.00437 700 0.80901 0.00427 725 0.80913 0.00420 750 0.80968 0.00415 775 0.80952 0.00408 800 0.81046 0.00399 825 0.81040 0.00393 850 0.81095 0.00386 875 0.81072 0.00381 900 0.81008 0.00378 925 0.80976 0.00373 950 0.80903 0.00368 975 0.80798 0.00366 1000 0.80706 0.00362 1025 0.80808 0.00358 1050 0.80790 0.00353 1075 0.80787 0.00350 1100 0.80756 0.00345 1125 0.80693 0.00340 1150 0.80632 0.00337 1175 0.80640 0.00333 1200 0.80612 0.00329 1225 0.80651 0.00327 1250 0.80657 0.00324 1275 0.80660 0.00319 1300 0.80623 0.00316 1325 0.80645 0.00311 1350 0.80608 0.00309 1375 0.80553 0.00306 1400 0.80461 0.00304 1425 0.80446 0.00301 1450 0.80434 0.00297 1475 0.80445 0.00295 1500 0.80435 0.00293 1525 0.80462 0.00290 1550 0.80521 0.00287 1575 0.80568 0.00284 1600 0.80523 0.00282 1625 0.80599 0.00280 1650 0.80595 0.00277 1675 0.80627 0.00275 1700 0.80648 0.00273 1725 0.80626 0.00271 1750 0.80665 0.00268 1775 0.80652 0.00266 1800 0.80647 0.00265 1825 0.80654 0.00264 1850 0.80586 0.00262 1875 0.80565 0.00261 1900 0.80575 0.00259 1925 0.80627 0.00257 1950 0.80633 0.00255 1975 0.80648 0.00253 2000 0.80657 0.00251 ============================================================================== [10] WHERE A DIFFERENT MODELLING CHOICE CHANGES THE ANSWER ============================================================================== Assumption 5 was that seeds are always present, so a dead population can return the moment conditions allow. Drop it: seed once at the cold end and never again. with re-seeding : regulated L 0.704 to 1.556, width 0.852, drift 10.421 C without re-seeding : regulated L 0.768 to 1.558, width 0.790, drift 10.908 C difference in width: 0.062 (7.3% of the re-seeded value) History dependence. Sweep the same grid downward from the hot end. descending: regulated L 0.622 to 1.212, width 0.590, drift 14.998 C max |T(ascending) - T(descending)| = 33.8202 C, at L = 1.498 hysteresis under this seeding rule : PRESENT, 33.82 C A third convention: seed every luminosity fresh from bare ground, with no memory of the previous step at all. fresh-seed: regulated L 0.704 to 1.212, width 0.508, drift 10.668 C against continuation width 0.852: the upper edge moves from L = 1.556 to L = 1.212, a loss of 40.4% of the regulated range. Reason: at high L a bare planet is already too hot for white daisies to germinate, so they can only be there if they were there already. Regulation at the hot end is inherited. Heat transport q, the softest parameter in the model: q (K^4) L_low L_high width drift(C) suppression 0.000e+00 0.770 1.140 0.370 28.953 1.00 2.000e+08 0.765 1.515 0.750 4.668 11.19 5.000e+08 0.755 1.520 0.765 5.644 9.47 1.000e+09 0.740 1.535 0.795 7.261 7.67 2.060e+09 0.705 1.555 0.850 10.235 5.87 4.000e+09 0.645 1.590 0.945 15.626 4.35 8.000e+09 1.020 1.600 0.580 8.423 4.26 q = 0 is a perfectly conducting planet: every patch sits at the planetary mean, black and white daisies feel identical temperatures, and the two species become one. Whatever survives there is the part of the result that owes nothing to having two colours. ============================================================================== [11] SUMMARY: CLUB VALUE BESIDE THE VALUE IT SHOULD MATCH ============================================================================== quantity club analytic difference area sum a_b+a_w+a_g, worst RK4 step 8.882e-15 0 8.882e-15 bare planet T vs Stefan-Boltzmann (K) 0.000e+00 0 0.000e+00 T_black + T_white in coexistence (K) 591.000187 591.000000 1.868e-04 bare ground in coexistence band 0.326530785 0.326530786 -1.732e-09 T_black in coexistence band (C) 27.338283 27.338283 3.209e-07 steady state |a_g beta - gamma|, interior 1.971e-04 0 1.971e-04 heat coefficient 4 T_opt^3 q' (K^4) 2.0643e+09 2.0600e+09 4.251e+06 equal-albedo control vs bare (K) 2.842e-13 0 2.842e-13 headline numbers for the article: regulated luminosity interval : 0.704 to 1.556 (width 0.852) residual temperature drift inside it : 10.421 C bare planet span over the same interval: 60.236 C suppression factor : 5.78 mean regulated temperature : 22.147 C growth optimum : 22.35 C bare planet span over the full sweep : 73.339 C tight band |T - T_opt| < 2 C : L 0.816 to 1.162 (width 0.346) drift inside the tight band : 3.961 C suppression inside the tight band : 6.65 largest single-step temperature jump : +32.294 C at L = 1.558 Monte Carlo regulated width : 0.8066 +- 0.0025 (SE, n = 2000) Monte Carlo residual drift : 10.0103 +- 0.0928 C Monte Carlo suppression factor : 6.403 +- 0.051 peak total daisy cover : 0.7000 at L = 1.480 ============================================================================== total wall clock: 416.0 s ==============================================================================