MODELING OF DIFFERENTIAL ROTATION IN RAPIDLY ROTATING SOLAR-TYPE STARS

We do not use the anelastic approximation here. The equations in an inertial frame can be expressed as

$$\begin{eqnarray} && \frac{\partial \rho _1}{\partial t}= -\frac{1}{r^2}\frac{\partial }{\partial r}(r^2v_r\rho _0) -\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta v_\theta \rho _0), \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial v_r}{\partial t}&=& -v_r\frac{\partial v_r}{\partial r} -\frac{v_\theta }{r}\frac{\partial v_r}{\partial \theta } +\frac{v_\theta ^2}{r} -\frac{1}{\rho _0} \left[ \rho _1 g+\frac{\partial p_1}{\partial r} \right]\nonumber\\ &&+\big(2\Omega _0\Omega _1+\Omega _1^2\big)r\sin ^2\theta +\frac{F_r}{\rho _0}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial v_\theta }{\partial t}&=& {-}v_r\frac{\partial v_\theta }{\partial r} -\frac{v_\theta }{r}\frac{\partial v_\theta }{\partial \theta }-\frac{v_rv_\theta }{r} -\frac{1}{\rho _0}\frac{1}{r}\frac{\partial p_1}{\partial \theta }\nonumber\\ &&+\big(2\Omega _0\Omega _1+\Omega _1^2\big)r\sin \theta \cos \theta +\frac{F_\theta }{\rho _0}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial \Omega _1}{\partial t}&=& {-}\frac{v_r}{r^2}\frac{\partial }{\partial r}[r^2(\Omega _0+\Omega _1)]\nonumber\\ &&-\frac{v_\theta }{r\sin ^2\theta }\frac{\partial }{\partial \theta } [\sin ^2\theta (\Omega _0+\Omega _1)] +\frac{F_\phi }{\rho _0 r\sin \theta }, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial s_1}{\partial t}&=& {-}v_r\frac{\partial s_1}{\partial r} -\frac{v_\theta }{r}\frac{\partial s_1}{\partial \theta } +v_r\frac{\gamma \delta }{H_p}\nonumber\\ &&+\frac{\gamma -1}{p_0}Q +\frac{1}{\rho _0 T_0}\mathrm{div}(\kappa _\mathrm{t}\rho _0T_0\mathrm{grad} s_1), \end{eqnarray} \tag{ 5 }$$

where Ω₀ is a constant value that represents the angular velocity of the rigidly rotating radiative zone. We set it as a parameter in Table 1. γ is the ratio of specific heats, with the value for an ideal gas being γ = 5/3. κ_t is the coefficient of turbulent thermal conductivity. δ = ∇ − ∇_ad represents superadiabaticity, where ∇ = d(ln T)/d(ln p) (see Section 2.2). g denotes gravitational acceleration. Following from this, the perturbation of pressure p₁ and pressure scale height H_p are expressed as

$$\begin{eqnarray} && p_1=p_0 \left(\gamma \frac{\rho _1}{\rho _0}+s_1 \right), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} && H_p=\frac{p_0}{\rho _0 g}. \end{eqnarray} \tag{ 7 }$$

s₁ is dimensionless entropy normalized by the specific heat capacity at constant volume c_v. Turbulent viscous force F follows from

$$\begin{eqnarray} &&F_r=\frac{1}{r^2}\frac{\partial }{\partial r}(r^2R_{rr}) +\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta R_{\theta r}) -\frac{R_{\theta \theta }+R_{\phi \phi }}{r},\nonumber\\ \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&F_\theta =\frac{1}{r^2}\frac{\partial }{\partial r}(r^2R_{r\theta }) +\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta R_{\theta \theta }) +\frac{R_{r\theta }-R_{\phi \phi }\cot \theta }{r},\nonumber\\ \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&F_\phi =\frac{1}{r^2}\frac{\partial }{\partial r}(r^2R_{r\phi }) +\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta R_{\theta \phi }) +\frac{R_{r\phi }+R_{\theta \phi }\cot \theta }{r},\nonumber\\ \end{eqnarray} \tag{ 10 }$$

with the Reynolds stress tensor

$$\begin{eqnarray} &&R_{ik}=\rho _0 \left[ \nu _\mathrm{tv} \left(E_{ik}-\frac{2}{3}\delta _{ik}\mathrm{div}\,{\bf v} \right) +\nu _\mathrm{tl}\Lambda _{ik} \right]. \end{eqnarray} \tag{ 11 }$$

Here ν_tv is the coefficient of turbulent viscosity and ν_tl is the coefficient of the Λ effect (Kitchatinov & Rüdiger 1995), a non-diffusive angular momentum transport caused by turbulence. ν_tv and ν_tl are expected to have the same value, since both effects are caused by turbulence, i.e., thermal-driven convection. We discuss this in more detail in Section 2.3. E_ik denotes the deformation tensor, which is given in spherical coordinates by

$$\begin{eqnarray} &&E_{rr}=2\frac{\partial v_r}{\partial r}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&E_{\theta \theta }=2\frac{1}{r}\frac{\partial v_\theta }{\partial \theta }+2\frac{v_r}{r}, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&E_{\phi \phi }=\frac{2}{r}(v_r+v_\theta \cot \theta), \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&E_{r\theta }=E_{\theta r}=r\frac{\partial }{\partial r} \left(\frac{v_\theta }{r}\right)+\frac{1}{r}\frac{\partial v_r}{\partial \theta }, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&E_{r\phi }=E_{\phi r}=r\sin \theta \frac{\partial \Omega _1}{\partial r}, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&E_{\theta \phi }=E_{\phi \theta }=\sin \theta \frac{\partial \Omega _1}{\partial \theta }. \end{eqnarray} \tag{ 17 }$$

An expression for the Λ effect (Λ_ik) is given later. The amount of energy that is converted by the Reynolds stress from kinematic energy to internal energy is given by

$$\begin{eqnarray} &&Q=\sum _{i,k}\frac{1}{2}E_{ik}R_{ik}. \end{eqnarray} \tag{ 18 }$$

Table 1. Significant Parameters of the Simplified Model

Case	Ω0	ν0v	Λ0	λ	δc
	(nHz)	(cm2 s−1)
1	1Ω☉ = 430	3 × 1012	1	15°	0
2	2Ω☉ = 860	3 × 1012	1	15°	0
3	4Ω☉ = 1720	3 × 1012	1	15°	0
4	8Ω☉ = 3440	3 × 1012	1	15°	0
5	16Ω☉ = 6880	3 × 1012	1	15°	0
6	1Ω☉ = 430	12 × 1012	1	15°	0
7	1Ω☉ = 430	6 × 1012	1	15°	0
8	1Ω☉ = 430	1.5 × 1012	1	15°	0
9	1Ω☉ = 430	3 × 1012	2	15°	0
10	1Ω☉ = 430	3 × 1012	0.5	15°	0
11	1Ω☉ = 430	3 × 1012	0.25	15°	0
12	1Ω☉ = 430	3 × 1012	1	25	0
13	1Ω☉ = 430	3 × 1012	1	15°	1 × 10−6
14	2Ω☉ = 860	3 × 1012	1	15°	1 × 10−6
15	4Ω☉ = 1720	3 × 1012	1	15°	1 × 10−6
16	8Ω☉ = 3440	3 × 1012	1	15°	1 × 10−6
17	16Ω☉ = 6880	3 × 1012	1	15°	1 × 10−6

Download table as:  ASCIITypeset image

We use an adiabatic hydrostatic stratification for the spherically symmetric reference state of ρ₀, p₀, and T₀. Gravitational acceleration is assumed to have ∼r⁻² dependence, since the radiative zone (r < 0.65 R_☉) has most of the solar mass. This is expressed as

$$\begin{eqnarray} && \rho _0(r)=\rho _\mathrm{bc} \left[ 1+\frac{\gamma -1}{\gamma }\frac{r_\mathrm{bc}}{H_\mathrm{bc}} \left(\frac{r_\mathrm{bc}}{r}-1\right) \right]^{1/(\gamma -1)}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} && p_0(r)=p_\mathrm{bc} \left[ 1+\frac{\gamma -1}{\gamma }\frac{r_\mathrm{bc}}{H_\mathrm{bc}} \left(\frac{r_\mathrm{bc}}{r}-1\right) \right]^{\gamma /(\gamma -1)}, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} && T_0(r)=T_\mathrm{bc} \left[ 1+\frac{\gamma -1}{\gamma }\frac{r_\mathrm{bc}}{H_\mathrm{bc}} \left(\frac{r_\mathrm{bc}}{r}-1\right) \right], \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} && g(r)=g_\mathrm{bc} \left(\frac{r}{r_\mathrm{bc}} \right)^{-2}, \end{eqnarray} \tag{ 22 }$$

where ρ_bc, p_bc, T_bc, H_bc = p_bc/(ρ_bcg_bc), and g_bc denote the values at the base of the convection zone r = r_bc of density, pressure, temperature, pressure scale height, and gravitational acceleration, respectively. In this study we use r_bc = 0.71 R_☉, with R_☉ representing the solar radius (R_☉ = 7 × 10¹⁰ cm). We adopt solar values ρ_bc = 0.2 g cm⁻³, p_bc = 6 × 10¹³ dyn cm⁻², T_bc = mp_bc/(k_Bρ_bc) ∼ 1.82 × 10⁶ K, and g_bc = 5.2 × 10⁴ cm s⁻², where k_B is the Boltzmann constant and m is the mean particle mass. Figure 1 shows the profiles of background density, pressure and temperature, and gravitational acceleration.

Zoom In Zoom Out Reset image size

Although the real Sun's stratification is not adiabatic in the convection zone, our reference state is valid, since the absolute value of superadiabaticity is small. In order to include the deviation from adiabatic stratification, we assume that superadiabaticity δ has the following profile:

$$\begin{eqnarray} &&\delta =\delta _\mathrm{conv}+\frac{1}{2}(\delta _\mathrm{os}-\delta _\mathrm{conv}) \left[ 1-\tanh \left(\frac{r-r_\mathrm{tran}}{d_\mathrm{tran}} \right) \right]. \end{eqnarray} \tag{ 23 }$$

Here δ_os and δ_conv denote the values of superadiabaticity in the overshoot region. r_tran and d_tran denote the position and the steepness of the transition toward the subadiabatically stratified overshoot region, respectively. Superadiabaticity in the convection zone is defined as

$$\begin{eqnarray} &&\delta _\mathrm{conv}=\delta _\mathrm{c}\frac{r-r_\mathrm{sub}}{r_\mathrm{max}-r_\mathrm{sub}}, \end{eqnarray} \tag{ 24 }$$

where r_max denotes the location of the upper boundary. We specify δ_os = −1.5 × 10⁻⁵, r_tran = 0.725 R_☉, r_sub = 0.8 R_☉, and d_tran = d_sub = 0.0125 R_☉ in our simulations. δ_c is considered as a free parameter. The entropy gradient can be expressed as

$$\begin{eqnarray} &&\frac{ds_0}{dr}=-\frac{\gamma \delta }{H_p}. \end{eqnarray} \tag{ 25 }$$

The third term of Equation (5), v_rγδ/H_p, includes the effect of deviations from adiabatic stratification. The term indicates that an upflow (downflow) can make negative (positive) entropy perturbations in the subadiabatically stratified layers (δ < 0).

We assume the coefficients of turbulent viscosity and thermal conductivity to be constant within the convection zone, and these smoothly connect with the values of the overshoot region. We assume that the diffusivities only depend on the radial coordinate:

$$\begin{eqnarray} && \nu _\mathrm{tv}=\nu _\mathrm{os}+\frac{\nu _\mathrm{0v}}{2} \left[ 1+\tanh \left(\frac{r-r_\mathrm{tran}+\Delta }{d_{\kappa \nu }} \right) \right]f_c(r), \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} && \nu _\mathrm{tl}=\frac{\nu _\mathrm{0l}}{2} \left[ 1+\tanh \left(\frac{r-r_\mathrm{tran}+\Delta }{d_{\kappa \nu }} \right) \right]f_c(r), \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} && \kappa _\mathrm{t}=\kappa _\mathrm{os}+\frac{\kappa _0}{2} \left[ 1+\tanh \left(\frac{r-r_\mathrm{tran}+\Delta }{d_{\kappa \nu }} \right) \right]f_c(r), \end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray} && f_c(r)=\frac{1}{2} \left[ 1+\tanh \left(\frac{r-r_\mathrm{bc}}{d_\mathrm{bc}} \right) \right], \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} && \Delta =d_{\kappa \nu }\tanh ^{-1}(2\alpha _{\kappa \nu }-1), \end{eqnarray} \tag{ 30 }$$

where ν_0v, ν_0l, and κ₀ are the values of the turbulent diffusivities within the convection zone, and ν_os and κ_os are the values in the overshoot region. We specify ν_0l = κ_0l = 3 × 10¹² cm² s⁻¹, ν_os = 6 × 10¹⁰ cm² s⁻¹, and κ_os = 6 × 10⁹ cm² s⁻¹, and we treat ν_0v as a parameter. α_κν specifies the values of the turbulent diffusivities at r = r_tran, i.e., ν_tv = ν_os + α_κνν_0v, ν_tl = α_κνν_0l, and κ_t = κ_os + α_κνκ₀ at r = r_tran. d_bc and d_κν are the widths of transition. We specify α_κν = 0.1, d_bc = 0.0125 R_☉, and d_κν = 0.025 R_☉. As already mentioned, the coefficients for turbulent viscosity and the Λ effect are different in our model from those of Rempel's (2005). There are two reasons for this. One is that we intend to investigate the influence of both effects on stellar differential rotation separately (see Section 4.2). The other reason is that the formation of a tachocline in a reasonable amount of time requires a finite value (though small) for the coefficient of turbulent viscosity even in the radiative zone, in which there is likely to be weak turbulence (Rempel 2005). Figure 2 shows the profiles of ν_tv, ν_tl, and κ_t.

Zoom In Zoom Out Reset image size

In this study, we adopt the non-diffusive part of the Reynolds stress, called the Λ effect. The Λ effect transports angular momentum and generates differential rotation. The Λ effect tensors are expressed as

$$\begin{eqnarray} && \Lambda _{r\phi }=\Lambda _{\phi r}=+L(r,\theta)\cos (\theta +\lambda), \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} && \Lambda _{\theta \phi }=\Lambda _{\phi \theta }=-L(r,\theta)\sin (\theta +\lambda), \end{eqnarray} \tag{ 32 }$$

where L(r, θ) is the amplitude of the Λ effect and λ is the inclination of the flux vector with respect to the rotational axis. We use for the amplitude of the Λ effect the expressions

$$\begin{equation} f(r,\theta)=\sin ^l\theta \cos \theta \tanh \left(\frac{r_\mathrm{max}-r}{d} \right), \end{equation} \tag{ 33 }$$

$$\begin{equation} L(r,\theta)=\Lambda _0\Omega _0\frac{f(r,\theta)}{\mathrm{max}|f(r,\theta)|}, \end{equation} \tag{ 34 }$$

where d = 0.025 R_☉. λ and Λ₀ are free parameters. The value of l needs to be equal to or larger than 2 to ensure regularity near the pole, so we set l = 2. The Λ effect does not depend on v_r, v_θ, or Ω₁, meaning it is a stationary effect. We emphasize that the Λ effect depends on stellar angular velocity Ω₀, since the Λ effect is generated by turbulence and Coriolis force. The more rapidly the star rotates, the more angular momentum the Λ effect can transport. The dependence of Λ₀ and λ on stellar angular velocity is discussed in Section 4.3.

Using the modified Lax–Wendroff scheme with total variation diminishing artificial viscosity (Davis 1984), we solve Equations (1)–(5) numerically for the northern hemisphere of the meridional plane in 0.65 R_☉ < r < 0.93 R_☉ and 0 < θ < π/2. We use a uniform resolution of 200 points in the radial direction and 400 points in the latitudinal direction in all of our simulations. Each simulation run is conducted until it reaches a stationary state. All the variables ρ₁, v_r, v_θ, Ω₁, and s₁ are equal to zero in the initial condition. At the top boundary (r = 0.93 R_☉) we adopt stress-free boundary conditions for v_r, v_θ, and Ω₁ and set the derivative of s₁ to zero:

$$\begin{equation} \frac{\partial v_r}{\partial r}=0, \end{equation} \tag{ 35 }$$

$$\begin{equation} \frac{\partial }{\partial r} \left(\frac{v_\theta }{r} \right)=0, \end{equation} \tag{ 36 }$$

$$\begin{equation} \frac{\partial \Omega _1}{\partial r}=0, \end{equation} \tag{ 37 }$$

$$\begin{equation} \frac{\partial s_1}{\partial r}=0. \end{equation} \tag{ 38 }$$

The boundary conditions for v_r, v_θ, and s₁ at the lower boundary (r = 0.65 R_☉) are the same as those at the top boundary. Differential rotation connects with the rigidly rotating core at the lower boundary, so we adopt Ω₁ = 0 there. At both radial boundaries, we set ρ₁ to make the right-hand side of Equation (2) equal zero. At the pole and the equator (θ = 0 and π/2) we use the symmetric boundary condition:

$$\begin{eqnarray} && \frac{\partial \rho _1}{\partial \theta }=0, \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} && \frac{\partial \Omega _1}{\partial \theta }=0, \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} && \frac{\partial v_r}{\partial \theta }=0, \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} && v_\theta =0, \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} && \frac{\partial s_1}{\partial \theta }=0. \end{eqnarray} \tag{ 43 }$$

Due to the low Mach number of the expected flows, a direct compressible simulation is problematic, so adopting the same technique as Rempel (2005), we reduce the speed of sound by multiplying the right-hand side of Equation (1) by 1/ζ². The equation of continuity is therefore replaced with

$$\begin{eqnarray} &&\frac{\partial \rho _1}{\partial t}+ \frac{1}{\zeta ^2}\mathrm{div}(\rho _0{\bf v})=0. \end{eqnarray} \tag{ 44 }$$

The speed of sound then becomes ζ times smaller than the original speed. We use ζ  =  200 in all our calculations. This technique can be used safely in our present study since we only discuss stationary states, so the factor ζ becomes unimportant. The validity of this technique is carefully discussed by Rempel (2005). We test our code by reproducing the results presented by Rempel (2005) and check the numerical convergence by runs with different grid spacings. After checking and cleaning up at every time step, conservation of total mass, total angular momentum, and total energy are maintained through the simulation runs.

In this section, we discuss the cases with angular velocities up to 16 times the solar value (represented by Ω_☉), placing an emphasis on the morphology of stellar differential rotation. Figure 3 shows the results of our calculations which correspond to cases 1–5 in Table 1. It is found that the larger the stellar angular velocity is, the more likely it is for differential rotation to be in the Taylor–Proudman state, in which the contour lines of the angular velocity are parallel to the rotational axis. To evaluate these results quantitatively, we define a parameter which denotes the morphology of differential rotation. We call it the non-Taylor–Proudman parameter (hereafter the NTP parameter), which is expressed as

$$\begin{eqnarray} P_\mathrm{ntp}&=&\frac{1}{R_\odot ^2\Omega _0^2}\int \frac{\partial \Omega _1^2}{\partial z}dV\nonumber\\ &=&\frac{1}{R_\odot ^2\Omega _0^2}\int \left(\cos \theta \frac{\partial }{\partial r} -\frac{\sin \theta }{r}\frac{\partial }{\partial \theta } \right) \Omega _1^2\, dV, \end{eqnarray} \tag{ 46 }$$

where Ω₀ is the angular velocity of the radiative zone. When the NTP parameter is zero, differential rotation is in the Taylor–Proudman state. Conversely, differential rotation is far from the Taylor–Proudman state with a large absolute value of the NTP parameter. The value of the NTP parameter with various stellar angular velocities is shown in Figure 4. The NTP monotonically decreases with increases in stellar angular velocity. These results indicate that with large stellar angular velocity values, differential rotation approaches the Taylor–Proudman state. These results are counterintuitive, however, since we do not expect differential rotation to approach the Taylor–Proudman state with increasing stellar angular velocity values, since the Λ effect, which is a driver of the deviation from the Taylor–Proudman state, is proportional to stellar angular velocity Ω₀. These are the most significant findings of this paper, so hereafter in this section we discuss these unexpected results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next discuss the temperature difference between the equator and the pole at the base of the convection zone (r = 0.71 R_☉). Since temperature is given as a function of entropy by

$$\begin{eqnarray} &&T_1=\frac{T_0}{\gamma } \left[s_1+(\gamma -1)\frac{p_1}{p_0} \right], \end{eqnarray} \tag{ 47 }$$

and it is easier to measure than entropy, we use it here for discussing the thermal structure of the simulation results in the convection zone. Further, although it is mentioned in Section 3 that the entropy gradient is crucial for breaking the Taylor–Proudman constraint, the temperature difference can be used as its proxy. Figure 5 shows the relationship between stellar angular velocity Ω₀ and temperature difference ΔT at r = 0.71 R_☉, where ΔT = max (T₁(r_bc, θ)) − min (T₁(r_bc, θ)). Although the temperature difference monotonously increases with larger stellar angular velocity values, it is not enough to make the rotational profile largely deviate from the Taylor–Proudman state. This can be explained by using the thermal wind equation, which is a steady state solution of Equation (45):

$$\begin{eqnarray} &&0= r\sin \theta \frac{\partial \Omega ^2}{\partial z} -\frac{g}{\gamma r}\frac{\partial s_1}{\partial \theta }. \end{eqnarray} \tag{ 48 }$$

The inertial term and the diffusion term are neglected here. This equation indicates that, for a given value of the NTP, we need an entropy gradient proportional to Ω²₀. However, our simulation results show that ΔT∝Ω^0.58₀, which means that as Ω₀ increases, the thermal driving force becomes insufficient to push differential rotation away from the Taylor–Proudman state. In other words, the latitudinal entropy gradient in rapidly rotating stars is so small that differential rotation stays close to the Taylor–Proudman state. In our model, meridional flow generates latitudinal entropy gradient at the base of the convection zone. It is conjectured that the insufficient thermal drive is due to a slow meridional flow.

Zoom In Zoom Out Reset image size

We next investigate the dependence of meridional flow on stellar angular velocity. Figure 6 shows the radial profile of latitudinal velocity v_θ at θ = 45°, using the results of cases 1, 2, and 9. In case 2, stellar angular velocity is twice that of case 1 (the solar value). In case 9, stellar angular velocity is equal to the solar value, and the amplitude of the Λ effect is two times the value in case 1. Figure 6 shows that meridional flow does not depend on stellar angular velocity, while it correlates with the Λ effect. Considering Equation (34), the Λ effect increases with larger values of stellar angular velocity, since the amplitude of the Λ effect is proportional to Ω₀. The reason why differential rotation in rapidly rotating stars is close to the Taylor–Proudman state is that meridional flow does not become fast with large stellar angular velocity values.

Zoom In Zoom Out Reset image size

We interpret the result that the speed of meridional flow does not depend on stellar angular velocity in our model as follows. With large values of stellar angular velocity, more angular momentum is transported by the Λ effect (note that the Λ effect is proportional to Ω₀ in Equation (34)), so meridional flow obtains more energy from differential rotation. The energy gain does not result in an increase in speed because of the associated enhancement of the Coriolis force, which bends the meridional flow in the longitudinal direction. Another explanation is possible in terms of angular momentum transport. The angular momentum fluxes from both meridional flow and the Reynolds stress (Λ effect) must be balanced in a steady state. The former is proportional to v_mΩ₀ and the latter is proportional to Ω₀, where v_m is the amplitude of meridional flow. Therefore, meridional flow does not depend on stellar angular velocity (Miesch 2005). Our results (Figure 6) indicate that with a larger stellar angular velocity (case 2), the above mechanism does not generate a fast meridional flow. However, this does not occur when only the Λ effect is large (case 9).

In this subsection, we discuss angular velocity difference ΔΩ at the surface and the relationship between our results and previous observations. We conduct numerical simulations to investigate the physical process which determines ΔΩ (cases 1 and 6–11). We define angular velocity difference as ΔΩ = max (Ω₁(r_max, θ)) − min (Ω₁(r_max, θ)).

ΔΩ is determined by two opposing effects, a smoothing effect from turbulent viscosity and a steepening effect from the Λ effect. In a stationary state these two effects cancel each other out. Latitudinal flux for turbulent viscosity and the Λ effect can be written as ρ₀ν_0vΔΩ/Δθ and ρ₀ν_0lΛ₀Ω₀, respectively. Because these two have approximately the same value, ΔΩ can be estimated as

$$\begin{eqnarray} &&\Delta \Omega \sim \frac{\nu _\mathrm{0l}}{\nu _\mathrm{0v}}\Lambda _0\Omega _0\Delta \theta, \end{eqnarray} \tag{ 49 }$$

where Δθ denotes the differential rotation region.

In order to confirm Equation (49), we conduct two sets of simulations, first varying the value of turbulent viscosity (ν_0v) and second the amplitude of the Λ effect (Λ₀). Note that the setting for turbulent viscosity does not reflect a real situation, since the coefficients of turbulent viscosity and the Λ effect should have a common value. Nonetheless, this is necessary for the purpose of our investigation. The simulation results are shown in Figures 7 and 8. We obtain ΔΩ∝ν^−0.88_0v and ΔΩ∝Λ^1.1₀, which are consistent with Equation (49).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the results of the dependence of ΔΩ on Ω₀ (cases 1–5). Asterisks denote the difference at the surface between the equator and the pole, squares show the difference between the equator and the colatitude θ = 45°, and triangles are the difference between the equator and the colatitude θ = 60°. The difference at low latitudes (squares and triangles) monotonically increases with stellar angular velocity. However, this is not the case for angular velocity difference between the equator and the pole (asterisk). As we discussed in Section 4.1, when stellar rotation velocity is large, the Taylor–Proudman state is achieved, meaning that the gradient of angular velocity at the surface concentrates in lower latitudes. Due to this concentration, Δθ becomes smaller in Equation (49) with larger values of Ω₀. Thus, ΔΩ₀ does not show an explicit dependence on Ω₀. At low latitudes, Δθ is fixed and the angular velocity difference increases with stellar angular velocity. We obtain ΔΩ∝Ω^0.43₀ (between the equator and the colatitude θ = 45°: squares) and ΔΩ∝Ω^0.55₀ (between the equator and the colatitude θ = 60°: triangles). This indicates that ΔΩ/Ω₀ decreases with stellar angular velocity. These results are consistent with previous stellar observations (Donahue et al. 1996; Reiners & Schmitt 2003; Barnes et al. 2005).

Zoom In Zoom Out Reset image size

In this section, we discuss the dependence of meridional flow and differential rotation on free parameters. The parameter set is shown in Table 1 (cases 12–17). At first we investigate the influence of the variation of the Λ effect. The Λ effect has two free parameters, i.e., amplitude Λ₀ and inclination angle λ (see Section 2.4). Amplitude is thought to become smaller with a larger stellar angular velocity, due to the saturation of the correlations such as 〈v'_rv'_ϕ〉 and 〈v'_θv'_ϕ〉, where v'_r, v'_θ, and v'_ϕ are the radial, latitudinal, and longitudinal component turbulent velocities, respectively. Figure 6 shows that meridional flow becomes slower with a smaller Λ₀, keeping the Ω₀ value constant (case 10). It is clear with the result of Section 4.1 that meridional flow becomes slow with a larger angular velocity when the variation of Λ₀ is included. Brown et al. (2008) reported this effect with their three-dimensional hydrodynamic calculation. When meridional flow is slow, the entropy gradient generated by the subadiabatic layer is small and differential rotation approaches the Taylor–Proudman state.

The inclination angle is thought to be small with large stellar angular velocity values, since the motion across the rotational axis is restricted (Kichatinov & Rüdiger 1993). In case 12, differential rotation with a small inclination angle (λ = 25) is calculated. Other parameters are the same as case 1. The radial distribution of meridional flow is shown in Figure 6. Meridional flow becomes faster with a smaller inclination angle. Because of the efficient angular momentum transport in the z-direction when the inclination angle is small, the second term on the right-hand side of Equation (45) is large. This generates a large ω_ϕ, i.e., fast meridional flow.

In summary, we found that rapid stellar rotation causes two opposing effects on the speed of meridional flow. The speed is reduced by the suppression of Λ₀, while it is enhanced by the angular momentum transport along the axial direction with a smaller λ. Although the results of the three-dimensional calculation suggest that meridional flow becomes slower with a larger stellar angular velocity, our model cannot draw a conclusion about the speed of meridional flow in rapidly rotating stars.

Next, we investigate the influence of superadiabaticity in the convection zone. In cases 13–17, superadiabaticity in the convection zone δ_c = 1 × 10⁻⁶. The differences of the NTP parameters with adiabatic and superadiabatic convection zones $(P_{\mathrm{ntp}(\delta _\mathrm{c}=0)}-P_{\mathrm{ntp}(\delta _\mathrm{c}=10^{-6})})/P_{\mathrm{ntp}(\delta _\mathrm{c}=0)}$ are shown in Figure 10. The NTP parameter values with a superadiabatic convection zone are smaller than those with an adiabatic convection zone, since meridional flow in the superadiabatic convection zone makes the entropy gradient small. This result is suggested by Rempel (2005). Note that the difference between the values of the NTP parameters with an adiabatic and those with a superadiabatic convection zone decreases as the stellar angular velocity increases, since the generation of entropy gradient by the subadiabtic layer becomes ineffective with a larger stellar angular velocity.

Zoom In Zoom Out Reset image size