text:
# Spectral Properties of Compressible Magnetohydrodynamic Turbulence from Numerical Simulations

## 1 Introduction

The molecular interstellar medium (ISM) is observed to be turbulent at all spatial scales – from the dense, small prestellar cores where the turbulent Mach number $`_s\sigma _v/c_s`$ is order-unity, to the low (mean) density, giant molecular cloud complexes (GMCs) where $`_s`$ is 10 or more (e.g. Blitz (1993)); here $`\sigma _v`$ is the velocity dispersion and $`c_s`$ is the sound speed. Such high Mach numbers imply a strongly compressible state on large spatial scales, which is manifest in the inhomogeneous structure of observed GMCs. Although magnetic field strengths are observationally difficult to obtain (see e.g. Crutcher (1999)), magnetizations at up to kinetic-equipartition levels are possible for GMCs that form in a flux-conservative way from the diffuse ISM. GMCs may thus be thought of as gaseous entities pervaded by compressible magnetohydrodynamic (MHD) turbulence. This turbulence is believed to have many important effects on cloud structure and evolution, but remains only incompletely understood.

In recent years, large-scale numerical MHD simulations have been introduced as a way to explore the fundamental properties of compressible MHD turbulence and hence interpret the dynamics of GMCs (see e.g. Vázquez-Semandeni et al. (2000)). An initial emphasis of these studies was to assess the dissipation rate of compressible MHD turbulence (e.g. Stone et al. (1998), hereafter Paper II; Mac Low et al. (1998); Padoan & Nordlund (1999)), with the general conclusion that magnetic fields do not appreciably lengthen the dissipation time compared to the ratio $`L/\sigma _v`$ associated in an unmagnetized medium with either the crossing time to form shocks (for compressible flow) or the turnover time of an eddy (for incompressible flow). Other studies have focused on characterizing various aspects of the density, velocity, and magnetic field structure, both intrinsic and as observable with spectral or continuum diagnostics (see e.g. Ostriker (2002) and references therein).

An important property of compressible MHD turbulence that has not previously been extensively analyzed is the power spectrum – namely, how the fluctuations of dynamical variables (especially velocity and the magnetic field) vary with spatial scale and angular direction. 111The work of Cho & Lazarian (2002), contemporaneous with our study, also analyzes some spectral aspects of compressible MHD turbulence. These spectral properties reflect how energy is transferred from input scales to smaller (dissipation) scales, and potentially also to larger scales. Spatial variations in heating may depend on the nature of the “forward” energy cascade (to small scales), while the viability of tapping internal sources of mechanical energy (such as stellar outflows) to provide large-scale turbulent support of a cloud depends on the efficiency of the “inverse” energy cascade. Determining what signatures the input scale leaves on the power spectrum may help identify the sources of turbulence. Because spectral properties depend on the mean magnetization of the medium, spectral characteristics inferred from molecular line data cubes or from polarized extinction/emission maps may be useful for diagnosing the magnetic field strength. One of the most well-know characteristics of turbulence in molecular clouds is the general power-law increase of linewidth with spatial scale, as initially described by Larson (1981). Analytical methods are currently being developed (e.g. Heyer & Schloerb (1997); Rosolowsky et al. (1999); Ossenkopf & Mac Low (2002); Lazarian & Pogosyan (2000); Ostriker et al. (2001)) that may enable relatively detailed measurements of the velocity power spectrum from molecular line observations.

With the above motivations in mind, in this paper we analyze a high-resolution (up to $`512^3`$) set of supersonic ($`_s=5`$) driven-turbulence, ideal MHD simulations to characterize their basic spectral properties. Using four different models with Alfvén speeds in the mean magnetic field from zero up to $`2\sigma _v`$, we focus in particular on how the spectral properties depend on the mean magnetization of the system.

The plan of this paper is as follows: We begin in §2 with a brief review of theoretical models and numerical results on power spectra in other regimes of turbulence (incompressible and/or unmagnetized), to place the current work in context. In §3 we describe the numerical methods used to construct our models and the parameter sets we have adopted. Section 4 presents the results of our spectral analyses. First, in §4.1, we examine the angle-averaged (one-dimensional) power spectra of the magnetic and velocity fields, for shear and compressive components separately and in combination. Section 4.2 then analyzes the global, scale-dependent, anisotropy introduced to the power spectral density (PSD) in the presence of a dynamically-strong mean magnetic field. Finally, in §5, we conclude with a discussion of our results and comparison to the results from other numerical simulations, theoretical model predictions, and current observational measures of turbulence.

## 2 Spectral Properties of Turbulence

Most commonly, “turbulence” refers to a complex state of randomly fluctuating velocity fields. Depending on the nature of the medium and the strength of the turbulence, there may also be accompanying fluctuations in the local magnetic field and in the density and temperature. The scale-dependent structure of the turbulence can be described in many ways, with perhaps the simplest characterization in terms of power spectra – i.e. the Fourier transforms of the fluctuating part of the fluid variables. In general, a power spectrum $`P(𝐤)`$ may depend on all three components of the wave vector $`𝐤`$. If the turbulence is isotropic, then $`P`$ depends only on the magnitude $`k(k_x^2+k_y^2+k_z^2)^{1/2}`$ of $`𝐤`$ (also written as $`k_r`$), and a one-dimensional spectrum averaged over polar and azimuthal angles fully characterizes $`P`$. When there is a preferred direction in the medium, as for the case with a strong, large-scale magnetic field, then turbulence may be anisotropic. For a single preferred direction $`\widehat{x}`$, $`P`$ depends only on the components of $`𝐤`$ parallel and perpendicular to $`\widehat{x}`$, and may be characterized by a two-dimensional spectrum $`P(k_{},k_{})`$ averaged over cylinders centered on the $`\widehat{x}`$ axis.

An equivalent way of characterizing the scale dependence of turbulence is in terms of the variance as a function of physical scale; “Larson’s-Law” empirical linewidth-size relations use this representation. For the case of isotropic turbulence with a power-law spectrum, the variance of the corresponding fluid variable will also increase as a power law as a function of the spatial scale over which it is averaged. For example, if $`v^2(𝐤)k^n`$, then the velocity dispersion increases as a function of physical size $`\mathrm{}`$ as $`\sigma _v(\mathrm{})\mathrm{}^{(n3)/2}`$. Larson (1981) originally reported a linewidth-size relation $`\sigma _v(\mathrm{})\mathrm{}^{0.38}`$; larger exponents $`0.50.6`$ have been reported for larger scales with more homogeneous data sets (e.g. Solomon et al. (1987); Brunt & Heyer (2002)). Because a number of issues remain unresolved in extracting the velocity power spectrum from molecular-line data cubes, at present we may regard observed turbulence in GMCs as broadly consistent with a range of possible spectral indices.

The most familiar model of turbulence is that of Kolmogorov (1941), conceived to describe the energy cascade from large to small (dissipation) scales in incompressible, unmagnetized flows. Local, isotropic interactions of eddies with velocities $`v(\mathrm{})`$ at physical scales $`\mathrm{}`$ have durations $`\mathrm{}/v(\mathrm{})`$. An assumption of a conservative energy cascade within an “inertial range” then implies $`v^3(\mathrm{})/\mathrm{}`$ is scale-independent, so that the velocity dispersion on the scale $`\mathrm{}`$ is $`v(\mathrm{})\mathrm{}^{1/3}`$. Assuming isotropy and converting to a power spectrum using $`_kd^3k^{}P(𝐤^{})v^2(\mathrm{})`$ for $`k=2\pi /\mathrm{}`$, the implied angle-averaged Kolmogorov power spectrum thus obeys the scaling $`P(k)k^{2/3}\times k^3k^{11/3}`$. Numerical simulations in the mildly compressive limit (decay models with $`<1`$, and driven models with $`1`$) support the velocity scaling predicted in the Kolmogorov model while showing that the compressive component of energy amounts to 5-10% of the total (Porter et al., 1994, 1998, 1999). Because the Kolmogorov model assumes an incompressible medium, it cannot be expected to apply at large scales in a molecular cloud, where the turbulence is strongly supersonic.

At an opposite extreme from Kolmogorov’s energy-conservative, incompressible turbulence cascade is the Burgers (1974) model of shock-dominated turbulence. In this and related models, the presence of shocks transfers energy immediately from large-scale motions to dissipative scales, rather than via a cascade through “eddies” at intermediate spatial scales. The spectrum $`P(k)k^4`$ associated with Burgers turbulence corresponds to the Fourier transform of a collection of step functions representing many individual shock profiles. Conversion from $`k`$-space to the spatial domain yields $`v(\mathrm{})\mathrm{}^{1/2}`$.

Neither the Kolmogorov nor the Burgers model incorporates magnetic fields, which may be important to the mode of energy transfer through scales – particularly if the total magnetic energy density nears that in the turbulent velocity field. Approaches to theoretical modeling of magnetized turbulence have focused on the incompressible case, in the limit in which fluctuations in the magnetic field are small compared to its mean value; an excellent current review is contained in Cho et al. (2002b). The physical picture in the early theory of Iroshnikov (1963) and Kraichnan (1965) (collectively, “IK”) is based on energy transfer to smaller scales initiated by collisions of counter-propagating Alfvén wave packets. An assumption of an energy-conservative cascade and isotropy (with the implication that many collisions are needed to dissolve a wavepacket) yields a spectrum $`P(k)k^{7/2}`$. However, consideration of the resonant conditions for three-wave interactions among Alfvén waves suggests that the cascade will in fact not be isotropic, but instead should preferentially transfer energy in the direction perpendicular to the mean magnetic field (Shebalin et al., 1983) (hereafter SMM83). Numerical simulations in 2D and 3D by SMM83 and Oughton et al. (1994) directly demonstrated that spectral anisotropy extended in the $`k_{}`$ direction is indeed present; this anisotropy increases toward smaller spatial scales.

A model of incompressible MHD turbulence that uses the colliding-wavepacket picture of IK, while replacing the unphysical assumption of isotropy with the primarily perpendicular cascade of SMM83, was proposed by Goldreich & Sridhar (1995, GS; see also Goldreich & Sridhar (1997)). GS introduce the concept of a “critically-balanced” cascade in which the nonlinear (“eddy turnover”) time of a wavepacket and its propagation time remain comparable, so that $`k_{}v_Ak_{}v(k_{})`$. Combining this relation with an energy-conservative near-perpendicular cascade ($`k_{}v^3(k_{})constant`$) yields a scale-dependent anisotropy law $`k_{}k_{}^{2/3}`$. Numerical simulations in the incompressible limit by Cho & Vishniac (2000), Maron & Goldreich (2001), and Cho et al. (2002a) have verified (via structure functions for coordinates aligned with the local $`𝐁`$-field) that the expected scale-dependent anisotropy is present. The power spectrum from these incompressible MHD simulations is similar to or slightly steeper than the Kolmogorov value (see also Müller & Biskamp (2000)).

To our knowledge, no analytic model of comparable conceptual simplicity to those of Kolmogorov, Burgers, or Goldreich & Sridhar in their respective regimes presently exists for compressible MHD turbulence. In practice, elements of the latter two theories may both apply: at large scales, velocities are strong enough to drive shocks, while at small scales (away from shocks), velocity perturbations are small compared to the sound speed and magnetic perturbations are small compared to the mean magnetic field. However, the overall result cannot be a simple hybrid for the system as a whole, because there is no clean separation of scales when shocks produced by large scale motions transfer energy immediately to the dissipation scale. Results from lower-resolution numerical work indicate relatively steep power spectra (Mac Low & Ossenkopf, 2000), while Cho & Lazarian (2002) have recently reported results showing anisotropy of the expected sense (see also earlier work in the weakly compressible limit by Matthaeus et al. (1996, 1998); Oughton et al. (1998)). The analysis of the present paper focuses on comparing turbulent spectra in strongly compressible models ($`_s=5`$) with varying background magnetic field strength, based on simulations with the highest resolutions ($`512^3`$) performed to date.

## 3 Models and Parameters

Our turbulent cloud models are obtained on a spatial grid of $`256^3`$ zones (models A, B, C, and D) or $`512^3`$ zones (models A2, C2). These models were studied in previous work (Stone et al. (1998), hereafter Paper II), which focused on the steady-state dissipation rate of turbulence. Models A-D here represent the same simulations with dissipation characteristics listed in Table 1 of Paper II, while models A2 and C2 are higher resolution counterparts.

Each model begins as a cubic box of dimension L, initially filled with a stationary, isothermal plasma of uniform density $`\rho _0=\overline{\rho }`$, and threaded by a uniform magnetic field $`𝐁`$ = $`B_0\widehat{x}`$. For each model, the initial magnetic field strength $`B_0`$ is parameterized by $`\beta c_s^2/v_{A,0}^2=c_s^2/(B_0^2/4\pi \overline{\rho })`$, being proportional to the ratio of the thermal to magnetic pressure. We consider models representing a strong field case ($`\beta =0.01`$; models A, A2), a moderate field case ($`\beta =0.1`$; model B), a weak field case ($`\beta =1`$; models C, C2), and a purely hydrodynamic case ($`\beta =\mathrm{}`$, model D).

We evolve each model in time using the ZEUS code (see Stone & Norman (1992a), Stone & Norman (1992b)) to integrate the compressible, ideal MHD equations

$$\frac{\rho }{t}=(\rho 𝐯),$$
(1)
$$\frac{𝐯}{t}+𝐯𝐯=\frac{P}{\rho }\frac{B^2}{8\pi \rho }+\frac{𝐁𝐁}{4\pi \rho },$$
(2)
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁).$$
(3)
We adopt an isothermal equation of state $`P=c_s^2\rho `$, where the sound speed $`c_s`$ is constant in space and time. We view this as an acceptable approximation for modeling most of the gas within a GMC (see discussion in Ostriker et al. (1999)). We do not include any explicit viscous or resistive terms other than shock-capturing artificial viscosity in the MHD equations, and thus dissipation is purely numerical.

Boundary conditions for all of our models are periodic. For a given model, $`\beta `$ remains constant in time due to the periodic boundary conditions, which maintain volume averaged $`\delta 𝐁`$=0. The initial uniformity of $`B_0`$ and $`\rho _0`$ also implies a spatially and temporally constant mass-to-flux ratio, modulo the numerical reconnection effect that may occur if magnetic flux elements of opposite signs are advected into a single cell.

The initial velocity field of each model is $`𝐯_\mathrm{𝟎}=0`$, with turbulent driving implemented by introducing velocity perturbations $`\delta 𝐯`$ at time intervals $`\mathrm{\Delta }t=0.001t_s`$ (where $`t_s=L/c_s`$ is the sound crossing time over the box). The perturbations follow a Gaussian random distribution with Fourier power spectrum

$$|\delta 𝐯(𝐤)|^2k^6\mathrm{exp}\left(8k/k_{peak}\right).$$
(4)
Peak driving occurs at $`k_{peak}L/2\pi =8`$ ($`\lambda _{peak}=L/8`$) for all models except model $`C_{k4}`$, for which $`k_{peak}L/2\pi =4`$ so that it may later be used in a rescaling comparison with model C2. For all of our models, the energy input rate $`\dot{E}`$ due to driving, given by $`\dot{E}/\rho _0L^2c_s^3=1000`$, is constant. The driving is incompressive, with the constraint imposed that $`\delta 𝐯`$=0. The perturbations are normalized such that no net momentum is added. This is accomplished by adding a spatially-uniform component of velocity to compensate for any net momentum associated with the sum of sinusoidal velocity increments, $`\rho \delta 𝐯`$=0 . The resulting saturated-state energies associated with each model were presented in Paper II and will be further analyzed later in §3.

Additional model parameters are provided in Table 1 (see also Paper II), including the turbulent Mach number $`_s\sigma _v/c_s`$ and the Alfvén Mach numbers $`_{A,0}\sigma _v/v_{A,0}`$, $`_{A,rms}\sigma _v/v_{A,rms}`$. Here the Alfvén speed $`v_{A,0}^2B_0^2/4\pi \overline{\rho }`$ associated with the mean magnetic field is distinguished from $`v_{A,rms}^2B^2/4\pi \overline{\rho }`$, calculated via a mass-weighted average of the squared Alfvén speed in each cell, $`|𝐁_0+\delta 𝐁|^2/4\pi \rho `$, over the box. The velocity dispersion is calculated by a mass-weighted average over the box, $`\sigma _v^2|𝐯|^2\rho /\overline{\rho }`$. The sonic Mach numbers are typically $`_s5`$, while values for $`_{A,0}`$ vary from $`0.5`$ for the strong field models up to 5 for the weak field models. The total turbulent energies listed in Table 1 are defined as $`E_K=(1/2)d^3x|𝐯|^2\rho `$ and $`\delta E_B=(1/8\pi )d^3x(|𝐁|^2B_0^2)`$, respectively. Dynamical time scales of interest include the flow crossing time $`t_f(L)=L/\sigma _v`$ and the Alfvén wave crossing time $`t_A(L)=L/v_{A,0}`$ over the box, which are given by $`t_s/_s`$ and $`t_s\beta ^{1/2}`$, respectively. Note that in Paper II, we frequently made use of the flow crossing time at the driving scale, $`t_f(\lambda _{peak})=t_f(L)/8=t_s/8_s`$. The values reported in Table 1 are for the saturated state at late times when the driving and dissipation have reached a balance. 222Values for the $`256^3`$ models are averaged over times $`t=0.25t_s`$ and $`0.30t_s`$. Due to computational expense, the $`512^3`$ models were not evolved for as many time steps. Values for model A2 and C2 are given at times $`t=0.10t_s`$ and $`t=0.15t_s`$ respectively. These times are well beyond the time $`t_f(\lambda _{peak})0.025t_s`$ needed to allow for the development of shocks arising from bulk flow motions.

For the purpose of relating the dimensionless units of our simulations to physical units, we give the sound crossing time as $`t_s=53(L/10\mathrm{p}\mathrm{c})(T/10\mathrm{K})^{1/2}\mathrm{Myr}`$, and the Alfvén wave crossing time as $`t_A=7.6(L/10\mathrm{p}\mathrm{c})(n_{H_2}/10^2\mathrm{cm}^3)^{1/2}(B_0/10\mu \mathrm{G})^1\mathrm{Myr}`$. The value of $`\beta `$ for a given model determines the corresponding physical value of the initial magnetic field strength as $`B_0=1.4\mu \mathrm{G}\beta ^{1/2}(T/10\mathrm{K})^{1/2}(n_{H_2}/10^2\mathrm{cm}^3)^{1/2}`$, while the sound speed is $`c_s=0.19\mathrm{kms}^1(T/10\mathrm{K})^{1/2}`$. The driving power input per volume is then given by $`\dot{E}/L^3=1.04\times 10^{25}\mathrm{ergs}\mathrm{cm}^3\mathrm{s}^1(n_{H_2}/10^2\mathrm{cm}^3)(L/10\mathrm{p}\mathrm{c})^1(c_s/0.2\mathrm{kms}^1)^3`$. Because the current simulations do not incorporate self-gravity, the length scale $`L`$ is arbitrary.

## 4 Results

In this section, we present the results of our spectral analyses. §4.1 focuses on angle-averaged spectral properties; we evaluate inertial-range slopes for various components of the power, and assess how these depend on $`B_0`$. In §4.2, we turn to the issue of characterizing anisotropy in the power spectra. Spectral analysis was performed on the models using a Danielson-Lanczos fast Fourier transform (FFT) and computing the one-sided power spectral density (PSD).

A question of some interest is the relative overall importance and scaling behavior of the compressive versus the shear modes of supersonic MHD turbulence. To compare these, we define separately the shear and the compressive component of the velocity PSD. These are computed in $`k`$-space, after the FFT has been performed, as

$$v_{shear}^2(\stackrel{}{𝐤})=|\widehat{𝐤}\times 𝐯(\stackrel{}{𝐤})|^\mathrm{𝟐}$$
(5)
$$v_{comp}^2(\stackrel{}{𝐤})=|\widehat{𝐤}𝐯(\stackrel{}{𝐤})|^\mathrm{𝟐}$$
(6)
We investigate the behavior of these separate components in both §4.1 and 4.2. In Table 1, we give the volume-integrated energies associated with these two components, for all models.

### 4.1 Angle-averaged power spectra

To study power spectra as a function of length scale only (i.e. angle-averaged), we bin in spherical shells in $`𝐤`$-space so as to generate PSD’s as a function of $`k=k_r(k_x^2+k_y^2+k_z^2)^{1/2}`$. The resulting functions are the shell-averaged (average spectral energy per mode) power spectrum $`P(k)`$ and can be related to the shell-integrated (total in a spherical shell) energy spectrum as $`dE(k)=4\pi k^2dkP(k)/2`$.

For model comparison, we will describe each of the resulting $`P(k)`$ curves qualitatively in terms of four distinct ranges. These are, in order of increasing $`k`$ (decreasing length): the environmental range, driving range, inertial range, and dissipative range. The environmental range extends from the largest spatial scales (the computational box) to the scales at which the forcing is applied. For saturated, driven turbulence, the $`P(k)`$ curve is nearly flat through these largest scales – i.e., there is near energy equipartition among modes. As $`k`$ increases, there is a turnover shortly after $`k>k_{peak}`$, where $`k_{peak}`$ defines the maximum in the forcing function (eq. 4), indicating the start of the inertial range. In the inertial range of $`k`$, the curve $`P(k)`$ is often fit well by a power-law, with slopes varying roughly from -3.3 to -4.9, depending on model parameters. For the present simulations, the inertial range typically spans from $`kL/2\pi `$ = 10 to 30 or 40, but is extended to 60 for the higher resolution models. In some cases, the choice of inertial range is somewhat subjective; the details of such cases are included in the discussion that follows. Finally, as $`k`$ increases still, the curve turns steeper once again at what we identify as the (purely numerical) dissipation range. Through the dissipation range, slopes tend to range between -6.5 and -7. The dissipation range is determined by purely numerical effects, as the code does not explicitly include dissipation terms except for shock-capturing artificial viscosity in the MHD equations. This generic form for $`P(k)`$ applies to both the kinetic and magnetic components of $`P(k)`$, for both strongly and weakly magnetized models, as seen e.g. in Figures 1a,b and Figures 2a,b, respectively.

The primary quantitative measure we evaluate from the $`P(k)`$ in our models is the slope in the inertial range. Table 2 shows the exponent $`n`$ of power-law fits, $`P(k)k^n`$, in the inertial range for the PSD of the (specific) kinetic energy $`P_\mathrm{K}(k)v^2(k)/2`$, the magnetic field energy $`P_B(k)\delta B^2(k)/8\pi \overline{\rho }`$, the combination $`P_{turb}(k)P_\mathrm{K}(k)+P_B(k)`$, the incompressive component of (specific) kinetic energy $`P_{shear}(k)v_{shear}^2(k)/2`$, and the compressive component of (specific) kinetic energy $`P_{comp}(k)v_{comp}^2(k)/2`$. The values presented in Table 2 for the $`256^3`$ models represent the analysis of an averaged data set obtained from two time snap-shots, t=0.25$`t_s`$ and 0.30$`t_s`$ (corresponding to 1.25$`t_f(L)`$ and 1.5$`t_f(L)`$ or 10$`t_f(\lambda _{peak})`$ and 12$`t_f(\lambda _{peak})`$). The values reported for the $`512^3`$ models A2 and C2 are at times $`t=0.10t_s`$ and $`t=0.15t_s`$, respectively.

While many spectra from the models exhibit an easily identifiable power-law range, giving typical errors in slope fits $`\pm 0.1`$, some important exceptions warrant detailed discussion. The PSD curves for density $`\rho `$ (not shown) have no discernible inertial range for any of the models. Curves of $`P_B(k)`$ typically have very little inertial range, and thus there is some ambiguity in choosing the boundaries of the range to fit; typically a range within $`kL/2\pi `$=10 to 30 was used. The values cited were generated by averaging several fits around the center of the most likely power-law, resulting in slope fit errors of $`\pm 0.2`$. In several cases for $`P_\mathrm{K}(k)`$ and $`P_{shear}(k)`$, particularly those of the weak field model, there is a slight “break” in the middle of what we have identified as the inertial range, with shallower slope at larger $`k`$. In these cases, two approaches were taken to fitting the inertial range: (1) an single average over the entire range and (2) fitting the two subranges on either side of the break separately. Results of the first method are cited in Table 2. Details of the second method are given in the discussion of the general results which follows.

One striking effect is that the slopes characterizing the inertial ranges of $`P_{turb}(k)`$, $`P_\mathrm{K}(k)`$, and $`P_{shear}(k)`$ all appear to scale with $`\beta `$, while those of the magnetic field and compressive velocity do not exhibit trends which are as definitive. Values for the PSD slopes of $`P_\mathrm{K}(k)`$ steepen as -4.0, -4.3, -4.7, -4.8 for $`256^3`$ models of increasing $`\beta `$ = 0.01, 0.1, 1.0, and $`\mathrm{}`$, respectively. Steepening of the same sense is also evident for models A2, C2, with $`P_\mathrm{K}(k)`$ slopes -3.8 and -4.3, respectively. The presence of a magnetic field thus appears to reduce intermediate-scale dissipation of energy, in that the inertial range slopes become flatter as one goes from models of no or small magnetic field to those with greater mean field strength.

Spectra $`P_B(k)`$ and $`P_{comp}(k)`$ are notably different from kinetic energy spectra. The slopes of $`P_B(k)`$ vary only between -3.6 and -3.3 as one moves from the strong (A2) to weak (C2) field high-resolution models (with similar behavior at lower resolution). As the mean field weakens, the inertial range of $`P_B`$ becomes progressively harder to identify, being almost entirely unrecognizable for the weak field models. The inertial range of $`P_{comp}(k)`$ is much broader than that of other variables, particularly for the strong field model where the power-law fits extend over nearly twice the range of the corresponding fits for $`P_{shear}(k)`$. However, for moderate and weak-field low resolution models, the extent of the inertial range of $`P_{comp}(k)`$ is reduced again such that an average of fits is needed and the resulting errors in the indices reported are $`\pm 0.2`$.

Notice that there are two main differences between the spectrum of low resolution models and those for their high resolution counterparts. First, the indices in Table 1 for low resolution models tend to be roughly 0.2 steeper than those of the high resolution models. In some cases, this is due to the difference in identification of the inertial ranges. For the low resolution model, the inertial ranges are of lesser extent, which tends to result in fit being placed either closer to the turnover from the driving range or the turnover into the dissipation range; both cases lead to an artificial steepening of the fit. Second, several of the curves for $`P_{shear}(k)`$ and $`P_\mathrm{K}`$(k), for both high and low resolution models, exhibit a slight “break” in the inertial range, as previously remarked. The consequence is that the average fit made through this break results in an apparent discrepancy between the indices reported for high and low resolution counterparts.

The results of fitting either side of the slight spectral break within the inertial range with separate power-laws reveals that in fact there is agreement between the high and low resolution models, with some interesting caveats. The resulting indices of $`P_\mathrm{K}(k)`$ are 4.0,4.3,(4.9:4.5), and 4.8 for models A, B, C, and D, respectively, with the notation $`(n_l:n_r)`$ giving the indices of fits to the left and right of the break (when one exists) with individual errors of $`\pm 0.1`$. Compare these to the indices found for high resolution counterparts A2 (4.0:3.5) and C2 (4.9:3.9), again with individual errors of $`\pm 0.1`$. The fits to the left of the break within the inertial range are in agreement. The case of $`P_{shear}(k)`$ is similar. Average and “left-right” indices $`n`$ and $`(n_l:n_r)`$ of fits to $`P_{shear}(k)`$ for low resolution models are 4.0, (4.5:4.0), (5.0:4.5), 4.9 for models A, B, C, and D, respectively. Again, comparing with the indices for the high resolution counterparts A2 (3.8:3.5) and C2 (5.0:4.0) reveals agreement to the left (smaller $`k`$) of the break within the inertial range. For $`P_B(k)`$, the two-part slopes are (3.5:3.7) for model A2, which may be compared with 3.6 for the “average” slope, and 3.5 for model A.

Note that when a break is apparent within the inertial range, it occurs well beyond the turnover that separates the driving range from the inertial range. Also, the portion of the spectra to the left of the break do not in any way mimic the input driving spectra. Although the current simulations do not have sufficient resolution to be definitive on this point, we speculate that the spectra to the right of the break (larger $`k`$) represent the true asymptotic regime, while the spectra to the left of the break (smaller $`k`$) represent a transition between driving and asymptotic regimes. The steeper slope in the transition range may reflect the envelope of wave distributions in the process of relaxing, via interactions, from the input spectrum to the asymptotic spectrum. This relaxation involves spreading in $`k`$-space accompanied by a decline in the peak amplitude.

We now consider two examples to examine more specific effects which the magnetic field strength has on the form of the PSD curves. Figure 1 illustrates the strong magnetic field case, with the spherically binned PSD for a high-resolution $`512^3`$ model overlayed with that of the lower resolution $`256^3`$ model (models A2 and A respectively, both at time $`t=0.10t_s`$), both of $`\beta `$=0.01.

In Figure 1, there is clear indication of an inertial range between $`kL/2\pi =10`$ and 60 for model A2, showing the cascade of power from larger scale to smaller, dissipative scale at $`kL/2\pi >60`$. The power-law fits to the inertial ranges of $`P_\mathrm{K}(k)`$, $`P_{shear}(k)`$, and $`P_B(k)`$ indicate scalings consistent with that of Kolmogorov turbulence ($`n=11/3`$) or shallower, while that of $`P_{comp}(k)`$ is significantly steeper. Notice from Figure 1 that the lower resolution model (A) has a more limited inertial range than that of the high-resolution model. The doubling of the resolution in effect doubles the extent of the inertial range, because of the greater separation between the driving scales and the dissipative scales associated with the computational grid.

This interpretation is supported by Figure 2, which illustrates the weak field model. Here we again plot $`P(k)`$ curves for a high-resolution model, with $`k_{peak}L/2\pi =8`$. We overlay these with the spectrum from a $`256^3`$ model which was driven with $`k_{peak}L/2\pi =4`$, rather than 8 as before. Here, both models are shown at $`t=0.15t_s`$. For the $`256^3`$ model, we have rescaled $`k`$ by a factor of 2, and rescaled the magnitude by a factor of 10, in order to compensate for the change in $`k_{peak}`$ and the difference in total energies of each model, respectively. Note that the inertial ranges are of the same extent and that the curves practically line up, because the ratio of peak driving scale to dissipation scale are identical.

While driving at lower wave numbers evidently provides a larger inertial range on the $`k>k_{peak}`$ side, our standard choice of $`k_{peak}L/2\pi =8`$ significantly larger than unity allows us to investigate inverse-cascade effects, that is, transfer of energy from the driving scale to the larger “environmental” scales. Evidence of the growth of power on large scales is presented in Figure 3, which shows time series of magnetic and kinetic power spectra from successive model snap-shots (times $`t=0.05t_s`$, $`0.10t_s`$, and $`0.30t_s`$). For both the strong-field and unmagnetized models, power at low $`k`$ increases over time. The unmagnetized and magnetized models have similar, nearly flat, $`P_\mathrm{K}(k)`$ at low $`k`$, whereas $`P_B(k)`$ rises slightly toward $`k=0`$ for the magnetized model. Thus, while indications of an inverse cascade are present, they may be principally due to compressibility.

Note that Figure 3 further indicates that a steady-state saturation has been reached in the spectra for values of $`k>k_{peak}`$ by time $`t=0.05t_s`$. Similarly in Paper II, Figure 1a shows that total energy saturation of forced models was also reached by time $`t=0.05t_s`$. The characteristic timescale for saturation of hypersonic ($`_s>1`$) turbulence is that which characterizes the bulk flow velocity at the driving scale $`\lambda _{peak}`$, namely the flow crossing time $`t_f(\lambda _{peak})=t_f(L)/8=0.025t_s(L)`$ for $`_s=5`$. This is because the relevant shock formation time scale is $`t_f(\lambda _{peak})`$, rather than the sound crossing time $`t_s`$ that would be needed for shock formation in a nonlinear wave steepening process.

How important is the compressive component of the velocity field in the dynamics of turbulence? At moderate Mach numbers, it depends on the ambient magnetic field strength. In Figure 4, we show an overlay of $`P_{shear}(k)`$ and $`P_{comp}(k)`$ for several models of differing field strength, illustrating the difference in magnitude and slope of the shear versus compressible components of the turbulent velocity spectra. Those shown in Figure 4 are the lower resolution $`256^3`$ models, averaged over times $`t=0.25t_s`$ and $`0.30t_s`$.

Table 2 provides for several comparisons between models. For the strong field case (models A2, A), the compressive energy is an order of magnitude smaller than the shear energy, and has a steeper slope (-3.7 vs. -4.3 for the $`512^3`$ models, or -4.0 vs. -4.6 for the $`256^3`$ models). As the magnetic field strength decreases, the difference in relative magnitudes and slopes of the shear and compressive power decreases until, for the pure hydrodynamic model, the difference in magnitude is about a factor of 3 and the slopes, -4.7 and -4.9, are almost indistinguishable, given the errors in average fitting of $`P_{comp}(k)`$ of model D. The magnitude difference in part simply reflects the fact that the forcing is incompressive. For weak field cases, the spectral slopes of the shear and compressive components are also very close (-4.4 and -4.2 for the $`512^3`$ model or -4.7 and -4.5 for the $`256^3`$ model), given the afore mentioned details of fitting averages.

The values we find for the total shear and compressive energies for each model, as listed in Table 1, indicate that $`E_{comp}`$ accounts for 24% of the total kinetic energy for the unmagnetized case (model D), 21% for model C, 19% for model B, and finally only 12% for the strong field case (model A). Comparison of models A, A2; C, C2 shows no dependence on numerical resolution of these ratios.

### 4.2 Directional dependence of the power spectra

We begin with a qualitative look at the turbulent anisotropy as seen in Figures 5, 6. These represent gray scale plots of $`(v/c_s)^2`$ and $`(\delta B)^2/(4\pi \overline{\rho }c_s^2)`$ from a 2D slice taken through the computational box (z=0) for the strong field case ($`512^3`$ model A2) at time $`t=0.10t_s`$. Notice that the structures are elongated along the direction of the mean magnetic field (horizontal axis). There are clearly stronger small scale variations in directions perpendicular to the mean field, suggesting more power in larger $`𝐤_{}`$ and in smaller $`𝐤_{}`$.

To study more quantitatively the global spectral anisotropy induced by the presence of the magnetic field, we have binned the PSD in annular shells, concentric about the axis of the ambient magnetic field direction $`\widehat{x}`$, so as to compare the PSD in the plane perpendicular to the field with that in the direction parallel to the field. The two-dimensional angle-averaged power spectra thus generated are functions of $`k_{}=(k_y^2+k_z^2)^{1/2}`$ and $`k_{}=k_x`$. The 2D contour maps of PSD($`k_{}`$,$`k_{}`$) provide a visual means by which to qualitatively characterize the effect of the magnetic field.

We have computed two-dimensional power spectra for all our models, and find that they become more and more anisotropically distributed with increasing mean magnetic field strength, with more of the power concentrated along the directions perpendicular to that of the ambient magnetic field.

Figure 7 shows contour plots of power spectra for the strong field case ($`\beta =0.01`$, model A2, $`t=0.10t_s`$). In addition to spectra of $`P_\mathrm{K}(k_{},k_{})`$ and $`P_B(k_{},k_{})`$, we also show the shear $`P_{shear}(k_{},k_{})`$ and compressive $`P_{comp}(k_{},k_{})`$ velocity power spectra separately. The solid contours represent the levels of constant power, from log(P($`k_{}`$,$`k_{}`$)) = -3 (smallest $`k`$) to -9 (largest $`k`$) in steps of 1. Three circular dashed curves at k=50, 100, and 150 are overlayed to illustrate a perfectly isotropic relation between $`k_{}`$ and $`k_{}`$ for comparison.

From Fig. 7, notice that the power remains concentrated at larger length scales along the magnetic field (small $`k_{}`$) direction, while extending to smaller length scales in the two perpendicular directions (large $`k_{}`$). There is also a clear indication that the anisotropy is dependent on scale size in the $`P_{comp}`$ spectra (Fig. 7d).

Compared to our results for the $`\beta =0.01`$ model, we find that the power distributions are dramatically more isotropic with a decrease in field strength. Figure 8 shows contour plots for the power in the weak field case ($`\beta =1.0`$, model C2, $`t=0.15t_s`$). While the magnetic field and total velocity field are clearly isotropic, there is a small but real anisotropy along a 45 degree angle between $`k_{}`$ and $`k_{}`$ for both the shear and compressive components of velocity. Note how the compressive field (Fig. 8d) has deficiency in power (extended contour) along the 45 degree line, while the shear field (Fig. 8c) has a surplus (contracted contour) along the same line; the total velocity obtained as the sum of these is, however, isotropic (Fig. 8a).

To further investigate direction-dependent behavior of the power spectra, we have taken slices through our 2D contour plots along $`k_{}=1`$ and $`k_{}=1`$ to generate 1D plots of $`P(k_{})`$ and $`P(k_{})`$, respectively. For all our models, the qualitative morphological features of the directional $`P(k_{})`$ and $`P(k_{})`$ spectra include driving, inertial, and dissipative ranges similar to the angle-averaged PSD curves already shown in Figures 1 and 2.

For the strong-field model, there are clear differences in the directional power that reflect the spectral anisotropy evident in Figure 7; we illustrate the results in Figure 9. From Figure 9, spectral anisotropy, as represented in the difference in magnitude of the overlayed $`P(k_{})`$ and $`P(k_{})`$ curves, appears strongest in the cases of $`P_\mathrm{K}`$ and $`P_{shear}`$, with the scale dependence of this anisotropy clearly evident in the differing slopes characterizing the inertial ranges. The magnitude of the power with variations in the direction parallel to $`𝐁`$ is comparable to that of the power with variations perpendicular to the field only for the largest length scales. As $`k`$ increases, $`P(k_{})`$ falls off much faster than $`P(k_{})`$, such that the anisotropy increases with decreasing length scale.

Within the driving range, the incompressible nature of the driving is clearly evident from the directional spectra in Figure 9, as seen in by the hump near $`k_{peak}`$ which appears strongly in the $`P(k_{})`$ curve for $`v_{shear}`$, but not at all in the curves for $`v_{comp}`$. Within the inertial range, both curves roughly follow a power law. It is difficult to identify a clear separation between the inertial and dissipative ranges for the $`P(k_{})`$ curves. Note that from Figure 9d, for the largest scales there is isotropy in the distribution of compressive velocities, but as $`k`$ increases, $`P_{comp}(k_{})`$ becomes larger than $`P_{comp}(k_{})`$.

To construct quantitative measures of the anisotropy and to see how these vary with changing magnetic field strength, we can compare indices of power-law fits through the inertial range of each of these $`P(k_{})`$ and $`P(k_{})`$ curves, as done before with the angle-averaged power spectra. Table 3 shows the index $`n`$ of the $`k^n`$ power-law fits for $`P_B`$, $`P_\mathrm{K}`$, $`P_{shear}`$, and $`P_{comp}`$ for all models; in Figure 9, these fits are compared to the directional spectra for model A2. Note that such inertial range fits are only able to show an anisotropy in the strong field models. The previously mentioned subtle variations in the weak field models $`P_{comp}`$ and $`P_{shear}`$ (Fig. 8) are not seen here.

There is a general trend that the $`P(k_{})`$ curves for the total and component velocity power spectra exhibit a steeper slope for the weak field than the moderate field case, and more so for the hydrodynamic case. This is especially seen for $`P_K(k_{})`$, where the indices of the power-law fits are 3.5, 4.3, 4.7, and 4.8 for the $`\beta =`$ 0.01, 0.1, 1.0, and $`\mathrm{}`$ models respectively. Note however that the reverse trend is evident for $`P_B(k_{})`$, with spectral indices of 4.5, 4.0, 3.6 for the $`\beta =`$ 0.01, 0.1, and 1.0 models respectively.

The strong field in models A and A2 is responsible for the anisotropy in the power spectra as reflected in the relative slopes along $`k_{}`$ and $`k_{}`$ directions. The steepest profiles for any component of our models are for the strong field model along the direction parallel to the field, while the slope in the perpendicular direction for the same model is shallowest. For $`\beta =0.01`$ (model A2), $`P_B(k_{})k_{}^{4.5}`$ while $`P_B(k_{})k_{}^{3.0}`$. For the same model, $`P_K(k_{})k_{}^{4.7}`$ and $`P_K(k_{})k_{}^{3.4}`$. This is indicative of the fact that increasingly small amounts of power are transfered into the $`k_{}`$ direction at large $`k`$.

As we shall discuss further in §5, the exponents of the power-law fits to $`P(k_{})`$ and $`P(k_{})`$ (for both $`P_\mathrm{K}`$ and $`P_B`$) can be related to a power-law scaling behavior between $`k_{}`$ and $`k_{}`$ that describes the spectral anisotropy (cf. GS). This can be quantified (see Cho & Vishniac (2000)) by asking, for a given power contour, what are the $`k_{}`$ and $`k_{}`$ intercepts? Starting with our 2D spectra, we compute a $`k_{}`$-$`k_{}`$ relation by performing an interpolation of the P($`k_{}`$) curve into the P($`k_{}`$) curve. The resulting function $`k_{}`$($`k_{}`$) provides information on how $`k_{}`$ and $`k_{}`$ are related across a given range of power contours. This scaling is not constant across all magnitudes of power, but is relatively constant within the inertial range. Figure 10 shows plots of $`k_{}`$($`k_{}`$) from the strong field model for $`P_B`$, $`P_\mathrm{K}`$, $`P_{shear}`$, and $`P_{comp}`$. In the inertial range, the relation for $`P_B`$, $`P_\mathrm{K}`$, and $`P_{shear}`$ is consistent with $`k_{}k_{}^{2/3}`$ or $`k_{}^{3/4}`$; the compressive velocity has quite different behavior, roughly $`k_{}k_{}^{0.8}`$. The scaling behavior for $`P_{comp}`$ holds over an extensive range, both at small $`k`$ and for $`k`$ within the expected dissipative range.

We conclude this section with a discussion of one possible application of this work to observations. Namely, we propose a method for constraining the value of the mean magnetic field strength within observed cloud cores by analyzing the degree of anisotropy in centroid maps in molecular tracers. Consider a cloud core which has a nonzero component of the mean magnetic field lying in the plane of the sky ($`xy`$-plane) and a relatively uniform density $`\rho `$ (for the later to be true, the velocities must be sonic or below). Given spatially uniform molecular excitation, the centroid of an emission line at a position $`(x,y)`$ in a projected map represents the mean value of the line-of-sight turbulent velocity, $`v_z(x,y)`$=$`v_z(x,y,z)𝑑z/𝑑z`$. If one takes the Fourier transform of the 2D velocity centroid map, the result $`v_z(k_x,k_y)`$ is identical to taking the $`k_z=0`$ plane from the Fourier transform $`v_z(k_x,k_y,k_z)`$ of the full 3D data cube. In Figure 11a and 11b, we compare the results of these two equivalent calculations, i.e. the FFT of the “centroid map” $`v_z(x,y)`$ and the $`k_z=0`$ slice taken from the FFT of the full 3D data cube of $`v_z`$ for the strong-field model A. The same anisotropy previously discussed, albeit with greater noise, is evident in projected maps. 333We note, however, that anisotropy in velocity centroid maps can be washed out if the emissivity is highly nonuniform, as expected for any region having supersonic flow. Direct tests on our $`=5`$ models show, for example, that density-weighted velocity centroid maps have isotropic power spectra.

In Figure 11c and 11d, we show that $`v_z(k_x,k_y)`$ from the weaker-field models (B,C) has significantly lower anisotropy than does model A. Taken together, these results suggest that anisotropy in the Fourier transform of velocity centroid maps can provide a lower bound on the mean magnetic field strength in an observed cloud. From Table 1, $`_{A,0}=v^2^{1/2}/v_{A,0}`$ in model A is 0.5, while $`_{A,0}`$ in model B is 1.5. Here $`v_{A,0}`$ is the Alfvén speed in the mean field, which for projections along $`\widehat{z}`$ is also the mean plane-of-the-sky field, $`v_{A,p}`$. Allowing for a factor of $`1/\sqrt{3}`$ reduction for $`v_z^2^{1/2}`$ compared to $`v^2^{1/2}`$, this implies that when $`v_z^2^{1/2}/v_{A,p}`$ is $`0.3`$, anisotropy would be clearly evident (unless washed out by other effects such as strong but isotropic density variations), whereas when $`v_z^2^{1/2}/v_{A,p}`$ is $`0.9`$ or larger, negligible anisotropy is expected.

## 5 Summary and discussion

In this paper we present spectral analysis of driven, highly compressive ($`_s5`$) MHD turbulence from a collection of 3D numerical simulations whose parameters span a range believed to apply within GMCs. We investigate models that were the subject of previous work on dissipation of turbulence (Paper II), as well as new high resolution counterparts. Each of our models evolves an isothermal plasma in time subject to an applied incompressive turbulent driving of the velocity field. We parameterize the strength of the mean magnetic field $`B_0`$ by $`\beta =c_s^2/v_{A,0}^2`$ and compare models of varying $`\beta `$. Our main findings, and their relation to other work, are outlined as follows:

1. The angle-averaged spectra for the total velocity and its shear and compressive components ($`P_\mathrm{K}`$, $`P_{shear}`$, $`P_{comp}`$) all exhibit distinctive environmental, driving, inertial, and dissipative ranges. In some of the spectra, there is evidence within the inertial range (i.e. for wavenumbers not directly populated by forcing) of a slight break. In these cases, the portion to the left of the break may represent a transition region, and the portion to the right an asymptotic regime. For strong-field models, the spectrum of $`P_B`$ also shows these distinct ranges. The spectra for $`P_B`$ shows no inertial range for the weaker field models, while spectra of $`\rho `$ show none for any model.

2. Signs of a nonzero but weak inverse cascade toward modal energy equipartition (i.e. flat $`P(k)`$) is evident for both magnetized and unmagnetized models at $`k<k_{peak}`$. This develops within half a flow crossing time over the scale of the box ($`0.1t_s`$, see Fig. 3), corresponding to four crossing times at the driving scale $`\lambda _{peak}`$. Observations of the turbulent spectra in clouds do not show such flattening towards large scales. Our work thus places serious constraints on the importance of turbulent driving at length scales $`\mathrm{}L`$ within GMCs, in comparison to inheritance of turbulence with an extended power spectrum from the process of cloud formation.

3. For the strong field case (model A2), the values of the index $`n`$ for the spectra of $`P_{turb}=(P_\mathrm{K}+P_B)`$, $`P_\mathrm{K}`$, and $`P_{shear}`$ are in the range $`3.73.8`$, while the value for $`P_B`$ is slightly less steep, $`3.5`$. The velocity spectra averaged over the inertial range are thus consistent with the strong MHD turbulence model of GS, which has $`n=11/3`$. In recent numerical work for incompressible MHD, Cho & Vishniac (2000) reported a 1D spectrum consistent with $`n=11/3`$, as did Müller & Biskamp (2000) and Biskamp & Müller (2000). 444Both groups note that numerical schemes employing hyperviscosity show spectral flattening at large $`k`$ from the so-called bottleneck effect; because ZEUS uses finite-differencing rather than a spectral method, our results do not show this bottleneck. Maron & Goldreich (2001) report a shallower law consistent with $`n=7/2`$ for their angle-averaged energy spectra; although they regard this as anomalous due to its inconsistency with GS’s theory, we find it intriguing that our magnetic power spectra are similar. Interestingly, when we fit just the highest-$`k`$ portion of the inertial range, we also find $`n_r=3.5`$ for $`P_K`$ and $`P_{shear}`$ for model A2. Cho & Lazarian (2002) also report a slope similar to $`11/3`$ for Alfvén modes in compressible MHD simulations, although their spectra are somewhat noisy because the inertial range falls at relatively low $`k`$ in their models.

4. For the weak field case (model C2), values $`n`$ of the spectral index for $`P_{turb}`$, $`P_\mathrm{K}`$, $`P_{shear}`$, and $`P_{comp}`$ are in the range $`4.04.4`$, steeper than $`11/3`$ and even exceeding the index $`4`$ associated with Burgers’s spectrum. Although these measurements may be affected in part by spectral steepening toward the driving range, the slopes to the right of the “break” in the inertial range (e.g. 3.9 for $`P_K`$) remain steeper than those in model A2 (e.g. 3.5 for $`P_K`$). Thus, a cloud permeated by a weak mean magnetic field – even if the total magnetic field strength is significant (cf. the total magnetic energy in Table 1) – would be expected to yield a distinctly steeper velocity spectrum than a cloud with a stronger mean field. $`P_B`$ for the weak-field model has a much shallower nominal spectral slope, $`3.3`$, but with large uncertainties due to fitting a limited inertial range.

5. As the mean magnetic field strength decreases, the slopes of the spectra for $`P_{turb}`$, $`P_\mathrm{K}`$, and $`P_{shear}`$ all become steeper, implying that a stronger magnetic field inhibits dissipation somewhat at intermediate length scales. For a weak-field system, velocities are super-Alfvénic (leading to shock dissipation) over a larger range of scales than in a stronger-field system, so that the near-incompressive limit for which the GS model applies would be shifted to larger $`k`$. More generally, compressional effects become progressively stronger for decreasing field strength: the ratio of compressive energy $`E_{comp}/E_\mathrm{K}`$=0.12. 0.19, 0.21, and 0.24 for models of $`\beta =0.01`$, 0.1, 1.0, and $`\mathrm{}`$, respectively.

In part, the lower proportion of compressive energy in the strong-field models must reflect the tendency for energy injected into the Alfvénic cascade (which is purely shear) to remain there (cf. GS, Maron & Goldreich (2001), Cho & Lazarian (2002)). Because losses from compressive components are more rapid than from shear components, however, the steady-state ratio $`E_{comp}/E_{total}`$ may differ from the relative losses through shocks vs. other dissipation (e.g. model A2 loses $`23\%`$ of the input power in shocks, while $`E_{comp}/E_{total}=0.07`$)

6. Two-dimensional power spectra $`P(k_{},k_{})`$ in our strong-field model indicate global anisotropy, with spectral energy preferentially aligned along $`k_{}`$ for $`P_B`$, $`P_\mathrm{K}`$, $`P_{shear}`$, and $`P_{comp}`$. Our weak-field models, however, do not show global anisotropy in their power spectra (but see below). These findings are consistent with the results of Matthaeus et al. (1998) and Oughton et al. (1998) that global spectral anisotropy (as measured using the mean anisotropy angle $`\theta _\omega `$) decreases as the ratio of fluctuating to mean magnetic field increases.

7. Enhancement of anisotropy at large $`k`$ is clearly evident in the two-dimensional power spectra $`P(k_{},k_{})`$, as well as in the cuts $`P(k_{})`$ and $`P(k_{})`$, for our strong-field model. The global spectral anisotropy of $`P_{turb}`$, $`P_B`$, and $`P_{shear}`$ in strong-field models is consistent with the relation $`k_{}k_{}^{2/3}`$ proposed by GS, while for $`P_{comp}`$ we find $`k_{}k_{}^{4/5}`$.

The global spectral anisotropy scaling we identify in our strong-field model has previously been shown to describe local spectral anisotropy for simulated incompressible MHD turbulence by Cho & Vishniac (2000) and Maron & Goldreich (2001). As those works emphasized, the anisotropy predicted by GS is for $`\widehat{k}`$ directions defined with respect to the local magnetic field. If the large-scale field is strong in the sense that the turbulence is sub-Alfvénic (as for our models A, A2), then the local directions of field lines do not depart strongly from $`\widehat{B}_0`$; for weak large-scale fields, on the other hand, local directions of $`\widehat{B}`$ are uncorrelated with $`\widehat{B}_0`$ (see Fig. 2 of Paper II). Thus, our finding of global spectral anisotropy only in our strong-field models is consistent with expectations, since only for models A, A2 does $`\widehat{B}_0`$ correlate well with the local field.

Discerning small-scale local anisotropy even when turbulence is super-Alfvénic on large scales is possible using second-order directional structure functions for coordinates defined with respect to the local magnetic field (Cho & Vishniac, 2000; Maron & Goldreich, 2001). Although we have not applied these specialized techniques here, our results for weak-field models are also consistent with the findings of Milano et al. (2001) that global spectral measures may be isotropic even when local anisotropy is present as identified by other methods.

8. The dependence of global velocity spectral anisotropy on the mean magnetic field strength suggests a potentially important observational diagnostic. Namely, observational evidence of turbulent anisotropy at a given spatial scale would imply that the value of the “mean field” Alfvén Mach number – i.e. based on the plane-of-the-sky $`|𝐁|`$ for that spatial scale – is less than unity. In molecular cloud core regions where the density varies weakly in space, Fourier transforms of velocity centroid maps correspond directly to slices through the 3D velocity power spectrum, so evidence of anisotropy is immediately accessible (see Fig. 11). In regions where the flow is supersonic, convolution of the anisotropic velocity power spectrum with the isotropic density power spectrum tends to wash out anisotropy in velocity centroid maps, but more complex analysis techniques using the full position-velocity data cubes may still be able to ascertain the degree of anisotropy in the velocity power spectrum itself. Further investigation will show if such spectral diagnostics offer a practical indirect means of constraining the magnetic flux through a cloud core, and hence in determining whether it is magnetically sub- or super-critical with respect to gravitational collapse.

9. Finally, we briefly comment on the relationship of our results to proposals presented in Boldyrev (2002) and Boldyrev et al. (2002) that the spectral slope for compressible MHD turbulence may be derived from a modified version of a She-Lévêque analysis that provides intermittency corrections to Kolmogorov scalings (see She & Lévêque (1994); Dubrulle (1994)). Those works suggest that for a range of fractal dimensions for the dissipative structures, the index of the power spectrum would lie in the range $`3.743.89`$. For our fits over the whole inertial range, we find indices for $`P_K`$ of 3.8 and 4.3, respectively, for models A2 and C2, and indices 3.7 and 4.4 for just the shear component of velocity $`P_{shear}`$. Fitting just the high-$`k`$ portion of the inertial range for model A2 (C2), we find $`n_r`$ of 3.5 (3.9) for $`P_K`$ and 3.5 (4.0) for $`P_{shear}`$. Thus, while our results for the strong-field model A2 are potentially consistent with the Boldyrev et al proposal, the spectra of our weak-field model C2 are steeper. While detailed discussion of the reason for this difference is beyond the scope of this paper, we note that the basic assumption of a conservative cascade (implicit in the assumption that the slope of the third-order structure function $`\zeta (3)=1`$) is more applicable in the strong-field model that the weak-field model. It remains an open challenge to develop an analytic model that incorporates both direct shock dissipation and multi-scale energy cascades in a single framework.

We are grateful to Charles Gammie for helping to initiate this project and to the referee for a valuable and careful report. This work was supported in part by NASA grants NAG53840 and NAG59167.


format: text

text:
# Self-consistent approach to deformation of intruder states in neutron-deficient Pb and Po

## Abstract

We present systematic calculations of the properties of $`11^{}`$ isomers in neutron-deficient 184-198Pb and 188-200Po. These states are based on the $`\pi (h_{9/2}i_{13/2})_{K^\pi =11^{}}`$ configuration. They are calculated in the framework of the Hartree-Fock-Bogoliubov method with a Skyrme interaction and density-dependent pairing force. The energies and deformations of the $`11^{}`$ states are compared to those of the intruder $`0^+`$ states in both Pb and Po isotopes. In the most neutron-deficient Po isotopes, the calculations predict, below a weakly oblate $`11^{}`$ state, another oblate $`11^{}`$ state which is even more deformed than the intruder oblate $`11^{}`$ state in their Pb isotones. The energies and quadrupole moments of the $`11^{}`$ isomers, corresponding to a weakly oblate nuclear shape, are in fair agreement with the available experimental data.

The neutron-deficient nuclei around the $`Z=82`$ shell closure constitute one of the most famous examples of shape coexistence Heyde83 ; Wood92 . In Pb isotopes, oblate, prolate and superdeformed configurations coexisting with a spherical ground state were predicted theoretically (see AFN90  and references therein) and identified experimentally Wood92 ; JK91  more than ten years ago. By now, low-lying $`0^+`$ states, which are a characteristic feature of shape coexistence, are known from 202Pb down to 184Pb VDCo ; BijnAnd ; Andrei00 ; AndreiEPJ99 . The interpretation of these low-lying $`0^+`$ states relies on proton excitations across the $`Z=82`$ spherical shell gap. In a shell model picture HeJo87 ; HeVI94 , the energy of these multi-particle-multi-hole (mp-mh) configurations is lowered by the quadrupole-quadrupole and pairing interactions. In mean-field models BeNa89 ; Naz93 , these structures are described by a state with an oblate (prolate) quadrupole deformation, whose low energy is due to the occupation of deformed intruder orbits, such as the down-sloping $`\pi h_{9/2}`$-orbital for an oblate state.

High-spin isomers constructed on these $`0^+`$ states are of particular interest since the measurement of their static electromagnetic moments provides direct information on both their single-particle structure and their deformation. Several isomeric states are known in Pb isotopes. In particular, $`11^{}`$ levels have been observed from 188Pb to 198Pb (see, e.g. Ref. JuHe01 ). They are not known in heavier Pb, most probably because they are not yrast. These states are interpreted as a two-quasiparticle (2qp) excitation based on the $`\pi (h_{9/2}i_{13/2})_{K^\pi =11^{}}`$ single-particle configuration as confirmed in 196Pb by the measurement of its $`g`$-factor PeHe87 . Experimental information on deformation is extracted from the collective bands based on the excited $`0^+`$ VDCo ; PeHe87 ; Andrei00 ; DrBy02  and $`11^{}`$ states ClWa93 ; DrBy98 . Recent data on static quadrupole moments for $`11^{}`$ states in 194,196Pb VyCh02 ; VyOr02  support the mean-field interpretation of these states as weakly deformed oblate configurations.

The experimental excitation energies of the $`11^{}`$ states and of the $`0_2^+`$ states in Pb isotopes are shown in the left part of Fig. 1. Both exhibit a similar behavior. The steady decrease of energy with decreasing neutron number indicates that, in a shell model picture, the residual quadrupole and pairing interactions, or, in mean-field models, the deformation effects, are maximal at the middle of the neutron sub-shell.

Similar coexisting structures based on intruder configurations are known experimentally above and below the $`Z=82`$ shell closure, in the neutron-deficient Hg, Pt and Po isotopes Wood92 ; JuHe01 . The energy systematics of the low-energy excited $`0^+`$ states in neutron-deficient 200-196Po BijnPRL ; JuHe01  strongly resembles that in Pb. The observation of bands based on these states seems to indicate a strong mixing between different shapes. The $`11^{}`$ isomers have been identified from the semi-magic 210Po down to 194Po JuHe01 . In contrast to the Pb case, their excitation energy is nearly constant as a function of the neutron number.

Previous theoretical studies of these isotopes have been mainly devoted to the description of the $`0^+`$ intruder states, in particular in Pb. A large variety of models has been used: a schematic spherical shell model VDCo ; HeJo87 ; DCDe00  and its symmetry-based boson truncations FoHe02 , the Strutinsky approach with a Woods-Saxon potential MaPa77 ; BeNa89 ; Naz93 , and several variants of self-consistent mean-field methods TaFl93 ; ChEg01 ; HeVa01 ; DuBe02 . A comprehensive study of the low-energy spectroscopy and ground-state properties of light Po isotopes has been presented within several schematic approaches in Ref. OrHe99  (see review JuHe01  for more references). On the other hand, the only systematic study of isomers in these neutron-deficient nuclei is based on cranked Nilsson-Strutinsky type calculations with a schematic Woods-Saxon potential and is devoted to 186-200Pb BeNa89 ; ClWa93 . The quadrupole moments of $`11^{}`$ intruder states in Pb and Po down to 194Pb and 198Po have also been analyzed recently within the particle-core coupling model VyOr02 . The deformation properties of the $`11^{}`$ state in 196Pb have been discussed in the framework of the tilted-axis cranking model VyCh02 . The magnitude of the $`11^{}`$ nuclear deformation is indeed a basic ingredient in this model for the description of magnetic rotational bands in both even-$`A`$ and odd-$`A`$ neutron-deficient Pb isotopes ClMa00 ; VyCh02 . However, up to now, no fully microscopic investigation of these isomeric configurations has been performed.

In this letter, we present the results of the first fully self-consistent description of the properties of $`11^{}`$ intruder isomers in 184-198Pb and in 188-200Po based on the Hartree-Fock-Bogoliubov (HFB) method. Our aim is to analyze the nature of these states and to study the evolution of their deformation as a function of neutron number. The simultaneous consideration of Pb and Po isotopes will allow also to determine the effect of the $`Z=82`$ spherical shell closure on the deformation trends. Since several Po isotopes have extremely shallow energy surfaces as a function of the axial quadrupole moment, we have also performed selected configuration-mixing calculations to validate the construction of isomers on an oblate configuration.

All our calculations have been performed with the same set of nucleon-nucleon effective interactions. In the mean-field channel, we use the Skyrme parameterization Sly4 Sly4  and in the pairing channel, a density-dependent zero-range force mainly active on the nuclear surface as introduced in Ref. RiBo99 . The HFB equations have been solved on a three-dimensional cubic mesh TeHe96  for reflection symmetric nuclear shapes with the method presented in Refs. GaBo94 ; TeHe96  and include an approximate particle-number projection with the Lipkin-Nogami prescription LN . The isomeric states are constructed fully self-consistently as $`2qp`$ excitations as in Ref. HeJa98 , selecting the quasiparticles whose dominant components in the HF basis are the $`9/2^{}`$ and $`13/2^+`$ particle states.

In Pb isotopes, we have taken for the energy and the deformation of the first excited $`0^+`$ states those of the oblate minimum of the deformation energy curves. Configuration-mixing calculations TaFl93 ; ChEg01 ; HeVa01 ; DuBe02  have shown that this approximation is well justified, even when the potential well of the oblate minimum is very shallow. The topology of the energy curves for Po isotopes is more complex, with several shallow minima on both prolate and oblate sides. The properties of the first excited $`0^+`$ state have therefore been determined by performing configuration-mixing calculations with respect to the axial quadrupole collective coordinate. To avoid the ambiguity resulting from the fact that the states obtained from the configuration-mixing calculations do not have the same mean-particle number, we have included in the configuration mixing a projection on the correct proton and neutron numbers.

This set of methods has already been applied with success to several other problems, such as the description of neutron-rich nuclei TeHe96 ; TeFl96 , of superdeformed rotational bands TeFl97 ; HeJa98  and of heavy nuclei DuBo01 .

Deformation energy curves
The deformation energy curves obtained for 194-184Pb and 198-188Po isotopes are shown in Figs. 2 and 3, respectively. These curves have been obtained by HF+BCS calculations with a constraint on the axial mass quadrupole moment $`qq_{20}=2z^2x^2y^2`$. Our formalism for configuration mixing calculations is indeed up to now limited to the mixing of HF+BCS wave functions and does not permit to use HFB ones. For this reason, we shall base our discussion of the deformation properties on the HF+BCS energy curves. We have checked that very similar curves are obtained within the constrained HFB method and that the trend in the $`0_2^+`$ excitation energy as a function of the neutron number is the same. In order to make a connection with simpler models, we have labelled orbitals by approximate Nilsson indices.

For all Pb isotopes (Fig. 2), the ground state corresponds to the spherical minimum. With decreasing neutron number, two structures appear, first as shoulders, then as secondary minima on the oblate (starting from $`A=192`$) and prolate sides (starting from $`A=188`$). The oblate minimum occurs around $`q=10`$ b which is the deformation corresponding to the crossing between the down-sloping $`9/2[514]`$ mean-field orbital and the up-going orbital $`1/2[440]`$. On the prolate side, the minimum at $`q20`$ b corresponds to the deformed shell gap due to the $`1/2[530]`$, $`1/2[541]`$ and $`3/2[532]`$ intruder orbitals and the rising of the $`1/2[400]`$, $`3/2[402]`$ and $`11/2[505]`$ orbitals. For 186Pb, the excitation energy of the prolate minimum significantly overestimates the experimental data. This discrepancy is irrelevant for the description of the isomeric states which are constructed on the oblate minimum. The energy of this prolate state is very sensitive to the details of the mean-field method. Several other mean-field calculations, based either on the HF+BCS method and other Skyrme interactions TaFl93 ; DuBe02 , or on the HFB method with the Gogny force ChEg01 , predict a prolate minimum lower in energy than the oblate one.

As we have explained in the introduction, for the Pb isotopes, we have approximated the energy of the oblate $`0^+`$ state by the energy of the oblate minima. These energies are compared with the experimental data in the left part of Fig. 1. The decrease of the HF+BCS energies with decreasing neutron number is in good agreement with the data, although the experimental energies are overestimated. This overestimation is even larger by approximately 300 keV in HFB calculations.

The deformation energy curves for Po isotopes (Fig. 3) exhibit a complicated coexistence of several shallow minima. To check whether these minima are stable against triaxial deformations, we have performed selected calculations including constraints on the triaxial quadrupole moment. These calculations show that the symmetric minima corresponding to $`q\pm 8`$ b ($`\beta \pm 0.1`$) obtained for all Po isotopes are separated by a very low barrier of about 100 keV in the $`\gamma `$ plane. They are therefore probably unstable against $`\gamma `$ vibrations. The shoulder that appears in 198Po for a deformation around $`15`$ b becomes a deep well for lighter isotopes which is the ground state for 192,190Po. In 188Po, the oblate minimum is nearly degenerate in energy with a prolate minimum at $`q19`$ b. The Strutinsky calculations of Refs. MaPa77 ; BeNa89 ; OrHe99  agree qualitatively with our results. It is clear from this discussion that a study of the low lying levels of Po isotopes must take into account the vibrational motion with respect to the axial quadrupole mode.

Properties of low-lying oblate $`0^+`$ states: configuration-mixing calculations for Po isotopes
We have performed configuration-mixing calculations with respect to the axial quadrupole moment for the Po isotopes within the generator coordinate method (GCM), as was presented in Refs. BoDo90 ; HeBo93 . The basis for the calculations is a set of wave functions $`\mathrm{\Phi }(q)`$ obtained by HF+BCS calculations with a constraint on the axial quadrupole moment. These wave functions $`\mathrm{\Phi }(q)`$ are projected on the correct number of protons and neutrons. Collective wave functions are then constructed as a linear superposition of the projected non-orthogonal basis functions. The weights corresponding to the different quadrupole moments are determined by diagonalizing a discretized version of the Hill and Wheeler equation HiWh53 ; BHR03 . This procedure provides as many eigenstates as discretization points of the quadrupole moment. The lowest eigenvalue corresponds to the ground state with a wave function spread over the axial quadrupole moment. As we shall see, depending on the structure of this wave function, it represents either a vibration within a well-defined well or a mixing of mean-field states with different shapes. In the same way, the other eigenvalues of the Hill and Wheeler equation represent excited states. Depending upon their shape, they can be classified into three main categories. Some states are mainly located in a well that is clearly separated from the ground state (as it is the case for a super-deformed well). Other states can be $`\beta `$ vibrations with respect to a lower state in the same potential well. Finally, some states can represent a mixing of states corresponding to different shapes. We have not included triaxial deformations in these configuration-mixing calculations. Some Po isotopes being soft against the $`\gamma `$ degree of freedom, these deformations may affect some of our results. However, their inclusion would considerably increase the computation time, and since our main purpose is to calculate the $`11^{}`$ isomers, we have limited the configuration mixing to axial shapes.

The energies of the first few states obtained by the diagonalization of the Hill and Wheeler equation are shown for selected Po isotopes in the upper part of Fig. 4, together with the deformation energy curves. The triangles representing the energy of the configuration-mixed states are positioned at their average quadrupole moment. The amplitudes of the collective eigenfunctions (which are related to the mixing coefficients by an integral transformation RS80 ) are plotted in the lower part of the figure for the three lowest states obtained for each isotope. The wave functions of these states are spread over a large range of quadrupole moments. In 200,198,196Po, the ground state and the first excited state are a mixing of an oblate and a prolate configuration, resulting in an average small oblate deformation. This mean quadrupole moment is close to the one of the first shallow oblate minimum of the deformation energy curve. In the three isotopes, the second excited $`0^+`$ state has an average quadrupole moment $`q10`$ b, corresponding to a weakly oblate shape. While it has only a slightly dominant weight in the ground-state wave function of 200,198,196Po, this oblate configuration becomes the ground state in 194,192,190Po, with still a spreading over a large range of quadrupole moments. In 190Po, it is nearly degenerate with a predominantly prolate configuration. These results agree well with the interpretation of the experimental data that the ground states of 200-196Po isotopes are mixed but nearly spherical, the ground state of 194Po is mixed, with a large oblate component and that of 192Po is predominantly oblate BijnAnd ; HeCo99 . Experiments also suggest a strong mixing between oblate and prolate shapes in the 190Po ground state (JuHe01 ; Po190 ).

In the right part of Fig. 1, we compare the three first eigenstates of the configuration-mixing calculation with the available experimental data. The states that have the largest oblate mean deformation are indicated by filled triangles. The experimental excited $`0^+`$ states are known in 200,198,196Po BijnPRL  and are interpreted as 4p-2h excitations in the spherical shell model (see Refs. BijnPRL ; HeCo99 ; JuHe01  and references therein) or as states with an oblate deformation in mean-field models. The experimental trend is well reproduced by the most oblate states, although it is not the first excited state. Note however that the experimental energies are intermediate between those of the two first excited states. This may indicate that the interaction that we use leads to a too large mixing between the deformed configurations. Such a result could be modified by the inclusion of triaxial deformations, but in a way which can hardly be anticipated.

From the deformation energy curves for the Po isotopes, the $`11^{}`$ states can be constructed on either of the two oblate minima. The configuration-mixing analysis clearly favors the slightly oblate minimum for the heavier isotopes. For the lighter isotopes (194-188Po), the mean deformation of the configuration-mixed ground state is intermediate between the two. For these nuclei, we shall therefore explore whether an $`11^{}`$ isomer can be constructed for both deformation regions.

Properties of the $`11^{}`$ isomers
The $`11^{}`$ states in the Pb and Po isotopes have been constructed as a $`2qp`$ excitation on an HFB vacuum. As in Ref. HeJa98 , the $`2qp`$ state is calculated fully self-consistently, which means that effects such as core polarisation due to the time-odd terms of the interaction and the conservation of the number of particles or changes of deformation are taken into account. Since the single-particle basis, which is formed by the eigenstates of the Hartree-Fock Hamiltonian, does not diagonalize the $`qp`$ basis, the $`qp`$ excitations cannot be defined exactly by single-particle states. We have always chosen the $`2qp`$ wave functions which have the largest overlaps with the lowest $`9/2^{}`$ and $`13/2^+`$ single-particle states. In practice, these overlaps are always larger than 0.95. As in previous calculations of $`1qp`$ and $`2qp`$ excited states HeJa98 ; DuBe02 , the isomeric state has been obtained axial in all cases, although triaxial deformations are allowed by the symmetries of our computer code.

The experimental and theoretical excitation energies of the $`11^{}`$ states are compared in Fig. 1 for the Pb and Po isotopes. In the Pb isotopes, the energy of the $`11^{}`$ states decreases with decreasing neutron number, both in experiment and theory. The experimental trends are quite well reproduced by our calculation, although the excitation energy is overestimated by 600 to 800 keV. This overestimation is most probably related to a better description in the HFB formalism of the vacuum (corresponding to the ground state) than of $`qp`$ excitations. Our calculation shows that the excitation energy of the $`11^{}`$ state is the lowest in 190,188Pb and increases in 186,184Pb for which there are no experimental data yet.

For the Po isotopes, the $`11^{}`$ states are known down to 194Po. The isomer has been identified in 192Po HeCo99 , but its energy could not be obtained from the data. In 200-196Po, a minimum corresponding to an $`11^{}`$ state could be constructed only for small oblate deformations. The theoretical excitation energies are in good agreement with the data. In the case of 194-188Po, $`2qp`$ states have been obtained for both oblate deformations. The isomers corresponding to smaller (larger) absolute values of the quadrupole moment are shown by filled (empty) diamonds in Fig. 1 (right part). The excitation energies of both solutions are close to each other. These states have to be further mixed by a configuration-mixing calculation. One should not expect a significant change in the excitation energy of the lowest mixed eigenstate, the ground state being lowered in a similar way by configuration mixing.

The energies of the $`11^{}`$ states are nearly constant as a function of the neutron number and the isomers occur at a lower energy than in Pb. The $`11^{}`$ energies do not follow the steep decreasing trend as the intruder $`0^+`$ energies. This suggests that the structure of the $`11^{}`$ isomers is rather similar for all $`Z=84`$ isotopes, at least down to 196Po. In lighter Po isotopes, the $`11^{}`$ states built on the most oblate minimum become lower in energy than the weakly oblate ones, and configuration-mixing calculations would be necessary to get more accurate predictions on the $`11^{}`$ energy and deformation.

In the terminology of phenomenological models, the constant energy of the $`11^{}`$ states in heavier Po isotopes, indicates that the self-consistency effects due to the $`2qp`$ excitation do not affect strongly the core. This is not the case for the Pb isotopes, in which the $`Z=82`$ magic number is broken by the $`2qp`$ excitation, inducing an additional core deformation, as also demonstrated by the trend in the intrinsic quadrupole moments.

The calculated intrinsic charge quadrupole moments of the $`11^{}`$ states in Pb and Po are plotted in Fig. 5. For comparison, there are also given two available data points for 194,196Pb extracted from the measured spectroscopic quadrupole moments in the assumption of the axial symmetry and $`K=11`$ VyCh02 . The theoretical results are in good agreement with the experimental values. The deformation parameters ($`\beta _2`$, $`\beta _4`$) defined as in Ref. CiDo96  vary for Pb isotopes from $`(0.15,0.014)`$ for 198Pb to $`(0.19,0.018)`$ for 184Pb. The quadrupole moments of the intruder isomers in the Pb-isotopes are significantly larger than those of the $`11^{}`$ states corresponding to the weakly oblate shape in the heaviest Po isotopes. This reflects the influence of core polarisation due to the breaking of the $`Z=82`$ core in the Pb isotopes.

In lighter Po isotopes, two candidates for an $`11^{}`$ isomer have to be mixed by a configuration-mixing calculation. It is difficult to predict the quadrupole moment of the mixed state which could have an intermediate value between those of the two mean-field states. From these results, one can infer a smooth change from a slightly deformed oblate shape of $`11^{}`$ states in the heaviest Po to a more deformed oblate shape in the lightest. Experimental observation of the rotational bands above these isomers or a direct measurement of the spectroscopic quadrupole moments would be helpful to clarify the situation. For the unambiguous $`11^{}`$ states in the heavier Po isotopes, the deformation parameters ($`\beta _2`$, $`\beta _4`$) vary from $`(0.114,0.007)`$ for 200Po to $`(0.117,0.007)`$ for 196Po.

In summary, we have, for the first time, applied the HFB method with a Skyrme interaction and a realistic pairing force to the study of the $`11^{}`$ isomers in a series of neutron-deficient Pb and Po nuclei. The isomers have been constructed self-consistently as $`2qp`$ excitations $`\pi (h_{9/2^{}}i_{13/2^+})_{11^{}}`$. Their main properties are in good agreement with the available experimental data. We obtain that the energy of $`11^{}`$ states in Pb isotopes reaches minimum around $`N=106108`$ and increases in lighter isotopes for which there are no data yet. For Po isotopes, the $`11^{}`$ energy is predicted to remain roughly constant as a function of neutron number. The calculations support the assumption of oblate deformation of these nuclei in their $`11^{}`$ states, with a steady increase towards neutron-deficient isotopes. The deformation is more pronounced in Pb than in Po, but in lighter Po isotopes the deformation of the $`11^{}`$ states might become of the same order or even larger than that in the Pb isotopes. We have also studied the properties of the low-lying $`0^+`$ states in Po isotopes within the configuration-mixing approach on the basis of HF+BCS wave functions. These calculations show that, in all cases, the lowest $`0^+`$ states do not have a well defined deformation, and oblate, spherical and prolate shapes are strongly mixed, both in the ground and the excited $`0^+`$ states.

###### Acknowledgements.

We acknowledge discussions with A. Andreyev, D. Balabanski, M. Huyse, K. Van de Vel and P. Van Duppen on experimental data and thank M. Bender for his help in the GCM calculations. N.A.S. thanks Kris Heyde for a careful reading of the manuscript and valuable comments. This work is supported by the DWTC grant IUAP $`\mathrm{\#}P5/07`$.