Convective dynamics with mixed temperature boundary conditions: why thermal relaxation matters and how to accelerate it

Astrophysical simulations of convection frequently impose different thermal boundary conditions at the top and the bottom of the domain in an effort to more accurately reflect the natural system being modeled. In this work, we study Rayleigh-B\'enard convection (RBC) with mixed (``FT'') temperature boundary conditions. We aim to understand how FT boundaries change the nature of the convective solution compared to the traditional choice of fixing the temperature at the top and bottom of the domain (``TT''). We find that simulations with FT boundaries experience a long thermal rundown while simulations with TT boundaries do not. In the relaxed state, the mean behavior of an FT simulation corresponds to an equivalent simulation with TT boundaries. The fast evolution of TT simulations can therefore be taken advantage of to rapidly relax FT simulations. We show that these findings carry over to more complicated problems through a brief examination of rotating convection. We furthermore find that FT boundaries introduce minor asymmetries into the flow fields, with the fixed flux boundary producing more extreme events than the fixed temperature boundary, but that these asymmetries do not appreciably affect the bulk flows. We conclude that thermal relaxation occurs in two stages: (1) changes to the experimental energy reservoir and (2) changes to the stratification of the experiment. Through the proper choice of boundary conditions or initial conditions, the first of these stages can be bypassed, and the second stage seems to be irrelevant for RBC.


I. INTRODUCTION
Convection is a crucial heat transport mechanism in the atmospheres and interiors of stars and planets. Numerical simulations are a commonly-used tool in studies of geophysical or astrophysical convection. These studies range from examinations of convection in the simplified Boussinesq approximation [1][2][3] to highly complex "dynamo simulations" which include magnetism and atmospheric density stratification [4,5]. Regardless of complexity, numerically simulated convection is fundamentally driven by some combination of imposed boundary conditions and internal heating profiles [6]. In studies of Boussinesq convection, the standard choice is to hold constant the temperature difference across the domain by fixing the temperature at the upper and lower boundaries. However, a common choice of thermal boundary conditions in astrophysical convection [7][8][9][10][11][12][13] is to fix the flux entering the domain through the bottom boundary and to fix the value of a thermodynamic quantity (e.g., temperature or entropy) at the top boundary. We are unaware of any study which has examined the consequences of imposing the "mixed" boundaries that are frequently favored in astrophysical convection studies.
In this work, we examine how the choice of using "mixed" boundaries affects the evolved nonlinear convective state in the simplest possible model: Rayleigh-Bénard convection (RBC) under the Boussinesq approximation. In RBC, temperature is the only thermodynamic quantity and throughout this work we will refer to the choice of fixing the flux at the bottom and temperature at the top as "FT" boundary conditions. We will refer to the common choice of fixing temperature at both boundaries as "TT" boundaries, and fixing the flux at both boundaries as "FF" boundaries 1 . It is generally assumed that FT and FF boundaries should fundamentally behave similarly [6], and the behavior of FF boundaries is well-known [14,15]. FF and TT boundaries exhibit the same scaling of convective heat transport (quantified by the Nusselt number, Nu) as a function of increased convective driving (quantified by the Rayleigh number, Ra) [15]. It is natural to assume that FT simulations should follow the same Nu vs. Ra scaling laws as FF and TT boundary conditions. However, FT boundaries introduce complexities into the convective solution which neither FF nor TT boundaries are exposed to. First, the evolved mean temperature of a simulation with FT boundaries differs from the initial mean temperature, and therefore the thermal reservoir of the convective system must evolve ("relax") over time [16]. Second, FT boundary conditions are fundamentally asymmetric, and it is unclear if these asymmetries affect the evolved convective solution.
In this paper, we investigate the thermal relaxation of and the asymmetries in RBC with FT boundary conditions. We also compare relaxed FT solutions to TT solutions. Our results suggest that the thermal relaxation of FT systems is very long compared to TT systems, in which it is nearly instantaneous. The thermal relaxation of FT simulations is analogous to a sweep through parameter space in which dynamics are sampled over a range of values of Ra. We find that this long thermal relaxation can be bypassed by using the results of TT simulations as initial conditions for FT simulations. Finally, FT boundaries create some asymmetries in the convective flows, but these asymmetries do not meaningfully change the mean convective state compared to TT simulations.
We present these findings as follows. In section II, we describe our simulation setup and numerical methods. In section III, we first describe our findings with respect to the time evolution of FT systems, and then describe the asymmetries in these systems. In section IV, we show that these findings carry over to a more complex system (rotating Rayleigh-Bénard convection) with some interesting implications. Finally, in section V, we summarize our findings and briefly describe the implications of this work for the field of astrophysical convection.

II. SIMULATION DETAILS
We study incompressible RBC under a freefall nondimensionalization; for details of this nondimensionalization, we refer readers to our previous work [16]. In section IV, we study convection in the presence of vertical global rotation [17], and include the Coriolis term in the momentum equation for generality. The Boussinesq equations of motion are where u = (u, v, w) is the velocity, T = T 0 (z) + T 1 (x, y, z, t) is the temperature (where T 0 is the initial profile and T 1 are the fluctuations around that profile), is the reduced kinematic pressure [16] which enforces the incompressibility constraint, and ω = ∇ × u is the vorticity. The dimensionless control parameters are the Rayleigh (Ra), Prandtl (Pr), and Ekman (Ek) numbers, defined respectively as where u ff is the freefall velocity, g is the gravity, α is the coefficient of thermal expansion, L z is the domain depth, ν and κ are respectively the viscous and thermal diffusivity, Ω is the global rotation frequency, and ∆ is the nondimensional temperature scale (defined below). These parameters set the freefall Reynolds (Re ff ) and Péclet (Pe ff ) numbers, and throughout this work we hold Pr = 1 so that Re ff = Pe ff . In non rotating RBC (section III), we set Ek = ∞.
The extent of our numerical domain vertically is z = [−0.5, 0.5] and horizontally is x, y = [−Γ/2, Γ/2], where Γ is the aspect ratio. The initial temperature profile, T 0 (z) = 0.5 − z, is unstable and linearly decreases from a value of 1 to 0 across the domain. The temperature scale, ∆, is set by the temperature jump across the domain (∆ = ∆T 0 = T 0 (z = 0.5)−T 0 (z = −0.5)) for TT boundaries or by the temperature gradient length scale (∆ = L z ∂ z T 0 ) for FT boundaries. We respectively define a temperature (Ra ∆T ) and a flux (Ra ∂zT ) Rayleigh number for these cases, We impose FT and TT boundary conditions respectively as In section III, we study two-dimensional (2D) convection where ∂ y = v = 0. For comparison with the literature, we specify Γ = 2 and these simulations employ no-slip, impenetrable boundaries, For this choice of boundary conditions, the critical value of the Rayleigh number is Ra ∂zT = 1296 for FT boundaries and Ra ∆T = 1708 for TT boundaries [6]. Due to this difference in onset, it is reasonable to expect important differences between FT and TT solutions at low supercriticalities. However, for the supercriticalities of O(10 5+ ) studied in Sec. III, we do not expect this difference in linear stability to be very important. The rotating cases in section IV employ stress-free, impenetrable boundaries, We follow previous work [18] and study three-dimensional (3D) tall, skinny boxes with Γ = 10λ c (Ek), where λ c (Ek) is the wavelength of convective onset at the specified value of Ek. For the cases studied here at Ek = 10 −6 , and for TT boundaries, λ c (10 −6 ) ≈ 4.81 × 10 −2 and the critical Rayleigh number is Ra ∆T ≈ 9.2 × 10 8 . We utilize the Dedalus 2 pseudospectral framework [19,20] to evolve Eqs. (1-3) forward in time. Our 2D simulations use an implicit-explicit (IMEX), third-order, four-stage Runge-Kutta timestepping scheme RK443; our 3D rotating simulations use the second-order, two-stage Runge-Kutta scheme RK222 [21]. The code used to run simulations and to create the figures in this work are available publicly online in a repository of supplemental materials [22] 3 . Variables are time-evolved on a dealiased Chebyshev (vertical) and Fourier (horizontal, periodic) domain in which the physical grid dimensions are 3/2 the size of the coefficient grid. We fill T 1 with random white noise whose magnitude is 10 −6 /Pe ff , and which is vertically tapered to zero at the boundaries. We filter this noise spectrum in coefficient space, such that only the lower 25% of the coefficients have power; this low-pass filter is used to avoid populating the highest wavenumbers with noise in order to improve the stability of our spectral timestepping methods.

A. Nondimensional Output Quantities
Throughout this work we will measure and report the evolved value of the Nusselt number (Nu). We define and measure Nu instantaneously as where represent a volume average ( A ≡ A dx dz/Γ in 2D and A ≡ A dx dy dz/Γ 2 in 3D for some quantity A), and ∆T = ∂ z T is the (negative) temperature difference between the top and bottom plate. In an evolved, statistically stationary state, Nu = 1 + Pe ff wT = ∂ z T | z={−0.5,0.5} when TT boundaries are employed, and Nu = (−∆T ) −1 when FT boundaries and a flux nondimensionalization are employed. This implies that Nu is the conversion between a temperature and flux nondimensionalization such that the equilibrated state of any convective solution is characterized by Ra ∆T and Ra ∂zT according to Ra ∂zT = Ra ∆T Nu. (11) Throughout this work, we will also measure the evolved Péclet number (Pe) and in section IV we will measure the Rossby number (Ro). These nondimensional quantities are defined as where |A| represents the magnitude of the vector A.

A. Time Evolution
In Fig. 1, we compare the time evolution of the temperature field of an FT simulation with Ra ∂zT = 4.83 ×10 10 to two TT simulations (with Ra ∆T = 10 10 and Ra ∆T = 10 9 , respectively). As shown in the top four panels, we see the expected convective roll solution in both TT simulations (top row) and at early and late times in the FT simulation (bottom row). Interestingly, we find highly asymmetrical dynamics at early times in the FT simulation (bottom left), in which the temperature anomaly in the cold plume is much greater than in the warm plume. This excess cold material slowly fills the domain and mixes, reducing the temperature difference between the top and bottom plates To first order, both cases have similar flow structures: a large convective cell and plumes which break apart into small turbulent eddies. However, in the FT case, the temperature anomaly of the cold plume is much larger than the hot plume, which does not appear on this color scale. (Right two panels) Dynamics in a TT case at Ra∆T = 10 9 and in the relaxed state of the previously pictured FT simulation with Ra∆T ≈ 10 9 . The relaxed FT simulation is visually indistinguishable from its comparable TT simulation. (Bottom four panels) Probability distribution functions (PDFs) of the full temperature field in each of the four dynamical panels pictured above. PDFs sample the dynamics in each simulation once every simulation freefall time unit over the course of 500 freefall time units. The black vertical line shows the median value, and the grey outline shows the 68% confidence interval, or where the cumulative distribution function (CDF)'s value ranges from 0.16 to 0.84. (Top row) In both TT simulations, the temperature field has a mean value at T = 0.5 and a symmetric distribution around that peak with maxima at the fixed values of the boundaries. (Bottom left) At early times in the FT simulation, the modal value of the PDF constantly moves left (towards the cold fixed-temperature boundary). (Bottom right) At late times, the temperature PDF from the cold fixed-temperature value (on the left) to the mode is indistinguishable from the TT PDF, but from the mode to the fixed-flux boundary there is a large tail characterized by low-probability, hot elements. from ∆T = −1 to ∆T = −Nu −1 in the relaxed state. In this relaxed state, the supply of warm fluid from the bottom plume and cold fluid from the top plume come into balance, and the FT dynamics (bottom right) are indistinguishable from TT dynamics (top right).
In the bottom four panels, we examine these temperature fields statistically by displaying their probability distribution functions (PDFs). To create these PDFs, we sample the full simulation temperature field once every simulation freefall time unit over the span of 500 time units. We interpolate the (unevenly spaced) vertical Chebyshev grid points onto an evenly spaced grid before histogramming the flow values into 200 bins and creating the PDFs.
We find that this statistical analysis of the simulations tells the same story as the dynamical images shown above. The temperature field in both of the TT simulations (top row) is dominated by the modal temperature of 0.5 in the bulk; a smaller fraction of the domain is filled with equal portions of hotter/colder material (mostly contained in the plumes), and the temperature field is rigidly bounded by the fixed-temperature boundary values. The story is more complex for the FT simulation. At early times (lower left), the FT simulation is characterized by two features: an extreme tail (to the left) that characterizes the cold plume at the upper boundary, and a migrating modal temperature that shifts from the right (hotter) to the left (cooler) as cold material mixes in the interior. At late times (lower right), the FT simulation's PDF is indistinguishable from the TT PDF between the cold (left), fixed-temperature boundary and the modal value. From the modal value towards warmer temperatures, we find that the fixed-flux (lower boundary) is capable of producing more extreme temperature events and results in a more extended PDF tail. This long tail is explored further in section III D.
In the left panels of Fig. 2, we examine the time evolution of scalar quantities from the FT simulation shown in Fig. 1 (orange lines) and compare it to the TT simulation with Ra ∆T = 10 9 (purples lines). Simulation time is shown in nondimensional freefall units on the x-axis; the latest time displayed for each simulation, t final , is subtracted for direct comparison of the relaxed states. Traces of Ra ∆T and Ra ∂zT are shown in the top-left panel. In the FT simulation, Ra ∆T relaxes to its final value over thousands of simulation time units, and this final value is the input value of the equivalent TT case. In comparison, Ra ∂zT for the TT case instantaneously reaches its final value, which is the input value for the FT simulation. This discrepancy in evolution timescales, where TT simulations evolve quickly and FT simulations evolve slowly, is also seen in the equilibration of Nu (middle panel) and Pe (bottom panel).
The right panels of Fig. 2 show that the relaxation of Ra ∆T in FT simulations is akin to a sweep through Ra ∆T parameter space. The orange (Ra ∂zT = 4.83 × 10 10 , as on the left) and yellow (Ra ∂zT = 2.61 × 10 9 ) lines show the evolution of FT simulations, and the arrows give the sense of time in the simulations. For comparison, we plot results from TT simulations (purple circles) and the reported results of ref. [24] (black crosses). The purple circles filled with orange and yellow circles are comparison TT simulations for the relaxed states of the FT simulations. The top-right panel is a scaling plot for Nu vs. Ra ∆T compensated by the best fit reported in ref. [15]. The bottom-right panel is a scaling plot of Pe vs Ra ∆T compensated by the expected scaling [2]. We find that FT simulations carry marginally more flux (higher Nu) and are more turbulent (higher Pe) than comparable TT simulations as they relax through this parameter space. These heightened values of Nu and Pe suggest that the dynamics do not immediately "forget" the higher-Ra ∆T state that they recently timestepped through on their way to achieving thermal relaxation.
FT simulation dynamics evolve slowly during thermal relaxation, and these images, PDFs, and traces demonstrate the importance of waiting for thermal relaxation to be achieved when conducting an FT simulation.

B. Evolved Structure
In Fig. 3, we compare the time-and horizontally-averaged profiles of the temperature and fluxes in the evolved FT and TT cases presented in Fig. 2. Time averages are taken over 500 nondimensional freefall time units for the TT case and over 2000 nondimensional freefall time units for the FT case, sampled once every 0.1 time units (the difference in averaging window is explained in Sec. III C 1). In the three left panels, we display profiles of (top) the mean temperature, (middle) the mean temperature in upflows (solid) and downflows (dashed), and (bottom) the convective enthalpy flux (F enth = wT , solid) and the conductive flux (F cond = −Pe −1 ff ∇T , dashed). Most of the interesting structure is in the boundary layers, located between the sides of the plots and the thin vertical black lines. Zoomed in plots of the bottom and top boundary layers are respectively shown in the middle and right columns. Inset panels show the percentage difference between the FT and TT solutions. In the flux panels (bottom row), we do not plot the percentage difference in the conductive flux, as this quantity is undefined in the bulk of the interior where that flux is zero. The conductive flux of the two cases agrees to within a few % in the boundary layers, and the largest deviations away from zero in the interior are O(0.01) in the plotted units.
We find good (∼ 1%) agreement between the FT and TT temperature profiles and enthalpy fluxes throughout the full depth of the domain, with slightly larger differences near the bottom boundary where the boundary conditions differ. When we split the temperature profile into upflows and downflows, we find that FT upflows/downflows are slightly warmer/cooler than their TT counterparts at the hot, bottom boundary. These differences are interesting, and are explored further in section III D and Fig. 5, but vanish in the interior and do not likely affect the convective dynamics appreciably. Achieving thermal relaxation in pure FT simulations is computationally costly for two reasons: (1) the turbulent dynamics at the large initial Ra ∆T require more spectral modes to resolve than the equilibrated state (compare the left and right dynamics in Fig. 1), and (2) thousands of freefall times must pass during relaxation (see Fig. 2). For example, for the cases displayed in the left panels of Fig. 2 with a modest Ra ∆T = 10 9 , the shown evolution of 10 4 time units of the FT simulation cost ∼ 4.5 × 10 5 cpu-hours, while the TT equivalent case cost only 5.6 × 10 4 cpu-hoursnearly an order of magnitude difference. In the previous sections, we have shown that relaxed FT and TT simulations are very similar. It should therefore be possible to use results from a TT simulation to quickly reach the relaxed state of a comparable FT simulation, saving up to an order of magnitude in computational cost.
We now briefly describe a procedure which takes advantage of the fast evolution of TT simulations to achieve equilibrated FT simulations. In short, the evolved flow fields in the TT simulation are properly re-nondimensionalized and used as initial conditions for an FT simulation. Throughout this work, we will refer to simulations conducted this way as "TT-to-FT" simulations. To achieve this, we perform these steps: 1. Run a TT simulation to its statistically stationary state (∼ 100+ freefall time units). Measure Nu in that state.
2. Re-nondimensionalize from ∆ = ∆T → ∂ z T , and from Ra ∆T →Ra ∂zT , as in Eqn. 11. In our freefall nondimensionalization, this means setting the velocities in the FT simulation to u FT = u TT / √ Nu, and setting the temperature field to T FT = T TT /Nu.  We show this procedure in practice in Fig. 4. In the left three panels, we display the temporal behavior of (top) Ra, (middle) the flux at the bottom boundary, and (bottom) the temperature difference between the top and bottom boundaries. We take the full evolution of the Ra ∆T = 10 9 TT simulation shown in Fig. 2, then change its boundary conditions to FT at Ra ∂zT = 4.83 × 10 10 and restart the simulation using the above procedure. The change from TT to FT boundaries occurs at the time denoted by the thin vertical line. Unlike in the FT case displayed in Fig. 2, there is no thermal rundown in the FT state, due to the rapid relaxation achieved during the TT portion of the simulation.
In the right four panels of Fig. 4, we compare PDFs of flow fields in this TT-to-FT simulation and the comparable FT simulation which we timestepped through thermal relaxation. Shown are PDFs of the temperature field (upper left), enstrophy (upper right), convective flux (lower left), and vertical velocity (lower right). In Table I, we display the first four moments of each of these distributions, where A is a flow quantity, P (A) is the PDF of A, µ is the mean, σ is the standard deviation, ∆A is the spacing between the discrete PDF bins, and i is the index of the bin. The PDFs are qualitatively similar, and there is generally good agreement between the moments of the PDFs. The remaining discrepancies between the PDFs are small and seem to primarily be due to the randomness of the flows in the time windows over which we sampled the simulations.
We note briefly that this TT-to-FT mechanism is only one of many ways of accelerating the thermal relaxation of an FT simulation. We discuss other mechanisms, and explore one in detail, in our previous work [16]. We note however that the TT-to-FT setup described here is likely the least complicated mechanism for achieving rapid relaxation in a simplified RBC setup that we are aware of. The successful degree with which this mechanism reproduces the evolved dynamics suggests that thermal relaxation occurs in two parts: 1. Changes to the simulation energy reservoir, and 2. Restratification of the experiment.
The thermal energy reservoir of TT simulations does not change between the initial and final state. The rapid relaxation of TT simulations therefore suggests that experimental restratification occurs rapidly in RBC. The long rundown of FT experiments on display in Fig. 2 is entirely due to the energy reservoir (the mean temperature) drifting over time. Put differently, the classic RBC setup for T 0 (z) is a bad choice of initial conditions for FT boundaries, and TT-to-FT simulations use TT dynamics to choose a more ideal set of initial conditions.

Changing Timescales
One surprising result of our FT simulations is that the nondimensional dynamical freefall timescale is a poor description of the evolved freefall timescale. As noted above in our TT-to-FT procedure, the velocities in an FT simulation are smaller than the velocities in a TT simulation by a factor of √ Nu. As a result, every simulation freefall time unit in a TT simulation samples a factor of √ Nu more dynamics than a nondimensional freefall time in an FT simulation. This is on display in the left panels of Fig 4, where the flow variations occur at a higher temporal frequency in the TT portion of the simulation than in the FT portion.
These results imply that the true nondimensional freefall time of an FT simulation is τ ff = √ Nu in simulation freefall time units. However, we note that Ra ∂zT is the input value of Ra for an FT simulations, and that Ra ∂zT /Ra ∆T = Nu. The thermal diffusion timescale, τ th , which scales like τ th ∝ Ra input , is also larger in an FT simulation than in a comparable TT simulation by a factor of √ Nu. One convective freefall timescale therefore occurs over the same fraction of a diffusion timescale regardless of boundary conditions or temperature nondimensionalization, τ ff /τ th = f (Nu). These findings suggest that comparing FT and TT dynamics may be more straightforward under a thermal diffusion nondimensionalization [6] than the freefall nondimensionalization we have used throughout this work.

D. Asymmetries induced by mixed boundaries
We now study in more detail the asymmetries introduced into a solution by FT boundaries. We run a TT and TT-to-FT simulation at Ra ∆T = 10 10 and Ra ∂zT = 9.51 × 10 11 , respectively. In Fig. 5, we examine the dynamical nature of the asymmetries which FT boundaries introduce into the simulation near the fixed-flux boundary. In the left panel, we plot PDFs of the temperature fields in comparable TT and FT simulations. These PDFs agree remarkably well near the mean and for cold temperatures (near the fixed temperature boundary), but diverge in the tail of the PDF for hot temperatures where T /∆T 0.8, where the boundary conditions differ. Interestingly, there are no temperature fluctuations which exceed the specified boundary values in the convective domain for TT simulations. However, the FT PDF has a much longer tail and the FT solution achieves fluid parcels which are hotter than the average bottom boundary value by more than 50%. In order to understand how this is possible, we examine a snapshot of the FT simulation's temperature anomaly in the middle panel. We have outlined a portion of a cold plume near the upper (fixed-temperature) boundary and a portion of a hot plume near the lower (fixed-flux) boundary, and these regions are magnified in the rightmost panels. The TT upper boundary suppresses temperature anomaly at the upper boundary and regulates the temperature minima which can be achieved. The fixed-flux lower boundary does no such suppression and allows for extreme temperature values to be achieved in the plume-launching area, thus allowing for the asymmetry in the tails of the temperature PDF.
We note briefly that these asymmetries do not seem to affect mean or volume-averaged quantities in these simulations appreciably (see the agreement between FT and TT in Figs. 2&3). However, the fact that fixed-flux boundaries produce a wider temperature distribution with more extreme values may be important in some astrophysical studies. We explore this further in the discussion in section V.

IV. ROTATING RAYLEIGH-BÉNARD CONVECTION
We now extend our study to a more complicated experiment: 3D rotating RBC with an Ekman number of 10 −6 . We study a TT case at Ra ∆T = 2.75 × 10 9 , and an FT case at Ra ∂zT = 2.1 × 10 10 (the supercriticality of the TT case is ∼ 3). These simulations employ stress free boundary conditions which allow for the generation of mean flows such as large scale vortices (LSV) [25].
In the left three panels of Fig. 6, we compare the time evolution of the FT and TT cases. The top left panel shows the evolution of Ra ∂zT and Ra ∆T . Even in the presence of strong rotation, the TT immediately equilibrates, but the FT case takes thousands of freefall times to achieve thermal relaxation. In the middle panel, we show the evolution of Ro; the evolved flows in both simulations exhibit rotationally constrained dynamics with Ro ≈ 0.1, but the flows in the FT simulation relax to this state from an initially unconstrained state (Ro ≈ 1). This implies that the thermal relaxation process can walk through the parameter space of flow balances (e.g., the balance between Inertial and Coriolis forces) in addition to the Ra ∆T parameter space. In the bottom panel, we display the evolution of Pe over time. Strangely, the peak value of Pe occurs a few hundred freefall times after the convective transient. After achieving this peak value, Pe monotonically decreases toward its relaxed state.
In the upper right panel of Fig. 6, we plot Nu vs. Ra for rotating simulations. Select TT cases are plotted as cyan circles with purple outlines (where the cyan color denotes the value of Ek = 10 −6 according to the color bar). The evolution of the FT case in the left panels is shown as a thick orange line with a cyan interior and the black arrows show the direction of time. The TT case that corresponds to the FT case is a purple star with a cyan interior. We have additionally included some literature data from numerical simulations (circles) and experiments (diamonds) as reported in the appendix tables of ref. [26]. These experiments were conducted in a cylindrical geometry at a different Pr, and are not meant to be one-to-one-comparable, but are meant to guide the eye to the nature of the parameter space of rotating convection. The solid black line is the best-fit line for rotationally unconstrained simulations with Ra ≥ 10 10 from ref. [26]. As expected, the scaling of Nu vs. Ra is steep in the rotationally constrained regime [3,27], which these simulations trace through. As in Fig. 2, the FT values of Nu are once again somewhat elevated above the comparable TT simulations.
In the bottom right three panels of Fig. 6, we plot the vertically integrated vertical vorticity in the simulation at three different times. In the left panel, a dominant LSV which is aligned with the global rotation dominates the simulation at early times. Over thousands of freefall times, this LSV evolves into a long-lived vortex pair, displayed in the middle panel. Finally, in the evolved state, this vortex pair solution begins to oscillate with domain-wide jets, such as those displayed in the right panel. We find that the TT solution shows this oscillatory behavior between vortex pairs and jets immediately and throughout the full 5000 freefall timescales of evolution that we simulated.
We suspect that the strange behavior of Pe in the bottom left panel can be explained by the evolution of the dominant flow structures over time. At early times, the initially large value of Ra ∆T in the FT case drives the displayed dominant LSV. This powerful driving injects energy into the LSV, causing Pe to grow. As Ra ∆T and convective driving decrease over time, the LSV saturates and then starts to wind down, leading to the "bump" in the Pe trace.
We once again find it important to briefly note the difference in computational cost between the FT and TT simulations conducted here. The TT simulation shown in the left panels of Fig. 6 only cost 2.6 × 10 4 cpu-hours. By comparison, the cost of the FT simulation shown in the same panels was two orders of magnitude larger -2.3 × 10 6 cpu-hours. The TT simulation's coefficient resolution was 128 3 . The FT simulation's initial resolution required to resolve the convective transient was 512 × 384 2 coefficients. We reduced the resolution to 256 × 384 2 after 100 freefall times, and then later to 128 × 384 2 after ∼ 3.3 × 10 3 freefall times. At each of these times, we found that lowering the horizontal coefficient resolution of the simulation did not reproduce the simulation solution with fidelity. This suggests that small scale turbulent velocity structures-which are injected by the vigorous transient and perhaps associated with the LSV-are long lived throughout the thermal evolution of the simulation.

V. CONCLUSIONS & DISCUSSION
In short, we find that FT simulations experience a long thermal relaxation which is not experienced by TT simulations and, to first order, FT boundaries do not introduce important asymmetries into the solution.  [18]. In the converged state, we see oscillatory behavior between this vortex pair behavior and jets (right panel). The TT case exhibits the oscillatory behavior between vortex pairs and jets throughout its whole evolution. The three vertical black lines in the left panels signify the times at which these snapshots are taken.
In this paper, we have studied the time evolution of Rayleigh-Bénard convection (RBC) under two different formulations of the thermal boundary conditions: "FT" boundaries, where the flux is fixed at the bottom and temperature is fixed at the top, and "TT" boundaries, where temperature is fixed at the top and bottom. Through studying this relaxation and the relaxed states of both simulations, we come to the following conclusions: 1. Thermal relaxation in RBC has two components: (a) changes in the energy reservoir and (b) changes in the stratification. We find that the long relaxation of FT simulations is due to changes in the energy reservoir; this reservoir is roughly constant in TT simulations due to the lack of evolution of the average domain temperature. The rapid evolution of our TT simulations suggests that RBC restratifies itself instantaneously.
2. Dynamical measurements taken during the thermal relaxation of an FT simulation may be misleading. Dynamics during the relaxation are more turbulent than in the evolved state, and exhibit evolving flow balances in the equation of motion (as quantified by e.g., the Rossby number).
4. Great computational expense achieving thermal relaxation in an FT simulation can be avoided by using the evolved state of a TT simulation as a "better" set of initial conditions for an FT simulation.
5. Despite minor asymmetries near the fixed-flux boundary, we find no meaningful difference between the mean state of FT and TT simulations.
We now describe some lessons that should be applied from this work to astrophysical convection, and comment on some open areas of research. Throughout this work, we have made the assumption that convection is only "interesting" in its final, fully equilibrated state. In nature, convection is not always in an equilibrium state. For example, in the late stages of the lifetimes of stars, some core burning regions have sufficiently short lifetimes that they likely do not come into thermal relaxation [28,29]. The use of FT boundaries or initial conditions that we have here considered to be "bad" choices may help in understanding these transient lifetime stages. However, for most convective studies where the lifetime of the natural convective system is much larger than its Kelvin-Helmholtz timescale, it is essential to study relaxed convection, and our results point towards the importance of either choosing good initial conditions (TT or TT-to-FT simulations) or running simulations to thermal relaxation.
One question which our study of RBC is not able to address is: how long does it take for a complex convective system to restratify? Our fully convective domains restratified instantaneously, but it is likely that mixed convectiveand-stably-stratified domains [11,13,30,31] should have regions that are not turbulently mixed by convection which could also have long relaxation timescales. It would be extremely helpful for future studies to examine relaxational timescales in systems where the energy reservoir is fixed, but where convection does not effectively mix the whole domain. Fortunately, clever techniques (e.g., as we explored in ref. [16]) can likely be used to rapidly restratify atmospheres in such simulations.
RBC is fundamentally symmetrical, but many natural convective processes occur in density-stratified domains in which the symmetries of the problem are fundamentally broken. In the present study, we observed that flux boundaries produce more extreme thermodynamic events than temperature boundaries. In studies of overshooting convection, it is possible that plumes produced by a flux boundary layer could launch further into a stable layer than plumes produced by a temperature boundary. Some authors have aimed to quantify the nature of overshooting plumes from a convective region into a stable region [11,31], and it is unclear if different choices of boundary conditions could change the observed distribution of overshooting plumes observed there.
Some of the most complex astrophysical convection experiments aim to understand self-consistently evolving magnetic dynamos in rotating, spherical, magnetohydrodynamical domains [32][33][34][35]. These dynamo simulations involve large numbers of timesteps through many freefall timescales in order to study the generation and evolution of magnetic fields and mean flows. We found in our FT rotating simulation that the unrelaxed state generated a mean flow (a LSV, Fig. 6) that was much more intense and large-scale than the eventual flows that developed in the relaxed state. If we had terminated our FT rotating simulation too early, we would not have seen the eventual destruction of this LSV or the later oscillatory behavior between jets and vortex pairs. Many dynamo simulations are performed in highly turbulent regimes at the cutting-edge of what is achievable using modern computational resources. As a result, timestepping through thousands of freefall timescales is not possible in these simulations. It is therefore crucial that dynamo simulations be set up in such a manner as to avoid large changes to the system's energy reservoir such as those that we observed and studied here. Some authors who study astrophysical convection [10,[35][36][37] employ FF boundary conditions, and our results here suggest that such a choice may be ideal in those complex simulations.
In conclusion, we note that our results here should provide astrophysical convection simulations with reason for optimism. Some problems that we encounter (e.g., long thermal rundown in FT simulations) can be completely avoided through a careful understanding of the numerical system being solved. Input quantities are the boundary conditions (BCs), Rayleigh number (Ra), coefficient resolution (nx×ny×nz, or horizontal × vertical), and the total simulation run time in freefall units t simulation and in cpu-hours. Output quantities are the Nusselt (Nu), Péclet (Pe), and Rossby (Ro) numbers. Reported values of Nu, Re, and Ro are the sample mean over the last 500 freefall time units. Reported uncertainties are the standard deviation of the sample mean; when the uncertainty is not reported, it is smaller than the number of reported digits. The "Nu comp" values are comparison Nu values reported in ref. [24]. Resolutions marked by a * show the initial, highest resolution utilized in the simulation. The 2D FT Ra = 4.83 × 10 10 simulation's resolution was changed to 1024 × 2048 about 500 freefall time units after transient. The rotating FT case's resolution was reduced to 384 2 × 256 about one hundred freefall time units after transient, and was further reduced to 384 2 × 128 about 3.3 × 10 3 freefall times after transient.