Revisiting the Sun's Strong Differential Rotation along Radial Lines

Current state-of-the-art models of the solar convection zone consist of solutions to the Navier-Stokes equations in rotating, 3D spherical shells. Such models are highly sensitive to the choice of boundary conditions. Here we present two suites of simulations differing only in their outer thermal boundary condition, which is either one of fixed-entropy or fixed-entropy-gradient. We find that the resulting differential rotation is markedly different between the two sets. The fixed-entropy-gradient simulations have strong differential rotation contrast and isocontours tilted along radial lines (in good agreement with the Sun's interior rotation revealed by helioseismology), while the fixed-entropy simulations have weaker contrast and contours tilted in the opposite sense. We examine in detail the force balances leading to the different rotation profiles and find that the poleward transport of heat by Busse columns plays a key role. We conclude that the Sun's strong differential rotation along radial lines may result from the solar emissivity being invariant with latitude (which is similar to the fixed-entropy-gradient condition in our models) and the poleward transport of heat by Busse columns. In future work on convection in the solar context, we strongly advise modelers to use a fixed-gradient outer boundary condition.


INTRODUCTION
Helioseismology has revealed in detail the internal rotation profile of the solar convection zone (CZ; e.g., Thompson et al. 2003;Howe et al. 2005), as shown in Figure 1. The most notable properties of the rotation rate are that the equator rotates significantly faster than the high-latitude regions and that the isorotation contours are tilted significantly with respect to the rotation axis, falling largely along radial lines. Furthermore, there are two shear layers at the top and bottom of the CZ: at the top, the contours bend toward the equator in a region known as the near-surface shear layer (NSSL), and at the bottom, the differential rotation in the CZ transitions to solid-body rotation, over a narrow boundary layer called the tachocline. Prior to helioseismic probing, most theoreticians had assumed that the differential rotation that is observed directly at the surface would imprint into the interior along isosurfaces parallel to the rotation axis, hence satisfying the Taylor-Proudman theorem. The helioseismic observations have clearly demonstrated that this theoretical supposition was wrong.
For the last several decades, global, 3D supercomputer simulations of hydrodynamic convection in rotating spherical shells have succeeded in achieving rotation profiles that are fast at the equator and slow at the poles. However, simulations generally have a weaker overall differential rotation contrast than that of the Sun. If the contrast is defined to be the difference in rotation rate between the equator and 60 • latitude, expressed as a percentage of the "frame" rotation rate, then for the Sun, this magnitude is ∼20%. Most simulations, on the other hand, have magnitudes of ∼10%, although there are some notable exceptions (e.g., Brown et al. 2010;Matilsky & Toomre 2020). To date, however, 1 loren.matilsky@colorado.edu  [1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009]. The rotation rate has been obtained using a regularized least-squares (RLS) inversion, which is sensitive only to the equatorially symmetric part of the rotation. The dashed lines are at a 25 • angle to the rotation axis and align with the rotation contours at mid-latitudes. Image credit: Howe et al. (2005), extended to include GONG data until 2009. there is no systematic physical explanation for why these exceptional simulations have high rotation contrast.
Simulations have also struggled to achieve rotation contours in the bulk of the CZ that are significantly tilted from the axis, as seen in Figure 1. With some exceptions (e.g., Miesch et al. 2006;Elliot et al. 2000), the simulations generally have rotation contours aligned with the rotation axis. In the case of Miesch et al. (2006), the contours were tilted systematically by imposing a latitudinal entropy gradient at the base of the CZ. By modifying the magnitude of this entropy gradient, Miesch et al. (2006) controlled the magnitude of the contour tilt. Imposing this boundary condition was inspired by thermal wind balance in the tachocline below imprinting on the CZ.
In this paper, we explore the role that the outer thermal boundary condition plays in the resulting differential rotation. We consider the two most commonly used options: fixed entropy (FE), in which the entropy at the outer boundary is fixed to a constant value, and fixedflux (FF), in which the radial entropy gradient at the outer boundary is fixed, therefore implying an outward conductive flux that is independent of latitude (see Anders et al. 2020 for a detailed analysis of these boundary conditions for simulations of Rayleigh-Bénard convection). The FF condition is more appropriate for the Sun, since the radiant flux from the solar photosphere does not appear to vary substantially with latitude. Rast et al. (2008) analyzed full-disk images from the Precision Solar Photometric Telescope at the Mauna Lea Solar Observatory and calculated non-magnetic contributions to solar photospheric intensity. In both continuum and Ca II K intensity distributions, only a ∼0.1-0.2% variation was observed, corresponding to a solar pole at most ∼2.5 K warmer (in terms of effective temperature) than the equator.
We do not address the dynamics of the near-surface nor tachocline shear layers in this work. In particular, our models have an impenetrable bottom boundary that does not allow for the convective overshoot of downflow plumes into the stable region that may play a role in the origin of the tachocline of shear. The dynamical maintenance of the NSSL is still an open question, as discussed in Hotta et al. (2015) and Matilsky et al. (2019), and models tend to only display signs of nearsurface shear if they have high density contrast across the CZ. To avoid high computational cost, the models in this work only have modest density contrast and do not attempt to achieve a near-surface shear layer. Solar-like differential rotation (fast equator and slow poles) in spherical-shell convection is thought to be due to Busse columns (also called "Taylor columns" and "banana cells") transporting angular momentum outward (e.g., Busse 2002;Jones et al. 2011). In this work, we show that Busse columns also transport heat poleward and equivalently drive a solar-like differential rotation through a thermal wind, or latitudinal pressure gradients. We find that the thermal wind in our FF cases drives stronger differential rotation magnitudes and achieves more significant tilt in the rotation contours than the corresponding FE cases. Elliot et al. (2000) noted this effect for one simulation, but did not explore the underlying mechanism. Here we present a detailed analysis of the manner by which the outer thermal boundary condition modifies the resulting differential rotation in our simulation suite.
In Section 2, we describe the parameter space explored by our simulation suite, as well as the mathematical details of the FF and FE boundary conditions. In Section 3, we describe the basic results of our experiment, focusing on the achieved differential rotation. In Section 4, we quantify the force balance achieved in our models, which, for the radial and latitudinal directions, consists of a thermal wind in spherical geometry. In Section 5, we examine the latitudinal transport of energy by Busse columns that is responsible for the thermal wind. In Section 6, we discuss how the effects of the thermal wind are modified by the outer thermal boundary condition. In Section 7, we discuss our simulation results in the context of the Sun.

NUMERICAL EXPERIMENT
We consider time-dependent, 3D simulations of a rotating, stratified spherical shell of fluid representative of the solar CZ. We use the open-source code Rayleigh 0.9.1 (Featherstone & Hindman 2016a;Matsui et al. 2016;Featherstone 2018), which solves the equations of hydrodynamics in spherical geometry. Our domain is a spherical shell with inner radius r i and outer radius r o . We describe this domain using the standard spherical coordinates (r, θ, φ) and corresponding unit vectors (ê r ,ê θ , e φ ). When convenient, we also use cylindrical coordinates (λ, φ, z) = (r sin θ, φ, r cos θ) and unit vectors (ê λ , e φ ,ê z ), where λ is the cylindrical radius and z is the axial coordinate perpendicular to the equatorial plane.
The thermodynamic reference state is chosen to be temporally steady and spherically symmetric with adiabatic stratification (see Jones et al. 2011 for a complete description). The density varies by a factor of ∼20 (three density scale heights) across the layer. We denote the pressure, density, temperature, and entropy by P , ρ, T , and S, respectively, using overbars to indicate the fixed reference state and the absence of overbars to indicate deviations from the reference state.
Rayleigh uses an anelastic approximation (e.g., Gough 1969;Gilman & Glatzmaier 1981;Jones et al. 2011), which removes sound waves from the system, making the maximum allowable timestep much larger since it is only limited by the flow velocity rather than the sound speed. The fluid equations are then given by (e.g., Featherstone & Hindman 2016a) and ρT ∂S ∂t is the fluid velocity in the rotating frame, c p is the constant-pressure specific heat, g is the gravitational acceleration due to a solar mass M located inside the inner radius of the spherical shell at its center, D = ρν[∇v + ∇v T − (2/3)(∇ · v)I] is the Newtonian viscous stress tensor, I is the identity tensor, ν is the kinematic viscosity, and κ is the thermal diffusivity. Because it is not computationally possible to resolve convection in the solar regime down to the turbulent microscale, ν and κ must be regarded as "eddy" diffusivities, which for simplicity we choose to be spatially constant, and equal so that the Prandtl number is unity. The internal heating function, which physically represents heating due to radiation, is chosen to have Table 1 Common input parameter values for all simulations r i 5.00×10 10 cm ro 6.59×10 10 cm cp 3.50×10 8 erg K −1 g −1 γ 1.67 ρ i 0.181 g cm −3 L 3.85 ×10 33 erg s −1 M 1.99 ×10 33 g Pr ≡ ν/κ 1.00 Note. -All values are displayed to 3 significant digits in c.g.s. units.
the fixed radial profile Q(r) = α[P (r) − P (r o )], with the normalization constant α chosen such that a solar luminosity L is forced through the domain (see Featherstone & Hindman 2016a).
The equation of state for the system is that of a perfect gas subject to small thermodynamic perturbations about the reference state: with γ being the ratio of specific heats. We adopt stress-free and impenetrable boundary conditions to conserve angular momentum and mass: where r i and r o refer to the inner and outer radii of the spherical shell, respectively. In all cases, the inner thermal boundary condition allows no flux of energy into the system through the bottom boundary: Input parameters common to all the simulations explored here are shown in Table 1.
2.1. Outer thermal boundary condition The main purpose of this work is to characterize the influence of the outer thermal boundary condition on the behavior of the resulting differential rotation. We consider models with different background rotation rates Ω 0 and diffusion values (ν = κ) and for each model analyze two sub-cases: and The solar luminosity that is injected into the system via internal heating is ultimately carried out through the outer surface via thermal conduction, which in our models arises from entropy gradients (see Equation (3)). For the fixed-entropy condition (7), the interior is initially heated (leading to S > 0 in the lower parts of the CZ) while the entropy at the outer surface is "pinned" to zero. This naturally establishes a thermal boundary layer (sharp entropy gradients ∂S/∂r < 0) just below the outer surface. The steepness of the gradient (i.e., the strength of the outward conductive loss of energy) is allowed to vary with latitude.
For the fixed-flux condition (8), the outer thermal boundary layer is present from the beginning of the simulation. The steepness of the entropy gradient (and thus the energy loss) at the outer surface is independent of latitude by construction, and is forced to have exactly the value needed to carry out a solar luminosity. The fixed-flux condition is thus more "solar-like," since in the Sun there is no observed latitudinal dependence of the emergent intensity, which is equal to the energy lost via radiative cooling at the photosphere.
The thermal conductive boundary layer stands in contrast to the real solar photosphere, in which radiative cooling removes the heat from a very thin (∼100 km) outer layer. The cooling drives very small temporal and spatial scales of motion compared to the deep interior (such as granulation and supergranulation), making its inclusion in global models difficult. Researchers have sought to address the problem in various ways. Nelson et al. (2018) implemented stochastic driving of convection by near-surface plumes designed to mimic the effects of supergranulation, finding that that the flow structures and transport properties were significantly altered in the deep CZ. Hotta et al. (2019) simulated the whole CZ with no rotation or magnetic field, coupling the global spherical shell to a Cartesian box that solved the equations of radiative transfer in the photosphere. They found that the near-surface motions had a weak influence on the deep interior. Regardless of its relevance to interior flow structures, correctly capturing the small-scale near-surface flows in global models is currently prohibitively expensive computationally. In order to explore a wider range of parameter space, we thus only consider the FE and FF conditions here.

SIMULATION RESULTS
We label simulations with a prefix that signifies the outer boundary condition ("FE" for Equation (7) and "FF" for Equation (8)), followed by the value of the diffusion constant ν = κ (in units of 10 12 cm 2 s −1 ), followed by the value of the rotation rate (in units of the solar Carrington value, Ω ≡ 2.87 × 10 −6 rad s −1 , or Ω /2π ≡ 456 nHz). For example, "case FE4-3" refers to a simulation with an FE outer boundary, for which ν = κ = 4 × 10 12 cm 2 s −1 throughout the domain, and Ω 0 = 3Ω . We note that all our models rotate at either two or three times the solar rate to avoid anti-solar differential rotation (fast poles, slow equator). The problem that highly turbulent simulations tend to be anti-solar has been called the "convective conundrum" (O'Mara et al. 2016). Anti-solar states can be avoided by either raising the rotation rate or lowering the luminosity, and we choose the former. Table 2 contains the values of the non-dimensional parameters, as well as the grid resolution, for each of the simulations considered in this work. Following the notation of Featherstone & Hindman (2016b), we parameterize the strength of the imposed driving in each simulation through a bulk "flux Rayleigh number," and the degree of rotational constraint through an Ekman number (imposed a priori), or a bulk Rossby number (calculated a posteriori), In the preceding equations, the length scale H is taken to be the shell depth r o − r i , the tildes refer to volume averages over the full spherical shell, and F refers to the the energy flux associated with the conduction and convection in equilibrium (see Featherstone & Hindman 2016b). The typical convective velocity amplitude v refers to the rms of the velocity with the longitudinally averaged part subtracted, the mean being taken over time and over the full volume of the shell (during the latter portion of run time for which there is statistical equilibrium-generally ∼3/4 of the total run time).
We quantify the magnitude of the overall differential rotation as the difference in the outer-surface rotation rate between the equator and 60 • -latitude, normalized by the frame rotation rate: where is the rotation rate of the fluid as a function of r and θ and the angular brackets denote a combined temporal and longitudinal average. From Table 2, the differential rotation fraction ∆Ω/Ω 0 is significantly greater for the FF cases-on average, by ∼40%). For comparison, the solar value of ∆Ω/Ω 0 is 0.20 (see Howe et al. 2000, Figure 1), where we have taken Ω 0 = Ω , the sidereal Carrington rate, and r o to lie just below the near-surface shear layer. Figure 2 shows how the differential rotation fraction scales with the reduced Rayleigh number, which takes into account the fact that in rotating models the critical Rayleigh number for the onset of convection is greater than in non-rotating models (Chandrasekhar 1961). The reduced Rayleigh number thus represents a better parameterization of the supercriticality of the system than simply the Rayleigh number. From Figure 2, not only do the FF cases have a stronger rotation contrast than the FE cases, but they also scale (linearly) nearly twice as fast with the reduced Rayleigh number. The rotation contrast in the FF cases does not quite reach the solar value, but progressively gets closer the higher the level of supercriticality.
To illustrate exactly where the "extra" rotation contrast in the FF cases is located, we plot the rotation rate at the outer surface for the cases rotating at 3Ω in Figure 3. Most of the additional contrast occurs at high latitudes, where the polar regions in the FF cases rotate significantly more slowly than in their FE counterparts. Additionally, the equator in the FF cases rotates slightly faster than in the FE cases. For all simulation pairs, the difference in contrast between the FE case and the FF  case is greater the smaller the value of the diffusion (or the higher the level of turbulence). Figure 4 shows contour plots of rotation rate in the upper meridional plane for the entire simulation suite. Clearly there is a striking difference between the tilts of the rotation contours in the FE and FF simulations.
In this paper, we define all rotation-contour tilt angles (or simply tilts) with respect to the rotation axis, with zero tilt corresponding to alignment of the contour with the rotation axis. We use the sign convention for tilt angle illustrated in Figure 5. Under this convention, the solar rotation contours have positive tilts at all latitudes in the bulk of the CZ (above the tachocline and below the NSSL, as shown in Figure 1). We thus define the tilt angle of a rotation contour at a given point in the meridional plane as which is consistent with the sign convention shown in Figure 5 for solar-like differential rotation, in which the contours further from the rotation axis correspond to a higher rotation rate.
Describing the solar rotation contours as "tilted along radial lines," as is often done, is technically misleading. Radial tilt implies a specific dependence of the contour  tilt angle with latitude, namely, contours that fan radially outward from the center of the Sun. In Figure 1, by contrast, the bulk-CZ tilts are roughly constant at ∼25 • for mid-latitudes, are smaller than ∼25 • for low latitudes (where radially-aligned tilts would be greater), and are greater than ∼25 • for high latitudes (where radiallyaligned tilts would be smaller). To avoid confusion, we will not use the term "radial tilt" and instead describe the rotation-contour tilt (in the Sun and in our simulations) simply as "positive" or "negative," using the sign convention illustrated in Figure 5.
In Figure 6, we show the values of the rotation-contour tilt angle at mid-depth for a subset of our models and for the Sun. The positive tilt for the FF cases is obvious, with the maximum tilt angle being about +15 • for the highest value of the diffusion (ν = 4 × 10 12 cm 2 s −1 ). This is still substantially lower than the solar value for contour tilt, which attains a maximum value of ∼25 • in   (14) shown as a function of latitude for the cases rotating at 3Ω (except case FE10-3) and for the Sun. The profiles are taken at the middle of the shell for our models and the middle of the CZ for the Sun. The northern and southern hemispheres have been averaged assuming odd symmetry for tilt angle. Dashed lines correspond to the FE cases and solid lines to the FF cases. For the solar tilt angle, we use the inversion from GONG data 1995-2004 as reported in Howe et al. (2005) and shown in Figure 1.
These results show that having an FF outer boundary condition and increasing the level of supercriticality may increase the rotation contrast toward the solar value, but it also decreases the magnitude of radial tilt. Nonetheless, the FF outer boundary is still significantly more desirable than the FE outer boundary in all cases, since the FE models with highest rotation contrast have the most negatively-tilted contours at high latitudes.

THERMAL WIND BALANCE
We find that to leading order, the longitudinally and temporally averaged force balance in the meridional directions r and θ (or λ and z) is dominated by the Coriolis, pressure, and buoyancy forces for each simulation in this work: Here, the angular brackets denote a combined temporal and longitudinal average. In the Earth's atmosphere, a "thermal wind" describes a situation in which geostrophic balance (pressure balancing the Coriolis force) holds in the horizontal directions and hydrostatic balance (pressure balancing gravity) holds in the vertical direction. Equation (15) thus represents the generalization of a thermal wind to the solar geometry, in which the vertical and horizontal extents of the flow structures are comparable (unlike in the Earth's atmosphere where the vertical extent is very thin). Furthermore, the flows in the solar geometry generally have a vertical (radial) component, unlike in a classical thermal wind for which the flows are purely horizontal.
The θ-component of Equation (15) may be rearranged to yield Ω(r, θ) ≈ Ω 0 + 1 Ω 0 ρr 2 sin 2θ which is a purely geostrophic equation, since the buoyancy force is radial. Figure 7 shows a representative example of geostrophic balance for the FF2-3/FE2-3 pair.  (16), which we denote by Ω P (rotation rate from the pressure), normalized by the frame rotation rate Ω 0 . The two hemispheres have been averaged assuming even symmetry about the equator.
Clearly Equation (16) is very well satisfied for both cases, with deviations from geostrophy being no more than 1 part in 10 3 in the bulk of the meridional plane and 1 part in 10 2 at isolated regions by the equator and pole. The same is true for all the cases considered in this work, indicating that the differential rotation profile in our simulations is almost completely determined by the pressure profile, and vice versa. The fact that the differential rotation magnitudes are ∼40% greater in the FF cases compared to the FE cases is thus a consequence of the greater latitudinal pressure gradients.
To assess why there are opposite signs of tilt for the rotation contours in the FF and FE simulations, we differentiate Equation (16) with respect to the axial coordinate z, using the radial component of Equation (15) to eliminate terms (or equivalently, take the φ-component of the curl of Equation (15)), yielding The tilt of the rotation contours is thus determined by the entropy distribution in the final thermodynamic state.
In Figure 8, we show the average profiles for entropy, pressure, and temperature in the meridional half-plane for the FF2-3/FE2-3 pair. Case FF2-3 (which is representative of all the FF cases in the simulation suite) displays a monotonically increasing entropy from equator to pole. Case FE2-3, on the other hand, has a nonmonotonic entropy profile: in most of the spherical shell (except just below the outer surface), the entropy from equator to pole increases up to ∼20 • latitude, then decreases. In both cases, the pressure and temperature deviations (normalized by the background reference state) are substantially greater (by a factor of ∼30 in the case of the pressure) than the entropy deviation. The profiles of pressure and temperature in the meridional plane thus tend to mirror one another, with high temperature regions corresponding to high pressure regions and vice versa (compare the last two columns of Figure 8).
The balance described by Equation (17) is shown for the representative simulation pair FF2-3/FE2-3 in Figure 9. There is good balance in the deep layers, although Figure 8. Temporally and longitudinally averaged entropy, pressure, and temperature deviations from the spherically symmetric mean in the meridional half-plane (averaged assuming even symmetry about the equator) for cases FF2-3 and FE2-3, normalized by the reference state profiles. The spherical mean · · · sph has been removed to show the variation from equator to pole. significant departures near the outer surface, which has been noted frequently in past work (e.g., Brun et al. 2011;Augustson et al. 2012;Hotta et al. 2015). Quantitatively, the error in Equation (17) (shown in the righthand column of Figure 9) is ∼10% in the bottom 80% of the layer and ∼50% in the upper 20% of the layer. For solar-like differential rotation (fast equator and slow poles), positively-tilted rotation contours (the FF cases) correspond to ∂Ω/∂z < 0, which arises from ∂S/∂θ < 0 at all latitudes, as in Figure 8(c). Similarly, the FE cases (which have contours tilted negatively at high-latitudes and positively at low latitudes) all have ∂S/∂θ > 0 at high latitudes and ∂S/∂θ < 0 at low latitudes, as in Figure 8(a). (16) and (17), a thermal wind in spherical geometry fundamentally consists of pressure and entropy differences in latitude. Poles that are high-pressure and high-entropy relative to lower latitudes (which we have shown lead to strong differential rotation with positively-tilted contours) are expected to be established by the preferentially poleward transport of energy. In our simulations, this transport arises from the action of the convective Busse-column rolls. These rolls manifest at convective onset as an overstable, low-frequency prograde wave (e.g., Unno et al. 1989) or, as it is called in the geophysics literature, a thermal Rossby wave. This wave consists of a series of convective rolls or Busse columns that gird the equator. Each roll is rotationally aligned and the sign of the vorticity alternates from roll to roll. Furthermore, each roll is in geostrophic balance; hence, the alternating sign of the vorticity corresponds to every other roll being a zone of high pressure, with lowpressure rolls interleaving the high-pressure rolls. Since the ends of the columns (at mid-latitudes) have neutral pressure, the pressure anomolies at the equator cause poleward axial flow in the high-pressure rolls and equatorward flow in the low-pressure rolls (e.g., Figure 1 in Gilman 1983). The resulting strong correlation between pressure and the direction of the latitudinal flow leads to a net poleward enthalpy flux through pressure work.

In light of Equations
The effect just described is easiest to illustrate for models that are barely supercricital. Here the profiles for the velocity and thermodynamic variables are dominated by the wavenumber associated with the most unstable mode. For the range of Ekman numbers spanned by our simulation suite, the resulting Busse columns are mostly localized in the outer half of the shell by radius and at low latitudes (see Jones et al. 2009 for a linear instability analysis of the problem). Figures 10(a,  b) show the instantaneous convective radial velocity and convective colatitudinal energy transport in the highly diffusive case FE10-3, which lies in the barely supercritical regime. Each upflow and downflow (pairs of which trace one Busse column roll) has an associated colatitudinal energy transport that is, on average, negative in the northern hemisphere and positive in the southern hemisphere, implying preferentially poleward energy transport.
The geostrophic nature of the Busse columns is illustrated in Figure 11, as is the resulting axial component of the flow. In the top row (case FE10-3), panel a shows that the Busse-column rolls alternate between high and low pressure. Panel b shows that the high-pressure rolls are each anticyclonic (have negative vorticity), while the low-pressure rolls are cyclonic. Finally, panel c shows that each high-pressure anomaly corresponds to poleward flow (v z > 0 in the northern hemisphere), while each low-pressure anomaly corresponds to equatorward flow (v z < 0). In the bottom row (the more supercritical case FF4-3), the Busse columns are less regularly spaced, but still largely alternate between anticyclonic regions of high pressure and cyclonic regions of low pressure (panels d, e). The axial flow associated with the Busse columns in case FF4-3 (panel f ) then leads to poleward energy transport through pressure work, just as in case FE10-3. Figures 10(c, d ) show the radial velocity and convective energy transport in the comparatively more turbulent case FF2-3. The flow structures are more intricate and fine-scale than in the barely supercritical regime, but the imprint of the most unstable mode remains: many Busse column rolls correspond to sites of negative colatitudinal energy transport in the northern hemisphere and positive transport in the southern hemisphere, again yielding preferentially poleward energy transport.

EFFECT OF OUTER THERMAL BOUNDARY CONDITION
Given that Busse columns direct energy poleward, spherical-shell convection simulations can achieve equilibrium only by forming conductive gradients that balances the poleward convective enthalpy flux. In general, such conductive transport can be achieved in two distinct ways. As the pole heats up and the equator-to-pole contrast increases, a latitudinal gradient will form that transports heat equatorward. Additionally, the increased temperature of the pole can lead to enhancement of the radial gradients in the outer thermal boundary layer, thus causing the poles to lose heat more efficiently (i.e., become superluminous). In the FF cases, the outer thermal boundary condition precludes the second of these options because the radial gradients are fixed. Hence, the FF models must rely solely on the development of a pole-to-equator conductive gradient. In the FE models, both types of gradients are possible. Therefore, the amount that the pole must be heated before equilibrium can be achieved is less for the FE models than it is for the FF models. The outer thermal boundary condition thus has a direct influence on the latitudinal contrast in the temperature, entropy, and pressure, with the FF boundary condition being conducive to strong contrast in all the thermodynamic variables. In the presence of thermal wind balance, the FF boundary condition thus leads to enhanced contrast in the differential rotation and positively-tilted isocontours in the rotation rate.
Mathematically, we illustrate the combined effects of the outer thermal boundary condition and latitudinal energy transport using the steady-state total energy equation for the fluid. Using Equations (1)-(3), this equation takes the form of a balance of fluxes: where is the temporally and longitudinally averaged total energy flux in the meridional plane and we have defined the averaged convective, conductive, radiative, viscous, and meridional-circulation fluxes through respectively. Note that ρT S + P = c p T , so the terms with T v and T v in the convective and meridionalcirculation fluxes represent the combined effects of heat advection and pressure work. Technically the flux due to transport of kinetic energy (proportional to v 2 v ) has convective terms (e.g., (v ) 2 v ) and meridionalcirculation terms (e.g., v 2 v ). For simplicity, we include all the kinetic-energy terms in the convective flux since they are in general small. We are interested in the total latitudinal transport of energy, and so we integrate the total flux in Equation (19) over conical surfaces at constant latitude: For the FF cases, there can be no net transport of energy in latitude due to the absence of conductive losses in the polar regions through the outer boundary. In other words, I θ (θ) ≡ 0. In the FE cases, by contrast, there is a net poleward energy transport because the pole is allowed to have enhanced conductive loss through the thermal boundary layer at the outer surface (i.e., the poles are superluminous). Thus, I θ (θ) will in general be negative in the northern hemisphere and positive in the southern hemisphere. Figure 12(a) shows the total latitudinal energy flux in case FF2-3 after the system has achieved statistical equilibrium. The total flux is very close to zero at all latitudes, indicating a well-equilibrated state. The dominant transport components are the convective flux, which transports energy preferentially poleward due to the Busse columns, and the conductive flux, which transports energy equatorward. The monotonic entropy gradient of Figure 8(d ), and by extension radially tilted contours in the FF cases, is thus seen to be a result of the response by conduction to the convective transport of energy to the poles. Figure 12(b) shows the total latitudinal energy transport in case FE2-3.
The poles are clearly superluminous-i.e., there is a net poleward energy transport due to the convection. For all the FE cases explored here, the energy loss at the poles is even greater than the heating by the convection; the conductive flux is thus forced to change sign at mid-latitudes (around ±20 • ), transporting energy poleward in concert with the Busse columns. This results in the non-monotonic entropy profile of Figure 8(a), leading to the tilts of the rotation contours being negative at high latitudes and positive at low latitudes, as seen in the top row of Figure  4.

DISCUSSION AND CONCLUSIONS
We have shown that the differential rotation achieved in global, 3D convection simulations is well-described by a thermal wind and highly sensitive to the outer thermal boundary condition. The FF boundary tends to yield more solar-like rotation profiles (strong contrast with positively-tilted contours), while the FE boundary yields weaker contrast and negatively-tilted contours. In light of these results, we now discuss the likelihood that the Sun's strong rotation contrast and positively-tilted contours arise from thermal wind balance in the deep interior and the fact that the radiative flux from the solar photosphere does not vary appreciably with latitude.
The first question is whether the force balance Equation (15), which should in general be true for low Rossby numbers, holds in the Sun. The interior solar Rossby number is currently unknown, but recent helioseismic estimates (Hanasoge et al. 2012;Greer 2015) give Ro 0.1 in the deep interior. Thus, it is likely that Equation (15) (and the derivative thermal wind Equations (16) and (17)) apply in the solar CZ, except perhaps in the layers just below the photosphere. Thermal wind balance in the Sun is also not prohibited by observational results. Back-of-the-envelope estimates (see the Appendix) show that the latitudinal thermodynamic contrasts required to drive a solar-like differential rotation through a thermal wind are no more than ∼1 part in 10 5 . At mid-depth, this corresponds to a latitudinal temperature contrast of ∼25 K, which is well within the constraints of helioseismology (e.g., Brun et al. 2010).
The second question is whether the Sun's Busse columns send energy preferentially poleward. In general, stellar convection transitions through a series of convective regimes as the supercriticality (measured by the reduced Rayleigh number R) increases (Hindman et al. 2020). Both the least supercritical case in our work (FE10-3, for which R ∼ 2) and the most supercritical cases (the pair FF2-3 and FE2-3, for which R ∼ 30) have a strong preference for poleward transport by Busse columns, suggesting that the poleward transport is a feature of the most unstable mode of convection that stays present as the flows get ever more complex. Of course, it is currently unknown in which regime the highly supercritical Sun (for which R ∼ 10 10 ) resides, and so future work must assess the transport properties of Busse columns in more supercritical regimes than the cases presented here. Thermal coupling to the tachocline (not considered in this work) may also play a role in the Sun's latitudinal energy transport. Rempel (2005) suggested in the context of mean field theory that latitudinal entropy gradients can arise in the radiative interior from meridional circulations and transmit to the CZ, most likely providing an additional source of poleward heat transport. Finally, it is an open question how the Sun might transport energy equatorward to maintain equilibrium. In our simulations, the net poleward transport of energy by Busse columns is a few per cent of the solar luminosity, which is counteracted almost entirely by conduction in the FF cases ( Figure 12). In the Sun, both the thermal and radiative diffusivities are ∼10 7 cm 2 s −1 at mid-CZ, which (if a thermal wind were operating with a temperature contrast of ∼25 K) would correspond to a latitudinal flux of ∼10 −7 L . It is unclear whether the Busse-column transport in the Sun could be this low. There clearly is a tendency for the net transport to go down in our more supercritical simulations (to 0.019L in case FF2-3 from 0.031L in case FF4-3) but we cannot predict how this relationship scales to a Sun-like supercriticality. It is also possible that the poleward transport by the solar Busse-columns is higher than the above estimate and the equatorward transport is accomplished by a different mechanism, such as magnetic fields (i.e., a Poynting flux) or flux from the meridional circulation.
On a more practical note, it is important to use an FF outer boundary condition in solar simulations for two reasons. First, maintaining a strong differential rotation is particularly relevant for dynamo models, since the dynamo cannot be excited without sufficient shear (e.g., Brown et al. 2010;Guerrero et al. 2016;Matilsky & Toomre 2020;Bice & Toomre 2020). Second, using an FE outer boundary condition leads to superluminous poles, which are directly at odds with solar observations. Figure 13 shows the conductive flux as a function of latitude at the top of the domain for the FE cases. For the most turbulent case FE2-3, the polar flux reaches a value in excess of the solar luminosity by about 20%. This is far greater than the observationally-constrained value of < 1% for the Sun (Rast et al. 2008).
We very much view this paper as a complement to Miesch et al. (2006). In that work, a systematic tilt of the rotation contours was achieved by imposing a latitudinal entropy gradient at the inner boundary. And indeed, for all our FF cases, there is a similar entropy gradient at the inner boundary, which is absent in the FE cases. In other words, Miesch et al. (2006) showed that it is possible to tilt the rotation contours by imposing a preferred geostrophic balance, and here we have shown how this preferred balance is naturally established as the result of poleward energy transport by Busse columns and the FF outer boundary condition.
We conclude that the observed strong solar differential rotation with positive tilt in the Sun may be a consequence of (1) the tendency of Busse columns to transport heat preferentially poleward, (2) the presence of a thermal wind operating in the deep solar interior, and (3) the invariance of the solar radiative flux with latitude at the photosphere. There are no doubt other processes not considered here that play a role, such as magnetism, coupling to the tachocline, and modified turbulent transport in the highly supercritical solar regime. Nevertheless, in the appropriate convective regime, the three effects above can lead to a differential rotation similar to the Sun's.