VPLanet: The Virtual Planet Simulator

TL;DR

VPLanet is a versatile, modular simulation software that models planetary system evolution over billions of years, integrating multiple physics modules to study habitability and system dynamics.

Contribution

It introduces a flexible, user-configurable framework for simulating diverse planetary and stellar phenomena without recompilation.

Findings

Validated by reproducing observations and past results

Supports large parameter sweeps for system analysis

Enables assessment of different physical processes' importance

Abstract

We describe a software package called VPLanet that simulates fundamental aspects of planetary system evolution over Gyr timescales, with a focus on investigating habitable worlds. In this initial release, eleven physics modules are included that model internal, atmospheric, rotational, orbital, stellar, and galactic processes. Many of these modules can be coupled simultaneously to simulate the evolution of terrestrial planets, gaseous planets, and stars. The code is validated by reproducing a selection of observations and past results. VPLanet is written in C and designed so that the user can choose the physics modules to apply to an individual object at runtime without recompiling, i.e., a single executable can simulate the diverse phenomena that are relevant to a wide range of planetary and stellar systems. This feature is enabled by matrices and vectors of function pointers that are dynamically allocated and populated based on user input. The speed and modularity of VPLanet enables large parameter sweeps and the versatility to add/remove physical phenomena to assess their importance. VPLanet is publicly available from a repository that contains extensive documentation, numerous examples, Python scripts for plotting and data management, and infrastructure for community input and future development.

Tables

Table 1. Table 1: BINARY Intermediate Values for Kepler-16

Parameter [units]	Leung & Lee (2013)	BINARY
$R_0$ [AU]	$0.7016$	$0.7016$
$n_k$ [yr-1]	$10.0823$	$10.0839$
$n_0/n_k$	$1.00702$	$1.00702$
${\kappa }_0/n_k$	$0.99224$	$0.99224$
${\nu }_0/n_k$	$1.02158$	$1.02158$
$C_0$	$0.000159$	$0.000160$
$C_1^0$	$-0.000282$	$-0.000282$
$C_2^0$	$-0.000589$	$-0.000589$
$C_3^0$	$-0.000049$	$-0.000049$
$C_1^{+}$	$0.000005$	$0.000005$
$C_2^{+}$	$-0.000033$	$-0.000033$
$C_3^{+}$	$-0.000006$	$-0.000006$
$C_1^{-}$	$0.035772$	$0.035776$
$C_2^{-}$	$0.002438$	$0.002440$
$C_3^{-}$	$0.000110$	$0.000110$

Table 2. Table 2: Properties of Selected Tidally Perturbed Worlds

Parameter	Gl 581	Gl 581 d	Jupiter	Ioa
$m_p$ ($M_{\oplus }$)	1.03 $\times {10}^5$	5.6	317.828	0.015
$r_p$ ($R_{\oplus }$)	31.6	1.6	11.209	0.286
$a$ ($R_{\oplus }$)	-	5123	-	66.13
$e$	-	0.0549	-	0.38
$P$ (d)	94.2	1	0.47
$ϵ$ (∘)	0	23.5	3.08a	0.0023
$r_g$	0.5	0.628	0.5	0.27
$k_2$	0.5	0.024059	0.3	1.5
$Q$	${10}^6$	100	${10}^5$	100

Table 3. Table 3: CoRoT-7 System Parameters

Parameter	Value
$M_{\star }$ [$M_{\odot }$]	0.93
$m_b$ [$M_{\oplus }$]	$8$
$m_c$ [$M_{\oplus }$]	$13.6$
$a_{b,\text{initial}}$ [AU]	$0.0188$
$a_{c,\text{initial}}$ [AU]	$0.0462$
$e_{b,\text{initial}}$	$0$
$e_{c,\text{initial}}$	$0.2$
$Q_b$	$100$
$Q_c$	$100$

Table 4. Table 4: Physical and orbital parameters for the system shown in Figs.  23 – 24 .

Parameter	Planet b	Planet c
$m_p$ ($M_{\oplus }$)	1	18
$r_p$ ($R_{\oplus }$)	1	1.5
$a$ (AU)	0.125	0.2246
$e$	0.1	0.1
$i$ (∘)	0.5	0.001
$\varpi$ (∘)	248.87	356.71
$\Omega$ (∘)	20.68	20
$ϵ$ (∘)	45	-
$p_A$ (∘)	0	-
$Q$	100	-
$k_2$	0.3	-

Table 5. Table 5: Disturbing function

Term
D0.1	$f_1+\left(h^2+k^2+h^{\prime 2}+k^{\prime 2}\right)f_2+\left(p^2+q^2+p^{\prime 2}+q^{\prime 2}\right)f_3$
	$+{\left(h^2+k^2\right)}^2{f}_4+\left(h^2+k^2\right)\left(h^{\prime 2}+k^{\prime 2}\right)f_5+{\left(h^{\prime 2}+k^{\prime 2}\right)}^2{f}_6$
	$+[(h^2+k^2)(p^2+q^2)+(h^{\prime 2}+k^{\prime 2})(p^2+q^2)$
	$+(h^2+k^2)(p^{\prime 2}+q^{\prime 2})+(h^{\prime 2}+k^{\prime 2})(p^{\prime 2}+q^{\prime 2})]f{}_7$
	$+\left[{\left(p^2+q^2\right)}^2+{\left(p^{\prime 2}+q^{\prime 2}\right)}^2\right]f_8+\left(p^2+q^2\right)\left(p^{\prime 2}+q^{\prime 2}\right)f_9$
D0.2	$(h{h}^{\prime }+k{k}^{\prime })[f_{10}+(h^2+k^2)f_{11}+(h^{\prime 2}+k^{\prime 2})f_{12}$
	$+(p^2+q^2+p^{\prime 2}+q^{\prime 2})f_{13}]$
D0.3	$(p{p}^{\prime }+q{q}^{\prime })[f_{14}+(h^2+k^2+h^{\prime 2}+k^{\prime 2})f_{15}$
	$+(p^2+q^2+p^{\prime 2}+q^{\prime 2})f_{16}]$
D0.4	$\left(h^2{h}^{\prime 2}-k^2{h}^{\prime 2}-h^2{k}^{\prime 2}+k^2{k}^{\prime 2}+4h{h}^{\prime }k{k}^{\prime }\right)f_{17}$
D0.5	$\left(h^2{p}^2-h^2{q}^2-k^2{p}^2+k^2{q}^2+4hkpq\right)f_{18}$
D0.6	$\left[h{h}^{\prime }\left(p^2-q^2\right)-k{k}^{\prime }\left(p^2-q^2\right)+2pq\left(h{k}^{\prime }+k{h}^{\prime }\right)\right]f_{19}$
D0.7	$\left(h^{\prime 2}{p}^2-h^{\prime 2}{q}^2-k^{\prime 2}{p}^2+k^{\prime 2}{q}^2+4h^{\prime }{k}^{\prime }pq\right)f_{20}$
D0.8	$\left(h^{2}p{p}^{\prime }-h^{2}q{q}^{\prime }-k^{2}p{p}^{\prime }+k^{2}q{q}^{\prime }+2hk{p}^{\prime }q+2hkp{q}^{\prime }\right)f_{21}$
D0.9	$\left[\left(h{h}^{\prime }+k{k}^{\prime }\right)\left(p{p}^{\prime }+q{q}^{\prime }\right)+\left(h{k}^{\prime }-k{h}^{\prime }\right)\left(p{q}^{\prime }-q{p}^{\prime }\right)\right]f_{22}$
D0.10	$\left[\left(h{h}^{\prime }+k{k}^{\prime }\right)\left(p{p}^{\prime }+q{q}^{\prime }\right)+\left(h{k}^{\prime }-k{h}^{\prime }\right)\left(q{p}^{\prime }-p{q}^{\prime }\right)\right]f_{23}$
D0.11	$\left[\left(h{h}^{\prime }-k{k}^{\prime }\right)\left(p{p}^{\prime }-q{q}^{\prime }\right)+\left(h{k}^{\prime }+k{h}^{\prime }\right)\left(p{q}^{\prime }+q{p}^{\prime }\right)\right]f_{24}$
D0.12	$\left(h^{\prime 2}p{p}^{\prime }-h^{\prime 2}q{q}^{\prime }-k^{\prime 2}p{p}^{\prime }+k^{\prime 2}q{q}^{\prime }+2h^{\prime }{k}^{\prime }{p}^{\prime }q+2h^{\prime }{k}^{\prime }p{q}^{\prime }\right)f_{25}$
D0.13	$\left(h^2{p}^{\prime 2}-h^2{q}^{\prime 2}-k^2{p}^{\prime 2}+k^2{q}^{\prime 2}+4hk{p}^{\prime }{q}^{\prime }\right)f_{18}$
D0.14	$\left[h{h}^{\prime }\left(p^{\prime 2}-q^{\prime 2}\right)-k{k}^{\prime }\left(p^{\prime 2}-q^{\prime 2}\right)+2p^{\prime }{q}^{\prime }\left(h{k}^{\prime }+k{h}^{\prime }\right)\right]f_{19}$
D0.15	$\left(h^{\prime 2}{p}^{\prime 2}-h^{\prime 2}{q}^{\prime 2}-k^{\prime 2}{p}^{\prime 2}+k^{\prime 2}{q}^{\prime 2}+4h^{\prime }{k}^{\prime }{p}^{\prime }{q}^{\prime }\right)f_{20}$
D0.16	$\left(p^2{p}^{\prime 2}-p^2{q}^{\prime 2}-q^2{p}^{\prime 2}+q^2{q}^{\prime 2}+4pq{p}^{\prime }{q}^{\prime }\right)f_{26}$

Table 6. Table 6: Stellar properties used in encounter model (reproduced from García-Sánchez et al., 2001 )

Stellar type	$\Delta M_V$	$v_{h,i}$ (km s-1)	${\sigma }_{\star i}$ (km s-1)	$n_{\star i}$ (10-3 pc-3)
B0	(-5.7, -0.2)	18.6	14.7	0.06
A0	(-0.2, 1.3)	17.1	19.7	0.27
A5	(1.3, 2.4)	13.7	23.7	0.44
F0	(2.4, 3.6)	17.1	29.1	1.42
F5	(3.6, 4.0)	17.1	36.2	0.64
G0	(4.0, 4.7)	26.4	37.4	1.52
G5	(4.7, 5.5)	23.9	39.2	2.34
K0	(5.5, 6.4)	19.8	34.1	2.68
K5	(6.4, 8.1)	25.0	43.4	5.26
M0	(8.1, 9.9)	17.3	42.7	8.72
M5	(9.9, 18)	23.3	41.8	41.55
WD	-	38.3	63.4	3.00
Giants	-	21.0	41.0	0.43

Table 7. Table 7: Parameters used in the EBM

Variable	Value	Units	Physical description
$C_L$	$1.55\times {10}^7$	J m-2 K-1	land heat capacity
$C_W$	$4.428\times {10}^6$	J m-2 K-1 m-1	ocean heat capacity per meter of depth
$m_d$	70	m	ocean mixing depth
$D$	0.58	W m-2 K-1	meridional heat diffusion coefficient
$\nu$	0.8		coefficient of land-ocean heat flux
$A$	203.3	W m-2	OLR parameter
$B$	2.09	W m-2 K-1	OLR parameter
${\alpha }_L$	0.363		albedo of land
${\alpha }_W$	0.263		albedo of water
${\alpha }_i$	0.6		albedo of ice
$f_L$	0.34		fraction of latitude cell occupied by land
$f_W$	0.66		fraction of latitude cell occupied by water

Table 8. Table 8: Parameters used in the ice sheet model

Variable	Value	Units	Physical description
$T_{freeze}$	273.15	K	freezing point of water
$L_h$	$3.34\times {10}^5$	J kg-1	latent heat of fusion of water
$r_{\text{snow}}$	$2.25\times {10}^{-5}$	kg m-2 s-1	snow/ice deposition rate
$A_{\text{ice}}$	$2.3\times {10}^{-24}$	Pa-3 s-1	deformability of ice
$n$	3		exponent of Glen’s flow law
${\rho }_i$	916.7	kg m-3	density of ice
${\rho }_s$	2390	kg m-3	density of saturated sediment
${\rho }_w$	1000	kg m-3	density of liquid water
$D_0$	$7.9\times {10}^{-7}$	s-1	reference sediment deformation rate
${\mu }_0$	$3\times {10}^9$	Pa s	reference sediment viscosity
$m$	1.25		exponent in sediment stress-strain relation
$h_s$	10	m	sediment depth
${\varphi }_s$	22	degrees	internal deformation angle of sediment
$T_b$	5000	years	bedrock depression/ rebound timescale
${\rho }_b$	3370	kg m-3	bedrock density

Table 9. Table 9: Initial Radiogenic Properties of Earth

Isotope	Half-life	Energy/decay	$Q_{rad,core}$(0)	$Q_{rad,man}$(0)	$Q_{rad,crust}$(0)
	(Gyr)	(${10}^{-12}$ J)	(TW)	(TW)	(TW)
26Al	0.717	0.642	0	0	0
40K	1.8178	0.213	33.86	36.16	13.89
232Th	20.202	6.834	0.145	6.52	3.83
235U	1.015	6.555	0.50	20.25	10.39
238U	6.452	8.283	0.12	11.67	5.16

Table 10. Table 10: ThermInt constants.

Symbol	Value	Units	Description
$A_{\nu }$	$3\times {10}^5$	J mol-1	Viscosity activation energy in (K6)
$A_{\mu }$	$2\times {10}^5$	J mol-1	Nominal shear modulus activation energy in (K7)
$A_{sol}$	$-1.160\times {10}^{-16}$	K$/$m3	Solidus coefficient in (K33)
$\alpha$	$3\times {10}^{-5}$	K-1	Thermal expansivity of mantle
${\alpha }_c$	$1\times {10}^{-5}$	K-1	Thermal expansivity of core
$B$	$2.5$	nd	Melt fraction coefficient in (K8)
$B_{sol}$	$1.708\times {10}^{-9}$	K$/$m2	Solidus coefficient in (K33), calibrated
$\beta$	$1/3$	nd	Convective cooling exponent in (K29)
${\beta }_{st}$	$1.71\times {10}^4$	GPa	Effective mantle stiffness
$c_m$	1265	J kg-1 K-1	Specific heat of mantle
$c_c$	840	J kg-1 K-1	Specific heat of core
$C_{sol}$	$-9.074\times {10}^{-3}$	K$/$m	Solidus coefficient in (K33), calibrated
$D$	2891	km	Mantle depth
$D_{Fe}$	7000	km	Iron solidus length scale
$D_N$	6340	km	Core adiabatic length scale
$D_{sol}$	$1.993\times {10}^4$	K	Solidus coefficient in (K33), calibrated
${\delta }_{ph}$	$6$	nd	Rheology phase coefficient in (K8, K9)
$E_G$	$3\times {10}^5$	J kg-1	Gravitational energy density release at ICB
${ϵ}_{UM}$	0.7	nd	Upper mantle adiabatic temperature drop
${ϵ}_{LM}$	1.3	nd	Lower mantle adiabatic temperature jump
${ϵ}_c$	0.8	nd	Average core to CMB adiabatic temperature drop
${\varphi }^{\ast }$	$0.8$	nd	Rheology phase coefficient in (K8, K9)
$g_{UM}$	$9.8$	${\text{m s}}^{-2}$	Upper mantle gravity
$g_{LM}$	$10.5$	${\text{m s}}^{-2}$	Lower mantle gravity
$g_c$	$10.5$	${\text{m s}}^{-2}$	CMB gravity
${\gamma }_c$	$1.3$	nd	Core Gruneisen parameter
${\gamma }_{dip}$	0.2	nd	Magnetic dipole intensity coefficient in (K39)
${\gamma }_{ph}$	$6$	nd	Rheology phase coefficient in (K8, K9)
$k_{UM}$	$4.2$	W m-1 K-1	Upper mantle thermal conductivity
$k_{LM}$	$10$	W m-1 K-1	Lower mantle thermal conductivity
$\kappa$	${10}^{-6}$	m2 s-1	Mantle thermal diffusivity
$L_{Fe}$	750	kJ kg-1	Latent heat of inner core crystallization
$L_{melt}$	320	kJ kg-1	Latent heat of mantle melting
$L_e$	$2.5\times {10}^{-8}$	W $\Omega$ K-1	Lorentz number
$M_m$	$4.06\times {10}^{24}$	kg	Mantle mass
$M_c$	$1.95\times {10}^{24}$	kg	Core mass
${\mu }_{ref}$	${10}^5$	Pa	Reference shear modulus in (K7)
${\mu }_0$	$4\pi \times {10}^{-7}$	H m-1	Magnetic permeability
${\nu }_{ref}$	$6\times {10}^7$	m2s-1	Reference viscosity
${\nu }_{LM}/{\nu }_{UM}$	$2$	nd	Viscosity jump from upper to lower mantle
$Q_{rad,0}$	$60$	TW	Initial mantle radiogenic heat flow (J07)
$R$	6371	km	Surface radius
$R_c$	3480	km	Core radius
$R_m$	4925	km	Radius to average mantle temperature $T_m$
$R{a}_c$	660	nd	Critical Rayleigh number
${\rho }_c$	11900	kg m-3	Core density
${\rho }_{ic}$	13000	kg m-3	Inner core density
${\rho }_m$	4800	kg m-3	Mantle density
${\rho }_{melt}$	2700	kg m-3	Mantle melt density
${\rho }_{solid}$	3300	kg m-3	Mantle upwelling solid density
$\Delta {\rho }_{\chi }$	700	kg m-3	Outer core compositional density difference
${\sigma }_c$	$10\times {10}^5$	S m-1	Core electrical conductivity
$T_{Fe,0}$	5600	K	Iron solidus coefficient in (K22)
$\xi$	$5\times {10}^{-4}$	nd	Rheology phase coefficient in (K8, K9)

Equations

$$T_{j}=\eta\times\min\Big{|}\frac{x_{j}}{dx_{j}/dt}\Big{|},$$

$$\Delta \alpha >\eta f_i{\left(\frac{d{f}_i}{d\alpha }\right)}^{-1},$$

$$\sin \Psi =\Vert \frac{(\mathbf{k}\times \mathbf{n})\times (\mathbf{s}\times \mathbf{n})}{|\mathbf{k}\times \mathbf{n}||\mathbf{s}\times \mathbf{n}|}\Vert ,$$

$$\cos \Psi =\frac{(\mathbf{k}\times \mathbf{n})\cdot (\mathbf{s}\times \mathbf{n})}{|\mathbf{k}\times \mathbf{n}||\mathbf{s}\times \mathbf{n}|}.$$

$$F_{\mathrm{EL}}=\frac{{ϵ}_{\mathrm{XUV}}{\mathcal{F}}_{\mathrm{XUV}}{R}_{\mathrm{p}}}{4G{M}_{\mathrm{p}}{K}_{\mathrm{tide}}{m}_{\mathrm{H}}},$$

$$R_{XUV}=\frac{R_p^2}{H\ln (p_{XUV}/p_s)+R_p}$$

$$F_{\mathrm{H}}=\{\begin{array}{ll}{F}_{\mathrm{EL}}& \text{if }{F}_{\mathrm{EL}}<F_{\mathrm{diff}}\\ F_{\mathrm{EL}}{\left(1+\frac{X_{\mathrm{O}}}{1-X_{\mathrm{O}}}\frac{m_{\mathrm{O}}}{m_{\mathrm{H}}}\frac{m_{\mathrm{c}}-m_{\mathrm{O}}}{m_{\mathrm{c}}-m_{\mathrm{H}}}\right)}^{-1}& \text{if }{F}_{\mathrm{EL}}\ge F_{\mathrm{diff}},\end{array}$$

$$F_{\mathrm{diff}}=\frac{(m_{\mathrm{O}}-m_{\mathrm{H}})(1-X_{\mathrm{O}})b_{\mathrm{diff}}g{m}_{\mathrm{H}}}{k_{\mathrm{boltz}}{T}_{\mathrm{flow}}}$$

$$m_{\mathrm{c}}=\frac{1+\frac{m_{\mathrm{O}}^2}{m_{\mathrm{H}}^2}\frac{X_{\mathrm{O}}}{1-X_{\mathrm{O}}}}{1+\frac{m_{\mathrm{O}}}{m_{\mathrm{H}}}\frac{X_{\mathrm{O}}}{1-X_{\mathrm{O}}}}{m}_{\mathrm{H}}+\frac{k_{\mathrm{boltz}}{T}_{\mathrm{flow}}{F}_{\mathrm{H}}^{\mathrm{ref}}}{\left(1+X_{\mathrm{O}}\left(\frac{m_{\mathrm{O}}}{m_{\mathrm{H}}}-1\right)\right)b_{\mathrm{diff}}g}$$

$$\begin{split}\Phi(R,\phi,z)&=\sum_{k=0}^{\infty}\Biggl{[}\Phi_{0k0}(R)-\frac{1}{2}\left(\frac{z}{R}\right)^{2}\Phi_{2k0}(R)+...\Biggr{]}\cos k(\phi-M_{b}-\varpi_{B})\\ &+e_{AB}\sum_{k=0}^{\infty}\Biggl{(}k\Biggl{[}\Phi_{0k0}(R)-\frac{1}{2}\left(\frac{z}{R}\right)^{2}\Phi_{2k0}(R)+...\Biggr{]}\\ &-\frac{1}{2}\Biggl{[}\Phi_{0k1}(R)-\frac{1}{2}\left(\frac{z}{R}\right)^{2}\Phi_{2k1}(R)+...\Biggr{]}\Biggr{)}\cos(k(\phi-\varpi_{B})-(k+1)M_{B})\\ &+e_{AB}\sum_{k=0}^{\infty}\Biggl{(}-k\left[\Phi_{0k0}(R)-\frac{1}{2}\left(\frac{z}{R}\right)^{2}\Phi_{2k0}(R)+...\right]\\ &-\frac{1}{2}\left[\Phi_{0k1}(R)-\frac{1}{2}\left(\frac{z}{R}\right)^{2}\Phi_{2k1}(R)+...\right]\Biggr{)}\cos(k(\phi-\varpi_{B})-(k-1)M_{B}),\end{split}$$

$$\Phi_{jk0}(R)=-\frac{2-\delta_{k0}}{2}\Biggl{[}(-1)^{k}\frac{m_{A}}{(m_{A}+m_{B})}b^{k}_{(j+1)/2}(\alpha_{A})+\frac{m_{B}}{(m_{A}+m_{B})}b^{k}_{(j+1)/2}(\alpha_{B})\Biggr{]}\frac{G(m_{A}+m_{B})}{R},$$

$$\Phi_{jk1}(R)=-\frac{2-\delta_{k0}}{2}\Biggl{[}(-1)^{k}\frac{m_{A}}{(m_{A}+m_{B})}\alpha_{A}Db^{k}_{(j+1)/2}(\alpha_{A})+\frac{m_{B}}{(m_{A}+m_{B})}\alpha_{B}Db^{k}_{(j+1)/2}(\alpha_{B})\Biggr{]}\frac{G(m_{A}+m_{B})}{R},$$

$${\alpha }_i=\frac{a_{AB}{m}_i}{R(m_A+m_B)},$$

$$\ddot{R}-R{\dot{\varphi }}^2=-\frac{\partial \Phi }{\partial R},R\ddot{\varphi }+2\dot{R}\dot{\varphi }=-\frac{1}{R}\frac{\partial \Phi }{\partial \varphi },\ddot{z}=-\frac{\partial \Phi }{\partial z}.$$

$$R=R_0+R_1(t),\varphi ={\varphi }_0(t)+{\varphi }_1(t),z=z_1(t),$$

$$n_k=\sqrt{G(m_A+m_B)/R_0^3}$$

$${\varphi }_0(t)=n_{0}t+{\psi }_0$$

$$n_0^2={\left[\frac{1}{R}\frac{d{\Phi }_{000}}{dR}\right]}_{R_0},$$

$$\begin{split}R&=R_{0}\Biggl{(}1-e_{\text{free}}\cos(\kappa_{0}t+\psi)-C_{0}\cos M_{\text{B}}-\sum_{k=1}^{\infty}\Biggl{[}C_{k}^{0}\cos k(\phi_{0}-M_{b}-\varpi_{B})\\ &+C_{k}^{+}\cos(k(\phi_{0}-\varpi_{B})-(k+1)M_{b})+C_{k}^{-}\cos(k(\phi_{0}-\varpi_{B})-(k-1)M_{b})\Biggr{]}\Biggr{)},\end{split}$$

$$C_{0}=-e_{AB}\left[\frac{d\Phi_{001}}{dR}\right]_{R_{0}}\Biggl{/}[R_{0}(k_{0}^{2}-n_{AB}^{2})],$$

$$C_{k}^{0}=\left[\frac{d\Phi_{0k0}}{dR}+\frac{2n\Phi_{0k0}}{R(n-n_{AB})}\right]_{R_{0}}\Biggl{/}(R_{0}[\kappa_{0}^{2}-k^{2}(n_{0}-n_{AB})^{2}]),$$

$$C_{k}^{\pm}=e_{AB}\left[\pm k\frac{d\Phi_{0k0}}{dR}-\frac{d\Phi_{0k1}}{2dR}+\frac{kn(\pm 2k\Phi_{0k0}-\Phi_{0k1})}{R(kn-(k\pm 1)n_{AB}}\right]_{R_{0}}\Biggl{/}(R_{0}[\kappa_{0}^{2}-(kn_{0}-(k\pm 1)n_{AB})^{2}]).$$

$${\kappa }_0^2={\left[R\frac{d{n}^2}{dR}+4n^2\right]}_{R_0}$$

$$\begin{split}\phi&=n_{0}t+\phi_{0}+\frac{2n_{0}}{\kappa_{0}}e_{\text{free}}\sin(\kappa_{0}t+\psi)+\frac{n_{0}}{n_{AB}}D_{0}\sin M_{b}\sum_{k=1}^{\infty}\Biggl{[}\frac{n_{0}}{k(n_{0}-n_{AB})}D_{k}^{0}\sin k(\phi_{0}-M_{b}-\varpi_{B})\\ &+\frac{n_{0}}{kn_{0}-(k+1)n_{AB}}D_{k}^{+}\sin(k(\phi_{0}-\varpi_{B})-(k+1)M_{b})+\frac{n_{0}}{kn_{0}-(k-1)n_{AB}}D_{k}^{-}\sin(k(\phi_{0}-\varpi)\\ &-(k-1)M_{b})\Biggr{]},\end{split}$$

$$D_0=2C_0,$$

$$D_k^0=2C_k^0-{\left[\frac{{\Phi }_{0k0}}{R^{2}n(n-n_{AB})}\right]}_{R0},$$

$$D_{k}^{\pm}=2C_{k}^{\pm}-e_{AB}\Biggl{[}\frac{k(\pm 2k\Phi_{0k0}-\Phi_{0k1})}{2R^{2}n(kn-(k\pm 1)n_{AB})}\Biggr{]}_{R_{0}}.$$

$$z=R_0{i}_{\text{free}}\cos ({\nu }_{0}t+\zeta ),$$

$${\nu }_0^2={\left[-\frac{{\Phi }_{200}}{R^2}\right]}_{R_0}.$$

$$h$$

Full text

1 Introduction

Exoplanetary systems display a diversity of morphologies, including a wide range of orbital architectures and planetary densities. This heterogeneity likely leads to a wide range of evolutionary trajectories that result in disparate planetary properties. As astronomers and astrobiologists probe these worlds to determine their atmospheric and surface properties, a comprehensive model of the physical effects that sculpt a planet can help prioritize targets and interpret observations. Here we describe a software package called VPLanet that self-consistently simulates many processes that influence the evolution of gaseous and terrestrial planets in a range of stellar systems. Our approach allows for coupling of simple models by simultaneously solving ordinary and partial differential equations (ODEs and PDEs) and explicit functions of time to track the evolution of, and feedbacks among, interior, atmospheric, stellar, orbital, and galactic processes. Below we describe a set of physical models (called modules) that simulates these phenomena, as well as their assumptions and limitations. We validate the code, including both individual modules and a subset of module combinations, by reproducing key observations of the Earth, Solar System bodies, and known stellar systems. Where observations are lacking or unavailable, we reproduce a selection of previously published results. The software is open source and includes documentation, examples, and the opportunity for community involvement.111https://github.com/VirtualPlanetaryLaboratory/vplanet

While VPLanet is designed to model an arbitrary planetary system, the primary motivation for creating VPLanet is to investigate the potential habitability of exoplanets with a single code that can capture feedbacks across the range of physical processes that affect a planet’s evolution. TRAPPIST-1 (Gillon et al., 2016, 2017) is a good example of why coupled processes are needed to understand planetary evolution. In that system, the host star dimmed during the pre-main sequence, orbital interactions between planets are strong, tidal heating and rotational braking are significant, and stellar activity can drive atmospheric mass loss. As discussed below, VPLanet combines a set of theoretical models to provide a first order approach to investigate these processes and the feedbacks among them.

We define habitable to mean a planet that supports liquid surface water. Since astrochemical and planet formation studies find that most small planets form with significant inventories of water and bioessential elements (e.g., van Dishoeck et al., 2014; Morbidelli et al., 2018), the most pertinent question for potentially habitable planets may be “do they still have water?” The presence of water is controlled by the planet’s interior, atmosphere, orbit, host star, and even the galaxy. Our approach does not presume a priori that any one process dominates the evolution because numerous processes can severely impact a planet’s potential to support liquid water, even in the habitable zone (HZ; Kasting et al., 1993; Kopparapu et al., 2013). For example, tidal heating may be strong enough to drive a runaway greenhouse (Barnes et al., 2013), the pre-main sequence evolution of low-mass stars may desiccate (remove water from) a planet (Luger & Barnes, 2015), and stable resonant orbital oscillations can drive extreme eccentricity cycles (Barnes et al., 2015). Furthermore, all these processes can operate in a given system, and hence a rigorous model of planetary system evolution, including habitability, should include as broad a range of physics as possible (Meadows & Barnes, 2018).

Coupled evolutionary models could be valuable tools for interpreting the environments of newly discovered Earth-sized planets with nearly Earth-like levels of incident stellar radiation (“instellation”) orbiting other stars. These planets may be capable of supporting liquid water, but their habitability is currently an open question given the plethora of physical processes at play, i.e., a planet’s presence in the star’s habitable zone does not imply the planet supports liquid water. Moreover, worlds like Proxima Centauri b (Anglada-Escudé et al., 2016) and TRAPPIST-1 c–g (Gillon et al., 2017) orbit bright enough host stars that their atmospheres will be probed with future ground- and space-based facilities (Meadows et al., 2018; Lincowski et al., 2018), which will provide constraints on otherwise unobservable surface habitability. This paper describes a benchmarked model of planetary system evolution that can be used to simulate the evolution of $\sim 1{M_{\oplus}}$, $\sim 1{R_{\oplus}}$ terrestrial planets with approximately Earth-like material properties and structure, with the goal of understanding their potential for habitability.

Though we are interested in Earth-like planets, the model is not restricted to Earth-radius and Earth-mass planets; rather, the model is intended to investigate a wide range of planetary systems. The phenomena described above are also important for uninhabitable worlds such as GJ 1132 b (Berta-Thompson et al., 2015) and 40 Eri b (Ma et al., 2018). VPLanet can be used to simulate many such planets (and moons) to infer their histories.

The individual modules of VPLanet can simulate many aspects of planetary evolution, but the ability to couple modules together facilitates novel investigations, as demonstrated in several previously-published studies. Deitrick et al. (2018a, b) combined the orbital, rotational, and climate modules to show that potentially habitable exoplanets can become globally glaciated if their orbital and rotational properties evolve rapidly and with large amplitude. Fleming et al. (2018) showed that the coupled stellar-tidal evolution of tight binary stars can lead to orbital evolution that ejects circumbinary planets. Lincowski et al. (2018) combined stellar evolution and atmospheric escape to track water loss and oxygen build-up on the TRAPPIST-1 planets. Finally, Fleming et al. (2019) simulated the coupled stellar-tidal evolution of binary stars to show how unresolved binaries impact gyrochronology age estimates of stellar populations.

VPLanet is a fast and flexible code that combines a host of semi-analytic models, which are all written in C. It can provide quick insight into an individual planetary system with a single simulation (such as calculating tidal heating or stellar evolution), or can perform parameter sweeps and generate ensembles of evolutionary tracks that can be compared to observations. Alternatively, VPLanet can be combined with machine learning algorithms to identify key parameters (Deitrick et al., 2018b). These capabilities can provide direct insight or can complement research with more sophisticated tools, such as 3-D global circulation models, which are too computationally expensive to explore vast parameter spaces. For example, VPLanet can be used to isolate the most important phenomena (see $\S$ 2) and explore parameter space, and then more complicated models can target interesting regions of that parameter space to provide further insight and observational predictions.

The modularity and flexibility of VPLanet is enabled by matrices and vectors of function pointers, in which individual elements represent addresses of functions. In this framework, a user can specify a range of physics to be simulated, i.e., select modules at runtime, and the software dynamically allocates the memory and collates the appropriate derivatives for integration. This approach allows users to trivially add or remove physics, e.g., tides or stellar evolution, and assess their relative importance in the system’s evolution. Moreover, this approach allows a “plug and play” development scheme in which new physics can be added and coupled with minimal effort.

The VPLanet code and its ecosystem have been designed to ensure transparency, accessibility, and reproducibility. For example, figures presented below that are derived from VPLanet output include a link (in the electronic version of the paper) to the location in the VPLanet repository that contains the input files and scripts that generate the figure, e.g., examples. The approximate run time for the simulation(s) necessary to generate the figure is also listed. The software is open source and freely available for all to use.

The objectives of VPLanet are to 1) simulate newly discovered exoplanets to assess the probability that they possess surface liquid water, and hence are viable targets for biosignature surveys, 2) model diverse planetary and stellar systems, regardless of potential habitability, to gain insights into their properties and history, and 3) enable transparent and open science that contributes to the search for life in the universe. This paper is organized as follows. In $\S$ 2 we describe the VPLanet algorithm and underlying software that enables flexibility in how the modules are connected. Readers interested only in the validation of the physics can skip to $\S$ 3 in which we begin the qualitative description of the 11 fundamental modules and demonstrate reproducibility of previously published results. Details describing the quantitative results and implementation are relegated to Appendices A–K. Briefly, the 11 modules, in alphabetical order, are as follows: simple thermal atmospheric escape models with AtmEsc ($\S$ 3, App. A), an analytic model of circumbinary planet orbital evolution with BINARY ($\S$ 4, App. B), 2nd and 4th order secular models of orbital evolution with DistOrb ($\S$ 5, App. C), a semi-analytic model for rotational axis evolution due to orbital evolution and the stellar torque with DistRot ($\S$ 6, App. D), an approximate model for tidal effects with EqTide ($\S$ 7, App. E), a model of Oort Cloud object orbits adopted to capture wide binary orbits that includes galactic migration with GalHabit, ($\S$ 8, App. F), an energy balance climate model with an explicit treatment of ice sheet growth and retreat with POISE ($\S$ 9, App. G), radiogenic heating throughout planetary interiors with RadHeat ($\S$ 10, App. H), an $N$-body orbital model with SpiNBody ($\S$ 11, App. I), stellar evolution, including the pre-main sequence, with STELLAR ($\S$ 12, App. J), and an internal thermal and magnetic evolution model that is calibrated to Earth and Venus with ThermInt ($\S$ 13, App. K). In $\S$ 14 we reproduce previous results that couple multiple modules. In $\S$ 15 we discuss the value of the coupled model and how it may be used to prioritize targets for life-detection observations. Appendix L is a list of all symbols arranged by module, and Appendix M describes VPLanet’s accessibility as well as customized tools that streamline its usage.

2 The VPLanet Algorithm

Models of planetary system evolution must be both comprehensive enough to simulate a planetary system with an arbitrary architecture, as well as flexible enough to ensure that only appropriate physics are applied to individual system members. The VPLanet approach is to engineer a software framework in which the user selects established models, and the executable assembles the appropriate subroutines to calculate the evolution. After checking for consistent input, the code simulates the entire system by solving the relevant equations simultaneously. This approach provides a simple interface and, more importantly, the opportunity to isolate important processes by easily turning on/off certain physics without needing to recompile. For example, the role of radiogenic heating in a planet’s evolution can be toggled by simply removing the module name (RadHeat) from the appropriate line in the input file.

The key software design that permits this flexibility is the use of arrays and matrices of function pointers, i.e., the elements contain memory addresses of functions. With this approach VPLanet reads in the modules requested by the user and assembles the appropriate governing equations into a multi-dimensional matrix in which one dimension corresponds to the bodies, the second to variables, and the third to processes that modify the variables, see Fig. 1. Parameters that are directly calculated are called “primary variables” and all other parameters are derived from them.

VPLanet is written in C, whose standard version includes the function pointer matrices described above. Moreover, C generates an executable that is usually as fast or faster than any other computer language, and speed of calculation is critical for the high-dimensional problem of planetary system evolution. Note that all modules are in C, but the support scripts, described in more detail in Appendix M, are written in python.

The modules described in the next 11 sections all solve ODEs or are explicit functions of time (with the exception of POISE, which solves PDEs in latitude and time). Thus, to couple the modules together and simulate diverse phenomena, VPLanet loops over the dimensions of the function pointer matrix and solves the equations simultaneously. For parameters influenced by multiple processes, VPLanet calculates the derivatives and then sums them to obtain the total derivative.

As an example, consider a system consisting of a star, a close-in planet (i.e., with an orbital period less than 10 days), and a more distant planetary companion. The star undergoes stellar evolution, the inner planet and the star experience mutual tidal effects, and the two planets gravitationally perturb each other. If the user applies the STELLAR and EqTide modules to the star, the EqTide, DistOrb, and DistRot modules to the inner planet, and the DistOrb and DistRot modules for the outer planet, then VPLanet will simulate the coupled stellar-orbital-tidal-rotational evolution of the 3 objects. In this case, the eccentricity of the inner planet is modified by tidal effects and gravitational perturbations from the other planet, so in the matrix the inner planet’s dimension for eccentricity processes would have 2 members for the two derivatives.

VPLanet integrates the system of equations using a fourth order Runge-Kutta method with dynamical timestepping. At each timestep, the values of all the derivatives are used to calculate a timescale for each process on each primary variable:

[TABLE]

where $T_{j}$ is the the timescale of the $j$th process, $\eta$ is a number less than 1 that is tuned to provide the desired accuracy, and $x_{j}$ is a primary variable, e.g., the star’s luminosity or a planet’s obliquity. All primary variables are updated over the same timescale $T_{j}$, i.e., the timestep is set by the fastest-changing variable. In many cases some variables are updated at a faster cadence than necessary for convergence, but this methodology ensures that all effects are properly modeled. Note that the timestep is computed from individual derivatives, not the summed derivatives, to ensure that all phenomena are accurately modeled. In the two-planet system described above, the derivatives for eccentricity from both tides and perturbations are calculated and included in Eq. 1. This approach is similar to the stand-alone code EqTide 222https://github.com/RoryBarnes/EqTide (Barnes, 2017), which served as the foundation for VPLanet.

The VPLanet operational flow chart is shown in Fig. 2. VPLanet first reads in the options set by the user in $\geq 2$ input files. The primary input file contains the top-level instructions, such as integration parameters (if any), units, and the list of the members of the system. Each system member has a “body file” which contains all the initial conditions, module-specific parameters, and output option selections for that object. After reading in all the options, they are vetted for completeness and inconsistencies in a process called “verify.” After passing this step, the physical state and evolution may be self-consistently calculated.

After verification, a log file can be written in which the initial state of the system is recorded in SI units, which are the units used in all internal VPLanet calculations, although we note that the user can specify both the input and output units in the primary input file. At this point, if the user requested an integration, it begins. The evolution is broken up into four parts: 1) calculation of “auxiliary properties,” i.e., parameters needed to calculate the derivatives of the primary variables, 2) one step (forward or backward) is taken, 3) the state of the system is written to a file if requested, 4) if necessary, changes to the matrix are implemented in “force behavior” (e.g., if a planet becomes tidally locked, its rotational frequency derivative due to tides is removed), and 5) checks for any threshold the user set to halt the integration are performed (e.g., the user may specify that the execution should be terminated once all water is lost from a planet). Once the integration is complete, the final conditions of the system can be recorded in the log file.

Note that in addition to coupling modules through simultaneous solutions of ODEs and PDEs, coupling also can occur in planetary interiors. For example, in a tidally heated terrestrial planet, the tidal power is a function of temperature (e.g., Driscoll & Barnes, 2015), but the temperature is also a function of the tidal power. VPLanet couples these connections during the Auxiliary Properties step in Fig. 2.

In the following sections, we briefly describe and validate each of the 11 modules that control the evolution of planetary systems in VPLanet. More detailed information about each module can be found in the Appendices.

3 Atmospheric Escape: AtmEsc

The erosion of a planet’s atmosphere due to extreme stellar radiation is among the biggest challenges to its habitability, particularly around low mass stars (e.g., Lissauer, 2007; Scalo et al., 2007; Luger & Barnes, 2015). The AtmEsc module models the escape of planetary atmospheres and their surface volatiles following simple parametric energy-limited and diffusion-limited prescriptions, which are discussed in detail in Appendix A. In order to validate our approach, in Figure 3 we present a reproduction of Figure 3 in Lopez & Fortney (2013) using AtmEsc. In that study, the authors used a coupled thermal evolution / photoevaporation model to explain the density dichotomy in the Kepler-36 system, which hosts two highly irradiated planets: a low-density mini-Neptune and a high-density super-Earth. The authors showed how, at fixed instellation, the initial core mass dictates the evolution of the gaseous envelope of a sub-giant planet: planets with high-mass cores hold on to their envelopes, while those with low-mass cores are more easily stripped by photoevaporation. As in Lopez & Fortney (2013), we plot the percentage of b’s gaseous envelope that is lost as a function of the core mass (top panel), initial envelope mass (center panel), and initial total mass (bottom panel) for different initial envelope mass fractions, assuming an escape efficiency $\epsilon_{XUV}=0.1$, a planet age of 5 Gyr, instellation one hundred times that received by Earth, and the XUV evolution model of Ribas et al. (2005) for a solar-mass star. We model the planet’s radius with the evolutionary tracks for super-Earths of Lopez et al. (2012) and Lopez & Fortney (2014).

We find, as those authors did, that the core mass displays the tightest correlation with the fraction of the envelope that is lost. Our values generally agree, although AtmEsc predicts slightly higher escape rates at lower core mass. For instance, in the top panel, AtmEsc finds that planets with core masses up to $\sim 3-4\mathrm{M}_{\oplus}$ are completely stripped for all values of the envelope fraction; in the top panel of Figure 3 in Lopez & Fortney (2013), only planets with masses less than $\sim 2-3\mathrm{M}_{\oplus}$ are completely evaporated. This result is likely due to the fact that the study of Lopez & Fortney (2013) employs a fully coupled thermal evolution/atmospheric escape model, in which the simulation tracks the evolution of the internal entropy and luminosity of the planet and any feedbacks with the escape process. In contrast, AtmEsc makes use of pre-computed radius grids as a function of mass and irradiation to calculate the atmospheric escape rate, and so is unable to capture feedback effects from the escape on the thermal evolution. However, at higher core mass the difference between the two studies is on the order of tens of percent, well within the uncertainties in parameters like $\epsilon_{XUV}$ and other observational constraints.

4 Circumbinary Planet Orbits: BINARY

The recent discovery of transiting circumbinary planets (CBPs) by Kepler provides intriguing laboratories to probe the orbital dynamics of such systems. Non-axisymmetric gravitational perturbations from the central binary force CBPs into oscillating orbits that display short and long-term non-Keplerian behavior. The BINARY module computes the orbit of a massless test particle on a circumbinary orbit using the analytic theory derived by Leung & Lee (2013). By generalizing the work of Lee & Peale (2006) to the case of an eccentric central binary orbit, Leung & Lee (2013) modeled the orbit of a massless CBP as a combination of the circular motion of the CBP’s guiding center with radial and vertical epicylic oscillations induced by non-axisymmetric components of the central binary’s gravitational potential. We discuss the full model and implementation details of the Leung & Lee (2013) theory in App. B.

Leung & Lee (2013) validated their analytic formalism against direct N-body simulations, and here we validate our implementation of their analytic theory by reproducing their Figure 4 that depicts the orbital evolution of the CBP Kepler-16 b (Doyle et al., 2011). For the initial conditions, we used the orbital parameters for both the binary and CBP given in Table 1 in Leung & Lee (2013) and set the CBP’s $e_{\text{free}}=0.03$ following Leung & Lee (2013). The results of our validation simulation are shown in Fig. 4 and are in excellent agreement with Figure 4 from Leung & Lee (2013). We further validate BINARY by comparing the intermediate quantities used in the Leung & Lee (2013) theory (see Appendix B) for the Kepler-16 system in Table 1. We find that BINARY perfectly reproduces nearly all of the intermediate quantities from Leung & Lee (2013). The maximum error between the BINARY values and Leung & Lee (2013) is $0.02\%$ on $n_{k}$, a negligible difference.

5 Approximate Orbital Evolution: DistOrb

An approximate solution to the orbital evolution of a planetary system can be derived from a quantity known as the disturbing function, the non-Keplerian component of the gravitational potential in a multi-body system. The disturbing function is most useful when written as a Fourier expansion in the orbital elements. First derived by Lagrange and Laplace (see, for example, chapter 7 of Murray & Dermott, 1999), this approach produces ODEs for the evolution of orbital parameters. Outside of mean motion resonances, an orbit-averaged, or secular, disturbing function can be used, which has the advantage that large time steps (hundreds of years) can be taken. The DistOrb module in VPLanet is based on the fourth order (in eccentricity and inclination) secular solution derived in Murray & Dermott (1999) and Ellis & Murray (2000). The theory is described in detail in Appendix C; here, we present results from the model.

Figures 5 and 6 show the orbital evolution of the inner and outer solar system planets, respectively, as calculated by DistOrb and HNBODY. The latter software package is an N-body code that calculates the gravitational evolution from first principles (Rauch & Hamilton, 2002). For the DistOrb runs, we used the fourth-order integration (a second order Lagrange-Laplace solution can also be utilized in the code; see Appendix C). Here we compare to HNBODY because that model also contains the general relativistic corrections, however, for the solar system these effects are small and HNBODY results appear almost identical to other N-body integrators.

The inclination evolution in DistOrb compares extremely well with HNBODY. The eccentricity evolution compares reasonably well for most of the planets; the largest error is in the amplitude of Mercury’s eccentricity variation. Figure 7 shows the absolute errors in the eccentricity and inclination evolution between the two models. The errors are largest for Mercury and Mars after $\sim 3$ Myr. The errors in $e$ and $i$ for Mercury, Mars, and Uranus grow in time. For Earth and Venus, the error also grow in time but remain smaller. The errors for Jupiter, Saturn, and Neptune are periodic and are explained easily by a slight mis-match in frequencies—this leads to a drift in the relative phase between the solutions, producing errors that are periodic and stable.

We do not expect to perfectly reproduce the N-body solution with this secular model, as the Solar System is affected by the proximity of Jupiter and Saturn to a 5:2 mean-motion resonance (Lovett, 1895). The other source of error is Mercury’s relatively large eccentricity ($e\sim 0.2$), which the fourth-order model does not handle as well as the direct N-body solution. For the Solar System, DistOrb performs as well as previous studies (Murray & Dermott, 1999). Since the orbital elements of the planets are known with a high degree of precision, the N-body solution is clearly desirable. However, for exoplanetary systems, for which the errors in $e$ and $i$ are quite large, DistOrb offers a computational advantage: the above simulation of eight planets runs in $\sim 200$ s on a modern CPU, about 10 times faster than HNBODY. This speed allows for broader explorations of parameter space. For systems which are near to mean-motion resonances or which have tight constraints on the orbital parameters (as in our Solar System), SpiNBody may be preferred over DistOrb.

6 Rotational Evolution from Orbits and the Stellar Torque: DistRot

There are a number of physical processes that affect the position of a planet’s spin axis. The module DistRot captures the physics of two processes: the torque acting on the equatorial bulge by the host star, and the motion of the planet’s orbital plane (e.g., the change in inclination). The model was derived in Kinoshita (1975) and Kinoshita (1977); see Appendix D for details.

Figure 6 shows the obliquity evolution for Earth and Mars using DistRot. We show the obliquity for Earth over the next Myr with and without the effect of the Moon, which compares extremely well with Laskar et al. (1993), Fig. 11. For the case including the Moon, the relative error between our solution and that of Laskar et al. (1993) is $\sim 0.6\%$ on average, with a maximum of $\sim 2\%$. For the case without, the relative error is $\sim 5\%$ on average, with a maximum $\sim 13\%$ (on the final point at 1 Myr). The larger error for the latter comes from a slight mis-match in the frequencies and the larger amplitude of the oscillation. Note that we do not directly include the effect of the Moon—here, the effect is mimicked by forcing Earth’s precession rate to the known value, $50.290966^{\prime\prime}/$yr (Laskar et al., 1993). This forced precession rate can be modified to recreate the effects of an arbitrary moon.

Mars is more challenging because its obliquity is known to be very sensitive to orbital frequencies (Ward & Rudy, 1991; Ward, 1992; Touma & Wisdom, 1993; Laskar et al., 2004). We show the obliquity evolution utilizing two methods for the orbital evolution over the last 10 Myr. The first couples DistRot directly to DistOrb. In this case, the obliquity evolution over the last $\sim 5$ Myr compares well with previous studies (e.g., Touma & Wisdom, 1993, Fig. 1), but it does not contain the well known shift to a higher obliquity state at $\sim 5$ Ma, owing to an imperfect representation of the orbital frequencies (see below). Still, we match the $\sim 20^{\circ}$ oscillation with periods of $\sim 120$ kyr, $\sim 1.25$ Myr, and $\sim 2.5$ Myr (Ward, 1992; Touma & Wisdom, 1993). In the second case, we use the orbital evolution from HNBODY as input into DistRot. This second case does produce the $5$ Ma obliquity shift, indicating that the problem lies with the accuracy of the orbital model, not with DistRot itself.

7 Tidal Effects: EqTide

Tidal effects modify a planet’s orbit, rotational properties, and internal power. The EqTide module employs the equilibrium tide models originally developed by Darwin (1880). We specifically use formulations called the constant-phase-lag (CPL) and constant-time-lag (CTL) models developed by Ferraz-Mello et al. (2008) and Leconte et al. (2010), respectively. See Appendix E for more details on these models.

In Fig. 9, we show the rotational evolution of the putative exoplanet Gl 581 d (Udry et al., 2007), which may in fact be an artifact and not actually a planet (Robertson et al., 2014). We nonetheless consider this example as it was examined by Heller et al. (2011) and thus provides a straightforward validation of EqTide. Physical and orbital parameters are listed in Table 2 (Udry et al., 2007). Fig. 9 is very similar to Fig. 6 in Heller et al. (2011) in which CPL is labeled “FM08” and CTL is “Lec10,” with differences less than 1%, most likely due to updates of fundamental constants (Prša et al., 2016).

In Fig. 10 we show the tidal heating surface flux of Io as a function of $e$ and $\varepsilon$. The current heat flux is 1.5 – 3 W/m*-2* (Veeder et al., 1994, 2012), and the physical and orbital parameters of Jupiter and Io are listed in Table 2. Io’s orbital eccentricity has been damped to 0.004, and is likely in a Cassini state, (see $\S$ 14.4 and Bills & Ray, 2000), with an obliquity of 0.0023∘, a displacement that remains below the detection threshold. The predicted heat flow of Io is about 2–3 times higher than observed, which has led some researchers to speculate that Io’s heat flow is not in equilibrium (Moore, 2003), although it is also likely that the equilibrium tide formalism does not accurately reflect Io’s tidal response. Nonetheless, VPLanet successfully reproduces Io’s surface energy flux.

8 Galactic Evolution: GalHabit

The GalHabit module accounts for two effects of the galactic environment on the orbits of binary star systems: the galactic tide and perturbations from passing stars. Such effects have been shown to impact the stability of planetary systems (Kaib et al., 2013) and so should be considered when studying planets in wide binary systems. We utilize the secular approach of Heisler & Tremaine (1986) for the galactic tide and model stellar encounters following the Monte-Carlo formulations in Heisler et al. (1987) and Rickman et al. (2008).

We model the evolution of an M dwarf orbiting the sun, with a mass of $0.12$ M*⊙*, $a=10000$ au, $e=0.7$, and $i=80^{\circ}$. This approach is similar to the simulation used in Figure 1 of Kaib et al. (2013). Two examples are shown: the evolution of the orbit under the galactic tide alone and the evolution with both the tide and random stellar encounters. We cannot reproduce the evolution shown in Kaib et al. (2013) exactly because our stellar encounters are drawn randomly. Nonetheless, we qualitatively recover the behavior of such a system. The timescale of tidal evolution is approximately 3 Gyr, similar to the simulation shown in Figure 1 of Kaib et al. (2013). The order of magnitude of the change in semi-major axis due to impulses is on the order of $100-1000$ au, also similar to Kaib et al. (2013).

We have tested the impulse approximation against an HNBODY simulation with the Bulirsch-Stoer integration method in the case above (errors are expected to increase with decreasing semi-major axis), for a total of 337,235 comparison simulations. The scale of the change in the periastron distance is on the order of 100–1000 AU, which agrees well with the errors found by using the impulse approximation by Rickman et al. (2005).

9 Climates of Habitable Planets: POISE

To model climate, VPLanet utilizes an energy balance model (EBM) appropriate for Earth-like atmospheres. POISE is based on North & Coakley (1979) who solved for the temperature and albedo as a function of latitude and day of year and forced by the incoming instellation. An additional component, a model for the dynamics of ice sheets on land, is based on Huybers & Tziperman (2008). The instellation is calculated directly from the planet’s orbital and rotational parameters, and the luminosity of the host star.

In Figure 12 we show the seasonal cycle for Earth as modeled by POISE. Here, we have included ice sheet growth on land, but sea ice is modeled simply as an increase in surface albedo when temperatures are below freezing. We have compared the outgoing longwave radiation and absorbed stellar fluxes to satellite data of Earth (Barkstrom et al., 1990); the full comparison is shown in Deitrick et al. (2018b). We find good agreement between the model and data for the OLR, with the largest errors ($\sim 20$%) occurring near the poles and the tropics. The error at these locations is due to the relatively coarse nature of the model, which is unable to capture the effects of atmospheric circulation and water vapor there. The error for the absorbed stellar radiation is similar over the equator and mid-latitudes ($\leq 20$%), though becomes much larger toward the poles. This is primarily a result of the simple parameterization of the albedo, which does not take into account the additional contribution from the atmosphere, and does not capture albedo variations with longitude. Overall, POISE performs well despite the model’s simplified nature. The ice mass balance shows the sum of annual accumulation and melting, which equals ice flow convergence for a stable ice sheet. Here, the negative values represent potential melting—this value is calculated even in the absence of ice. POISE can simulate an ice-covered, ice-free, or partially ice-covered planet, and variability among all three.

The hemispheric asymmetry (most easily seen in the left-hand panels of Fig. 12) is a result of Earth’s eccentricity and the relative angle between perihelion and equinox. Since perihelion occurs near the beginning of the calendar year, southern summer is more intense than northern, and a greater amount of ice melt can occur.

10 Radiogenic Heating: RadHeat

The major sources of radiogenic heating are included in the module RadHeat, specifically 235U, 238U, 232Th, 40K, and 26Al. On Earth, and presumably other planets, these isotopes are distributed through planetary cores, mantles, and crusts. 26Al is relatively short-lived, but can provide enormous power on planets that form within a few Myr. The current amount of radiogenic power in Earth is poorly constrained (see, e.g., Araki et al., 2005), and VPLanet assumes the current total power is 24.3 TW (Turcotte & Schubert, 2002). The evolution of these reservoirs on Earth, described in more details in $\S$ 13 and App. K is shown in Fig. 13. For more information on this module, consult Appendix H.

11 Accurate Orbital Evolution: SpiNBody

The VPLanet N-body model, SpiNBody, directly calculates the gravitational acceleration between massive objects. As this calculation is from first principles, SpiNBody is valid for any configuration. Figure 14 reproduces Fig. 13 from Barnes et al. (2015). This system consists of a solar-type star, an Earth-mass planet in the HZ, and outer Neptune-mass planet in a 3:1 resonance. The eccentricity occasionally surpasses 0.9999, but the system is stable for 10 Gyr. (Note that the N-body model assumes point masses, so the planet and star don’t merge in this example, which is chosen to validate an extreme case.) The bottom right panel of Fig. 14 shows that energy and angular momentum are conserved to within 1 part in $10^{9}$, similar to MERCURY (Chambers, 1999), and validates this VPLanet module. A simulation of the Solar System is also presented in Appendix I.

12 Stellar Evolution: STELLAR

Evolving stellar parameters shape the dynamics and habitability of stellar and exoplanetary systems. For example, exoplanets orbiting in the habitable zone of late M-dwarfs likely experienced an extended runaway greenhouse during the host star’s superluminous pre-main sequence phase, potentially driving extreme water loss (e.g., Luger & Barnes, 2015). The combination of evolving stellar radii and magnetic braking, the long-term removal of stellar angular momentum arising from the coupling of the stellar wind with the surface magnetic field, dictate the stellar angular momentum budget, molding observed stellar rotation period distributions as a function of stellar age and mass (e.g., Skumanich, 1972; McQuillan et al., 2014; Matt et al., 2015). Moreover, in stellar binaries, coupled stellar-tidal evolution depends sensitively on the evolving stellar radii and rotation state, driving orbital circularization on the pre-main sequence (e.g., Zahn, 1989), potentially destablizing any circumbinary planets they may harbor (e.g., Fleming et al., 2018), and can strongly impact the stellar rotation period evolution (e.g., Fleming et al., 2019).

VPLanet’s stellar evolution module, STELLAR, tracks the evolution of the fundamental stellar parameters of low-mass (M_{\star}~{}\hbox to0.0pt{<\hss}{\lower 4.30554pt\hbox{\sim}}1.4\mathrm{M_{\odot}}) stars, including a star’s radius, radius of gyration, $r_{g}$, effective temperature, luminosity, XUV luminosity, and rotation rate. STELLAR models stellar evolution via a bicubic interpolation of the Baraffe et al. (2015) models of solar metallicity stars over mass and time. STELLAR computes the XUV luminosity according to the product of the luminosity and Eq. (J3) for stars that drive atmospheric escape and water loss (see $\S$ 3; Luger & Barnes, 2015). Furthermore, STELLAR allows the user to model the long-term angular momentum evolution of low-mass stars using one of three magnetic braking models (Reiners & Mohanty, 2012; Repetto & Nelemans, 2014; Matt et al., 2015, 2019). Here, we demonstrate the modeling capabilities of STELLAR, while in Appendix J, we describe the numerical and theoretical details of the module.

We demonstrate the evolution predicted by STELLAR in Figure 15 that depicts the evolution of the stellar radius, $r_{g}$, temperature, luminosity, and XUV luminosity, the latter computed using Eq. (J3) and assuming $f_{sat}=10^{-3}$, $t_{sat}=1$ Gyr, and $\beta_{XUV}=-1.23$ (Ribas et al., 2005), all as functions of time for several stars ranging in mass from late M dwarfs ($M_{\star}=0.08\mathrm{M_{\odot}}$) to late F dwarfs ($M_{\star}=1.3\mathrm{M_{\odot}}$). Our tracks agree with present-day solar values and display the well-known extended pre-main sequence phase of M dwarfs (e.g., Luger & Barnes, 2015).

13 Geophysical Evolution: ThermInt

The thermal and magnetic evolution of the interior of rocky Earth-like planets is modeled in ThermInt by solving the coupled heat balance in the mantle and core. Parameterized heat flow scalings and material properties appropriate for Earth are assumed. A nominal Earth-model is calibrated to approximately reproduce the main constraints on the thermal and magnetic evolution: present-day mantle potential temperature of 1630 K, mantle surface heat flow of 38 TW, upwelling mantle melt fraction of 7%, inner core radius of 1221 km, and a continuous core magnetic field. The radiogenic abundances assumed in the mantle are based on bulk silicate Earth models (Arevalo et al., 2013; Jaupart et al., 2015) and 3 TW of K40 in the core today to maintain a continuous dynamo (Driscoll & Bercovici, 2014). Using the default Earth values, ThermInt + RadHeat produce an “Earth interior model” with the following values at 4.5 Gyr: $Q_{surf}=33.4(\pm 5)$ TW, $Q_{cmb}=13.5(\pm 4)$ TW, $T_{UM}=1587(-34,+161)$ K, $T_{CMB}=4000(\pm 200)$ K, $R_{ic}=1224$ km, $M=80$ ZAm2, and inner core nucleation at $3.97$ Gyr ($0.53$ Ga). The error bars in Figure 16 reflect the range of values in Jaupart et al. (2015). For the nominal Earth interior model, Figure 16 shows the thermal evolution of the mantle and core temperatures, heat flows, thermal boundary layer thicknesses, mantle viscosities, upwelling mantle melt fraction, and mantle melt mass flux over time. Figure 17 shows the thermal and magnetic evolution of the core in terms of inner core radius, core buoyancy fluxes, geodynamo magnetic moment, and magnetopause radius over time. There are no error bars in Fig. 17 because $R_{ic}$ is from the “preliminary reference Earth model (PREM Dziewonski & Anderson, 1981) and the other quantities (core buoyancy fluxes, magnetic moment, and magnetopause radius) are calibrated to give Earth-like values. In both figures, the evolution concludes at the modern Earth and all parameters are consistent with observations.

14 Multi-Module Applications

In this section we present results in which the previous modules are coupled together to reproduce previously published results.

14.1 Milankovitch Cycles

In this subsection we demonstrate VPLanet’s ability to reproduce Milankovitch cycles on Earth. This section utilizes DistOrb, DistRot, and POISE. Note that in order to reproduce the effect of Earth’s moon on Earth’s obliquity, we force the precession rate to be $50.290966^{\prime\prime}$ year*-1* (Laskar et al., 1993). This choice does not perfectly match the dynamics of the Earth-Moon-Sun system, but it is close enough to replicate the physics of the ice age cycles. The surface properties are shown in Figure 18 (see Huybers & Tziperman, 2008, Figure 4 for comparison), for a 200,000 year window. The ice sheets in the northern hemisphere high latitude region grow and retreat as the obliquity, eccentricity, and climate-precession-parameter, or CPP ($e\sin{(\varpi+\psi)}$), vary. The ice deposition rate is less than that used by Huybers & Tziperman (2008) and so the ice accumulation per year is slightly smaller. The ice ablation occurs primarily at the ice edge (around latitude $60^{\circ}$) and is slightly smaller than Huybers & Tziperman (2008), peaking at $\sim 2.1$ m yr*-1*, compared to their $\sim 3$ m yr*-1*.

We note that in this framework, the net growth and retreat of ice sheets is highly sensitive to the tuneable ice deposition rate, $r_{\text{snow}}$. With $r_{\text{snow}}$ too low, we do not build up ice caps on Earth at all. With $r_{\text{snow}}$ too high, the ice sheets grow so large that they become insensitive to orbital forcing. With $r_{\text{snow}}\sim 2.25\times 10^{-5}$ kg m*-2* s*-1*, we roughly reproduce the Earth’s ice age cycles at $\sim 40,000$ years and $\sim 100,000$ years over a 10 million year simulation.

Note that Huybers & Tziperman (2008) used a deposition (precipitation) rate of 1 m per year, which is indeed roughly the global average for Earth. However, precipitations rates are diminished in cold regions, which suggests that a diminished value for ice deposition is justified. Additionally, we find in VPLanet that when a value of 1 m per year is used, the ice sheets become too thick and are insensitive to orbital forcing. Note as well that Huybers & Tziperman (2008)’s climate (EBM) components are coded on different grids and utilize different parameterizations for radiative transfer.

There are a number of differences between our reproduction of Milankovitch cycles and those of Huybers & Tziperman (2008). Most notably, our ice sheets tend to persist for longer periods of time, taking up to three obliquity cycles to fully retreat. As previously stated, we also require a lower ice deposition (snowing) rate than Huybers & Tziperman (2008) in order to ensure a response from the ice sheets to the orbital forcing. We attribute these differences primarily to the difference in energy balance models used for the atmosphere. For example, our model has a single-layer atmosphere with a parameterization of the OLR tuned to Earth, while Huybers & Tziperman (2008) used a multi-layer atmosphere with a simple radiative transfer scheme. Further, while the Huybers & Tziperman (2008) model contained only land, our model has both land and water which cover a fixed fraction of the surface. The primary effect of having an ocean in this model is to change the effective heat capacity of the surface. This dampens the seasonal cycle and affects the ice sheet growth and retreat. Thus, our seasonal cycle is somewhat muted compared to theirs, and our ice sheets do not grow and retreat as dramatically on orbital time scales. Ultimately, our ice age cycles are more similar to the longer late-Pleistocene cycles than to $\sim 40,000$ year cycles of the early-Pleistocene.

Even though we cannot perfectly match the results of Huybers & Tziperman (2008), we find the comparison acceptable. Both models make approximations to a number of physical processes and thus have numerous parameters that have to be tuned to reproduce the desired behavior. Thus, it is no surprise that we do not reproduce their results precisely. Furthermore, despite all of the crude assumptions made in the energy balance model (the parameterization of radiative transfer, the reduction of the sphere to a single dimension, etc), we are nevertheless able to produce the ice age frequencies at $\sim 40$ kyr and $\sim 100$ kyr, suggesting that our simple model captures the basic physical processes involved. Finally, the details of Earth’s ice ages (such as the dominance of the $100$ kyr cycle in the late-Pleistocene), remain difficult to capture and explain, even with sophisticated models (Abe-Ouchi et al., 2013; Sánchez Goñi et al., 2019).

14.2 Evolution of Tight Binary Stars

Stars form large and contract until the central pressure become large enough for fusion. For stars in tight binary systems, large tidal torques rapidly circularize the orbit early in the system’s history. In a classic study of the evolution of short-period binary stars, Zahn & Bouchet (1989) found that orbits of binaries with an orbital period of less than about 8 days tidally circularized before they reach the Zero Age Main Sequence (ZAMS), consistent with contemporaneous observations. More recent surveys have found that circularization extends to binary orbital periods of about 10 days (e.g., Meibom & Mathieu, 2005; Lurie et al., 2017). Here we reproduce the Zahn & Bouchet result by coupling the stellar evolution models of Baraffe et al. (2015), incorporated in the STELLAR module, and the equilibrium tide CPL model via EqTide. See Sections 7 and 12 and their corresponding Appendices, E and J, respectively, for a more in-depth discussion of those models.

Adopting the initial conditions used to produce Figure 1 from Zahn & Bouchet (1989), we model an equal mass $1{M_{\odot}}-1{M_{\odot}}$ binary with initial orbital eccentricity $e=0.3$ and an orbital period of 5 days. The initial stellar rotation rate to mean motion ratio is set at $\Omega/n=3$ in line with the estimates from Zahn & Bouchet (1989), who assumed conservation of angular momentum during the stellar accretion phase. For our tidal model, we set $Q=1.25\times 10^{5}$ and $k_{2}=0.5$, both reasonable values for stars given the wide range of assumed values in the literature (e.g., Barnes et al., 2013; Fleming et al., 2018). The parameters $Q$ and $\tau$ can, and probably do, vary as a function of the forcing frequency (see e.g., Penev et al., 2018), but VPLanet does not (yet) include this complication.

The results of the simulation are depicted in Figure 19. Our results are in good agreement with Zahn & Bouchet (1989), Fig. 1. In this case, a quantitative comparison is unwarranted as the Zahn & Bouchet (1989) model used older stellar evolution models, so we do not expect the results to be a perfect match. After an initial increase in orbital eccentricity, the binary circularizes within the first ${\sim}10^{6}$ years before the ZAMS, in agreement with Zahn & Bouchet (1989). The transition between increasing and decreasing eccentricity occurs when $e=\sqrt{1/19}$ at the $\Omega/n=1.5$ transition, as expected from the CPL model (see Section E.1). The orbital period peaks at $t=10^{5}$ years and then decreases as the orbit circularizes. One difference between the two model predictions is that Zahn & Bouchet (1989) find an increase in $\Omega/n$ from unity to over 2 near $10^{6}$ years, before the stars tidally lock again after about $10^{9}$ years. In our model, the stars remain tidally locked after $10^{5}$ years as we force stars with rotation periods close to the orbital period to remain tidally locked to prevent numerical instabilities (see the appendix of Barnes et al. (2013) for a more in-depth discussion of this numerical necessity). This tidal locking formalism can still model complex physical interactions near the tidally locked state (e.g. subsynchronous rotation, Fleming et al., 2019), but it can easily be disabled.

14.3 Interiors of Tidally Heated Planets

In this multi-module application we model the gravitational tidal dissipation in the interior of an Earth-like planet and its orbit. The modules used are ThermInt, RadHeat, and EqTide. The tidal dissipation equations used in this application is the “orbit-only” model from Driscoll & Barnes (2015), see also Appendix E.3, and the dissipation efficiency depends on the temperature of the mantle. To reproduce the results of Driscoll & Barnes (2015) the dissipation efficiency in the orbital equations ($\dot{e}$ and $\dot{a}$) is approximated by $Im(k_{2})\approx k_{2}/\mathcal{Q}$, where $\mathcal{Q}=\eta\omega/\mu$ is the Maxwell tidal efficiency.

This example reproduces the results of Driscoll & Barnes (2015) for three tidally evolving planets orbiting a $0.1$ solar-mass star to within a few percent. (Note that some of the underlying physics in their interior model has been updated in ThermInt, so we do not expect an exact match.) Figures 20 and 21 compare the thermal, magnetic, and orbital evolution of three Earth-like planets each with the same initial eccentricity of 0.5 and initial Semi-major axes of $0.01$, $0.02$, $0.05$ all in AU. Figure 20 reproduces Driscoll & Barnes (2015), Fig. 5 and Fig. 21 reproduces Driscoll & Barnes (2015), Fig. 4. The main features that are reproduced from Driscoll & Barnes (2015) are that the tidal power peaks at $T_{UM}=1800$ K, the orbits circularize around $10^{-4}$ Gyr for $e_{0}=0.01$, $10^{-2}$ Gyr for $e_{0}=0.02$, and $>10$ Gyr for $e_{0}=0.05$, the inner core nucleates around 2 Gyr, and mantle melt mass flux goes to zero around 5 Gyr.

14.4 Tidal Damping in Multi-Planet Systems

In multi-planet systems in which one or more planets is close enough to the host star for tides to damp the orbit, orbital and rotational properties can reach an approximately fixed state as tides remove energy and angular momentum from planetary orbits and rotations. In this subsection, we reproduce the damping of orbital evolution into the “fixed point solution” (Wu & Goldreich, 2002; Zhang & Hamilton, 2008) and the damping of rotational cycles into a Cassini state (Colombo & Shapiro, 1966; Ward & Hamilton, 2004; Brasser et al., 2014; Deitrick et al., 2018a).

One difficulty arises from the fact that the semi-major axis decays, and its evolution is not accounted for in DistOrb, which ignores terms involving the mean longitude. The functions, $f_{i}$, in the disturbing function (Table 5) are computationally expensive. In the secular approximation used in DistOrb, the semi-major axes of the planets do not change, thus the $f_{i}$ do not change. So when DistOrb is used without EqTide, we calculate all of the $f_{i}$ values at the start of the simulation and store them in an array. This accelerates the model by over a factor of 100 compared to recalculating $f_{i}$ every time step.

However, the tidal forces in EqTide do change the semi-major axes, and so when the two models are coupled, we must recalculate $f_{i}$. Rather than recalculate every time step, we additionally store the derivatives, $df_{i}/d\alpha$, and the value of $\alpha$ (the semi-major axis ratio for each pair of planets) at which $f_{i}$ was calculated. Every time step, then, we can calculate the change in $\alpha$, and recalculate $f_{i}$ only when

[TABLE]

where $\eta$ is a user-defined tolerance factor. The smaller $\eta$ is, the more accurate the simulation will be, at the expense of computation time.

14.4.1 Apsidal Locking

Multi-planet systems in which one or more planets experience strong tidal damping of $e$ can reach a so-called fixed point state in which the angular momentum exchange between the planets ends and the longitudes of periastron circulate with identical frequency (Wu & Goldreich, 2002; Rodríguez et al., 2011). In Fig. 22, we use VPLanet’s DistOrb and EqTide modules to reproduce the evolution of $e$ and $\varpi$ for the CoRoT-7 system examined by Rodríguez et al. (2011). This figure should be compared to Figs. 2–3 from Rodríguez et al. (2011). Note that our model is purely secular whereas Rodríguez et al. (2011) directly integrated the equations of motion, including resonant effects, so our results may slightly differ. We adopted the initial conditions from Rodríguez et al. (2011) and list them in Table 3. Our VPLanet simulations qualitatively reproduce the evolution examined in the original work by Rodríguez et al. (2011). Within the first few Myrs, the inner planet, CoRoT-7 b, experiences large eccentricity oscillations that drive its inward tidal migration. As these eccentricity oscillations damp towards 0, both planets enter the fixed point state where the differences between their longitudes of periastron damp to 0, after which their $\varpi$’s circulate together with the same frequency. The tidal eccentricity damping for both planets occurs rapidly within the first 10 Myr, similar to the damping time of about 7 Myr found in Rodríguez et al. (2011).

14.4.2 Cassini States

In this section we consider the damping of a planet’s obliquity into a Cassini state (Colombo & Shapiro, 1966; Ward & Hamilton, 2004; Winn & Holman, 2005; Deitrick et al., 2018a). In a two-planet system, such a configuration occurs when a planet’s orbital and rotational angular momentum vectors remain coplanar with the system’s total angular momentum vector. Such an alignment could occur by chance, but it is most likely for a world to reach a Cassini state if its obliquity experiences both damping and excitation, i.e., a damped driven configuration. In that case, the obliquity reaches a non-zero equilibrium value, such as the $\sim 6^{\circ}$ obliquity of the moon, as first noted by Giovanni Cassini himself.

The physics and mathematics of Cassini states have been discussed at length in the literature, but, briefly, any given three (or more) body system will have Cassini states available to its members. Each member has up to 4 possible Cassini states available, but only up to 2 can be stable. One or more separatrices exist in the phase space, and if damping is present, the stable Cassini states represent attractors.

A Cassini state can be quantitatively identified using the following relations from Ward & Hamilton (2004), see also (Deitrick et al., 2018a):

[TABLE]

where $\mathbf{k}$, $\mathbf{n}$, and $\mathbf{s}$ are the vectors associated with the perpendicular to the appropriate reference plane (the invariable plane or Laplace plane, for example), the angular momentum of the body’s orbit, and the angular momentum of the body’s rotation. Alternatively, the complimentary relation can be used:

[TABLE]

A Cassini state occurs when $\cos\Psi$ and/or $\sin\Psi$ oscillate about 1, 0, or -1, with the equilibrium value depending on the particular Cassini state.

To demonstrate evolution into a Cassini state, we construct a simulation based on Figure 2 of Winn & Holman (2005). The planetary system parameters are listed in Table LABEL:tab:cassini, and with a stellar mass of $1{M_{\odot}}$. We use the EqTide, DistOrb, and DistRot modules to perform this experiment. The evolution of the system is shown in Fig. 23 and the obliquity settles into an equilibrium value of $\sim 59^{\circ}$. In this case, we include first order GR corrections, see $\S$ C.

Fig. 24 shows the phase space of this configuration, as well as the evolution of this test system. Cassini states 1 and 2 are stable, but state 4 is a saddle point and hence is unstable. The gray curves show lines of constant Hamiltonian (Equation 5 in Winn & Holman 2005) and the black curve shows the separatrix between states 1 and 2. In this case, the system is attracted to Cassini state 2 after $\sim 500$ kyr.

The results here for Cassini states should be regarded as approximate. The addition of triaxiality, due to rigidly-supported features or tidal forces, modifies the precession constant and becomes particularly important in spin-orbit resonance (Goldreich & Soter, 1966; Peale, 1969; Hubbard & Anderson, 1978; Jankowski et al., 1989; Bills, 2005; Baland et al., 2016). This changes the precise locations and evolution of the states. VPLanet does not yet include the scaling for triaxiality due to tides or the resonant terms in the obliquity evolution. Importantly, several studies have shown that state 2 may be destabilized by these effects or other physical processes (Gladman et al., 1996; Levrard et al., 2007; Fabrycky et al., 2007).

14.5 Evolution of Venus

14.5.1 Interior

The interior of Venus is modeled using ThermInt + RadHeat, where Venus is assumed to have the same core mass fraction ($0.32$), radiogenic budget in the mantle ($20$ TW today) and core ($3$ TW today), and mantle and core melting curves. The difference between the Venus and Earth models is Venus’ mantle is assumed to be in a stagnant lid which reduces the mantle heat flow. The main constraint on the interior is to ensure no dynamo today. To achieve this the bulk mantle viscosity and activation energy are assumed to be $2\times 10^{9}$ m2 s*-1* and $3.5\times 10^{5}$ J mol*-1*. We note that this model is similar to that in Driscoll & Bercovici (2014), with the main difference here being that the viscosity depends on melt fraction.

Figures 25 and 26 show the thermal evolution of the mantle and core. Due to the isolation of the stagnant lid the core heats up over time and the mantle cools very little, which causes thermal convection and dynamo action to cease in the Venusian core around $0.5$ Gyr. Decreasing the mantle viscosity would cause the dynamo to die later or not at all. These results are very similar to those of Driscoll & Bercovici (2014), Fig. 6, but are not an exact match as our model has been updated slightly.

14.5.2 Atmospheric Loss

Venus is widely believed to have once had a substantial surface water inventory (comparable to that of Earth) that was subsequently lost to both thermal and non-thermal escape processes (Watson et al., 1981; Kasting et al., 1984; Kasting, 1988; Chassefière, 1996a, b; Gillmann et al., 2009). However, the total initial amount of water and the rate at which hydrogen escaped are extremely uncertain, as measurements of D/H fractionation (e.g., Donahue et al., 1982) only place a lower limit on the total amount lost.

Figure 27 shows the evolution of Venus’ water content due to hydrodynamic escape in the first ${\sim}200$ Myr following its formation, assuming the XUV evolution law of Ribas et al. (2005) and the escape efficiency model of Bolmont et al. (2012). The results in the left panel of Figure 27 are consistent with estimates that Venus may have lost on the order of one to a few terrestrial oceans of water in the first several hundred Myr. The dashed vertical line at $t=280$ Myr is the timescale for the loss of 1 TO predicted by Watson et al. (1981). As can be seen from the figure, our model predicts that nearly 5 TO can be lost in that amount of time. This discrepancy is due to two reasons. First, Watson et al. (1981) assumed a constant value of $\epsilon_{\mathrm{XUV}}\mathcal{F}_{\mathrm{XUV}}\approx 2\,\mathrm{erg/cm^{2}/s}$ for Venus, which they computed from estimates of the current XUV flux at Earth and an efficiency $\epsilon_{\mathrm{XUV}}=0.15$. Since then, studies have shown that the XUV flux from the Sun was about two orders of magnitude higher during the first 100 Myr (Ribas et al., 2005), resulting in a much shorter timescale for ocean loss. Second, Watson et al. (1981) did not account for the hydrodynamic drag of oxygen, which strongly damps the net rate of ocean loss during the first several tens of Myr (Luger & Barnes, 2015). Together, these effects lead to a timescale for the loss of 1 TO that is approximately 3 times shorter: about 100 Myr.

For reference, the right panel in the Figure shows the amount of photolytically produced oxygen that is retained by the planet, ranging from a few to several hundred bars. The vast majority of this oxygen would have gone into oxidation of the surface.

14.6 Water Loss During the Pre-Main Sequence

We recover the result that planets near the inner edge of the HZ of low mass M dwarfs can lose one to several oceans of water and produce hundreds to a few thousand bars of atmospheric O2. In Figure 28 we reproduce Figure 7 in Luger & Barnes (2015), showing the amount of water lost and the amount of atmospheric oxygen that builds up for a water-rich Earth-mass planet orbiting an M dwarf. Our water loss estimates are somewhat lower than those in Luger & Barnes (2015) (by up to a factor of 2), primarily because of the lower escape efficiency predicted by the model of Bolmont et al. (2017), which we use here, see Appendix A. We also account for the increasing mixing ratio of oxygen as water is lost, which acts to slow the escape of hydrogen.

15 Discussion and Conclusions

In the previous sections we described the VPLanet algorithm ($\S$ 2) and how individual modules ($\S\S$ 3–13) and module combinations ($\S$ 14) reproduce various previous results. Moreover, VPLanet has already been used for in novel investigations (Deitrick et al., 2018a, b; Fleming et al., 2018; Lincowski et al., 2018; Fleming et al., 2019), demonstrating that this approach can can provide new insight into planetary system evolution and planetary habitability. Furthermore, these insights can be tested; for example, Fleming et al. (2018) coupled STELLAR and EqTide to derive a mechanism that removes circumbinary planets orbiting tight binaries, an hypothesis that can be falsified by upcoming TESS observations — should it discover such planets, then the model is incorrect.

While the previous sections showed the coupling of many modules, not all module couplings have been tested yet. This situation is partly due to some modules being incompatible, e.g., DistOrb and SpiNBody, but also due to the sheer number of combinations that are possible. With this first release, 20 module combinations have been validated against observations and previous results, see $\S\S$ 3–14. Future research will explore more combinations, but this first version of the code includes a large number of processes that affect planetary system evolution. Future versions will include new physics and additional couplings between modules.

VPLanet can simulate a wide range of planetary systems, but it is still an incomplete model of planetary evolution. The relatively simple modules have important limitations and caveats, which are discussed at length in the previous sections and appendices. Users should consult these sections prior to performing simulations to ensure that they are not pushing the models into unrealistic regions of parameter space. Furthermore, we also urge caution when coupling module combinations not explicitly validated here as their stability and/or accuracy cannot be guaranteed.

The modularity of VPLanet and the spread of modules in this first release provide a framework to build ever more sophisticated models. Future versions could be tailored to particularly interesting systems such as TRAPPIST-1, or particular observations such as planetary spectra. For example, a magma ocean module could be created based on the Schaefer et al. (2016) model for GJ 1132 b, and could be combined with tidal heating and a range of radiogenic heating by including the RadHeat and EqTide modules. Or stellar flaring could be added to provide more realistic simulations of atmospheric mass loss.

The discovery of life beyond the Solar System is challenging, in part because resources are scarce and planets are complicated systems. VPLanet’s flexibility and speed permits parameter sweeps that can help allocate those resources efficiently, be they telescopes or computer time for more sophisticated, i.e., computationally expensive, software packages. While numerous models and codes have been created to simulate planetary evolution, we are aware of none that is as broad and flexible as VPLanet. This paper has described not just its physics modules, but also a novel software design that facilitates interdisciplinary science: the function pointer matrix (see $\S$ 2). Furthermore the open source nature of the code, extensive documentation, and code integrity checks (see Appendix M) ensure transparency and reproducible results. These software engineering practices combined with the rigorous validations described in $\S\S$ 3–14 ensure that VPLanet is a reliable platform for the study of planetary system evolution and planetary habitability.

This work was supported by the NASA Virtual Planetary Laboratory Team which is funded under NASA Astrobiology Institute Cooperative Agreement Number NNA13AA93A, and Grant Number 80NSSC18K0829. Additional support was provided by NASA grants NNX15AN35G, and 13-13NAI7_0024. DPF is supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program - Grant 80NSSC17K0482. This work also benefited from participation in the NASA Nexus for Exoplanet Systems Science (NExSS) research coordination network. We thank an anonymous referee whose comments greatly improved the quality of this manuscript. We are also grateful for stimulating conversations with Brian Jackson, Héctor Martinez-Rodriguez, Terry Hurford, Ludmila Carone, Juliette Becker, John Ahlers, Quadry Chance, and Nathan Kaib.

Appendix A The AtmEsc Module

The escape of a planet’s atmosphere to space is an extremely complex process. The rate at which a gaseous species escapes from a planet strongly depends on factors including, but not limited to, the magnetic properties of the planet and the host star, the space weather the planet is exposed to, the wavelength-dependent irradiation of the planet’s atmosphere, as well as the temperature-pressure profile of the atmosphere and its detailed composition, down to the abundance of trace gases that can act as coolants. Decades of work on solar system bodies have enabled the measurement and modeling of the escape fluxes from the Earth, Mars, and Venus using complex hydrodynamic and kinetic models (Hunten, 1973; Watson et al., 1981; Donahue et al., 1982; Kasting & Pollack, 1983; Hunten et al., 1987; Zahnle et al., 1988; Chassefière & Leblanc, 2004). For extrasolar planets, however, the situation is drastically different. Even for the most well-studied exoplanets, little is known at present about their bulk properties other than their radii, their instellations, and occasionally their masses. Some constraints have been placed on the bulk atmospheric composition of some hot exoplanets via transit transmission spectroscopy, but even in the most favorable cases, little is known other than the presence or absence of a large hydrogen/helium envelope (e.g., Nortmann et al., 2018; Allart et al., 2019) or loose constraints on the presence of simple molecules such as CO2 and H2O (e.g., Line et al., 2014; MacDonald & Madhusudhan, 2019). Stellar activity measurements can yield information about the space weather that some of these exoplanets are exposed to, but the measurement of an exoplanet’s magnetic properties is yet to be made (e.g., Driscoll & Olson, 2011; Lynch et al., 2018). On the observation front, hydrogen escape fluxes have been inferred for only a few large, hot exoplanets from Lyman-alpha absorption measurements (e.g., Odert et al., 2019).

However, while precious little is known about the atmospheric escape process from most (individual) exoplanets, recent studies have leveraged the statistical information from the ensemble of all known exoplanets to infer trends in atmospheric escape as a function of planet size and irradiation (Lopez & Rice, 2016; Owen & Wu, 2017). These studies show that the distribution of radii of hot exoplanets discovered by the Kepler mission are well explained, on average, by a surprisingly simple atmospheric escape model, introduced by Watson et al. (1981) and based on investigations of the solar wind by Parker (1964). In what is commonly referred to as an energy-limited model, the escape from a planetary atmosphere is driven by the supply of energy to the upper atmosphere by stellar extreme ultraviolet (XUV; 1–1000Å) photons, which are absorbed by hydrogen atoms and converted into kinetic energy. In the simplest form of the model, a fixed fraction $\epsilon_{\mathrm{XUV}}$ of the incoming XUV energy goes into driving the escape (Watson et al., 1981; Erkaev et al., 2007; Lammer et al., 2013; Volkov & Johnson, 2013; Johnson et al., 2013). For a hydrogen-dominated atmosphere, the energy-limited particle escape rate $F_{\mathrm{EL}}$ is obtained by equating the energy provided by XUV photons to the energy required to lift the atmosphere out of the gravitational potential well:

[TABLE]

where $\mathcal{F}_{\mathrm{XUV}}$ is the XUV energy flux, $M_{\mathrm{p}}$ is the mass of the planet, $R_{\mathrm{p}}$ is the planet radius, $\epsilon_{\mathrm{XUV}}\approx 0.1$ is the XUV absorption efficiency, and $K_{\mathrm{tide}}$ is a tidal correction term of order unity (Erkaev et al., 2007). The total escape rate is this quantity integrated over the surface area of the planet, whose effective radius to incoming XUV energy is $R_{\mathrm{XUV}}\approx R_{p}$. For terrestrial planets, we compute this quantity as

[TABLE]

(Lehmer & Catling, 2017), where $H$ is the atmospheric scale height, $p_{XUV}$ is the pressure at the effective XUV absorption level, and $p_{s}$ is the pressure at the surface.

In the absence of detailed information about the properties of an exoplanet that can control or modulate the atmospheric escape rate, we implement this simple model for atmospheric escape in VPLanet, with a few modifications to explicitly model the escape rate from potentially habitable terrestrial planets. Our model closely follows that of Luger et al. (2015) and Luger & Barnes (2015). Here we briefly discuss the principal equations and slight modifications to the models presented in those papers.

In VPLanet, we model atmospheric escape from two basic types of atmospheres: hydrogen-dominated atmospheres, such as that of an Earth or super-Earth with a thin primordial hydrogen/helium envelope, and water vapor-dominated atmospheres, such as that of a terrestrial planet in a runaway greenhouse. In the former case, we compute the escape in the energy-limited regime, Equation (A1), and assume that the hydrogen envelope must fully escape before any other volatiles can be lost to space, given the expected large diffusive separation between light H atoms and other atmospheric constituents. If the envelope is not lost by the time the star reaches the main sequence, we shut off the escape process to account for the transition to ballistic escape predicted by Owen & Mohanty (2016). We model the planet’s radius with the evolutionary tracks for super-Earths of Lopez et al. (2012) and Lopez & Fortney (2014). If an exoplanet loses its H/He envelope, we compute its solid radius using the Sotin et al. (2007) mass-radius relation. The XUV flux is computed from stellar evolution tracks (Appendix J) and the XUV absorption efficiency parameter $\epsilon_{\mathrm{XUV}}$ is a tunable constant.

In the case of a terrestrial planet with no hydrogen/helium envelope, we assume atmospheric escape only takes place if the total flux incident on the planet exceeds the runaway greenhouse threshold, computed from the equations in Kopparapu et al. (2013). Typically, the fluxes experienced by planets in or near the habitable zone are not high enough to drive the hydrodynamic escape of the high mean molecular weight bulk atmosphere. Although water vapor can be photolyzed by stellar ultraviolet photons, liberating hydrogen atoms that can go on to escape hydrodynamically, on Earth this process is strongly inhibited by the stratospheric cold trap, which prevents water molecules from reaching the upper atmosphere. However, during a runaway greenhouse, the surface temperature exceeds the temperature of the critical point of water (647 K) and the surface oceans fully evaporate, leading to an upper atmosphere that is dominated by water vapor (e.g., Kasting, 1988). As in Luger & Barnes (2015), we use the energy-limited formalism to compute the loss rate of a planet’s surface water via escape of hydrogen to space, with modifications to allow for the hydrodynamic drag of oxygen by the escaping hydrogen atoms. The total hydrogen particle escape rate is (Luger & Barnes, 2015):

[TABLE]

where

[TABLE]

is the diffusion-limited flux of oxygen atoms through a background atmosphere of hydrogen and

[TABLE]

is the crossover mass, the largest particle mass that can be dragged upward in the flow. In the expressions above, $m_{\mathrm{H}}$ and $m_{\mathrm{O}}$ are the masses of the hydrogen and oxygen atoms, respectively, $k_{\mathrm{boltz}}$ is the Boltzmann constant, $T_{\mathrm{flow}}$ is the temperature of the hydrodynamic flow, set here to 400 K (Hunten et al., 1987; Chassefière, 1996b), $b_{\mathrm{diff}}=4.8\times 10^{17}(T_{\mathrm{flow}}/\mathrm{K})^{0.75}\ \mathrm{cm^{-1}s^{-1}}$ (Zahnle & Kasting, 1986) is the binary diffusion coefficient for the two species, $g$ is the acceleration of gravity, and $X_{\mathrm{O}}$ is the oxygen molar mixing ratio at the base of the flow, equal to $\frac{1}{3}$ when the upper atmosphere is water vapor-dominated. As in Tian (2015) and Schaefer et al. (2016), we account for the increasing mixing ratio of oxygen at the base of the hydrodynamic flow, which slows the escape of hydrogen. Tian (2015) finds that, as oxygen becomes the dominant species in the upper atmosphere, the Hunten et al. (1987) formalism predicts that an oxygen-dominated flow can rapidly lead to the loss of all O2 from planets around M dwarfs. However, hydrodynamic oxygen-dominated escape requires exospheric temperatures $\sim m_{\mathrm{O}}/m_{\mathrm{H}}=16$ times higher than that for a hydrogen-dominated flow, which is probably unrealistic for most planets. Following the prescription of Schaefer et al. (2016), we therefore shut off oxygen escape once its mixing ratio exceeds $X_{\mathrm{O}}=0.6$ (corresponding to an equal number of O2 and H2O molecules at the base of the flow), switching to the diffusion-limited escape rate of hydrogen. Finally, users can either choose a constant value for $\epsilon_{\mathrm{XUV}}$ or model it as a function of the incoming XUV flux as in Bolmont et al. (2017). In Bolmont et al. (2017), the authors modeled atmospheric loss with a set of 1D radiation-hydrodynamic simulations that allowed them to calculate the XUV escape efficiency, $\epsilon_{\mathrm{XUV}}$. In Fig. 29, we show a subset of the range of $\epsilon_{\mathrm{XUV}}$ values calculated in Bolmont et al. (2017) for some values of stellar XUV flux received by the planet. As STELLAR can calculate the star’s changing XUV flux over time, the XUV escape efficiency will change over time if set by the Bolmont et al. (2017) model.

Currently, AtmEsc does not model the absorption of oxygen by surface sinks, although users can run the code in two limiting cases: efficient surface sinks, corresponding to (say) a reducing magma ocean that immediately absorbs any photolytically produced oxygen; and inefficient surface sinks, corresponding to (say) a fully oxidized surface, leading to atmospheric buildup of O2 over time. Upcoming modifications to AtmEsc will couple it to the geochemical evolution of the planet’s mantle in order to more realistically compute the rate of oxygen buildup in a hydrodynamically escaping atmosphere.

As a word of caution, it is important to reiterate that the energy-limited formalism we adopt in AtmEsc is a very approximate description of the escape of an atmosphere to space. The heating of the upper atmosphere that drives hydrodynamic escape is strongly wavelength dependent and varies with both the composition and the temperature structure of the atmosphere, which we do not model. Moreover, line cooling mechanisms such as recombination radiation scale non-linearly with the incident flux. Non-thermal escape processes, such as those controlled by magnetic fields, flares, and/or coronal mass ejections, lead to further departures from the simple one-dimensional energy-limited escape rate. Nevertheless, as we argued above, several studies show that for small planets the escape rate does indeed scale with the stellar XUV flux and inversely with the gravitational potential energy of the gas (e.g., Lopez et al., 2012; Lammer et al., 2013; Owen & Wu, 2013, 2017) and that $\epsilon_{\mathrm{xuv}}\approx 0.1$ is a reasonable median value that predicts the correct escape fluxes within a factor of a few. Since presently we have little information about the atmospheric structure of exoplanets, we choose to employ the energy-limited approximation and fold all of our uncertainty regarding the physics of the escape process into the XUV escape efficiency $\epsilon_{\mathrm{xuv}}$. Future versions of VPLanet can include diverse models, such as radiation-recombination limited escape (e.g., Murray-Clay et al., 2009), to more accurately track atmospheric evolution.

Appendix B The BINARY Module

Here we describe the analytic theory for circumbinary orbits of test particles derived by Leung & Lee (2013). We adopt a cylindrical coordinate system centered on the binary center of mass, the system barycenter in this case, and consider the test particles to be massless circumbinary planets (CBPs). Assuming that the binary orbit lies in the $x-y$ plane, and that the orbit of the CBP is nearly coplanar, Leung & Lee (2013) approximate the gravitational potential felt by the CBP due to the binary at the position $(R,\phi,z)$ as

[TABLE]

for integer $k$, binary mean anomaly $M_{b}$, binary orbital eccentricity $e_{AB}$, binary longitude of the periapse $\varpi_{B}$, and R is the radial distance from the CBP to the binary center of mass. This expression, correct to first order in $e_{AB}$ as indicated by the ellipses, contains the two gravitational potential components from the stars that are in general not axisymmetric. The two components, $\Phi_{jk0}(R)$ and $\Phi_{jk1}(R)$, are given by the following expressions

[TABLE]

and

[TABLE]

where $b^{k}_{(j+1)/2}$ is a Laplace coefficient, $D=\partial/\partial\alpha$, $m_{A}$ and $m_{B}$ are the masses of the primary and secondary star, respectively, $G$ is the Universal Gravitational constant, and $\delta_{k0}$ is the Kroeneker delta function. For a CBP located at cylindrical position $R$, $\alpha_{A}$ and $\alpha_{B}$ are the normalized semi-major axis of the CBP relative to each star, given by

[TABLE]

where $a_{AB}$ is the binary orbital semi-major axis and the index $i$ is $A$ for the primary and $B$ for the secondary star, respectively.

Given the approximation for the binary gravitational potential in Eq. (B1), $\Phi$, the equations that govern the motion of the CBP in cylindrical coordinates are given by

[TABLE]

Following Lee & Peale (2006), Leung & Lee (2013) approximate the orbit of the CBP as small epicyclic deviations from the circular motion of the guiding center via

[TABLE]

where $R_{0}$ is the cylindrical radius of the CBP guiding center in the plane of the binary orbit, [math] subscripts denote the position of the guiding center, $R_{1}$ is the small, time-dependent epicyclic radial deviation from the guiding center radius.

The Keplerian mean motion at the radius of the CBP’s guiding center is

[TABLE]

and the phase angle of the circular motion of the CBP’s guiding center is

[TABLE]

for time $t$ and arbitrary constant phase offset, $\psi_{0}$.

The CBP mean motion is given by

[TABLE]

where this expression is evaluated at the radius of the guiding center, $R=R_{0}$.

Given these assumptions and definitions, Leung & Lee (2013) solve the equations of motion given in Eq. (B5) to obtain the radial position, $R$, of the CBP relative to the binary center of mass, the position of the CBP above or below the orbital plane of the binary, $z$, and the phase angle of the CBP, $\phi$, over the course of its orbit as an analytic function of time. The CBP’s radial position, $R$, is given by the following expression

[TABLE]

where $e_{free}$ is the CBP’s free eccentricity, a free parameter of the model, and $\psi$ is an arbitrary phase offset. Although the summations in Eq. (B10) extend to $\infty$, we follow Leung & Lee (2013) and truncate the summation at $k=3$ for computational speed with minimal loss of accuracy.

The variables $C_{0}$, $C_{k}^{0}$, and $C_{k}^{\pm}$ represent the approximate radial amplitudes for the forced oscillations due to the non-axisymmetric components of the binary gravitational potential and are given by the following

[TABLE]

[TABLE]

and

[TABLE]

In these expressions, $k_{0}$ is the epicyclic frequency and is given by

[TABLE]

for CBP mean motion, $n$ given by Eq. (B9), and is evaluated at CBP radial position $R=R_{0}$.

The CBP’s cylindrical phase angle relative to the binary’s center of mass is

[TABLE]

where $\varphi$ is an arbitrary phase offset and the variables $D_{0}$, $D_{k}^{0}$, and $D_{k}^{\pm}$ are given by

[TABLE]

[TABLE]

and

[TABLE]

The cylindrical position of the CBP above or below the plane of the binary, $z$, is decoupled from the epicyclic motion of the CBP orbital radius and phase angle and is simply

[TABLE]

where $\zeta$ is an arbitrary phase offset and $i_{\text{free}}$ is the free inclination, a free parameter of the model. The vertical frequency $\nu_{0}$ is given by

[TABLE]

Appendix C The DistOrb Module

Our model for the orbital evolution, called DistOrb (for “Disturbing function Orbits”), uses the literal, 4th order disturbing function developed in Murray & Dermott (1999) and Ellis & Murray (2000). We use only the secular terms, meaning that the rapidly varying terms that depend on the mean longitudes of the planets are ignored on the assumption that these terms will average to zero over long timescales. This assumption is valid as long as no planets are in the proximity of mean-motion resonances.

There are two solution methods in DistOrb: the first is a direct Runge-Kutta integration of the fourth-order equations of motion; the second is the Laplace-Lagrange eigenvalue solution, which reduces the accuracy in the disturbing function to 2nd order, but returns a solution that is explicit in time, and thus provides a solution in much less computation time and is stable by assumption. The 4th order solution tends to produce results that match better with N-body; however, the Laplace-Lagrange solution can be a powerful predictive tool because it provides secular frequencies directly. One example of the use of such frequencies is in the prediction of Cassini states, as shown in $\S$ 14.4.2 (see also Ward & Hamilton, 2004; Brasser et al., 2014; Deitrick et al., 2018a).

The equations of motion are Lagrange’s equations (see Murray & Dermott, 1999). In the secular approximation the equations for semi-major axis and mean longitude, and any disturbing function derivative with respect to these variables, are ignored. Additionally, to avoid singularities in the equations for the longitudes of pericenter and ascending node, which occur at zero eccentricity and inclination, respectively, we rewrite Lagrange’s equations and the disturbing function in terms of the variables (a form of Poincaré coordinates):

[TABLE]

where $e$ is the orbital eccentricity, $i$ is the inclination, $\Omega$ is the longitude of ascending node, and $\varpi=\omega+\Omega$ is the longitude of periastron (see Figure 30).

Lagrange’s equations for secular theory are then:

[TABLE]

where $\mathcal{R}$ is the disturbing function, and $a$, $n$, and $e$ are the semi-major axis, mean motion, and eccentricity, respectively. See Berger & Loutre (1991) for the complete set of Lagrange’s equations in $h$, $k$, $p$, and $q$, including mean-motion (i.e., resonant) effects.

General relativity is known to affect the apsidal precession (associated with eccentricity) of planetary orbits, so we include a correction to Equations (C5) and (C6). Following Laskar (1986), the apsidal corrections are:

[TABLE]

where

[TABLE]

and $c$ is the speed of light.

For systems well away from mean-motion resonances, DistOrb does a reasonable job at modeling the orbital evolution, up to eccentricities of $\sim 0.4$ or mutual inclinations of $\sim 35^{\circ}$, largely capturing the frequencies and amplitudes of oscillations to within $5-10\%$. Note, however, that the accuracy is difficult to quantify in general because of the highly non-linear nature of the orbital problem. See Deitrick et al. (2018a) for further comparisons. However, as noted in Murray & Dermott (1999), one of the series expansions used in deriving the disturbing function diverges at $e=0.6627434$. If the eccentricity of a planet ever exceeds this value, DistOrb immediately halts. In this case, the user will need to switch to an N-body integrator such as SpiNBody, or another N-body code, which can still be coupled to DistRot and POISE if desired.

Because planetary systems have a large number of parameters (masses, number of planets, eccentricities, inclinations, etc.), the full limitations of DistOrb are difficult to map. Thus we advise users to always compare a small selection of cases with an N-body model, to understand DistOrb’s applicability to the desired planetary system.

Finally we present, for the sake of completeness, the disturbing function as used in DistOrb, in the variables $h,k,p,$ and $q$ in Table 5. These were originally derived by Ellis & Murray (2000); we have simply applied coordinate transformations and calculated derivatives with respect to the new coordinates. This disturbing function, in its original form, can also be seen in Murray & Dermott (1999). We will not restate the semi-major axis functions, $f_{1},f_{2},f_{3}$ and so on, in this paper, as they are taken directly from Table B.3 of Murray & Dermott (1999).

The secular disturbing function, for any pair of planets, is:

[TABLE]

for the inner body, and

[TABLE]

for the outer body. Here, $a^{\prime}$ is the semi-major axis of the exterior planet and the mass factors are $\mu=\kappa^{2}m$ and $\mu^{\prime}=\kappa^{2}m^{\prime}$, where $m$ is the mass of the interior planet, $m^{\prime}$ is the mass of the exterior planet, and $\kappa$ is Gauss’ gravitational constant. Finally,

[TABLE]

where the terms D0.1, D0.2, and so on, are given in Table 5.

Appendix D The DistRot Module

The rotational axis model, DistRot (for “Disturbing function Rotation”), is derived from the Hamiltonian for rigid body motion introduced by Kinoshita (1975, 1977) and has been used extensively (e.g., Laskar, 1986; Laskar et al., 1993; Armstrong et al., 2004, 2014). In the absence of large satellites (such as the Moon), the equations of motion for a rigid planet are:

[TABLE]

where $\psi$ is the “precession angle” (see Figure 30 and the following paragraph), $\varepsilon$ is the obliquity, and,

[TABLE]

Note the sign error in Equation (8) of Armstrong et al. (2014), corrected in our Equation (D5). Our Equation (D3) does not contain the lunar constants present in Equation (24) of Laskar (1986). Here, $\varepsilon$ represents the obliquity, $p$ and $q$ are the inclination variables from Eqs. (C3) and (C4), $\dot{p}$ and $\dot{q}$ are their time derivatives (Equations [C7] and [C8]), $\kappa$ is the Gaussian gravitational constant, $M_{\star}$ is the mass of the host star in solar units, $\nu$ is the rotation frequency of the planet in rad day*-1*, $CM^{-1}r^{-2}$ is the specific polar moment of inertia of the planet, and $J_{2}$ is the gravitational quadrupole of the (oblate) planet. The final term in Equation (D1), $p_{g}$, accounts for precession due to general relativity and is equal to $\delta_{R}/2$ (Barker & O’Connell, 1970), where $\delta_{R}$ is the apsidal precession rate and is given by Equation (C11).

The symbol $\psi$ refers to the precession angle, defined as $\psi=\Lambda-\Omega$, where $\Lambda$ is the angle between the vernal point $\vernal$, the position of the Sun/host star at the planet’s northern spring equinox, and the location of the ascending node, $\Omega$, measured from some reference direction $\vernal_{0}$ (often taken to be the direction of the vernal point at some reference date, hence the use of the symbol $\vernal$333The vernal point occurred in the constellation Aries during Ptolemy’s time; thus it is also called the “first point of Aries” and the “ram’s horn” symbol is used.), see Figure 30. The convention of defining $\vernal$ as the location of the Sun at northern spring equinox is sensible for the solar system since we observe from Earth’s surface; however, it is a confusing definition to use for exoplanets, for which the direction of the rotation axis is unknown. We adhere to the convention for the sake of consistency with prior literature. In the coming decades, it may become possible to determine the obliquity and orientation of an exoplanet’s spin axis; in that event, care should be taken in determining the initial $\psi$ for obliquity modeling. One can equivalently refer to $\vernal$ as the position of the planet at its northern spring equinox $\pm 180^{\circ}$. The relevant quantity for determining the instellation, however, is the angle between periastron and the spring equinox, $\Delta^{*}=\omega+\Lambda+180^{\circ}=\varpi+\psi+180^{\circ}$ (See Section G).

An additional complication for obliquity evolution is, of course, the presence of a large moon. We do not include the component of the Kinoshita (1975) model that accounts for the lunar torque because the coefficients used are specific to the Earth-Moon-Sun three body problem and were calculated from the Moon’s orbital evolution (and are therefore not easily generalized). However, in DistRot we can approximate the effect of the Moon by forcing the precession rate, Eq. (D3), to be equal to the observed terrestrial value. The effect of this is shown in Figure 8.

Equation (D1) for the precession angle contains a singularity at $\varepsilon=0$. To avoid numerical instability, we instead recast Equations (D1) and (D2) in terms of the rectangular coordinates:

[TABLE]

The third coordinate, $\chi$, is necessary to preserve sign information when the obliquity crosses 90*∘* (see Laskar et al., 1993). The equations of motion for these variables are then:

[TABLE]

The value of $J_{2}$ is a function of the planet’s rotation rate and density structure. In hydrostatic equilibrium the planet’s gravitational quadrupole moment scales with the rotation rate squared (Cook, 1980; Hubbard, 1984). This result simply comes from force balance: equipotential surfaces will be those for which the combined gravitional and centrifugal force is constant. Then in hydrostatic equilibrium, the physical surface will be an equipotential surface. The rapid rotators of the Solar System (Earth, Mars, and the giants) appear to be close to hydrostatic equilibrium, in that the $J_{2}$ is dominated by the rotational bulge, rather than topographic or other non-uniform features. For the slow rotators, such as Venus, the contribution to $J_{2}$ from rotation is small compared with the contribution from topographic features (Yoder, 1995), indicating a departure from hydrostatic balance. In DistRot, we give the option to assume that planets are in hydrostatic equilibrium, and scale $J_{2}$ according to:

[TABLE]

where $\nu$ is the rotation rate, $r$ is the mean planetary radius, and $M$ is its mass. Earth’s measured values of $J_{2\oplus}=1.08265\times 10^{-3}$, $\nu_{\oplus}=7.292115\times 10^{-5}$ radians s*-1*, $R_{\oplus}=6.3781\times 10^{6}$ m, and $M_{\oplus}=5.972186\times 10^{24}$ kg are used for reference ($J_{2}$ from Cook 1980, $R_{\oplus}$ and $M_{\oplus}$ from Prša et al. 2016). This approach is identical to the method used by Brasser et al. (2014). Like that study, we assume that $J_{2}$ cannot go below the measured value of Venus from Yoder (1995), which would ordinarily occur at rotation periods $\gtrsim 13$ days. The assumption behind this scaling is then is that there is a limit to hydrostatic equilibrium for a partially rigid body. Alternatively, the dynamical ellipticity, $E_{D}=J_{2}M^{-1}r^{-2}$, and the polar moment of inertia (see below) can be independently set, allowing for departures from hydrostatic equilibrium.

Note that the current version of VPLanet does not yet take into account the additional contribution to $J_{2}$ from tidal distortion (see Hubbard & Anderson, 1978).

The specific polar moment of inertia of a planet, $CM^{-1}r^{-2}$, is always between 0.2 and 0.4. The $CM^{-1}r^{-2}$ value of the Earth, for example, is $0.33$, (Cook, 1980).

Finally, we note that DistRot can also ingest output from previous orbital simulations, such as from N-body models, and compute the rotational axis through finite differencing. This functionality thus allows for coupling to external orbital models, if desired.

Appendix E The EqTide Module

The tidal model we use is commonly called the equilibrium tide model and was first conceived by George Darwin, grandson of Charles (Darwin, 1880). This model assumes the gravitational potential of the tide raiser on an unperturbed spherical surface can be expressed as the sum of Legendre polynomials (i.e., surface waves) and that the elongated equilibrium shape of the perturbed body is slightly misaligned with respect to the line that connects the two centers of mass. This misalignment is due to dissipative processes within the deformed body and leads to a secular evolution of the orbit, as well as the spin angular momenta of the two bodies. Furthermore, the bodies are assumed to respond to the time-varying tidal potential as though they are damped, driven harmonic oscillators. As described below, this approach leads to a set of six coupled, non-linear differential equations, but note that the model is linear in the sense that there is no coupling between the surface waves which sum to the equilibrium shape. A substantial body of research is devoted to tidal theory (e.g., Darwin, 1880; Goldreich & Soter, 1966; Hut, 1981; Ferraz-Mello et al., 2008; Wisdom, 2008; Efroimsky & Williams, 2009; Leconte et al., 2010), and the reader is referred to these studies for a more complete description of the derivations and nuances of equilibrium tide theory.

E.1 Constant-Phase-Lag Model

In the CPL model of tidal evolution, the angle between the line connecting the centers of mass and the tidal bulge is constant. This approach is commonly utilized in planetary studies (e.g., Goldreich & Soter, 1966; Greenberg, 2009) and the evolution is described by the following equations:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The quantity $Z^{\prime}_{i}$ is

[TABLE]

where $k_{2,i}$ are the Love numbers of order 2, and $Q_{i}$ are the tidal quality factors. The parameter $\xi_{i}$ is

[TABLE]

where $i$ and $j$ refer to the two bodies, and $r_{g}$ is the “radius of gyration,” i.e., the moment of inertia is $M(r_{g}R)^{2}$. The signs of the phase lags are

[TABLE]

with sgn($x$) the sign of any physical quantity $x$, i.e., sgn($x)~{}=~{}+1,-1$ or 0.

The tidal heating of the $i$th body is due to the transformation of rotational and/or orbital energy into frictional heating. The heating from the change in the orbital energy is

[TABLE]

and that from the change in rotational energy is

[TABLE]

The total heat in the $i$th body is therefore

[TABLE]

The rate of evolution and amount of heating are set by three free parameters: $Q,k_{2}$, and $r_{g}$. $\dot{E}_{\mathrm{tide},i}^{\mathrm{CPL}}$ can be plugged into Eq. (K1) for worlds that are experiencing tidal heating.

Goldreich (1966) suggested that the equilibrium rotation period for both bodies is

[TABLE]

Murray & Dermott (1999) present a derivation of this expression, which assumes the rotation rate may take a continuum of values. However, the CPL model described above only permits 4 “tidal waves”, and hence does not contain the resolution to permit this continuum. Instead the rotation rate is synchronous up to $e=\sqrt{1/19}~{}\approx~{}0.23$ and then jumps instantaneously to $\omega/n=1.5$, i.e., a 3:2 spin-orbit frequency ratio. The next phase jump occurs at $e=\sqrt{2/19}\approx 0.32$ but is not present in the 2nd order CPL model. Therefore the evolution at larger $e$ predicted by the CPL model may not be qualitatively correct. We urge caution when interpreting CPL results above $e=0.32$.

E.2 The Constant-Time-Lag Model

The constant-time-lag (CTL) model assumes that the time interval between the passage of the perturber and the tidal bulge is constant. This assumption allows the tidal response to be continuous over a wide range of frequencies, unlike the CPL model. But, if the phase lag is a function of the forcing frequency, then the system is no longer analogous to a damped, driven harmonic oscillator. Therefore, this model should only be used over a narrow range of frequencies, see Greenberg (2009). However, this model’s use is widespread, especially at high $e$, so we use it to evaluate tidal effects as well. Compared to CPL, this model predicts larger tidal heating and evolution rates at high $e$ due to the inclusion of higher order terms. Therefore, the CPL and CTL models probably bracket the actual evolution.

The evolution is described by the following equations:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

The tidal heating of the $i$th body is therefore

[TABLE]

$\dot{E}_{\mathrm{tide},i}^{\mathrm{CTL}}$ can be plugged into Eq. (K1) for worlds that are experiencing tidal heating.

It can also be shown that the equilibrium rotation period for both bodies is

[TABLE]

which for low $e$ and $\psi=0$ reduces to

[TABLE]

There is no general conversion between $Q_{\mathrm{p}}$ and $\tau_{\mathrm{p}}$. Only if $e~{}=~{}0$ and $\psi_{\mathrm{p}}~{}=~{}0$, when merely a single tidal lag angle $\varepsilon_{\mathrm{p}}$ exists, then

[TABLE]

as long as $n-\omega_{\mathrm{p}}$ remains unchanged. Hence, a dissipation value for an Earth-like planet of $Q_{\mathrm{p}}~{}=~{}12$ (Williams et al., 1978) is not necessarily equivalent to a tidal time lag of 638 s (Lambeck, 1977), so the results for the tidal evolution will intrinsically differ among the CPL and the CTL model.

E.3 The Orbit-Only Model

The calculations of $\S$ 14.3 use the orbital evolution model of Driscoll & Barnes (2015), which only considered the orbital effects and ignored dissipation in the star. We refer to this case as the “orbit-only model.” In this case the planet’s semi-major axis $a$ and eccentricity $e$ evolution rates are (Goldreich, 1966; Jackson et al., 2009; Ferraz-Mello et al., 2008)

[TABLE]

and

[TABLE]

The mean motion can be replaced with $n^{2}=GM_{*}/a^{3}$, and after rearrangement we obtain

[TABLE]

The differential equations for thermal evolution (K3, K4) and orbital evolution (E22, E21) are solved simultaneously to compute coupled thermal-orbital evolutions. The tidally heated Earth model example uses the approximation $Im(k_{2})=k_{2}/\mathcal{Q}$, where $\mathcal{Q}=\eta\omega/\mu$, in Eqs.  (E21–E22).

Appendix F The GalHabit Module

The module GalHabit (“Galactic Habitability”) is designed to account for the effects of galactic scale processes on binary stars. The main processes are the galactic tide, stellar migration, and stellar encounters, which, as shown by Kaib et al. (2013), can lead to the destabilization of planetary systems.

F.1 Galactic tides

The galactic tide is a differential force on gravitationally-bound, widely separated objects, such as wide binary star systems or Oort cloud comets orbiting the sun. This force arises because of the variation in the gravitational force with distance from the midplane and distance from the galactic center. At the sun’s galactocentric distance ($\sim$ 8 kpc), the galactic gravitational potential is nearly axisymmetric (disk-like), where the density is highest in the mid-plane and gradually decreases with some scale height in the $Z-$direction, i.e., perpendicular to the plane of the disk (Heisler & Tremaine, 1986). However, the galactic bulge and halo do provide a spheroidal component to the galactic potential that becomes more important close to the galactic center and far from the mid-plane (Kordopatis et al., 2015). GalHabit does not at present account for such departures from axisymmetry, and so should be used with caution when modeling systems in such regions.

The effects of the axisymmetric galactic tide on a system may be described by the time-averaged (or secular) Hamiltonian (Heisler & Tremaine, 1986):

[TABLE]

where $\mu=G(M_{\text{c}}+M)$, $L=\sqrt{\mu a}$, $J=\sqrt{\mu a(1-e^{2})}$, $J_{z}=J\cos{i}$, $\rho_{0}$ is the local galactic density, $M_{c}$ is the central mass, $M$ the mass of the orbiter, $a$ is its semi-major axis, $e$ is its eccentricity, $i$ is its inclination of the orbit with respect to the galactic mid-plane, and $\omega$ is its argument of periastron. The variables $L$, $J$, and $J_{z}$ are thus canonical momenta associated with the energy of the orbit (in momentum units), the total angular momentum, and the $Z-$component of the angular momentum, respectively. From Hamilton’s equations we know that, since the canonical angles associated with $L$ and $J_{z}$ do not appear explicitly in $H_{\text{av}}$, $L$ and $J_{z}$ are perfectly conserved. Since $\omega$, which is associated with $J$ does appear, however, $J$ is free to vary with time. The eccentricity and inclination are thus coupled through the definitions of $J$ and $J_{z}$, and as one changes under the influence of the tide, the other must compensate, resulting in “Lidov-Kozai-like” cycling of $e$ and $i$.

The Hamiltonian above results in four equations of motion (Hamilton’s equations):

[TABLE]

where $\Omega$ is the longitude of ascending node and $l$ is the mean anomaly. Since $L$ is constant and we are modeling the orbit-averaged (secular) evolution of the system, we can disregard Equation (F5) — this choice ultimately affects none of the other variables.

These equations are modeled in this form and the other orbital elements ($e$ and $i$) can be calculated from $J$ and $J_{z}$. Since in future versions of VPLanet we intend to include the additional dynamics of triple star systems, we desire a form of these equations that is more easily generalized to different coordinate systems. Cartesian coordinate systems are ideal as rotational and translational transformations are easily applied. In our case, we utilize the Cartesian components of the angular momentum vector (per unit mass), $\vec{J}=(J_{x},J_{y},J_{z})$, and the eccentricity vector, $\vec{e}=(e_{x},e_{y},e_{z})$. They are defined as:

[TABLE]

The derivatives are calculated via chain-rule from Hamilton’s equations above and written in terms of orbital elements (when convenient), resulting in

[TABLE]

and

[TABLE]

where, for the $x-$components,

[TABLE]

for the $y-$components,

[TABLE]

and for the $z-$components,

[TABLE]

Note that Equations (F12), (F13), (F18), (F21), and (F24) contain singularities at either $e=0$ or $i=0$; however, the offending terms cancel when the algebra in Equations (F12 - F17) is fully carried out. We perform this algebra as necessary to remove the singularities, which are likely to trigger numerical instabilities or unphysical behavior. VPLanet contains the resulting equations.

Our integration variables for the galactic tide have thus grown from three ($J$, $\Omega$, and $\omega$) to five ($J_{x}$, $J_{y}$, $e_{x}$, $e_{y}$, and $e_{z}$); this is the trade-off of using this more flexible coordinate system. In practice, the tidal model is still very fast (simulations of billions of years, with tides alone, take seconds on a CPU), and thus this is an acceptable exchange for future flexibility.

For output, the orbital elements need to be calculated from the Cartesian vectors. The eccentricity is simply $e=|\vec{e}|$, while $i$ and $\Omega$ can be calculated from $J=|\vec{J}|$, $J_{x}$, $J_{y}$, and $J_{z}$ (Equations [F6]— [F8]). For $\omega$, Mardling & Lin (2002) give

[TABLE]

where $\vec{n}=\hat{Z}\times\vec{J}$ and $\hat{Z}$ is the unit vector perpendicular to the galactic plane. Note that $\vec{n}$ contains the full vector $\vec{J}$, not its unit vector $\hat{J}$, and that in the equations above $\hat{n}=\vec{n}/|\vec{n}|$. The semi-major axis, $a$, is constant under the galactic tide alone; however, additional forces (see next Section) lead to changes in $a$. To recalculate $a$, we can first calculate $e$ and then utilize the definition of $J$.

The local density around the host system is calculated using the mass model of Kordopatis et al. (2015), which accounts for disk stars, gas, dark matter, and stars in the galactic bulge. Technically, the latter two components have spheroidal distributions and thus do not fit the axisymmetric assumption of the tidal model. The time-averaging of the galactic Hamiltonian to produce Eq. (F1) requires axisymmetry (Heisler & Tremaine, 1986). Generally, within the disk of the Milky Way, these components are small and axisymmetry is a decent approximation. However, close to the galactic center or at high $Z$, this approximation may not hold. In that case, the unaveraged Hamiltonian (or a reformulation of $H_{\text{av}}$) for the galactic tide may be needed.

F.2 Encounters with passing stars

The second force supplied by the galactic environment is the gravitational perturbation of wide orbits by passing stars. Close encounters between the solar system and nearby stars are thought to affect the Oort cloud and thus to contribute to the appearance of long-period comets (Duncan et al., 1987; Heisler et al., 1987; Kenyon & Bromley, 2004; Rickman et al., 2008; Collins & Sari, 2010). Similarly, such close encounters will affect the orbits of widely separated binary stars (Kaib et al., 2013), potentially affecting orbiting planets. We use the techniques developed by Heisler et al. (1987) and Rickman et al. (2008) to account for these perturbations. The process is described as follows.

At the start of a simulation, we calculate the total stellar encounter rate and the time until the next encounter. The encounter rate is based on the local stellar properties determined by García-Sánchez et al. (2001), see Table 6. Given an “encounter radius” (i.e., the distance at which the system and a passing star are considered to encounter each other), $R_{\text{enc}}$, the rate of encounters is given by

[TABLE]

where $n_{\star i}$ is the number density of each stellar type and $v_{i}=\sqrt{v_{h,i}^{2}+\sigma_{\star i}^{2}}$ (García-Sánchez et al., 2001). Note that $R_{\text{enc}}$ is a fixed distance at which encounters are considered to begin (the size of the encounter box, say), not the impact parameter (i.e. the distance of closest approach). The encounter velocity, $v_{i}$, is a function of both the stellar type’s velocity dispersion, $\sigma_{\star i}$, and the host system’s peculiar (or “apex”) velocity relative to each stellar type, $v_{h,i}$. Using the numbers in García-Sánchez et al. (2001) and $R_{\text{enc}}=1$ pc yields a total encounter rate (summed over 13 stellar types) of 10.5 Myr*-1* for the Sun in its current location. The density of stars increases toward the galactic center, while the velocity dispersion increases over time in galactic simulations (Minchev et al., 2012; Roškar et al., 2012). For systems placed at different galactocentric distances, we scale the number density according to the mass model of Kordopatis et al. (2015). The velocity dispersion is scaled as $\sqrt{t}$, where $t$ is the simulation time. The next encounter time is then

[TABLE]

where $\xi$ is a random number on the interval $(0,1]$ and is the time of closest approach to the host system’s primary star (see below).

When the simulation time exceeds the time of the next stellar encounter, the initial position, mass, and velocity of the passing star are selected randomly. The position is given by two random numbers that give the angles $\theta_{\star}$ and $\phi_{\star}$—these are the co-latitudinal and longitudinal position of the passing star on a sphere of radius $R_{\text{enc}}$, centered on the primary component of the host system. The passing star’s initial position, $\vec{R}_{\star}=(x_{\star},y_{\star},z_{\star})$, is then

[TABLE]

To obtain the stellar mass, we need to build a distribution of the stellar encounter frequencies for each stellar type and utilize rejection sampling. Similar to Eq. (F28), the encounter frequencies are given for each stellar type by

[TABLE]

with values again given by García-Sánchez et al. (2001). The distribution of stellar types extends across a magnitude range of 23.7 (the bin edges are given in Table 6) for the main sequence. Giant and white dwarfs are given bin sizes of $\Delta M_{V,i}=1$ each, resulting in a total magnitude range of 25.7. We now draw a random magnitude, $\xi_{M_{V}}$, on the interval $[-7.7,18]$, with $-7.7<M_{V}<-6.7$ assigned to giants and $-6.7<M_{V}<-5.7$ assigned to white dwarfs (for selection purposes only), and calculate the relative frequency of an encounter as $f_{\text{rel}}=f_{\text{enc},i}/\Delta M_{V,i}$ for the corresponding stellar type. We next draw another random number, $\xi_{\text{samp}}$, on the interval $[0,f_{\text{enc}}]$, where $f_{\text{enc}}$ is the total encounter rate and sets the threshold for acceptance. If $f_{\text{rel}}<\xi_{\text{samp}}$, we then redraw $\xi_{M_{V}}$ and $\xi_{\text{samp}}$ until $f_{\text{rel}}>\xi_{\text{samp}}$. Once a value for the magnitude is accepted, the mass of the passing star is calculated using the formulae in Reid et al. (2002) for the main sequence. White dwarfs are assigned a mass of 0.9 $M_{\odot}$ and giants a mass of 4 $M_{\odot}$. The white dwarf and giant masses are archetypal values from Allen’s Astrophysical Quantities (Cox, 2000). Following Rickman et al. (2005, 2008) and Kaib et al. (2013), we have used these discrete values as a crude first approximation. Future versions of VPLanet will sample their masses from a distribution.

The three components of the velocity of the star are selected from a normal distribution with width $\sigma_{\star i}/\sqrt{3}$. This velocity is relative to the local standard of rest, and thus we must account for the host system’s apex velocity. The velocity between the host and the passing star is then

[TABLE]

where $\vec{v}_{\star}$ is the passing star’s velocity as randomly selected and $\vec{v}_{h}$ is the host system’s apex velocity relative to the selected star’s type. To account for the additional effect of the host system’s velocity on the encounter rate (and selection process), we reject and redraw the velocity if a random number on $(0,1]$ is greater than $|\vec{v}_{h,\star}|/v_{\text{max}}$, where $v_{\text{max}}=v_{h}+3\sigma_{\star i}$ (see Rickman et al., 2008). We similarly reject and redraw $\vec{v}_{h,\star}$ if its radial component $(\vec{v}_{h,\star})_{r}$ is positive, in which case the passing star would be moving away from the host system.

Our next step is to advance the mean anomaly, $l$, of the orbiting body in the host system (whether it is a comet or secondary star). This value is calculated from $l(t)=l(t-\Delta t)+n\Delta t$, where $n=\sqrt{\mu/a^{3}}=\mu^{2}/L^{3}$. The mean motion, $n$, is merely the first term in Eq. (F5). The model is not fully self-consistent in that the second term in Eq. (F5) (the effect of the tide on the mean anomaly) is not included (indeed, this is what makes the model “secular”), however, this term is generally very small compared to $n$ (Heisler & Tremaine, 1986).

Next, we calculate the impact parameters of the passing star. There are two that interest us: the closest approach to the primary (or central) star of the host system, $\vec{b}_{1}$, and the closest approach to the orbiter, $\vec{b}_{2}$, which do not occur at the same instant in time. The time of the encounter, $t_{\text{next}}$, is taken to be the time of closest approach to the primary. The impact parameters are calculated using simple kinematics assuming that the position of each body in the host system is fixed at the time of the encounter, $t_{\text{next}}$, and that the passing star’s velocity is constant, i.e., the encounter timescale is much shorter than the orbital timescale. This approach translates to

[TABLE]

where $\vec{R}_{2}$ is the position of orbiting body of the host system at time $t_{\text{next}}$ and

[TABLE]

Finally, we apply the perturbation using the impulse approximation (Remy & Mignard, 1985). This takes the form of an instantaneous change in the orbiter’s velocity, given by

[TABLE]

where $m_{\star}$ is the mass of the passing star and $v_{\star,2}$ is the relative speed between the orbiter and the passing star. We then recalculate the osculating elements based on the new velocity and finally recalculate $\vec{e}$ and $\vec{J}$. In summary, the path of the passing star is as follows: the star begins at some random location a distance of $R_{\text{enc}}$ away from the host, moves toward the host and companion until it reaches the distances $b_{1}$ and $b_{2}$, where it interacts gravitationally with each body, then proceeds to leave the encounter area.

We then randomly select the time until the next encounter, $t_{\text{next}}$, using the procedure previously described, and then continue integration of the tidal model, etc., until $t_{\text{next}}$ is exceeded.

F.3 Radial migration

Radial migration of the host system is treated as a sudden change in the galactocentric position, since galactic simulations show the process is typically very rapid (Roškar et al., 2008). In this case, the encounter rate is recalculated at the time of migration, as is the local mass density, $\rho_{0}$. Currently, a system is only able to migrate once per simulation.

Appendix G The POISE Module

The climate model, POISE (Planetary Orbit-Influenced Simple EBM), is a one-dimensional EBM based on North & Coakley (1979), with a number of modifications, foremost of which is the inclusion of a model of ice sheet growth, melting, and flow. The model is one-dimensional in $x=\sin{\phi}$, where $\phi$ is the latitude. In this fashion, latitude cells of size $dx$ do not have equal width in latitude, but are equal in area. The general energy balance equation is

[TABLE]

where $C(x)$ is the heat capacity of the surface at location $x$, $T$ is the surface temperature, $t$ is time, $D$ is the coefficient of heat diffusion between latitudes (due to atmospheric circulation), $I(x,t)$ is the outgoing long-wave radiation (OLR) to space (i.e., the thermal infrared flux), $S(x,t)$ is the incident instellation (stellar flux), and $\alpha$ is the planetary albedo and represents the fraction of the instellation that is reflected back into space. Note that POISE assumes a constant atmospheric mass, and so should not be coupled to AtmEsc.

Though the model lacks a true longitudinal dimension, each latitude is divided into a land portion and a water portion. The land and water have distinct heat capacities and albedos, and heat is allowed to flow between the two regions. The energy balance equation can then be separated into two equations (North & Coakley, 1979; Deitrick et al., 2018b):

[TABLE]

where we have employed the co-latitudinal component of the spherical Laplacian, $\nabla^{2}$ (the radial and longitudinal/azimuthal components vanish). The effective heat capacity of the ocean is $C_{W}^{eff}=m_{d}C_{W}$, where $m_{d}$ is an adjustable parameter representing the mixing depth of the ocean. The parameter $\nu$ is used to adjust the land-ocean heat transfer to reasonable values, and $f_{L}$ and $f_{W}$ are the fractions of each latitude cell that are land and ocean, respectively.

The instellation received as a function of latitude, $\phi$, and declination of the host star, $\delta$, is calculated using the formulae of Berger (1978). Declination, $\delta$, varies over the course of the planet’s orbit for nonzero obliquity. For Earth, for example, $\delta\approx 23.5^{\circ}$ at the northern summer solstice, $\delta=0^{\circ}$ at the equinoxes, and $\delta\approx-23.5^{\circ}$ at the northern winter solstice. Because $\delta$ is a function of time (or, equivalently, orbital position), the instellation varies, and gives rise to the seasons (again, assuming the obliquity is nonzero). For latitudes and times where there is no sunrise (e.g., polar darkness during winter):

[TABLE]

while for latitudes and times where there is no sunset:

[TABLE]

and for latitudes with a normal day/night cycle:

[TABLE]

Here, $S_{\star}$ is the solar/stellar constant, $\rho$ is the distance between the planet and host star normalized by the semi-major axis (i.e., $\rho=r/a$), and $H_{0}$ is the hour angle of the star at sunrise and sunset, and is defined as:

[TABLE]

The declination of the host star with respect to the planet’s celestial equator is a simple function of its obliquity $\varepsilon$ and its true longitude $\theta$:

[TABLE]

See also Laskar et al. (1993) for a comprehensive derivation. For these formulae to apply, the true longitude should be defined as $\theta=f+\Delta^{*}$, where $f$ is the true anomaly (the angular position of the planet with respect to its periastron) and $\Delta^{*}$ is the angle between periastron and the planet’s position at its northern spring equinox, given by

[TABLE]

Above, $\varpi$ is the longitude of periastron, and $\psi$ is the precession angle. Note that we add $180^{\circ}$ because of the convention of defining $\psi$ based on the vernal point, $\vernal$, which is the position of the sun at the time of the northern spring equinox.

A point of clarification is in order: EBMs (at least, the models employed in VPLanet) can be either seasonal or annual. The EBM component of POISE is a seasonal model—the variations in the instellation throughout the year/orbit are resolved and the temperature of the surface at each latitude varies in response, according to the leading terms in Equations (G2) and (G3). In an annual model, the instellation at each latitude is averaged over the year, and the energy balance equation, Eq. (G1), is forced into “steady state” by setting $\partial T/\partial t$ equal to zero (either numerically or analytically). By “steady state”, we mean that the surface conditions (temperature and albedo) come to final values and remain there. Seasonal EBMs, on the other hand, can be in a stable balance, in that the orbit-averaged surface conditions remain the same from year to year, but the surface conditions vary throughout the year.

The planetary albedo is a function of surface type (land or water), temperature, and zenith angle. For land grid cells, the albedo is:

[TABLE]

while for water grid cells it is:

[TABLE]

where $Z$ is the zenith angle of the sun at noon and $P_{2}(x)=1/2(3x^{2}-1)$ (the second Legendre polynomial). This last quantity is used to approximate the additional reflectivity seen at shallow incidence angles, e.g., at high latitudes on Earth. The zenith angle at each latitude is given by

[TABLE]

The albedos, $\alpha_{L}$, $\alpha_{W}$ (see Table 7), not accounting for zenith angle effects, are chosen to match Earth data (North & Coakley, 1979) and account, over the large scale, for clouds, various surface types, and water waves. Additionally, the factor of $0.08$ in Equations (G10) and (G11) is chosen to reproduce the albedo distribution in North & Coakley (1979). The ice albedo, $\alpha_{i}$, is a single value that does not depend on zenith angle due to the fact that ice tends to occur at high zenith angle, so that the zenith angle is essentially already accounted for in the choice of $\alpha_{i}$. Equation (G10) indicates that when there is ice on land ($M_{\text{ice}}>0$), or the temperature is below freezing, the land takes on the albedo of ice. Though there are multiple conditionals governing the albedo of the land, in practice the temperature condition is only used when ice sheets are turned off in the model, since ice begins to accumulate at $T=0^{\circ}$ C, so that ice is always present when $T<-2^{\circ}$ C. Equation (G11) indicates a simpler relationship for the albedo over the oceans: when it is above freezing, the albedo is that of water (accounting also for zenith angle effects); when it is below freezing, the albedo is that of ice, $\alpha_{i}$.

The land fraction and water fraction are constant across all latitudes, which is roughly equivalent to a single continent that extends from pole to pole. The current version of POISE does not contain other geographies because we have not yet developed a consistent method to handle the flow of ice when land fraction varies.

North & Coakley (1979) utilized a linearization of the OLR with temperature:

[TABLE]

where, for Earth, $A=203.3$ W m*-2* and $B=2.09$ W m*-2* ∘C-1, and $T$ is the surface temperature in *∘*C. This linearization is a good fit to the observations of Earth (Warren & Schneider, 1979). A side benefit is that it allows the coupled set of equations to be formulated as a matrix problem that can be solved using an implicit Euler scheme (Press et al., 1987) with the following form:

[TABLE]

where $T_{n}$ is a vector containing the current surface temperatures, $T_{n+1}$ is a vector representing the temperatures to be calculated, and $C$, $A$, $S$, and $\alpha$ are vectors containing the heat capacities, OLR offsets (Eq. [G13]), instellation at each latitude, and albedos, respectively. The matrix $\mathscr{M}$ contains all of the information on the left-hand sides of Eqs. (G2) and (G3) related to temperature. The time-step, $\Delta t$, is chosen so that conditions do not change significantly between steps, resulting in typically 60 to 80 time-steps per orbit. The new temperature values can then be calculated by taking the dot-product of $\mathscr{M}^{-1}$ with the right-hand side of Eq. (G14). The large time step allowed by this integration scheme greatly speeds the climate model, permitting simulations for millions of years in hours of computation time.

We model ice accumulation and ablation in a similar fashion to Armstrong et al. (2014). Ice accumulates on land at a constant rate, $r_{\text{snow}}$, when temperatures are below 0*∘* C. Melting/ablation occurs when ice is present and temperatures are above 0*∘* C, according to the formula:

[TABLE]

where $M_{\text{ice}}$ is the surface mass density of ice, $\sigma=5.67\times 10^{-8}$ W m*-2* K*-4* is the Stefan-Boltzmann constant, $L_{h}$ is latent heat of fusion of ice, $3.34\times 10^{5}$ J kg*-1* and $T_{\text{freeze}}=273.15$ K. The factor $\xi$ is used to scale the ice ablation. This can be done to reproduce values similar to Earth ($\sim$3 mm ∘C-1 day*-1* (see Braithwaite & Zhang, 2000; Lefebre et al., 2002; Huybers & Tziperman, 2008) or to experiment with other melting rates.

The ice sheets flow across the surface via deformation and sliding at the base. We use the formulation from Huybers & Tziperman (2008) to model the changes in ice height due to these effects. Bedrock depression is important in this model (despite the fact that we have only one atmospheric layer and thus do not resolve elevation-based effects), because the flow rate is affected. The ice flow (via Huybers & Tziperman, 2008) is:

[TABLE]

where $h$ is the height of the ice, $H$ is the height of the bedrock (always negative or zero, in this case), $A_{\text{ice}}$ represents the deformability of the ice, $\rho_{i}$ is the density of ice, $g$ is the acceleration due to gravity, and $n$ is the exponent in Glen’s flow law (Glen, 1958), where $n=3$. The ice height and ice surface mass density, $M_{\text{ice}}$ are simply related via $M_{\text{ice}}=\rho_{i}h$. The first term inside the derivative represents the ice deformation; the second term is the sliding of the ice at the base. The latitudinal coordinate, $y$, is related to the radius of the planet and the latitude, $y=R\phi$, thus $\Delta y=R\Delta x(1-x^{2})^{-1/2}$. Finally, $u_{b}$, the ice velocity across the sediment, is:

[TABLE]

as described by Jenson et al. (1996). The constant $D_{0}$ represents a reference deformation rate for the sediment, $\mu_{0}$ is the reference viscosity of the sediment, $h_{s}$ is the depth of the sediment, and $m=1.25$. The shear stress from the ice on the sediment is:

[TABLE]

The shear strength increases linearly with depth, with a constant rate of increase, $b_{\text{sed}}$, given by the material properties of the soil:

[TABLE]

where $\rho_{s}$ and $\rho_{w}$ are the density of the sediment and water, respectively, and $\phi_{s}$ is the internal deformation angle of the sediment. We adopt the same numerical values as Huybers & Tziperman (2008) for all parameters related to ice and sediment (see Table 8), with a few exceptions. We use a value of $A_{\text{ice}}$ (ice deformability) that is consistent with ice at 270 K (Paterson, 1994), and a value of $r_{\text{snow}}$ (the precipitation rate) that best reproduces Milankovitch cycles on Earth (see Section 14.1). Note also that the value of $D_{0}$ in Table A2 of Huybers & Tziperman (2008) appears to be improperly converted for the units listed (the correct value, from Jenson et al. (1996), is listed in the text, however). With Eqs. (G17) and (G18), Eq. (G16) can be treated numerically as a diffusion equation, with the form:

[TABLE]

where,

[TABLE]

and $D_{\text{ice}}$ is evaluated at each time-step, at every boundary to provide mass continuity. We solve the diffusion equation numerically using a Crank-Nicolson scheme (Crank et al., 1947).

The bedrock depresses and rebounds locally in response to the changing weight of ice above, always seeking isostatic equilibrium. The equation governing the bedrock height, $H$, is

[TABLE]

where $T_{b}$ is a characteristic relaxation time scale, $H_{eq}=0$ is the ice-free equilibrium height, and $\rho_{b}$ is the bedrock density (Clark & Pollard, 1998; Huybers & Tziperman, 2008). We again adopt the values used by Huybers & Tziperman (2008) (see Table 8).

Because of the longer time-scales (years) associated with the ice sheets, the growth/melting and ice-flow equations are run asynchronously in POISE. First, the EBM, Eq. (G1), is run for a user set number of orbital periods (for example, 3-5), and ice accumulation and ablation is tracked over this time frame, but ice-flow (Equation G16) is ignored. Integrating the EBM without ice flow for a few orbits is generally acceptable as long as the total time is lower than the required time step for the ice flow, which is typically around 5-10 years. The annually-averaged ice accumulation/ablation is then calculated from this time-frame and passed to the ice-flow time-step, which can be much longer (years). The EBM is then re-run periodically to update accumulation and ablation and ensure that conditions vary smoothly and continuously.

To clarify, the hierarchy of models and their time-steps (when POISE is part of a VPLanet calculation) is as follows:

The EBM (shortest time-step): run for a duration of several orbital periods with time-steps on the order of days. The model is then rerun at the end of every orbital/obliquity time-step and at user-set intervals throughout the ice-flow model. 2. 2.

The ice-flow model (middle time-step): run at the end of every orbital time-step (with time-steps of a few orbital periods), immediately after the EBM finishes. The duration of the model will follow one of two scenarios:

1. (a)

If the orbital/obliquity time-step is sufficiently long, the EBM is rerun at user-set intervals, then the ice-flow model continues. The ice-flow model and the EBM thus alternate back-and-forth until the end of the orbit/obliquity time-step. 2. (b)

If the orbital/obliquity time-step is shorter than the user-set interval, the ice-flow model simply runs until the end of the orbital time-step. 3. 3.

The rest of VPLanet (longest time-step). The time-steps are set by the fastest changing variable amongst those parameters, see $\S$ 2.

This approach is shown schematically in Figure 31. The user-set interval discussed above must be considered carefully. The assumption is that annually-averaged climate conditions like surface temperature and albedo do not change much during the time span over which the ice-flow model runs. For the results in Section 14.1, we choose a value that ensures that the ice-flow does not run so long that it changes the albedo over more than one grid point without updating the temperature and ice balance (growth/ablation) via the EBM.

The initial conditions for the EBM are as follows: The first time the EBM is run, the planet has zero ice mass on land, the temperature on both land and water is set by the function

[TABLE]

where $\phi$ is the latitude. This choice gives the planet a mean temperature of $\sim 14^{\circ}$ C, ranging from $\sim 28^{\circ}$ C in the tropics to $\sim-13^{\circ}$ at the poles, or a “warm start” condition. The initial albedo of the surface is calculated from the initial temperatures. In order to adjust the initial conditions to better suit non-Earth-like initial conditions (i.e., different insolation parameters), we then perform a “spin-up” phase, running the EBM iteratively until the mean temperature between iterations changes by $<0.1^{\circ}$ C, without running the orbit, obliquity, or ice-flow models, to bring the seasonal EBM close to equilibrium at the actual stellar flux the planet receives and its actual initial obliquity. This spin up phase generally takes a few seconds of computation time. The threshold value of $0.1^{\circ}$ was chosen so that the EBM is close to equilibrium with the stellar forcing without the spin up phase running longer than 10-20 seconds, however, a different value may be more appropriate for different scenarios. In that case, users should adjust this value and/or the initial temperature distribution to suit the nature of the planet being modeled. There does exist the possibility that, for certain values of this error threshold, the solution will oscillate indefinitely between two solutions and the spin up phase will never end. We advise users to be mindful of whether the combination of initial conditions and spin up process produces sensible results. Then, every time the EBM is rerun (at the user-set interval or the end of the orbit/obliquity time-step), the initial conditions are taken from the previous EBM run (temperature distribution) and the end of the ice-flow run (albedo, ice mass).

Appendix H The RadHeat Module

Radiogenic heat production in the Earth is generated primarily by the decay of ${}^{238}U$, ${}^{235}U$, ${}^{232}Th$, and ${}^{40}K$ in the crust, mantle, and core. Other species, such as 26Al, are also potential heat sources on exoplanets. The radiogenic power produced by species $i$ in reservoir $j$ is

[TABLE]

where $\lambda_{i,1/2}=\ln 2/\tau_{i,1/2}$, $\tau_{i,1/2}$ is the half-life, $t$ is time, $Q_{i.j}(0)$ is the initial heat production at $t=0$, and $j=$ core, mantle, crust. The initial radiogenic heat can be input as a power, number of atoms (typically $\sim 10^{42}$), or mass within each reservoir. Table  9 shows our default Earth values.

Appendix I The SpiNBody Module

In addition to the DistOrb method described in Section C, orbital evolution can also be modelled using a classic N-body simulation. The VPLanet module SpiNBody simulates the orbital evolution of a system by calculating the gravitational forces on each body using Newtonian gravity. The gravitational force on an orbital body $i$ by $j$ additional bodies is given by:

[TABLE]

where $G$ is Newton’s gravitational constant, $m_{i}$ is the mass of body $i$, and $r_{i,j}$ is the distance between body $i$ and body $j$. Thus, the barycentric position and velocity derivatives are given by:

[TABLE]

[TABLE]

Because this is an N-body model, it is several orders of magnitude slower than DistOrb. Unlike DistOrb, however, SpiNBody can simulate systems of large eccentricity and mutual inclination, as well as systems with orbital resonances. Note, because SpiNBody uses Cartesian coordinates and derivatives, it cannot currently couple to modules (e.g., DistOrb and EqTide) that calculate orbital evolution based on the derivatives of the osculating elements. Fig. 32 shows the evolution of Earth’s orbit due to the Sun and other 7 planets for both VPLanet in black, and HNBODY’s 4th order Runge-Kutta method (i.e., the same method as VPLanet) in red. At the end of the simulations, the difference in eccentricity is about 2%, but is a function of timestep.

Appendix J The STELLAR Module

The STELLAR module simulates the evolution of the stellar radius, radius of gyration ($r_{g}$), effective temperature, luminosity, XUV luminosity, and rotation rate over time. STELLAR evolves the first four quantities using a bicubic spline interpolation over mass and time of the evolutionary tracks of Baraffe et al. (2015) for solar-metallicity stars. STELLAR tracks the XUV luminosity of stars using the empirical broken power law model of Ribas et al. (2005),

[TABLE]

where $L_{\mathrm{bol}}$ is the total (bolometric) stellar luminosity, i.e., from the Baraffe et al. (2015) grids, $\beta_{\mathrm{XUV}}$ is the power law exponent (equal to $-1.23$ for Sun-like stars; Ribas et al., 2005), and $f_{\mathrm{sat}}$ is the initial ratio of XUV to bolometric luminosity. Prior to $t=t_{\mathrm{sat}}$, the XUV luminosity is said to be “saturated,” as observations show that the ratio $L_{\mathrm{XUV}}/L_{\mathrm{bol}}$ remains relatively constant at early times (e.g., Wright et al., 2011). For G and K dwarfs, $t_{\mathrm{sat}}\approx 100\mathrm{Myr}$ and $f_{\mathrm{sat}}\approx 10^{-3}-10^{-4}$ (Jackson et al., 2012). Because of poor constraints on the ages of field M dwarfs, the saturation timescale for these stars is not known, although it is likely longer than that for K dwarfs. Following Luger & Barnes (2015), we adopt $t_{\mathrm{sat}}=1\,\mathrm{Gyr}$ and $f_{\mathrm{sat}}=10^{-3}$ as default values in STELLAR.

In addition to tracking the evolution of fundamental stellar properties, STELLAR models the angular momentum evolution of stars undergoing stellar evolution, e.g., radius contraction, and angular momentum loss due to magnetic braking. VPLanet includes three different prescriptions from Reiners & Mohanty (2012), Repetto & Nelemans (2014), and Matt et al. (2015), with the latter including the corrections of Matt et al. (2019).

The Reiners & Mohanty (2012) magnetic braking model, which links angular momentum loss with stellar magnetic field strength, is given by

[TABLE]

where C=2.66\times 10^{3}\text{ (gm{}^{5}$cm${}^{-10}$s${}^{3}$)${}^{1/3}}, \omega_{crit}=8.56\times 10^{-6}\text{ s{}^{-1}} for $M>0.35$ M*⊙*, and \omega_{crit}=1.82\times 10^{-6}\text{ s{}^{-1}} for $M\leq 0.35M_{\odot}$. In this model, the dichotomy between unsaturated ($\omega<\omega_{crit}$) and saturated ($\omega\geq\omega_{crit}$) assumes that for fast enough rotation rates, i.e., the stars in the “saturated” regime, the stellar magnetic field strength is constant, whereas in the unsaturated regime, the stellar magnetic field strength depends on the rotation rate. The Repetto & Nelemans (2014) model is simply based off of the observed spin-down of Sun-like stars by Skumanich (1972) and is given by

[TABLE]

where $\gamma=5\times 10^{-25}\text{ s m}^{-2}$. STELLAR also allows the user to select the Matt et al. (2015) magnetic braking model that successfully reproduced observed features in the rotation period distribution of Kepler field stars. The Matt et al. (2015) model is based on scaling relations with the stellar mass, radius, and Rossby number, the ratio of the convective turnover timescale to the rotation period, and is given by the following equations:

[TABLE]

where $\tau_{cz}$ is the convective turnover timescale, $\tau_{cz\odot}$ is the solar convective turnover timescale, 12.9 days, and

[TABLE]

is the inverse critical Rossby number for saturation in solar units. Matt et al. (2015) finds that saturation occurs when $Ro\leq Ro_{\odot}/\chi$, where $\chi=10$. The normalized torque is given by

[TABLE]

We compute $\tau_{cz}$ using the fit presented in Cranmer & Saar (2011), given by

[TABLE]

for stellar effective temperature $T_{eff}$.

We model the total change in a star’s rotation rate, $\omega$ due to radius and $r_{g}$ evolution, assuming conservation of angular momentum, and magnetic braking according to

[TABLE]

for the stellar moment of inertia, $I=mr_{g}^{2}R^{2}$. We compute numerical time derivatives of the stellar radius and $r_{g}$ using our interpolation of the Baraffe et al. (2015) grids and $\dot{J}_{\star}$ is the angular momentum loss due to magnetic braking given by one of the above models. In Fig. 33, we plot the rotation period evolution for late M-dwarfs and Sun-like stars subject to Eq. (J10) for the three magnetic braking models implemented in VPLanet.

In Fig. 34, we compare our STELLAR magnetic braking implementation to the Matt et al. (2015) model by reproducing their Fig. 3. We plot the rotation period distribution of a synthetic cluster of single stars with ages ${\sim}4$ Gyr from the publicly-available VPLanet STELLAR simulations of Fleming et al. (2019)444Simulation data are available at https://github.com/dflemin3/sync/tree/master/Data and is given in the MagneticBraking VPLanet example directory (see the Fig. 34 caption). in Fig. 34 and overplot the Kepler stellar rotation period distribution measured by McQuillan et al. (2014). We find that our 4 Gyr-old synthetic cluster matches the upper envelope of the Kepler distrubtion, in good agreement with Matt et al. (2015) and validating our implementation.

Appendix K The ThermInt Module

The thermal history of the interior is determined by the time evolution of the average mantle $T_{m}$ and average core $T_{c}$ temperatures. The energetic state of the core governs inner core crystallization and the dynamo.

K.1 Governing Differential Equations

The conservation of energy in the mantle is

[TABLE]

where $Q_{surf}$ is the total mantle surface heat flow, $Q_{conv}$ is heat conducted through the lithospheric thermal boundary layer that is supplied by mantle convection, $Q_{melt}$ is heat loss due to the eruption of upwelling mantle melt at the surface, $Q_{i,man}$ is heat generated by radioactive decay of species $i$ in the mantle, $Q_{cmb}$ is heat lost from the core across the core-mantle boundary (CMB), $Q_{man}$ is the secular heat lost from the mantle, $Q_{tidal}$ is heat generated in the mantle by tidal dissipation, and $Q_{L,man}$ is latent heat released by the solidification of the mantle. Crustal heat sources have been excluded because they do not contribute to the mantle heat budget. Note that heat can be released from the mantle in two ways: via conduction through the upper mantle thermal boundary layer ($Q_{conv}$) and by melt eruption ($Q_{melt}$). Detailed expressions for heat flows and temperature profiles as functions of mantle and core properties are given below.

Similarly, the thermal evolution of the core is governed by the conservation of energy in the core,

[TABLE]

where $Q_{core}$ is core secular cooling, $Q_{i,core}$ is radiogenic heat production of species $i$ in the core, and heat released by the solidification of the inner core is $Q_{icb}=\dot{M}_{ic}(L_{icb}+E_{icb})$, where $\dot{M}_{ic}$ is the change in inner core mass $M_{ic}$, and $L_{icb}$ and $E_{icb}$ are the latent and gravitational energy released per unit mass at the inner-core boundary (ICB).

Thermal evolution equations for the average mantle $T_{m}$ and core $T_{c}$ temperatures are derived by using the secular cooling equation $Q_{j}=-c_{j}M_{i}\dot{T}_{j}$, where $c$ is specific heat and $j$ refers to either mantle or core, in Eqs. (K1) and (K2). Solving for $\dot{T}_{m}$ and $\dot{T}_{c}$ gives the mantle and core thermal evolution equations,

[TABLE]

[TABLE]

where the denominator of (K4) is the sum of core specific heat and heat released by the inner core growth, $A_{ic}$ is inner core surface area, $\rho_{ic}$ is inner core density, $\eta_{c}$ is a constant that relates average core temperature to CMB temperature, $dR_{ic}/dT_{cmb}$ is the rate of inner core growth as a function of CMB temperature, and $L_{Fe}$ and $E_{G}$ are the latent and gravitational energy released at the ICB per unit mass. Table LABEL:tab:thermint lists the values for these parameters, which are taken from McKenzie (1984), McKenzie & Bickle (1988), Hirschmann (2000), Elkins-Tanton (2008), Barnes et al. (2013), and Jaupart et al. (2015).

An additional conservation of energy can be added for the cooling of the lithosphere for a stagnant lid planet. Alternatively, the stagnant lid heat flow can be approximated by

[TABLE]

where $\epsilon_{stag}\approx 1/25$ (Driscoll & Bercovici, 2014).

K.2 Mantle Rheology

Effective mantle viscosity follows an Arrhenius Law form,

[TABLE]

where $\nu=\eta/\rho_{m}$ is kinematic viscosity, $\rho_{m}$ is mantle density, $\nu_{ref}$ is a reference viscosity, $E_{\nu}$ is the viscosity activation energy, $R_{g}$ is the gas constant, $T$ is either upper or lower mantle temperature, and $\epsilon_{phase}$ accounts for the weakening effect of a solid to liquid phase change (see Table LABEL:tab:thermint for a list of constants). Shear modulus, which is required to calculate tidal heating (see App. K.3), is similarly described,

[TABLE]

This model predicts the rapid drop in shear modulus with melt fraction demonstrated experimentally by Jackson et al. (2004). The reference shear modulus $\mu_{ref}=10^{5}$ Pa and effective stiffness $\beta_{st}=1.71\times 10^{4}$ GPa are calibrated by $k_{2}=0.3$ and $Q$= 100 for the present-day mantle.

The influence of melt fraction $\phi$ on viscosity is modeled following the parameterization of Costa et al. (2009),

[TABLE]

and

[TABLE]

where $\Phi=\phi/\phi^{*}$, “erf” is the error function, and $\phi^{*}$, $\xi$, $\gamma_{ph}$, and $\delta_{ph}$ are empirical constants (Table LABEL:tab:thermint). The nominal value of $A_{\mu}=2\times 10^{5}$ J mol*-1* produces a dissipation peak when the melt fraction is about $50\%$.

The mid-mantle and lower mantle viscosities can be related to the upper mantle viscosity by

[TABLE]

where $f_{\nu,man}$ and $f_{\nu,LM}$ are constant coefficients. The mantle Rayleigh number is

[TABLE]

where $\alpha$ is thermal expansivity, $D_{man}=R-R_{cmb}$ is mantle shell thickness, $\kappa$ is thermal diffusivity, and mantle convective temperature jump is,

[TABLE]

where $T_{UM}$ and $\Delta T_{UM}$ are the temperatures at the bottom of and across the upper mantle thermal boundary layer; $T_{LM}$ and $\Delta T_{LM}$ are the temperatures at the top of and across the lower mantle thermal boundary layer; $T_{cmb}$ is CMB temperature.

K.3 Tidal Heating

Tidal heating associated with gravitational tides is assumed to occur by visco-elastic dissipation in the mantle. The power dissipated by tidal strain in a planet in synchronous rotation is (Segatz et al., 1988; Driscoll & Barnes, 2015),

[TABLE]

where $G$ is the gravitational constant, $M_{*}$ is the central mass, $R_{p}$ is planet radius, $\omega$ is orbital frequency, $e$ is orbital eccentricity, $a$ is orbital semi-major axis, and $\mbox{Im}(k_{2})$ is the imaginary part of the complex second order love number $k_{2}$. If planetary rotation is synchronous, then the tidal frequency is equal to the mean motion $\omega=n=\sqrt{GM_{*}/a^{3}}$, and the tidal power reduces to

[TABLE]

This expression for tidal dissipation is the product of three physical components: (1) tidal efficiency ($-$Im$(k_{2})$), (2) star-planet size ($M_{*}^{5/2}R_{p}$), and (3) orbit ($e^{2}/a^{15/2}$).

The one-dimensional dissipation model in Eq. (K14) assumes a homogeneous body with uniform stiffness and viscosity. We assume that dissipation is dominated by the material properties at the base of the upper mantle thermal boundary layer where viscosity is expected to be at a minimum in the mantle. To derive the dissipation efficiency ($-$Im$(k_{2})$) we first define the Love number,

[TABLE]

where $\mu$ is shear modulus and $\beta_{st}$ is effective stiffness. Writing shear modulus as a complex number and using the constitutive relation for a Maxwell body, one can derive the dissipation efficiency in Eq. (K14) as

[TABLE]

where $\eta=\nu\rho$ is dynamic viscosity (Henning et al., 2009). We note that this model does not involve a tidal $\mathcal{Q}$ factor, rather the rheological response of the mantle is described entirely by $\mbox{Im}(k_{2})$. For comparison with other models, one can compute the standard tidal $\mathcal{Q}$ factor of the Maxwell model as

[TABLE]

The common approximation is then -Im$(k_{2})\approx k_{2}/Q$.

K.4 Geotherm

The mantle temperature profile is assumed to be adiabatic everywhere except in the thermal boundary layers where it is conductive. The adiabatic temperature profile in the well mixed region of the mantle is approximated to be linear in radius, which is a good approximation considering that mantle thickness $D=2891$ km is significantly less than the adiabatic scale height $H=c_{p}/\alpha g\approx 12650$ km,

[TABLE]

where the solid adiabatic gradient is $\gamma_{ad}\approx 0.5~{}\mbox{K/km}$. In the thermal boundary layers the conductive temperature solutions,

[TABLE]

replace the adiabat. Thermal boundary layer temperature jumps are $\Delta T_{UM}=T_{UM}-T_{g}$ and $\Delta T_{LM}=T_{cmb}-T_{LM}$, where $T_{g}$ is ground temperature, and thermal boundary layer depth is $\delta$ (Driscoll & Bercovici, 2014).

The core temperature profile is assumed to be adiabatic throughout the entire core, and the thermal boundary layers within the core are ignored. This approximation is valid because the low viscosity and high thermal conductivity of liquid iron produce very small thermal boundary layers that are insignificant. The core adiabatic profile is approximated by

[TABLE]

where $D_{N}\approx 6340$ km is an adiabatic length scale (Labrosse et al., 2001). The iron solidus is approximated by Lindemann’s Law,

[TABLE]

where $T_{Fe,0}=5600$ K, $\gamma_{c}$ is the core Gruneisen parameter, $D_{Fe}=7000$ km is a constant length scale (Labrosse et al., 2001), and $\Delta T_{\chi}$ is a liquidus depression due to the presence of a light element, e.g., potassium. If $\Delta T_{\chi}=0$, then the inner core radius can be written as a function of CMB temperature by equating Eqs. (K21) and (K22) and solving for $r=R_{ic}$, giving

[TABLE]

where $T_{cmb}$ is CMB temperature, $T_{Fe,cen}$ is the liquidus at the center of the core, and $T_{Fe,CMB}$ is the liquidus at the CMB. The time derivative of Eq. (K23) gives the rate of inner core growth as a function of core cooling rate $\dot{T}_{c}$ as,

[TABLE]

Assuming a constant inner core density $\rho_{ic}$, inner core mass growth rate is approximated by

[TABLE]

where $A_{ic}$ is the area of the inner core, and the change in inner core radius with CMB temperature can be computed by the temperature derivative of Eq. (K23), giving

[TABLE]

K.5 Mantle and Core Heat Flows

In this subsection we define the remaining heat flows that appear in the mantle and core energy balance, Eqs. (K1)—(K2). The convective cooling of the mantle $Q_{conv}$ is proportional to the temperature gradient in the upper mantle thermal boundary layer,

[TABLE]

where $A$ is surface area and $k_{UM}$ is upper mantle thermal conductivity. $Q_{conv}$ is written in terms of $T_{m}$ and the thermal boundary layer thickness $\delta_{UM}$ by requiring that the Rayleigh number of the boundary layer $Ra_{UM}$ be equal to the critical Rayleigh number for thermal convection $Ra_{c}\approx 660$ (Howard, 1966; Solomatov, 1995; Sotin & Labrosse, 1999; Driscoll & Bercovici, 2014). Solving for the critical thermal boundary layer thickness from $Ra_{UM}=Ra_{c}$ gives,

[TABLE]

Using this relationship in Eq. (K27) gives,

[TABLE]

where the thermal boundary layer temperature jump $\Delta T_{UM}$ has been replaced by $\Delta T_{UM}\approx\epsilon_{UM}\Delta T_{m}$, $\epsilon_{UM}=\exp(-(R_{UM}-R_{m})\alpha g/c_{p})\approx 0.7$, which is the adiabatic temperature decrease from the average mantle temperature to the bottom of the upper mantle thermal boundary layer, $\Delta T_{m}=T_{m}-T_{g}$, and the mantle cooling exponent is $\beta=1/3$ (Driscoll & Bercovici, 2014).

Similar to the mantle convective heat flow, the CMB heat flow is,

[TABLE]

where $A_{c}$ is core surface area and $k_{LM}$ is lower mantle thermal conductivity. The lower mantle and CMB temperatures, $T_{LM}$ and $T_{cmb}$, are extrapolations along the mantle and core adiabats: $T_{LM}=\epsilon_{LM}T_{m}$ and $T_{cmb}=\epsilon_{c}T_{c}$, where $\epsilon_{LM}=\exp(-(R_{LM}-R_{m})\alpha g/c_{p})\approx 1.3$ and $\epsilon_{c}\approx 0.8$. Analogous to Eq. (K28), the lower mantle thermal boundary layer thickness $\delta_{LM}$ is derived by assuming the boundary layer Rayleigh number is critical, giving

[TABLE]

Core secular cooling is

[TABLE]

where $M_{c}$ is core mass, $c_{c}$ is core specific heat, and $\dot{T}_{c}$ is the rate of change of the average core temperature $T_{c}$.

K.6 Mantle Melting

The mantle solidus is approximated by a third-order polynomial (Elkins-Tanton, 2008),

[TABLE]

where the coefficients are constants (see Table LABEL:tab:thermint). This solidus is calibrated to fit the following constraints: solidus temperature of $1450$ K at the surface, solidus temperature of $4150$ K at the CMB (Andrault et al., 2011), and present-day upwelling melt fraction of $f_{melt}=8\%$. The liquidus is assumed to be hotter by a constant offset $\Delta T_{liq}=500$ K, so $T_{liq}(r)=T_{sol}(r)+\Delta T_{liq}$.

Mantle melt heat loss (or advective heat flow) is modeled as,

[TABLE]

where $\epsilon_{erupt}=0.2$ is the assumed efficiency of magma eruption to Earth’s surface, $\dot{M}_{melt}$ is melt mass flux (see below), $L_{melt}$ is latent heat of the melt, $c_{m}$ is specific heat of the melt, and $\Delta T_{melt}$ is the excess temperature of the melt at the surface (see below). The latent heat released by mantle melting does not contribute to the cooling of the mantle. This formulation of heat loss is similar to the “heat pipe” mechanism invoked for Io (O’Reilly & Davies, 1981; Moore, 2003), where melt eruption is a significant source of heat loss. We note that this mechanism is more important for stagnant lid planets, where the normal conductive heat flow is lower (Driscoll & Bercovici, 2014).

The melt mass flux $\dot{M}_{melt}$ is the product of the upwelling solid mass flux times the melt mass fraction $f_{melt}$,

[TABLE]

where solid density is $\rho_{solid}$, volumetric upwelling rate is $\dot{V}_{up}=1.16\kappa A_{p}/\delta_{UM}$, $z_{UM}=R-\delta_{UM}$, and melt fraction is

[TABLE]

This model predicts a ridge melt production of $\dot{M}_{melt}=2.4\times 10^{6}~{}\mbox{kg~{}s}^{-1}$ for $\delta_{UM}=80$ km and $f_{melt}=0.1$, similar to present-day global melt production estimates (Cogné & Humler, 2004).

We define the magma ocean as the region of the mantle with temperature exceeding the liquidus. Given the geotherm in Eqs. (K18) and (K20) and the liquidus $T_{liq}(r)$ similar to Eq. (K33), the mantle will mainly freeze from the bottom of the convecting mantle up because the liquidus gradient is steeper than the adiabat (e.g., Elkins-Tanton, 2012). However, if the core is hot enough, a second melt region exists in the lower mantle boundary layer, where the temperature gradient exceeds the liquidus and the mantle freezes towards the CMB.

Latent heat released from the solidification of the mantle is

[TABLE]

where $L_{melt}$ is the latent heat released per kg and $\dot{M}_{sol}$ is the solid mantle growth rate. The growth rate is calculated assuming a uniform mantle density $\rho_{m}$ so that $\dot{M}_{sol}=\rho_{m}\dot{V}_{sol}$, where $\dot{V}_{sol}=-\dot{V}_{liq}$. The rate of change of the liquid volume of the mantle is approximated by

[TABLE]

where $\dot{T}_{m}$ is the mantle secular cooling rate and $dV_{liq}/dT_{m}$ is linearly approximated by $8\times 10^{17}~{}\mbox{m}^{3}\mbox{K}^{-1}$, which is the change in liquid volume from a $90\%$ liquid to a completely solid mantle. This approximation implies that the latent heat released due to mantle solidification is linearly proportional to the mantle secular cooling rate, and the ratio of the latent heat flow to the mantle secular cooling heat flow is $Q_{L,man}/Q_{sec,m}\approx 0.24$. For example, a mantle solidification time of 100 Myr corresponds to an average latent heat release of $Q_{L,man}\approx 400$ TW over that time.

K.7 Core Dynamo

Given the thermal cooling rate of the core, the magnetic dipole moment $\mathcal{M}$ is estimated from the empirical scaling law,

[TABLE]

where $\gamma_{d}$ is the saturation constant for fast rotating dipolar dynamos, $\mu_{0}=4\pi\times 10^{-7}\mbox{H~{}m}^{-1}$ is magnetic permeability, $D_{c}=R_{c}-R_{ic}$ is the dynamo region shell thickness, $R_{c}$ and $R_{ic}$ are outer and inner core radii, respectively, and $F_{c}$ is the core buoyancy flux (Olson & Christensen, 2006). We assume that the field is dipolar, ignoring the complicating influences of shell thickness and heterogeneous boundary conditions. In this formulation a positive buoyancy flux implies dynamo action, which is a reasonable approximation when the net buoyancy flux is large, but may overestimate the field strength at low flux. The total core buoyancy flux $F_{c}$ is the sum of thermal and compositional buoyancy fluxes,

[TABLE]

where the thermal and compositional buoyancy fluxes are

[TABLE]

where the subscript $c$ refers to bulk core properties, core convective heat flux is $q_{c,conv}=q_{cmb}-q_{c,ad}$, gravity at the ICB is approximated by $g_{ic}=g_{c}R_{ic}/R_{c}$, and the outer core compositional density difference is $\Delta\rho_{\chi}=\rho_{c}-\rho_{\chi}$ with $\rho_{\chi}$ the light element density. For simplicity, the expression for light element buoyancy Eq. (K42) ignores buoyancy due to latent heat release at the ICB because it is a factor of $3.8$ less than buoyancy of the light elements.

The isentropic core heat flux at the CMB, proportional to the gradient of Eq. (K21), is

[TABLE]

where core thermal conductivity is approximated by the Wiedemann-Franz law,

[TABLE]

and electrical conductivity is $\sigma_{c}$ and $L_{c}$ is the Lorentz number. For typical values of high pressure-temperature iron, $\sigma_{c}=10\times 10^{5}$ $\Omega^{-1}m^{-1}$ (Pozzo et al., 2012; Gomi et al., 2013), $L_{c}=2.5\times 10^{-8}~{}\mbox{W}\Omega\mbox{K}^{-1}$, and $T_{cmb}=4000$ K, the core thermal conductivity is $k_{c}=100~{}\mbox{Wm}^{-1}\mbox{K}^{-1}$.

Appendix L List of Symbols

In this appendix we present the definitions of all symbols used in this manuscript. Those that are general are listed first, followed by those that are specific to individual models.

\twocolumngrid

L.1 General Symbols

$a$ = semi-major axis

$e$ = eccentricity

$G$ = Newton’s gravitational constant

$j$ = index

$L$ = luminosity

$M$ = mass

$M_{*}$ = stellar mass

$M_{p}$ = planetary mass

$M_{j}$ = mass of body $j$

$n$ = mean motion (orbital frequency)

$\Omega$ = rotational frequency

$\Omega_{j}$ = rotational frequency of body $j$

$P$ = period

$\psi$ = obliquity

$\psi_{j}$ = obliquity of body $j$

$R$ = body radius

$R_{p}$ = planetary radius

$r_{g}$ = radius of gyration

$T_{eff}$ = effective temperature

$t$ = time

L.2 AtmEsc Symbols

XUV = X-ray + UV luminosity

$b_{\mathrm{diff}}$ = binary diffusion coefficient

$\epsilon_{\mathrm{XUV}}$ = XUV absorption efficiency

$F_{\mathrm{diff}}$ = diffusion-limited particle escape flux

$F_{\mathrm{EL}}$ = energy-limited particle escape flux

$F_{\mathrm{H}}$ = hydrogen particle escape flux

$\mathcal{F}_{\mathrm{XUV}}$ = XUV flux

$k_{\mathrm{boltz}}$ = Boltzmann constant

$K_{\mathrm{tide}}$ = tidal enhancement factor

$m_{\mathrm{c}}$ = crossover mass

$m_{\mathrm{H}}$ = hydrogen atom mass

$m_{\mathrm{O}}$ = oxygen atom mass

$T_{\mathrm{flow}}$ = temperature of hydrodynamic flow

$X_{\mathrm{O}}$ = oxygen mixing ratio

L.3 BINARY Symbols

CBP = circumbinary planet

$m_{i}$ = mass of $i^{th}$ star

$a_{AB}$ = binary semi-major axis

$e_{AB}$ = binary eccentricity

$M_{b}$ = binary mean anomaly

$\varpi_{B}$ = binary longitude of the periapse

$\alpha_{i}$ = normalized CBP orbital distance from the $i^{th}$ star

$\delta_{k0}$ = Kroenker delta function

$b^{k}_{(j+1)/2}$ = Laplace coefficient

$\Phi_{jk0}$ = 0th component of binary gravitational potential

$\Phi_{jk1}$ = 1st component of binary gravitational potential

$R$ = radial distance from binary center of mass

$R_{0}$ = CBP radial guiding center

$R_{1}$ = time-dependent epicyclic radial deviation from the guiding center radius

$\psi_{0}$ = arbitrary phase offset in the $R$ direction

$\varphi$ = arbitrary phase offset in the $\phi$ direction

$\xi$ = arbitrary phase offset in the $z$ direction

$n_{k}$ = CBP Keplerian mean motion

$n_{0}$ = CBP mean motion

$n_{AB}$ = binary mean motion

$e_{\textrm{free}}$ = binary free eccentricity

$i_{\textrm{free}}$ = binary free inclination

$\kappa_{0}$ = Radial CBP epicyclic frequency

$\nu_{0}$ = Vertical CBP epicyclic frequency

L.4 DistOrb Symbols

$e$ = orbital eccentricity

$\varpi$ = longitude of pericenter

$i$ = orbital inclination

$\Omega$ = longitude of ascending node

$h$ = first Poincaré coordinate for eccentricity

$k$ = second Poincaré coordinate for eccentricity

$p$ = first Poincaré coordinate for inclination

$q$ = second Poincaré coordinate for inclination

$\mathcal{R}$ = disturbing function

$n$ = mean motion

$\delta_{R}$ = post-Newtonian correction factor

$c$ = speed of light

$\mu$ = mass factor of inner planet in a pair

$\mu^{\prime}$ = mass factor of outer planet in a pair

$\kappa$ = Gaussian gravitational constant

$\mathcal{R_{D}}$ = direct terms disturbing function

$D0.i$ = term $i$ of direct disturbing function

$f_{i}$ = semi-major axis functions in disturbing function

L.5 DistRot Symbols

$\psi$ = precession angle

$\varepsilon$ = obliquity of planet’s spin axis

$R(\varepsilon)$ = precession rate due to stellar torque

$p$ = first Poincaré coordinate for inclination

$q$ = second Poincaré coordinate for inclination

$A(p,q)$ = precession/obliquity term from orbital plane motion

$B(p,q)$ = precession/obliquity term from orbital plane motion

$\Gamma(p,q)$ = precession/obliquity term from orbital plane motion

$p_{g}$ = geodetic (relativistic) precession rate

$\kappa$ = Gaussian gravitational constant

$M_{\star}$ = mass of host star

$M$ = mass of planet

$J_{2}$ = gravitational quadrupole moment due to planetary oblateness

$a$ = semi-major axis

$\nu$ = rotation rate of planet

$C$ = polar moment of inertia of planet

$r$ = equatorial radius of planet

$S_{0}$ = correction term for precession due to eccentric orbit

$\Lambda$ = angle between vernal point and ascending node

$\Omega$ = longitude of ascending node

$\vernal$ = vernal point (location of sun on northern spring equinox)

$\xi$ = first rectangular coordinate for obliquity and precession

$\zeta$ = second rectangular coordinate for obliquity and precession

$\chi$ = third rectangular coordinate for obliquity and precession

$J_{2\oplus}$ = gravitational quadrupole moment of Earth

$\nu_{\oplus}$ = rotation rate of Earth

$R_{\oplus}$ = equatorial radius of Earth

$M_{\oplus}$ = mass of Earth

L.6 EqTide Symbols

$\varepsilon$ = sign of the phase lag

$k_{2}$ = Love number of degree 2

$P_{eq}$ = equilibrium spin period

$Q$ = tidal quality factor

$\xi$ = constant in EqTide calculations

$Z$ = constant in EqTide calculations

L.7 GalHabit Symbols

$Z$ = height above galactic midplane

$H_{\text{av}}$ = time-averaged Hamiltonian for galactic tide

$\mu$ = mass factor

$G$ = Newtonian gravitational constant

$\rho_{0}$ = density of stars

$M_{c}$ = mass of central star

$M$ = mass of orbiter

$a$ = semi-major axis

$L$ = canonical momentum associated with orbital energy

$J$ = canonical momentum associated with total orbital angular momentum

$J_{z}$ = canonical momentum associated with $Z$ angular momentum

$\omega$ = argument of periastron

$i$ = inclination of system with respect to galactic plane

$\Omega$ = longitude of ascending node

$l$ = mean anomaly

$\vec{J}$ = angular momentum vector

$J_{x}$ = x component of angular momentum

$J_{y}$ = y component of angular momentum

$\vec{e}$ = eccentricity vector

$e_{x}$ = x component of eccentricity vector

$e_{y}$ = y component of eccentricity vector

$e_{z}$ = z component of eccentricity vector

$\hat{e}$ = unit vector with direction of eccentricity vector

$\hat{J}$ = unit vector with direction of angular momentum vector

$\hat{n}$ = unit vector used in calculation of $\omega$

$\vec{n}$ = vector used in calculation of $\omega$

$\hat{Z}$ = unit vector perpendicular to galactic plane

$f_{\text{enc}}$ = frequency of stellar encounters

$R_{\text{enc}}$ = radius at which stellar encounters begin

$v_{i}$ = velocity of typical encounter with stellar type $i$

$n_{\star,i}$ = number density of stellar type $i$

$v_{h,i}$ = average velocity of system relative to stellar type $i$

$\sigma_{\star i}$ = velocity dispersion of stellar type $i$

$t_{\text{next}}$ = time of next encounter

$\xi$ = random number used in $t_{\text{next}}$ calculation

$\vec{R}_{\star}$ = starting position of passing star

$x_{\star}$ = x component of passing star’s starting position

$y_{\star}$ = y component of passing star’s starting position

$z_{\star}$ = z component of passing star’s starting position

$x_{\star}$ = x component of passing star’s starting position

$\theta_{\star}$ = co-latitudinal angle of passing star’s starting position

$\phi_{\star}$ = longitudinal angle of passing star’s starting position

$f_{\text{enc},i}$ = frequency of stellar encounters of stellar type $i$

$M_{V}$ = magnitude of passing star

$\Delta M_{V,i}$ = magnitude bin size of stellar type $i$

$f_{\text{rel}}$ = frequency of stellar encounters of stellar type $i$ per magnitude

$\xi_{M_{V}}$ = randomly drawn magnitude of passing star

$\xi_{\text{samp}}$ = random number used in rejection sampling for magnitude

$\vec{v}_{h,\star}$ = relative velocity vector between system and passing star

$\vec{v}_{\star}$ = passing star’s velocity vector

$\vec{v}_{h}$ = velocity of system relative to local standard of rest for stellar type $i$

$v_{\text{max}}$ = threshold velocity used for passing star velocity selection

$n$ = mean anomaly of orbiter

$\vec{b}_{1}$ = impact parameter of passing star relative to host system’s primary star

$\vec{b}_{2}$ = impact parameter of passing star relative to host system’s orbiter

$\vec{R}_{2}$ = position of orbiter relative to primary

$\Delta t_{1}$ = time difference until closest approach of passing star to primary

$\Delta t_{2}$ = time difference until closest approach of passing star to orbiter

$\Delta\vec{v}_{2}$ = instantaneous change in velocity of orbiter due to passing star

$m_{\star}$ = mass of passing star

L.8 POISE Symbols

OLR = outgoing longwave radiation

$\phi$ = latitude

$x$ = sine of latitude

$C$ = heat capacity of atmosphere

$T$ = temperature at/near surface

$D$ = coefficient of heat diffusion

$I$ = OLR

$S$ = instellation

$\alpha$ = albedo

$C_{L}$ = heat capacity over land

$C_{W}^{eff}$ = heat capacity over water

$T_{L}$ = temperature over land

$T_{W}$ = temperature over water

$f_{L}$ = land fraction of grid cell

$f_{W}$ = water fraction of grid cell

$\nu$ = land-ocean heat transport parameter

$C_{W}$ = head capacity over water per meter of depth

$m_{d}$ = depth of ocean mixing layer

$\delta$ = declination of sun/host star

$S_{\star}$ = flux at planet’s distance from host star

$\rho$ = planet-star distance, in semi-major axis units

$H_{0}$ = hour angle of host star at sunrise/sunset

$\varepsilon$ = obliquity

$\theta$ = true longitude of planet

$f$ = true anomaly of planet

$\Delta^{*}$ = angle between periastron and planet’s position at northern spring equinox

$\varpi$ = longitude of periastron

$\psi$ = precession angle

$\vernal$ = vernal point

$\alpha_{L}$ = albedo over land

$\alpha_{W}$ = albedo over water

$P_{2}$ = second Legendre polynomial

$M_{\text{ice}}$ = mass of ice on land

$Z$ = zenith angle of sun/star at noon

$A$ = OLR coefficient

$B$ = OLR coefficient

$\mathscr{M}$ = matrix used in energy balance model integration

$\xi$ = ice ablation tuning parameter

$\sigma$ = Stefan-Boltzmann constant

$L_{h}$ = latent heat of fusion for water

$T_{\text{freeze}}$ = freezing temperature of water

$h$ = height of ice sheet above ground

$H$ = height of bedrock

$y$ = latitudinal coordinate

$R$ = radius of planet

$A_{\text{ice}}$ = deformability of ice

$\rho_{i}$ = density of ice

$g$ = gravitational acceleration at the surface

$n$ = Glen’s flow law exponent

$u_{b}$ = flow speed of ice at the base

$D_{0}$ = reference deformation rate of sediment

$a_{\text{sed}}$ = shear stress from ice on sediment

$b_{\text{sed}}$ = rate of increase of shear strength with depth of sediment

$\mu_{0}$ = reference viscosity of sediment

$m$ = coefficient used in basal flow calculation

$h_{s}$ = depth of sediment

$\rho_{s}$ = density of sediment

$\rho_{w}$ = density of water

$\phi_{s}$ = internal deformation angle of sediment

$D_{\text{ice}}$ = diffusion coefficient for ice flow

$r_{\text{snow}}$ = rate of snow accumulation

$T_{b}$ = relaxation time scale of bedrock

$H_{eq}$ = ice-free equilibrium height of bedrock

$\rho_{b}$ = bedrock density

$T_{0}$ = initial temperature distribution of surface

L.9 RadHeat Symbols

$Q_{i,j}$ = radiogenic heat production of species $i$ in reservoir $j$

$\lambda_{i,1/2}$ = radiogenic decay constant of species $i$ = $\ln 2/\tau_{i,1/2}$

$\tau_{i,1/2}$ = radiogenic halflife of species $i$

L.10 SpiNBody Symbols

$F_{i}$ = Newtontian gravitational force on orbital body $i$

$m_{i}$ = Mass of orbital body $i$

$r_{i,j}$ = Distance between $i$ and $j$

$\dot{v_{i}}$ = Net gravitational acceleration of body $i$

$\dot{x_{i}}$ = Net velocity of body $i$

L.11 STELLAR Symbols

$\beta_{\mathrm{sat}}$ = XUV power law index

$L_{\mathrm{bol}}$ = Stellar bolometric (total) luminosity

$L_{\mathrm{XUV}}$ = Stellar XUV luminosity

$f_{\mathrm{sat}}$ = Saturation ratio of $L_{\mathrm{XUV}}$ to $L_{\mathrm{bol}}$

$t_{\mathrm{sat}}$ = XUV saturation timescale

$r_{g}$ = Stellar radius of gyration

L.12 ThermInt Symbols

CMB = core mantle boundary

ICB = inner core boundary

$Q_{surf}$ = total mantle surface heat flow

$Q_{conv}$ = heat conducted through the lithospheric thermal boundary layer that is supplied by mantle convection

$Q_{melt}$ = heat loss due to the eruption of upwelling mantle melt at the surface

$Q_{cmb}$ = heat lost from the core conducted across the CMB

$Q_{man}$ = secular (sensible) heat lost from the mantle

$Q_{tidal}$ = heat generated in the mantle by tidal dissipation

$Q_{L,man}$ = latent heat released by the solidification of the mantle

$Q_{core}$ = core secular cooling

$Q_{i,core}$ = radiogenic heat production of species $i$ in the core

$Q_{icb}$ = heat released by the solidification of the inner core

$M_{ic}$ = inner core mass

$L_{icb}$ = latent energy released per unit mass at the ICB

$E_{icb}$ = gravitational energy released per unit mass at the ICB

$T_{g}$ = ground temperature

$T_{UM}$ = upper mantle temperature

$T_{m}$ = average mantle temperature

$T_{LM}$ = lower mantle temperature

$T_{CMB}$ = core-mantle boundary temperature

$T_{c}$ = average core temperature

$\delta_{UM}$ = thermal boundary layer thickness of upper mantle

$\delta_{LM}$ = thermal boundary layer thickness of lower mantle

$\epsilon_{erupt}$ = mantle melt mass extrusive eruption fraction

$A$ = Surface area of planet

$A_{CMB}$ = CMB area

$A_{ICB}$ = ICB area

$A_{\nu}$ = Viscosity activation energy in (K6)

$A_{\mu}$ = Nominal shear modulus activation energy in (K7)

$A_{sol}$ = Solidus coefficient in (K33)

$\alpha$ = Thermal expansivity of mantle

$\alpha_{c}$ = Thermal expansivity of core

$B$ = Melt fraction coefficient in (K8)

$B_{sol}$ = Solidus coefficient in (K33)

$\beta$ = Convective cooling exponent in (K29)

$\beta_{st}$ = Effective mantle stiffness

$c_{m}=$ = Specific heat of mantle

$c_{c}$ = Specific heat of core

$C_{sol}$ = Solidus coefficient in (K33), calibrated

$D$ = Mantle depth

$D_{Fe}$ = Iron solidus length scale

$D_{N}$ = Core adiabatic length scale

$D_{sol}$ = Solidus coefficient in (K33), calibrated

$\delta_{ph}$ = Rheology phase coefficient in (K8, K9)

$E_{G}$ = Gravitational energy density release at ICB

$\epsilon_{UM}$ = Upper mantle adiabatic temperature drop

$\epsilon_{LM}$ = Lower mantle adiabatic temperature jump

$\epsilon_{c}$ = Average core to CMB adiabatic temperature drop

$\epsilon_{phase}$ = mantle rheology melt reduction factor

$F_{c}$ = total core buoyancy flux

$F_{th}$ = core thermal buoyancy flux

$F_{\chi}$ = core compositional buoyancy flux

$\phi^{*}$ = Rheology phase coefficient in (K8, K9)

$g_{UM}$ = Upper mantle gravity

$g_{LM}$ = Lower mantle gravity

$g_{c}$ = CMB gravity

$\gamma_{c}$ = Core Gruneisen parameter

$\gamma_{dip}$ = Magnetic dipole intensity coefficient in (K39)

$\gamma_{ph}$ = Rheology phase coefficient in (K8, K9)

$\gamma_{ad}$ = solid mantle adiabatic temperature gradient

$k_{2}$ = real second order Love number in (K15)

$Im(k_{2})$ = imaginary second order Love number in (K16)

$k_{UM}$ = Upper mantle thermal conductivity

$k_{LM}$ = Lower mantle thermal conductivity

$\kappa$ = Mantle thermal diffusivity

$k_{c}$ = core thermal conductivitiy

$L_{Fe}$ = Latent heat of inner core crystallization

$L_{melt}$ = Latent heat of mantle melting

$L_{e}$ = Lorentz number

$M_{m}$ = Mantle mass

$M_{sol}$ = mass of solid mantle

$V_{liq}$ = volume of liquid mantle

$M_{c}$= Core mass

$M_{ic}$= Inner core mass

$\mathcal{M}$ = core dynamo dipolar magnetic moment

$\mu_{ref}$ = Reference shear modulus in (K7)

$\mu_{0}$ = Magnetic permeability

$\nu_{ref}$ = Reference viscosity

$f_{\nu,UM,LM}$ = Viscosity jump from upper to lower mantle

$q_{c,conv}$ = core convective heat flux

$q_{c,ad}$ = core adiabatic heat flux

$\mathcal{Q}$ = tidal quality factor

$R$ = Surface radius

$R_{c}$ = Core radius

$R_{ic}$ = Inner core radius

$R_{m}$ = Radius to average mantle temperature $T_{m}$

$Ra_{c}$ = Critical Rayleigh number

$\rho_{c}$ = Core density

$\rho_{ic}$ = Inner core density

$\rho_{m}$ = Mantle density

$\rho_{melt}$ = Mantle melt density

$\rho_{solid}$ = Mantle upwelling solid density

$\Delta\rho_{\chi}$ = Outer core compositional density difference

$\sigma_{c}$ = Core electrical conductivity

$T_{Fe,0}$ = Iron solidus coefficient in (K22)

$T_{Fe,cmb}$ = Iron solidus temperature at CMB

$T_{Fe,cen}$ = Iron solidus temperature at center of core

$\xi$ = Rheology phase coefficient in (K8, K9)

\onecolumngrid

Appendix M VPLanet Accessibility and Support Software

In this final appendix, we describe our approaches for accessibility and reproducibility. In general, we adhere to the recommendations of the National Acadamies 2018 report Open Science by Design555https://www.nap.edu/catalog/25116/open-science-by-design-realizing-a-vision-for-21st-century. This study identified numerous challenges in science as software and datasets grow very large. VPLanet brings together numerous models from disparate fields of science and therefore risks becoming inscrutable to many professional scientists. Moreover, VPLanet can produce large amounts of high-dimensional data that can be difficult for the user to digest. More philosophically, given that the challenge of detecting life beyond the Solar System is daunting and the potential societal impact of such a discovery is huge, it is critical that the scientific process that culminates in such an announcement be as transparent as possible. To address these issues, VPLanet has been designed and released to maximize usability and clarity, including code testing and data management. The next subsections describe Python scripts (requiring Python 3.x) to explore parameter space, visualize VPLanet output, efficiently store that output, and the git repository that hosts VPLanet.

M.1 Parameter Sweeps: VSPACE

VSPACE is a Python code that conveniently generates initial conditions for VPLanet simulations but does not run VPLanet itself. This option is ideal when the user wishes to run a large number of simulations, perhaps on a supercomputer. Initial conditions can be generated in a random (Monte-Carlo) fashion or over a regular grid. The user creates a template plain text file that contains a list of all desired input files, with a list of input options to be changed, iterated over, or randomly drawn. VSPACE operates on this text file, determines which files are to be copied to new directories, which input parameters are to be varied, and then creates a new directory that contains (as individual directories) all the desired simulations. When used in random mode, the user sets the total number of simulations and the parameters to be sampled along with details about the distributions. Current distribution types are uniform, Gaussian/normal, uniform in sine, or uniform in cosine. When used in grid mode, the user does not need to specify the total number of simulations. One simply sets all of the parameters to be varied, along with the spacing type (linear or logarithmic), and spacing size, or the total range and number of grid points. VSPACE automatically determines the number of simulations and identifies every permutation of the desired set of input parameters using the itertools library. There is no limit to the number of parameters that can be iterated over in a single set of simulations, thus the user should be cautious about generating large numbers of files. The examples/IoHeat directory in the VPLanet repository contains a VSPACE example.

M.2 Plotting: VPLot

The VPLot package is recommended for interacting with and visualizing the output of VPLanet simulations. VPLot is coded in Python and is installable via the pip command or from source from its own GitHub repository666https://github.com/VirtualPlanetaryLaboratory/vplot. VPLot reads .log and .forward files generated by a VPLanet run in the current working directory and provides an object-oriented interface to the simulation output, with convenient plotting functions. VPLot can also be run from the command line, providing a quick graphical view of the simulation output. For more information on VPLot, please see its online documentation. Note that all plots in this document that display VPLanet output have utilized VPLot.

M.3 Data Management: BigPlanet

BigPlanet is an accompanying Python package for storing the output of a large number of VPLanet simulations. For large numbers of VPLanet simulations, i.e., a suite over a large grid of initial conditions, the total size of the output for all simulations can exceed what can be stored in memory, necessitating new solutions. BigPlanet solves this issue by using the HDF5 data storage format777https://support.hdfgroup.org/HDF5/ for efficiently storing full time series simulation outputs for easy access. BigPlanet leverages the hierarchical structure of the HD5F format by storing each simulation as a group, each simulation body as a subgroup of the simulation, and each variable as a subgroup of its parent body. BigPlanet wraps the HDF5 file interface in easy-to-use python code, using the H5PY888https://www.h5py.org/ python package, to allow users to easily read and manipulate VPLanet simulation outputs from storage for in-memory analysis. BigPlanet is extensively documented, and includes several example scripts and Jupyter notebooks that demonstrate its typical use cases.

M.4 The VPLanet Repository

The VPLanet repository contains the source code, example scripts, and documentation to facilitate its use. The VPLanet team continuously checks that all new updates reproduce past results through continuous integration999https://travis-ci.org/, and periodically searches for memory issues using valgrind101010http://valgrind.org/ and addresssanitizer111111https://github.com/google/sanitizers/wiki/AddressSanitizer. Thus, the master branch in this git repository is stable, and the results presented in this manuscript can all be obtained with it. For those who prefer not use git, gzipped tarballs are also available. The figures in this document were generated with version 1.0.

Community input and additions are welcome. Pull requests should be issued to the “dev” branch and will be reviewed before being merged with the master branch. As new features and modules are added to the software suite, this document will be updated and available in the Manual sub-directory of the VPLanet repository.