text:
# The Electron Spectrum in 3C 279 and the Observed Emission Spectrum

## 1 Introduction

Blazars are active galactic nuclei characterized by high energy emission and short time variability (e.g., Krolik, 1999). Recent observations of blazars show that they are powerful sources of high energy emission (Ulrich, Maraschi, & Urry, 1997, for review), which is explained by nonthermal emission by particles in relativistic jets. The emission of blazars in the GeV energy range was observed by EGRET onboard the Compton Gamma Ray Observatory (Thompson et al., 1995; Mukherjee et al., 1997; Hartman et al., 2001). These observations show that a large fraction of emission power from blazars is in $`\gamma `$-ray band (see Ghisellini et al. 1998 and Collmar 2001 for reviews).

Multiwavelength observations of 3C 279 (redshift $`z=0.538`$) were performed several times (Maraschi et al., 1994; Wehrle et al., 1998; Hartman et al., 2001). Based on these results, $`\gamma `$-rays from 3C 279 are found to be explained by inverse Compton scattering of external UV radiation (e.g., Inoue & Takahara, 1996). The sources of seed soft photons (e.g., accretion disks (Dermer & Schlickeiser, 1993), broad-line regions (Sikora, Begelman, & Rees, 1994), and beamed radiation reflected by broad-line region clouds (Ghisellini & Madau, 1996)) have been proposed by many authors. Recently, Hartman et al. (2001) assumed soft photons from both the accretion disk and the broad-line region in order to model the emission spectra from 3C 279. Ballo et al. (2002), on the other hand, applied a broad-line region model to 3C 279.

Most emission models have assumed that electrons obey a power law (e.g., Hartman et al., 2001) or a broken power law (e.g., Inoue & Takahara, 1996; Ballo et al., 2002). However, an alternative electron spectrum is possible as a result of the Klein-Nishina (K-N) effects; the smaller cooling rate due to the K-N effects makes the electron spectrum harder (Dermer & Atoyan, 2002), and the emissivity of scattered photons is smaller (Georganopoulos, Kirk, & Mastichiadis, 2001). Thus, the electron spectrum must be solved self-consistently, including the cooling rate in the K-N regime.

In previous work, we solved the kinetic equations of electrons and photons simultaneously in a relativistically moving blob (Li & Kusunose 2000; Kino, Takahara, & Kusunose 2002, hereafter KTK). Here we include the inverse Compton scattering of external soft photons (external Compton scattering, hereafter EC) and obtain the physical quantities for 3C 279 such as the bulk Lorentz factor, the energy densities of electrons and magnetic fields, etc.

The variability of blazars on small timescales is often observed (e.g., Wehrle et al., 1998). The time variation of fluxes is not only useful for obtaining constraints on the size of the emission region but is also important for studying the emission and particle acceleration mechanisms (Mastichiadis & Kirk, 1997; Kusunose, Takahara, & Li, 2000; Li & Kusunose, 2000; Böttcher and Chiang, 2002). Flares observed in X-ray and $`\gamma `$-ray regimes might be explained by the internal collisions of plasma blobs in a jet and the subsequent cooling. Without considering the details of particle acceleration after the collision, it is possible to obtain the emission spectrum by assuming that the time evolution occurs quasi-steadily. In this Letter we present a model to explain the flare of 3C 279 observed in 1996 February.

In §2, kinetic equations for electrons and photons are described. Numerical results are given in §3 for a steady state. Results for time variability are shown in §4. Finally, §5 is devoted to the summary of our results and discussion.

## 2 Model

The formulation used in this Letter is the same as in Li & Kusunose (2000) and KTK. We assume that a spherical blob with radius $`R`$ moves with a relativistic speed with Lorentz factor $`\mathrm{\Gamma }`$. Accelerated electrons obeying a power law are injected into the blob uniformly, and they are cooled by synchrotron emission and inverse Compton scattering.

### 2.1 Kinetic Equations

The equation describing the time evolution of the electron number spectrum in the blob is given by

$$\frac{n_e(\gamma )}{t}=\frac{}{\gamma }\left[\left(\frac{d\gamma }{dt}\right)_{\mathrm{loss}}n_e(\gamma )\right]\frac{n_e(\gamma )}{t_{e,\mathrm{esc}}}+q(\gamma ),$$
(1)
where $`\gamma `$ is the Lorentz factor of electrons and $`n_e(\gamma )`$ is the number density of electrons per unit $`\gamma `$. We assume that nonthermal electrons with a power law are injected at a rate:

$$q(\gamma )=q_0\gamma ^p\mathrm{exp}(\gamma /\gamma _{\mathrm{max}})\text{for}\gamma \gamma _{\mathrm{min}},$$
(2)
where $`q_0`$ is the normalization. The injection is assumed to continue during a simulation with the rate given above. The energy-loss rate of electrons is denoted by $`(d\gamma /dt)_{\mathrm{loss}}`$, which is due to synchrotron radiation and inverse Compton scattering. The synchrotron-loss rate is calculated by the formulation given by Robinson & Melrose (1984) for mildly relativistic electrons and by Crusius & Schlickeiser (1986) for relativistic electrons. The cooling rate of inverse Compton scattering is calculated by Coppi & Blandford (1990); although the EC process is anisotropic, we use the formulation for isotropic scattering. The value of the electron-escape timescale, $`t_{e,\mathrm{esc}}`$, is set to be $`3R/c`$, assuming that the escape of electrons from the jet is likely due to advection. We use a method developed by Xiao, Yabe, & Ito (1996) for numerical calculations.

To obtain the isotropic photon field in the blob frame, i.e., except the EC component, we solve the following equation:

$$\frac{n_{\mathrm{ph}}(\epsilon )}{t}=\dot{n}_\mathrm{C}(\epsilon )+\dot{n}_{\mathrm{em}}(\epsilon )\dot{n}_{\mathrm{abs}}(\epsilon )\frac{n_{\mathrm{ph}}(\epsilon )}{t_{\gamma ,\mathrm{esc}}},$$
(3)
where $`n_{\mathrm{ph}}(\epsilon )`$ is the photon number density per unit energy $`\epsilon `$; Compton scattering and synchrotron emission are denoted by $`\dot{n}_\mathrm{C}(\epsilon )`$ and $`\dot{n}_{\mathrm{em}}(\epsilon )`$, respectively; $`\dot{n}_{\mathrm{abs}}(\epsilon )`$ is absorption due to the self-absorption of synchrotron process and pair production; and $`n_{\mathrm{ph}}(\epsilon )/t_{\gamma ,\mathrm{esc}}`$ denotes the escape of photons (we set $`t_{\gamma ,\mathrm{esc}}=R/c`$). The EC component is separately calculated following Georganopoulos et al. (2001; see also Dermer & Schlickeiser 2002) because the external photons are highly anisotropic in the blob frame. We assume that in the rest frame of the central black hole, the external radiation is isotropic with energy density $`u_{\mathrm{ext}}`$ and monochromatic with energy $`\epsilon _{\mathrm{ext}}`$. In the blob frame, the comoving energy density of the soft photons is $`u_{\mathrm{ext}}^{}=u_{\mathrm{ext}}\mathrm{\Gamma }^2(1+\beta _\mathrm{\Gamma }^2/3)`$, where $`\beta _\mathrm{\Gamma }=(1\mathrm{\Gamma }^2)^{1/2}`$ (Dermer & Schlickeiser, 2002). On the other hand, the energy of the external photons is $`\epsilon _{\mathrm{ext}}^{}\mathrm{\Gamma }\epsilon _{\mathrm{ext}}`$.

The observed radiation is boosted by the Doppler effect characterized by the beaming factor $`𝒟=[\mathrm{\Gamma }(1\beta _\mathrm{\Gamma }\mathrm{cos}\theta _{\mathrm{obs}})]^1`$, where $`\theta _{\mathrm{obs}}`$ is the angle between the jet propagation and the line of sight (Blandford & Königl, 1979). In this Letter, we assume $`𝒟=\mathrm{\Gamma }`$. The parameters to be determined are $`\mathrm{\Gamma }`$, $`\gamma _{\mathrm{min}}`$, $`\gamma _{\mathrm{max}}`$, $`p`$, $`q_0`$, $`B`$, $`R`$, $`\epsilon _{\mathrm{ext}}`$, and $`u_{\mathrm{ext}}`$.

The comoving quantities are transformed back into the observer’s frame depending on the beaming factor and the redshift. The observed photon energy $`\epsilon _{\mathrm{obs}}=\epsilon 𝒟/(1+z)`$, and the time duration $`dt_{\mathrm{obs}}=dt(1+z)/𝒟`$.

## 3 Numerical Model of a Steady State

We use the observational data in Hartman et al. (2001). In particular, the spectral energy distribution (SED) observed between 1996 January 16 and 30 (P5a in their paper) is used for our model of a steady state. 3C 279 is in a preflare stage in this period and is not necessarily in a quiescent sate. However, we fitted the data assuming that the spectrum is approximated by a steady state solution of the kinetic equations.

In Figure 1, our model is plotted with the observed data. We obtain the following values for the parameters: $`\mathrm{\Gamma }=25`$, $`B=0.3`$ G, $`R=7\times 10^{16}`$ cm, $`p=1.8`$, $`\gamma _{\mathrm{min}}=4`$, $`\gamma _{\mathrm{max}}=1.8\times 10^3`$, the injection rate $`2.3\times 10^4\text{cm}^3\text{s}^1`$, $`\epsilon _{\mathrm{ext}}=50`$ eV, and $`u_{\mathrm{ext}}=2\times 10^4\text{ergs}\text{cm}^3`$. The spectrum below $`10^{15}`$ Hz is produced by synchrotron radiation, although the radio emission is probably from the outer regions of the jet. The radiation between $`10^{16}`$ and $`10^{19}`$ Hz is from synchrotron-self-Compton scattering, and the $`\gamma `$-rays above $`10^{20}`$ Hz are from the inverse Compton scattering of external soft photons. As shown in Figure 1 by the dashed line, the photons scattered by electrons with $`\gamma _{\mathrm{min}}`$ appear around $`10^{20}`$ Hz. Thus the flux around $`10^{20}`$ Hz is sensitive to the value of $`\gamma _{\mathrm{min}}`$. However, because the data in this regime are mostly upper limits, the constraint on $`\gamma _{\mathrm{min}}`$ is weak.

Based on our numerical calculations, it is found that $`u_{\mathrm{kin}}/u_B=36`$, where $`u_{\mathrm{kin}}`$ is the energy density of nonthermal electrons and $`u_B`$ is the magnetic energy density. We also find that the Poynting power is given by $`L_{\mathrm{Poy}}=\pi R^2cu_B𝒟^2=1.03\times 10^{45}\text{ergs}\text{s}^1`$ and that the kinetic power is given by $`L_{\mathrm{kin}}=\pi R^2cu_{\mathrm{kin}}𝒟^2=3.69\times 10^{46}\text{ergs}\text{s}^1`$. The dominance of the kinetic energy of nonthermal electrons over that of magnetic fields is the same as that for TeV blazars such as Mrk 421 and Mrk 501 (KTK).

In Figure 2, the electron spectrum for the SED shown in Figure 1 is presented. It is found that the electron spectrum is rather flat in the $`\gamma `$-$`\gamma ^2n_e(\gamma )`$ plot. This is because the cooling of high-energy electrons is inefficient, owing to the K-N effects (Blumenthal, 1971). Note that the external photon energy in the blob frame is $`\epsilon ^{}\mathrm{\Gamma }\epsilon _{\mathrm{ext}}`$ and that the K-N effects become effective for electrons with $`\gamma _{\mathrm{K}\mathrm{N}}m_ec^2/\epsilon ^{}400`$, which is less than $`\gamma _{\mathrm{max}}`$. The spectral shape is contrary to the models often employed to fit the observed data of blazars, in which $`n_e(\gamma )\gamma ^p`$ below a break energy $`\gamma _{\mathrm{br}}`$ and $`n_e(\gamma )\gamma ^{(p+1)}`$ above $`\gamma _{\mathrm{br}}`$; the value of $`\gamma _{\mathrm{br}}`$ is determined by $`t_{\mathrm{cool}}=t_{e,\mathrm{esc}}`$, where $`t_{\mathrm{cool}}`$ is the electron cooling time. Nominally, $`\gamma _{\mathrm{br}}`$ is about 30. As seen in Figure 2, the actual spectral shape is much different from this conventional expectation. The flat $`\gamma `$-ray spectrum in $`10^{21}10^{24}`$ Hz in the $`\nu `$-$`F_\nu `$ plot is the result of these K-N effects.

The optical depth for photon absorption against electron-positron pair production is calculated (Gould & Schréder, 1967). The optical depth is much smaller than unity and $`\gamma `$-rays escape without absorption.

## 4 Flare

The time variability of emission spectra from blazars is quite common. 3C 279 also shows time variation from radio to $`\gamma `$-rays (e.g., Wehrle et al., 1998). To understand the jet formation mechanisms and particle acceleration processes, it is important to understand the mechanism that causes the time variability. Sikora et al. (2001) and Spada et al. (2001) show how flare light curves behave according to the shock-in-jet models, which assume that blobs with different speeds collide and that the kinetic energy of the bulk motion is dissipated. Ballo et al. (2002) interpreted the time variability of 3C 279 as being the result of various values of $`\mathrm{\Gamma }`$, although they assumed a broken power-law spectrum of electrons without solving the kinetic equation.

We apply the internal shock model, by which Takahara et al. (2003) will explain the time variability of Mrk 421, to a flare observed from 3C 279. Without considering the details of the collision process, we simply assume that a coalesced blob has a bulk Lorentz factor $`\mathrm{\Gamma }_f`$ and that various physical quantities scale with $`\mathrm{\Gamma }`$. We examine whether or not the flares are explained by a change in $`\mathrm{\Gamma }`$, keeping a constant opening angle for the jet cone. The distance of the collision location from the base of the jet, $`d`$, scales proportional to $`\mathrm{\Gamma }^2`$, while the lateral size of the observed region behaves as $`R\mathrm{\Gamma }`$, which is roughly equal to the shell thickness in the comoving frame, because the angle from the line of sight is limited by $`\mathrm{\Gamma }^1`$. The soft photon energy density changes as $`u_{\mathrm{ext}}\mathrm{\Gamma }^2d^2\mathrm{\Gamma }^2`$. We assume that the kinetic and Poynting powers scale proportional to $`\mathrm{\Gamma }^2`$ so that the total number flux remains constant. Thus, the injection rate in the comoving frame, $`q`$, scales proportional to $`d^2R^1\mathrm{\Gamma }^5`$, and the magnetic field scales proportional to $`\mathrm{\Gamma }^2`$. We also assume $`\gamma _{\mathrm{min}}`$ and $`\gamma _{\mathrm{max}}`$ do not depend on $`\mathrm{\Gamma }`$. The value of $`\epsilon _{\mathrm{ext}}^{}`$ scales proportional to $`\mathrm{\Gamma }`$ because of the Doppler effects. Using this scaling, we calculate a flare spectrum. Assuming that the spectrum shown in Figure 1 is in a quiescent state, the spectrum observed between 1996 January 30 and February 6 is fitted. Here the bulk Lorentz factor in Figure 1 is denoted as $`\mathrm{\Gamma }_q`$, and that of a shocked blob is given by $`\mathrm{\Gamma }_f`$. When $`\mathrm{\Gamma }_f=1.6\mathrm{\Gamma }_q=40`$, the $`\gamma `$-ray spectrum is fitted fairly well by our model (Figure 3), although there are some minor discrepancies: The slope in the optical-UV regime is slightly different; the flux between $`10^{19}`$ and $`10^{20}`$ Hz is underestimated; and the flux between $`10^{20}`$ and $`10^{21}`$ Hz is overestimated.

## 5 Summary and Discussion

We used the kinetic equations of photons and electrons to calculate emission spectra from 3C 279, assuming that the emission is from a relativistically moving blob almost along the line of sight. Our numerical solution shows that the $`\gamma `$-rays are produced by inverse Compton scattering of external soft photons, which is consistent with previous work (e.g., Inoue & Takahara, 1996; Hartman et al., 2001; Ballo et al., 2002). However, we found that the electron spectrum is not a broken power law with a steeper spectrum above a break energy. Although electrons are cooled before escape (except those in the lowest energy range), the cooling efficiency for high-energy electrons decreases because of the K-N effects. As a result, a flat GeV $`\gamma `$-ray spectrum is obtained in the $`\nu `$-$`\nu F_\nu `$ plot. It should be noted that there is an alternative model of GeV emission that assumes different soft photon sources, i.e., the accretion disk and broad line regions (Hartman et al., 2001). In their model, because of the different temperatures of soft photons, a broad GeV emission is formed by inverse Compton scattering.

Our result of $`\epsilon _{\mathrm{ext}}=50`$ eV seems to be a little high compared with the conventional value of around 10 eV. However, the spectra of ionizing radiation that form the broad line clouds are not directly constrained observationally but should extend to energy high enough to multiply ionize atoms of heavy elements. From a theoretical point of view, even the emission from the accretion disk is known to deviate from the blackbody because of scattering effects. Thus, our choice may be suitable. Recently, it was suggested that IR radiation from dust might be important, if the energy dissipation occurs mainly far away from the central region (Błażejowski et al., 2000; Sikora et al., 2002); although we have not examined this case, the K-N effects will be weaker for this component.

The kinetic energy density of nonthermal electrons in the blob is an order of magnitude larger than the magnetic energy density, which is the same result as for TeV blazars (KTK). As a result, the particles transfer more energy in the jet from the central region to the outer region than the magnetic fields do.

We demonstrated that the flare observed between 1996 January 30 and February 6 is explained fairly well by the internal shock model with simple scaling laws described by a change in the bulk Lorentz factor. It is remarkable that the flare is fitted with those simple assumptions.

It was pointed out by Sikora & Madejski (2000) that the bulk motion of cold electrons/positrons in the jet scatters the external photons, resulting in observed emission peaking at energy $`\epsilon _{\mathrm{BC}}\mathrm{\Gamma }^2\epsilon _{\mathrm{ext}}`$. Our model with $`\mathrm{\Gamma }=25`$ and $`\epsilon _{\mathrm{ext}}=50`$ eV implies $`\epsilon _{\mathrm{BC}}31`$ keV ($`\nu _{\mathrm{BC}}7.5\times 10^{18}`$ Hz). The observed luminosity of bulk Compton (BC) scattering is estimated as

$$\begin{array}{cc}\hfill L_{\mathrm{BC}}& \mathrm{\Gamma }^2\left(\frac{4}{3}c\sigma _\mathrm{T}u_{\mathrm{ext}}\mathrm{\Gamma }^2\right)n_c𝑑V\hfill \\ & 3.0\times 10^{39}\left(\frac{\mathrm{\Gamma }}{25}\right)^4\frac{u_{\mathrm{ext}}}{2\times 10^4\text{ergs}\text{cm}^3}\left(\frac{R}{7\times 10^{16}\text{cm}}\right)^3n_c\text{ergs}\text{s}^1,\hfill \end{array}$$
(4)
where $`\sigma _\mathrm{T}`$ is the Thomson cross section, $`n_c`$ is the number density of cold electrons/positrons, and the integration is done over the blob volume. For 3C 279 with luminosity distance $`d_L=2.4\times 10^3`$ Mpc, where the Hubble constant $`H_0=75\text{km}\text{s}^1\text{Mpc}^1`$ is assumed, $`\nu F_\nu `$ at $`\epsilon _{\mathrm{BC}}`$ is $`L_{\mathrm{BC}}/(4\pi d_L^2)7\times 10^{14}n_c/n_e`$ $`\text{ergs }\text{cm}^2\text{s}^1`$, where $`n_e=1.65\times 10^4\text{cm}^3`$ is the number density of nonthermal electrons shown in Figure 2. The value of $`n_c/n_e`$ is unknown, and unless it is about 100 or more, the BC spectrum is not observable; note that the observed value is $`10^{11}\text{ergs}\text{cm}^2\text{s}^1`$ at $`\epsilon _{\mathrm{BC}}31`$ keV. Regions between the base of the jet and the blob may contribute to BC scattering because more cold electrons may exist there than in the blob. Let us assume that the bulk Lorentz factor is constant and that the number density of the cold electrons behaves as $`n_c\mathrm{}^2`$ for $`\mathrm{}>\mathrm{}_0`$, where $`\mathrm{}`$ is the distance from the central region and $`\mathrm{}_0`$ is the critical distance where the bulk Lorentz factor saturates or the jet becomes optically thin. If we assume that $`u_{\mathrm{ext}}`$ scales proportional to $`\mathrm{}^2`$, the predicted value of $`L_{\mathrm{BC}}`$ is $`d/\mathrm{}_0`$ times larger than that of equation (4), noting that $`n_cdV`$ is measured in the comoving frame of the blob. Even if there are no cold electrons is the blob, in the inner-jet regions there should exist a corresponding number of cold electrons to relativistic electrons in the blob. Thus, we obtain the constraint of $`\mathrm{}_0>0.01d`$. If $`n_c/n_e`$ in the blob is more than 100, too much BC emission is predicted, which probably means that $`u_{\mathrm{ext}}`$ does not increase much with a decrease in $`\mathrm{}`$. In this case some modifications of our model are needed.

This work has been partially supported by Scientific Research Grants (M.K. and F.T.: 13440061; M.K.:15037210; F.T.: 14079205 and 14340066) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.


format: text

text:
# A search for H₂ around pre-main-sequence stars

## 1 Introduction

Circumstellar disks play a vital part in both star and planet formation, and if we are to understand these processes then it seems we must understand the composition and physics of disks first. The discovery that disks can survive well into the main sequence lifetime of a star (after planet formation ended in our own system) has opened up an exciting window on the planet formation process (Zuckerman 2001). This means that we can directly probe the environment of planet formation, even though planets themselves may be undetectable.

The dominant constituent of the disks around pre-main-sequence stars is $`\mathrm{H}_2`$. However, observing emission from $`\mathrm{H}_2`$ which represents the bulk of the warm ($``$200K) gas in the disks is not straightfoward. Near-infrared (Bary et al. 2002), UV emission (Valenti et al. 2000) or absorption observations (Roberge et al. 2001) trace only the hottest gas or are dependent on the line of sight through the cirumstellar disk. Mid-infrared observations of the pure rotational emission lines of $`\mathrm{H}_2`$ trace the bulk of the warm gas, but are difficult from ground-based sites, due to telluric absorption and the high mid-infrared background emission. Thus, many previous studies of the gas component of such disks has relied on observations of tracer molecules, particularly CO (Mannings & Sargent 2000; Duvert et al. 2000). While these observations have been successful in tracing the velocity structure in disks, it seems that CO is a poor indicator of the total gas content. This is because CO can be destroyed by UV photons or frozen onto dust grains, leading to a severe underestimate of the total amount of gas. Most estimates of disk mass are based on sub-mm observations of the dust continuum (Beckwith et al. 1990; Sylvester et al. 2001), but these rely on assumptions about the composition and size distribution of the dust grains, and also on the gas/dust mass ratio. More importantly, we know that the dust and gas do not evolve in the same way – disks around main-sequence stars are composed mainly of dust, while disks around pre-main-sequence stars are composed mainly of gas. It is important to be able to measure gas properties independently from the dust properties using direct observations of $`\mathrm{H}_2`$.

Recently, the Infrared Space Observatory (ISO) has been used to search for $`\mathrm{H}_2`$ pure-rotational emission lines from pre-main-sequence stars, and also young main-sequence stars which are known to have dust disks. Thi et al. (2001b) report detections of the $`\mathrm{H}_2`$ $`J=31`$ transition ($`\lambda =17.035`$ $`\mu \mathrm{m}`$) and $`\mathrm{H}_2`$ $`J=20`$ transition ($`\lambda =28.221`$ $`\mu \mathrm{m}`$). $`\mathrm{H}_2`$ masses derived from these observations imply large masses of $`\mathrm{H}_2`$ in the disks around main-sequence targets (Thi et al. 2001a), reversing the previous belief that dust is the main component. However, as the ISO beam is large (14 $`\times `$ 27 arcsec) compared with the disks and the observations have moderate spectral resolution ($`R2000`$), these results cannot distinguish circumstellar disk emission from that of a foreground cloud or extended source.

Follow up observations of both the pure rotational emission lines and $`\mathrm{H}_2`$ in absorption have provided strong evidence that the warm $`\mathrm{H}_2`$ is extended in sources selected from the ISO studies. Using the the echelle spectrometer, TEXES, Richter et al. (2002) obtained observations of six sources, three of which have reported ISO $`\mathrm{H}_2`$ detections (Thi et al. 2001b): HD163296, AB Aur and GG Tau. They obtained upper limits on the line fluxes from five sources, with a possible 2 sigma detection of the $`\mathrm{H}_2`$ $`J=31`$ S(0) line in AB Aur and conclude that the majority of the $`\mathrm{H}_2`$ emission detection by ISO is extended on scales of 5 arcsec or more. The far-UV spectrum of Beta Pictoris presented by Lecavelier des Etangs et al. (2001) shows a complete absence of $`\mathrm{H}_2`$ absorption lines, contrary to expectations if large quantities of $`\mathrm{H}_2`$ are present within the edge-on dust disk.

We present the results of our attempts to confirm the ISO results using MICHELLE, the mid-infrared echelle spectrometer on the ground-based United Kingdom Infrared Telescope (UKIRT). We targeted two pre-main-sequence stars, CQ Tau (distance 100 pc) and AB Aur (distance 142 pc), both of which have good ISO detections and one of which (AB Aur) was also observed by Richter et al. (2002). In both cases, the dust disk is expected to be unresolved spatially by MICHELLE.

## 2 Observations and data reduction

Mid-IR echelle spectroscopy of AB Aur and CQ Tau was carried out on 2001 December 29 and 30 using MICHELLE (Glasse et al. 1997) on UKIRT. AB Aur and CQ Tau were observed in the $`\mathrm{H}_2`$ $`J=31`$ transition ($`\lambda =17.035`$ $`\mu \mathrm{m}`$), providing direct comparison with the ISO observations. AB Aur was also observed in the $`\mathrm{H}_2`$ $`J=42`$ transition ($`\lambda =12.279`$ $`\mu \mathrm{m}`$). The $`\mathrm{H}_2`$ $`J=20`$ transition ($`\lambda =28.221`$ $`\mu \mathrm{m}`$) was used in the ISO observations, but is not readily accessible from the ground due to the structure of the atmospheric transmission at 28 $`\mu \mathrm{m}`$. Details of the observations are shown in Table 1.

MICHELLE is a mid-IR imager and spectrometer which uses a 320$`\times `$240-pixel Si:As array and operates between 8 and 25 microns. For these observations the echelle grating was used, to provide maximum sensitivity to the $`\mathrm{H}_2`$ lines which are expected to be unresolved. AB Aur was observed with a 2-pixel slit width; CQ Tau with a 4-pixel slit width. The resolving power was $`R=\mathrm{15\hspace{0.17em}200}`$ at 12 $`\mu \mathrm{m}`$ and $`R=\mathrm{18\hspace{0.17em}600}`$ ($`R=\mathrm{9\hspace{0.17em}300}`$) at 17 $`\mu \mathrm{m}`$ for a two (four) pixel wide slit; as a result the $`\mathrm{H}_2`$ lines are well separated from the bright sky lines. A single order from the grating is selected using blocking filters. The pixel field of view varies with wavelength due to anamorphic demagnification. The values are given in Table 1.

The point-like sources were all observed by nodding along the 70 arcsec long slit, which was orientated north-south on the sky. The separation of the two beams was 15 arcsec and the dwell-time in each beam was 40 seconds.

Data reduction was done using software written in the idl programing language. Pairs of frames from different nod positions were subtracted, giving a set of difference frames each containing a positive and a negative version of the source. In each frame, the positive and negative versions were combined to give a single 2D spectrum. Before the set of 2D spectra were combined into a final image, the dataset was despiked using an iterative sigma cut. To do this, all of the 2D spectra were compared, and datum were flagged as bad if they were further than a specified number of standard deviations from the mean of their pixel. This sigma cut was done three times with threshold settings of 6, 4 and 3 sigma.

Once a final 2D spectrum had been produced, 1D spectra were extracted in two ways: a wide extraction which summed all the emission in a 6.5 arcsec region centered on the source, and an optimal extraction of the continuum. The optimal extraction method takes a weighted mean of the pixels around the source, giving more weight where the continuum is brighter. This method gives the lowest possible noise on the continuum level, but if the $`\mathrm{H}_2`$ emission region is more extended than the continuum emission then it will underestimate the $`\mathrm{H}_2`$ line flux. For both extraction methods a sky subtraction was done first, estimating the residual sky signal in each column using pixels between 3.5 and 7 arcsec from the source.

Standard stars were observed and reduced in the the same way as the science targets. The final spectrum for each object was divided through by a standard spectrum to remove wavelength dependant variations in atmospheric and instrumental transmission. The standard stars were observed using the same rows on the detector as the sources, and division by the standard effectively flat-fields the detector. For the 17 $`\mu \mathrm{m}`$ observations of AB Aur the standard star observation was saturated, so we were unable to make this correction. This causes variation in the continuum level, but would not prevent us from detecting a narrow emission line.

Wavelength calibration was obtained by measuring the positions of sky lines, and comparing the positions to a model sky spectrum. We estimate a 1 sigma error on the wavelength calibration is about one pixel ($`3\times 10^4`$ $`\mu \mathrm{m}`$), and is largely due to the difficulty in centroiding the broad sky lines. The spectra were shifted to the object’s rest frame, accounting for the radial velocity of the earth and object’s heliocentric radial velocity (+8.9 km s-1 for AB Aur and +22.6 km s-1 for CQ Tau (Barbier-Brossat & Figon 1999)). Flux calibration was based on the observed continuum level and published continuum flux measurements from Thi et al. (2001b) for both the 17 $`\mu \mathrm{m}`$ spectra and Richter et al. (2002) for the AB Aur 12 $`\mu \mathrm{m}`$ spectrum.

## 3 Results

Figures 1 and 2 show the final spectra produced using an optimal extraction. At the location of sky lines the noise increases dramatically, so these regions have been masked out (i.e. they were not used when estimating the continuum level and the overall noise). On the 17 $`\mu \mathrm{m}`$ figures a gaussian with the total line flux measured by Thi et al. (2001b) is shown, to allow comparison between the ISO results and our spectra. As the ISO observations do not resolve the line, we have assumed an intrinsic line width of 20 km s-1 (see Section 4 for discussion of this value), and an instrumental width of 16 km s-1 for AB Aur and 32 km s-1 for CQ Tau, giving a total width of 26 km s-1 and 38 km s-1 respectivley. The MICHELLE sensitvities should have produced a clear detection of the line, if all of the flux in the $`14\times 27`$ arcsec ISO beam were contained within our extracted beams, as expected for emission from an unresolved disk.

To estimate the line flux, all of the flux within a 40 km s-1 wide bin centered on the systematic velocity was summed, after subtracting the continuum level. Fluxes and upper limits for all extraction methods are shown in Table 2. There is no evidence of $`\mathrm{H}_2`$ emission in either of the 17 $`\mu \mathrm{m}`$ spectra. The 12 $`\mu \mathrm{m}`$ optimally extracted spectrum gives flux of $`4.8\pm 2.1`$ ergs s-1 cm-2, but as this is only a 2.3 sigma result should not be regarded as a definite detection. However, this result is consistent with the tentative detection by Richter et al. (2002), who measure $`(2.0\pm 1.0)\times 10^{14}`$ ergs s-1 cm-2 with a width of 10.5 km s-1.

## 4 Discussion

Despite the enhanced sensitivity to narrow emission lines from a point-like source, these new spectra of AB Aur and CQ Tau show no convincing detections of either the 0-0 S(1) or 0-0 S(2) lines. In order for the 0-0 S(1) $`\mathrm{H}_2`$ emission to be detectable by ISO but not by MICHELLE, the emission would have to be either spectrally broad or spatially extended. We consider these two explanations in turn.

The spectral resolving power of the ISO observations is $`\mathrm{R}(=\frac{\lambda }{\delta \lambda })=2400`$, and the $`\mathrm{H}_2`$ lines Thi et al. (2001b) observed are unresolved. In calculating our upper limits, we assume that the lines had an intrinsic width of 20 km s-1. The dominant line broadening mechanisms for gas in a disk are thermal broadening, and bulk motion due to the disk’s rotation. Assuming the gas has a temperature of 200 K (as estimated by Thi et al. (2001b)), the thermal broadening would only be 1.5 km s-1. The amount of rotational broadening depends on the distance from the star, and the inclination of the disk. If we assume a distance of 2AU (a reasonable choice, given the gas temperature of 200 K), this gives a width of 15 km s-1, assuming an inclination of $`45^{}`$. Therefore, 20 km s-1 is a conservative estimate of the line width. In fact, the line FWHM would have to be over 100 km s-1 to make the ISO and MICHELLE measurements consistent.

The possibility of extended emission, due to excitation of a remnant cloud of molecular gas by UV photons or by shocks, is plausible. Assuming that the ISO measured flux is uniformly distributed over the $`14\times 27`$ arcsec aperture, the fluxes are consistent with UV-excited emission strengths predicted by Burton et al. (1992) for moderate UV field strengths ($`10^310^4`$ times that of the interstellar medium) and moderate densities ($`n_{\mathrm{H}_2}10^310^4\mathrm{cm}^3`$). Over the MICHELLE beam used to obtain the spectra shown in Figure 1, the line flux would be an order of magnitude fainter than the sensitivity limit obtained. Moreover, the sky subtraction region (3.5-7 arcsec radius) means that our observations are insensitive to emission at radius $`>`$3.5 arcsec. The hypothesis of extended emission is strengthened for AB Aur since a map of 12CO 3-2 shows an extended molecular cloud around the central source, with a size of at least 30 arcsec (Thi et al. 2001b). For CQ Tau, no 12CO 3-2 line was detected, but this could be reconciled with extended $`\mathrm{H}_2`$ emission if the CO has been destroyed by photodissociation (Dent et al. 1995).

In interpreting our results, an important question concerns the optical depth of the dust disk. If the dust is optically thick at the wavelengths we observed and the gas and dust are at the same temperature, then the gas would be equally likely to absorb or emit a photon, so no emission line would be produced. Line emission will only occur if there are gaps in the dust disk, or if there is a hot atmosphere above the disk. Though these situations are quite plausible, only a small fraction of the gas in the disk would be observed in this case, and any calculation of the gas mass would underestimate the true mass. It is therefore essential that we know the optical depth of the disk in order to correctly interpret our results.

The optical depth of the disks can be estimated from radiative transfer models based on the SED of the disks. Natta et al. (2001) made models of the disk SEDs for both AB Aur and CQ Tau, and we can use their fitted parameters to determine whether the disks are optically thick at 12 and 17 $`\mu \mathrm{m}`$. Based on their estimates of disk mass and disk radius, and using the usual assumption that the dust opacity at 1.3 mm is $`\kappa _{1.3\mathrm{mm}}=0.01`$ cm2 g-1, the average optical depth of the disks at 1.3 mm is 0.03 for AB Aur and 0.24 for CQ Tau. Then, using a power law dust opacity index ($`\kappa \lambda ^\beta `$, where $`\beta `$ was estimated to be 2.0 for AB Aur and 1.0 for CQ Tau in Natta et al. 2001), the optical depth can be estimated at shorter wavelengths. This indicates that the disks become optically thick at wavelengths shorter than 200 $`\mu \mathrm{m}`$ for AB Aur and 300 $`\mu \mathrm{m}`$ for CQ Tau. Both disks must therefore have a high optical depth at 12 and 17 $`\mu \mathrm{m}`$.

This result has two important implications. Firstly, this means that we cannot use our upper limits on line flux to calculate upper limits on the mass of $`\mathrm{H}_2`$ in the disk, as most of the gas in the disk is hidden by the optically thick dust. The $`\mathrm{H}_2`$ line at 12 $`\mu \mathrm{m}`$ (if confirmed) must arrise from a hot atmosphere above the disk, or from gaps in the dust disk. Secondly, if the the lines detected by ISO do arise in the disk, then only a small fraction of the gas was observed and the total gas mass must be much higher than was estimated.

In conclusion, we have searched for pure rotational $`\mathrm{H}_2`$ emission lines from the circumstellar disks of AB Aur and CQ Tau. We tentatively detect emission from AB Aur, but find no evidence for emission from CQ Tau. The likely reason for this non-detection is that the disks are optically thick in the mid-IR. Our upper limits on line flux are significantly smaller than previous line flux estimates based on ISO observations, which suggests that the emission detected by ISO is extended on scales of at least 6 arcsec (or $`>`$ 60AU).

Mapping the emission over the area covered of the ISO beam is one approach to confirming these conclusions. Detecting the faint extended emission will require a considerable gain in sensitivity over these observations, such as is anticipated using MICHELLE on the Gemini North 8-m telescope.

The United Kingdom Infrared Telescope is operated by the Joint Astronomy Centre on behalf of the UK Particle Physics and Astronomy Research Council. We wish to thank the staff at UKIRT for enabling these observations. We also thank Mark Wyatt and Alistair Glasse for helpful discussions of this work.


text:
# Layer charge instability in unbalanced bilayer systems in the quantum Hall regime

## Abstract

Measurements in GaAs hole bilayers with unequal layer densities reveal a pronounced magneto-resistance hysteresis at the magnetic field positions where either the majority or minority layer is at Landau level filling factor one. At a fixed field in the hysteretic regions, the resistance exhibits an unusual time dependence, consisting of random, bidirectional jumps followed by slow relaxations. These anomalies are apparently caused by instabilities in the charge distribution of the two layers.

Hysteretic phenomena are widespread in nature. They are common magnetic materials, and often indicate a non-equilibrium situation associated with a phase transition and the presence of domains berttoti . Recently, hysteresis has also been reported in various two-dimensional (2D) carrier systems in semiconductor structures at low temperatures and high magnetic fields cho ; kronmueller ; piazza ; eom ; etienne ; smet . In these cases, magneto-resistance ($`\rho _{xx}`$) hysteresis appears in the quantum Hall (QH) regime when two Landau levels (LLs) with opposite spin are brought into coincidence. While the 2D systems studied have been notably different, the common thread in these experiments is that there is a magnetic transition involving the carrier spin tomas .

Here we present hysteretic $`\rho _{xx}`$ data in 2D bilayer systems in the QH regime. The hysteresis in these systems has a different origin and is caused by a non-equilibrium charge distribution in the two layers. We studied the magneto-transport coefficients of GaAs bilayer hole systems with unequal layer densities. When the interlayer tunneling is sufficiently small, $`\rho _{xx}`$ of the bilayer system exhibits a pronounced hysteresis at perpendicular magnetic field ($`B`$) positions close to where either the majority or minority layer is at LL filling factor one. Most remarkable is the time dependence of $`\rho _{xx}`$ at a fixed field in the hysteretic regime, when the two layers are closely spaced. As a function of time, $`\rho _{xx}`$ exhibits large, random, sudden jumps toward higher and lower values, followed by a slow decay in the opposite direction. The data may signal an instability in the charge distribution of the two layers, i.e., an instability associated with the pseudospin (layer), rather than spin, degree of freedom.

We studied nine GaAs bilayer hole samples from six different wafers, all grown on GaAs (311)A substrates and modulation doped with Si. In all samples, the holes are confined to two 15nm-wide GaAs quantum wells which are separated by AlAs or AlAs/AlGaAs barriers with thickness $`7.5W200`$nm. The rather thick barrier combined with the large effective mass of GaAs 2D holes mass  reduces considerably the tunneling between the two layers tunneling . As grown, the samples have layer densities of $`7\times 10^{10}`$ cm-2, and low temperature ($`T`$) mobilities of $`35`$ m2/Vs. Metallic top and bottom gates were added to control the densities in the layers. We studied several types of devices, including 2.5$`\times `$2.5mm square samples and ones with patterned Hall bars; in these samples the ohmic contacts contact both layers. One sample was fabricated using a selective depletion scheme eisenstein-APL  that allows probing the transport characteristics of individual layers. The measurements were performed in a dilution refrigerator down to $`T=20`$mK.

Data of Fig. 1 highlight some of the results of our study. In (a) we show a set of traces where $`\rho _{xx}`$ was measured for a sample with $`W=11`$nm as $`B`$ was ramped up or down. For the traces of Fig. 1(a) the total bilayer density ($`p_{tot}`$) is kept constant at $`5.5\times 10^{10}`$ cm-2 while charge is transferred from one layer to another using back- and front-gate biases. We define the charge transfer from one layer to another as $`\mathrm{\Delta }p=(p_Bp_T)/2`$, where $`p_B`$ and $`p_T`$ are the densities of bottom and top layers, respectively note2 . At a given value of $`B`$, we define the filling factor, $`\nu `$, of the bilayer system as the ratio between $`p_{tot}`$ and the LL degeneracy, $`eB/h`$. We also introduce the filling factors for top and bottom layers, $`\nu _T`$ and $`\nu _B`$ respectively, as the ratio between the layer density and $`eB/h`$.

The data of Fig. 1(a) show that when the bilayer system is balanced ($`\mathrm{\Delta }p`$=0; top trace) $`\rho _{xx}`$ is independent of the direction $`B`$ is ramped. However, as soon as the system is imbalanced ($`|\mathrm{\Delta }p|>0`$) a strong hysteresis develops in $`\rho _{xx}`$. For values of $`\mathrm{\Delta }p<4.7\times 10^9`$ cm-2, $`\rho _{xx}`$ displays hysteretic behavior in two field ranges near $`\nu =2`$, one near $`\nu _T=1`$ and another near $`\nu _B=1`$. When $`\mathrm{\Delta }p4.7\times 10^9`$ cm-2, the hysteresis exists only near $`\nu _B=1`$. For sufficiently large $`\mathrm{\Delta }p`$, no hysteresis is observed. The amplitude of the hysteresis also decreases as $`T`$ is increased (data not shown) and vanishes completely above $`T`$230mK, roughly independent of $`\mathrm{\Delta }p`$.

To probe the contribution of the spin degree of freedom to the hysteresis in our bilayer systems, we performed measurements in tilted magnetic fields on a sample very similar to the one shown in Fig. 1(a). In this experiment the direction of the field was kept at an angle $`\theta `$ with respect to the normal to the plane of the bilayer system. For $`\theta `$ ranging from $`0^{}`$ to $`80^{}`$, corresponding to a six-fold increase of the total field (and therefore of the Zeeman energy) in the hysteretic region, the position of the hysteresis in perpendicular magnetic field did not change at all. If the hysteresis were caused by an instability associated with the spin degree of freedom, one would expect etienne ; smet  that the applied parallel field would change the position and magnitude of the hysteresis. These results rule out spin as being responsible for the hysteresis.

To better understand the origin of the observed hysteresis, we fabricated another sample from a different wafer, also with $`W=11`$nm, using a selective depletion scheme eisenstein-APL , and aimed to independently probe each layer of the bilayer system. The data are shown in the right panel of Fig. 1. In (b) we plot $`\rho _{xx}`$ for the bilayer system, that is, when the ohmic contacts are connected to both layers, for both up and down $`B`$-sweeps. In (c) and (d) we show $`\rho _{xx}`$ traces for the top and bottom layers, measured separately, but at the same pair of layer densities as in panel (b) traces. Two features of these data are noteworthy. First, the traces of (b) exhibit hysteresis in two ranges of $`B`$. The hysteresis between 1.45 and 1.8T matches well the position of $`\nu _T=1`$ QH state of the top layer as seen in (c), while the hysteresis located between 2.5 and 3.2T overlaps the QH state of the bottom layer (see (d)). This observation confirms that the hysteresis in $`\rho _{xx}`$ of the bilayer system takes place when one of the layers is at filling factor one. Second, data of (c) and (d) show that each individual layer exhibits hysteresis when the other layer is at filling factor one, i.e., $`\rho _{xx}`$ of the top layer exhibits hysteresis when $`\nu _B=1`$, and vise versa.

The observation of hysteretic magneto-resistance in unbalanced bilayer systems has precedence. Zhu et al. zhu  reported hysteresis in 2D electron systems with a parallel conducting layer. In their case, the parallel layer was a parasitic, low-mobility, doping layer at a distance of 75 to 95 nm away from the high-mobility 2D electrons. Similar to our data, they observed hysteresis in $`\rho _{xx}`$ when the layer containing the high-mobility 2D electrons was in a QH state. They also presented a simple model to explain the observed hysteresis. In their model, the hysteresis comes about because of a non-equilibrium charge distribution in the layers. As $`B`$ is swept, thanks to the Landau quantization, the Fermi levels of both layers oscillate. These oscillations lead to temporary imbalances between the chemical potentials of the two layers. The potential imbalance is particularly abrupt and large when the high-mobility 2D electrons enter a QH state as their Fermi level jumps by a significant amount, equal to the separation between the adjacent LLs. With increasing time, of course, the Fermi levels of the two layers have to come to equilibrium since the latter are shorted together via the ohmic contacts. But this equilibration can take a long time in the QH effect regime: it has to take place via the layers’ edges and the ohmic contacts since the (bulk) states in the center of the 2D layer that is in the QH state are localized and the layer sheet conductivity is very small. As a result, $`\rho _{xx}`$, which is recorded as $`B`$ is swept at a finite rate, can show a hysteretic behavior. Consistent with their model, Zhu et al. found that when the $`B`$-sweep is interrupted in the hysteretic region and $`\rho _{xx}`$ is monitored as a function of time, it decays approximately exponentially toward an equilibrium value. Moreover, they found the time constant of the decay in reasonable agreement with estimates based on the parameters of their experiment.

The hysteresis in our samples resembles what Zhu et al. observe and likely has a similar origin. Our data of Fig. 1 in fact explicitly show that the hysteresis happens when one of the layers is in a QH state, and it is the resistivity of the layer which is not in the QH state that is hysteretic. The interpretation that the hysteresis indicates a charge transfer between the two layers is also supported by our data: the hysteresis observed in $`\rho _{xx}`$ of either the top or the bottom layer measured separately appears much like a horizontal shift of the trace in $`B`$, as if the layer had slightly different density when ramping $`B`$ up or down \[see, e.g., Fig. 1(c)\]. The time evolution of $`\rho _{xx}`$ in our samples, however, is qualitatively different from the observations of Zhu et al., and points to a very unusual relaxation. In the remainder of the paper, we describe this evolution and speculate on its possible origin.

We have studied a number of samples with varying $`p_{tot}`$, $`\mathrm{\Delta }p`$, and barrier width $`7.5W200`$nm. Samples with $`W=7.5`$nm do not show any hysteresis; this is likely because the interlayer tunneling is sufficiently large so that the two layers stay in equilibrium during the $`B`$-sweep. We observe hysteresis for all samples with $`11W200`$nm, but the time evolution of $`\rho _{xx}`$ critically depends on $`W`$. Examples are shown in Figs. 2 and 3. For samples with $`W=11`$nm (Fig. 2), the time evolution is simply wild! It displays sudden jumps in $`\rho _{xx}`$, followed by a slow relaxation after each jump. Note that $`\rho _{xx}`$ jumps toward both higher and lower values, reminiscent of a bistability, although the jumps do not happen between fixed values of $`\rho _{xx}`$. In between jumps, $`\rho _{xx}`$ follows a slow relaxation, in the opposite direction of the jump, that can be fitted well by a double exponential \[see Fig. 2(b) inset\]. It is noteworthy that even when measured over days (up to 2.5$`\times 10^5`$s), we did not observe any tendency towards a settling of the jumps. We wish to emphasize that, outside the hysteretic region, $`\rho _{xx}`$ is independent of time to within less than 0.3%.

We have attempted to quantify the characteristics of this evolution by two parameters: average frequency and amplitude of the jumps. Our $`T`$-dependence measurements show that, at a fixed $`B`$, the average jump amplitude decreases as $`T`$ increases. Not surprisingly, the jumps are no longer visible above the temperature where the hysteresis vanishes. On the other hand, the average frequency appears to be independent of $`T`$. Also, the jump frequency and amplitude are independent of the magnitude of the sample current, as long as the current is kept sufficiently small ($`10`$ nA) jouleheating . Interestingly, the jumps and decays appear to continue even when the current is completely turned off note3 .

The sample with $`W=30`$nm shows a behavior qualitatively similar to the one with $`W=11`$nm, although both the frequency and size of $`\rho _{xx}`$ jumps are smaller \[Fig. 3(a)\]. Data for the $`W=50`$nm sample, however, are qualitatively different \[(Fig. 3(b)\]: there are by far fewer $`\rho _{xx}`$ jumps (typically one jump every few $`10^3`$ s), and $`\rho _{xx}`$ appears to decay with time. Finally, for the $`W=200`$nm sample \[Fig. 3(c)\], we typically observe a simple decay with time jumps-in-200nm . Interestingly, for the samples of Figs. 3(b) and (c), $`\rho _{xx}`$ continues to decay with an ever increasing time constant.

The time evolutions we observe for the sample with $`W=200`$nm barrier \[Figs. 3(c)\] is qualitatively similar to the observation of Zhu et al. zhu . The time dependences for the samples with barrier widths $`W=11`$ and 30nm, however, are very unusual and cannot be understood in a simple model where the bilayer system slowly and steadily relaxes to an equilibrium state. In these samples, $`\rho _{xx}`$ displays sudden jumps that do not have any tendency to settle, at least over a time scale of days. We do not know the origin of these time evolutions. The sudden jumps in $`\rho _{xx}`$ in our bilayers with small $`W`$ bear some resemblance to the so-called Barkhausen jumps, which are observed in magnetic materials berttoti . The Barkhausen jumps occur when the magnetic system finds a lower energy state available and one or several domains change orientation. A tantalizing speculation is that in the bilayer systems with close layer separation, the interlayer interaction acts as an opposing force to the charge transfer caused by the Fermi level difference dlb  . In this scenario, the two opposing mechanisms may mediate the creation of a complicated layer charge density pattern or pseudospin domains.

Another possibility is that the observed jumps are not intrinsic to the sample, but rather triggered by external sources (e.g. electromagnetic noise). If so, it is a puzzle why the jumps are much more frequent in the bilayer samples with smaller layer separation and are seen only in the hysteretic region.

In summary we report an unusual time dependence associated with hysteretic magneto-resistance in GaAs bilayer holes with close layer separation. The resisitivity exhibits sudden jumps with time, possibly caused by a layer charge instability.

We thank DOE and NSF for support and D. Haldane for helpful discussions.


text:
# Dynamo action in turbulent flows

## 1 Introduction

The dynamo problem is mainly related to the understanding of the processes of magnetic field generation due to the motions of conducting fluids, the so-called dynamo action (Cowling Cowling1934 (1934), Moffatt Moffatt1978 (1978), Parker Parker1979 (1979)). One often divides dynamo theory into two theoretical regimes: First, the kinematic dynamo problem in which magnetic fields and flows are decoupled. The velocity of the flow is prescribed and the evolution of the magnetic field is governed by the induction equation (here in dimensionless form):

$$\frac{𝐁}{t}=\times (𝐮\times 𝐁)+\frac{1}{\mathrm{Re}_\mathrm{m}}^2𝐁,$$
(1)
where B is the magnetic field and u is the prescribed velocity field. The quantity $`\mathrm{Re}_\mathrm{m}=u\mathrm{}/\eta `$ is the magnetic Reynolds number and it is often huge in astrophysical systems, where $`u`$ and $`\mathrm{}`$ are characteristic velocity and length scales and $`\eta `$ is the magnetic diffusivity. A flow acts as a fast dynamo if the growth rate does not tend to zero in the limit of infinite magnetic Reynolds number. The solar dynamo is an astrophysical example of fast dynamo action since it operates on the convective time scale of the fluid which is fast compared to that of the Ohmic diffusion. In the fast kinematic dynamo problem one seeks velocity fields that lead to the exponential amplification of the magnetic field through induction Eq. (1). A powerful mechanism for growth of the magnetic field is the stretching of the magnetic field lines (cf. Archontis et al. Archontis+ea03 (2003)), and is primarily achieved by chaotic flows. The kinematic approach is valid when the magnetic field is weak and there are no further dynamical effects. In the second, fully dynamical dynamo regime, the velocity field is not prescribed and the exponential growth of the magnetic field saturates when the field becomes strong enough to modify the flow sufficiently through the feedback by the Lorentz force. The saturation state of the dynamo then emerges as a self consistent solution to the equations of nonlinear magneto-hydrodynamics. In a fully self-consistent dynamical experiment the forcing would have internal causes, e.g. of convective origin—here we are concerned with the case of an externally prescribed forcing function.

The nonlinear properties of fast dynamos have received considerable attention and a variety of both analytical and numerical studies have provided a valuable insight into the nature of the dynamo action and saturation (Nordlund et al. Nordlund+ea92 (1992); Brandenburg et al. Brandenburg+ea95 (1995); Cattaneo et al. Cattaneo+ea96 (1996); Zienicke et al. Zienicke+ea98 (1998); Brummell et al. Brummell+ea98 (1998)). However, a difficulty that arises in numerical simulations is that often, fast dynamo flows become hydrodynamically unstable at fluid Reynolds number (Re$`=u\mathrm{}/\nu `$, $`\nu `$ being the viscosity) greater than a critical value (Podvigina & Pouquet Podvigina+Pouquet1994 (1994)) and spatio-temporal turbulence appears to be nascent for higher values of Re. As a result the desired flow velocity is modified even in the kinematic regime and the problem contains the complexity of that of a turbulent dynamo. Thus, the choice of the forcing function that drives the flow is crucial. We use a forcing function with an amplitude that is allowed to evolve with time, keeping the average kinetic energy approximately constant through both laminar and turbulent phases. The growth rates and the saturation level of such turbulent dynamos are examined with respect to the magnetic and fluid Reynolds numbers. Also of great interest is the understanding of the processes at work in the kinematic (both laminar and turbulent phase of the flow) and in the saturated regime. We present an analysis of the Lorentz work and Joule dissipation, which shows that they are in close balance and that dynamo-work occurs mainly in the weak part of the magnetic field where dissipation is insignificant.

The paper is organized as follows. Section 2 contains the equations, the numerical method and a description of the driving of the flow that has been used in the simulation. Results for turbulent dynamo action are presented in Section 3. The structures of the velocity field and the magnetic field in physical space are discussed in Section 4. Whether or not the growth rates and the saturation level of the dynamo depend on the magnetic or fluid Reynolds number is examined in Section 5 and Section 6. A qualitative understanding of the nature of the dynamo is discussed in Section 7, using results from a kinematic dynamo experiment. Section 8 contains the overall conclusions of the present numerical simulations.

## 2 The simulation

### 2.1 The equations

The compressible MHD equations are solved numerically in a periodic computational domain, with periodicity of $`2\pi `$ in all three directions. Apart from the magnetic induction equation Eq. (1), the additional equations solved are the following:

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\rho 𝐮,`$ (2)
$`{\displaystyle \frac{(\rho 𝐮)}{t}}`$ $`=`$ $`P+𝐣\times 𝐁+𝐟(\rho \mathrm{𝐮𝐮}),`$ (3)
$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`(e𝐮)P𝐮+\mathrm{Q}_\mathrm{v}+\mathrm{Q}_\mathrm{J}+\mathrm{Q}_{\mathrm{cool}}`$ (4)
where $`\rho `$ is the fluid density, $`P`$ the pressure, $`𝐣`$ is the electric current density, $`e`$ is the internal energy and $`𝐟`$ is an external forcing term. $`\mathrm{Q}_\mathrm{v}`$ and $`\mathrm{Q}_\mathrm{J}`$ are the viscous and Joule dissipation respectively. In the experiments, a Newtonian cooling term, $`\mathrm{Q}_{\mathrm{cool}}=(TT_0)/\tau _{\mathrm{cool}}`$ is used, where $`T`$ is the temperature and $`\tau _{\mathrm{cool}}=5`$. The cooling term only works as a sink of the thermal energy, in order to balance the heat generation by dissipative terms in the energy equation and the actual magnitude of the cooling time is of little importance with no effect on the results of our simulation.

The above equations are solved numerically on a staggered mesh. The time stepping is performed by a third order predictor-corrector method (Hyman Hyman1979 (1979)). The derivatives and interpolations are of 6th and 5th order respectively and the numerical scheme conserves $`𝐁=0`$ exactly. Numerical solutions are obtained on a grid of a maximum of $`180^3`$ points, using a modified version of the code by Galsgaard, Nordlund and others (Galsgaard & Nordlund Galsgaard+Nordlund1997 (1997); Nordlund et al. Nordlund+ea92 (1992)). The initial magnetic field is chosen to be a weak random perturbation with an amplitude of $`10^5`$ in non-dimensional units.

### 2.2 Flow considerations

There is a considerable body of work in the literature that deals with fast dynamo action in simple, steady and three-dimensional flows. In our simulation the velocity is chosen to have vigorous (fast) dynamo properties in the kinematic regime, where the Lorentz force is negligible. One such class of flows, which are well known candidates for fast dynamo action is the ABC flows. The form of the velocity is given by

$`𝐮_{\mathrm{ABC}}`$ $`=`$ $`A(0,\mathrm{sin}kx,\mathrm{cos}kx)+B(\mathrm{cos}ky,0,\mathrm{sin}ky)`$
$`+`$ $`C(\mathrm{sin}kz,\mathrm{cos}kz,0)`$
This periodic flow is the sum of three steady Beltrami waves, parameterized by $`A`$, $`B`$ and $`C`$ and has the property $`\times 𝐮𝐮`$ (velocity is parallel to vorticity). If one or more of the constant coefficients ($`A`$, $`B`$ and $`C`$) is zero, the flow is integrable and is not a fast dynamo. If all three are non-zero the flow is non-integrable (Dombre et al. Dombre+ea86 (1986)) and contains a mixture of chaotic regions and regular islands. The special case $`A=B=C=1`$ was introduced by Childress (Childress+ea70 (1970)) as a model for kinematic dynamo action. Arnold et al. (Arnold+ea83 (1983)) noted that a steady three-dimensional flow with infinite conductivity and chaotic streamlines favors the growth of magnetic fields. Their results were followed by the work of Galloway et al. (Galloway+ea84 (1984)) and Moffatt et al. (Moffatt+ea85 (1985)) who studied dynamo action in ABC flows with finite conductivity. A detailed description of the amplification process responsible for the kinematic dynamo action, when the wavenumber of the flow $`k`$ is equal to one or higher, has been provided by Dorch (Dorch2000 (2000)) and Archontis et al. (Archontis+ea03 (2003)). The hydrodynamic stability of the ABC flow has been studied by Podvigina & Pouquet (Podvigina+Pouquet1994 (1994)) who found that for values of fluid Reynolds number above a critical one ($`\mathrm{Re}_\mathrm{c}=13`$) the ABC flow destabilizes, first to time dependent but still smooth states, and then to a turbulent state. In our simulation, the definition of the Reynolds numbers is the same as in Galanti et al. (Galanti+ea92 (1992)).

In the following experiments the initial flow is taken to be an ABC flow with $`A=0.9`$, $`B=1`$, $`C=1.1`$ and wavenumber equal to unity. The reason for choosing values of $`A`$, $`B`$, and $`C`$ slightly different from unity, is that the $`A=B=C`$ case is too special, because of its exact symmetry, and thus turbulence takes longer to develop. The purpose of our work is to study dynamo action in the turbulent regime of driven ABC flows since astrophysical dynamos typically occur in turbulent environments with very high fluid and magnetic Reynolds numbers. However, the behavior of turbulent dynamos in the high Re, $`\mathrm{Re}_\mathrm{m}`$ limit may be expected to be to some extent generic. An increase in the Reynolds numbers involves smaller and smaller length scales where dissipation actually occurs. However, the rate of dissipation is determined by the large scale developments of the flow and is independent of the magnitude of the viscosity once the Reynolds numbers are large. Thus, there is good hope that also dynamo properties are asymptotically independent of the Reynolds numbers. That hope is indeed borne out by the present numerical experiments, as is discussed in more details below.

## 3 Results

As a starting point, the forcing function which is imposed is such that the ABC flow Eq. (2.2) appears. The forcing necessary to accomplish this must have the following form:

$$𝐟=\frac{1}{\mathrm{R}_\mathrm{e}}^2𝐮_{\mathrm{ABC}}$$
(6)
However, in order to ensure that the kinetic energy is maintained at a value close to the initial value, the amplitude of the forcing is controlled by a differential equation, which compensates for the increased dissipation in the turbulent phase by increasing the forcing when the kinetic energy tends to drop below the nominal value. The amplitude factor in front of the driving force is given by

$$K=\frac{\mathrm{E}_0}{\mathrm{E}_{\mathrm{kin}}}e^L,$$
(7)
where $`L`$ is determined from

$$\frac{dL}{dt}=\frac{1}{\tau _{\mathrm{t}.\mathrm{o}.}}ln\frac{\mathrm{E}_0}{\mathrm{E}_{\mathrm{kin}}},$$
(8)
with $`L(t=0)=0`$. The differential equation evolves the amplitude factor $`e^L`$ on the turn over time scale, $`\tau _{\mathrm{t}.\mathrm{o}.}`$, if the kinetic energy tends to deviate from the nominal value $`\mathrm{E}_0`$. The immediate factor $`\mathrm{E}_0/\mathrm{E}_{\mathrm{kin}}`$ helps to reduce the time delay that results from the integral nature of the exponential factor, and the delayed response of the kinetic energy on changes in the driving. The kinetic energy control is an essential feature, in that otherwise the amplitude of the velocity in the turbulent phase would become arbitrarily small, at increasingly large fluid Reynolds number, since the forcing in the laminar phase needs to scale with $`1/\mathrm{Re}`$. It is indeed necessary to maintain the kinetic energy in order to maintain the same turn over time and actual fluid Reynolds number, since the formal fluid Reynolds number is defined in terms of the RMS velocity of the laminar ABC flow.

Figure 1 shows the evolution of the kinetic energy in time for the first 100 time units. Two phases may be identified. The laminar phase which extends up to $`t=30`$ and the turbulent regime for $`t>30`$. Note that the kinetic energy drops by up to about 40% when the flow enters the turbulent phase. After a short adjustment time, the forcing control returns the kinetic energy to the nominal level, where it is then maintained, except for minor fluctuations.

The magnetic energy for the same time period during the laminar phase increases exponentially and the growth rate is about 0.15. When the flow enters the turbulent phase ($`t>30`$) the magnetic energy drops slightly and later ($`t=40`$) continues to grow exponentially with a similar growth rate as in the laminar phase (Fig. 2).

An analysis of the external work in the above two distinct phases shows that the average work in the turbulent regime is much larger than in the laminar case (Fig. 3). On the first hand, the asymptotic level of the average work ($`0.25`$) is independent of Re and incidentally corresponds closely to the level that would be obtained for the critical value ($`\mathrm{Re}_\mathrm{c}`$). On the other hand, the average work in the laminar regime is proportional to $`1/\mathrm{Re}`$, and can become arbitrarily small for high enough Reynolds numbers.

This illustrates the point made above, about the forcing control; without it, the amplitude of the velocity in the turbulent phase would drop very much and the kinetic energy could not be maintained in order to have constant fluid Reynolds number during the simulation.

The enhancement of the driving force above the laminar value is shown in Fig. 4. The factor in Eq. (7) is by definition $`=1`$ during the laminar phase of the flow but increases during the turbulent phase to maintain the kinetic energy. The asymptotic level of the factor is about 10 for $`\mathrm{Re}=100`$.

Figure 5 shows the viscous dissipation as a function of time for the laminar and turbulent phase. During the laminar phase the dissipation is constant and is balanced by the external work. The viscous dissipation increases when the flow becomes turbulent. It reaches an asymptotic level when the turbulent velocity is maintained at the same level as that of the laminar flow.

Figure 6 compares the Joule dissipation to the work done against the Lorentz force and also shows the temporal evolution of the magnetic energy for the kinematic regime of the dynamo.

The work done against the Lorentz force is leading the Joule dissipation and is ultimately responsible for the increase in the magnetic energy. There is a close balance between them with a small positive difference during the increase of the magnetic energy and small negative difference when the magnetic energy decreases, see Fig. 6a, Fig. 6b. The fluctuations of the Lorentz work lead similar fluctuations of the Joule dissipation when the magnetic energy increases or decreases. Their difference in the saturated state of the dynamo fluctuates around zero, as expected. The work done by the fluid on the magnetic field is, on average, converted into Joule heat and none is left to increase the magnetic energy.

The accuracy in the numerical method is demonstrated in Fig. 6c, where the $`d\mathrm{E}_{\mathrm{mag}}/dt`$ is overplotted (symbols) on top of the difference $`\mathrm{W}_\mathrm{L}\mathrm{Q}_\mathrm{J}`$. These are averages over the periodic box, for which the flux term in the magnetic energy equation vanish. The magnetic energy equation is

$$\frac{\mathrm{E}_{\mathrm{mag}}}{t}=𝐅_{\mathrm{poynt}}𝐮𝐟_\mathrm{L}\mathrm{Q}_\mathrm{J}.$$
(9)
The source of magnetic energy is the work done against the Lorentz force and the sink is the Joule dissipation. Magnetic energy is transported by the Poynting flux, $`𝐅_{\mathrm{poynt}}=𝐄\times 𝐁`$. According to Eq. (9), $`d\mathrm{E}_{\mathrm{mag}}/dt`$ should be identically equal to $`\mathrm{W}_\mathrm{L}\mathrm{Q}_\mathrm{J}`$, and the two are indeed very nearly the same.

Note that the ratio of the Lorentz work to the Joule dissipation in Fig. 6d remains close to unity both during intervals of time when the magnetic energy is growing and during intervals of time when the magnetic energy is decaying.

## 4 Magnetic and velocity patterns

During the laminar phase the flow has points with two-dimensional stable manifolds (so called $`\alpha `$-type stagnation points), where the magnetic field is concentrated along magnetic “flux cigars”. There are also points with two-dimensional unstable manifolds, similar to the $`\beta `$-type stagnation points of the A=B=C=1 flow, where magnetic flux sheets are formed (see Dorch Dorch2000 (2000)).

Figure 7 (left panel) shows that flux cigars form an inclined triangle around a $`\beta `$-type stagnation point (located at the center of the box). Another one, which is located at the upper left corner, is pointing towards the plane of the $`\beta `$-type point while the flow is still laminar. The right panel of the same figure shows low velocity and high magnetic field isosurfaces at $`t=30`$. The velocity field comes into a turbulent phase and the regions of low velocity are disturbed, but the triangle formed by the flux cigars is still visible. The low velocity isosurfaces no longer contain points of exact stagnation, but still correspond to regions of vigorous stretching (see Archontis et al. Archontis+ea03 (2003)), which is why we choose to visualize them.

Figure 8 is a visualization of the same velocity and magnetic field strength for $`t=80`$. The velocity field has a more turbulent configuration and the regions with velocities less that $`4\%`$ of the peak have disappeared. The magnetic flux cigars also disappear and the magnetic energy is confined into curved tubes and weaker sheets.

## 5 Saturation of the dynamo

Two experiments were performed to show the dependence of the saturation level of the turbulent dynamo on the magnetic and fluid Reynolds number. Cases with $`\mathrm{Re}_\mathrm{m}=\mathrm{Re}=200`$ and $`\mathrm{Re}_\mathrm{m}=\mathrm{Re}=400`$ are studied and the temporal evolution of the magnetic energy is followed in the kinematic and saturated regime.

Figure 9 shows the growth of the magnetic energy in the beginning of the kinematic regime for $`\mathrm{Re}_\mathrm{m}=\mathrm{Re}=200`$ and $`\mathrm{Re}_\mathrm{m}=\mathrm{Re}=400`$. The numerical resolution used were $`80^3`$ and $`160^3`$, respectively. After an initial transient phase the growth rate for the second experiment is very similar to the growth rate of the first one. The difference in the magnetic energy amounts to a factor of 2.7. Thus if the magnetic energy for the $`160^3`$ experiment is divided with 2.7, it falls right on top of the $`80^3`$ experiment with $`\mathrm{Re}_\mathrm{m}=100`$ (see Fig. 9).

Keeping the same values of magnetic and fluid Reynolds number, the temporal evolution of the magnetic energy is examined in the saturated regime. The magnetic energy eventually stops growing exponentially and saturates. A result of considerable significance is that for higher values of the fluid Reynolds number, but constant magnetic Prandtl number $`(\mathrm{Pr}_\mathrm{m}=\nu /\eta )`$, the saturation level of the turbulent dynamo increases, rather than decreases (Vainshtein & Cattaneo Vainshtein+Cattaneo1992 (1992)).

Figure 10 shows the temporal evolution of the magnetic energy over a small epoch of the non-linear regime, illustrating the effect of increased numerical resolution. There is an initial transient, as the larger resolution allows field concentrations to collapse further. Note that there are then sections of time where the growth rates and the decay rates are quite similar in the $`\mathrm{Re}_\mathrm{m}=400`$ and $`\mathrm{Re}_\mathrm{m}=200`$ cases. The sections after the largest maximum in the two cases, and the rise towards the next maximum are two examples. This might be fortuitous, but more likely it corresponds to parallel evolution of similar magnetic structures.

The overall level is slightly higher, and the evolution is slightly delayed in the high $`\mathrm{Re}`$ case, because flux structures that were barely resolved at $`80^3`$ are allowed to collapse in the $`160^3`$ case. This is indeed also what happens in the kinematic case, but here the size of most structures is controlled by the Lorentz force, and not by diffusion, so the effect is only marginal. For the same reason, the effect may be expected to vanish altogether as the magnetic Reynolds number goes to infinity. In addition, the size of the magnetic structures in the kinematic regime scales as $`\mathrm{Re}_\mathrm{m}^{1/2}`$. In the saturation regime the size is again determined by the balance between the advection and dissipation of the magnetic energy and presumably changes following the same scaling with no crucial dependence on the value of $`\mathrm{Re}`$.

After some time, the saturation level for $`\mathrm{Re}_\mathrm{m}=400`$, indeed fluctuates around a level that is only marginally larger than for $`\mathrm{Re}_\mathrm{m}=200`$. This is small compared to the factor 2.7 for the kinematic regime of the dynamo.

One also notices that there are more fluctuations in the higher resolution case (e.g. the extra bump on the way down after the largest maximum), as a result of the small scales resolved when going to higher $`\mathrm{Re}_\mathrm{m}`$. That trend is expected to continue with increasing magnetic Reynolds number.

## 6 Magnetic Prandtl number and growth rate

In Section 5 we examined the growth rate and the saturation level of a fast dynamo when the magnetic Prandtl number is equal to one. Equally important is the study of the degree of their dependence when $`\mathrm{Pr}_\mathrm{m}`$ is varied. Thus, additional experiments were performed where the forcing again is such that it maintains the kinetic energy at a value close to the initial one, apart from small oscillations.

Not surprisingly, the growth rates during the laminar phase ($`t<30`$) of the flow are found to be the same (Figures 11, 12). By construction, the velocity is exactly the same in the laminar phase, independent of Re, and the effect of increasing $`\mathrm{Re}_\mathrm{m}`$ is simply, as in the comparison of the $`\mathrm{Re}_\mathrm{m}=200`$ and $`\mathrm{Re}_\mathrm{m}=400`$ experiments, to give a different transient at the beginning of the simulation.

Less trivially, the average growth rates in the turbulent regime also turn out to be similar with no crucial dependence on the magnetic Prandtl number for high $`\mathrm{Re}`$. This is illustrated in Figure 12, where the difference in the magnetic energy during the turbulent phase ($`t>30`$) is very small between the two experiments with $`\mathrm{Re}=200`$ and $`\mathrm{Re}=300`$. This may also be seen as a consequence of the generic properties of high Reynolds number turbulence. More specifically, by one of Kolmogorov’s main assumptions, energy dissipation (and hence the forcing) is independent of the fluid viscosity once the Reynolds number is high enough. Thus, and by analogy with the previous experiments, we expect that the only important change in the magnetic field is the development of more fine structure with increasing magnetic Reynolds number. The growth rate is determined by the large scale properties of the flow, and in particular by the turn over time, which by construction does not change between these experiments.

Finally, the above experiments reinforce the conclusion from Section 5, that the saturation level of the dynamo does not depend crucially on the value of $`\mathrm{Pr}_\mathrm{m}`$, but only increases slightly with $`\mathrm{Re}_\mathrm{m}`$ as it is shown in Figure 11.

## 7 Nature of the dynamo

The results of the above non-linear experiments are similar to the results of previous studies of turbulent hydromagnetic convection (Nordlund et al. Nordlund+ea92 (1992)) and both results indicate that during the dynamo action there is a close balance between the work done against the Lorentz force and Joule dissipation with a relatively small positive (negative) difference when the magnetic energy increases (decreases). In other words, work and dissipation are in near balance, with a relatively small difference responsible for the average growth of the magnetic energy. One may rightly ask why this is so; is the small positive difference then just fortuitous, and might the difference not just as well be negative for a slightly different flow?

To address this question we make use of a kinematic dynamo experiment, where the analysis is simpler than in a turbulent case. Turbulence makes the analysis of the Lorentz work difficult—there is much time dependent noise, which corresponds to waves averaging out over time. The work is much less noisy when the velocity is prescribed and/or has not yet entered the turbulent phase, and we therefore analyze the results of an experiment that deals with kinematic dynamo action. The velocity field is chosen to be the $`A=B=C=1`$ flow with $`k=1`$ and the magnetic Reynolds number is $`\mathrm{Re}_\mathrm{m}=300`$.

The convergence and the diffusion rates at the points of strong flux concentration are of the order of $`\sqrt{2}`$, which is much larger than the dynamo growth rate (which is about 20 times smaller, e.g. see Childress & Gilbert STF (1995)), so the advection and diffusion are in close balance in the neighborhood of the “flux cigars”. During the exponential growth of the magnetic energy the size of the flux cigars does not change dramatically, but the magnetic energy increases gradually.

A detailed analysis for the latter case confirms that this is indeed the case, and that the bulk part of the Joule dissipation and the work done against the Lorentz force are in detailed balance. This is understandable as a balance between kinetic energy being converted into magnetic energy, that is then immediately converted into heat through Joule dissipation.

Significant net positive work (defined by the difference $`\mathrm{W}_\mathrm{L}\mathrm{Q}_\mathrm{J}`$) occurs instead in regions where the magnetic field is weak, and the dissipation is very weak as well. In these regions, which cover 90% of the volume in this example, there is almost “pure work”; the Joule dissipation is much smaller than the work, and the net work is almost identical to the net work in the whole box.

In other words: The dissipation is strong only in the neighborhood of the $`\alpha `$-type stagnation points, but there it is balanced almost exactly by advection/convergence.

Table (1) shows the estimated values of the average work and dissipation over the $`89\%`$ of the computational volume. On the one hand, the dissipation is very weak (less than $`0.12\%`$ of the maximum dissipation level) and on the other hand, almost all the net work occurs in those $`89\%`$ of the volume. Moreover, the largest amount of dissipation and work occur in the remaining $`11\%`$ of the volume and they are in almost perfect detailed balance.

The above results lead to the conclusion that dynamo action occurs primarily in regions where the field is weak. The net work is much higher than the dissipation and as a result the field is amplified by dynamo action. Most of the net work occurs when the weak field is bent and stretched. The magnetic energy is then distributed over the volume through the Poynting flux and the field lines are dragged out and pile up against the local flux cigars where balance occurs between stretching and diffusion. The weak field is almost perfectly advected with the fluid and this advection of the weak field would be possible even in the saturated regime where the flow still has a good grip on the weakest field, although the strongest field is now exerting a strong resistance against motion/stretching, through the Lorentz force.

The above picture of the dynamo process is most likely related to the observation that the growth rates found for the kinematic ABC dynamos are always consistent with the large scale turn over times, even though the ratio of the magnetic energy to the rate of Joule dissipation ($`\mathrm{E}_{\mathrm{mag}}/\mathrm{Q}_\mathrm{J}`$), which is the time over which the magnetic energy is destroyed and replenished, is much shorter.

In conclusion, the work done on the weak part of the field is responsible for the increase in the magnetic energy by bending and stretching the magnetic field lines. The exponential amplification dominates the Joule dissipation which is weak through almost all of the volume. The detailed balance between the Lorentz work and the dissipation occurs only around the points where strong magnetic field structures are formed. Thus, a dynamo appears to be a process where amplification of the weak field plays the dominant rôle in the growth of the magnetic energy.

## 8 Conclusions

In this paper turbulent dynamo action has been studied and some new aspects of kinematic dynamo action have been illustrated. It was found to be important to keep the kinetic energy level of the same order through the simulation, by automatically adjusting the amplitude of the driving. Using relatively high values of $`\mathrm{Re}_\mathrm{m}`$ and $`\mathrm{Re}`$ we find that the growth rate of the magnetic energy amplification apparently depends on the stretching ability of the flow rather than on the values of the above parameters. This is obvious even in the turbulent phase, for different values of the magnetic Prandtl number, where the dissipation is increased. The level at which the magnetic energy saturates increases slightly with $`\mathrm{Re}_\mathrm{m}`$ but the transient in the saturated regime is expected to become insignificant at sufficiently high $`\mathrm{Re}_\mathrm{m}`$.

An important result is the existence of a near balance between the Lorentz work and the Joule dissipation, which is apparent in both the laminar and the turbulent regime of the forced ABC flow. Remarkably, the balance originates primarily from small regions where strong magnetic flux structures are concentrated. The nature of the dynamo is such that the net growth of magnetic energy comes about through stretching and folding of the weak magnetic field in the rest of the volume, where dissipation is very weak.

The aim of this paper was not to study dynamo action on a global scale (including aspects such as differential rotation, stratification and convection) but rather to study the different stages of dynamos produced by a family of hydrodynamically unstable flows. However, the physical processes and the dynamo mechanisms studied here are likely to be of similar importance in more realistic stellar dynamo settings.

###### Acknowledgements.

VA thanks the EU-TMR for support through a Marie Curie Fellowship. SBFD was supported by the Danish Natural Science Research Council. The work of ÅN was supported in part by the Danish National Research Foundation, through its establishment of the Theoretical Astrophysics Center. Access to computational resources granted by the Danish Center for Scientific Computing is gratefully acknowledged.