text:
# Hall effect in the marginal Fermi liquid regime of high-𝑇_𝑐 superconductors

## I Introduction

Soon after the discovery of high-temperature superconductivity in the cuprate compounds, it was found that all the transport properties for compositions near those for the maximum $`T_c`$ have anomalous temperature and/or frequency dependence, in contrast to what is expected for Landau Fermi liquids. Most of the anomalies, in resistivity and optical conductivity, in the frequency dependence of Raman intensity, in the tunneling spectra, and in the temperature dependence of the copper nuclear relaxation could be understood by the “Marginal Fermi Liquid” (MFL) phenomenology.mfl  Within MFL, scale-invariant fluctuations in both in the charge and magnetic sectors are assumed to have the form,

$`\mathrm{Im}\chi (𝐪,\omega ,T)`$ $`=`$ $`N_0(\omega /T),\omega T`$ (1)
$`=`$ $`N_0(\mathrm{sgn}\omega ),T\omega \omega _c,`$
characteristic of a quantum critical point. $`\omega `$ is the frequency and $`𝐪`$ the momentum of the fluctuations, $`N_0`$ is the density of energy states per unit volume and $`\omega _c`$ is a high-frequency cutoff. At long wavelengths this form is observed directly in Raman scatteringSlakey . A principal predictionmfl  from this hypothesis is that at low energies ($`\omega T`$), the inelastic part of the single-particle relaxation rate has the MFL form $`\lambda T`$ with coefficient $`\lambda `$ having negligible momentum dependence either along or perpendicular to the Fermi surface. This form has been confirmed in angle-resolved photoemission (ARPES) experimentsarpes  and leads directly to the observed linear temperature dependence of the resistivity.

However, the temperature dependence of anomalies in the normal state Hall effect and magnetoresistance could not be understood by the MFL hypothesis alone. For example, for a situation near optimal doping, where the resistivity is linear in temperature, the expectation is that the cotangent of the Hall angle ($`\sigma ^{xx}/\sigma ^{xy}`$) should also be linear. Experiment shows it to be more nearly quadratic. This left open the possibility that essential new physics near optimum doping may not be captured by the MFL scenario. In this paper we show that the temperature dependence of the Hall angle may be understood quantitatively by a proper application of transport theory using the measured single-particle relaxation rates.

The single-particle self energy $`\mathrm{\Sigma }(𝐤,\omega ,T)`$ is measurable in angle-resolved photoemission (ARPES) experiments.arpes ; pnas  As stated above, within MFL, the prediction is that the inelastic part of $`\mathrm{\Sigma }`$ is $`𝐤`$-independent and of the MFL form,mfl  proportional to $`\omega `$ for $`\omega T`$ and to $`T`$ for $`T\omega `$. ARPES experiments arpes  do find the inelastic part $`\mathrm{\Sigma }`$ of this form but find in addition an elastic part which varies in magnitude around the Fermi surface. Thus, as explained in detail earlier,pnas  on the Fermi surface ($`\omega =0`$), the experimentally-measured self energy consists of the MFL part $`\lambda T`$ and an anisotropic $`T`$-independent elastic part:

$$\mathrm{Im}\mathrm{\Sigma }(\stackrel{}{k},T)=\lambda T+\gamma (\widehat{k}),$$
(2)
According to experiment, $`\gamma (\widehat{k})`$ increases by about a factor 4 to 5 going from the $`(\pi ,\pi )`$ to the $`(\pi ,0)`$ direction along the Fermi surface.pnas  A crucial point is that even at its minimum value, it is more than an order of magnitude larger than the transport scattering rate due to impurities obtained by extrapolating the normal state resistivity to $`T=0`$. As argued in Ref. 4, such behavior arises if $`\gamma (\widehat{k})`$ comes from small-angle impurity scattering. We discuss this point further in Sec. VII.

In a previous communication,prl  we described, for high-$`T_c`$ superconductors, how Eq. (1) can be used in a Boltzmann equation analysis to account for observed anomaliesong ; exp ; hwang ; ando ; raffy ; drew  in the Hall effect. We performed a calculation (page 4655 of Ref. 1) using a simple Fermi surface in order to give an example of how a new contribution corr  could dominate the conventional result and thus account for experimental observations. However it was pointed out to us by V. Yakovenko vy  that we erred in the form we chose to parameterize $`\gamma (\widehat{k})`$.err  The purposes of this paper are to give a more complete solution of the Boltzmann equation, correct the above error, and to present a model calculation which illustrates how, with a proper parameterization and a reasonable choice of parameters, the presence of the anisotropic impurity scattering together with other anisotropies can account for the observed behavior of the Hall angle.ong ; exp ; hwang ; ando ; raffy ; drew  Further, we generalize our results to the complex Hall angle at finite frequencies which has been recently measured.matt

A point about the experimental results needs stressing. In many presentations of the data, for example Refs. 9,10, single temperature non-integer power-law fits to the Hall number or the Hall angle have been made, with varying success, while integer power law fits ong  do not usually faithfully represent the data. There is no physical reason to expect integer power law behavior and our analysis shows quite generally that the solution of the Boltzmann equation can give a sum of contributions with different temperature dependences. In particular, it is usually foundando ; raffy  that a temperature power law with a power less than 2 can fit the data for the cotangent of the Hall angle. As we show here, this is indistinguishable from the ratio of terms each of which is a sum of contributions with different temperature dependences.

Our results have two significant conclusions. First, anisotropies, especially in the scattering rate, introduce corrections to the calculation of magnetotransport properties so that conventional results such as $`R_H=(nec)^1`$, or $`\mathrm{tan}\theta _H=\omega _c\tau _{tr}`$ are not in general valid.hlu1  Second, the fact that the Hall effect in the cuprates can be understood within transport theory, with only the measured single-particle relaxation rate, the measured resistivity and Fermi surface quantities as inputs implies that no new physics is involved other than what leads to the behavior of those quantities. We conclude that the MFL, which gives the experimentally observed temperature dependence of the relaxation rate, is sufficient to account for all the the anomalous magnetotransport near optimal doping. This point of view is reinforced by recent frequency-dependent complex Hall effect measurementsmatt  whose results are consistent with the analysis presented here, as discussed in Sec. VI.

## II Transport Equation

At long wavelengths and low frequencies the derivation of Boltzmann equation relies only on conservation laws and is valid whether the system is a Fermi liquid or not. Therefore we start with the linearized Boltzmann equation for the deviation $`g(𝐤,t)`$ of the momentum distribution function from its equilibrium value $`f(𝐤)`$ in the presence of uniform static electric and magnetic fields,

$$\frac{g_𝐤}{t}+e𝐄𝐯_𝐤\frac{f}{ϵ_𝐤}+\frac{e}{\mathrm{}c}(𝐯_𝐤\times 𝐁)_𝐤g_𝐤=C_𝐤$$
(3)
Here $`C(𝐤)`$ is the collision operator. The solution to this equation for the stationary case ($`\omega =0`$) is given by KSV

$$g(𝐤)=e\mathrm{}\underset{𝐤^{}}{}A_{𝐤,𝐤^{}}^1\left[𝐄𝐯_𝐤^{}\left(\frac{f}{ϵ_𝐤^{}}\right)\right],$$
(4)
where

$$A_{𝐤,𝐤^{}}=\mathrm{}\left[\frac{1}{\tau (𝐤)}+\frac{e}{\mathrm{}c}𝐯_𝐤\times 𝐁_𝐤\right]\delta _{𝐤,𝐤^{}}C_{𝐤,𝐤^{}}.$$
(5)
Here $`C_{𝐤,𝐤^{}}`$ is the “scattering-in” term in the collision operator for the Boltzmann equation and $`\mathrm{}/\tau (𝐤)=_𝐤^{}C_{𝐤,𝐤^{}}`$, is the “scattering-out” term equal to the single-particle relaxation rate. It is evident that the distribution $`g(𝐤)`$ is not only determined by the energy of the state $`𝐤`$ but also by the anisotropy of the scattering. The distribution is depleted in directions of large net scattering and augmented in directions of small net scattering.

We calculate the conductivities using the single-particle scattering rate of Eq. (2). The conductivity tensor is

$$\sigma ^{\mu \nu }=\frac{e^2\mathrm{}}{\mathrm{\Omega }}\underset{𝐤,𝐤^{}}{}v_{\mu ,𝐤}A_{𝐤,𝐤^{}}^1v_{\nu ,𝐤^{}}\left(\frac{f}{ϵ_𝐤^{}}\right),$$
(6)
where $`\mathrm{\Omega }`$ is the sample (taken to be a plane) area. We expand $`A^1`$ from Eq. (4) in powers of $`𝐁`$: KSV

$`A^1=𝒯`$ $``$ $`(e/c)𝒯(𝐯\times 𝐁)𝒯`$ (7)
$`+`$ $`(e/c)^2𝒯(𝐯\times 𝐁)𝒯(𝐯\times 𝐁)𝒯,`$
where (from Eq. (2.12) of Ref. 17)

$$𝒯_{𝐤,𝐤^{}}=\frac{1}{\mathrm{}}[\tau _𝐤\delta _{𝐤,𝐤^{}}+\underset{𝐤^{\prime \prime }}{}\tau _𝐤C_{𝐤,𝐤^{\prime \prime }}𝒯_{𝐤^{\prime \prime },𝐤^{}}]$$
(8)
The first term in $`A^1`$ gives the longitudinal conductivity, the second the Hall conductivity $`\sigma ^{xy}`$ and the magnetoconductivity may be calculated from the third term.

As discussed above, the $`\tau _𝐤`$ is the MFL scattering plus the angle dependent impurity piece:

$$1/\tau _𝐤=1/\tau _M+1/\tau _i(\widehat{k}).$$
(9)
Our approach to solving Eq. (7) is to take advantage of the properties of $`1/\tau _M`$ and $`1/\tau _i`$. The kernel $`C`$ in Eq. (8) comprises the vertex corrections for the various scattering mechanisms. Within MFL, the $`1/\tau _M`$ is an $`s`$-wave scattering process so the MFL interaction does not appear in $`C`$. That leaves

$$C_{𝐤,𝐤^{}}=2\pi \delta (ϵ_𝐤ϵ_𝐤^{})|u_i(\theta ,\theta ^{})|^2,$$
(10)
where, e.g., $`\theta `$ is the angle of $`\widehat{k}`$ along the Fermi surface. We have argued prl ; pnas  that the impurity scattering matrix element $`u_i(\theta ,\theta ^{})`$ involves scattering through small angles only. This property allows an expansion of the RHS of Eq. (8) in powers of a “small angle scattering parameter”, $`\theta _c`$ which can itself be a function of $`\theta `$. We illustrate such an expansion for the calculation of $`1/\tau _i`$:

$`{\displaystyle \frac{1}{\tau _i}}(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{𝐤^{}}{}}|u_i(\theta ,\theta ^{})|^2\delta (ϵ^{}ϵ)`$ (11)
$``$ $`\theta _cU(\theta ,\theta )N(\theta )+{\displaystyle \frac{\theta _c^3}{24}}\left[{\displaystyle \frac{d^2}{d\theta ^2}}U(\theta ,\theta ^{})N(\theta ^{})\right]_{\theta ^{}=\theta }.`$
Here we have used

$$\underset{𝐤}{}=\frac{1}{2\pi }𝑑ϵN(\theta )𝑑\theta ,N(\theta )=\frac{\mathrm{\Omega }}{2\pi \mathrm{}}\frac{dk_t/d\theta }{v(\theta )}.$$
(12)
$`N(\theta )`$ is the density of states per unit energy at the Fermi surface at angle $`\theta `$ and $`v(\theta )`$ is the Fermi velocity. $`dk_t`$ is an infinitesimal taken tangent to the Fermi surface at the angle $`\theta `$. We have abbreviated $`|u_i(\theta ,\theta ^{})|^2`$ by $`U(\theta ,\theta ^{})`$. We shall describe below how experiments show that $`\theta _c`$ is indeed sufficiently small to justify such an expansion.

## III Longitudinal Conductivity

We begin by determining the longitudinal conductivity within the model. We first give the solution for a general Fermi surface and then specialize to a particular example to illustrate the results.

From Eqs. (5,6), we have

$$\sigma ^{xx}=\frac{e^2\mathrm{}}{\mathrm{\Omega }}\underset{\mathrm{𝐤𝐤}^{}}{}v_𝐤^x𝒯_{𝐤,𝐤^{}}v_𝐤^{}^x(\frac{f}{ϵ^{}})$$
(13)
Define the vector L as

$$𝐋_𝐤=\underset{𝐤^{}}{}𝒯(𝐤,𝐤^{})(f/ϵ^{})𝐯_𝐤^{},$$
(14)
so that

$$\sigma ^{xx}=\frac{e^2\mathrm{}}{\mathrm{\Omega }}\underset{𝐤}{}v_𝐤^xL_𝐤^x.$$
(15)
Using Eq. (8) in Eq. (14), we find the integral equation for $`𝐋`$:

$$𝐋_𝐤=\frac{\tau _𝐤}{\mathrm{}}\left[(\frac{f}{ϵ})v_𝐤^x+\underset{𝐤^{\prime \prime }}{}C_{𝐤,𝐤^{\prime \prime }}𝐋_{𝐤^{\prime \prime }}\right].$$
(16)
As discussed earlier, $`C`$ contains only the small angle scattering so that

$$𝐋_𝐤=\frac{\tau _𝐤}{\mathrm{}}\left[(\frac{f}{ϵ})𝐯_𝐤+𝑑\theta ^{}N(\theta ^{})U(\theta ,\theta ^{})𝐋_𝐤^{}\right],$$
(17)
which, since $`𝐋`$ is restricted to the Fermi surface, can be rewritten as

$`{\displaystyle \frac{\mathrm{}}{\tau (\theta )}}𝐋(\theta )`$ $`=`$ $`({\displaystyle \frac{f}{ϵ}})𝐯(\theta )+{\displaystyle \frac{\mathrm{}}{\tau _i(\theta )}}𝐋(\theta )`$ (18)
$`+`$ $`{\displaystyle 𝑑\theta ^{}N(\theta ^{})U(\theta ,\theta ^{})[𝐋(\theta ^{})𝐋(\theta )]},`$
where $`1/\tau _i=1/\tau 1/\tau _M`$ is given in Eq. (11). Given the small-angle scattering restriction on $`U(\theta ,\theta ^{})`$, we expand the difference $`𝐋(\theta ^{})𝐋(\theta )`$ and find the differential equation (primes indicate derivatives with respect to $`\theta `$)

$$𝐌(\theta )=𝐯(\theta )+u(\theta )𝐌^{\prime \prime }(\theta )+2u^{}(\theta )𝐌^{}(\theta ),$$
(19)
where

$$𝐌(\theta )=\frac{\mathrm{}}{\tau _M(f/ϵ)}𝐋(\theta ),$$
(20)
and

$`u(\theta )`$ $`=`$ $`(\theta _c^3/24)(\tau _M/\mathrm{})N(\theta )U(\theta ,\theta ),`$
$`u^{}(\theta )`$ $`=`$ $`(\theta _c^3/24)(\tau _M/\mathrm{})[{\displaystyle \frac{d}{d\theta ^{}}}N(\theta ^{})U(\theta ,\theta ^{})]_{\theta ^{}\theta }.`$ (21)
This differential equation can be solved for the components of $`𝐌`$, hence $`𝐋`$, once the $`\theta `$-dependences of $`𝐯,u`$ and $`u^{}`$ are known. This equation, which is basic to the evaluation of both the longitudinal and the Hall conductivities is completely equivalent to Eq. (10) of Ref. 5. Although $`u`$ contains a small parameter an iterative solution of Eq. (19), as in our previous work,prl  is not valid. We solve Eq. (19) exactly in this paper.

As we discussed previously,prl ; err  the details of the anisotropies of the quantities entering the transport coefficients determine the magnitudes of the various contributions, in particular to $`\sigma ^{xy}`$ in the presence of a magnetic field. The shape of the Fermi surface is especially important since it determines the size of the contribution to $`\sigma ^{xy}`$ from the isotropic marginal Fermi liquid scattering rate $`1/\tau _M`$. In what follows, we take the anisotropy of the impurity scattering from ARPES data and we assume a simple form for the Fermi surface velocity.

A schematic of the Fermi surface is shown in Fig. 1, where the angular variable $`\theta `$ is shown. A general Fermi surface respecting the square symmetry of the CuO2 layers must have Fermi velocities in the first quadrant of the form

$`v_x`$ $`=`$ $`{\displaystyle \underset{n}{}}v_n\mathrm{sin}[(2n+1)\theta ]`$
$`v_y`$ $`=`$ $`{\displaystyle \underset{n}{}}v_n(1)^n\mathrm{cos}[(2n+1)\theta ].`$ (22)
For our calculations we take a form for $`𝐯(\theta )`$ which is the simplest extension beyond the circular Fermi surface (similar to that shown in Fig. 1) consistent with Eq. (22). In the first quadrant,

$`v_x`$ $`=`$ $`v_0(\mathrm{sin}\theta +\rho \mathrm{sin}3\theta )`$
$`v_y`$ $`=`$ $`v_0(\mathrm{cos}\theta \rho \mathrm{cos}3\theta )`$ (23)
For this choice for $`𝐯(\theta )`$, the density of states of Eq. (12) in all quadrants is

$$N(\theta )=\frac{\mathrm{\Omega }}{2\pi \mathrm{}}\frac{k_0(1\rho )^{1/4}}{v_0(1\rho \mathrm{cos}4\theta )^{5/4}},$$
(24)
where $`k_0`$ is the value of $`k_F(\theta )`$ at $`\theta =0`$.

Thus, according to Eq. (15), $`\sigma ^{xx}`$ is given by:

$$\sigma ^{xx}=\frac{e^2\tau _M(1\rho )^{1/4}}{\pi ^2\mathrm{}}\frac{k_0}{v_0}_0^{\pi /2}𝑑\theta \frac{v_x(\theta )M_x(\theta )}{(1\rho \mathrm{cos}4\theta )^{5/4}}.$$
(25)
## IV Hall Conductivity

For a magnetic field perpendicular to the $`xy`$ plane, $`𝐯\times 𝐁=[\mathrm{\Omega }/2\pi \mathrm{}N(\theta )]B/\theta `$ and from Eqs. (6,7),

$$\sigma ^{xy}=\frac{e^3B}{2\pi c}\underset{𝐤,𝐤^{},𝐤1}{}v_𝐤^x\frac{𝒯(𝐤,𝐤_1)}{N(\theta _1)}\frac{d𝒯(𝐤_1,𝐤^{})}{d\theta _1}(\frac{f}{ϵ^{}})v_𝐤^{}^y.$$
(26)
We rewrite this using Eq. (14):

$$\sigma ^{xy}=\frac{e^3B}{2\pi c}\underset{𝐤,𝐤1}{}v_𝐤^x\frac{𝒯(𝐤,𝐤_1)}{N(\theta _1)}\frac{d}{d\theta _1}L_{𝐤1}^y.$$
(27)
We define

$$𝐊_𝐤=\underset{𝐤^{}}{}𝒯(𝐤,𝐤^{})\frac{1}{N(\theta ^{})}\frac{d}{d\theta ^{}}𝐋_𝐤^{},$$
(28)
so that

$$\sigma ^{xy}=\frac{e^3B}{2\pi c}\underset{𝐤}{}v_𝐤^xK_𝐤^y.$$
(29)
Using Eq. (16) in Eq. (31), we get an integral equation for $`𝐊`$:

$`𝐊_𝐤`$ $`=`$ $`{\displaystyle \frac{\tau _𝐤}{\mathrm{}}}\left\{{\displaystyle \frac{1}{N(\theta )}}𝐋_{𝐤}^{}{}_{}{}^{}+{\displaystyle \underset{𝐤^{\prime \prime }}{}}C_{𝐤,𝐤^{\prime \prime }}𝐊_{𝐤^{\prime \prime }}\right\}`$ (30)
$`=`$ $`{\displaystyle \frac{\tau _𝐤}{\mathrm{}}}\left\{{\displaystyle \frac{𝐋_{𝐤}^{}{}_{}{}^{}}{N(\theta )}}+{\displaystyle 𝑑\theta ^{}N(\theta ^{})U(\theta ,\theta ^{})𝐊(\theta ^{})}\right\}.`$
We proceed exactly as in the analysis for $`𝐋`$, Eqs. (16-21). For the required $`K_y`$, we find the differential equation

$$Z(\theta )=M_y^{}(\theta )(1\rho \mathrm{cos}4\theta )^{5/4}+u(\theta )Z^{\prime \prime }(\theta )+2u^{}(\theta )Z^{}(\theta )$$
(31)
where $`Z(\theta )`$ determines $`K_y(\theta )`$ as

$$K_y(\theta )=\left(\frac{\tau _M}{\mathrm{}}\right)^2\frac{2\pi }{\mathrm{\Omega }}\frac{\mathrm{}v_0}{k_0}\frac{1}{(1\rho )^{1/4}}(\frac{f}{ϵ})Z(\theta ).$$
(32)
and $`M_y`$ is obtained as the solution of Eq. (19). As in that case, an iterative solution of Eq. (31) is not valid. Combining all our results into Eq. (29), we find

$$\sigma ^{xy}=\frac{e^3B}{\pi ^2c}\left(\frac{\tau _M}{\mathrm{}}\right)^2_0^{\pi /2}\frac{d\theta }{(1\rho \mathrm{cos}4\theta )^{5/4}}v_x(\theta )Z(\theta )$$
(33)
When there is no small-angle scattering, $`M_y=v_y`$ and $`Z(\theta )=v_0(\mathrm{sin}\theta 3\rho \mathrm{sin}3\theta )(1\rho \mathrm{cos}4\theta )^{5/4}`$. In that case,

$$\sigma ^{xy}_0^{\pi /2}𝑑\theta (\mathrm{sin}^2\theta 3\rho ^2\mathrm{sin}^2\theta )$$
(34)
which vanishes when $`\rho =1/\sqrt{3}`$. This is the familiar result of the vanishing of $`\sigma ^{xy}`$ when the Fermi surface consists of equal portions of positive and negative curvature. One sees that the anisotropic small angle contributions, which have a leading contribution proportional to $`\tau _M^3`$, i.e. to $`1/T^3`$, do not vanish there.

## V Choice of Parameters and Comparison with Experiments

We now describe the evaluation of the conductivities using experimental data from ARPES and from the longitudinal transport in zero magnetic field. The strategy is as follows. The anisotropies in the problem will be determined from the anisotropy of the impurity scattering $`1/\tau _i(\theta )=2\gamma (\widehat{k})/\mathrm{}`$ as determined from ARPES. For simplicity, we shall let the density of states $`N(\theta )`$ be responsible for the anisotropy of $`1/\tau _i`$, and hence for the quantity $`u(\theta )`$, see Eq. (21). The measured anisotropy of $`1/\tau _i`$ from the $`\pi ,\pi `$ direction to the $`\pi ,0`$ direction in the Brillouin zone then determines the velocity anisotropy parameter $`\rho `$ from Eqs. (11,24). These parameterizations give the correct behavior at the edges of the Brillouin zone. The small angle parameter $`\theta _c`$ can then be found from the “residual resistance” ratio ($`RR`$) of the resistivity at $`T_c`$ to the extrapolated value at $`T=0`$. That is sufficient to determine the $`T`$-dependence of $`\mathrm{cot}\theta _H=\sigma ^{xx}/\sigma ^{xy}`$. The magnitude of $`\mathrm{cot}\theta _H`$ is then fixed by the effective mass $`\mathrm{}k_0/v_0`$, which we leave as an adjustable parameter. It will be seen that the effective mass has a quite reasonable value.

The conductivities are given by Eqs. (25,33). The quantity $`u`$ which appears in the differential equations for $`𝐌`$ and $`Z`$ is defined as $`u(\theta ;T)=(\theta _c^2/24)\tau _M/\tau _i(\theta )=(\tau _M/\mathrm{})(\theta _c^3/24)N(\theta )U(\theta ,\theta )`$. The measured anisotropy pnas  of $`\gamma (\widehat{k})=\mathrm{}/2\tau _i`$ is about 4.5. We assign this to $`N(\theta )`$, thus determining $`\rho `$ from Eq . (24) as $`\rho =0.54`$. Therefore, we have

$`1/\tau _i(\theta )`$ $`=`$ $`(1/\tau _0)(1\rho \mathrm{cos}4\theta )^{5/4},`$
$`u(\theta ;T)`$ $`=`$ $`u_0(T)(1\rho \mathrm{cos}4\theta )^{5/4},`$
$`\rho `$ $`=`$ $`0.54.`$ (35)
The quantity $`u_0(T)=\tau _M(T)(\theta _c^2/24)/\tau _0`$. $`\tau _M(T)`$ and $`\tau _0`$ may be determined from ARPES as described in Ref. 3. From ARPES, $`\mathrm{}/\tau _i0.24`$ eV at $`\theta =\pi /8`$. So $`\mathrm{}/\tau _0=0.24`$ eV. From the same source, $`\mathrm{}/\tau _M=0.015(T/100)`$ eV, where $`T`$ is in Kelvin. This determines $`u_0`$ as $`u_0(T)=16(100/T)(\theta _c^2/24)`$. The remaining parameter $`u_0`$, that is $`\theta _c`$, is determined from the resistance ratio $`RR`$ as follows.

The longitudinal resistivity in zero magnetic field $`\rho ^{xx}`$ is given by the inverse of $`\sigma ^{xx}`$ from Eq. (25). Unfortunately we do not have an analytic solution of the differential equation Eq. (19) for $`M_x`$, although the coefficients are now known from Eq. (35). Numerical integrations of Eqs. (19,25) are performed for different values of $`u_0(T)=(1600/T)(\theta _c^2/24)`$. For a given value of $`\theta _c`$, $`\rho ^{xx}`$ is almost precisely linear in $`T`$ above $`T=100`$ K. See Fig. 2. The $`T=0`$ intercept of $`\rho ^{xx}`$, $`\rho _0`$, is obtained by extrapolating the found high temperature linear dependence. The $`RR`$ is defined as $`\rho ^{xx}(T=100)/\rho _0`$, and is typically about 8, say. This is obtained for $`u_0(T).022(100/T)`$ so that the characteristic small angle parameter is indeed small, $`\theta _c^2/240.0014`$.

All parameters are now fixed with the exception of $`m^{}=\mathrm{}k_0/v_0`$. The numerical integrations of Eqs. (19,31,33) for $`M_y,Z`$ and $`\sigma ^{xy}`$ respectively can be carried out for different values of $`T`$ and combined with the previous result for $`\sigma ^{xx}(T)`$ to give $`\mathrm{cot}\theta _H(T)=\sigma ^{xx}(T)/\sigma ^{xy}(T)`$:

$$\mathrm{cot}\theta _H(T)=\frac{0.82}{\omega _c\tau _M}\frac{m^{}}{m}\frac{𝑑\theta n(\theta )v_x(\theta )M_x(\theta )}{𝑑\theta n(\theta )v_x(\theta )Z(\theta )},$$
(36)
where we have used $`\rho =0.54`$, $`\omega _c`$ is the cyclotron frequency with the bare mass and $`n(\theta )=(10.54\mathrm{cos}4\theta )^{5/4}`$.

We now compare our result with the experiment of Chien, et al.ong  They measured the Hall angle as a function of temperature for various concentrations of Zn impurities in near-optimally doped YBCO. Their data at $`B=8`$ T ($`\omega _c=1.28\times 10^{12}`$ s-1) is reproduced in Fig. 3.

We examine the data for zero Zn concentration. In Fig. 4, we have replotted the $`x=0`$ points from Fig. 3 (large dots) and the theoretical result (small dots) of Eq. (36) with $`m^{}/m=1.5`$, so adjusted to give the best agreement. It is important that $`m^{}`$ turns out to be reasonable. We emphasize that the parameters we picked were determined from ARPES experiments on BiSCO, while the Hall data is on YBCO. Thus, our intent here is only to show that the scattering rates characteristic of high-$`T_c`$ superconductors in the normal state near optimal doping can by themselves account for the temperature dependence of the Hall angle, not to demonstrate quantitative agreement with a particular experiment. It is seen that the experimental temperature dependence is fairly well reproduced.

Notice that the experimental curves have a negative curvature on a $`T^2`$ plot. While this might suggest a power law less than 2 and several investigators have tried such a fit,ando ; raffy  our picture indicates that this behavior is a consequence of the fact that the conductivities are each the sum of two terms with approximate $`1/T^n`$ and $`1/T^{n+1}`$ ($`\sigma ^{xx}:n=1,\sigma ^{xy}:n=2`$) dependences.

The data of Ref. 6 and our theory both suggest that in the absence of in-plane impurities such as Zn, the extrapolated $`\mathrm{cot}\theta _H`$ at $`T0`$ is zero. In-plane impurities will in general not give scattering restricted to small angles. Indeed, their influence is evident in the measured resistivity. We could include their effects by assuming them to be isotropic scatterers which add a constant $`1/\tau _z`$ to the isotropic linear in temperature MFL rate $`1/\tau _M`$. This has the effect of an almost parallel upward shift of $`\mathrm{cot}\theta _H`$ so that it has a zero temperature intercept proportional to $`1/\tau _z`$, i.e. to the in-plane impurity concentration. We illustrate this qualitatively using some data from Ref. 9. In Fig. 5, we show the data (large dots) for Bi2Sr2-xLaxCuO6 with $`x=0.44`$ (optimal doping). It is seen that the data does not extrapolate with a zero intercept. We interpret this as signalling the presence of in-plane scatterers. The theoretical curve (small dots) uses similar parameters as in Fig. 4 but with an in-plane impurity contribution added to $`1/\tau _M`$.

## VI Complex Hall conductivity

In this section, we comment on the recent ac Hall effect results of Grayson, et al.matt  For low frequencies and $`T>T_c`$, we can replace $`1\tau _M1/\tau _{tr}`$ by $`1/\tau _{tr}i\omega `$. This follows immediately from Eqs. (3-5): the distribution function $`g(𝐤,\omega )`$ and therefore the transport properties are obtained from the $`\omega =0`$ results by this replacement.

From Fig. 2 it is seen that the resistivity is almost precisely linear in $`T`$ in the normal state: $`\rho ^{xx}=a+bT`$, where $`a`$ is the (small) residual resistivity due to the small-angle impurity scattering. Thus, $`a/b`$ is about 10, say. We have also found that $`\sigma ^{xy}(T)`$ is, in the temperature range of interest, very closely of the form $`c/T^2+d/T^3`$, with $`c/d0.01`$K-1. Therefore, neglecting quantities of order 0.1, our prediction for the Hall angle is

$`\theta _H(\omega ,T)`$ $``$ $`\mathrm{tan}\theta _H=\sigma ^{xy}\rho ^{xx}`$ (37)
$``$ $`A\left[{\displaystyle \frac{1}{(1/\tau _{tr}i\omega )^2}}+{\displaystyle \frac{c}{d}}{\displaystyle \frac{1}{(1/\tau _{tr}i\omega )}}\right],`$
where $`A`$ and $`c/d`$ are constant in temperature. The latter ratio is now in units of seconds; for our example parameters it is about $`2\times 10^{14}`$s. We emphasize that a single relaxation rate $`1/\tau _{tr}`$ enters all our expressions. The constant $`A`$ in Eq. (37) corresponds to the quantity $`\omega _H\mathrm{\Omega }_p`$ of Ref. 15. In that reference, Fig. 4 shows that $`A`$ is indeed temperature and frequency independent. The $`c/d`$ term is the conventional term. In Ref. 15, it was shown that by itself it cannot account for the experimental data. Just as the dc data shows deviations from $`T^2`$ behavior for $`\mathrm{cot}\theta _H`$, we expect that the ac data should be fit by a combination of the two terms in Eq. (37). At higher frequencies or temperatures, the conventional term should finally dominate. The complex Hall angle measurements convincingly demonstrate that just one ineleastic transport rate determines all the frequency and temperature dependence of the transport properties of the high-$`T_c`$ superconductors in the marginal Fermi liquid region, i.e. in their normal state near optimum doping.

When calculating the ac conductivity, the limits of applicability of the Boltzmann equation should be kept in mind. The results are only valid for $`\omega \tau 1`$. In this regime the conservation laws (in terms of bare particles) completely determine transport and the Boltzmann equation deals with them properly. In a microscopic theory, the effects of small-angle scattering calculated here appear as corrections to the coupling of the external magnetic field to the carriers. At high enough frequency these vertex corrections must vanish. The crossover frequency and the behavior of transport properties in the intermediate regime can only be determined by a microscopic calculation. Experimental results at high frequencies do depart from the predictions of Eq. (37).

## VII Concluding Remarks

Here we discuss sundry issues related to the theory and calculations presented in this paper. We emphasize that this work pertains only to the Hall effect in optimally or overdoped samples. For the latter, the linear-in-temperature scattering rate $`1/\tau _M`$ should be replaced in the theory by a rate with the $`T`$-dependence of the observed resistivity. In underdoped cuprates, additional physical considerations appear to determine the transport and equilibrium properties.

### VII.1 Magnetoresistance

Earlier,prl  we gave a plausibility argument that the observed magnetoresistance would follow quantitatively from the same considerations which lead to the observed Hall effect. We leave it as a future exercise to calculate this quantity directly from the solution of the Boltzmann equation as an extension of the present theory.

### VII.2 Small-angle scattering, Fermi surface geometry, bilayer splitting …

We have modeled the Fermi velocity through Eq. (23), which with $`\rho 0.54`$ is in itself enough to give the observed angular variation of the elastic part of the single-particle relaxation rate. As we have seen, this also gives agreement with the observed temperature dependence of the Hall angle. This value of $`\rho `$ is 10% away from the value at which the customary Hall coefficient would be zero. This seems to us reasonable, since in a model with nearest neighbor Cu-O hopping alone, the Hall effect is zero at half-filling and with O-O hopping included, the density for zero Hall effect shifts substantially towards hole doping. The optimal doping composition lies in the range 15-20% hole doping. In this connection, it is important to remember that we are using parameters obtained from experiments on BSCCO to fit Hall data on YBCO. Therefore the precise numbers are not too meaningful; we only argue that we have shown that the small-angle scattering scenario gives the qualitative features of the experiments.

In our numerical analysis, we have chosen the simplest parameterization of the Fermi velocity to show the plausibility of the small-angle effect in relation to the experiments. Our parameterization should not be expected to represent the Fermi surface adequately. For the actual Fermi surface several terms in Eq. (22) for the Fermi velocity would have to be included with coefficients smaller than $`\rho `$. We note however that the successive terms give compensatingly larger coefficients since the higher harmonics produce larger derivative terms in Eq. (36) for $`\sigma _{xy}`$. For a truer Fermi surface the advantage of analytical calculations is lost. These remarks are of academic interest at the present time since the actual Fermi surface is not accurately known.

More relevant is the fact that the Fermi surface of Bi-2212, the ARPES data for which we have used here, has two sheets coming from the bilayer splitting. The bilayer splitting varies as a function of angle, being largest in the $`\pi ,0`$ direction. This arises because of the geometry of the interlayer orbitals.ok  This angular variation is essentially the same as that of the observed anisotropic elastic scattering. The fact that the bilayer splitting is not resolved in the experimental data we have used indicates the presence of an elastic scattering mechanism which couples the layers. We surmise that for interlayer or intralayer scattering due to impurities between the planes, the same orbitals are involved and therefore the same angular dependence arises in scattering. In a single-layer model as treated here, this effect is parameterized through the choice of the angle-dependence of the Fermi velocity as in Eq. (23).

### VII.3 Comments on related work

As already remarked, although the analytical expressions derived in our earlier workprl  are correct, we erred, as pointed out by V. Yakovenko, in our choice of the parameterization in which they were evaluated. This error has been remedied here. Our more complete analysis with a simpler parameterization in fact gives results consistent with experiment. We have shown that not only do we get the right temperature dependence but that with a reasonable effective mass, the absolute value of the cotangent of the Hall angle is obtained within 10% of the data from the same parameterization.

In Ref. 19, R. Hlubina argues that the new term discovered by us, if adequate to explain the Hall angle, would make an unacceptable contribution to the longitudinal resistivity. Since our parameters are actually determined by the measured resistivity, this criticism is invalid. It is difficult to make a direct comparison with the calculations of Ref. 19 since a quite different parameterization of the Fermi surface is used there. However, we can note that the parameterization family used in Hlubina’s evaluation does not include ours and gives a much larger conventional contribution to $`\sigma ^{xy}`$ than the one we have used.

A recent paper of Carter and Schofieldsch  addresses the main point of our original paper.prl  Their work consists of two parts. One is analytical, the other is a numerical solution of equations we derived in Ref. 5. We cannot comment on the adequacy of the numerical work, but Carter and Schofield are indeed correct that the customary term (with a different temperature dependence than the new term) is zero only for a particular choice of Fermi surface. This is of course well-known: the traditional Hall angle is zero only for “particle-hole symmetry” defined in terms of the curvature of the Fermi surface. However, as we have shown here, it is not necessary that the customary term be zero (and all the temperature dependence come from our new term) to get good agreement with experiment. It is sufficient that the customary term is about a factor of four smaller than the new term at $`T100K`$. In fact, the experimental data for $`\mathrm{cot}\theta _H`$, when plotted against $`T^2`$, usually shows a slight downward curvature. This point has been discussed in Sec. I and is adequately demonstrated in the experimental data shown in Fig3. 3-5.

The analytical calculation of Ref. 20, expressed in their Eqs. (7-10), is not sufficiently general to contain the anisotropies required to show the effect we have derived. Their calculation for $`\sigma ^{xy}`$ is equivalent to retaining only the impurity scattering contribution to the conventional term, not the new contribution which we have derived. See footnote 22 of Ref. 20 for a comment on this. Actually, the analytical calculation of Ref. 20 is based on a circular Fermi surface and we agree that such a parameterization never leads to a large enough effect.

## VIII Acknowledgments

We thank V. Yakovenko, D. Drew and M. Grayson for helpful discussions of our analysis and A.J. Schofield for informative communications. N.P. Ong and Y. Ando generously provided experimental data. EA was supported in part by NSF grant DMR99-76665. The authors acknowledge the hospitality of the Aspen Center for Physics where part of this work was carried out.


format: text

text:
# A Free, Fast, Simple and Efficient TVD MHD Code

## 1  Introduction

Astrophysical fluids in which the magnetic field plays an important role are common in nature. As a few examples, consider magnetized interstellar gas, accretion disks, molecular clouds, and jets. With the advent of high speed computers and ever improving MHD codes, considerable theoretical progress has been made through numerical simulation of otherwise intractable problems.

A major challenge to solving flux conservative systems of equations, such as ideal fluids and MHD, is the spontaneous development of shock discontinuities. Finite differencing across discontinuities leads to divergences and instabilities. Modern codes implement various aspects of “flux limiters” (Harten 1983) to achieve stability near shocks and second order accuracy away from shocks. Recently, several shock capturing methods which solve the MHD equations in flux conservative form with upwind finite differencing have been developed. Enforcing the $`\mathbf{}𝑩=0`$ constraint is key to the accuracy of these codes Near shock fronts, derivatives are ill-defined, and the divergence constraint can be maximally violated. Evans & Hawley (1988) first noted that this “constrained transport” (CT) can easily be enforced to machine precision by (1) defining the magnetic field at cell faces instead of centers, and (2) using the same EMF, computed on cell corners, to update the magnetic flux through each neighboring face. Using CT and various shock capturing schemes, several groups (see e.g., Tóth (2000) for a review of different methods) have now produced robust, efficient MHD codes.

The detailed algorithm used for the finite differencing varies between the different groups. We do not give an exhaustive review of the literature but rather compare only to widely used codes, or those similar to ours. The Zeus code (Stone & Norman 1992) partially updates certain fluid and magnetic field quantities along Alfven, but not magnetosonic, characteristics to avoid short lengthscale instabilities in shear Alfven waves. Ryu et al. (1998) use Harten’s TVD method, which evolves the fluid along all the characteristics by constructing the linearized eigenvectors. Common to these two methods is the need to first compute the EMF over the whole grid, then perform a spatial averaging of the EMF, and then update the magnetic field.

In this paper we implement the divergence constraint in a slightly different way than previous investigators. We show that individual pieces of the EMF can be used in advection-constraint steps, without the need to store the computed EMF’s over the whole grid. This gives us a sizeable savings in memory than if the EMF’s were stored. Furthermore, by using Jin and Xin’s (1995) “symmetric” 111 This method is called symmetric since it decomposes each fluid quantity into left and right moving parts, each of which can be differenced in an upwind fashion. method of computing TVD fluxes, we can reduce the operations count relative to codes which manifestly evolve the fluid along characteristics.

In section 2 we review the MHD equations. In section 3 we describe our numerical method. Tests of the code are presented in section 4. Section 5 contains a discussion of the merits and drawbacks of the code. Section 6 contains the conclusions. We briefly review the Jin & Xin (1995) method for solving one dimensional advection equations in an appendix.

## 2  Equations

The MHD equations expressing conservation of mass, momentum and energy, as well as magnetic flux freezing are (Landau & Lifshitz 1984)

$$_t\rho +\mathbf{}\left(\rho 𝒗\right)=0$$
(1)
$$_t\left(\rho 𝒗\right)+\mathbf{}\left(\rho 𝒗𝒗+P_{}𝜹𝒃𝒃\right)=\rho 𝒂$$
(2)
$$_te+\mathbf{}\left[(e+P_{})𝒗𝒃𝒃𝒗\right]=\rho 𝒗𝒂$$
(3)
$$_t𝒃=\mathbf{}\times \left(𝒗\times 𝒃\right)$$
(4)
$$\mathbf{}𝒃=0$$
(5)
$$P_{}=p+\frac{b^2}{2}$$
(6)
$$e=\frac{\rho v^2}{2}+\frac{p}{\gamma 1}+\frac{b^2}{2}.$$
(7)
Here $`\rho `$ and $`e`$ are the mass and (total) energy densities, $`𝒗`$ is the velocity, $`P_{}`$ is the total pressure, $`p`$ is the gas pressure, $`𝒃=𝑩/\sqrt{4\pi }`$ is the magnetic field in terms of $`\sqrt{4\pi }`$, $`𝒂`$ is an externally imposed acceleration, and $`𝜹`$ is the Kronecker delta symbol. In eq.7 we have used an ideal gas equation of state with internal energy $`\epsilon =p/(\gamma 1)`$, where $`\gamma `$ is the ratio of specific heats. The infinite conductivity limit has been used so that the electric field is $`𝑬=𝒗\times 𝑩/c`$. The electric force has been ignored since it is assumed that charge separation is negligible on the scales of interest.

## 3  Numerical method

First we describe the update of the magnetic field in two dimensional advection- constraint steps. We then briefly review the update of the fluid variables along one dimension. Next we discuss how operator splitting and Runge-Kutta can be used to make the code second order accurate in space and time. Finally, we discuss boundary conditions, fine tunings of the code, and the parallel implementation.

### 3.1  solution of the induction equation in 2D

We use operator splitting to reduce the problem into a series of smaller decoupled equations. Alternating the order of operators in the correct fashion allows one to achieve net second order accuracy. In this prescriptions, we hold the fluid variables fixed to update the magnetic field. The magnetic field is defined on cell faces (see fig.1) in order to satisfy the $`\mathbf{}𝒃=0`$ to machine precision. Let the cell centers be denoted by $`(i,j,k)(x_i,y_j,z_k)`$, and faces by $`(i\pm 1/2,j,k)`$, $`(i,j\pm 1/2,k)`$, and $`(i,j,k\pm 1/2)`$, etc. For convenience, let the cells have unit width. The magnetic field is then stored in arrays

$`bx(i,j,k)`$ $`=`$ $`b_{i1/2,j,k}^x`$
$`by(i,j,k)`$ $`=`$ $`b_{i,j1/2,k}^y`$
$`bz(i,j,k)`$ $`=`$ $`b_{i,j,k1/2}^z.`$ (8)
The flux out of cell $`(i,j,k)`$ is then

$`[{\displaystyle d^3x\mathbf{}𝒃}]_{ijk}`$ $`=`$ $`b_{i+1/2,j,k}^xb_{i1/2,j,k}^x`$ (9)
$`+`$ $`b_{i,j+1/2,k}^yb_{i,j1/2,k}^y`$
$`+`$ $`b_{i,j,k+1/2}^zb_{i,j,k1/2}^z.`$
Since the magnetic field is defined on cell faces, where the magnetic flux in a cell is evaluated, it is possible to enforce $`\mathbf{}𝒃=0`$ to machine precision (Evans & Hawley 1988). If we defined the magnetic field at some other location, the divergence could only be kept to zero to truncation error arising from interpolation to the faces. The truncation errors are necessarily large near shock discontinuites. Next we describe a method to evolve the field which preserves $`\mathbf{}𝒃=0`$, if it is so initially.

When the induction equation is written out in spatial components, it is apparent that the terms involving $`𝓔𝒗\times 𝒃`$ come in six pairs. For instance, the terms involving $`v_yb_x`$ are

$`_tb_x+_y(v_yb_x)`$ $`=`$ $`0`$
$`_tb_y`$ $`=`$ $`_x(v_yb_x).`$ (10)
The first equation is just the advection of $`b_x`$ along the $`y`$ direction, the second equation is a constraint which enforces $`\mathbf{}𝒃=0`$. The key point to note here is that to enforce $`\mathbf{}𝒃=0`$ we must use the same EMF computed in the advection step during the constraint step; otherwise $`\mathbf{}𝒃=0`$ will only be zero up to truncation error. We accomplish this by finding a second order accurate, upwind EMF $`v_yb_x`$ for the advection step to update $`b_x`$, and then immediately use this same EMF for the constraint step to update $`b_y`$. Ryu et al. (1998) first store the EMF’s over the entire 3D grid, then average the EMF’s, and then update the field. We construct the EMF using Jin and Xin’s (1995) TVD method, which is described in the appendix. Note that the velocity $`v_{ijk}^y`$ must be interpolated to the same position as the magnetic field $`b_{i1/2,j,k}^x`$ with second order accuracy.

Jin and Xin’s (1995) symmetric method introduces a “flux freezing speed” $`c`$, which must be greater than or equal to the maximum speed at which information can travel. Since we are holding the fluid variables fixed, the flux freezing speed for the advection-constraint equation is just $`c=|v_y|`$.

### 3.2  solution of the fluid equations in 1D

Now we briefly describe the fluid update. A more complete discussion is given in Trac & Pen (2003). The magnetic field is held fixed, and interpolated to grid centers with second order accuracy. Let $`𝒖=(u_1,u_2,u_3,u_4,u_5)=(\rho ,\rho v_x,\rho v_y,\rho v_z,e)`$ represent the volume averaged quantities positioned at the center of each cell. For advection along the $`x`$ direction, the Euler, continuity, and energy equations can be written in flux conservative form as

$`_t𝒖+_x𝑭`$ $`=`$ $`0`$ (11)
where the flux vector is given by

$$𝑭=\left(\begin{array}{c}\rho v_x\hfill \\ \rho v_x^2+P_{}b_x^2\hfill \\ \rho v_xv_yb_xb_y\hfill \\ \rho v_xv_zb_xb_z\hfill \\ (e+P_{})v_xb_x𝒃𝒗\hfill \end{array}\right)$$

and the pressure is determined by $`p=(\gamma 1)(e\rho v^2/2b^2/2)`$. We hold the magnetic field fixed during the fluid update, and interpolate $`𝒃`$ to cell centers for second order accuracy.

Eq.11 can be solved by symmetric TVD, described in the appendix. The flux freezing speed is taken to be $`c=|v_x|+(\gamma p/\rho +b^2/\rho )^{1/2}`$, which is the maximum speed information can travel 222 We have taken the maximum speed of the fast magnetosonic wave over all directions..

### 3.3  extension to 3D

Let the fluid update step for a time $`\mathrm{\Delta }t`$ along $`x`$ be denoted by “fluidx”, and update of $`b_x`$ along $`y`$ by “bxalongy”. Operator splitting requires us to apply each operator first in forward, and then in reverse order to advance by two timesteps. We implemented two versions. In one, we advance forward using the sequence of operations: fluidx, byalongx, bzalongx, fluidy, bxalongy, bzalongy, fluidz, bxalongz, byalongz, and then the reverse byalongz, bxalongz, fluidz, bzalongy, bxalongy, fluidy, bzalongx, byalongx, fluidx. A second implementation, used in the public version of this code, is to transpose the fluid variables and spatial dimensions (see, e.g. Press et al. (1996) pg.984). This is only easily done if two of the dimensions are equal. Transposing has the benefit of high efficiency on cache based computers, where we only need to read data in column order. It is convenient to implement the 3D code by using a single routine for the fluid update along the $`x`$ direction and a single routine for advection of $`b_y`$ along $`x`$, and constraint of $`b_x`$ along y. The order of the spatial indices is transposed to take account of the other directions. The advantage of having one subroutine is that it can be heavily optimized.

The time step is set by the fastest speed at which information travels over the grid. Since the fluid update is more restrictive, the time step is set to be

$`\mathrm{\Delta }t`$ $`=`$ $`\mathrm{cfl}\left[max(|v_x|,|v_y|,|v_z|)+(\gamma p/\rho +b^2)^{1/2}\right]^1`$ (12)
where cfl $`\mathrm{¡}\mathrm{}1`$ is generally set to cfl $`0.7`$ for stability.

### 3.4  boundary conditions

The standard boundary conditions such as periodic, continuous, or reflecting can be easily enforced by specifying values of the variables in “ghost zones” adjacent to the physical grid. These ghost zones are needed when interpolating $`𝒗`$ to cell faces, $`𝒃`$ to cell centers, and in the one-dimensional advection routines.

We have implemented the boundary conditions in two different ways. The first method is to pad the grid with a large number ($`610`$) of extra cells at each boundary. Both the on-grid and off-grid variables are evolved in time, but so many extra cells are used that the boundary cells only need be updated once per double time step. This method is useful for parallel implementations in which buffer zones are used to represent a small number of cells in adjacent regions. The second method to implement the boundary conditions requires that one write specific routines for interpolation or derivatives which specify the off grid values. This method requires less computation, and is preferable for serial applications. We find it convenient to evolve “extra” values of the magnetic field variables. That is, we evolve $`b_{i1/2,j,k}^x`$ for $`j=1,\mathrm{},n_y`$ and $`k=1,\mathrm{},n_z`$ but $`i=1,\mathrm{},n_x+1`$. This is useful for three reasons. First, this allows $`\mathbf{}𝒃`$ to be computed over all cells. Second, to update $`b_{n_x,j+1/2,k}^y`$ the fluxes $`v_yb_x`$ are needed at in the boundary cells (see fig.1) with $`i=n_x+1/2`$. Third, $`𝒃`$ can be interpolated to cell centers without the need to specify off grid values.

### 3.5  fine tunings of the code

Since the TVD limiters are nonlinear, sinusoidal waveforms can tend to become “clipped”, or boxy-looking. We find that these nonlinear distortions can be minimized by using constant flux freezing speed, set to be the maximum along that advection line. Additional stability can be gained by multiplying the flux freezing speed by a constant multiplicative factor, although this increases the number of time steps needed and makes the code more diffusive. Empirically we find that smoothing the velocity field which advects the magnetic field can lead to less damping of the slow mode.

In production runs, we have found an occasional failure of the code when the Courant condition is pushed too close to the limit. In the operator split approach, the time step is fixed at the beginning of a double time step, and determined from the Courant condition at the beginning. During the time step, this condition may change, leading to an instability if it exceeds the initial constraint. Our solution has been to be sufficiently conservative using a choice of cfl$`0.7`$. A more efficient procedure would be to measure the change in the Courant condition during the sweeps, and use this as an indicator in subsequent time steps. And should a given sweep step be instable, one can always break it into two substeps.

### 3.6  Parallel implementation

We have implemented a fully distributed version in MPI. After a full set of operators in one dimension, we update the buffer zones. For hydrodynamics, only 3 buffer cells are required. The magnetic field requires interpolation, and we use 16 buffer cells for magnetized simulations. A full three dimensional domain decomposition is implemented, where we update the buffers in the appropriate direction after each dimensional operator. Since only large faces are communicated, latency of communication is negligible, but signficant bandwidth is required to move the buffer zones. The communication is performed asynchronously, and computations proceed during the communication stage. Within each node, OpenMP is used to utilize multiple processors in a node without the overhead of buffer cells and communications.

We tested the parallel implementation on the CITA McKenzie beowulf cluster. The main cluster has 256 nodes of dual Intel Pentium-4 Xeon processors running at 2.4 Ghz, 1 GB of RAM, dual gigabit ethernet, and 80 GB of disk. The networking consists of bristles of 16 machines with one gigabit port connected to a switch. The second gigabit port is used to interconnect the bristles in a cubical layout. The nominal bi-section bandwidth is 128 Gbit/sec. For cubical problems, the usable bandwidth is higher when the domain decomposition is matched to the cube as it is for our runs. To minimize communications overhead, we mapped the computational grid layout to coincide with the physical network interconnect. The largest problem that we have been able to run in memory is $`1400^3`$ grid zones, which takes about 40 seconds per double time step. The fine grained OpenMP parallelism within each nodes allows the code to benefit from the intra-node hyperthread speedup. The code also fully vectorizes for the SSE2 parallel execution units. Due to the large number of buffer cells required, about 1/3 of the computation and memory are used by these buffers. The operation count of the van Leer limiter relaxing TVD algorithm is 33 floating point operations per variable per time step. The flux computation takes an additional 7 operations averaged per variable, for a total count of 40. Each double time step consists of 6 sweeps of 8 variables, or about 2000 floating point operations. Our execution speed corresponds to a sustained rate of over 200 Gflop on the cluster, which is about 5% of theoretical peak speed of 4.8 Tflop in single precision.

## 4  Tests

In this section we present tests of the code on MHD waves and shocks.

For all the tests we use $`\gamma =5/3`$ and set the box size $`L`$ to be equal in all directions. Period boundary conditions were used for the wave tests, and continuous boundary conditions (all variables continuous across the boundary) in the shock. The wave tests are two-dimensional while the shock tests are along one dimension. The van Leer limiter (see, e.g. Trac & Pen (2003)) and a constant freezing speed and cfl=0.7 (see section 3.5) were used throughout.

### 4.1  torsional Alfven Waves

Torsional Alfven waves are exact nonlinear solutions of the compressible MHD equations. In the absence of any perturbations or noise, they should propagate without steepening, making them a good test for numerical codes. 333 In the presence of any infinitesimal noise, Alfven waves are unstable to decay into three other waves (Goldstein 1978; Derby 1978). For large amplitude waves the growth time becomes comparable to the Alfven wave period. Hence care must be used in applying this test to very nonlinear waves as noise arising from truncation error or the nonlinear flux limiters may grow exponentially.

We perform tests for four different resolutions $`n_y=n_z=16,32,64,128`$, where $`n_{y,z}`$ are the number of grid points in the $`y`$ and $`z`$ directions. Different fluid pressures are used corresponding to low and high $`\beta =c_s^2/b^2=0.1,10`$. The exact solution we input to the code is $`\rho =1`$, $`e=p/(\gamma 1)+0.5+A^2`$,

$`b_x`$ $`=`$ $`A\mathrm{cos}(k(y+z)\omega t)`$ (13)
$`b_y`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(1+A\mathrm{sin}(k(y+z)\omega t)\right)`$ (14)
$`b_z`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(1A\mathrm{sin}(k(y+z)\omega t)\right)`$ (15)
$`v_x`$ $`=`$ $`A\mathrm{cos}(k(y+z)\omega t)`$ (16)
$`v_y`$ $`=`$ $`{\displaystyle \frac{A}{\sqrt{2}}}\mathrm{sin}(k(y+z)\omega t)`$ (17)
$`v_z`$ $`=`$ $`{\displaystyle \frac{A}{\sqrt{2}}}\mathrm{sin}(k(y+z)\omega t).`$ (18)
The wavenumber and frequency are $`k=2\pi /L`$ and $`\omega =\sqrt{2}k`$ for the lowest order mode. We set the ratio of wave to background field to be $`A=0.1`$ for the low gas pressure test, and $`A=1`$ for the high gas pressure test. As long as the thermal energy is larger than the kinetic energy of the wave, we can also test very non-linear waves 444For sufficiently small $`\beta 1`$ and large amplitude $`A1`$, we found the code is susceptible to a short lengthscale instability. The timescale over which this instability develops depends strongly on $`\beta `$ and $`A`$, and also the number of grid points. It is unclear to us whether or not this is the physical decay instability of Alfven waves due to perturbations seeded by the truncation error. .

The waves were propagated for one wave period. The result is read along the $`z`$ axis and plotted against the exact solution in fig.2. For clarity, we plotted two periods of the wave (only one period was simulated). The second order convergence is apparent from the figures: as one doubles the number of grid cells, the errors goes down by a factor of 4.

Even the significantly non-linear solutions are well behaved and also converge quadratically.

### 4.2 magnetosonic waves

We have tested magnetosonic waves in the linear regime $`\delta pp`$ for $`\beta =0.1`$ and $`10`$. The background magnetic field is $`𝒃=𝒆_z`$, background pressure is $`p=\beta `$, and the wavevector is $`𝐤=k(𝐞_y+𝐞_z)`$, where $`k=2\pi /L`$ and $`L=L_y=L_z`$. The exact solution we plug into the code is

$`\delta p`$ $`=`$ $`pA\mathrm{cos}(k(y+z)\omega t)`$
$`\rho `$ $`=`$ $`1+{\displaystyle \frac{\delta p}{c_s^2}}`$
$`\delta 𝒗`$ $`=`$ $`{\displaystyle \frac{\omega }{\omega ^2\omega _b^2}}\left(𝒌{\displaystyle \frac{\omega _b^2}{\omega ^2}}\widehat{𝒃}\widehat{𝒃}𝒌\right){\displaystyle \frac{\delta p}{\rho }}`$
$`\delta 𝒃`$ $`=`$ $`{\displaystyle \frac{k^2}{\omega ^2\omega _b^2}}\left(𝒃\widehat{𝒌}\widehat{𝒌}𝒃\right){\displaystyle \frac{\delta p}{\rho }}`$ (19)
where $`c_s^2=\gamma p/\rho `$, $`\omega _b^2=k^2b^2/\rho `$, and $`A1`$ is the wave amplitude. The fast and slow mode frequencies are given by $`\omega _{f,s}^2=.5(\omega _b^2+\omega _s^2)\pm 0.5\left((\omega _b^2+\omega _s^2)^24\omega _s^2\omega _a^2\right)^{1/2}`$ where $`\omega _s^2=c_s^2k^2`$ and $`\omega _a^2=(𝒌𝒃)^2/\rho `$. We evolved the waves for one period. The results for the fast and slow waves are shown in figures 3 and 4, respectively. The slow wave is subject to substantial diffusion as compared to the fast wave since its frequency is so much lower for these extreme values of $`\beta `$.

### 4.3  the MHD Riemann problem

These 1d tests of the code involve a shock tube along the x-axis, as in fig.2a,b of Ryu et al. (1998). We used continuous boundary conditions and 1024 grid points for both tests.

The initial conditions for fig.5 are $`(\rho ,v_x,v_y,v_z,p,b_x,b_y,b_z)=(1,10,0,0,20,5/\sqrt{4\pi },5/\sqrt{4\pi },0)`$ for the left side and $`(1,10,0,0,1,5/\sqrt{4\pi },5/\sqrt{4\pi },0)`$ for the right side. The code is run for a time $`0.08L`$. The result agrees well with fig.2a of Ryu et al. (1998). The following features can be seen. The steep discontinuities at $`x0.1`$ and $`x0.85`$ are fast shock fronts where the incoming flow converts its kinetic energy into thermal energy and compresses the transverse field $`b_y`$. As new matter falls on, this shock is regenerated and maintains it’s steep profile as it moves outward. At $`x0.6`$ and $`x0.5`$ are a slow shock and slow rarefaction, respectively. The slow shock again compresses the fluid but decreases the transverse field. At $`x0.55`$ the two phases of the initial gas configuration with different entropies form a contact discontinuity. Pressure, magnetic field and velocity are continuous while density and thermal energy experience a discontinuity. This discontuity moves rightward across the grid, and the TVD advection of such discontinuities results in some smearing or diffusion of the structure. No physical mechanism steepens this contact discontinuity once it smears, and a slow numerical diffusion is visible in this, and all generic TVD codes which do not introduce explicit contact steepeners. There are no significant oscillations. Both our and Ryu et al. (1998)’s solution have a slight overshoot in some variables in the first postshock cell, but this effect does not persist onto subsequent cells.

The initial condition for fig.6 has velocity and magnetic field components in all directions, and hence exhibits additional structures such as rotational discontinuities. The values are $`(\rho ,v_x,v_y,v_z,p,b_x,b_y,b_z)`$ = $`(1.08,1.2,0.01,0.5,0.95,2/\sqrt{4\pi },3.6/\sqrt{4\pi },2/\sqrt{4\pi })`$ on the left side and $`(1,0,0,0,1,2/\sqrt{4\pi },4/\sqrt{4\pi },2/\sqrt{4\pi })`$ on the right hand side. The code is run for a time $`0.2L`$. The results again agree with Ryu et al. (1998). The following features may be seen: fast shocks at $`x0.3`$ and $`0.9`$, rotational discontinuity at $`x0.53`$ right next to a slow shock at $`x0.55`$, contact discontinuity at $`x0.6`$, slow shock and rotational discontinuity at $`0.68`$ and $`0.70`$ respectively.

The shock tests have two basic types of structures: self-steepening shock fronts and non-evolutionary contact and rotational discontinuities. For the active shock fronts, fig.5 and 6 show the shocks extend over leading to a nominal shock resolution of around 2 grid cells. The internal contact discontinuities arise from the discontinuities in the initial conditions, and diffuse numerically as they advect over the grid. At discontinuities, the solution is non differentiable which is in general a challenge to numerical schemes. In the TVD approach, the scheme drops to first order accuracy, with some associated diffusivity. This is discussed in more detail in Trac & Pen (2003). The shock fronts are self steepening, so the first order diffusivity is less noticable.

## 5  Discussion of the Merits and Drawbacks of Our Numerical Scheme

The code described in this paper differs from its predecessors ( TVD method for shock capturing, enforcement of $`\mathbf{}𝒃=0`$ to machine precision ) in two main respects. First, we solve the induction equation in 2D advection-constraint steps without storing intermediate fluxes over the entire grid. Storage of the fluxes would require a $`3n^3`$ array, nearly half the memory used for the basic variables. Second, we use Jin and Xin’s (1995) method to implement TVD, which requires only $`30`$ floating point operations per grid cell per time step per variable. The benefits of these two methods are high resolution per grid cell, low operations count, and simplicity of coding (the public version of the 3D code is only 400 lines long.)

The price of separating the fluid and magnetic field updates into two steps, rather than updating all variables at once, is that the coupling between the velocity and magnetic field may be relatively weak compared to other methods. It was this point which led Stone & Norman (1992) to use an update along Alfven characteristics. However, our tests on both linear and nonlinear waves, as well as shocks, seem to indicate that in most circumstances the code performs rather well and no instabilities arise.

There are regimes in which we know the current code to be inaccurate or unstable. The generic setting is one where the characteristic families have very different velocities. These can occur at low $`\beta `$ for large amplitude Alfven waves, or highly supersonic flow with weak embedded shocks. For such regimes, customized modifications to the algorithm may be necessary (for the high mach number regime see, e.g. Trac and Pen 2003, in preparation; for the low $`\beta `$ case see Turner, Stone, Krolik, and Sano, in preparation). Also, since Jin and Xin’s (1995) method does not explicitly evolve the fluid variables along characteristics, our code may be more diffusive for low frequency waves when the ratio of fast and slow wave speeds is large (either large or small $`\beta `$.)

## 6 Conclusions

We have presented the algorithm and tests for a simple and robust MHD code which incorporates all features of modern high resolution shock capturing. It is second order accurate away from extrema, requires no memory overhead beyond storing the fluid variables, optimizes easily to many computer architectures, and offers simplicity in the coding. We have tested this code on linear and nonlinear MHD waves as well as shocks. The single processor code can be freely downloaded at

```
 http://www.cita.utoronto.ca/~pen/MHD

```

. We have also implemented a parallel version, which scales well on very large commodity beowulf clusters.

We thank Neal Turner for several useful discussions. Phil Arras is an NSF Astronomy and Astrophysics Postdoctoral Fellow. Computing resources were provided by the Canada Foundation for Innovation. SW is supported by the Taiwan NSC.

## Appendix A  advection in one dimension

Here we review Jin and Xin’s (1995) solution of the advection equation. We will focus on a scalar equation, but extension to a vector equation is straightforward. A more detailed discussion as well as code has recently been published by Trac & Pen (2003).

The advection equation for a quantity $`u`$ with flux $`f`$ is

$`_tu+_xf=0.`$ (A1)
Jin and Xin’s (1995) symmetric method is to define the new variable $`w=f/c`$ and equations for $`u`$ and $`w`$

$`_tu+_x(cw)`$ $`=`$ $`0`$
$`_tw+_x(cu)`$ $`=`$ $`0.`$ (A2)
These equations can be written in terms of left and right moving variables by defining $`u_r=(u+w)/2`$ and $`u_l=(uw)/2`$. These variables satisfy the equations

$`_tu_r+_x(cu_r)`$ $`=`$ $`0`$
$`_tu_l_x(cu_l)`$ $`=`$ $`0,`$ (A3)
which describe information propagating to the right and left, respectively.

To solve eq.A3 over a full time step with second order accuracy, we first advance $`u_r`$ and $`u_l`$ over a half time step using the first order upwind donor cell formula. These values are then used to construct a second order accurate upwind flux using any of the known nonlinear TVD limiters such as minmod, Van Leer, or superbee. Finally, given the updated values for $`u_r`$ and $`u_l`$, we reconstruct $`u=u_r+u_l`$.

For stability, the value of the flux freezing speed $`c`$ must be chosen larger than the speed at which information propagates. As discussed in the text, we set $`c=|v|`$ when advecting the magnetic field, and $`c=\mathrm{cfl}(|v|+(\gamma p/\rho +b^2/\rho )^{1/2})^1`$.

How can one relate TVD to the “artificial viscosity” schemes? These schemes add in a nonlinear viscosity term in order to prevent instabilities, as well as damp away oscillations which may occur near disconintuities. However, this viscosity tends to prevent the formation of discontinuities on scales of order one cell, severely degrading the resolution of the simulation. TVD may be viewed as a strongly nonlinear flux limiter which adds just enough diffusion to prevent numerical instabilities. TVD can often capture shocks in only one or two cells. Away from discontinuities, maxima or minima, TVD is second order in space, but at a maxima it is only first order.


text:
# A Turbulent Origin for Flocculent Spiral Structure in Galaxies: II. Observations and Models of M331footnote 11footnote 1Based on data obtained at the Canada-France-Hawaii Telescope

## 1 Introduction

The flocculent spiral structure in several nearby galaxies was recently shown to have a power spectrum for azimuthal scans that resembles the power spectrum of HI emission from the Large Magellanic Clouds (Elmegreen, Elmegreen, & Leitner 2003; hereafter Paper I). This power spectrum has a long-wavelength part that falls as $`1/k`$ for wavenumber $`k`$, an intermediate part that falls approximately as $`k^{5/3}`$, and a short wave part that falls as $`k^1`$ again. Individual stars contribute most to the short waves, and the brightest stars can dominate this part of the power spectrum if they are not removed. The distance scale for the $`k^{5/3}`$ part typically extends up to several hundred parsecs, depending on galaxy and galacto-centric radius. For NGC 5055, which is a large flocculent galaxy, the $`k^{5/3}`$ part extends up to 1 kpc.

Aside from the short-wave contamination from individual bright stars, these optical power spectra closely resemble the power spectrum of HI emission from the LMC (Elmegreen, Kim, & Staveley-Smith 2001) and dust absorption from the nuclear regions of two Sa-type galaxies (Elmegreen, Elmegreen & Eberwein 2002). This similarity led us to suggest that young stars follow the turbulent gas as they form, and that the outer scale for this turbulence is comparable to the disk thickness or the inverse Jeans length in the interstellar medium. The implication of this is not only that gravitational instabilities form flocculent arms, which was well-known before from numerical simulations (Sellwood & Carlberg 1984), but also that these instabilities generate much of the turbulence in the interstellar medium, which both structures the gas and causes the stars to form in fractal patterns (see also Huber & Pfenniger 2001; Wada, Meurer, & Norman 2002).

Paper I reviewed the observations of power spectra in local gas and dust (e.g., Crovisier & Dickey 1983; Gautier, et al. 1992; Green 1993; Armstrong et al. 1995; Stützki et al. 1998; Schlegel, Finkbeiner, & Davis 1998; Deshpande, Dwarakanath, & Goss 2000; Dickey et al. 2001) and in whole galaxies (Stanimirovic, et al. 1999; Stanimirovic et al. 2000; Elmegreen, Kim, & Staveley-Smith 2001). It also discussed the possible links between gaseous density structures and the structures that come from velocity in a turbulent medium (e.g., Lazarian & Pogosyan 2000; Goldman 2000; Stanimirovic & Lazarian 2001; Lazarian et al. 2001; Lithwick & Goldreich 2001; Cho & Lazarian 2003). Here we model the stellar light distribution in two galaxies to illustrate that star formation follows this turbulent gas and acquires a similar power spectrum. The nature of interstellar turbulence and star formation are not understood well enough to explain the results. Several processes are possible, including the formation of clouds and stellar complexes in moving sub-clouds that act like passive scalars in a large-scale turbulent flow (Goldman 2000; Boldyrev, Nordlund, & Padoan 2002; Paper I), the formation of clouds and stars in turbulence-compressed regions (Klessen, Heitsch, & Mac Low 2000; Ossenkopf, Klessen, & Heitsch 2001; Klessen 2001; Padoan et al. 2001a; Padoan & Nordlund 2002) and the formation of stars in shocks that are driven by other stars in a turbulent medium (Elmegreen 2002a). All of these processes give about the same power law structure for density and star formation.

The advent of large-scale digital CCD images measuring several thousand pixels on a side has made power spectrum analyses of galactic structure meaningful. Large images are necessary to get enough dynamic range to see the power law part of a power spectrum if there is one. In this respect, the optical survey of galaxies by the CFH12K camera of the Canada-France-Hawaii telescope (CFHT) is ideal. Here we analyze the CFHT image of M33 (Cuillandre, Lequeux, & Loinard 1999), which has an original pixel size of 0.206 arcsec and a binned size of 0.412 arcsec for the full galaxy mosaic. The maximum azimuthal circumference for the part we use is 10000 px, or 4120 arcsec, which gives scales covering 4 orders of magnitude. We take power spectra in both the azimuthal direction at equally spaced radii and along the spiral direction, following the arm and interarm regions.

We also model this power spectrum and the power spectrum of NGC 5055, a flocculent galaxy from Paper I, with fake-galaxy images of a fractal Brownian motion continuum plus discrete stars that follow this continuum.

## 2 Power Spectra of the Galaxy M33

Figure 1 shows two images of M33 with azimuthal and spiral lines at which scans of intensity were obtained and converted to 1D power spectra. The azimuthal scans are circles in the galaxy disk plane and the spiral scans go through both arm and interarm regions at a constant pitch angle (24.5 degrees). This pitch angle fits the strong arm in the south using deprojected coordinates (see Sandage and Humphreys 1980 for pitch angle fits to each arm). Power spectra were obtained for intensity scans along these directions, with each scan one pixel wide to maximize k-space resolution. For the azimuthal scans, 8 adjacent pixel-wide power spectra were averaged together to make the final power spectrum at each radius. These 8 adjacent scans for the 10 selected radii are shown in the figure as dark pixels. For the spiral scans, 21 adjacent pixel scans were used separately for each power spectrum, and then these 21 power spectra were averaged together to give the resultant power spectrum along each of the 8 spiral or interspiral cuts.

Figure 2 has the azimuthal intensity scan from the ellipse that is fourth out from the center of the R-band image, and it shows the power spectrum of this scan multiplied by $`k^{5/3}`$. The original scan is on the left and a version with the brightest stars removed is on the right. The power spectrum is the sum of the squares of the sine and cosine Fourier transforms, plotted in log-log as a function of the wavenumber, $`k`$. The abscissa in the figure is the wave number normalized so that the right-most point, with a value of 1, has $`k=1/2`$ pixels<sup>-1</sup>. The distance corresponding to any normalized $`k`$ value is $`2/k`$ pixels. For our images one pixel is 1.68 parsecs, assuming a distance to M33 of 0.84 Mpc (Freedman, Wilson & Madore 1991).

The power spectra in Figure 2 and other figures here have been multiplied by $`k^{5/3}`$ in order to flatten what is normally a $`k^{5/3}`$ spectrum for incompressible Kolmogorov turbulence and make this part easier to see. There are many possible explanations for a power spectrum that is approximately $`k^{5/3}`$, as reviewed in detail in Paper I. It is probably not the result of motions in a perfectly incompressible fluid, as for the Kolmogorov problem, but some degree of incompressibility is still possible (Goldman 2000; Boldyrev, Nordlund, & Padoan 2002; Lithwick & Goldreich 2001). Even if it is highly compressible, the power spectrum of density structure is about the same (see Lazarian & Pogosyan 2000; Cho & Lazarian 2003).

Bright stars are clipped from the intensity scans in 3 steps. First we find the average intensity level of a clipped version of the scan, clipped below zero and at a high enough level to remove the brightest and saturated stars, all of which are foreground. Then we find the running boxcar average 31 pixels wide of the same scan, clipped a second time at some intermediate height above the average to remove the pedestals of the brightest stars, but not clipped low enough to remove the stars in M33. Finally, we clip the scan at a level above the running average that is 3 times the Gaussian $`\sigma `$ for the rms noise in this final clipped version; this last step is done with several iterations. The first step is necessary to account for the exponential disk; the second to include large-scale variations from star-formation regions and spiral waves, and the third to reduce the impact of stars in the power spectrum.

The power spectrum has a slightly steeper rise in Figure 2 when bright stars are included, and it has a deeper dip at high frequency. M33 has many foreground stars and bright point sources inside the galaxy so star removal is difficult. Thus there is a residual high-$`k`$ dip in all of the results shown here. The rising part could also be a bit too steep as a result of residual stars. This translates into a non-normalized power spectrum that is too shallow in its decline. The galaxies in Paper I had fewer foreground stars and shallower rising parts in their $`k^{5/3}`$normalized power spectra. One of them, NGC 5055, will be modelled in the next section.

Power spectra of M33 for ten radii in three passbands are shown in Figure 3, and power spectra for the spiral scans in R band, along with the intensity profiles averaged over 21 adjacent single-pixel wide scans (as discussed above), are shown in Figure 4. The power spectra are all similar regardless of color, radius, arm, or interarm position. After multiplication by $`k^{5/3}`$ as in these figures, there is a slowly rising part on the left, a flat part in the middle-right, and a dip and second rising part on the far right. The density wave in M33 appears in the left-most few points in the azimuthal power spectra (Fig. 3), which turn up by a factor of $`10`$ compared to the extrapolated curves at low $`k`$. The left-most point corresponds to a one-arm spiral, and the second-left most point corresponds to a two-arm spiral. There is no such systematic turn up in the spiral arm and interarm power spectrum (Fig 4).

The similarity of the power spectra for the different passbands may be surprising considering the galaxy looks smoother in red colors. However the power spectrum only gives a relative measure of power at each spatial frequency, whereas the visual impression of smoothness is based mostly on the ratio of high frequency emission to average brightness. Slight differences in the 3 passbands are visible at high $`k`$, where the prominence of the stellar dip decreases toward the red.

The power spectra and intensity scans of a clipped H$`\alpha `$ image of M33 taken with the same instrument are shown in Figure 5. The dashed line has a slope of 1, which is similar to the slope of the $`k^{5/3}`$-normalized power spectrum of the optical emission in Figure 3. There are many small H$`\alpha `$ sources that act like stars in the intensity scans, giving dips in the power spectra at high frequency as for the optical images.

The top axes in Figures 3, 4, and 5 give the distance scale in parsecs that corresponds to $`2/k`$ on the abscissa. The outer scale for the horizontal part of the power spectrum is $`2/0.02100`$ pixels, which corresponds to 170 parsecs. This scale is about the same as what we found for galaxies in Paper I, although some aspects of the overall shape result from pixelation, as shown in the next section, and are therefore distance-dependent.

## 3 Model Power Spectra

Fractal Brownian motion models were made by filling half of a $`640\times 640\times 640`$ cube in $`k`$-space with random complex numbers having real and imaginary values between 0 and 1, and then multiplying these values by $`k^{\beta /2}`$ for $`k=\left(k_x^2+k_y^2+k_z^2\right)^{1/2}`$. The cube is only half full because of the symmetries in the Fourier transforms. The inverse Fourier transform of this cube, from complex to real numbers, gives a three-dimensional fractal with positive and negative numbers having a Gaussian distribution of values. This fractal is then exponentiated to give another fractal, now with a log-normal distribution of all-positive values. This last step is done to mimic models of turbulence which produce a log-normal distribution of density (e.g., Wada & Norman 2001). An example of a fractal cube made this way is in Figure 2 of Elmegreen (2002b). Three dimensional power spectra of the model reproduce the input power spectrum slope $`\beta `$ (which is twice the input power of $`k`$ because the power spectrum is made from the square of the Fourier transform), and one-dimensional power spectra of the model have a slope that is shallower by 2. We denote the slope of the 1D model power spectrum by $`\alpha `$ and confirmed experimentally that $`\alpha =\beta 2`$ for $`\beta `$ greater than 2. Input $`\beta <2`$ give power spectra indistinguishable from noise ($`\alpha =0`$). Other fractal Brownian motion models of interstellar structure and a more detailed discussion of $`\beta `$ are in Stützki et al. (1998).

We varied the intrinsic power $`\beta /2`$ from $`0.33`$ to $`+2.66`$ and plotted the power spectrum for each model using density scans along a strip and along a circle in the midplane. We also integrated the cube through one dimension with a Gaussian weighting factor centered on the midplane, to simulate the average Gaussian profile of the interstellar medium viewed through a disk face-on. These power spectra should represent the observations of a fractal gas, and in fact models of the HI emission from the LMC were fitted in this way (Elmegreen, Kim & Staveley-Smith 2001).

The optical power spectra presented here are for starlight, and models of this require an extra step. To simulate stars, we use small points of “emission” with Gaussian profiles having a $`\sigma `$ of 2 pixels and an intensity equal to a random number between 0 and 30; these values were chosen to match the width and depth of the high-$`k`$ dip. We consider three cases: (1) an image made only of stars with a random distribution in space; (2) an image made only of stars with a fractal distribution in space, and (3) an image made of the fractal continuum in addition to stars with the same fractal distribution in space. The fractal distribution of stellar points is done by going to each cell in the fractal model cube and choosing a random number between 0 and 1. If the density in the cube at this point (differenced from the minimum density and divided by the maximum density for normalization) exceeds some factor times the random number, then we place a star there with the assumed Gaussian profile. If the density is less than this value, then we do not place a star there. If adjacent points have stars, then we add all the stellar profiles together. The multiplicative factor for the (0,1) random number determines how many stars there are. A value of 0.4 was chosen to give intensity profiles that looked reasonably dense in resolved stellar sources. In case (3), we chose stars in the same way but added this starlight to the fractal continuum in each cell.

The most reasonable case is (3). In a typical galaxy there will be both isolated bright stars and a haze of blended faint stars. Case (3) assumes that both follow the fractal distribution of the gas as a result of moderately recent star formation. The underlying smooth distribution of old stars will not contribute much to the high spatial frequencies but will appear mostly in the lowest few $`k`$ values, like the spiral density wave. These can be ignored in both the modelled and the observed power spectra.

Figures 6, 7, and 8 show the results for three values of $`\alpha `$ which bracket our observations. For the best-fit $`\alpha =0.66`$ case in Figure 7, the power spectrum and intensity scan from the 4th ellipse out in the R band image is shown at the top. The left panel shows power spectra multiplied by $`k^{5/3}`$, the middle panel shows conventional power spectra with no multiplier, and the right panel shows the corresponding intensity scans. The power spectrum below the R band observation is for a strip from left to right in the fractal cube. Below that is the power spectrum from a circle at the maximum radius of the model cube, giving a circumference 1700 pixels long. Next comes the power spectrum from continuum plus stars in a fractal pattern, and below that is the power spectrum from the pure stellar image with a fractal pattern. At the bottom is the spectrum for a pure-star image with a random distribution of stellar positions. This pure-star power spectrum is the same in Figures 6, 7, and 8 because it does not depend on $`\alpha `$.

The best match between the star+continuum power spectrum and the observation of M33 is for $`\alpha 0.66`$. The $`\alpha =0`$ case has a power spectrum that is slightly too steep in its rising part on the left (meaning that in a conventional plot, not normalized with multiplication by $`k^{5/3}`$, the power spectrum is too shallow in its decline here). The $`\alpha =1.33`$ case has a power spectrum that is too slowly rising in this part (too steeply falling in a conventional plot). Generally, the star+continuum case is a better match than the fractal-star case, which has too large a dip at high frequency from the intrinsic stellar profiles. The random star case is a poor fit, having a power spectrum in the normalized plot with an increase as $`k^{5/3}`$, which corresponds to a flat power spectrum in a conventional plot (as seen at the bottom of the middle panel of Fig. 7).

The azimuthal power spectrum from the fifth ellipse out in the I-band image of NGC 5055 (Paper I; data from HST archives) is reproduced in Figure 9, along with several models. This galaxy has fewer foreground stars and so the dip at high frequency is very slight. The two models at the bottom are pure stars in a fractal pattern. The bottom one has a star brightness that is a random number multiplied by 5 and the next one up has a star brightness that is a random number multiplied by 20. Both star models have Gaussian $`\sigma =2`$ px. The power spectra are identical although the intensity traces on the right show these differences clearly. Next up from the bottom are three models with a fractal distribution of stars and the associated continuum. They each have a star brightness equal to a random number times 5 and the stars have $`\sigma =2`$ px again. This is only one-sixth the star brightness used for the M33 models, and the high-frequency dip is smaller too, in proportion. The three models differ in the slope of the intrinsic power spectrum, decreasing from $`\alpha =1.70`$ to $`1.50`$ to $`1.33`$ with increasing height in the figure. The power spectra show this difference at low frequency where the rising part gets slightly steeper as $`\alpha `$ decreases. The middle one, with $`\alpha =1.50`$, is the best fit to NGC 5055. The power spectrum of the pure continuum with no stars in this case is shown next to the observation of NGC 5055. It differs from the star+continuum case only slightly at high frequency. We conclude from these models that NGC 5055 has relatively faint point sources and an intrinsic power spectrum for the optical emission with a slope of $`1.5\pm 0.2`$.

The $`k^{1.5\pm 0.2}`$ power spectrum fit to NGC 5055 is clearly steeper than the $`k^{0.66\pm 0.66}`$ fit to M33. The difference could be the result of residual stars and unresolved clusters in the M33 galaxy, although NGC 5055, which is 9 times more distant, has 4 times better resolution from HST so the pixel scale in parsecs (3.5 pc) is only a factor of 2 greater for NGC 5055 than for M33. Residual foreground stars could be more important. M33 has many more foreground stars than NGC 5055 because of its large angular size. Foreground stars are randomly positioned and their power spectrum varies like $`k^0`$, as in Figures 6, 7, and 8. Contamination from numerous random stars will flatten the spectrum.

Figure 10 shows a fit to M33 that starts with the underlying galaxy model for NGC 5055 and adds random foreground stars with an amplitude 6 times larger than in the NGC 5055 model itself and a stellar width of $`\sigma =2`$ px. At the bottom is the power spectrum and intensity scan for the foreground stars. Next is the observation of NGC 5055, from Figure 9. The third curves up are the power spectrum and intensity scan for the sum of the best-fit model of NGC 5055 (the $`\alpha =1.5`$ model in Fig. 9) and the foreground stars. The top curves are the observations of M33. The model of M33 with foreground stars has the same intrinsic galactic properties as NGC 5055, namely a superposition of galactic stars and a continuum that both follow a fractal pattern giving a 1D power spectrum slope of $`1.5`$. The foreground stars make the slope appear shallower, $`0.66`$, like the observed slope in M33.

## 4 Conclusions

The power spectra of optical emission in galaxies show a characteristic structure that can be modelled by a superposition of a fractal unresolved brightness distribution and a corresponding fractal distribution of point-like sources. The intrinsic power spectrum of these distributions has a slope of $`0.66\pm 0.66`$ for M33 and a slope of $`1.5\pm 0.2`$ for NGC 5055. The relatively flat slope for M33 is probably the result of foreground stars, which are difficult to remove from the observations. A model of M33 with foreground stars superposed on the NGC 5055 model demonstrates this point. The power spectrum slope for the foreground-corrected version of M33, and the slopes for NGC 5055 and the other galaxies in Paper I are all comparable to the one-dimensional power spectrum from turbulence. We infer from this that young stars and clusters are distributed in a fractal pattern that is generated by a turbulent gas.

Acknowledgements: B.G.E is supported by NSF grant AST-0205097.