text:
# Magnetic Fields in Massive Stars. II. The Buoyant Rise of Magnetic Flux Tubes Through the Radiative Interior

## 1 Introduction

Although magnetic fields and related activity are believed to be nearly ubiquitous among stars of solar and lower mass, it is less certain that magnetism is a common characteristic of stars more massive than the Sun. To date, definitive detections of fields in stars with masses $`\mathrm{}>1.5M_{}`$ have, for the most part, been made for objects having anomalous chemical abundances (e.g., the chemically peculiar A and B stars; Landstreet 1992; Donati 1998). Recently, however, observations of cyclic variability in the properties of winds from luminous OB stars have been interpreted as evidence for the presence of large-scale magnetic fields in the surface layers and atmospheres of these objects (see, e.g., Kaper & Henrichs 1994; Kaper et al. 1997; Prinja 1998). These inferences have been bolstered by the unambiguous measurement of a weak ($`360`$ G) field in the chemically normal B1 IIIe star $`\beta `$ Cephei (Donati et al. 2001). These results suggest that magnetic fields of moderate strength might be more prevalent among hot stars than had previously been thought (see Charbonneau & MacGregor 2001, hereafter Paper I, for a more complete discussion).

At the present time, the origin of magnetism in massive stars is not well understood. It is unclear whether the inferred fields are fossil in nature, or are instead the product of dynamo activity in the stellar interior (see, e.g., Parker 1979). The latter explanation, if appropriate, raises additional interesting questions concerning the site of magnetic field generation inside hot stars. For stars having spectral types O and B, convection occurs in the innermost portion of the core, not in the outer envelope as in low-mass stars. Because convection is thought to be necessary in order that field regeneration via the so-called $`\alpha `$-effect take place, it follows that a hot star dynamo should be located deep within the stellar interior. Within the context of mean-field electrodynamics, the properties of kinematic core dynamo models for early-type stars have been investigated by Levy & Rose (1974) and Schüssler & Pähler (1978), and, more recently, in Paper I.

If the magnetic field of a hot star is produced by dynamo action in the convective core, then a mechanism for transporting the field to the stellar surface must be identified. The finite electrical conductivity of the envelope leads to the outward diffusion of any fields contained therein, but only over an extended period of time. Estimates indicate that for stars more massive than a few solar masses, the resistive diffusion time across the radiative interior exceeds the main sequence lifetime $`t_{MS}`$ (Schüssler & Pähler 1978). Another possibility is that dynamo fields are advected from the core to the surface by rotation-induced meridional circulation (see Paper I). For a star of mass $`M_{}`$ and equatorial rotation speed $`v_{rot}`$, a rough estimate for the circulation time $`t_c`$ as a fraction of $`t_{MS}`$ is $`(t_c/t_{MS})9\times 10^2(M_{}/M_{})^{0.8}(v_{rot}/1\mathrm{km}\mathrm{s}^1)^2`$ (Kippenhahn & Weigert 1994). In the case of a star with $`M_{}=10M_{}`$ and $`v_{rot}=150`$ km s<sup>-1</sup>, $`(t_c/t_{MS})0.25`$, indicating that even for relatively rapidly rotating stars, a time $`t_{MS}`$ is required to bring the field produced by the dynamo to the surface by this mechanism. In addition, the results of Paper I suggest that the circulatory flow generated by extremely rapid rotation is capable of interfering with the operation of a core dynamo.

Alternatively, in the Sun, magnetic flux emerges at the photosphere in the form of fibrils or flux tubes. The field is thought to assume this form in or near the dynamo domain at the bottom of the convective zone. The subsequent rise of a flux tube to the solar surface is driven by buoyancy, a consequence of the reduced density inside a tube that is in mechanical and thermal equilibrium with the surrounding, adiabatically stratified, field-free gas. If the formation of flux tubes is a process that is not specific to fields in the Sun, then a dynamo-generated field at the base of the extended radiative envelope of a hot star might also develop fibril structure. In this case, the buoyant force might likewise enable flux from deep in the interior to reach the surface in a time that is shorter than evolutionary time scales. An important distinction, however, is that unlike the the solar convection zone, the stratification within the radiative interior of a hot star is sub-adiabatic. The motion of buoyant magnetic elements in a thermodynamic environment like that in the envelope of an upper main sequence star has been studied by Gurm & Wentzel (1967) and Moss (1989).

In the present paper, we adopt the premise that hydromagnetic dynamo activity takes place inside the convective core of a rotating massive star, in a manner similar to that described in Paper I. We furthermore presume that the fields so-produced naturally form into discrete, toroidal flux tubes that remain in total pressure equilibrium with the external, unmagnetized stellar interior at all times. Using this conceptual framework in concert with the physical model described in §2, we determine the time-dependent position, velocity, and thermodynamic properties of a buoyant flux ring that begins its outward motion from a specified location on the core boundary. Results pertaining to the dynamics and rise times of rings with a variety of initial field strengths and cross-sectional radii are presented in §3. Our findings regarding the dynamical behavior of magnetic rings in a radiative environment and the efficacy of buoyancy as a flux transport mechanism in hot stars are summarized in §4.

## 2 Model

To investigate the buoyant transport of magnetic flux through the radiative interior of a hot star, we consider a spherical star of mass $`M_{}`$ and radius $`R_{}`$ that contains a central, convective core of mass $`M_c(<M_{})`$ and radius $`R_c(<R_{})`$. The star is assumed to rotate uniformly with angular velocity $`\mathrm{\Omega }`$ around an axis that coincides with the polar axis of a spherical coordinate system $`(r,\theta ,\varphi )`$. In the radiative region external to the core, we study the dynamics of an isolated concentration of azimuthally directed magnetic field (i.e., a flux tube) that at any given time takes the form of a circular ring, symmetric about the stellar rotation axis. A similar configuration has been utilized by Choudhuri & Gilman (1987) to examine the influence of rotation on the rise of flux tubes in the solar convection zone. We assume that the cross-sectional profile of the tube is a circle whose radius $`a`$ is much smaller than either the radius of curvature of the ring or the local scale height of the surrounding, unmagnetized, radiative envelope. The assumed geometry for the rising flux tube is illustrated in Figure 1. These approximations are here used to study the dynamics for the rise of magnetic flux to the stellar surface, and for comparisons of rise times with the main sequence lifetime of the star.

As it is described above, the magnetic ring satisfies the criteria for applicability of the thin flux tube approximation (e.g., Defouw 1976; Spruit 1981), and its internal properties can be reasonably taken to be uniform throughout. The ring is assumed to contain material with mass density $`\rho `$, so that its total mass is $`m=\rho V`$, where $`V=2\pi r\mathrm{sin}\theta \pi a^2`$ is the ring volume. It follows from conservation of mass that if a ring with cross-sectional radius $`a_0`$ is released from position $`(r_0,\theta _0)`$ at time $`t=0`$, its density at some $`t(>0)`$ is

$$\rho =\rho _0\left(\frac{a_0}{a}\right)^2\left(\frac{r_0\mathrm{sin}\theta _0}{r\mathrm{sin}\theta }\right),$$
(1)
where $`a`$ and $`(r,\theta )`$ are the tube radius and position at the later time. Likewise, the flux ring is presumed to be untwisted and to contain a toroidal magnetic field of the form $`𝐁=B𝐞_\varphi `$ with $`B`$ a constant. Conservation of magnetic flux then implies that $`B`$ is related to its initial value $`B_0`$ through the relation

$$B=B_0\left(\frac{a_0}{a}\right)^2.$$
(2)
A further consequence of the thin flux tube approximation is that because the time required for a fast magnetosonic wave to transit the tube cross section is short, the total pressure inside the ring can be assumed to instantaneously equilibrate with that of the surrounding stellar material. The quantitative expression of this lateral mechanical equilibrium condition is

$$p_e=p+\frac{B^2}{8\pi },$$
(3)
where $`p=\rho kT/\mu `$ is the gas pressure and the subscript $`e`$ refers to conditions in the external, field-free medium on the periphery of the tube. Note that by using the ideal gas law for $`p`$ and $`p_e`$, we are neglecting the radiative component of the total pressure. This omission restricts our analysis to stars less massive than about $`10M_{}`$, for which radiation pressure is a minor contributor to the support of the deep stellar interior. In addition, we assume that the mean molecular weight $`\mu `$ has the same constant value throughout the star. Our model thus describes conditions within a young, chemically homogeneous, main sequence star; it does not apply at later evolutionary stages when a flux tube that originates in the vicinity of the core-envelope interface may contain chemically processed material.

In a frame of reference that rotates with the stellar angular velocity, the uniform flux ring moves coherently in response to the axisymmetric forces applied to it; that is, each mass element within it has the same velocity

$$𝐮=u_r𝐞_r+u_\theta 𝐞_\theta +u_\varphi 𝐞_\varphi =\frac{\mathrm{d}r}{\mathrm{d}t}𝐞_r+r\frac{\mathrm{d}\theta }{\mathrm{d}t}𝐞_\theta +r\mathrm{sin}\theta \frac{\mathrm{d}\varphi }{\mathrm{d}t}𝐞_\varphi .$$
(4)
The velocity components $`(u_r,u_\theta ,u_\varphi )`$ and ring position $`(r,\theta )`$ can be obtained as functions of time by integration of the ring equation of motion,

$$\frac{\mathrm{d}𝐮}{\mathrm{d}t}+2𝛀\times 𝐮=\left(\frac{\rho _e\rho }{\rho _e+\rho }\right)\left[\frac{GM(r)}{r^2}𝐞_r\mathrm{\Omega }^2r(\mathrm{sin}^2\theta 𝐞_r+\mathrm{sin}\theta \mathrm{cos}\theta 𝐞_\theta )\right]$$

$$\frac{B^2}{4\pi r(\rho _e+\rho )}(𝐞_r+\mathrm{cot}\theta 𝐞_\theta )C_D\left(\frac{\rho _e}{\rho _e+\rho }\right)\frac{|𝐮_{}|𝐮_{}}{\pi a},$$
(5)
where $`M(r)`$ is the internal mass distribution of the star, $`C_D`$ is the drag coefficient, and $`𝐮_{}u_r𝐞_r+u_\theta 𝐞_\theta `$ is the transverse velocity. The individual components of the inertial and Coriolis terms on the left-hand side of equation (5) have the explicit forms

$$\left(\frac{\mathrm{d}𝐮}{\mathrm{d}t}+2𝛀\times 𝐮\right)_r=\frac{\mathrm{d}u_r}{\mathrm{d}t}\frac{(u_\theta ^2+u_\varphi ^2)}{r}2\mathrm{\Omega }u_\varphi \mathrm{sin}\theta ,$$
()
$$\left(\frac{\mathrm{d}𝐮}{\mathrm{d}t}+2𝛀\times 𝐮\right)_\theta =\frac{\mathrm{d}u_\theta }{\mathrm{d}t}+\frac{u_ru_\theta }{r}\frac{u_\varphi ^2}{r}\mathrm{cot}\theta 2\mathrm{\Omega }u_\varphi \mathrm{cos}\theta ,$$
()
$$\left(\frac{\mathrm{d}𝐮}{\mathrm{d}t}+2𝛀\times 𝐮\right)_\varphi =\frac{\mathrm{d}u_\varphi }{\mathrm{d}t}+\frac{u_ru_\varphi }{r}+\frac{u_\theta u_\varphi }{r}\mathrm{cot}\theta +2\mathrm{\Omega }(u_r\mathrm{sin}\theta +u_\theta \mathrm{cos}\theta ).$$
()
The equation of motion (5) and its many variations have been used to treat the dynamics of buoyant flux tubes in the Sun’s convective envelope by numerous authors, including Choudhuri & Gilman (1987), Choudhuri (1989), Moreno-Insertis, Schüssler, & Ferriz-Mas (1992), Cheng (1992), Fan, Fisher, & DeLuca (1993), Ferriz-Mas & Schüssler (1993, 1994), Fan, Fisher, & McClymont (1994), Caligari, Moreno-Insertis, & Schüssler (1995), and Fan & Fisher (1996), among others. The three terms on the right-hand side of equation (5) describe, respectively, the accelerations produced by: ($`i`$) the buoyant force, including the effect of the centrifugal reduction of the local gravitational acceleration; ($`ii`$) the magnetic tension force, arising in response to outward, stretching motion of the ring in the plane perpendicular to the rotation axis of the star; and ($`iii`$), the aerodynamic drag force, derived from that exerted on a straight, circular cylinder immersed in a steady flow with upstream velocity directed perpendicular to the axis of symmetry (Goldstein 1938). The appearance of the factor $`(\rho _e+\rho )^1`$ rather than $`\rho ^1`$ in the expressions for each of the accelerations reflects the fact that we have made provision for the so-called virtual (or enhanced) inertia of the ring, the hydrodynamical resistance to acceleration experienced by a rigid body that moves through a fluid (see, e.g., Batchelor 1967). Considerable discussion has been devoted to the reasons for including (or not including) this effect in the equation describing the balance of forces acting on a flux tube (see, e.g., Moreno-Insertis, Schüssler, & Ferriz-Mas, 1996; Fan & Fisher 1996, and references therein). As will become clear in subsequent sections, the overall picture of flux ring dynamics that emerges from solutions of equation (5) depends little on the presence or absence of virtual inertia in any of the terms.

The difference between the nature of the energy balance that prevails within the envelope of a massive star and that within the convection zone of the Sun leads to significant differences between the dynamics of flux tubes located in the two regions. The outer portion of the solar interior is adiabatically stratified, a consequence of efficient convective energy transport therein. Because of this, an adiabatic flux tube that is initially buoyant will remain so during the course of its rise toward the surface. Alternatively, the envelope of a hot star is sub-adiabatically stratified, a consequence of the radiative equilibrium conditions that exist throughout it. As a result, the density deficit of an initially buoyant, adiabatic tube will diminish as it rises, causing the upward buoyant acceleration to decrease as well. When adiabatic cooling causes the temperature contrast of the tube relative to its surroundings to attain the value $`(T_eT)/T_e=(1+\beta )^1`$ where $`\beta =(8\pi p/B^2)`$, $`\rho =\rho _e`$, and the buoyancy of the tube is reduced to zero.

The discussion of the preceding paragraph indicates that within the envelope of a massive star, the non-adiabatic, radiative interaction between a flux ring and its environment will play an important role in determining the magnitude of the buoyant acceleration the ring experiences, and thus, the time required for it to emerge at the surface of the star (see, e.g., Parker 1975; Moreno-Insertis 1983; Fan & Fisher 1996). In the present model, the heating associated with the diffusion of radiation from the external medium into the tube is described by the thermal energy equation (see Appendix A),

$$\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{32}{3}\frac{(\gamma 1)\sigma T_e^4}{\kappa _e\rho _epa^2}\frac{\mathrm{\Delta }T}{T_e}q,$$
(7)
where $`\mathrm{\Delta }T=(T_eT)`$, and $`S`$ is the entropy-like quantity

$$S=\mathrm{ln}\left[\left(\frac{p}{p_0}\right)\left(\frac{\rho _0}{\rho }\right)^\gamma \right].$$
(8)
In equations (7) and (8), $`\gamma `$ is the ratio of specific heats, $`\sigma `$ is the Stefan-Boltzmann constant, and $`\kappa _e`$ is the Rosseland mean opacity in the material surrounding the tube. In the absence of radiative heating (i.e., when $`q=0`$), equation (7) indicates that $`S`$ is constant, and the tube behaves adiabatically.

A general question concerning the flux tube approach is how far we can consider the flux tube to be thin in the thermal sense. We assume the heating that arises from the inward diffusion of radiation affects the tube interior uniformly. In reality, the fact that the heating time scale and the tube radius are both finite implies a differential response by the tube to the heat input. In particular, if the time scale for the dynamical adjustment of the tube ($``$ the sound crossing time) is much shorter than the diffusion heating time scale, then the temperature in the outermost portion of the tube can become elevated relative to the center, and the resultant enhancement of the buoyant acceleration can cause a deformation of the tube. Unless checked by some additional confining influence (e.g., twisted flux tube fields), this tendency might ultimately lead to the destruction of the flux concentration. While quantitative treatment of this effect is beyond the scope of the present exploratory model calculations, its potential occurence underscores the desirability of performing 2D MHD simulations of flux tubes in a radiative enviroment.

We use the model consisting of equations (1)-(8) to determine the motion of an initially buoyant flux ring in the radiative envelope of a hot star. The distributions of pressure ($`p_e`$), density ($`\rho _e`$), and temperature ($`T_e`$) are taken from a spherical, non-rotating stellar interior model for $`M_{}=9M_{}`$, kindly calculated for us by S. Jackson. As described in Appendix B, in order to facilitate numerical computations, we have derived an approximate, analytic representation of this model by assuming that $`p_e`$, $`\rho _e`$, and $`T_e`$ are related polytropically. For the purpose of these simulations, the star is assumed to rotate rigidly with a prescribed angular velocity $`\mathrm{\Omega }`$ that is sufficiently slow that departures from sphericity can be neglected.

At time $`t=0`$, a toroidal flux ring with specified values of the cross-sectional radius $`a_0`$ and plasma beta $`\beta _0(=8\pi p_0/B_0^2)`$ is released from the outer boundary of the convective core ($`r_0=R_c`$) at co-latitude $`\theta _0`$ (latitude $`\lambda _0`$). We suppose that the ring is initially in thermal equilibrium with the ambient stellar material, $`T_0=T_{e0}`$. Application of the mechanical equilibrium condition (3) then yields $`p_0=p_{e0}\beta _0/(1+\beta _0)`$ and $`\rho _0=\rho _{e0}\beta _0/(1+\beta _0)`$, so that $`(\rho _{e0}\rho _0)>0`$ and the ring is buoyant. The initial field strength $`B_0`$ and Alfvén speed $`u_{A0}`$ in the ring are given in terms of $`\beta _0`$ and the external properties at the starting position according to $`B_0=[8\pi p_{e0}/(1+\beta _0)]^{1/2}`$ and $`u_{A0}=[2p_{e0}/(\beta _0\rho _{e0})]^{1/2}`$.

The velocity components ($`u_r,u_\theta ,u_\varphi `$) are obtained as functions of $`t`$ by numerical integration of the three components of equation (5); the position ($`r,\theta `$) of the axisymmetric ring follows by integrating the relations $`\mathrm{d}r/\mathrm{d}t=u_r`$, $`\mathrm{d}\theta /\mathrm{d}t=u_\theta /r`$ given in equation (4). In addition, we simultaneously integrate the thermal energy equation (7) to obtain $`S`$, which enables us to update the tube properties at each time step. In particular, from the definition given in equation (8), the pressure $`p`$ within the tube can be expressed as

$$p=p_0e^S\left(\frac{\rho }{\rho _0}\right)^\gamma .$$
(9)
Similarly, if equation (1) is rewritten in the form

$$\left(\frac{a}{a_0}\right)^2=\left(\frac{\rho _0}{\rho }\right)\left(\frac{r_0\mathrm{sin}\theta _0}{r\mathrm{sin}\theta }\right),$$
(10)
then substitution in equation (2) yields

$$B=B_0\left(\frac{\rho }{\rho _0}\right)\left(\frac{r\mathrm{sin}\theta }{r_0\mathrm{sin}\theta _0}\right).$$
(11)
This result, together with equation (9), permits the pressure equilibrium condition (3) to be recast as

$$e^S\left(\frac{\rho }{\rho _0}\right)^\gamma +\frac{1}{\beta _0}\left(\frac{r\mathrm{sin}\theta }{r_0\mathrm{sin}\theta _0}\right)^2\left(\frac{\rho }{\rho _0}\right)^2\left(1+\frac{1}{\beta _0}\right)\left(\frac{p_e}{p_{e0}}\right)=0,$$
(12)
which is solved to obtain $`\rho /\rho _0`$, given the values of $`S`$ and $`p_e`$ that prevail at a particular time $`t`$ and location $`(r,\theta )`$. Equations (9), (10), and (11) are then used to derive $`p`$, $`a`$, and $`B`$, and $`T`$ is determined via the relation $`T/T_0=(p/p_0)(\rho _0/\rho )`$.

## 3 Results

We have used the thin flux tube physical model and solution method outlined in the preceding section to investigate the dynamical behavior of magnetic flux rings with $`10^4\beta _010^7`$ and $`10^4(a_0/h_{e0})10^2`$, where $`h_{e0}`$ ($`=2.65\times 10^{10}`$ cm; see Appendix B) is the pressure scale height in the stellar interior at the starting radius $`r_0=R_c`$. The maximum and minimum values adopted for the quantity $`\beta _0`$ correspond to toroidal fields with strengths $`B_0=1.82\times 10^5`$ G and $`5.77\times 10^6`$ G, respectively. In the computational results described below, the parameters $`\gamma `$, $`C_D`$, and $`\mathrm{\Omega }`$ were assigned the values, $`\gamma =5/3`$, $`C_D=1`$, and $`\mathrm{\Omega }=5.88\times 10^5`$ s<sup>-1</sup>, the latter value representing a surface equatorial rotation speed of about 150 km s<sup>-1</sup>. We have examined the motion of buoyant flux rings following release from the core boundary at initial latitudes in the range $`0^{}\lambda _045^{}`$ (i.e., $`\pi /4\theta _0\pi /2`$). Because the dynamics of rings in the equatorial plane ($`\lambda _0=0^{}`$) differ from those of rings with $`\lambda _0>0^{}`$, we treat the two cases separately in the ensuing discussion.

### 3.1 Flux Rings in the Equatorial Plane

Inspection of the equation of motion (5) reveals that the meridional component of each force acting on a flux ring that is initially at rest ($`𝐮=0`$) in the equatorial plane ($`\theta =\pi /2`$) vanishes identically. The subsequent motion of such a ring therefore remains in this plane, and is governed by just the $`r`$\- and $`\varphi `$-components of equation (5). Some of the properties of configurations of this kind are illustrated in Figures 2 and 3, which contain results pertaining to rings with $`\beta _0=10^4`$. The pairs of panels that make up the three rows of Figure 2 show, respectively, the time evolution of: ($`i`$) the radial position of the ring, $`(\mathrm{\Delta }r/r_0)=(rr_0)/r_0`$; ($`ii`$) the velocity components $`(u_r,u_\theta ,u_\varphi )`$ in units of the Alfvén speed $`u_{A0}(=7.49\times 10^5`$ cm s<sup>-1</sup>) in the ring at $`t=0`$; and ($`iii`$), the accelerations produced by each of the forces acting on the ring, in units of $`u_{A0}^2/R_{}`$. In each panel, time is measured in units of the Alfvén time $`t_A`$, defined as $`t_AR_{}/u_{A0}(=3.42\times 10^5`$ s $`4`$ days).

The dynamical evolution of a flux ring with initial cross-sectional radius $`a_0=10^4h_{e0}`$ is shown in panels (A), (C), and (E) of Figure 2. Note that in this case, the outward expansion of the ring ceases after $`tt_A`$, its radial position remaining fixed at later times. In the course of evolving toward this apparent equilibrium configuration, the ring experiences a brief, initial period of accelerated motion in the $`+𝐞_r`$-direction, followed by radial deceleration to a state of rest in which $`u_r=0`$ but $`u_\varphi 0`$. Examination of the ring force balance indicates that the initial expansion is driven (as expected) by buoyancy, while being opposed by the slightly smaller, inward-directed, magnetic tension force. The aerodynamic drag force, although small in magnitude, contributes to the damping of the oscillatory motions that are induced by the interplay between the larger buoyant and magnetic forces.

The flux ring is assumed to corotate with the stellar interior at $`t=0`$, so that its initial azimuthal velocity is $`u_{\varphi 0}=0`$ in the reference frame that rotates with angular velocity $`\mathrm{\Omega }`$. For $`t>0`$, the azimuthal motion of the ring is such that the angular momentum of the material contained within it is conserved. The buoyancy-driven increase in the radial position of the ring is therefore accompanied by a decrease in the rotation rate of the ring material. Hence, in the frame of reference that corotates with the star, the azimuthal velocity of the ring is in the $`𝐞_\varphi `$-direction, as can be seen in panel (C). According to equation (6$`a`$), this azimuthal motion produces an inward-directed Coriolis force that grows with increasing $`\mathrm{\Delta }r`$ until its magnitude, combined with that of the magnetic tension force, is sufficient to balance the buoyant force and prevent further radial expansion of the ring. A similar behavior has been observed in connection with models of equatorial magnetic flux rings located in the radiative layers just below the base of the solar convection zone (Moreno-Insertis, Schüssler, and Ferriz-Mas 1992; Ferriz-Mas 1996).

In panels (A) and (C) of Figure 3, we show how the thermodynamic state of a ring with $`\beta _0=10^4`$ and $`a_0=10^4h_{e0}`$ evolves over time. Panel (A) shows the normalized differences between the pressure, density, and temperature in the tube and in the material surrounding it (i.e., $`\mathrm{\Delta }p/p_e`$, $`\mathrm{\Delta }\rho /\rho _e`$, and $`\mathrm{\Delta }T/T_e`$, respectively, where for any quantity $`Q`$, $`\mathrm{\Delta }QQ_eQ`$) as functions of time. Beginning from a state in which $`T_0=T_{e0}`$, the ring cools slightly relative to its surroundings during the earliest and most rapid portion of its limited buoyant ascent. As its outward progress stalls, radiative heating acts to restore the condition of thermal equilibrium between the ring and its external environment. In panel (C), we display the time evolution of the quantities $`S`$ and $`q`$ (see equations  and ), along with the corresponding histories of the instantaneous heating and rise times, defined as $`t_Hq^1`$ (see Appendix A) and $`t_R|h_e/u_r|`$, respectively. Note that $`S`$ increases as the tube is radiatively heated, and approaches a constant value as $`u_r0`$, $`TT_e`$, and both $`t_H`$ and $`t_R`$ become large.

To further explore the nature of the equilibrium seen in Figures 2 and 3 for flux rings in the equatorial plane, we write the $`\varphi `$-component of the equation of motion (5) in the form

$$r\frac{\mathrm{d}u_\varphi }{\mathrm{d}t}+u_ru_\varphi +2u_r\mathrm{\Omega }r=\frac{\mathrm{d}}{\mathrm{d}t}\left(ru_\varphi +\mathrm{\Omega }r^2\right)=0,$$
(13)
where we have assumed $`u_\theta =0`$ for $`\theta =\pi /2`$, and have used the fact that $`u_r=\mathrm{d}r/\mathrm{d}t`$. Equation (13) implies that the specific angular momentum $`(ru_\varphi +\mathrm{\Omega }r^2)`$ has a constant value which, after recalling the initial condition $`u_{\varphi 0}=0`$, is readily established to be $`\mathrm{\Omega }r_0^2`$. This quantity can then be used to evaluate the azimuthal velocity component of the ring, yielding the result

$$u_\varphi =\left(\frac{\mathrm{\Omega }r_0^2}{r}\right)\mathrm{\Omega }r.$$
(14)
Similar considerations applied to the $`r`$-component of equation (5) lead to the equilibrium condition

$$\left(\frac{\rho _e\rho }{\rho _e+\rho }\right)\left[\frac{GM(r)}{r^2}\mathrm{\Omega }^2r\right]+\frac{u_\varphi ^2}{r}+2\mathrm{\Omega }u_\varphi \frac{B^2}{4\pi r(\rho _e+\rho )}=0,$$
(15)
where we have set $`u_r=(\mathrm{d}u_r/\mathrm{d}t)=0`$. Equation (15) expresses the fact that the radial equilibrium is characterized by a balance between the outward-directed buoyancy and centrifugal forces and the inward-directed Coriolis and magnetic tension forces. Making the substitution $`(r/r_0)=1+(\mathrm{\Delta }r/r_0)`$ and retaining only the lowest order terms in $`(\mathrm{\Delta }r/r_0)`$ ($`1`$), we obtain

$$\frac{\mathrm{\Delta }r}{r_0}\frac{1}{8(\mathrm{\Omega }r_0)^2}\left(\frac{GM_c}{\beta _0r_0}u_{A0}^2\right)=\frac{1}{8\beta _0(\mathrm{\Omega }r_0)^2}\left(\frac{GM_c}{r_0}\frac{2p_{e0}}{\rho _{e0}}\right),$$
(16)
where equation (14) and the approximation $`\rho _e\rho [1+(1/\beta _0)]`$ have been used in simplifying equation (15).

As is evident from equation (16), for a given stellar model, the equilibrium position of an equatorial flux ring depends only on the value of the parameter $`\beta _0`$; in the case of the present model, $`(\mathrm{\Delta }r/r_0)6\times 10^4(\beta _0/10^4)^1`$. This dependence has been confirmed by examination of solutions corresponding to a variety of input parameter values, including the particular example shown in panels (B), (D), and (F) of Figure 2, and in panels (B) and (D) of Figure 3. For this solution, $`\beta _0=10^4`$ as in the case of the ring considered above, but $`a_0=10^3h_{e0}`$. Because this initial cross-sectional radius is a factor of 10 larger than the previous value, the radiative heating rate $`q`$ and the acceleration produced by the aerodynamic drag force are smaller by factors of 100 and 10, respectively (see equations  and ). As can be seen in the relevant portions of Figures 2 and 3, a consequence of these reductions is that the ring is subject to vigorous buoyancy oscillations in the radial direction. Oscillations of a similar kind in toroidal flux tubes inside non-rotating stars have been studied by Spruit & van Ballegooijen (1982), and inside rotating stars by van Ballegooijen (1983), Moreno-Insertis, Schüssler, & Ferriz-Mas (1992), and Ferriz-Mas & Schüssler (1993). In the present example, the enhanced tendency toward oscillatory motion is directly attributable to modifications in the buoyant acceleration of the ring; these changes are themselves a result of the altered thermodynamic state of the ring material, produced by the diminished rate of radiative heating. The build-up in the magnitude of the Coriolis force, together with the action of the drag force, cause the amplitude of these oscillations to decrease over time, leaving the ring in the same equilibrium configuration that obtained in the case of the $`a_0=10^4h_{e0}`$ ring studied above.

These results suggest that the equilibrium of a flux ring located in the equatorial plane is stable with respect to small amplitude, radial displacements from the position given by equation (16). This conjecture has been verified by a series of numerical experiments in which an additional, radially directed force of specified amplitude is applied to the ring for a brief time interval after it has assumed its equilibrium position. In all cases, the ring is initially displaced in radius, but is quickly restored to its original location following cessation of the imposed forcing. Such behavior is not observed in response to the application of latitudinal forcing of the same type. In all of these cases, a small displacement out of the equatorial plane is accompanied by the development of a $`\theta `$-component of the magnetic tension force that accelerates the ring toward the pole. This is similar to the dynamical evolution of the so-called poleward-slip instability of equatorial flux rings in the radiative layers beneath the solar convection zone (see, e.g., Moreno-Insertis, Schüssler, & Ferriz-Mas 1992, and references therein). For the conditions considered in the present paper, the perturbed force balance is such that a ring initially contained within the equatorial plane moves over time to higher latitudes and larger radii. Once out of the equatorial plane, the dynamics of such a ring is found to be identical to that of one with $`\theta _0<\pi /2`$, and is discussed in §3.2 below.

### 3.2 Flux Rings Outside of the Equatorial Plane

We now consider the motion of toroidal magnetic flux tubes having initial positions $`r_0=R_c`$ with $`0<\theta _0<\pi /2`$. We focus on the properties of rings with $`\theta _0=\pi /3`$ (i.e., initial latitude $`\lambda _0=30^{}`$) since the behavior in this case is representative of that displayed by all rings having starting locations along the core-envelope interface, out of the equatorial plane. In Figure 4, we show the trajectory that a ring with $`\beta _0=10^4`$ and $`a_0=10^4h_{e0}`$ follows during the time interval between its release and $`t=10^3t_A10.85`$ years. As can be seen in panel (A), the edge of the ring traces a path that is nearly parallel to the stellar rotation axis over this period, extending from $`r_0=0.232R_{}`$ and $`\lambda _0=30^{}`$ to $`r=0.64R_{}`$ and $`\lambda =72^{}`$. In panel (B), we have magnified the horizontal scale in order to make the meridional motion of the ring more apparent. There it can be seen that for a time $`t5t_A`$ after the start of the calculation, the distance $`r\mathrm{sin}\theta `$ of the ring from the stellar rotation axis decreases as the ring moves toward higher latitudes along the periphery of the convective core. At later times, this moment arm length approaches its initial size, while the height $`r\mathrm{cos}\theta `$ of the ring above the equatorial plane steadily increases.

Additional information pertaining to the evolving dynamical and thermodynamical properties of the flux ring under consideration is presented in Figure 5. Inspection of panels (A) and (B) reveals that after a brief, initial period of dynamical readjustment, the ring moves continuously in the direction of higher latitudes (i.e., the $`𝐞_\theta `$-direction) and larger radii. This is unlike the behavior of rings with $`\lambda _0=0^{}`$ which, in the absence of suitable perturbations, do not depart from the equatorial plane as they evolve toward a final equilibrium state with $`u_r=u_\theta =0`$ (see §3.1). In the present case, although its upward progress slows considerably at later times, the ring never attains an equilibrium in which the forces acting in the $`𝐞_r`$\- and $`𝐞_\theta `$-directions balance to produce a state of rest.

Details concerning the radial and meridional dynamics of the flux ring are provided in panels (C) and (D) of Figure 5. When the ring is released at $`t=0`$, it experiences a net radial acceleration that is $`>0`$, resulting from the fact that the outward, buoyant force is larger than the inward, radial component of the tension force. In the $`𝐞_\theta `$-direction, the tension force is initially unopposed by any other force components, leading to a net meridional acceleration that is $`<0`$. The ring is therefore pulled toward the pole, and the horizontal distance $`r\mathrm{sin}\theta `$ between it and the rotation axis decreases (see Figure 4). This change in position has an important consequence for the subsequent motion of the ring. In order to conserve angular momentum, the azimuthal velocity $`u_\varphi `$ increases as the moment arm of the ring material becomes smaller (see panel \[B\]). Associated with this spin-up are substantial increases in the magnitudes of the $`r`$\- and $`\theta `$-components of the Coriolis force, both of which exert considerable influence over the dynamical evolution of the ring.

Using a procedure similar to that adopted in the case of an equatorial flux ring, the azimuthal component of the equation of motion (5) for a ring at higher latitudes can, after some manipulation, be written as

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(r\mathrm{sin}\theta u_\varphi +\mathrm{\Omega }r^2\mathrm{sin}^2\theta \right)=0,$$
(17)
from which it immediately follows that

$$u_\varphi =\left(\frac{\mathrm{\Omega }r_0^2\mathrm{sin}^2\theta _0}{r\mathrm{sin}\theta }\right)\mathrm{\Omega }r\mathrm{sin}\theta .$$
(18)
Note that in this case, because the initial radial and meridional motion of the ring corresponds to decreasing $`r\mathrm{sin}\theta `$, the sense of the azimuthal motion is such that $`u_\varphi >0`$. Hence, the accelerations produced by the Coriolis force in the $`𝐞_r`$\- and $`𝐞_\theta `$-directions, $`2\mathrm{\Omega }u_\varphi \mathrm{sin}\theta `$ and $`2\mathrm{\Omega }u_\varphi \mathrm{cos}\theta `$, respectively, are both $`>0`$. As a result, the tension-induced, poleward movement of the ring at early times begins to slow when the sum of the $`\theta `$-components of the Coriolis and drag forces becomes large enough to change the sign of the net meridional acceleration, making it $`>0`$. Likewise, the acceleration supplied by buoyancy and the radial Coriolis force component causes $`u_r`$ to grow until $`u_r=u_\theta \mathrm{cot}\theta `$, at which point contraction of the ring toward the rotation axis ceases; at later times, $`u_r>u_\theta \mathrm{cot}\theta `$, and the ring expands. The resulting increase in $`r\mathrm{sin}\theta `$ leads to a gradual reduction in $`u_\varphi `$, and leaves the ring in a state of near balance between the Coriolis and tension forces in both the $`𝐞_r`$\- and $`𝐞_\theta `$-directions.

The lack of equilibria for flux rings that originate at latitudes $`\lambda _0>0^{}`$ can be understood by examining the $`r`$\- and $`\theta `$-components of the equation of motion (5). Assuming that a stationary state with $`u_r=u_\theta =0`$ exists, the $`r`$-component can be written in a form reminiscent of the equilibrium condition (15) for flux rings in the equatorial plane,

$$\left(\frac{\rho _e\rho }{\rho _e+\rho }\right)\left[\frac{GM(r)}{r^2}\mathrm{\Omega }^2r\mathrm{sin}^2\theta \right]+\frac{u_\varphi ^2}{r}+2\mathrm{\Omega }u_\varphi \mathrm{sin}\theta \frac{B^2}{4\pi r(\rho _e+\rho )}=0,$$
(19)
while the $`\theta `$-component becomes

$$\left(\frac{\rho _e\rho }{\rho _e+\rho }\right)\mathrm{\Omega }^2r\mathrm{sin}^2\theta +\frac{u_\varphi ^2}{r}+2\mathrm{\Omega }u_\varphi \mathrm{sin}\theta +\frac{B^2}{4\pi r(\rho _e+\rho )}=0.$$
(20)
Note that were it not for the appearance of the gravitational acceleration in the radial component, equations (19) and (20) would be identical. This correspondence implies that the conditions for force balance in the $`𝐞_r`$\- and $`𝐞_\theta `$-directions can be simultaneously satisfied only if $`\rho =\rho _e`$. However, as can be seen in panel (E) of Figure 5, $`\rho <\rho _e`$ throughout the time interval covered by the calculation, indicating that a stationary equilibrium is not possible in this particular case.

In general, all magnetic flux rings with $`\lambda _0>0^{}`$ remain slightly cooler and less dense than the surrounding material during most of their ascent, a consequence (in part) of the total pressure balance condition given by equation (3). On the basis of equations (19) and (20) then, stationary equilibria are ruled out for toroidal flux tubes of this kind. Because $`T<T_e`$ (see panel \[E\] of Figure 5), a given tube is heated by the radiation that diffuses into it from the external medium. In response to this heating and in order to comply with the pressure equilibrium condition, the tube expands, becoming less dense as it does so. Buoyancy then carries the tube outward to a somewhat larger radius. Hence, apart from a brief period at the start of the motion when the buoyant acceleration is established by the initial conditions, the rate of rise of the tube is controlled by the rate at which it is radiatively heated. This behavior is evident in panel (F) of Figure 5, wherein it can be seen that for much of the time interval depicted, the heating and rise times vary in concert. Such a situation is in contrast to the dynamical and thermodynamical evolution observed previously for flux rings in the equatorial plane. The stationary equilibrium of an equatorial flux ring is characterized by $`T=T_e`$, a vanishing radiative heating rate, and a balance between buoyancy and the sum of the Coriolis and tension forces (see Figures  and ).

The considerations of the preceding paragraph can be used to derive an approximate expression for the radial rise speed of a buoyant flux ring. Interpreting the derivative on the left-hand side of equation (7) as $`(\mathrm{d}S/\mathrm{d}t)=u_r(\mathrm{d}S/\mathrm{d}r)`$, the sought-after component of the tube velocity is

$$u_r=\frac{q}{(\mathrm{d}S/\mathrm{d}r)},$$
(21)
where $`q`$ and $`S`$ are defined in equations (7) and (8). To evaluate $`(\mathrm{d}S/\mathrm{d}r)`$, we assume that since $`\beta 1`$, $`pp_e`$ and $`\rho \rho _e`$; then,

$$S\mathrm{ln}\left[\frac{p_{e0}}{p_0}\left(\frac{\rho _0}{\rho _{e0}}\right)^\gamma \right]+\left(\frac{\alpha \gamma }{\alpha }\right)\mathrm{ln}\left(\frac{p_e}{p_{e0}}\right),$$
(22)
where we have made use of the polytrope relation in the form $`(\rho _e/\rho _{e0})=(p_e/p_{e0})^{1/\alpha }`$. Thus,

$$\frac{\mathrm{d}S}{\mathrm{d}r}\left(\frac{\alpha \gamma }{\alpha }\right)\frac{\mathrm{d}\mathrm{ln}p_e}{\mathrm{d}r}=\left(\frac{\gamma \alpha }{\alpha }\right)\frac{1}{h_e},$$
(23)
and

$$u_r\left(\frac{\alpha }{\gamma \alpha }\right)qh_e,$$
(24)
where $`h_e=p_e/(\rho _eg)`$ is the pressure scale height. In Figure 6, we compare the approximation given in equation (24) with the radial velocity component obtained by numerical solution of the full equation of motion for $`\beta _0=10^4`$, $`a_0=10^4h_{e0}`$, and $`\lambda _0=30^{}`$. At times later than $`tt_A`$, the two results are indistinguishable, indicating that the rate of ascent is indeed regulated by the radiative heating of the ring. In fact, using the definitions of the heating and rise time scales, equation (24) implies that $`t_Rt_H`$, in accordance with the behavior seen in panel (F) of Figure 5.

Figures 7 and 8 contain a summary of results from the simulated rise of a toroidal magnetic flux tube released at latitude $`\lambda _0=30^{}`$ with $`\beta _0=10^4`$ and $`a_0=10^3h_{e0}`$. As was discussed in §3.1, the larger cross-sectional radius of this tube implies that the deceleration due to aerodynamic drag ($`a^1`$) and the radiative heating rate ($`a^2`$) are both reduced from the magnitudes these quantities have when $`a_0=10^4h_{e0}`$. Similar to the equatorial flux rings that were the subject of that section, in the present case, these reductions are responsible for engendering oscillatory motion of the tube during the early stages of its ascent toward the surface. Such behavior can be seen in Figure 7, in which (as in Figure 4) two views of the path taken by the rising tube are given. In panel (A), it is apparent that the larger tube is less buoyant; during the time interval $`\mathrm{\Delta }t=10^3t_A`$ between the first and last points on the trajectory, the radial position of the ring only increases from $`r_0=0.232R_{}`$ to $`r=0.275R_{}`$. On the expanded scale of panel (B), it is seen that for a time $`0.25t_A`$ following the start of the calculation, the ring moves steadily toward the pole, and thereafter executes latitudinal oscillations that decrease in amplitude as the distance of the ring from the equatorial plane increases.

Inspection of the relevant panels of Figure 8 indicates that as in the case of the $`a_0=10^4h_{e0}`$ tube considered earlier in this section, the initial poleward motion is driven by the magnetic tension force. However, in the present case, the tube acquires a higher meridional velocity and approaches closer to the rotation axis during this movement, a consequence of the fact that the deceleration due to drag is smaller. Also, the change in volume arising from the decrease in $`r\mathrm{sin}\theta `$, together with the diminshed efficiency of radiative heating, are sufficient to ultimately make $`\mathrm{\Delta }\rho <0`$, and to thereby change the sign of the buoyant acceleration (see panels \[C\] and \[E\]). Note that the the spin-up associated with the overall contraction of the ring causes both components of the Coriolis acceleration to increase substantially; the oscillatory behavior seen in Figure 7 results from the interplay between this acceleration and that produced by the magnetic tension force.

### 3.3 Rise Times

A primary focus of the present investigation is the buoyant transport of magnetic flux from the assumed site of its generation in the convective core to the stellar surface. Of particular interest is an estimate of the time required for this transport to take place. From the definition of $`u_r`$ given in equation (4), it follows that

$$t_r(r)=\underset{r_0}{\overset{r}{}}\frac{\mathrm{d}r}{u_r},$$
(25)
is the time it takes a flux ring to travel outward from $`r_0`$ to any radius $`r>r_0`$.

We have computed $`t_r(r)`$ for toroidal flux tubes having the properties listed in Table 1. For each case, the radial component of the tube velocity was determined by numerical integration of the full equation of motion (5) for a time interval $`0tt_1`$, where the value of $`t_1`$ for a given solution was chosen to be between $`10^3t_A`$ and $`3\times 10^4t_A`$, depending upon the magnitude of $`u_r`$ for that solution. For $`t>t_1`$, $`u_r`$ was evaluated using an analytic approximation to equation (24), derived in the following way.

According to equation (24), $`u_rqh_e`$ at times late enough that any transient behavior stemming from the adjustment of the tube to the imposed initial conditions has disappeared. The heating rate $`q`$ depends directly on the temperature difference $`\mathrm{\Delta }T`$ between the tube and its surroundings; using the pressure equilibrium condition (3), it can be shown that $`(\mathrm{\Delta }T/T_e)B^2/(8\pi p_e)`$, assuming $`\rho \rho _e`$. Inserting this result in the definition of $`q`$ (see Appendix A), we find that the dependence of $`u_r`$ on the physical properties of the tube and the external medium is given by (see also Gurm & Wentzel 1967; Parker 1975)

$$u_r\frac{32}{3}\left[\frac{\alpha (\gamma 1)}{\gamma \alpha }\right]\left(\frac{\sigma T_e^4}{\kappa _e\rho _e^2g}\right)\left(\frac{B^2}{8\pi pa^2}\right).$$
(26)
To further simplify this expression, note that as evidenced by Figures 4 and 7, much of the ascent toward the surface takes place with $`r\mathrm{sin}\theta `$ constant (i.e., $`u_ru_\theta \mathrm{cot}\theta `$). Therefore, from equations (1) and (2), $`\rho (\rho _e)a^2`$ and $`Ba^2\rho _e`$, so that $`u_r`$ varies as

$$u_r\frac{T_e^3}{\kappa _eg},$$
(27)
where we have approximated $`pa^2TT_e`$. To illustrate the validity of equation (27), we have evaluated the constant of proportionality using the properties of solution 1 (see Table 1) at the time $`t=50t_A`$. The rise speed so-derived is shown as a function of time for $`t50t_A`$ by the dotted line in Figure 6.

The results of our rise time computations are depicted in graphical form in Figure 9. There we show the time $`t_r(r)`$ required to reach a given radius $`r`$ within the interval $`0.232R_{}r0.960R_{}`$, for each of the solutions listed in Table 1. It is apparent that rings with smaller values of $`\beta _0`$ and $`a_0`$ generally traverse the radiative envelope in less time than do rings that are less strongly magnetized and have larger cross-sectional radii (see also Moss 1989). Based on the discussion of the preceding subsection, rings having smaller values of $`\beta _0`$ and $`a_0`$ are initially more buoyant, and are more likely to remain so because of the comparatively shorter time it takes to heat them by inward diffusion of radiation. With the exception of the rings corresponding to solutions 7, 9, and 10, the time needed to arrive at the stellar surface is less than the estimated main sequence lifetime of the 9 $`M_{}`$ model, $`t_{MS}10^{10}(M_{}/M_{})(L_{}/L_{})2.4\times 10^7`$ years. However, note that even the more slowly rising tubes of solutions 7, 9, and 10 attain radii in the range (0.92-0.95) $`R_{}`$ within a time $`<t_{MS}`$. In fact, all of the solutions shown in Figure 9 spend more time crossing the outer 10% of the stellar radius than they do traveling from $`r=r_0`$ to $`r=0.9R_{}`$. From equation (27), it can be seen that since $`T_e`$ decreases and $`\kappa _e`$ increases as $`r`$ approaches $`R_{}`$ in the outer envelope, $`u_r`$ becomes small and the rise time increases accordingly.

Examination of Table 1 and Figure 9 also reveals that the omission of virtual inertia and aerodynamic drag from the equation of motion has little impact on ring dynamics and rise times (see solutions 1, 3, and 4). This is because the flux rings considered herein do not experience significant, impulsive accelerations, and never achieve large rise speeds for prolonged periods of time. Instead, most of the ascent takes place at the slow, quasi-steady rate given by equation (24) (or the approximation ). This behavior accounts for another notable characteristic of the results shown in Figure 9, namely, that the rise time profiles of several solutions are nearly identical to one another. From equation (26), it is readily seen that the dependence of the rise speed on the intrinsic properties of a given flux ring is $`u_r(a^2\beta )^1`$. Consideration of the entries in Table 1 indicates that, as expected, solutions having a common profile of $`t_r(r)`$ in Figure 9 are characterized by the same value of the parameter combination $`(a_0^2\beta _0)^1`$.

Finally, we note that with the preceding developments, it is possible to derive an approximate expression for the rise time $`t_r(r)`$. To do this, we first use the equation of hydrostatic equilibrium in the form $`\mathrm{d}\rho _e/\mathrm{d}r=(g\rho _e)/c_e^2`$, where $`c_e`$ is the external sound speed, to convert the integration variable in equation (25) from $`r`$ to $`\rho _e`$; this procedure yields

$$t_r(r)=\underset{\rho _{e0}}{\overset{\rho _e}{}}d\rho _e\left(\frac{c_e^2}{\rho _egu_r}\right)\left(\frac{c_{e0}^2}{\rho _{e0}g_0u_{r0}}\right)\underset{\rho _{e0}}{\overset{\rho _e}{}}d\rho _e\left(\frac{\rho _e}{\rho _{e0}}\right)^{(\alpha 2)}\left(\frac{T_{e0}}{T_e}\right)^3\left(\frac{\kappa _e}{\kappa _{e0}}\right).$$
(28)
To obtain the second, approximate equality in equation (26), we have replaced $`u_r`$ by the scaling relation (27) and used the fact that $`c_e^2=c_{e0}^2(\rho _e/\rho _{e0})^{(\alpha 1)}`$ where $`c_{e0}^2=(\alpha p_{e0}/\rho _{e0})`$. The integrand in equation (26) can be simplified by assuming that the opacity follows Kramer’s law, $`\kappa _e\rho _eT_e^{3.5}`$. This assumption, together with the polytrope relation $`T_e\rho _e^{(\alpha 1)}`$ (see Appendix B), allows us to express the integrand as a function of $`\rho _e`$ only; performing the integration we obtain

$$t_r(r)=\left(\frac{2}{11\alpha 13}\right)\left(\frac{c_{e0}^2}{g_0u_{r0}}\right)\left[\left(\frac{\rho _e}{\rho _{e0}}\right)^{\frac{(11\alpha 13)}{2}}1\right].$$
(29)
If $`u_{r0}`$ is determined by evaluating equation (26) at $`r_0`$, the rise time estimate (27) becomes

$$t_r(r)=3.0835\times 10^3\left[\left(\frac{a_0}{h_{e0}}\right)^2\beta _0\right]\left[\left(\frac{\rho _e}{\rho _{e0}}\right)^{0.9349}1\right]\mathrm{years},$$
(30)
where the interior model of Appendix B has been used to fix the values of all physical quantities at the core-envelope interface. As an example, the dashed line in Figure 9 represents the rise time profile obtained from equation (30) for $`[(a_0/h_{e0})^2\beta _0]=10^2`$

## 4 Conclusions and Discussion

The foregoing examination of the physics of magnetic flux rings in the radiative envelope of a uniformly rotating hot star has revealed a range of possible dynamical behaviors, and has enabled us to estimate the efficiency of buoyancy as a means of transporting dynamo-generated fields from the core to the surface. Rings located in the equatorial plane can attain a stationary equilibrium in which $`u_r=0,u_\varphi 0`$, and the combined Coriolis and tension forces balance buoyancy. This state is stable against infinitesimal displacements of the ring in the radial direction, but unstable when the perturbations are meridionally directed, causing the ring position to shift from latitude $`0^{}`$. Apart from transient episodes of oscillatory behavior at the outset of their motion, rings released from latitudes other than $`0^{}`$ move to larger radii along paths that are very nearly parallel to the stellar rotation axis. During most of the ascent, the rate of rise of such a ring is controlled by the rate at which it is radiatively heated, with a near balance prevailing between the accelerations produced by the Coriolis and magnetic tension forces in the $`r`$\- and $`\theta `$-directions.

Strongly magnetized rings with smaller cross-sectional radii experience larger initial buoyant accelerations and have shorter radiative heating time scales. As a result, they attain higher rise speeds and require less time to traverse the envelope than do rings with weaker fields and/or larger cross-sectional radii (see also Moss 1989). A majority of the ring models enumerated in Table 1 reach the stellar surface in a time that is less than the main sequence lifetime of the 9 $`M_{}`$ star, the only exceptions being solutions with large values of $`\beta _0`$ and/or $`a_0`$. If the field strengths and flux tube sizes considered herein are realistic, this result suggests that the field produced by a core dynamo could manifest itself in the surface layers of a hot star at a relatively early stage of main sequence evolution. In light of the discussion of §3, the strength of the field in a buoyant flux ring that reaches the photosphere of the 9 $`M_{}`$ stellar model is estimated to be $`B_{}B_0(\rho _e/\rho _{e0})9.44(\beta _0/10^4)^{1/2}`$ G.

The preceding conclusions are drawn from results obtained using a simple model for a thin, isolated, untwisted flux ring immersed in a non-evolving, rigidly rotating stellar interior. Given the exploratory character of this investigation, a more detailed treatment of flux tube structure and dynamics is neither warranted nor practical. However, several potentially important effects that have been neglected in the present analysis can be included within the context of the basic model of §2. Among these are large-scale internal circulatory flows, the development (over time) of gradients in composition and angular velocity, the influence of radiation pressure on tubes inside more massive stars, mass loss, and the consequences of magnetic interference with radiative energy transport (e.g., thermal shadows and small-scale flows; Parker 1984).

Of the effects noted above, rotationally driven, meridional circulation may have the most significant impact on the results described in §3. As noted in that section, each of the ring models listed in Table 1 spends considerably more time traveling through the thin shell $`0.90R_{}r0.96R_{}`$ than in ascending from the starting radius to $`r=0.90R_{}`$ (see Figure 9). In solution 1, for example, the former displacement requires approximately $`1.29\times 10^5`$ years, while the latter occurs in only $`1.38\times 10^3`$ years, a difference in travel time of nearly a factor of 100. We point out that the inclusion of meridional circulation may decrease the time needed to transport magnetic flux through the layers just below the photosphere. Estimates of the circulation speed inside upper main sequence stars suggest that while such flows proceed quite slowly at great depth, they can become fast in the lower-density region near the surface (see, e.g., Kippenhahn & Weigert 1994). Hence, flux tube advection by a circulatory flow may act to supplement buoyant transport in the outermost portion of the stellar interior, thereby shortening the time between the production of fields in the core and their emergence in the atmosphere. The combined effects of buoyancy and advection will be investigated in the next paper of this series.

We are grateful to Yuhong Fan, M. Maheswaran and Paul Charbonneau for helpful conversations, and to Steven Jackson for providing the interior model used in the calculations, and to Wendy Mukluk for the rising flux tube figure. We thank the referee for a useful question regarding tube heating and for suggesting additional references. JPC is also grateful for a visiting scientist position at HAO.

## Appendix A  Derivation of the Ring Energy Equation

To derive the energy equation for a thin, axisymmetric magnetic flux ring, we begin with the first law of thermodynamics, written in the form,

$$\frac{\mathrm{d}}{\mathrm{d}t}+p\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{\rho }\right)=\frac{\mathrm{d}𝒬}{\mathrm{d}t},$$
(A1)
where $`=(p/\rho )/(\gamma 1)`$ is the specific internal energy of the gas, and $`(\mathrm{d}𝒬/\mathrm{d}t)`$ is the rate at which heat is input, per unit mass of ring material. Let $`F_e`$ be the radiative flux at the surface of the ring; in the case of interest to us here, $`T<T_e`$, and this flux is directed into the ring from the external medium. We assume that the ring is sufficiently opaque that all of the flux incident upon it is thermalized within it. For a ring of mass $`M`$, the specific radiative heating rate is then

$$\frac{\mathrm{d}𝒬}{\mathrm{d}t}=\frac{AF_e}{M}=\frac{2F_e}{\rho a},$$
(A2)
where $`A=2\pi a2\pi r\mathrm{sin}\theta `$ is the surface area of the ring. If $`T_e`$ is the local temperature gradient in the surrounding stellar material and $`𝐧`$ is a unit vector normal to the ring surface, the diffusive flux of radiation into the ring is approximately

$$F_e=\frac{16}{3}\frac{\sigma T_e^3}{\kappa _e\rho _e}(𝐧T_e)\frac{16}{3}\frac{\sigma T_e^3}{\kappa _e\rho _e}\frac{\mathrm{\Delta }T}{a},$$
(A3)
where $`\mathrm{\Delta }T=(T_eT)`$, and we have estimated $`|T_e|(\mathrm{\Delta }T/a)`$.

Explicit calculation of the derivatives appearing on the left-hand side of equation (A1) leads to

$$\frac{\mathrm{d}}{\mathrm{d}t}+p\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{\rho }\right)=\frac{p}{(\gamma 1)\rho }\frac{\mathrm{d}S}{\mathrm{d}t},$$
(A4)
where the quantity

$$S=\mathrm{ln}\left[\left(\frac{p}{p_0}\right)\left(\frac{\rho _0}{\rho }\right)^\gamma \right],$$
(A5)
has been referenced to the initial thermodynamic conditions in the ring. Substitution of the results given in equations (A2), (A3), and (A4) into equation (A1) yields the ring energy equation in the form used herein,

$$\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{32}{3}\frac{(\gamma 1)\sigma T_e^4}{\kappa _e\rho _epa^2}\frac{\mathrm{\Delta }T}{T_e}q.$$
(A6)
The time scale $`t_H`$ for radiative heating of the ring can be expressed in terms of the rate $`q`$ as $`t_H[/(\mathrm{d}𝒬/\mathrm{d}t)]=q^1`$.

## Appendix B  Approximate Stellar Structure

The distributions of $`p_e`$, $`\rho _e`$, and $`T_e`$ used in the calculations described in §3 are derived from a model for a spherical, non-rotating star with mass $`M_{}=9M_{}`$ and radius $`R_{}=3.678R_{}`$. The star is chemically homogeneous with composition $`X=0.7131,Z=0.0169`$, and contains a central, convective core of mass $`M_c=0.308M_{}`$ and radius $`R_c=0.232R_{}`$. The ambient physical conditions at the initial radial position $`r_0(=R_c)`$ of the flux ring are $`p_{e0}=1.3250\times 10^{16}`$ dyne cm<sup>-2</sup>, $`\rho _{e0}=4.7239`$ g cm<sup>-3</sup>, and $`T_{e0}=1.9976\times 10^7`$ K. In order to enhance the precision of computations involving ring properties that differ from those of the external medium by very small amounts, it is helpful to use an analytic approximation to the background stellar model. We discuss the derivation of one such representation in this appendix.

Unlike the convective envelope of the Sun, the radiative envelope of the adopted stellar model contains a significant fraction of the mass of the star, $`1(M_c/M_{})0.70`$. Consequently, the gravitational acceleration $`g(r)`$ in the inner portion of the envelope does not vary simply as $`r^2`$, but instead depends on the detailed distribution of mass. We approximate $`g(r)`$ within the region between the convective core and the stellar surface by using the following piecewise continuous fit to the computed interior model:

$$\begin{array}{cccc}g(z)& =& g_1\mathrm{exp}\left[\left(\frac{z\xi _1}{\mathrm{\Delta }_1}\right)^2\right],& z_0zz_1,\\ & =& g_2\mathrm{exp}\left[\left(\frac{z\xi _2}{\mathrm{\Delta }_2}\right)\right],& z_1zz_2,\\ & =& \frac{g_3}{z^2},& z_2z1,\end{array}$$
(B1)
where $`zr/R_{}`$ and

$$g_1=1.06\times 10^5\mathrm{cm}\mathrm{s}^2,g_2=8.90\times 10^4\mathrm{cm}\mathrm{s}^2,g_3=1.80\times 10^4\mathrm{cm}\mathrm{s}^2,$$

$$\xi _1=0.27,\mathrm{\Delta }_1=0.31,\xi _2=0.40,\mathrm{\Delta }_2=0.34,$$

$$z_0=0.232,z_1=0.422,z_2=0.641.$$

We assume that $`p_e`$ and $`\rho _e`$ satisfy a polytrope relation of the form $`p_e\rho _e^\alpha `$ where the index $`\alpha `$ has a constant value. This assumption, together with the approximation for $`g`$ given in equation ($`B1`$), enables us to integrate the equation of hydrostatic equilibrium in order to obtain the structural properties of the stellar radiative interior. The results of this procedure are:

$$\frac{\rho _e}{\rho _{e0}}=\left\{1\left(\frac{\alpha 1}{\alpha }\right)\left(\frac{\rho _{e0}g_1R_{}}{p_{e0}}\right)\left(\frac{\sqrt{\pi }\mathrm{\Delta }_1}{2}\right)\left[\mathrm{\Phi }\left(\frac{z\xi _1}{\mathrm{\Delta }_1}\right)\mathrm{\Phi }\left(\frac{z_0\xi _1}{\mathrm{\Delta }_1}\right)\right]\right\}^{\frac{1}{(\alpha 1)}},$$
(B2)
in $`z_0zz_1`$, where $`\mathrm{\Phi }(x)=(2/\sqrt{\pi })\underset{0}{\overset{x}{}}dte^{t^2}`$ is the error function;

$$\frac{\rho _e}{\rho _{e1}}=\left\{1\left(\frac{\alpha 1}{\alpha }\right)\left(\frac{\rho _{e1}g_2R_{}}{p_{e1}}\right)\left[\mathrm{exp}\left(\frac{\xi _2z_1}{\mathrm{\Delta }_2}\right)\mathrm{exp}\left(\frac{\xi _2z}{\mathrm{\Delta }_2}\right)\right]\right\}^{\frac{1}{(\alpha 1)}},$$
(B3)
in $`z_1zz_2`$, where $`\rho _{e1}=\rho _e(z_1)`$ and $`p_{e1}=p_e(z_1)`$; and,

$$\frac{\rho _e}{\rho _{e2}}=\left[1\left(\frac{\alpha 1}{\alpha }\right)\left(\frac{\rho _{e2}g_3R_{}}{p_{e2}}\right)\left(\frac{1}{z_2}\frac{1}{z}\right)\right]^{\frac{1}{(\alpha 1)}},$$
(B4)
in $`z_2z1`$, where $`\rho _{e2}=\rho _e(z_2)`$ and $`p_{e2}=p_e(z_2)`$. By virtue of the polytrope relation, the pressure and temperature in each of the three intervals are given by

$$\frac{p_e}{p_{ei}}=\left(\frac{\rho _e}{\rho _{ei}}\right)^\alpha ,\frac{T_e}{T_{ei}}=\left(\frac{\rho _e}{\rho _{ei}}\right)^{\alpha 1},$$
(B5)
where $`i=0,1,2`$.

A polytrope with $`\alpha =1.3518`$ yields a reasonable fit to the properties of the computed stellar interior model. Note that a polytropic representation of this kind is not completely consistent with the detailed model, a consequence of the fact that the simple assumed relation between $`p_e`$ and $`\rho _e`$ provides only an average description of the entire stellar envelope. A more accurate fit can be obtained if the index $`\alpha `$ is allowed to have a distinct value in each of the three intervals delineated above (e.g., $`\alpha _1=1.400`$, $`\alpha _2=1.330`$, and $`\alpha _3=1.277`$). However, the single index description is adequate for the purpose of the present computations.


format: text

text:
# References

On the $`q`$-analogues of the Zassenhaus formula
for disentangling exponential operators

R. SridharE-mail: sridhar@imsc.res.in and R. JagannathanE-mail: jagan@imsc.res.in

The Institute of Mathematical Sciences
C.I.T. Campus, Tharamani, Chennai, TN 600113, India

To our dear friend Srinivasa Rao with admiration, affection, and best wishes
on his sixtieth birthday

Abstract

Katriel, Rasetti and Solomon introduced a $`q`$-analogue of the Zassenhaus formula written as $`e_q^{(A+B)}`$ $`=`$ $`e_q^Ae_q^Be_q^{c_2}e_q^{c_3}e_q^{c_4}e_q^{c_5}\mathrm{}`$, where $`A`$ and $`B`$ are two generally noncommuting operators and $`e_q^z`$ is the Jackson $`q`$-exponential, and derived the expressions for $`c_2`$, $`c_3`$ and $`c_4`$. It is shown that one can also write $`e_q^{(A+B)}`$ $`=`$ $`e_q^Ae_q^Be_{q^2}^{𝒞_2}e_{q^3}^{𝒞_3}e_{q^4}^{𝒞_4}e_{q^5}^{𝒞_5}\mathrm{}`$. Explicit expressions for $`𝒞_2`$, $`𝒞_3`$ and $`𝒞_4`$ are given.

1. Introduction

As is well known, the Baker-Campbell-Hausdorff (BCH) formula

$`e^Ae^B`$ $`=`$ $`e^{A+B+\frac{1}{2}[A,B]+\frac{1}{12}([A,[A,B]]+[[A,B],B])+\mathrm{}}`$ (1)
helps express the product of two noncommuting exponential operators as a single exponential operator in which the exponent is, in general, an infinite series in terms of repeated commutators and several hundred terms of the series have been calculated using the computer. The dual of the BCH formula is the Zassenhaus formula

$`e^{A+B}`$ $`=`$ $`e^Ae^Be^{\frac{1}{2}[A,B]}e^{\frac{1}{6}[A,[A,B]]\frac{1}{3}[[A,B],B]}\mathrm{}`$ (2)
which helps disentangle an exponential operator into a product of, in general, an infinite series of exponential operators involving repeated commutators (see  for details). The BCH and the Zassenhaus formulas have several applications (see, for example, (-)).

With the advent of $`q`$-deformed algebraic structures in physics (see, for example,  and references therein) there has been a growing interest in $`q`$-generalization of several results of classical analysis. Katriel and Solomon  first obtained the $`q`$-analogue of the BCH formula and later Katriel, Rasetti and Solomon  proposed a $`q`$-analogue of the Zassenhaus formula. Let us briefly recall their results.

The Jackson $`q`$-exponential is defined by

$`e_q^x`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{[n]!}}`$ (3)
where

$`[n]`$ $`=`$ $`{\displaystyle \frac{1q^n}{1q}}`$ (4)
is the Heine basic number and

$`[n]!`$ $`=`$ $`[n][n1][n2]\mathrm{}[1][0]!=1`$ (5)
such that

$`\underset{q1}{lim}[n]`$ $`=`$ $`n\underset{q1}{lim}e_q^x=e^x.`$ (6)
In the following $`[n]`$ will refer to the $`q`$-deformed $`n`$ defined by (4) corresponding to the base $`q`$. If the base is different then it will be indicated explicitly; for example, $`[n]_{q^k}`$ will mean $`(1q^{kn})/(1q^k)`$. The $`q`$-exponential function $`e_q^{\alpha x}`$ is the eigenfunction of the Jackson $`q`$-differential operator

$`D_qf(x)`$ $`=`$ $`{\displaystyle \frac{f(x)f(qx)}{(1q)x}}`$ (7)
such that

$`D_qe_q^{\alpha x}`$ $`=`$ $`\alpha e_q^{\alpha x}.`$ (8)
The $`q`$-commutator is defined by

$`[X,Y]_q`$ $`=`$ $`XYqYX`$ (9)
and it obeys the $`q`$-antisymmetry property

$`[Y,X]_q`$ $`=`$ $`q[X,Y]_{q^1}.`$ (10)
The $`q`$-BCH formula found in  reads

$`e_q^Ae_q^B`$ $`=`$ $`e_q^{A+B+\frac{q}{[2]}[A,B]_{q^1}+\frac{q^2}{[2]3!}\left([A,[A,B]_q]_{q^1}+[[A,B]_q,B]_{q^1}\right)+\mathrm{}}.`$ (11)
Katriel and Solomon  have given the explicit expressions for the terms involving up to 4-tuple $`q`$-commutator. In the limit $`q`$ $``$ $`1`$ the above $`q`$-BCH formula (11) is seen to agree with the classical BCH formula (1).

When $`A`$ and $`B`$ satisfy the relation $`AB`$ $`=`$ $`q^1BA`$ it is found that $`[A,B]_{q^1}`$, $`[A,[A,B]_q]_{q^1}`$, $`[[A,B]_q,B]_{q^1}`$, and all the higher $`q`$-commutators vanish leading to the result of Schützenberger  and Cigler  (see also, ) :

$`e_q^Ae_q^B`$ $`=`$ $`e_q^{A+B}\text{if}AB=q^1BA.`$ (12)
The $`q`$-analogue of the Zassenhaus formula proposed in  reads

$`e_q^{(A+B)}=e_q^Ae_q^Be_q^{c_2}e_q^{c_3}e_q^{c_4}e_q^{c_5}\mathrm{}`$ (13)
where

$`c_2`$ $`=`$ $`[B,A]_q/[2],`$
$`c_3`$ $`=`$ $`([[B,A]_q,A]_{q^2}/[3]!)+([[B,A]_q,B]_q/[3]),`$
$`c_4`$ $`=`$ $`([[[B,A]_q,A]_{q^2},A]_{q^3}/[4]!)+([[[B,A]_q,B]_{q^2},B]_{q^3}/[2][4])`$ (14)
$`+([[[B,A]_q,A]_{q^2},B]_q/[2][4])+(q[[B,A]_q,[B,A]_q]_q/[2]^2[4]),`$
$`\mathrm{}.`$
In this article we shall see that it is possible to have a $`q`$-Zassenhaus formula written also as

$`e_q^{(A+B)}=e_q^Ae_q^Be_{q^2}^{𝒞_2}e_{q^3}^{𝒞_3}e_{q^4}^{𝒞_4}e_{q^5}^{𝒞_5}\mathrm{}`$ (15)
and give the explicit expressions for the first few terms.

2. The Zassenhaus formula

A standard procedure to obtain the Zassenhaus formula is as follows. Set

$`e^{x(A+B)}`$ $`=`$ $`e^{xA}e^{xB}e^{x^2C_2}e^{x^3C_3}\mathrm{}.`$ (16)
Differentiating both sides of (16) with respect to $`x`$ and multiplying it from the right by

$`e^{x(A+B)}`$ $`=`$ $`\mathrm{}e^{x^3C_3}e^{x^2C_2}e^{xB}e^{xA}`$ (17)
one obtains

$`A+B`$ $`=`$ $`A+e^{xA}Be^{xA}+e^{xA}e^{xB}(2xC_2)e^{xB}e^{xA}`$ (18)
$`+e^{xA}e^{xB}e^{x^2C_2}(3x^2C_3)e^{x^2C_2}e^{xB}e^{xA}+\mathrm{}.`$
The expressions $`e^{xA}Be^{xA}`$, $`e^xe^{xB}(2xC_2)e^{xB}e^{xA}`$, etc., are expanded again using the formula

$`e^{xA}Be^{xA}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{n!}}\{A^n,B\}`$ (19)
where the multiple commutator bracket $`\{A^n,B\}`$ is defined by

$`\{A^{n+1},B\}`$ $`=`$ $`[A,\{A^n,B\}]\{A^0,B\}=B.`$ (20)
Then, equation (18) becomes

$`0`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{n!}}\{A^n,B\}`$ (21)
$`+2x{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^{m+n}}{m!n!}}\{A^m,B^n,C_2\}`$
$`+3x^2{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^{l+m+2n}}{l!m!n!}}\{A^l,B^m,C_2^n,C_3\}+\mathrm{}`$
with

$`\{A^0,B^n,C_2\}`$ $`=`$ $`\{B^n,C_2\}`$
$`\{A^{m+1},B^n,C_2\}`$ $`=`$ $`[A,\{A^m,B^n,C_2\}]\mathrm{}.`$ (22)
Equating the coefficients of $`x^n`$ to zero in (21) one obtains

$`C_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}[A,B]`$
$`C_3`$ $`=`$ $`{\displaystyle \frac{1}{6}}[A,[A,B]]{\displaystyle \frac{1}{3}}[[A,B],B]\mathrm{}.`$ (23)
There are also alternative methods of derivation of the above result (see, for example,  for more details). For our purpose of deriving a $`q`$-analogue of the Zassenhaus formula we shall follow Karplus and Schwinger  (see Appendix I where $`\mathrm{exp}(A+B)`$ is expanded in powers of $`B`$ up to the second term using a method we adopt here).

3. A $`q`$-analogue of the Zassenhaus formula

Let

$`F(x)`$ $`=`$ $`e_q^{x(A+B)}F(0)=I`$ (24)
and write

$`e_q^{x(A+B)}`$ $`=`$ $`e_q^{xA}G(x)G(0)=I.`$ (25)
On $`q`$-differentiating with respect to $`x`$

$`D_qF(x)`$ $`=`$ $`(A+B)e_q^{x(A+B)}.`$ (26)
Note that $`q`$-differentiation obeys the $`q`$-Leibniz rule

$`D_q(fg)`$ $`=`$ $`(D_qf)g+f(qx)(D_qg)=(D_qf)g(qx)+f(D_qg).`$ (27)
Now, $`q`$-differentiating the right hand side of (25) using the $`q`$-Leibniz rule we have

$`(A+B)e_q^{x(A+B)}`$ $`=`$ $`Ae_q^{xA}G(x)+e_q^{qxA}D_qG(x)`$ (28)
$`=`$ $`Ae_q^{x(A+B)}+e_q^{qxA}D_qG(x).`$
Thus,

$`Be_q^{x(A+B)}`$ $`=`$ $`Be_q^{xA}G(x)=e_q^{qxA}D_qG(x).`$ (29)
Or

$`D_qG(x)`$ $`=`$ $`(e_q^{qxA})^1Be_q^{xA}G(x).`$ (30)
Recall that the inverse of $`q`$-differentiation is $`q`$-integration (see, for example, ), defined by

$`{\displaystyle _0^x}d_q\xi f(\xi )`$ $`=`$ $`(q1)x{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}q^nf(q^nx)`$ (31)
such that

$`D_q\left({\displaystyle _0^x}d_q\xi f(\xi )\right)`$ $`=`$ $`f(x).`$ (32)
Then, the formal solution of the $`q`$-differential equation (30) for $`G(x)`$ with the initial condition $`G(0)`$ $`=`$ $`I`$, is given by

$`G(x)`$ $`=`$ $`I+{\displaystyle _0^x}d_q\xi (e_q^{q\xi A})^1Be_q^{\xi A}G(\xi )`$ (33)
$`=`$ $`I+{\displaystyle _0^x}d_q\xi (e_q^{q\xi A})^1Be_q^{\xi A}`$
$`+{\displaystyle _0^x}d_q\xi _1(e_q^{q\xi _1A})^1Be_q^{\xi _1A}{\displaystyle _0^{\xi _1}}d_q\xi _2(e_q^{q\xi _2A})^1Be_q^{\xi _2A}`$
$`+\mathrm{}.`$
Using the well known result

$`e_q^xe_{q^1}^x`$ $`=`$ $`1`$ (34)
we can rewrite (33) as

$`G(x)`$ $`=`$ $`I+{\displaystyle _0^x}d_q\xi e_{q^1}^{q\xi A}Be_q^{\xi A}`$ (35)
$`+{\displaystyle _0^x}d_q\xi _1e_{q^1}^{q\xi _1A}Be_q^{\xi _1A}{\displaystyle _0^{\xi _1}}d_q\xi _2e_{q^1}^{q\xi _2A}Be_q^{\xi _2A}`$
$`+\mathrm{}.`$
We now require a $`q`$-analogue of the classical formula (19). By straightforward expansion and regrouping of terms in powers of $`x`$ we have

$`(e_q^{qxA})^1Be_q^{xA}`$ $`=`$ $`e_{q^1}^{qxA}Be_q^{xA}`$ (36)
$`=`$ $`B+x[B,A]_q+{\displaystyle \frac{x^2}{[2]!}}[[B,A]_q,A]_{q^2}`$
$`+{\displaystyle \frac{x^3}{[3]!}}[[[B,A]_q,A]_{q^2},A]_{q^3}+\mathrm{}.`$
Let us write this equation as

$`e_{q^1}^{qxA}Be_q^{xA}`$ $`=`$ $`B+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{[n]!}}X_n`$ (37)
with the definition

$`X_n`$ $`=`$ $`[\mathrm{}[[[B,A]_q,A]_{q^2},A]_{q^3}\mathrm{}],A]_{q^n}n=1,2,3,\mathrm{}.`$ (38)
Now, using (37) and (38) in (35), we get

$`G(x)`$ $`=`$ $`I+{\displaystyle _0^x}d_q\xi \left\{B+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\xi ^n}{[n]!}}X_n\right\}`$ (39)
$`+{\displaystyle _0^x}d_q\xi _1\left\{B+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\xi _1^n}{[n]!}}X_n\right\}{\displaystyle _0^{\xi _1}}d_q\xi _2\left\{B+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\xi _2^n}{[n]!}}X_n\right\}`$
$`+\mathrm{}.`$
Collecting the first few terms of the resulting series in powers of $`x`$ we have

$`G(x)`$ $`=`$ $`\left\{I+xB+{\displaystyle \frac{x^2B^2}{[2]!}}+{\displaystyle \frac{x^3B^3}{[3]!}}+\mathrm{}\right\}`$ (40)
$`+\left\{{\displaystyle \frac{x^2X_1}{[2]!}}+{\displaystyle \frac{x^4X_1^2}{[1][2]![4]}}+{\displaystyle \frac{x^6X_1^3}{[1]^2[2]![4][6]}}+\mathrm{}\right\}`$
$`+\left\{x^3\left({\displaystyle \frac{X_2}{[3]!}}+{\displaystyle \frac{BX_1}{[3]!}}+{\displaystyle \frac{X_1B}{[1][3]}}\right)+x^6\left(\mathrm{}\right)+\mathrm{}\right\}`$
$`+\mathrm{}.`$
Realizing that the terms in the first bracket sum to $`e_q^{xB}`$ let us rewrite (40) as

$`G(x)`$ $`=`$ $`e_q^{xB}(I+e_{q^1}^{xB}\{{\displaystyle \frac{x^2X_1}{[2]!}}+{\displaystyle \frac{x^4X_1^2}{[1][2]![4]}}+{\displaystyle \frac{x^6X_1^3}{[1]^2[2]![4][6]}}+\mathrm{}\}`$
$`+e_{q^1}^{xB}\{x^3({\displaystyle \frac{X_2}{[3]!}}+{\displaystyle \frac{BX_1}{[3]!}}+{\displaystyle \frac{X_1B}{[1][3]}})+x^6\left(\mathrm{}\right)+\mathrm{}\}+\mathrm{}).`$
Substituting the series expansion $`e_{q^1}^{xB}`$, using the relation

$`[n]_{q^1}`$ $`=`$ $`q^{1n}[n]_q`$ (42)
and after some straightforward algebra one can rewrite (LABEL:gx2) further as

$`G(x)`$ $`=`$ $`e_q^{xB}e_{q^2}^{x^2X_1/[2]}(I+e_{q^2}^{x^2X_1/[2]}\{xB+{\displaystyle \frac{qx^2B^2}{[2]}}{\displaystyle \frac{q^2x^3B^3}{[3]}}+\mathrm{}\}\mathrm{}`$ (43)
$`+e_{q^2}^{x^2X_1/[2]}\{x^3({\displaystyle \frac{X_2}{[3]!}}+{\displaystyle \frac{BX_1}{[3]!}}+{\displaystyle \frac{X_1B}{[3]}})+\mathrm{}\}+\mathrm{}).`$
Now, substituting the series expression for $`e_{q^2}^{x^2X_1/[2]}`$ and simplifying, one recognizes that one can pull out a factor $`e_{q^3}^{x^3\{([X_1,B]_q/[3])+(X_2/[3]!)\}}`$ from the above expression and write

$`G(x)`$ $`=`$ $`e_q^{xB}e_{q^2}^{x^2[B,A]_q/[2]}e_{q^3}^{x^3\{([[B,A]_q,B]_q/[3])+([[B,A]_q,A]_{q^2}/[3]!)\}}\mathrm{}`$
$`=`$ $`e_q^{xB}e_{q^2}^{qx^2[A,B]_{q^1}/[2]}e_{q^3}^{x^3\{(q^3[A,[A,B]_{q^1}]_{q^2}/[3]!)(q[[A,B]_{q^1},B]_q/[3])\}}\mathrm{}.`$
This shows that the general expression for $`G(x)`$ can be assumed to be

$`G(x)`$ $`=`$ $`e_q^{xB}e_{q^2}^{x^2𝒞_2}e_{q^3}^{x^3𝒞_3}e_{q^4}^{x^4𝒞_4}\mathrm{}.`$
The crucial point here is to note that the exponential factors in $`G(x)`$ have bases $`q`$, $`q^2`$, $`q^3`$, $`q^4`$, $`\mathrm{}`$, respectively, unlike in the formula (13) introduced by Katriel, Rasetti and Solomon . Thus, having recognized a new general form of $`G(x)`$ we can use the comparison method to determine $`𝒞_2`$, $`𝒞_3`$, $`𝒞_4`$, $`\mathrm{}`$, in (LABEL:ccc). To this end, we write (25) as

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{[n]!}}(A+B)^n`$ $`=`$ $`{\displaystyle \underset{j,k,l,m,\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^{j+k+2l+3m+\mathrm{}}}{[j]_q![k]_q![l]_{q^2}![m]_{q^3}!}}A^jB^k𝒞_2^l𝒞_3^m\mathrm{}`$ (46)
and compare the coefficients of equal powers of $`x`$. This leads to coupled equations for $`𝒞_2`$, $`𝒞_3`$, $`𝒞_4`$, $`\mathrm{}`$, which can be solved in terms of $`A`$ and $`B`$ after straightforward algebra. Thus, solving the first few equations we confirm the forms of $`𝒞_2`$ and $`𝒞_3`$ already obtained above (compare (LABEL:eee) and (LABEL:ccc)) and find that

$`𝒞_4`$ $`=`$ $`{\displaystyle \frac{1}{[4]!}}(q^6[A,[A,[A,B]_{q^1}]_{q^2}]_{q^3}`$ (47)
$`+q^3[3][[A,[A,B]_{q^1}]_{q^2},B]_qq[3][[[A,B]_{q^1},B]_q,B]_{q^2}).`$
Now, using (25), (LABEL:eee), (LABEL:ccc) and (47), we get the $`q`$-Zassenhaus formula as

$`e_q^{(A+B)}=e_q^Ae_q^Be_{q^2}^{q[A,B]_{q^1}/[2]}e_{q^3}^{\{(q^3[A,[A,B]_{q^1}]_{q^2}/[3]!)(q[[A,B]_{q^1},B]_q/[3])\}}`$
$`\times e_{q^4}^{\frac{1}{[4]!}\left(q^6[A,[A,[A,B]_{q^1}]_{q^2}]_{q^3}+q^3[3][[A,[A,B]_{q^1}]_{q^2},B]_qq[3][[[A,B]_{q^1},B]_q,B]_{q^2}\right)}`$
(48)
up to the first five terms. In the limit $`q`$ $``$ $`1`$ it is found that $`𝒞_2`$, $`𝒞_3`$ and $`𝒞_4`$ become the expressions for the classical Zassenhaus formula (see, for example, (7)-(9) of ). When $`AB`$ $`=`$ $`q^1BA`$ it is found that $`𝒞_2`$ $`=`$ $`𝒞_3`$ $`=`$ $`𝒞_4`$ $`=`$ $`\mathrm{}`$ $`=`$ $`0`$ leading to the Schützenberger-Cigler result (12).

In the literature on $`q`$-series there are two other definitions of the $`q`$-exponential:

$`e_q(x)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{(1q)(1q^2)\mathrm{}(1q^n)}}`$ (49)
$`E_q(x)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{q^{n(n1)/2}x^n}{(1q)(1q^2)\mathrm{}(1q^n)}}.`$ (50)
such that

$`\underset{q1}{lim}e_q((1q)x)`$ $`=`$ $`e^x`$ (51)
$`\underset{q1}{lim}E_q((1q)x)`$ $`=`$ $`e^x`$ (52)
$`e_q(x)E_q(x)`$ $`=`$ $`1.`$ (53)
Let us note that the corresponding Zassenhaus formulae can be found in straightforward ways by using the relations:

$`e_q(x)`$ $`=`$ $`e_q^{x/(1q)}`$ (54)
$`E_q(x)`$ $`=`$ $`e_{q^1}^{x/(1q)}.`$ (55)
Thus, we find, using (48) and (54),

$`e_q(A+B)`$ $`=`$ $`e_q^{(A+B)/(1q)}`$ (56)
$`=`$ $`e_q(A)e_q(B)e_{q^2}(q[A,B]_{q^1}/(1q))\mathrm{}`$
and

$`e_q(A+B)`$ $`=`$ $`e_q(A)e_q(B)\text{if}AB=q^1BA.`$ (57)
Similarly,

$`E_q(A+B)`$ $`=`$ $`e_{q^1}^{(A+B)/(1q)}`$ (58)
$`=`$ $`E_q(A)E_q(B)E_{q^2}([A,B]_q/(1q))\mathrm{}`$
and

$`E_q(A+B)`$ $`=`$ $`E_q(A)E_q(B)\text{if}AB=qBA.`$ (59)
4. Conclusion

To summarize, it is found that the classical Zassenhaus formula

$`e^{A+B}`$ $`=`$ $`e^Ae^Be^{\frac{1}{2}[A,B]}e^{\frac{1}{6}[A,[A,B]]\frac{1}{3}[[A,B],B]}`$ (60)
$`\times e^{\frac{1}{4!}\left([A,[A,[A,B]]]+3[[A,[A,B]],B]3[[[A,B],B],B]\right)}\mathrm{}`$
has a $`q`$-analogue given by

$`e_q^{(A+B)}=e_q^Ae_q^Be_{q^2}^{q[A,B]_{q^1}/[2]}e_{q^3}^{\{(q^3[A,[A,B]_{q^1}]_{q^2}/[3]!)(q[[A,B]_{q^1},B]_q/[3])\}}`$
$`\times e_{q^4}^{\frac{1}{[4]!}\left(q^6[A,[A,[A,B]_{q^1}]_{q^2}]_{q^3}+q^3[3][[A,[A,B]_{q^1}]_{q^2},B]_qq[3][[[A,B]_{q^1},B]_q,B]_{q^2}\right)}`$
$`\times \mathrm{}`$ (61)
where the $`q`$-exponential is defined by

$`e_q^x`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{[n]!}}.`$ (62)
Thus, we have shown that while Katriel, Rasetti and Solomon  have proposed a $`q`$-analogue of the Zassenhaus formula in the form

$`e_q^{(A+B)}=e_q^Ae_q^Be_q^{c_2}e_q^{c_3}e_q^{c_4}e_q^{c_5}\mathrm{},`$ (63)
it is possible to have a $`q`$-Zassenhaus formula written also as

$`e_q^{(A+B)}=e_q^Ae_q^Be_{q^2}^{𝒞_2}e_{q^3}^{𝒞_3}e_{q^4}^{𝒞_4}e_{q^5}^{𝒞_5}\mathrm{}.`$ (64)
We have also explicitly found the first few terms of the disentanglement formula (64).

It should be noted that once a $`q`$-Zassenhaus formula is given for the $`q`$-exponential defined by (62) for the other two common definitions of the $`q`$-exponential (49, 50) found in the literature the corresponding $`q`$-Zassenhaus formulae can be found using the relations (54, 55).


text:
# Nuclear-polarization effect to the hyperfine structure in heavy multicharged ions

## Abstract

We have investigated the correction to the hyperfine structure of heavy multicharged ions, which is connected with the nuclear-polarization effect caused by the unpaired bound electron. Numerical calculations are performed for hydrogenlike ions taking into account the dominant collective nuclear excitations. The correction defines the ultimate limit of precision in accurate theoretical predictions of the hyperfine-structure splittings.

During the last decade a significant progress has been achieved in investigations of the hypefine structure (hfs) in heavy multicharged ions. To date, accurate measurements of hfs splittings were performed for a number of elements by using optical spectroscopy at the experimental storage ring at GSI in Darmstadt and at the Super-EBIT at the Lawrence Livermore National Laboratory . One of the major purposes of these investigations consists in testing nontrivial effects of bound-state QED in intense nuclear fields with an accuracy on the level of a few percent. Experiments on the hfs splitting allow one also to probe internal nuclear structure, in particular, the magnetic moment distribution within the nucleus (the so-called Bohr-Weisskopf effect). Though a rough estimate of the latter effect can be obtained in the framework of the single-nucleon model, its accurate calculation with the use of a microscopic nuclear theory still has to be performed.

Two attempts to calculate the extended magnetization distribution taking into account nuclear many-body corrections should be mentioned here. The first approach is based on the dynamic correlation model , which allows to include the excited-core configurations in the nuclear wave function. Numerical calculations performed in the framework of this model provide accurate values for the nuclear magnetic moments. However, no satisfactory agreement with experimental results for the hfs splittings has been achieved . The second approach, which also takes into account the nuclear core polarization by the unpaired nucleon, is based on the theory of finite Fermi systems . In this case, the correction due to the nuclear magnetization distribution calculated for hydrogenlike bismuth turned out to be too small compared with the corresponding value deduced from the experiment.

At present it appears to be unlikely to obtain an accurate theoretical value for the Bohr-Weisskopf correction. However, one can eliminate it to a large extend in combined measurements of the hfs splittings in H- and Li-like ions . For the ground state of 209Bi80+, two independent calculations have lead to similar predictions, $`797.1(2)`$ meV  and $`797.15(13)`$ meV , respectively. While both theoretical results agree with the experimental value of $`820(26)`$ meV , further efforts to measure the predicted splitting with higher precision have not been successfull.

In the present paper we evaluate a correction, which has not been previously considered in calculations of the hyperfine structure. The correction is connected with the nuclear-polarization effect caused by the unpaired bound electron. While the corresponding contribution is relatively small in the case of hydrogenlike ions, it turns out to be non-negligible for accurate theoretical predictions of the hfs splittings in Li-like heavy ions. Since the uncertainty of our calculation is comparable with the magnitude of the nuclear-polarization effect itself, the latter cannot be completely eliminated by extracting the Bohr-Weisskopf contribution in accordance with Shabaev’s idea . The effect under consideration sets a natural limit up to which one can test bound-state QED, even if a specific difference of the hfs splittings of H- and Li-like ions is introduced .

In the following, we shall consider a hydrogenlike ion with a nonzero-spin nucleus, so that the total angular momentum of an atom $`𝑭`$ is defined by the coupling of the nuclear spin $`𝑰`$ with the total angular momentum of an electron $`𝒋`$. The hfs energy shift due to the magnetic-dipole interaction is given by

$$\mathrm{\Delta }E_\mu (F)=\frac{1}{2}\left[F(F+1)I(I+1)j(j+1)\right]A_{n\kappa }.$$
(1)
Here $`A_{n\kappa }`$ is the hfs constant, which depends on the electron state characterized by the standard set of quantum numbers. For electron states with $`j=1/2`$, the hfs splitting between the levels with $`F=I+1/2`$ and $`F=I1/2`$ is just $`\mathrm{\Delta }E_{n\kappa }=(I+1/2)A_{n\kappa }`$. Since the nuclear size is rather small with respect to the radius of the electron orbit, the hyperfine structure can be fairly understood in the framework of the external-field approximation. This allows to treat the magnetic field of a nucleus as a perturbing potential in calculations of the hfs splitting.

Employing the Dirac equation for the electron in the external field of an infinitely heavy pointlike nucleus, one obtains ($`\mathrm{}=c=1`$)

$$A_{n\kappa }^\mathrm{D}=\alpha (\alpha Z)^3\frac{m^2}{m_p}\frac{g_I\kappa }{j(j+1)}\frac{[2\kappa (n_r+\gamma )N]}{N^4\gamma (4\gamma ^21)},$$
(2)
where $`\alpha =e^2`$ is the fine-structure constant $`(e>0)`$, $`g_I`$ is the nuclear $`g`$ factor, $`n_r=n|\kappa |`$ is the radial quantum number, $`n`$ is the principal quantum number, $`\kappa =(j+1/2)(1)^{j+l+1/2}`$, $`\gamma =\sqrt{\kappa ^2(\alpha Z)^2}`$, $`N=\sqrt{(n_r+\gamma )^2+(\alpha Z)^2}`$, and $`m`$ and $`m_p`$ are the electron and proton masses, respectively. Because of various nuclear and QED effects, the experimental value for the hfs constant deviates from its Dirac prediction (2). To describe the extended nucleus, one usually employs the following parametrization

$$A_{n\kappa }^{\mathrm{ext}}=A_{n\kappa }^\mathrm{D}(1\delta _{n\kappa })(1ϵ_{n\kappa }),$$
(3)
where the corrections $`\delta _{n\kappa }`$ and $`ϵ_{n\kappa }`$ account for the nuclear charge and magnetic moment distributions within the nucleus, respectively . In addition to the approximation (3), the radiative corrections should be taken into account. All leading QED effects for the hyperfine structure have been independently calculated by different groups and the numerical results are consistent . For Li-like ions, one can develop a perturbation theory with respect to the parameter $`1/Z`$, which accounts for corrections arising from the electron-electron interaction . The recoil correction, which is due to the finite nuclear mass, is negligibly small for heavy ions.

Here we consider a correction $`\mathrm{\Delta }A_{n\kappa }`$ to the hfs constant (3) due to the nuclear-polarization effect, which is caused by the bound electron. More precisely, a core-polarization part of the effect is considered only, which is due to collective nuclear excitations. The corresponding single-nucleon contributions should be generally considered beyond the external-field approximation. However, they are assumed to be completely negligible. To describe nuclear polarization, we adopt a relativistic field theoretical approach, which incorporates the many-body theory for virtual nuclear excitations within bound-state QED for atomic electrons . This approach has been successfully applied in calculations of nuclear-polarization effects to the Lamb shift  and to the bound-electron $`g`$ factor . To some extent, the present formulae are quite similar to those derived in Ref. . The correction under consideration may be represented by the sum of contributions, which are referred to as the irreducible, the reducible, and the vertex parts.

The irreducible contribution to the hfs constant can be written in terms of a multipole decomposition as follows:

$$\mathrm{\Delta }A_{n\kappa }^{\mathrm{irr}}=\frac{\alpha }{2\pi }\frac{eg_I\mu _N\kappa }{j(j+1)}\underset{L0}{}B(EL)\underset{n_1,\kappa _1}{}\left[C_{j\frac{1}{2}L0}^{j_1\frac{1}{2}}\right]^2\frac{n\kappa |F_L|n_1\kappa _1n_1\kappa _1|F_L|\overline{n\kappa }}{\epsilon _{n\kappa }\epsilon _{n_1\kappa _1}\mathrm{sgn}(\epsilon _{n_1\kappa _1})\omega _L},$$
(4)
where $`\mu _N=e/(2m_p)`$ is the nuclear magneton, $`\epsilon _{n\kappa }`$ is the one-electron energy, $`\omega _L=E_LE_0`$ are the nuclear excitation energies with respect to the ground-state energy $`E_0`$ of the nucleus and $`B(EL)=B(EL;0L)`$ are the corresponding reduced electric transition probabilities. The sum over $`n_1`$ runs over the entire Dirac spectrum, while the sum over $`\kappa _1`$ is restricted to those intermediate states, where $`l+l_1+L`$ is even. A two-component radial vector $`r|n\kappa `$ is determined by

$$r|n\kappa =\left(\genfrac{}{}{0pt}{}{P_{n\kappa }(r)}{Q_{n\kappa }(r)}\right),$$
(5)
where $`P_{n\kappa }(r)=rg_{n\kappa }(r)`$ and $`Q_{n\kappa }(r)=rf_{n\kappa }(r)`$, with $`g_{n\kappa }(r)`$ and $`f_{n\kappa }(r)`$ being the upper and lower components of the Dirac wave function, respectively. The radial shape parametrizing the nuclear transitions is carried by the functions

$$F_L(r)=\frac{4\pi }{(2L+1)R_0^L}\left[\frac{r^L}{R_0^{L+1}}\mathrm{\Theta }(R_0r)+\frac{R_0^L}{r^{L+1}}\mathrm{\Theta }(rR_0)\right]$$
(6)
in the case of multipole excitations with $`L1`$ and

$$F_0(r)=\frac{5\sqrt{\pi }}{2R_0^3}\left[1\left(\frac{r}{R_0}\right)^2\right]\mathrm{\Theta }(R_0r)$$
(7)
for monopole excitations, respectively. Here $`R_0`$ is an average radius assigned to the nucleus in its ground state. In Eq. (4), the matrix element is given by

$$a|F_L|b=\underset{0}{\overset{\mathrm{}}{}}𝑑rF_L(r)\left[P_a(r)P_b(r)+Q_a(r)Q_b(r)\right].$$
(8)
The perturbed vector $`r|\overline{n\kappa }`$, which follows as

$$r|\overline{n\kappa }=\underset{n^{}}{\overset{n^{}n}{}}\frac{n^{}\kappa |\sigma _xr^2|n\kappa }{\epsilon _{n\kappa }\epsilon _{n^{}\kappa }}r|n^{}\kappa ,$$
(9)
can be evaluated analytically by means of the generalized virial relations for the Dirac equation  (see also Refs. ). In Eq. (9), $`\sigma _x`$ is the Pauli matrix.

The reducible contribution reads

$$\mathrm{\Delta }A_{n\kappa }^{\mathrm{red}}=\frac{\alpha }{4\pi }A_{n\kappa }^{\mathrm{ext}}\underset{L0}{}B(EL)\underset{n_1,\kappa _1}{}\left[C_{j\frac{1}{2}L0}^{j_1\frac{1}{2}}\right]^2\frac{n\kappa |F_L|n_1\kappa _1^2}{[\epsilon _{n\kappa }\epsilon _{n_1\kappa _1}\mathrm{sgn}(\epsilon _{n_1\kappa _1})\omega _L]^2},$$
(10)
where $`A_{n\kappa }^{\mathrm{ext}}`$ is the hfs constant given by Eq. (3). The sum $`l+l_1+L`$ again should be even.

The nuclear core-polarization correction to the hfs constant due to the vertex part is conveniently represented as the sum of a pole term

$`\mathrm{\Delta }A_{n\kappa }^{\mathrm{pol}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{eg_I\mu _N\kappa }{\sqrt{j(j+1)(2j+1)}}}{\displaystyle \underset{L0}{}}B(EL){\displaystyle \underset{n_1,\kappa _1}{}}{\displaystyle \frac{(2j_1+1)^{3/2}}{\sqrt{j_1(j_1+1)}}}\left[C_{j\frac{1}{2}L0}^{j_1\frac{1}{2}}\right]^2\left\{{\displaystyle \genfrac{}{}{0pt}{}{j_1}{j}}{\displaystyle \genfrac{}{}{0pt}{}{j_1}{j}}{\displaystyle \genfrac{}{}{0pt}{}{1}{L}}\right\}`$ (11)
$`\times `$ $`{\displaystyle \frac{n_1\kappa _1|\sigma _xr^2|n_1\kappa _1n_1\kappa _1|F_L|n\kappa ^2}{[\epsilon _{n\kappa }\epsilon _{n_1\kappa _1}\mathrm{sgn}(\epsilon _{n_1\kappa _1})\omega _L]^2}}`$ (12)
and of a residual term

$`\mathrm{\Delta }A_{n\kappa }^{\mathrm{res}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{\pi }}{\displaystyle \frac{\sqrt{2}eg_I\mu _N\kappa }{\sqrt{j(j+1)(2j+1)}}}{\displaystyle \underset{L0}{}}B(EL){\displaystyle \stackrel{}{}_{n_1,n_2}}{\displaystyle \underset{\kappa _1,\kappa _2}{}}\sqrt{2j_2+1}C_{j\frac{1}{2}L0}^{j_1\frac{1}{2}}C_{j\frac{1}{2}L0}^{j_2\frac{1}{2}}C_{j_2\frac{1}{2}11}^{j_1\frac{1}{2}}`$ (13)
$`\times `$ $`\left\{{\displaystyle \genfrac{}{}{0pt}{}{j_1}{j}}{\displaystyle \genfrac{}{}{0pt}{}{j_2}{j}}{\displaystyle \genfrac{}{}{0pt}{}{1}{L}}\right\}{\displaystyle \frac{n\kappa |F_L|n_1\kappa _1n_2\kappa _2|F_L|n\kappa }{\epsilon _{n\kappa }\epsilon _{n_2\kappa _2}\mathrm{sgn}(\epsilon _{n_2\kappa _2})\omega _L}}{\displaystyle \frac{n_1\kappa _1|\sigma _xr^2|n_2\kappa _2}{\epsilon _{n_1\kappa _1}\epsilon _{n_2\kappa _2}}},`$ (14)
respectively. Here $`\mathrm{\Delta }A_{n\kappa }^{\mathrm{pol}}`$ accounts for the terms with $`n_1=n_2`$ and $`\kappa _1=\kappa _2`$ in the sums over intermediate states. The prime in the sum in Eq. (14) indicates that $`\epsilon _{n_1\kappa _1}\epsilon _{n_2\kappa _2}`$ when $`\kappa _1=\kappa _2`$, i.e., the pole contribution is supposed to be omitted. In Eqs. (12) and (14), the value $`l+l_1+L`$ has to be even. A second condition in Eq. (14) is that the sum $`l_1+l_2`$ should be even as well. The total nuclear core-polarization contribution to the hfs constant is determined by the sum of all contributions given by Eqs. (4), (10), (12), and (14).

In Table I, we present numerical results for some hydrogenlike ions, which are of particular experimental interest. The calculations were performed taking into account a finite set of dominant collective nuclear excitations. To estimate the nuclear parameters, $`\omega _L`$ and $`B(EL)`$, in the case of nearly spherical nucleus of $`{}_{83}{}^{}{}_{}{}^{209}`$Bi, we employed experimental data corresponding to the low-lying vibrational levels in neighbouring even-even isotope of $`{}_{82}{}^{}{}_{}{}^{208}`$Pb, which were deduced from nuclear Coulomb excitation. In the case of giant resonances, we utilized phenomenological energy-weighted sum rules , which are assumed to be concentrated in single resonant states. In the present calculations, contributions due to monopole, dipole, quadrupole, and octupole giant resonances have been taken into account. The infinite summations over the entire Dirac spectrum were performed by the finite basis set method. Basis functions are generated via B splines including nuclear-size effects .

Concluding, we have evaluated a correction to the hyperfine structure in heavy multicharged ions, which is connected with the nuclear core-polarization effect caused by the unpaired bound electron. The correction is exhausted over distances of the order of the nuclear size and it is enhanced due to a singular behavior of the hyperfine-interaction operator. The uncertainty of our calculation can be as large as the nuclear-polarization effect itself. It yet determines the natural limitation for testing higher-order QED corrections in future experiments aiming for accurate hfs measurements. In the case that the experimental value for the ground-state hfs splitting is used to eliminate the Bohr-Weisskopf effect , the latter cannot be separated from the effect we have considered. Accordingly, the utmost precision for theoretical hfs predictions in lithiumlike heavy ions is determined by the nuclear-polarization correction to the ground state in hydrogenlike ions. In particular, for the ground state in 209Bi82+, the nuclear-polarization effect contributes on the level of about $`0.05`$ meV. This implies that although numerical calculations of bound-state QED corrections provide sufficiently stable results, the conservative estimate of uncertainties for the hfs splitting in 209Bi80+ quoted in Ref.  appears to be more realistic rather then the one predicted in Ref. . It is also worth noting that because of a similar scaling dependence of nuclear and radiative corrections upon the principal quantum number, significant cancellations of almost all corrections, except for the screened QED contribution, occur in a specific difference of the hfs splittings in H- and Li-like heavy ions . In this case, the nuclear-polarization effect might be most relevant for determining the uncertainties of accurate theoretical predictions.

###### Acknowledgements.

A.N. is grateful for financial support from RFBR (Grant No. 01-02-17246) and from the Alexander von Humboldt Foundation. G.S. and G.P acknowledge financial support from BMBF, DFG, and GSI.

|  | $`I^\pi `$ | $`\mu `$ (nm) | | $`\mathrm{\Delta }A_{1s}`$ (meV) | | | | | | $`\mathrm{\Delta }E_{1s}`$ (meV) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 159Tb64+ | $`3/2^+`$ | 2 | $`.014`$ | 0 | $`.25(2)^\mathrm{a}`$ | 0 | $`.26(2)^\mathrm{b}`$ | 0 | $`.51(2)^\mathrm{c}`$ | $`0.10(1)`$ |
| 165Ho66+ | $`7/2^{}`$ | 4 | $`.132`$ | 0 | $`.38(3)`$ | 0 | $`.27(2)`$ | 0 | $`.31(2)`$ | $`0.12(1)`$ |
| 175Lu70+ | $`7/2^+`$ | 2 | $`.2327`$ | 0 | $`.13(2)`$ | 0 | $`.22(2)`$ | 0 | $`.35(2)`$ | $`0.14(1)`$ |
| 187Re74+ | $`5/2^+`$ | 3 | $`.2197`$ | 0 | $`.25(2)`$ | 0 | $`.66(2)`$ | 0 | $`.91(2)`$ | $`0.27(1)`$ |
| 203Tl80+ | $`1/2^+`$ | 1 | $`.62226`$ | 0 | $`.16(2)`$ | 0 | $`.29(1)`$ | 0 | $`.31(1)`$ | $`0.31(1)`$ |
| 209Bi82+ | $`9/2^{}`$ | 4 | $`.1106`$ | 0 | $`.81(3)`$ | 0 | $`.10(1)`$ | 0 | $`.11(1)`$ | $`0.55(1)`$ |
| 235U91+ | $`7/2^{}`$ | $`0`$ | $`.38`$ | $`0`$ | $`.36(2)`$ | $`0`$ | $`.29(2)`$ | $`0`$ | $`.65(2)`$ | $`0.26(1)`$ |