text:
# Fracture and second-order phase transitions.

## Abstract

Using the global fiber bundle model as a tractable scheme of progressive fracture in heterogeneous materials, we define the branching ratio in avalanches as a suitable order parameter to clarify the order of the phase transition occurring at the collapse of the system. The model is analyzed using a probabilistic approach suited to smooth fluctuations. The branching ratio shows a behavior analogous to the magnetization in known magnetic systems with 2nd-order phase transitions. We obtain a universal critical exponent $`\beta 0.5`$ independent of the probability distribution used to assign the strengths of individual fibers.

preprint: HEP/123-qed

The interest in the fracture processes of heterogeneous media has increased in the last years . In the lab, a disordered material subjected to an increasing external load can be studied by measuring the acoustic emissions before the global rupture. It has been shown  that this intense precursory activity in the form of bursts of different microscopic sizes follows a well-defined power law. In despite of the many efforts and successes that have been recently achieved, the question of whether rupture exhibits the properties of a first-order or a second-order phase transition remains under discussion as well as what is the order parameter that indicates the type of transition.

From the theory side, the understanding of fracture in heterogeneous materials has progressed due to the use of lattice models and large scale simulations . In this field, it is important to use models able to describe the complexity of the rupture process, nevertheless, they should be simple enough as to permit analytical insights. To this class of models belong the well-known Fiber Bundle Models (FBM) widely used since their introduction more than forty years ago . In static FBM, a set of fibers (elements) is located on a supporting lattice and one assigns to its elements a random strength threshold sampled from a probability distribution. The set is loaded and fibers break when their loads exceed their threshold values. In the Equal Load Sharing (or global) FBM, which is the simplest scheme one can adopt to make the problem analytically tractable, one assumes that the load carried by failed elements is equally distributed among the surviving elements of the system.

In the present Letter, we explore the criticality of fracture in the global FBM using a novel probabilistic approach devised to smooth fluctuations. The scaling relations obtained and the behavior of the order parameter point out that the fracture of a fiber set with long range interaction undergoes a continuous phase transition.

Let us first recall the basic ingredients of the global FBM and how one proceeds in numerical simulations. The system under consideration is a set of $`N_0`$ elements located in a supporting lattice each one having at the initial state a zero load and a fixed strength threshold value sampled randomly from a probability distribution $`P(\sigma )`$. The system is then subjected to an external force $`F`$ such that each element increases its load in the same amount, $`\sigma `$. This individual stress (or load), $`\sigma `$, will act as the control parameter. The process of driving is done quasi-statically; i.e., the external force is increased at a sufficiently slow rate as to produce a single breaking event when the stress on the weakest element equals its threshold value. Then, the increase of $`F`$ stops and the load of the broken element is equally transferred. This implies that the load on any element is given by $`\sigma =F/n_s(F)`$, where $`n_s`$ is the number of surviving elements for a given $`F`$. The rupture of an element may induce secondary failures which in turn may trigger more failures and so on. This process of induced failure at constant external load, termed an avalanche, stops when all surviving elements carry a load lower than their thresholds. The system is then loaded again and the process is repeated until the final catastrophic avalanche provokes the total rupture of the material, which occurs at a critical load $`\sigma _c`$ that depends on the probability distribution from where the individual strengths were drawn, as well as on the system size. The FBM have been recently used in self-organized criticality (SOC), a theoretical framework widely used for the study of avalanche phenomena in disordered systems. It has been shown using these models that systems with plastic behavior can reach a SOC state just before the global rupture . A second case of self-organization with power law distributions in several quantities corresponds to the situation in which the fracture process coexists with a healing process .

In numerical simulations, the cycle of complete breakdown of the model is performed many times in order to average out the effect of fluctuations and obtain mean values. As we are interested in studying the behavior of the system as the critical point, or point of final collapse, is approached, it is of most importance to find a simple method able to capture the evolution of the system avoiding as much as possible the fluctuations appearing in numerical simulations. To introduce our probabilistic strategy, let there be a large set of $`N_0`$ elements. Suppose that each element carries a given load $`\sigma `$, which is zero at the initial state. The strength of each element is drawn from a probability distribution $`P(\sigma )`$. Different probability distributions can be considered. In materials science the Weibull distribution is widely used,

$$P(\sigma )=1e^{(\frac{\sigma }{\sigma _0})^\rho },$$
(1)
$`\rho `$ being the so-called Weibull index, which controls the degree of disorder in the system (the bigger the Weibull index, the smaller the disorder), and $`\sigma _0`$ is a load of reference. In the following we will assume $`\sigma _0=1`$, and therefore the loads are dimensionless. At this point, it is worth noting that the results and the formulas derived in the following hold for a wide class of probability distributions. We use here the Weibull distribution for definiteness but results have also been obtained for other distributions. Equation (1) represents the probability that an element fails under the individual load $`\sigma `$. Now, consider the case in which an element drawn from Eq. (1) supports a load $`\sigma _1`$ but breaks under a new load $`\sigma _2`$. The probability that this happens is given by

$$p(\sigma _1,\sigma _2)=\frac{P(\sigma _2)P(\sigma _1)}{1P(\sigma _1)}=1e^{(\sigma _2^\rho \sigma _1^\rho )}.$$
(2)
So, the probability $`q(\sigma _1,\sigma _2)`$ that an element that has survived to the load $`\sigma _1`$ also survives to the load $`\sigma _2`$ will be given by $`q(\sigma _1,\sigma _2)=1p(\sigma _1,\sigma _2)=e^{(\sigma _2^\rho \sigma _1^\rho )}`$.

To mimic the quasi-static increase in load on the system we impose the condition that under an external force $`F`$, the next breaking event consists of one single failure. Let suppose that after the latest avalanche, there are $`N_k`$ surviving elements each one bearing a load $`\sigma _k`$. The new individual load $`\sigma _l`$ needed to provoke the failure of just one more element is given by the solution of $`N_k1=N_kq(\sigma _k,\sigma _l)`$. Thus,

$$\sigma _l=\left[\sigma _k^\rho \mathrm{ln}\left(1\frac{1}{N_k}\right)\right]^{\frac{1}{\rho }},$$
(3)
where in Eq. (3) $`N_k=N_0`$ and $`\sigma _k=0`$ at the initial state. Elevating the external force up to the $`N_k\sigma _l`$ level, statistically speaking, one element breaks. As we are dealing with an equal load sharing set, the choice of the broken element is irrelevant because all the elements are equivalent. Once the first element fails, the redistribution of its stress takes place. This may induce other failures until the end of the avalanche.

Now, how many elements will survive to the situation in which $`n_1`$ elements with load $`\sigma _1`$ fail in an avalanche step? The new load on the intact $`N_1n_1=N_2`$ elements is $`\sigma _2=\frac{N_1\sigma _1}{N_2}`$. So, the number $`N_3`$ of elements that survive to the new load can be expressed as

$$N_3=N_2q(\sigma _1,\frac{N_1}{N_2}\sigma _1)=N_2q(\sigma _1,\sigma _2).$$
(4)
As a consequence of applying Eq. (4), $`N_2N_3`$ elements break and the new total number of intact fibers will support a bigger load $`\sigma _3`$. The avalanche may continue and Eq. (4) is applied again for the set of $`N_3`$ surviving elements. The iterative process will stop when no new element fails, which occurs when the right-hand side is equal to the left-hand side in Eq. (4). The general form of Eq. (4) is

$$N_{j+1}=N_jq(\sigma _{j1},\sigma _j),$$
(5)
with the conservation condition for the total load in the system during an avalanche

$$N_j\sigma _j=N_{j1}\sigma _{j1}$$
(6)
and the condition

$$N_j=N_{j+1}$$
(7)
which determines the end of the avalanche. The dynamics of the system is determined by Eq. (3), (5), (6). The size of an avalanche is given by the number of elements that break between two successive steps of external loading. The critical load, defined as the load needed to provoke the total collapse of the system, is equal to the load on the intact elements just before the final catastrophic avalanche. Note that in this probabilistic approach, in contrast to Monte Carlo simulations, we only need to store the information concerning the loads on the intact elements, that is, the details of the threshold list is omitted.

In the probabilistic strategy, we can proceed in two different ways in order to determine when an avalanche ends, to which we will refer as the continuous and the discrete cases. For the continuous case, the number $`N_{j+1}`$ of surviving elements is considered a real number. Strictly speaking, this means that the condition (7) is never fulfilled. So, the condition (7) is replaced in numerical calculations by using a factor $`\nu 1`$ that determines the end of an avalanche, i.e., if $`N_jN_{j+1}\nu `$ the avalanche stops; otherwise it continues. In the discrete case, $`N_{j+1}`$ is considered to be a whole number, so that after each iteration of Eq. (5), $`N_{j+1}`$ has to be rounded up. This is done comparing the remainder of $`N_{j+1}`$, $`\lambda `$, with a random number $`\alpha `$ uniformly distributed in the interval $`[0,1[`$. Thus, if $`\alpha \lambda `$, $`N_{j+1}`$ is equal to its whole part, otherwise, $`N_{j+1}`$ is equal to its whole part plus one. Next, we check whether the condition (7) is satisfied for the rounded value of $`N_{j+1}`$ or if a new iteration of Eq. (5) has to be performed. The continuous approach has the advantage that the fluctuations are ruled out, whereas for the discrete case the results are similar to those obtained by Monte Carlo simulations where it is necessary to average over many realizations in order to get accurate mean values. Remember that in this model the central limit theorem applies .

In Fig. 1 we have depicted the fraction of broken elements versus $`\sigma `$, for the continuous case and for four individual Monte Carlo simulations with a Weibull index $`\rho =2`$ and $`N_0=5000`$. No averaging has been done because our aim is just to illustrate the scatter of the results. As can be seen, the continuous probabilistic model gives a smooth curve, and provides an accurate value for the critical load $`\sigma _c`$, which analytically is given by $`\sigma _c=(\rho e)^{1/\rho }`$  in the limit of infinite $`N_0`$.

Now, we proceed to explore the behavior of some quantities as the critical point is reached. The results shown below have been obtained using the continuous approach ($`\rho =2`$). In Fig. 2 we show the interesting scaling relation for the average avalanche size. It turns out that the avalanche size near to the critical point diverges with an exponent $`\gamma =\frac{1}{2}`$ as $`s(\sigma _c\sigma )^\gamma `$. A similar behavior, through a mapping of a fuse network model to the global fiber bundle model used here, has been recently reported . We have also obtained the same scaling function for the derivative of the number of broken fibers with respect to the load on the system. The rate $`dN/d\sigma `$ diverges as $`(\sigma _c\sigma )^{\frac{1}{2}}`$, thus qualifying a critical mean field behavior as was already shown in Ref. . In Ref. , a similar scaling behavior is addressed for the derivative of the strain carried by the fibers with respect to the driving force.

Another way to shed light on the critical behavior of this type of system is to define a branching ratio $`\zeta `$ for each avalanche. This magnitude represents the probability to trigger future breaking events given an initial individual failure , and is related to the number of broken fibers by

$$\zeta =\frac{<z>1}{<z>}.$$
(8)
The above relation can be obtained by thinking of the evolution of fracture as a kind of branching process . In this process, each node gives rise to a number $`n`$ of new branches in the next time step. The average number $`<n>`$ of new branches is called the branching ratio. Let us denote by $`n_t`$ the number of branches at a given step $`t`$ of the branching process, and by $`t_{max}`$ the total number of time steps before it dies. Then, $`\zeta =1\frac{n_0}{_{t=0}^{t_{max}}n_t}`$. As, $`n_0=1`$, $`\zeta =1\frac{1}{n_{tot}}`$ where $`n_{tot}`$ is the total number of nodes developed in the branching process. For a fracture process, $`n_{tot}`$ is equal to the average number of failure events. So, Eq. (8) defines the branching ratio. We represent by $`<z>`$ the average number of elements that fail in one avalanche, which is a function of $`\sigma `$ and coincides with $`s`$. This analogy between fracture and branching processes has been previously used to study the criticality in the process of fragmentation of Hg drops . The branching ratio will then act as the order parameter. It takes the value 1 when the system is critical thereby representing a measure of the distance of the system from the critical state . We have numerically computed $`\zeta `$ by means of the continuous method. The results obtained for a system of $`N_0=50000`$ elements and several values of $`\rho `$ have been plotted in Fig. 3. It can be seen in this figure that $`1\zeta `$ approaches zero, in all cases, in a continuous fashion as the critical load is reached. Near the critical point, the relation $`1\zeta (\sigma _c\sigma )^\beta `$, with $`\beta =\frac{1}{2}`$ applies. This is the exponent of the order parameter that we should expect from a mean-field approach. Note the similarity of the figure with those obtained for the magnetization in known magnetic systems with second-order phase transitions. In the figure, the values of $`\zeta `$ are collected for all the avalanches except for the last, that which provokes the collapse of the system. The result that at $`\sigma _c`$ $`\zeta 1`$, is consistent with the previous result that the avalanche size diverges at the critical point. On the other hand, the branching ratio does not depend on the size of the system for large systems, in contrast with previous results in other fracturing systems .

The suggestion of Ref. is that fracture can be seen as a first-order phase transition close to a spinodal-like instability . There, by simulating models of electric breakdown and fracture, the authors present numerical and theoretical evidence of several scaling relations and of a discrete jump in some macroscopic properties. Here, we have obtained the same scaling relation for the rate of fiber failures as the critical point is reached, and for the avalanche sizes, which also diverge at that point. Our numerical results also fit the mean field result $`\gamma =\frac{1}{2}`$. It is true that the fraction of unbroken fibers just before the global rupture has a discontinuity; but from our point of view that is not enough as to set the conclusion that fracture can be describe as a first-order phase transition, since the concepts related to spinodal nucleation are not sufficiently well established in driven disordered systems.

Our alternative point of view has been to consider the branching ratio as an appropriate order parameter. According to the results obtained, the branching ratio goes continuously from zero to one. Note, additionally, that what changes discontinuously at $`\sigma _c`$ is the rate of change of $`\zeta `$ rather than $`\zeta `$ itself, which is in the essence of a continuous phase change. Therefore, the behavior of the branching ratio implies that the system undergoes a second-order phase transition as claimed in other analysis of fracture models . Our results suggest that fracture in heterogeneous systems with long range interaction can be described as a phase transition of the second-order type. Besides, it is important to bear in mind that in fiber models with local interactions, the order parameter $`\zeta `$, has a discontinuous jump typical of first order phase transitions. Finally, it is worth recalling that fracture of real materials is a process based on elasticity, and elasticity is a long-range phenomenon. In this respect, the global load-sharing model we have explored here could be a better analogy to real fracture than the local one.

Y.M thanks A. Vespignani and H. J. Jensen for useful and stimulating discussions held at the Abdus Salam International Centre for Theoretical Physics where this work began. Y. M. is supported by the AECI. This work was supported in part by the Spanish DGICYT under Project PB98-1594.


format: text

text:
# The Effects of Rotation on the Evolution of Rising Omega-loops in a Stratified Model Convection Zone

## 1 Introduction

Active regions represent the largest observable concentrations of magnetic flux on the Sun. The magnetic field within active regions exhibits a bipolar structure, which suggests that they are the tops of large $`\mathrm{\Omega }`$-shaped loops which have risen through the convection zone and have emerged through the photosphere. On average, active region bipoles are oriented nearly parallel to the East-West direction (Hale’s Polarity Law). This suggests that the subsurface field geometry is approximately toroidal. Since active regions appear to be the tops of large, rising $`\mathrm{\Omega }`$-loops, and since Hale’s law persists for years in a given solar cycle, we can infer that the toroidal flux must exist well below the solar surface, and that the magnetic flux must reside in a region where it is relatively free from disruption by buoyant instabilities or convective motions (Spruit & van Ballegooijen, 1982; van Ballegooijen, 1982; Ferriz-Mas & Schüssler, 1993, 1995). However, if the magnetic flux is embedded in the radiative zone (well below the convective zone) it would be so stable that it could not emerge through the surface on the time scale of a solar cycle (Parker, 1979; van Ballegooijen, 1982). This implies that the toroidal field likely resides in the “convective overshoot region” — a thin, slightly convectively stable layer located at the interface of the radiative and convective zones. This layer also appears to coincide with the “tachocline” (Kosovichev, 1996; Corbard et al., 1999), where the solid body rotational behavior of the radiative interior transitions into the differential rotation behavior observed in the convection zone. Thus, the toroidal field seems to be generated and stored in the overshoot layer. (Gilman, Morrow, & DeLuca, 1989; DeLuca & Gilman, 1991; Parker, 1993; MacGregor & Charbonneau, 1997; Durney, 1997; Dikpati & Charbonneau, 1999).

We now have a basic picture of how active regions might form: A portion of the toroidal flux layer eventually succumbs to one of several possible instabilities (eg. Caligari et al. 1995; Fan & Fisher 1996; Fan 2000; Wissink et al. 2000b), and forms a magnetically buoyant tube that begins to rise through the convection zone toward the photosphere. The top of the emerging flux loop is then observed at the photosphere as a bipolar magnetic region. The simplest way to model this picture is with the “thin flux tube” approximation (see Defouw 1976; Roberts & Webb 1978; Spruit 1981). In this simple model, it is assumed that magnetic flux tubes move through a field-free plasma, that the cross-section of the tube is much smaller than all other relevant length scales of the problem, and that pressure balance is always maintained across the tube. Given these assumptions, one can derive an equation of motion for a one-dimensional tube moving through a three dimensional medium. This model has had great success in explaining many of the observed properties of active regions, including (but not limited to) the variation and scatter of active region tilt as a function of solar latitude (D’Silva & Howard, 1993; D’Silva & Choudhuri, 1993; Fan et al., 1993; Caligari et al., 1995; Fisher et al., 1995; Longcope & Fisher, 1996; Fan & Fisher, 1996), the asymmetries in the orientation of active regions just after emergence (van Driel-Gesztelyi & Petrovay, 1990; Moreno-Insertis et al., 1994), and the latitudinal variation and fluctuations of field line twist from active region to active region (Longcope et al., 1998).

However, from the start, the thin flux tube approximation assumes that the magnetic field will retain its “tube-like” identity; that is, it will not fragment or disperse as it evolves. This assumption is called into question by recent two-dimensional MHD simulations which show that unless flux tubes possess a critical amount of initial field line twist along the tube — and observations of emerging active regions suggest that many do not (Pevtsov et al., 1995; Longcope et al., 1998, 1999) — the tubes will fragment (break apart) before they are able to emerge through the surface (Linton et al., 1996; Moreno-Insertis & Emonet, 1996; Emonet & Moreno-Insertis, 1998; Fan et al., 1998). Recently, Abbett et al. (2000) used three-dimensional MHD simulations to demonstrate that the degree of cohesion at the apex of an emerging $`\mathrm{\Omega }`$-loop depends strongly on the three-dimensional geometry of the loop. If the curvature of the $`\mathrm{\Omega }`$-loop at its apex is relatively large (ie. if the loop is tall and narrow), then the apex will emerge through the photosphere more cohesively than if the curvature is small (ie. if the $`\mathrm{\Omega }`$-loop is short and wide). Dorch & Nordlund (1998) and Dorch et al. (1999) also performed detailed three-dimensional MHD simulations of buoyant magnetic flux tubes. All of these 3-D studies still require a significant amount of initial field line twist in order to suppress the eventual fragmentation of the tube near the loop apex. Is there a way out this quandary which Abbett et al. (2000) refer to as “Longcope’s paradox”?

In this paper, we investigate how solar rotation affects the dynamic evolution of buoyant magnetic flux tubes. We find that the Coriolis force plays a dynamically important role during the tube’s ascent, and acts to maintain the tube’s cohesion even in the absence of initial field line twist, possibly resolving Longcope’s paradox. In Section 2 we briefly describe our computational methodology, discuss how we include the Coriolis force in our computational model, and define the relevant physical parameters of our problem. In Section 3 we analyze the results of our simulations, and develop a very simple model that effectively characterizes the simulation results pertaining to flux tube cohesion. In this section we also comment on possible observational consequences of the simulations. Finally, we present our summary and conclusions in Section 4.

## 2 Method

We use the methods described in Abbett et al. (2000) and Fan et al. (1999) to solve the three-dimensional equations of MHD in the anelastic approximation (Ogura & Phillips, 1962; Gough, 1969; Lantz & Fan, 1999). The anelastic equations result from an expansion of the equations of compressible MHD about a zeroth order reference state (which we take to be a field-free, adiabatically stratified polytrope of index $`m=1.5`$). In this approximation, the higher order time derivative of the perturbed density is neglected, effectively filtering out fast-moving acoustic waves. This is computationally advantageous, as it allows for much larger timesteps than would be possible in a fully compressible calculation. This approximation is valid in sub-surface regions where the acoustic Mach number is small ($`Mv/c_\mathrm{s}1`$), and where the local sound speed ($`c_\mathrm{s}`$) greatly exceeds the Alfvén speed ($`v_A`$) of the plasma (Glatzmeier, 1984). A detailed theoretical discussion of the anelastic formulation can be found in Lantz & Fan 1999 (and references therein).

The anelastic equations in a non-rotating reference frame are given by equations - in Abbett et al. (2000) (hereafter AFF). To consider the effects of solar rotation, we adjust the anelastic momentum equation to include non-inertial terms (we note, here, that due to a typesetting error in the final version of AFF, a cross-product symbol was omitted in the Lorentz force term of equation ). To consider the effects of solar rotation, we adjust the anelastic momentum equation to include non-inertial terms. In this study, we neglect centrifugal effects and concentrate only on the effects of the Coriolis force due to solid-body rotation. In this approximation, the anelastic momentum equation (in the rotating reference frame) becomes:

$$\rho _0\left(\frac{𝐯}{t}+𝐯\mathbf{}𝐯\right)=\mathbf{}p_1+\rho _1𝐠2\rho _0\left(𝛀\times 𝐯\right)+\frac{1}{4\pi }\left(\mathbf{}\times 𝐁\right)\times 𝐁+\mathbf{}𝚷.$$
(1)
Here, $`\rho _0`$ refers to the density stratification of the zeroth order reference state, while $`\rho _1`$, $`p_1`$, $`𝐯`$, and $`𝐁`$ represent the density, pressure, velocity and magnetic field perturbations. The viscous stress tensor is given by $`\mathrm{\Pi }_{ij}\mu (v_i/x_j+v_j/x_i2/3(\mathbf{}𝐯)\delta _{ij})`$, where $`\mu `$ is the coefficient of viscosity and $`\delta _{ij}`$ is the Kronecker delta function. $`𝛀`$ denotes the angular velocity about the Sun’s axis of rotation. Both the acceleration due to gravity, $`𝐠\mathbf{=}\mathbf{}\text{g}\widehat{𝒛}`$, and the coefficient of viscosity are assumed constant in our calculations.

We solve the non-dimensional form of the anelastic equations in the modified, local *f*-plane approximation of Brummell et al. (1996). This approach is similar to the standard *f*-plane approximation of geophysical fluid dynamics (Pedlosky, 1987), except that the transverse components of the Coriolis force within the rectangular domain are retained. Since $`𝛀`$ is assumed constant throughout the Cartesian box (positioned at a given latitude), we can apply simple, periodic boundary conditions in the horizontal directions, and stress-free, non-penetrating conditions at the upper and lower boundaries. Further simplification is possible by noting that the momentum density (in the anelastic approximation) and magnetic field vectors are both divergence-free. This allows us to define scalar potentials such that $`𝐁\mathbf{}\mathbf{}\mathbf{\times }\mathbf{}\mathbf{\times }𝓑`$ $`\widehat{𝒛}`$ $`\mathbf{+}\mathbf{}\mathbf{\times }𝓙`$ $`\widehat{𝒛}`$  and
$`𝝆_\mathrm{𝟎}𝐯\mathbf{}\mathbf{}\mathbf{\times }\mathbf{}\mathbf{\times }𝓦`$ $`\widehat{𝒛}`$ $`\mathbf{+}\mathbf{}\mathbf{\times }𝓩`$ $`\widehat{𝒛}`$ . Along with the other dependent variables of the system, these potentials are decomposed spectrally in the horizontal directions:

$$𝓠\mathbf{=}\underset{\mathrm{𝐦𝐧}}{\mathbf{}}\stackrel{\mathbf{~}}{𝓠}_{\mathrm{𝐦𝐧}}\mathbf{(}𝒛\mathbf{)}𝒆^{\mathrm{𝟐}𝝅𝒊\mathbf{\left(}𝒇_𝐦𝒙\mathbf{+}𝒈_𝐧𝒚\mathbf{\right)}}\mathbf{.}$$
(2)
Here, $`𝓠`$ and $`\stackrel{\mathbf{~}}{𝓠}_{\mathrm{𝐦𝐧}}`$ refer to a given dependent variable or scalar potential and its associated Fourier coefficient. $`𝒇_𝐦`$ and $`𝒈_𝐧`$ are defined as $`𝐦\mathbf{/}𝓛_𝒙`$ and $`𝐧\mathbf{/}𝓛_𝒚`$ respectively, where $`𝓛_𝒙`$ and $`𝓛_𝒚`$ denote the size of the domain in the $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ directions. The integers $`𝐦`$ and $`𝐧`$ range from $`\mathbf{}𝐍_𝒙\mathbf{/}\mathrm{𝟐}\mathbf{+}\mathrm{𝟏}`$ to $`𝐍_𝒙\mathbf{/}\mathrm{𝟐}`$ and $`\mathbf{}𝐍_𝒚\mathbf{/}\mathrm{𝟐}\mathbf{+}\mathrm{𝟏}`$ to $`𝐍_𝒚\mathbf{/}\mathrm{𝟐}`$. $`𝐍_𝒙`$ and $`𝐍_𝒚`$ refer to the number of zones in the $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ directions. The Fourier coefficients are then discretized with respect to the vertical direction (defined to be parallel to $`𝐠`$), and vertical derivatives are approximated by fourth-order, centered finite differences. Using the semi-implicit method of operator splitting, we apply the second-order Adams-Bashforth scheme to the advection terms, and the second-order Crank-Nicholson scheme to the diffusion terms (see Press et al. 1986 for a description of these methods).

A detailed, in-depth description of the techniques we employ to obtain a solution to the anelastic equations in a non-rotating reference frame can be found in Appendix A of Fan et al. (1999) (hereafter FZLF). Only minor modifications are necessary to extend that formalism to apply in the case of a rotating reference frame in the modified *f*-plane approximation. The equation for the time evolution of $`\stackrel{\mathbf{~}}{𝓩}_{\mathrm{𝐦𝐧}}`$, the Fourier coefficients of the spectral decomposition of $`𝓩`$, is obtained by evaluating the $`\widehat{𝒛}`$-component of the curl of the momentum equation. Thus, the addition of the Coriolis term to equation (1) results in an additional non-inertial term to equation (A16) of FZLF:

$$\mathbf{\left(}\frac{\mathbf{}\stackrel{\mathbf{~}}{𝓩}_{\mathrm{𝐦𝐧}}}{\mathbf{}𝒕}\mathbf{\right)}_{\mathrm{𝐜𝐨𝐫}}\mathbf{=}\mathbf{\hspace{0.33em}2}𝛀_𝒛\frac{\mathbf{}\stackrel{\mathbf{~}}{𝓦}_{\mathrm{𝐦𝐧}}}{\mathbf{}𝒛}\mathbf{+}\mathrm{𝟒}𝝅𝒊\stackrel{\mathbf{~}}{𝓦}_{\mathrm{𝐦𝐧}}\mathbf{\left(}𝒇_𝐦𝛀_𝒙\mathbf{+}𝒈_𝐧𝛀_𝒚\mathbf{\right)}\mathbf{.}$$
(3)
Here, the $`𝛀_𝒊`$’s represent the Cartesian components of the angular velocity vector, and $`\stackrel{\mathbf{~}}{𝓦}_{\mathrm{𝐦𝐧}}`$ denotes the Fourier coefficients of the scalar potential $`𝓦`$. Equation (A18) of FZLF results from taking the $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ components of the curl of the momentum equation as it appears in equation (A4). Thus, we must add another term to equation (A18) to reflect the addition of the Coriolis force to equation (1):

$$\mathbf{\left(}\frac{\mathbf{}\stackrel{\mathbf{~}}{𝝍}_{\mathrm{𝐦𝐧}}}{\mathbf{}𝒕}\mathbf{\right)}_{\mathrm{𝐜𝐨𝐫}}\mathbf{=}\mathbf{}\mathrm{𝟐}𝛀_𝒛\frac{\mathbf{}}{\mathbf{}𝒛}\mathbf{\left(}\frac{\stackrel{\mathbf{~}}{𝓩}_{\mathrm{𝐦𝐧}}}{𝝆_\mathrm{𝟎}}\mathbf{\right)}\mathbf{+}\frac{\mathrm{𝟒}𝝅𝒊}{𝝆_\mathrm{𝟎}𝒉_\mathrm{𝟎}}\stackrel{\mathbf{~}}{𝓦}_{\mathrm{𝐦𝐧}}\mathbf{\left(}𝒇_𝐦𝛀_𝒚\mathbf{}𝒈_𝐧𝛀_𝒙\mathbf{\right)}\mathbf{}\frac{\mathrm{𝟒}𝝅𝒊}{𝝆_\mathrm{𝟎}}\stackrel{\mathbf{~}}{𝓩}_{\mathrm{𝐦𝐧}}\mathbf{\left(}𝒇_𝐦𝛀_𝒙\mathbf{+}𝒈_𝐧𝛀_𝒚\mathbf{\right)}\mathbf{.}$$
(4)
In the above equation, $`\stackrel{\mathbf{~}}{𝝍}_{\mathrm{𝐦𝐧}}\mathbf{}\mathbf{}𝒊\mathbf{(}\mathrm{𝟐}𝝅\mathbf{)}^\mathbf{}\mathrm{𝟏}\mathbf{(}𝒇_𝐦^\mathrm{𝟐}\mathbf{+}𝒈_𝐧^\mathrm{𝟐}\mathbf{)}^\mathbf{}\mathrm{𝟏}\mathbf{(}𝒇_𝐦\stackrel{\mathbf{~}}{𝝎}_{𝒚\mathrm{𝐦𝐧}}\mathbf{}𝒈_𝐧\stackrel{\mathbf{~}}{𝝎}_{𝒙\mathrm{𝐦𝐧}}\mathbf{)}`$ is the scalar variable that corresponds to $`\stackrel{\mathbf{~}}{𝛀}_{\mathrm{𝐦𝐧}}`$ in equation (A18) of FZLF (we change notation so as to avoid confusion with the components of the angular velocity vector); $`\stackrel{\mathbf{~}}{𝝎}_{𝒙\mathrm{𝐦𝐧}}`$ and $`\stackrel{\mathbf{~}}{𝝎}_{𝒚\mathrm{𝐦𝐧}}`$ are the Fourier coefficients of the $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ components of the vorticity; and $`𝒉_\mathrm{𝟎}\mathbf{}\mathbf{}𝝆_\mathrm{𝟎}\mathbf{(}\mathbf{}𝝆_\mathrm{𝟎}\mathbf{/}\mathbf{}𝒛\mathbf{)}^\mathbf{}\mathrm{𝟏}`$ is the density scale height. The last necessary modification is to equations (A23) and (A24) of FZLF where we add the contribution from the horizontal average of the Coriolis term in the momentum equation: $`\mathbf{(}\mathbf{}\mathbf{(}𝝆_\mathrm{𝟎}\overline{𝒗}_𝒙\mathbf{)}\mathbf{/}\mathbf{}𝒕\mathbf{)}_{\mathrm{𝐜𝐨𝐫}}\mathbf{=}\mathrm{𝟐}𝝆_\mathrm{𝟎}𝛀_𝒛\overline{𝒗}_𝒚`$ and $`\mathbf{(}\mathbf{}\mathbf{(}𝝆_\mathrm{𝟎}\overline{𝒗}_𝒚\mathbf{)}\mathbf{/}\mathbf{}𝒕\mathbf{)}_{\mathrm{𝐜𝐨𝐫}}\mathbf{=}\mathbf{}\mathrm{𝟐}𝝆_\mathrm{𝟎}𝛀_𝒛\overline{𝒗}_𝒙`$. Here, $`\overline{𝒗}_𝒙`$ and $`\overline{𝒗}_𝒚`$ denote the horizontal averages of the $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ components of the velocity (note that $`\overline{𝒗}_𝒛`$ is zero).

We now turn our attention to the non-dimensional formulation of the anelastic equations as implemented in the code. As was described in AFF, many of the dimension-less physical variables and fundamental physical parameters are expressed in terms of the magnetic field strength along the axis of the initial flux tube, $`𝑩_\mathrm{𝟎}`$. With the addition of the Coriolis force to the momentum equation, we introduce another fundamental dimension-less quantity, $`𝑹_𝒃\mathbf{}𝒗_𝑨\mathbf{/}\mathbf{(}\mathrm{𝟐}𝛀𝑯_𝒓\mathbf{)}\mathbf{}\mathrm{𝟏}\mathbf{/}\mathbf{(}\mathrm{𝟐}𝛀^{\mathbf{}}\mathbf{)}`$, which we refer to as the “magnetic Rossby number” (a quantity also discussed by Schüssler & Solanki 1992). Here, $`𝛀^{\mathbf{}}`$ is defined as the non-dimensional rotation rate given in terms of the ratio of the characteristic velocity scale to the characteristic length scale of our problem; namely, the pressure scale height at the base of the computational domain ($`𝑯_𝒓`$), and the Alfvén speed along the axis of the initial tube ($`𝒗_𝑨\mathbf{}𝑩_\mathrm{𝟎}\mathbf{/}\mathbf{(}\mathrm{𝟐}\sqrt{𝝅𝝆_𝒓}\mathbf{)}`$, where the characteristic density, $`𝝆_𝒓`$, is taken to be the density at the lower boundary). The magnetic Rossby number is a measure of the relative importance of the Coriolis force on the dynamic evolution of our buoyant magnetic flux tube. Rotational effects become important when $`𝑹_𝒃\mathbf{<}\mathrm{𝟏}`$, and are less significant when $`𝑹_𝒃\mathbf{>}\mathrm{𝟏}`$. We therefore wish to investigate three representative cases: $`𝑹_𝒃\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$, $`𝑹_𝒃\mathbf{=}\mathrm{𝟏}`$, and $`𝑹_𝒃\mathbf{=}\mathrm{𝟐}`$. Of course, the solar rotation rate is not a free parameter — estimates of $`𝛀`$ in different parts of the convection zone can be obtained from helioseismologic measurements. Thus, varying the magnetic Rossby number to investigate the effect of solar rotation on a buoyant magnetic flux tube amounts to varying the initial axial field strength of the tube.

The specification of $`𝑩_\mathrm{𝟎}`$ impacts the scaling of two of the non-dimensional diffusive parameters: the Reynolds number ($`𝑹_𝒆\mathbf{}\mathbf{[}𝝆\mathbf{]}\mathbf{[}𝒍\mathbf{]}\mathbf{[}𝒗\mathbf{]}𝝁^\mathbf{}\mathrm{𝟏}`$), and the magnetic Reynolds number ($`𝑹_𝒎\mathbf{}\mathbf{[}𝒍\mathbf{]}\mathbf{[}𝒗\mathbf{]}𝜼^\mathbf{}\mathrm{𝟏}`$). Again, the characteristic length, velocity, and density scales ($`\mathbf{[}𝒍\mathbf{]}`$, $`\mathbf{[}𝒗\mathbf{]}`$, and $`\mathbf{[}𝝆\mathbf{]}`$) are given in terms of $`𝑯_𝒓`$, $`𝒗_𝑨`$, and $`𝝆_𝒓`$ respectively. For this study, we hold the coefficients of viscosity and magnetic diffusion ($`𝝁`$ and $`𝜼`$ respectively) constant between runs (where possible), and vary $`𝑩_\mathrm{𝟎}`$ such that we obtain our three representative values of the magnetic Rossby number. In this manner, we focus our attention on the effects of solar rotation on the dynamic evolution of a rising $`𝛀`$-loop. We note however, that as $`𝑩_\mathrm{𝟎}`$ becomes large, it is increasingly difficult to hold the coefficient of viscosity constant for a given grid resolution. This is because as the Reynolds number increases, features develop at the forefront of the rising tube and in the trailing eddies that are on the order of a few grid zones in size. If this occurs, the code becomes unstable, and the grid resolution must be increased to alleviate the problem. Table 1 lists the values of the relevant physical parameters for each of seven simulations. For the cases where $`𝑩_\mathrm{𝟎}\mathbf{}\mathrm{𝟏𝟎}^\mathrm{𝟓}`$G, values of $`𝝁`$ and $`𝜼`$ were increased to avoid instabilities. We otherwise would have to double the grid resolution in the tube cross-section to take full advantage of the FFT algorithms employed in the code, and such a large increase is unnecessary for the purposes of this study.

Our labeling convention (see column one of Table 1) consists of four letters which signify both the strength of $`𝑩_\mathrm{𝟎}`$, and the latitude of the computational box. For example, LFHL refers to the run which has an initial magnetic flux tube with a relatively *L*ow axial *F*ield strength positioned at a *H*igh *L*atitude, while HFLL refers to the run with an initial tube of relatively *H*igh axial *F*ield strength positioned at *L*ow *L*atitude — and so on. Run SS0 refers to the simulation where solar rotation was neglected (this notation was chosen in order to be consistent with the naming convention of AFF). In each case, we begin with a static, horizontal, cylindrical magnetic flux tube positioned near the bottom of a simulation box that spans 5.147 pressure scale heights vertically. The initial field configuration is given by $`𝐁\mathbf{=}𝑩_𝒙\mathbf{(}𝒓\mathbf{)}\widehat{𝒙}`$, where $`𝑩_𝒙\mathbf{(}𝒓\mathbf{)}\mathbf{}𝑩_\mathrm{𝟎}𝒆^{\mathbf{}𝒓^\mathrm{𝟐}\mathbf{/}𝒂^\mathrm{𝟐}}`$ and $`𝒓`$ is the radial distance to the central axis in the tube’s cross-section. The size of the initial tube, $`𝒂`$, is taken to be $`\mathbf{0.1}𝑯_𝒓`$ in each of the simulations. We choose not to include any initial field line twist about the central axis of the tube, since we wish to determine what role the Coriolis force *alone* plays in the overall fragmentation of the tube as it rises toward the photosphere. As was done in AFF, we apply an artificial entropy perturbation to the midsection of the tube at time $`𝒕\mathbf{=}\mathrm{𝟎}`$. This causes the tube to rise as an $`𝛀`$-loop without having to take the computationally expensive step of evolving each case self-consistently from a state of initial force balance (as was done in Caligari et al. 1995; Fan & Fisher 1996; Fan 2000).

## 3 Results

Each simulation begins with an untwisted cylindrical magnetic flux tube embedded near the base of our model convection zone. An entropy perturbation is applied to the flux tube such that in the absence of rotation, the tube will rise as an $`𝛀`$-loop, and the apex will approach the photospheric boundary in a highly “fragmented” state. That is, a majority of the magnetic flux will be concentrated away from the tube’s central axis along vortex pairs that are shed during the tube’s ascent. Following AFF, we consider a section of the $`𝛀`$-loop to be fragmented if the ratio of the magnetic field weighted second moments of position along the Frenet binormal and normal directions of the tube exceed $`\mathbf{1.5}`$ (see Section 3 of AFF for details on how tube fragmentation is determined, and see Figures 1 and 3 in AFF for visual examples of a fragmented and non-fragmented $`𝛀`$-loop). We expect that solar rotation will play a dynamically important role during the tube’s ascent, since the loop’s rise time is on the order of a rotation period for reasonable initial field strengths ($`\mathrm{𝟑}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$G $`\mathbf{}𝑩_\mathrm{𝟎}\mathbf{}\mathrm{𝟏}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟓}`$G, see Fisher et al. 2000).

We find that the impact of the Coriolis force on a rising flux tube is dramatic, as demonstrated in Figure 1. In this figure, we compare the magnetic field strength in the apex cross-section of two $`𝛀`$-loops: one which rises at $`\mathrm{𝟏𝟓}^{\mathbf{}}`$ latitude in the presence of solar rotation ($`𝛀\mathbf{}\mathbf{13.9}^{\mathbf{}}`$/day) with an initial axial field strength of $`𝑩_\mathrm{𝟎}\mathbf{=}\mathbf{2.67}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$G (corresponding to $`𝑹_𝒃\mathbf{=}\mathbf{0.5}`$), and an identical loop which rises in the absence of rotation. Two effects are immediately apparent. First, the $`𝛀`$-loop that rises through the rotating atmosphere has yet to attain the height of the loop that rises through the non-rotating atmosphere; and second, the loop subjected to the effects of the Coriolis force retains its cohesion, and does not fragment. Figure 2 is a volume rendering of a rising magnetic flux tube in a run where $`𝑹_𝒃\mathbf{=}\mathrm{𝟐}`$ and $`𝑩_\mathrm{𝟎}\mathbf{=}\mathbf{1.0}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟓}`$G. The coordinate axes are chosen such that $`\widehat{𝒚}`$ points northward, and $`\widehat{𝒛}`$ points radially outward. Again, we see that the tube retains its cohesion during its ascent. Figure 2 also demonstrates that in the presence of rotation, the $`𝛀`$-loop is somewhat asymmetric about the apex. That is, the trailing portion of the loop forms a steeper angle with respect to the $`\widehat{𝒙}`$ direction than the leading side of the loop.

A qualitative understanding of these effects can be obtained by examining Figure 3, which shows $`\mathbf{|}𝐁\mathbf{|}`$, the transverse components of the Coriolis force, and the flow field in the apex cross-section during the initial stages of run LFLL (an earlier snapshot of the same run that is displayed in the left frame of Figure 1). In this plot (and all others, unless otherwise stated) the physical quantities are displayed in naturalized units: the magnetic field strength is given in terms of the strength of the field along the axis of tube at $`𝒕\mathbf{=}\mathrm{𝟎}`$ ($`𝑩_\mathrm{𝟎}`$), the velocities are given in terms of the Alfvén speed along the axis of the initial tube ($`𝒗_𝑨`$), and the unit of time is the ratio of the pressure scale height at the base of the computational domain to the characteristic velocity, $`𝑯_𝒓\mathbf{/}𝒗_𝑨`$. For run LFLL, shown in Figure 3, $`𝑩_\mathrm{𝟎}\mathbf{=}\mathbf{2.67}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$G, $`𝒗_𝑨\mathbf{=}\mathbf{1.68}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$ cm s$`^\mathbf{}\mathrm{𝟏}`$ (for $`𝝆_𝒓\mathbf{=}\mathbf{0.2}`$ g cm$`^\mathbf{}\mathrm{𝟑}`$), and $`𝒕\mathbf{=}\mathbf{1.5}`$ corresponds to $`\mathbf{6.18}`$ days. The axes in the figure are given in terms of computational zones; in this case, each individual zone is $`\mathbf{0.017}𝑯_𝒓`$, or approximately $`\mathrm{𝟏}`$ Mm.

In the right-hand frame of Figure 3 we see that a strong, retrograde axial flow has developed as a result of the Coriolis force acting on rising material within the buoyant loop. This is in sharp contrast to the simulations in AFF, where solar rotation is neglected and the initial flux tube has little or no field line twist. In these cases, the strong axial flow is absent, and counter-rotating vortex pairs can quickly form on either side of the tube’s central axis (see also Wissink et al. 2000a; Dorch & Nordlund 1998). As was shown in AFF and Longcope et al. (1996), the non-vertical component of the hydrodynamic “lift” force per unit volume generated by flows in the vortex tubes acts to push them apart, ultimately resulting in fragmentation near the apex of the $`𝛀`$-loop. However, in the presence of rotation, we find that the circulation (and thus the degree of fragmentation of the loop) is substantially reduced. Roughly speaking, this follows from the fact that the Coriolis force can do no work (it always acts in a direction perpendicular to the flow); it simply redistributes kinetic energy between transverse and axial flows. The initial increase in the axial flow at the apex of the loop gives rise to a secondary effect: a transverse component of the Coriolis force normal to the angular velocity vector. At low latitudes, this component acts primarily in the $`\mathbf{}`$$`\widehat{𝒛}`$ direction, and opposes the buoyant rise of the loop; however, it also forces the top of the loop to drift northward over time.

In simulations where the effects of rotation are neglected, we find that material flows along the rising tube away from the apex. In the presence of rotation, this type of “draining” flow manifests itself as an asymmetry in the axial flow induced by the Coriolis force: the retrograde flow is stronger in the trailing side of the loop than it is in the leading side. As the Coriolis force acts upon these asymmetric axial flows, and as the buoyant, magnetized plasma attempts to preserve its angular momentum, the morphological asymmetry predicted by Moreno-Insertis et al. (1994) and apparent in Figure 2 quickly develops.

Figure 4 shows the transverse components of the Coriolis force, and the flow field in the apex cross-section of run MFLL at $`𝒕\mathbf{=}\mathbf{5.5}`$ (the same time, in naturalized units, as the cross-sections of Figure 1). In this run, $`𝑹_𝒃\mathbf{=}\mathrm{𝟏}`$, and the domain is again centered at $`\mathrm{𝟏𝟓}^{\mathbf{}}`$ latitude. This implies that $`𝑩_\mathrm{𝟎}\mathbf{=}\mathbf{5.35}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$G, $`𝒗_𝑨\mathbf{=}\mathbf{3.37}\mathbf{\times }\mathrm{𝟏𝟎}^\mathrm{𝟒}`$cm s$`^\mathbf{}\mathrm{𝟏}`$, and the unit of time is 2.06 days ($`𝒕\mathbf{=}\mathbf{5.5}`$ thus corresponds to $`\mathbf{11.3}`$ days). In an earlier snapshot of run MFLL ($`𝒕\mathbf{=}\mathbf{1.5}`$), the morphology and flow patterns are very similar to those shown for run LFLL in Figure 3. However, as the runs progress, differences in secondary and tertiary flows become more pronounced, and the loops begin to evolve in a dissimilar manner. This is evident when comparing the first frame of Figure 4 with the frames of Figure 1, which highlight the striking differences in morphology that ultimately develop between simulations with different values of the magnetic Rossby number. These figures show that as the initial axial magnetic field strength of the tube is increased, the tendency for the tube to fragment as it rises toward the photosphere is also increased.

### 3.1 An Approximate Model

To better understand the relationship between magnetic Rossby number, latitude of flux emergence, and loop fragmentation, we approximate the forcing terms of the momentum equation \[equation (1)\] by constant, order-of-magnitude estimates evaluated for our initial conditions, and valid for the initial stages of evolution. We then have a simplified expression for the velocity field that admits to an analytic solution:

$$\frac{𝒅𝐯}{𝒅𝒕}\mathbf{=}\mathbf{}\mathrm{𝟐}\mathbf{\left(}𝛀\mathbf{\times }𝐯\mathbf{\right)}\mathbf{+}𝐐\mathbf{.}$$
(5)
It is important to recognize that this simple approximation will not capture the factor of two reduction in acceleration that occurs in the non-rotating case because of “added inertia” effects (see Longcope et al. 1996). Nevertheless, it does capture the essence of the dynamic effects introduced by the Coriolis force. Let $`\mathbf{(}𝒒_𝒙\mathbf{,}𝒒_𝒚\mathbf{,}𝒒_𝒛\mathbf{)}`$ represent the Cartesian components of the constant driving term $`𝐐`$, and let $`𝜽`$ represent solar latitude. Then the general solution to the above equation is given by

$`𝒗_𝒙\mathbf{(}𝒕\mathbf{)}`$ $`\mathbf{=}`$ $`{\displaystyle \frac{𝒒_𝒙}{\mathrm{𝟐}𝛀}}\mathrm{𝐬𝐢𝐧}\mathrm{𝟐}𝛀𝒕\mathbf{}𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}\mathrm{𝐜𝐨𝐬}\mathrm{𝟐}𝛀𝒕\mathbf{+}𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}`$ (6)
$`𝒗_𝒚\mathbf{(}𝒕\mathbf{)}`$ $`\mathbf{=}`$ $`{\displaystyle \frac{𝒒_𝒙}{\mathrm{𝟐}𝛀}}\mathrm{𝐬𝐢𝐧}𝜽\mathrm{𝐜𝐨𝐬}\mathrm{𝟐}𝛀𝒕\mathbf{+}𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}\mathrm{𝐬𝐢𝐧}𝜽\mathrm{𝐬𝐢𝐧}\mathrm{𝟐}𝛀𝒕`$ (7)
$`\mathbf{+}\mathbf{(}𝒒_𝒚\mathbf{}\mathrm{𝟐}𝛀𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}\mathrm{𝐬𝐢𝐧}𝜽\mathbf{)}𝒕\mathbf{}{\displaystyle \frac{𝒒_𝒙}{\mathrm{𝟐}𝛀}}\mathrm{𝐬𝐢𝐧}𝜽`$
$`𝒗_𝒛\mathbf{(}𝒕\mathbf{)}`$ $`\mathbf{=}`$ $`\mathbf{}{\displaystyle \frac{𝒒_𝒙}{\mathrm{𝟐}𝛀}}\mathrm{𝐜𝐨𝐬}𝜽\mathrm{𝐜𝐨𝐬}\mathrm{𝟐}𝛀𝒕\mathbf{}𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}\mathrm{𝐜𝐨𝐬}𝜽\mathrm{𝐬𝐢𝐧}\mathrm{𝟐}𝛀𝒕`$ (8)
$`\mathbf{+}\mathbf{(}𝒒_𝒛\mathbf{+}\mathrm{𝟐}𝛀𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}\mathrm{𝐜𝐨𝐬}𝜽\mathbf{)}𝒕\mathbf{+}{\displaystyle \frac{𝒒_𝒙}{\mathrm{𝟐}𝛀}}\mathrm{𝐜𝐨𝐬}𝜽\mathbf{,}`$
where

$`𝒗_𝒙^{\mathbf{(}𝑷\mathbf{)}}`$ $`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}𝛀}}\mathbf{(}𝒒_𝒚\mathrm{𝐬𝐢𝐧}𝜽\mathbf{}𝒒_𝒛\mathrm{𝐜𝐨𝐬}𝜽\mathbf{)}\mathbf{.}`$ (9)
The total circulation over a cross-sectional slice of the flux tube is given by $`𝚪\mathbf{=}\mathbf{}_𝑺`$$`𝝎`$$`\mathbf{}\widehat{𝐧}_𝑺𝒅𝑺\mathbf{=}\mathbf{}_𝑪𝐯\mathbf{}𝒅`$$`\mathbf{}`$ , where $`𝑪`$ is the closed circuit bounding the cross-sectional surface $`𝑺`$, $`\widehat{𝐧}_𝑺`$ is the outward normal to $`𝑺`$, and $`𝝎`$$`\mathbf{=}`$$`\mathbf{}`$$`\mathbf{\times }𝐯`$ is the vorticity vector. As the tube begins to rise, oppositely directed circulation patterns form on either side of a line defined by the Frenet normal. The total circulation of these regions (at the loop apex) is given by $`𝚪\mathbf{=}\mathbf{}_𝑪𝐯\mathbf{}𝒅`$$`\mathbf{}`$$`\mathbf{}\mathbf{}_\mathbf{}_𝒂𝒅\mathbf{}_𝒂𝒗_𝒛\mathbf{}𝒂\mathbf{<}𝒗_𝒛\mathbf{>}`$, if we assume that the primary contribution to $`𝚪`$ stems from $`𝐯\mathbf{}𝒅`$$`\mathbf{}`$ across the diameter of the tube, $`𝒂`$. Since $`𝒗_𝒛`$ is roughly constant across the interior of the tube, we can approximate the circulation (and thus the forces acting to pull the tube apart) with our analytic solution for the velocity. Initially, the motion of the tube is determined by the buoyancy force, which acts in the $`\widehat{𝒛}`$ direction. Thus, we assume that $`𝒒_𝒙`$ and $`𝒒_𝒚`$ are zero, and express the circulation as:

$$𝚪\mathbf{(}𝒕\mathbf{)}\mathbf{=}\frac{𝒂𝒗_{𝑨}^{}{}_{}{}^{\mathrm{𝟐}}}{\mathrm{𝟐}𝑯_𝒓}\mathbf{\left(}𝒕\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}𝜽\mathbf{+}\frac{\mathrm{𝐬𝐢𝐧}\mathbf{(}\mathrm{𝟐}𝛀𝒕\mathbf{)}}{\mathrm{𝟐}𝛀}\mathrm{𝐜𝐨𝐬}^\mathrm{𝟐}𝜽\mathbf{\right)}\mathbf{.}$$
(10)
We plot our simple approximation of $`𝚪\mathbf{(}𝒕\mathbf{)}`$ for four representative values of $`𝑹_𝒃`$ in Figure 5.

In order to compare our approximation with the MHD solution, we must first obtain $`𝚪\mathbf{(}𝒕\mathbf{)}`$ directly from the simulation data. To accomplish this, we follow the methodology described in Section 3 of AFF, and first calculate the path of the magnetic flux tube through the background plasma. The position of the loop (and its fragments, if applicable) is specified by the magnetic field weighted first moments of position along a series of two-dimensional slices that closely correspond to the tube’s cross-sectional plane. Once the path of the tube is determined (see Section 3 of AFF), we construct the Frenet basis vectors, and define two regions on either side of a line defined by the Frenet normal where the oppositely directed circulation patterns will develop. From the simulation results, we then compute the total circulation across each of these regions. Figure 6 shows $`𝚪\mathbf{(}𝒕\mathbf{)}`$ for each simulation. One can compare these results with the simple approximation shown in Figure 5. It is clear that the general pattern of time dependence of the circulation is successfully reproduced by our simple model, though there are significant differences in the details.

The oscillations in $`𝚪\mathbf{(}𝒕\mathbf{)}`$ are inertia waves (see Batchelor 1967) with a period of $`𝒕\mathbf{=}\mathrm{𝟐}𝝅𝑹_𝒃`$ (here, $`𝒕`$ is expressed in normalized units of $`𝑯_𝒓\mathbf{/}𝒗_𝑨`$). Thus, the tendency for the loops to fragment is significantly reduced after the circulation begins to fall off at $`𝒕\mathbf{=}𝝅𝑹_𝒃\mathbf{/}\mathrm{𝟐}`$. The overall evolution of the loop depends strongly on the dynamics during the initial portion of its rise, since the early differences in secondary and tertiary flows induced by the Coriolis force have a significant impact on the eventual evolution of the tube. We find that if the tube is able to rise several of its diameters before $`𝒕`$ reaches $`𝝅𝑹_𝒃\mathbf{/}\mathrm{𝟐}`$, then the resulting $`𝛀`$-loop will ascend toward the photosphere in a fragmented state — otherwise, the tube evolves cohesively.

We can infer from our approximate model that, on average, the rise speed of the tube slowly increases over time. This prediction is confirmed by the simulations, which (for sufficiently high choices of $`𝑹_𝒆`$ and $`𝑹_𝒎`$) show the initial slow, linear increase in the rise speed of the tube, modulated by inertial oscillations. This behavior persists in the simulations until $`𝒕`$ becomes comparable to the diffusive timescale. If this occurs, magnetic diffusion begins to (artificially) reduce the buoyancy of the tube, eventually preventing it from approaching the upper boundary (this occurred in runs LFLL and LFHL after $`𝒕\mathbf{}\mathbf{34.0}`$). Of course, this problem can be mitigated by reducing the coefficient of magnetic diffusivity in the simulations, but this requires a significant increase in grid resolution (to avoid numerical instabilities); and since our analysis concentrates on the early stages of loop evolution, such a computationally expensive step is unwarranted.

### 3.2 Comparison with Two-dimensional Models

To investigate the role that three-dimensional geometry plays in the tube’s evolution, we compare our runs with simulations of buoyant, axially symmetric magnetic flux tubes moving in two dimensions. Figure 7 shows the tube cross-section of a 2-D run whose parameters are that of the 3-D run MFLL shown in the first frame of Figure 4. This 2-D case can be thought of as an infinitely-long, buoyant, axially symmetric $`𝛀`$-loop of zero apex curvature. A comparison of the apex cross-section shown in Figure 4 with the cross-section of Figure 7 reveals that a 2-D geometric constraint imposed on the solution (even in the early stages of loop evolution) significantly alters the dynamic evolution of the tube. This suggests that in the presence of solar rotation, the geometry of the $`𝛀`$-loop still plays an important role in determining the morphology of the emerging magnetic field — a result consistent with the findings of AFF for similar runs where solar rotation was neglected.

In 2-D, the Coriolis force still acts to suppress the fragmentation of the tube. However, the hydrodynamic “lift” force that acts to fragment the tube is overestimated due to the tube’s effectively infinite axial length scale. It is the binormal component of this force that acts to pull the tube apart (see Longcope et al. 1996; Fan et al. 1998; AFF), thus the overall tendency for the tube to retain its cohesion during its ascent is reduced as compared to the 3-D cases. In addition, the amount that the tube is deflected toward the pole is significantly overestimated in a 2-D geometry. This is because in 3-D, the component of the Coriolis force that acts to push the tube poleward acts only on the portion of the $`𝛀`$-loop where there are strong, mainly axial flows (eg. near the apex). In general, our calculations support the conclusions of Wissink et al. (2000a) who used 2-D simulations to demonstrate that in the presence of the Coriolis force, flux tubes remained somewhat more cohesive during their ascent. We note, however, that in 3-D, this effect is more dramatic.

### 3.3 Morphology of the Emerging Magnetic Field

There are several reasons why we must be cautious when making predictions about the behavior of emerging magnetic flux based on data obtained from our simulations. First, the anelastic approximation becomes marginal as the flux tube approaches the photospheric boundary, where densities are lower, and the acoustic Mach number of the flow approaches unity. Second, we intended for these simulations to investigate the effects of the Coriolis force alone on the evolution of our $`𝛀`$-loop. Thus, we did not include any field line twist along our initial tube, nor did we consider the effects of a tube rising through a dynamically convecting background state. Third, our models do not account for radiative losses that occur at the surface. However, given these caveats, we feel that we can still make some general predictions regarding expected magnetic field and velocity flow patterns if the magnetic flux emerges through a relatively uncluttered region on the solar disk.

Figure 8 shows the radial component of the magnetic field along horizontal slices taken underneath the apex of the rising $`𝛀`$-loops of runs HFLL, and SS0. This figure emphasizes the differences in morphology of emerging magnetic field between a simulation that includes the effects of solar rotation (run HFLL), and one that does not (run SS0). In the figures, outwardly directed field is displayed in white, and inwardly directed field is displayed in black. Superimposed over the each grey-scale image is the transverse velocity field. One of the most striking effects of the Coriolis force is the strong flow that rapidly develops in both the magnetized and unmagnetized plasma directed from the leading toward the trailing polarity (the leading polarity moves in the direction of solar rotation as the loop continues to rise). This flow is simply the result of the Coriolis force acting on the rising fluid just below the apex of each loop. This “LTT” (*L*eading *T*o *T*railing polarity) flow is evident in all of our simulations except in cases where solar rotation is neglected.

The Coriolis force manifests itself in several other interesting ways. First, since the trailing edge of the loop emerges at a steeper angle with respect to the surface than the leading edge (an effect most pronounced in the runs where $`𝑹_𝒃`$ is small), the leading polarity becomes “spread-out” in the east-west direction (along $`𝒙`$). Second, in the cases where $`𝑹_𝒃`$ is small, we observe that the leading polarity is positioned slightly closer to the equator than the trailing polarity — behavior that is roughly consistent with the “Joy’s Law” result from thin flux tube models (see D’Silva & Choudhuri 1993; Fan et al. 1994; Schüssler et al. 1994; Caligari et al. 1995; Fan & Fisher 1996; Fisher et al. 2000). Third, the north-south dispersion of the emerging flux is greatly reduced when solar rotation is included in the simulations. This reflects the increased cohesion of the flux rope during its ascent through our model convection zone.

## 4 Summary and Conclusions

AFF found that the degree of fragmentation of a rising $`𝛀`$-loop depends upon the 3-D geometry of the loop — the greater the apex curvature, the lesser the degree of fragmentation for a fixed amount of initial field-line twist. We are now able to extend this analysis to a rotating model convection zone, and investigate how the degree of apex fragmentation depends on the axial field strength of the initial magnetic flux tube. Assuming a given rotation rate, we find that:

1. A magnetic flux tube that rises buoyantly toward the photosphere in a rotating convection zone rises more slowly, and is better able to retain its cohesion than an identical flux tube that rises through a non-rotating convection zone — even in the absence of initial field line twist. If the tube is sufficiently buoyant so that it is able to rise a distance of several tube diameters before $`𝒕\mathbf{=}𝝅\mathbf{/}\mathrm{𝟒}𝛀`$ (where $`𝛀`$ is an estimate of the relevant solar rotation rate at a given latitude), then the tube will fragment during its rise toward the photosphere.
2. In the absence of forces due to convective motions, the stronger the initial axial field strength of the tube, the more fragmentation occurs as the tube rises.
3. An approximate model is able to predict the initial evolution of the circulation, and hence can effectively characterize the formation and subsequent interaction of the oppositely-directed vortex pairs of a fragmented $`𝛀`$-loop.
4. One must be careful not to over-interpret results from two-dimensional simulations of buoyant magnetic flux tubes. In a two-dimensional geometry, hydrodynamic forces due to vortex interaction are overestimated. The amount of poleward deflection due to the Coriolis force is also overestimated.
5. If magnetic flux emerges as an $`𝛀`$-loop, the Coriolis force will induce a relatively strong flow that is directed from the region of leading polarity to the region of trailing polarity. This flow is present in both the magnetized and unmagnetized plasma.

This work was funded by NSF grants AST 98-19727 and ATM 98-96316, by NASA grant NAG5-8468, and by ONR grant FDN00173-00-1-G901-03/03. The computations described here were supported by the National Center for Atmospheric Research under grant ATM 98-96316. Further computational support was provided by the National Computational Science Alliance. We would like to thank the authors of the FFTW package (Frigo & Johnson, 1997) for making their code publicly available.


text:
# Correlated hopping of electrons: Effect on the Brinkman-Rice transition and the stability of metallic ferromagnetism

## Abstract

We study the Hubbard model with bond-charge interaction (“correlated hopping”) in terms of the Gutzwiller wave function. We show how to express the Gutzwiller expectation value of the bond-charge interaction in terms of the correlated momentum-space occupation. This relation is valid in all spatial dimensions. We find that in infinite dimensions, where the Gutzwiller approximation becomes exact, the bond-charge interaction lowers the critical Hubbard interaction for the Brinkman-Rice metal-insulator transition. The bond-charge interaction also favors ferromagnetic transitions, especially if the density of states is not symmetric and has a large spectral weight below the Fermi energy.

The microscopic origin of magnetic ordering in systems like transition metals, transition-metal oxides, and high-temperature superconductors is intricate, since it is due to correlations between the electrons. The simplest model to attempt a description of such systems is the single-band Hubbard model

$`\widehat{H}_{\text{Hubbard}}`$ $`=`$ $`{\displaystyle \underset{ij\sigma }{}}t_{ij}\widehat{c}_{i\sigma }^+\widehat{c}_{j\sigma }^{}+U{\displaystyle \underset{i}{}}\widehat{n}_i\widehat{n}_i,`$ (1)
where hats indicate operators. This model describes the competition between kinetic and potential energy which is at the heart of the quantum-mechanical correlation problem. The Hubbard interaction $`U`$ represents the Coulomb repulsion of electrons in the same orbital at a given lattice site. It is given by the matrix element $`U`$ $`=`$ $`ii|V(𝒓𝒓^{\mathbf{}})|ii`$ of the Coulomb potential, and is typically on the order of a few eV. The matrix elements involving neighboring lattice sites $`i`$ and $`j`$ are generally smaller than $`U`$, but may not be negligibly small. One of them is the bond-charge interaction $`X_{ij}`$ $`=`$ $`ii|`$ $`V(𝒓𝒓^{\mathbf{}})`$ $`|ij`$, which is typically on the order of 0.1-1 eV and hence is comparable in magnitude to the tight-binding hopping amplitude $`t_{ij}`$. It describes a density-dependent hopping of the electrons

$`\widehat{H}_X`$ $`=`$ $`{\displaystyle \underset{ij\sigma }{}}X_{ij}\widehat{c}_{i\sigma }^+\widehat{c}_{j\sigma }^{}(\widehat{n}_{i\overline{\sigma }}+\widehat{n}_{j\overline{\sigma }}),`$ (2)
which only contributes if the lattice site from or onto which an electron with spin $`\sigma `$ is hopping is occupied by an electron with spin $`\overline{\sigma }`$ (“correlated hopping”). The effect of $`\widehat{H}_X`$ competes strongly with both the kinetic energy and the Hubbard interaction. Correlated hopping of spin-$`\sigma `$ electrons between two sites is enhanced if spin-$`\overline{\sigma }`$ electrons are present, but this in turn will cost the latter kinetic energy, as well as Coulomb energy for the double occupations. Moreover, in a band picture the coupling of densities and kinetic energy can lead to a band narrowing which lowers the amount of energy that is necessary for a ferromagnetic spin polarization. These mechanisms explain why, in principle, $`\widehat{H}_X`$ may play an important role in the stabilization of ferromagnetism and the localization of electrons.

The model $`\widehat{H}`$ $`=`$ $`\widehat{H}_{\text{Hubbard}}`$ $`+`$ $`\widehat{H}_X`$ cannot be solved exactly, and only mean-field and finite-size diagonalization results are available. Recently, Schiller showed how to incorporate $`\widehat{H}_X`$ into the framework of dynamical mean-field theory, but no numerical results have yet been obtained for $`\widehat{H}`$.

One of the standard tools to approach the correlated electron problem is the Gutzwiller wave function

$$|\mathrm{\Psi }_\text{G}=\underset{i}{}[1(1g)\widehat{n}_i\widehat{n}_i]|\mathrm{\Phi }_0,$$
(3)
where $`g`$ is a variational parameter ($`0g1`$), and the starting wave function $`|\mathrm{\Phi }_0`$ is a product state of spin-up and spin-down Fermi seas. By construction both $`|\mathrm{\Phi }_0`$ and $`|\mathrm{\Psi }_\text{G}`$ are translationally invariant and have a fixed particle density $`n`$ $`=`$ $`n_{}`$ $`+`$ $`n_{}`$ and magnetization $`m`$ $`=`$ $`n_{}`$ $``$ $`n_{}`$. The uncorrelated case $`U=X_{ij}=0`$ corresponds to $`g=1`$, while $`U=\mathrm{}`$ forbids any doubly occupied sites and thus corresponds to $`g=0`$. For $`m`$ $``$ $`0`$ and $`g0`$ the wave function $`|\mathrm{\Psi }_\text{G}`$ describes an itinerant ferromagnetic state. Starting wave functions $`|\mathrm{\Phi }_0`$ with other broken symmetries can also be considered. Here, however, we will only consider paramagnetism and ferromagnetism.

Using the Gutzwiller wave function one may, in principle, calculate expectation values of any operator $`\widehat{A}`$ as $`\widehat{A}_\text{G}`$ $`=`$ $`\mathrm{\Psi }_\text{G}|\widehat{A}|\mathrm{\Psi }_\text{G}/\mathrm{\Psi }_\text{G}|\mathrm{\Psi }_\text{G}`$. The energy expectation value $`E`$ $`=`$ $`\widehat{H}_\text{G}`$, when optimized with respect to $`g`$, is an upper bound for the exact ground-state energy of $`\widehat{H}`$ by the variational principle. The variational energy $`E`$ can be written as

$$E=\widehat{H}_\text{G}=\underset{\text{k}\sigma }{}ϵ_\text{k}n_{\text{k}\sigma }^{}+\underset{\text{k}\sigma }{}2\xi _\text{k}\text{Re}(x_{\text{k}\sigma }^{})+Ud,$$
(4)
where the $`k`$-space occupation $`n_{𝒌\sigma }^{}`$, bond-charge occupation $`x_{𝒌\sigma }^{}`$, and double occupation $`d`$ are defined by ($`L`$ is number of lattice sites)

$`n_{𝒌\sigma }^{}`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{ij}{}}e^{i𝒌(𝑹_i𝑹_j)}\widehat{c}_{i\sigma }^+\widehat{c}_{j\sigma }^{}_\text{G},`$ (5)
$`x_{𝒌\sigma }^{}`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{ij}{}}e^{i𝒌(𝑹_i𝑹_j)}\widehat{n}_{i\overline{\sigma }}\widehat{c}_{i\sigma }^+\widehat{c}_{j\sigma }^{}_\text{G},`$ (6)
$`d`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{i}{}}\widehat{n}_i\widehat{n}_i_\text{G}.`$ (7)
They each depend on $`g`$, $`n`$, and $`m`$. Here $`ϵ_𝒌`$ and $`\xi _𝒌`$ are the Fourier transforms of $`t_{ij}`$ and $`X_{ij}`$, respectively. By convention, $`t_{ii}`$ $`=`$ $`X_{ii}`$ $`=`$ $`0`$.

In general the expectation values in Eqs. (5)-(7) are not independent of one another. In particular, for the Gutzwiller wave function the $`k`$-space and bond-charge occupation are closely related. Using the techniques of Ref.  we obtain, for all Bravais lattices in arbitrary dimensions $`D`$,

$$x_{\text{k}\sigma }^{}=\frac{[(1+g)n_{\text{k}\sigma }^0g]n_{\text{k}\sigma }^{}n_{\text{k}\sigma }^0}{1g}+n_{\overline{\sigma }}n_{\text{k}\sigma }^0d,$$
(8)
with $`n_{𝒌\sigma }^0`$ $`=`$ $`n_{𝒌\sigma }^{}|_{g=1}`$ as the uncorrelated Fermi function. It should be noted that although $`x_{𝒌\sigma }^{}`$ and $`n_{𝒌\sigma }^{}`$ are linearly related, the bond-charge energy and the kinetic energy will generally not be, since the regions of the Brillouin zone inside ($`n_{𝒌\sigma }^0=1`$) and outside ($`n_{𝒌\sigma }^0=0`$) of the Fermi sea contribute differently.

Gutzwiller approximation.— Now we turn to evaluate the variational energy $`E`$ within the Gutzwiller approximation, which is known to yield the exact evaluation of expectation values in terms of the Gutzwiller wave function in the limit of infinite spatial dimensions ($`D`$ $``$ $`\mathrm{})`$. It describes a Fermi liquid with piecewise constant $`k`$-space occupation,

$$n_{\text{k}\sigma }^{}=n_\sigma +(n_{\text{k}\sigma }^0n_\sigma )q_\sigma ,$$
(9)
where the discontinuity at the Fermi surface is given by

$$q_\sigma =\frac{[\sqrt{(n_\sigma d)(1n+d)}+\sqrt{(n_{\overline{\sigma }}d)d}]^2}{n_\sigma (1n_\sigma )}.$$
(10)
The variational parameter $`g`$ is related to the double occupation $`d`$ by

$$g^2=\frac{(n_{}d)(n_{}d)}{(1n+d)d},$$
(11)
and it is convenient to use the latter as variational parameter. For the variational energy we obtain

$$E=E^0\underset{\sigma }{}(1q_\sigma )(ϵ_{0\sigma }+\xi _{0\sigma })+U_{\text{eff}}(dn_{}n_{}),$$
(12)
where $`ϵ_{0\sigma }`$ $`=`$ $`(1/L)_𝒌ϵ_𝒌n_{𝒌\sigma }^0`$ is the uncorrelated kinetic energy, and similarly $`\xi _{0\sigma }`$ $`=`$ $`(1/L)_𝒌\xi _𝒌n_{𝒌\sigma }^0`$. The effective Hubbard interaction $`U_{\text{eff}}`$ and the uncorrelated variational energy $`E^0`$ ($`=`$ $`E|_{g=1}`$) in Eq. (12) are given by

$`U_{\text{eff}}`$ $`=`$ $`U+{\displaystyle \underset{\sigma }{}}{\displaystyle \frac{(12n_\sigma )}{(1n_\sigma )n_\sigma }}\xi _{0\sigma },`$ (13)
$`E^0`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}(ϵ_{0\sigma }+2n_{\overline{\sigma }}\xi _{0\sigma })+Un_{}n_{}.`$ (14)
Within the Gutzwiller approximation the bond-charge interaction thus leaves the form of the variational energy unchanged, but enters into the effective kinetic energy and effective Hubbard interaction via $`\xi _{0\sigma }`$. The effect on the kinetic energy can be interpreted as a spin-dependent band narrowing or widening, which is also found in the Hartree-Fock approximation of the correlated hopping term and can lead to a stabilization of ferromagnetism. However the Gutzwiller approximation reveals two distinct effects that cannot be resolved in ordinary Hartree-Fock theory, where the suppression of double occupancies can only be achieved by spin polarizing the system. On the one hand, the kinetic energy increases if $`\xi _{0\sigma }`$ $`>`$ $`0`$, i. e., the effective hopping becomes smaller. This corresponds to an effective narrowing of the band, which lowers the amount of energy that must be expended for a ferromagnetic spin polarization. Furthermore the bond-charge interaction contributes to $`U_{\text{eff}}`$, which suppresses double occupancies already at lower values of the bare Hubbard interaction $`U`$. The balance between these two effects is determined by the optimal variational parameter $`d`$.

Brinkman-Rice metal-insulator transition.— We now discuss the effect of $`\widehat{H}_X`$ on the Brinkman-Rice transition that occurs in the Gutzwiller approximation at half-filling ($`n=1`$). For convenience we define the strength of the bond-charge interaction $`X`$ by

$$\xi _{0\sigma }=Xϵ_{0\sigma }.$$
(15)
Note that $`\xi _{0\sigma }`$ will remain proportional to $`ϵ_{0\sigma }`$ according to Eq. (15) for all densities if $`t_{ij}`$ and $`X_{ij}`$ have the same range, e. g., if they are nonzero only for nearest-neighbor sites. In the limit of $`D`$ $``$ $`\mathrm{}`$ both $`t_{ij}`$ and $`X_{ij}`$ must both be scaled as $`1/\sqrt{Z_{ij}}`$ ($`Z_{ij}`$ is the number of neighbors $`ij`$), which is compatible with Eq. (15). The dispersion $`ϵ_𝒌`$ enters only through the density of states (DOS) $`N(ϵ)`$ (which determines $`ϵ_{0\sigma }`$), as expected in dimension $`D=\mathrm{}`$. We will consider several densities of states below.

With the above definition of $`X`$, the Gutzwiller approximation energy for $`n=1`$ and $`m=0`$ simplifies to $`E`$ $`=`$ $`8d(12d)(1X)ϵ_0+Ud`$, where $`ϵ_0`$ $``$ $`_\sigma ϵ_{0\sigma }`$ $`<`$ $`0`$. Optimization with respect to $`d`$ yields a critical value for $`U`$,

$$U_\text{c}(X)=8|ϵ_0|(1X),$$
(16)
above which the localized state with $`d=0`$ is lowest in energy. Hence the Brinkman-Rice transition is moved to lower $`U`$ for $`X>0`$, i. e., the bond-charge interaction favors localization. (Only $`U0`$ and $`X1`$ will be considered from now on.) We find that the $`U`$ dependence of the double occupation $`d`$, the discontinuity of the $`k`$-space occupation $`q`$ ($``$ $`q_\sigma `$), and the energy $`E`$ is formally the same as in the original Brinkman-Rice theory for $`X`$ $`=`$ $`0`$, i. e.,

$`d`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(1{\displaystyle \frac{U}{U_\text{c}}}\right),`$ (17)
$`q`$ $`=`$ $`1{\displaystyle \frac{U^2}{U_\text{c}^2}},`$ (18)
$`E`$ $`=`$ $`{\displaystyle \frac{U_\text{c}}{8}}\left(1{\displaystyle \frac{U}{U_\text{c}}}^2\right),`$ (19)
except that $`U_\text{c}`$ now depends on $`X`$ \[Eq. (16)\]. We note that the simultaneous vanishing of $`q`$ and double occupation $`d`$ at a finite value of $`U`$ is characteristic of the Brinkman-Rice transition, in contrast to the numerical solution of the Hubbard model, where $`d`$ remains nonzero across the transition.

Ferromagnetic transition.— The instability of the paramagnetic state toward ferromagnetism can be determined from the bulk susceptibility $`\chi `$. For half-filling, we obtain

$$\frac{1}{\chi }=\frac{q}{2N(ϵ_F)}\left[(1X)\left(1\frac{pU(U+2U_\text{c})}{(U+U_\text{c})^2}\right)+\frac{rXU_\text{c}}{U+U_\text{c}}\right],$$
(20)
where we have introduced the dimensionless parameters

$`p`$ $`=`$ $`4N(ϵ_F)|ϵ_0|,`$ (21)
$`r`$ $`=`$ $`4N(ϵ_F)ϵ_F.`$ (22)
The Fermi energy $`ϵ_F`$ in Eq. (22) represents an absolute scale since the first moment of the DOS is fixed at zero (due to $`t_{ii}=0`$). There are two factors in $`\chi `$ that can diverge: either $`q0`$, i. e., the effective band mass $`m^{}/m=q^1`$ diverges at $`U_\text{c}(X)`$, indicating a localization transition, or the Stoner-type factor in square brackets in Eq. (20) vanishes at

$$U_{\text{fm}}(X)=\frac{\sqrt{r^2X^24p(1p)(1X)^2}rX}{2(1p)}+X1,$$
(23)
signaling an instability toward ferromagnetism. The latter instability precedes the localization transition whenever $`p>p_{\text{fm}}`$, where

$$p_{\text{fm}}=\frac{4}{3}\left(1+\frac{rX}{2(1X)}\right).$$
(24)
These results reduce to the known values $`p_{\text{fm}}`$ $`=`$ $`\frac{4}{3}`$ and $`U_{\text{fm}}`$ $`=`$ $`(\sqrt{p/(p1)}1)`$ $`U_\text{c}`$ for $`X=0`$.

Let us first consider the effect of the bond-charge interaction in the case of a symmetric DOS, $`N(ϵ)`$ $`=`$ $`N(ϵ)`$, which results if hopping takes place only between different sublattices of a bipartite lattice. The Fermi energy at half-filling is then $`ϵ_F`$ $`=`$ $`0`$, hence $`r=0`$. In this case $`p_{\text{fm}}`$ $`=`$ $`\frac{4}{3}`$, the same criterion as for $`X=0`$. On the other hand, for an asymmetric DOS ferromagnetism is favored by $`X>0`$ if $`r<0`$, i. e., $`ϵ_F<0`$. This is the case if the Fermi energy is below the center of mass of the DOS, which means that there is large spectral weight below the Fermi energy. The tendency toward ferromagnetism in such a situation was already proposed long ago.

We now consider nearest-neighbor hopping $`t_{ij}`$ $`=`$ $`t^{}/\sqrt{Z}`$ and bond-charge interaction $`X_{ij}`$ $`=`$ $`X`$ $`t^{}/\sqrt{Z}`$ on several infinite-dimensional lattices. Their densities of states and phase diagram are shown in Fig. 1. The Gaussian DOS of the hypercubic lattice, $`N(ϵ)`$ $`=`$ $`\mathrm{exp}(ϵ^2/2)/\sqrt{2\pi }`$, has a parameter $`p`$ $`=`$ $`1/\pi `$, and the semielliptic DOS of the Bethe lattice, $`N(ϵ)`$ $`=`$ $`\sqrt{4ϵ^2}/2\pi `$, has $`p=32/3\pi ^3`$, and both have $`r=0`$ due to particle-hole symmetry. Since $`p<\frac{4}{3}`$ in both cases, the metal-insulator transition at $`U_\text{c}`$, drawn as a solid line in the phase diagram in Fig. 1, will mask the ferromagnetic phase. The variational phase diagram for $`\widehat{H}_{\text{Hubbard}}`$ (i. e., $`X`$ $`=`$ $`0`$) on the hypercubic lattice was calculated by Fazekas et al., who predicted the ferromagnetic and antiferromagnetic phases to coexist as the system phase separates and to preempt the metal-insulator transition. Here we consider only homogeneous ferromagnetic phases, thus allowing the metal-insulator transition to take place, and do not attempt to distinguish between paramagnetic and ferromagnetic insulators at half-filling, which are degenerate in energy ($`E=0`$).

A more complicated scenario arises if the lattice system is not particle-hole symmetric, so that the DOS is asymmetric and thus $`r0`$. The generalized infinite-dimensional fcc lattice, with hopping scaled as $`t_{ij}`$ $`=`$ $`1/\sqrt{2D(D1)}`$, has a DOS (Ref. )

$$N(ϵ)=\frac{\mathrm{exp}((1+\sqrt{2}ϵ)/2)}{\sqrt{\pi (1+\sqrt{2}ϵ)}},$$
(25)
showing a square-root singularity at the lower band-edge. The Hubbard model (1) on this lattice has been studied numerically by Ulmke within the dynamical mean-field theory, who found ferromagnetism at low enough temperatures and band-filling. For half-filling one has $`ϵ_F`$ $`=`$ $`0.3854`$, $`p`$ $`=`$ $`1.6157`$, and $`r=1.0272`$. Thus ferromagnetism occurs in this case already for $`X=0`$, but the critical $`U`$ is lowered by the presence of $`X>0`$. Hence a ferromagnetic phase is found for $`U_{\text{fm}}(X)`$ $`<`$ $`U`$ $`<`$ $`U_\text{c}(X)`$.

Finally, for the class of densities of states with $`p`$ $`<`$ $`\frac{4}{3}`$ ferromagnetism is absent for $`X=0`$ for all $`U`$, and is only enabled by switching on the bond-charge interaction $`X>0`$. It is useful to consider a model DOS,

$$N(ϵ)=\frac{1+\sqrt{1a^2}}{2\pi }\frac{\sqrt{4ϵ^2}}{2+aϵ},$$
(26)
where the tunable parameter $`a`$ $`=`$ $`1`$ $`\mathrm{}`$ $`1`$ determines the distribution of spectral weight: for $`a=0`$ this DOS reduces to the Bethe DOS, whereas for $`a=1`$ there is a square-root singularity at the lower band-edge similar to the fcc DOS. For $`a=1`$ we have $`p`$ $`=`$ $`1.1353`$, $`r`$ $`=`$ $`0.5006`$, while for $`a=0.9`$ the parameters are $`p`$ $`=`$ $`1.1008`$ and $`r`$ $`=`$ $`0.2821`$. Thus in these cases only a metal-insulator transition is found for $`X=0`$, but for large enough $`X`$ a ferromagnetic phase is predicted, as shown in Fig. 1.

Away from half-filling.— Since the metal-insulator transition takes place only at half-filling, metallic ferromagnetism occupies a larger part of the phase diagram for $`n`$ $``$ $`1`$. Figure 2 shows results for the Bethe lattice with particle density $`n=0.9`$. Whereas for half-filling the strong-coupling phase was insulating, now there is metallic behavior for all couplings, with ferromagnetism setting in for large $`U`$ when $`X`$ is small, and moving to small $`U`$ when $`X`$ becomes large. Compared to Hartree-Fock theory the Gutzwiller approximation predicts a much reduced region of stability of ferromagnetism, which is due to the correlated nature of the Gutzwiller wave function. While in Hartree-Fock theory, owing to the lack of correlations, double occupation can be reduced only through a global spin polarization of the system, the Gutzwiller wave function describes a paramagnetic state with reduced double occupation controlled by the variational parameter $`g`$. Furthermore, the Hartree-Fock prediction of a phase boundary in a range where either $`U`$ or $`X`$ is comparable with the hopping amplitude (see Fig. 2), is not consistent with the weak-coupling nature of this approximation, i. e., self-consistent perturbation theory to first order in $`U`$ and $`X`$. Hence we expect that the Gutzwiller wave function in general provides a quantitatively better estimate than Hartree-Fock theory.

Figure 3 shows the $`U`$ vs $`n`$ phase diagram for the model DOS of Eq. (26) with $`a`$ $`=`$ $`0.9`$ for several values of $`X`$. Already for $`X=0`$ the Gutzwiller theory predicts a large region of ferromagnetic ground states, in qualitative agreement with numerical results. The bond-charge interaction again leads to a further stabilization of ferromagnetism.

Conclusion.— We found that within Gutzwiller’s approach the bond-charge interaction can enhance the instability towards ferromagnetism both at and away from half-filling. This effect is particularly strong when the uncorrelated DOS is asymmetric and there is large spectral weight below the Fermi energy. This provides further support for the conclusion that such a situation is favorable for ferromagnetism. At half-filling, the presence of the bond-charge interaction leads to a metal-insulator transition at lower values than in the standard Brinkman-Rice scenario, since it tends to immobilize the electrons. Although the Gutzwiller theory can be expected to be reliable only at small to intermediate couplings, it represents a major improvement over Hartree-Fock theory, which for example cannot describe a nonmagnetic localization transition.

In conclusion the bond-charge interaction leads to a subtle competition between paramagnetism, ferromagnetism, and localization. Of course, a variational method is not capable of proving the actual stability of a phase. It can only provide estimates for the occurrence of instabilities. Nevertheless, since the Gutzwiller theory treats kinetic and interaction effects nonperturbatively on the same footing, it provides additional insight into the physical mechanism behind these instabilities.

This work was supported in part by the Sonderforschungsbereich 484 of the Deutsche Forschungsgemeinschaft. We would like to thank K. Held for providing us with the Hartree-Fock data of Fig. 2.