text:
# Modeling eruptions of Karymsky volcano

## I Introduction

A wide variety of types of volcanic activity exists: from devastating explosions separated by calm periods of many years or even centuries to a steady continuous outpouring of magma. In this paper we study mechanisms that give rise to a Strombolian-type volcanic activity, which is somewhere between these two extremes and characterized by regularly repeated explosions (10-400 per day) followed by relatively calm periods. Often in the course of a longer explosion, the gas and ash emission exhibits audible and visible modulations with a rather robust period of about 1 second. These modulations, often found on the seismograms of Strombolian-type eruptions, received a special name, “chugging” (Benoit and McNutt, 1997), because they resemble a periodic noise produced by a steam engine. An example of a volcano exhibiting such activity, we consider Karymsky volcano on Kamchatka peninsula in the Far East of Russia.

We believe that these Strombolian-type eruption patterns are caused by peculiarities of the motion of the magma column, and propose a model that describes the magma motion in a volcanic conduit as a creep flow of viscoelastic compressible medium with shear-stress-dependent boundary conditions. We show that the dynamics of the model indeed exhibits 2 levels of quasiperiodic behavior and resembles the eruption patterns of the Strombolian-type volcanos both on long and short time scales. The physical transparency of the suggested model allows us to express the observed time scales using the material properties of the magma.

Our paper is organized in the following way: first, we briefly describe geological aspects of the current eruption of Karymsky volcano. Then, after referring to the existing explanations of Strombolian-type activity, the model is formally introduced. We derive simple analytical expressions for a dormant time, a duration of explosion , and a typical period of tremor, and then present a numerical solution to a system of dynamical equations describing the motion of the magma column. The paper is concluded with a discussion of results and possible directions for further studies.

## II The Karymsky Eruption of 1996 - 2000

Karymsky volcano (Fig. 1), located in the central part of the East Kamchatka volcanic belt, is one of the most active volcanoes in the far-eastern Russia. It is a typical andesite stratovolcano composed of a $`77007800`$-year old caldera $`5`$ km in diameter (Ivanov et al., 1991) and a more recent cone which started growing about $`5300`$ years ago and at present (September 2000) is $`700`$ m high. The cone is composed of accumulated lava and pyroclastic materials, and the summit of the cone is capped by a double crater. The active part of the crater, formed during the most recent eruption, widened from 90 m in 1996 to 190 m in 2000 (Fig. 2). The absolute height of the volcano is 1549 m.

Karymsky activity usually consists of long effusive-explosive eruptions, with the eruption previous to the one being studied taking place in 1970-1982. The current eruption at Karymsky volcano began in January 1996 and apparently entered a lull phase in December 2000. Each eruption consists mostly of the discrete quasiperiodic gas and ash bursts sending plumes 100-1000 m above the crater rim (Fig. 3).

Ash and steam clouds have been observed extending downwind for 10-50 km, rarely 100-200 km.

A lava flow field, produced by the very viscous magma, runs down the south-western slope of the cone and reaches the length of $`1.4`$ km and the width $`100200`$ m. Even during the periods of highest activity, the motion of the lava flows was very slow (few meters per day) and was noticeable only by periodic rock-falls, happening every 5-15 minutes.

During particularly strong explosions, the eruptive column collapsed causing pyroclastic flows. In the lower part of the eruptive column, hot volcanic bombs of the sizes up to 2, rarely 5 m, were frequently observed. The bombs were dense and often covered with a bread-crust-like structure formed during cooling. On the surface of some bombs, noticeable extrusive tracks were found; they were scratched on the very viscous moving magma by the protruding parts of the magma channel.

Total volumes of pyroclastic and effusive materials deposited on the cone of Karymsky during the eruption of 1996-2000 is estimated to be 0.0274 and 0.0229 km3, respectively. Eruptive products, which include magma, volcanic bombs, and ash, predominantly consist of andesite, usually of black and dark grey color. The most common mineral in the eruptive products is plagioclase (20-25 %). Crystals of plagioclase can reach a size of 3-5 mm. The concentration of olivine and pyroxene is less than 1%, their crystal sizes are usually 2-3 mm. A chemical content of the andesite, averaged over 6 samples from the 1996 - 2000 eruption, is listed below: $`SiO_2`$ – 61.54, $`TiO_2`$ – 0.86, $`Al_2O_3`$ – 16.67, $`Fe_2O_3`$ – 2.43, $`FeO`$ – 5.12, $`MnO`$ – 0.11, $`MgO`$ – 1.99, $`CaO`$ – 5.31, $`Na_2O`$ – 3.69, $`K_2O`$ – 1.58, $`P_2O_5`$ – 0.25.

The feeding system of Karymsky volcano was studied by various methods, including gravity surveys, aeromagnetic surveys, photogrammetry, seismology and geodesy (Zubin, 1971; Maguskin, 1982; Shirokov, 1988). Despite slight discrepancies in the estimates of the depth of the upper border of magma chamber and its size, these data agree that underneath the Karymsky volcano in the close proximity to the surface there exists a spherical magma chamber. The upper border of this chamber is in few kilometers below the sea level. The diameter of the magma chamber is estimated to be between 1.5 and 7 km. A magma channel of a diameter of 100-200 m leads from the magma chamber to the top crater.

## III dual periodicity of Karymsky eruption

A distinctive feature of Karymsky eruptions is its rather robust periodicity which is observed on two timescales. An explosion, and a following quiet period, define the first period of activity which varied between 3 and 20 minutes and sometimes extending up to an hour. An explosion itself can start either abruptly or gradually (Fig. 4) and proceed according to one of the following scenarios:

during shorter explosions, the activity decays quickly and monotonously (in less than 10 seconds), while the longer explosions ($`2040`$ s, rarely 1 minute) may contain periods of relatively high and low activity. Often during longer ($`>20`$s) explosions, the intensity of the acoustic and seismic signals emitted by the volcano exhibits chugging, i.e. modulations with a typical period of $`1`$ second. The chugging defines the second, shorter, period in the cyclic activity of the volcano. A typical chugging event consists of 10-20 cycles (Fig. 4). The explosion is followed by a quiet period characterized by a complete lack of any activity in the crater and usually lasting much longer than the explosions themselves. Similar temporal patterns have been also observed during previous Karymsky eruptions. (Tokarev and Firstov, 1967; Farberov et al., 1983).

The systematic seismic and acoustic studies of Karymsky volcano were conducted in 1997-1999 by 3 Russian-US expeditions. The seismograms presented in Fig. 4 were recorded in 1998 by a three-component wide-band seismometers CMG 40-T, but the most detailed information about a fine structure of the chugging signal was obtained using a set of infrasound microphones. For a more complete description of the seismic and infrasound recordings at Karymsky the reader is referred to (Johnson et al., 1998; Johnson, 2000; Johnson and Lees, 2000).

## IV Existing models of eruption periodicity

Currently there is no single well-established view on the mechanism causing periodic patterns in the eruptions of Strombolian type. Below we mention several possible explanations of origins of the periodicity. Lees and Bolton (1998) compared the eruption process to the steam exhaust from a hydrodynamic system with a nonlinearly controlled exhaust valve resembling a pressure cooker. Manga (1996) suggested that the periodicity is caused by a segregation of ascending bubbles into waves, as in a pint of Guiness beer. Julian (1994) proposed a mechanism of nonlinear hydrodynamic oscillations in the ascending magma flow. Jaupart and Vergniolle (1989) suggested that oscillations are caused by periodical collapse of a bubble foam trapped in pockets of the magma chamber. Recently, there appeared a paper by Delinger and Hoblit (1999) where an ad hoc hysteresis model was developed to describe a periodic behavior of silicic volcanos. However, self-consistent explanation for both long- and short-time periodicity and extended calm periods in the activity of Strombolian-type volcanoes is still missing. Our model, based on viscoelastic hydrodynamical description of the motion of the magma column with shear stress-dependent boundary conditions, accounts for both timescales.

## V The model

Let us first qualitatively discuss the motion of magma in the conduit. Since the temperature and pressure in the magma column decrease with the height, the viscosity of magma increases rapidly in the upper part of the conduit. Inspecting eruption products ejected to the surface, it is natural to deduce that at the very top of the column the magma is almost solid. However, deeper parts of the magma column are still hot and much less viscous. Given this significant difference in rheological properties of the upper and lower sections of the magma column, we consider the moving magma as consisting of two distinct parts: long lower part (LP) and short upper part (UP). We assume that the viscosity of the magma in the LP is negligible. Due to the large amount of dissolved gases, magma in the LP is compressible. We also assume that there is a constant supply of fresh magma to the bottom of the LP. Unlike the LP, magma in the UP is cold and viscous; in addition to that, as any media near the liquid-solid transition point, it possesses a certain degree of elasticity. Because of the small length of the UP compared to the LP ($`100`$ m vs. $`5`$ km), we can neglect compressibility of the UP. In summary, we consider a long cylindrical tube filled with a non-viscous compressible media, fed to the bottom with a constant velocity and a short viscoelastic plug (UP) on the top. The boundary conditions for the moving UP are controlled by a local value of a shear stress in the vicinity of the cylinder wall. When the boundary shear stress is low, the plug “sticks” to the cylinder wall and the boundary velocity is zero. As the shear stress increases, the magma near the wall undergoes phase transition (melts) to a much less viscous liquid state, and the boundary layer of the plug slips along the walls of the channel. This shear-induced phase transition happens when the viscoelastic state with the non-zero shear modulus becomes thermodynamically less stable than a pure liquid state with zero shear modulus. Phenomenologically, the transition between stick and slip boundary conditions for the very viscous body is similar to the effect of dry friction, where a motion begins only when a driving force exceeds some threshold value (Persson, 1999, Ch. 8).

Qualitatively, a cycle of the system evolution can be described as follows: when the pressure in the compressible LP of the column is low, the shear stress $`\sigma `$ in the plug is much smaller than the threshold phase transition value $`\sigma _{slip}`$, and the velocity of the UP is much less than the feeding velocity $`V_0`$. The difference in feeding and plug velocities results in the growth of the pressure of the compressible magma in LP which, in its turn, increases the shear stress in the plug. This stage corresponds to the dormant period between explosions. When the shear stress in the boundary layer of the UP exceeds the critical value $`\sigma _{slip}`$, the stick-slip phase transition occurs and the plug begins to move. This corresponds to the appearance of visible signs of a volcanic explosion.

Depending on the parameters of the model, the further evolution of the system proceeds according to one of the following two scenarios. For certain conditions, a relaxation of the accumulated in the UP viscoelastic deformations can result in a low, or even reverse, velocity of the boundary layers of the UP relative to the walls of the conduit. If these conditions persist for sufficient time, they can cause a drop of the shear stress below the slip-stick transition critical value $`\sigma _{stick}`$, which in turn can give rise to the reverse, slip-stick transition. However, the stick state will be very short-lived, since the accumulated excessive pressure in the LP and oscillatory motion of UP will result in the immediate increase of the boundary shear stress and a new stick-slip transition. Length of a period of such coupled viscoelastic-stick-slip process is controlled by the density, the elastic modulus, and the characteristic time of the phase transition. In the second scenario, sticking does not occur during short-periodic viscoelastic oscillations of the UP, and the period of these oscillations is roughly equal to that of a freely oscillating elastic membrane with mixed boundary conditions. Either of these oscillation processes can cause the short-period modulation of the activity of the volcano.

During the active phase, the velocity of the UP is much higher than the feeding velocity $`V_0`$ which results in an expansion of the LP and a decrease of pressure under the UP. When the pressure under the UP drops so that the shear stress in the UP is below the slip-stick threshold value $`\sigma _{stick}`$, a “long-term” sticking of the plug occurs, and the system enters a new quiet period of evolution. Because of the hysteresis ($`\sigma _{slip}<\sigma _{stick}`$) associated with the first order nature of the stick-slip transition, it takes some time to build enough pressure for a new explosion to begin. This cycle corresponds to the second, longer period (10-30 minutes) in the dynamics of eruption.

Although simple in nature, this model qualitatively explains both scales of periodicity of Karymsky volcano and does not contradict to the observations. Below a formal analysis of the dynamics of the model is presented.

Let us consider a cylinder (LP) of the radius $`R`$ and the length $`L`$, filled with an ideal non-viscous compressible medium feeded from the bottom with a constant velocity $`V_0`$. The top of LP ends with a short viscoelastic plug (UP) of length $`l`$, $`lL`$ (Fig. 5).

The velocity field in the UP, of density $`\rho `$ and driven by pressure $`P`$ and stress tensor $`\sigma _{ij}`$, is described by the Navier-Stokes equation,

$$\rho \frac{v_j}{t}=\frac{P}{x_j}+\frac{\sigma _{ij}}{x_i}.$$
(1)
Because of the high viscosity and the low velocity of magma in the UP, the nonlinear (convective) term is omitted. To account for the viscoelastic properties of the UP we use the simple linear Maxwell model (see, for example, Shore et al.), with a characteristic memory time $`\tau _m`$,

$$\tau _m\frac{\sigma _{ij}}{t}=\sigma _{ij}+\eta e_{ij}.$$
(2)
Here $`\eta `$ is the dynamic viscosity, and $`e_{ij}v_i/x_j+v_j/x_i`$ is the shear rate tensor. Incompressibility condition for the magma in the UP reads:

$$\frac{v_i}{x_i}=0.$$
(3)
For simplicity, we assume that the motion of the plug is axially symmetric and select the $`Z`$ direction of our cylindrical coordinate system parallel the axis of the cylinder directed along the magma flow. Also, we also consider the pressure gradient constant through the UP and being equal to $`\frac{P_0}{l}`$, where $`P_0`$ is a pressure at the bottom of the UP. Using (2,3), the Navier-Stokes equation (1) for $`V_z`$ ($`Z`$-component of the velocity of the plug averaged over the plug length) can be transformed to:

$$\rho \frac{V_z(r,t)}{t}=\frac{P_0(t)}{l}+\frac{\eta }{\tau _m}^t\mathrm{exp}\left(\frac{tt^{}}{\tau _m}\right)\mathrm{\Delta }_rV_z(r,t^{})𝑑t^{}.$$
(4)
Here $`\mathrm{\Delta }_r^2/r^2+1/r/r`$ is the radial part of the Laplacian operator in the cylindrical coordinates.

From the conservation of mass of the compressible media in the LP, it follows for the pressure $`P_0`$ on the bottom surface of the UP,

$$V(t)V_0=\beta L\frac{dP_0(t)}{dt},V(t)2_0^RV_z(r,t)r𝑑r/R^2,$$
(5)
where $`\beta 1/v(v/p)_T`$ is the compressibility (isothermal) of the media in the LP, and $`V(t)`$ is the velocity of the UP $`V_z(r,t)`$, averaged over the perpendicular crossection.

To complete the description of the dynamics of the system, the equations (4,5) must be complemented by boundary conditions for the velocity on the conduit walls, $`V_z(r=R,t)`$. These boundary conditions are controlled by the kinetics of the first-order phase transition driven by the value of the shear stress $`\overline{\sigma }(t)\sigma _{rz}(t,R)`$ near the wall and are written in a general mixed form:

$$V_z(R,t)=R\psi (t)\frac{V_z(r,t)}{r}|_{r=R}.$$
(6)
The quantity $`R\psi (t)`$ is usually called a “slipping length”. It shows at what distance inside the wall the linearly extrapolated velocity becomes equal to zero. For $`\psi (t)=0`$ the boundary conditions for the velocity are “stick”, i.e. $`V_z(R,t)=0`$; for $`\psi (t)>0`$ the boundary layers of the UP “slip” along the walls of the cylinder with some finite velocity, $`V_z(R,t)>0`$.

The time evolution of the dimensionless slipping length $`\psi (t)`$ is determined from the time-dependent Ginzburg-Landau equation, which describes relaxation of a system with a first-order phase transition:

$$\tau _{GL}\frac{d\psi }{dt}=\frac{dF}{d\psi },$$
(7)
where $`\tau _{GL}`$ sets a phase transition timescale. The free energy $`F(\psi )`$ is a usual double-well potential with a tilt depending on the value of the shear stress $`\overline{\sigma }(t)`$ on the boundary of the UP. We choose a double-parabola potential which has fixed positions of both minima, $`\psi =0`$ and $`\psi =\psi _0`$:

$$F(\psi )=\{\begin{array}{ccc}\psi ^2,& \psi <\psi ^{}& \\ (\psi \psi _0)^2h,& \psi \psi ^{}& \end{array}$$
(8)
Here $`h`$ is a tilt parameter; matching point $`\psi ^{}`$ is determined from continuity requirement, $`\psi ^{}=(\psi _0^2h)/2\psi _0`$. To determine $`h`$ and to illustrate how Eqs. (7,8) work, let us consider a typical stick-slip cycle. We go around a hysteresis curve shown in Fig. 6.

When the shear stress $`\overline{\sigma }(t)`$ is small and growing (line I), the boundary condition is “stick” and stationary, $`\psi =\frac{dF}{d\psi }=0`$. As the shear stress exceeds the critical value $`\sigma _{slip}`$, $`\psi (t)`$ starts moving towards $`\psi _0`$ with the rate $`2\psi _0/\tau _{GL}`$ (line II). After $`p`$ reaches $`\psi _0`$, the boundary condition becomes “slip” and stationary again, and the shear stress starts decaying (line III). After $`\overline{\sigma }(t)`$ drops below $`\sigma _{stick}`$, $`\psi `$ relaxes to 0 (line IV). The state I corresponds to the dormant period, the state III corresponds to the explosion without intermediate sticking; the transition intervals II, IV are usually very short. The free energy tilt parameter $`h`$ is linearly connected to $`\overline{\sigma }`$:

$$h=\psi _0^2\left[\frac{2\overline{\sigma }}{\sigma _{slip}\sigma _{stick}}\frac{\sigma _{slip}+\sigma _{stick}}{\sigma _{slip}\sigma _{stick}}\right].$$
(9)
Similarly to the Eq. (4), in the framework of the Maxwell model the shear stress $`\overline{\sigma }(t)`$ is expressed as a convolution of a shear rate $`\frac{V_z(r,t^{})}{r}`$ with an exponential memory function $`\mathrm{exp}(\frac{tt^{}}{\tau _m})`$

$$\overline{\sigma }(t)=\frac{\eta }{\tau _m}^t\mathrm{exp}(\frac{tt^{}}{\tau _m})\frac{V_z(r,t^{})}{r}|_{r=R}dt^{}$$
(10)
The equations (410) completely define the dynamics of our model. A very similar set of equation was derived by Shore et al. (1997) for the physically equivalent model of polymer extrusion. Because of the strong non-linearity of the Eq. (7), an analytic solution of these equations is impossible. We refer readers to Shore et al., 1997, for a detailed description of a linear stability analysis of this system. Contrary to a stable flow, a dynamical regime corresponding to volcanic activity is strongly unstable and nonlinear. However, large separation of characteristic timescales allows us to give reasonable analytical estimates of important temporal features of the model.

Let us first look at the shortest timescale which corresponds to small oscillations of the UP when it slips along the conduit walls without intermediate sticking. For this timescale the pressure term in (4) is considered to be constant. Dynamics of the UP is physically equivalent to the oscillation of an elastic circle membrane with the mixed boundary conditions (6). Only the lowest harmonic is considered,

$$V_z(r,t)=V_0e^{i\omega t}J_0\left(\frac{r}{R}\lambda \right),$$
(11)
where $`\lambda `$, a number of order of one, depends on the slipping length $`\psi `$: $`J_0(\lambda )=J_1(\lambda )\lambda \psi `$. Here $`J_0`$ and $`J_1`$ are the zero- and first-order Bessel functions. After plugging Eq. (11) into Eq. (4) and disregarding the constant pressure term, we obtain the following dispersion relation:

$$\omega =\frac{i\pm \sqrt{4\frac{\lambda ^2\eta \tau _m}{R^2\rho }1}}{2\tau _m}.$$
(12)
This yields for a period $`T_{osc}`$ of not too overdamped UP oscillations ($`T_{osc}\tau _m`$):

$$T_{osc}\frac{2\pi R}{\lambda }\sqrt{\frac{\rho \tau _m}{\eta }}.$$
(13)
Taking into account a commonly used relation between the Maxwell time and shear modulus $`G`$, $`\tau _m\eta /G`$ (see, for example, Webb (1997)), we observe that Eq. (13) reduces to a usual expressions for a period of a small free oscillations of an elastic membrane,

$$T_{osc}\frac{2\pi R}{\lambda }\sqrt{\frac{\rho }{G}}.$$
(14)
Now we go to much longer timescales, $`T\{T_{osc},\tau _m\}`$ and look first at the steady state UP boundary shear stress $`\overline{\sigma }=\frac{V_z(r)}{r}|_{r=R}`$. Since this timescale is much greater than the elastic memory time $`\tau _m`$, the viscoelastic term in (4) is relaxed to a usual viscous term, $`\eta \mathrm{\Delta }_rV_z(r,t)`$. For a steady state, the averaged over the crossection of the conduit velocity $`V(t)`$ must be equal to $`V_0`$. For an arbitrary slipping length $`\psi `$, a straightforward calculation yields:

$$\overline{\sigma }=\frac{4V_0\eta }{R}\frac{1}{1+4\psi }.$$
(15)
It immediately follows that in order for stick-slip repeating cycles to happen, the “stick” steady state shear stress ($`\overline{\sigma }`$ for $`\psi =0`$) should be greater than $`\sigma _{slip}`$ and “slip” steady state shear stress ($`\overline{\sigma }`$ for $`\psi =\psi _0`$) should be less than $`\sigma _{stick}`$. These conditions imply for the feeding velocity $`V_0`$:

$$\frac{\sigma _{slip}R}{4}<V_0<\frac{\sigma _{stick}R(1+4\psi _0)}{4}.$$
(16)
Besides the steady state properties, we can also evaluate characteristic times between the stick-slip and the slip-stick phase transitions, i.e. duration of the dormant and explosive periods. Assuming that these times are much longer than both $`T_{osc}`$ and $`\tau _m`$, we again disregard all memory effects in the viscous term in Eq. (4) and consider the flow inertialess by omitting the left-hand-side term $`\frac{V_z(r,t)}{t}`$. After averaging the remaining part of Eq. (4) over the crossection of the flow and plugging it in (5), we obtain for arbitrary slipping length $`\psi `$:

$$^t\frac{V_oV(t^{})}{\beta Ll}𝑑t^{}\frac{8\eta V(t)}{R^2(1+4\psi )}=0.$$
(17)
Combining the solution of this equation,

$$V(t)=V_o+Ce^{\frac{t}{t_0}},t_0\frac{8\eta \beta Ll}{R^2(1+4\psi )},$$
(18)
with expressions for stick and slip velocity that follow from (15), we obtain for dormant ($`T_d`$) and explosive ($`T_e`$) times:

$$T_d=\frac{8\eta \beta Ll}{R^2}\mathrm{ln}\frac{4\eta V_0\sigma _{stick}R}{4\eta V_0\sigma _{slip}R}$$
(19)
$$T_e=\frac{8\eta \beta Ll}{R^2(1+4\psi _0)}\mathrm{ln}\frac{\sigma _{slip}R(1+4\psi _0)4\eta V_0}{\sigma _{stick}R(1+4\psi _0)4\eta V_0}.$$
(20)
## VI NUMERICAL RESULTS

To check our analytical predictions and get an overall view of the system dynamics, we solved the Eqs. (410) numerically. We used a simple Eulerian finite-difference scheme with a uniform radial grid of 50-100 points. All the integrals are approximated by trapezoid formula. Given that time step is sufficiently small to avoid von Neumann-type instabilities, this simplest difference scheme proved to be sufficiently robust.

We chose the numerical parameters to be in general accordance with the literature (Murase and McBirney, 1973;, Borgia and Linneman, 1990; Griffiths and Fink, 1993;, Bagdassarov and Dorfman, 1994; Webb, 1997; Bagdassarov and Dorfman, 1998) and with our estimates (14,16,19,20). The following values were used:

* Density $`\rho =2400kg/m^3`$
* Dynamic viscosity $`\eta =10^9Pas`$
* Conduit raidus $`R=50m`$
* UP length $`l=100m`$
* LP length $`L=5km`$
* Compressibility $`\beta =\mathrm{4\hspace{0.17em}10}^{10}Pa^1`$
* Maxwell time $`\tau _m=0.4s`$
* Ginzburg-Landau time $`\tau _{GL}=10^3s`$
* Average velocity $`V_0=2\times 10^3m/s`$
* Slipping shear stress $`\sigma _{slip}=10^5Pa`$
* Sticking shear stress $`\sigma _{stick}=\mathrm{3\hspace{0.17em}10}^4Pa`$
* Dimensionless slipping length $`\psi =10`$

With this choice of numerical parameters a no-stick explosion takes place, during which the UP moves along the walls of the conduit with the constant slipping length $`\psi R`$ until it finally sticks to the wall at the end of the explosion. For such explosions a choice of phase transition time $`\tau _{GL}`$ does not affect the dynamics, given that it is shorter than explosion time, $`\tau _{GL}T_e`$.

Plots of averaged over crossection velocity $`V(t)`$ and pressure gradient vs. time are presented in Fig. 7. The simulation results for dormant and explosion times, as well for short-term oscillation period, $`T_d^s530`$ s, $`T_e^s21`$s, $`T_{osc}^s0.69`$ s are in good agreement with our simple analytical estimates, $`T_d=497`$ s, $`T_e=20`$ s, and $`T_{osc}=0.69`$ s. To calculate $`T_{osc}`$ we have to solve the equation $`J_0(\lambda )=J_1(\lambda )\lambda \psi `$, which for $`\psi =10`$ yields $`\lambda 0.43`$. For the stick boundary conditions ($`\psi =0`$), $`\lambda 2.3`$; which gives for the period of oscillations at the end of explosion cycle when UP sticks to the walls $`T_{osc}^{}0.13`$, which is again in very good agreement with the simulation data.

Another possible scenario for an eruption is when during explosion cycle the UP sticks to the walls numerous times. It happenes when the shear stress at the wall, modulated by the free membrane-like oscillations of the UP, drops below the sticking value $`\sigma _{stick}`$, and the slip-stick phase transition is sufficiently fast, $`\tau _{GL}T_{osc}`$. To simulate this regime, we increase the value of $`\sigma _{stick}`$ from $`\mathrm{3\hspace{0.17em}10}^4`$ to $`\mathrm{8.5\hspace{0.17em}10}^4`$ Pa and keep all other parameters constant. The long-time dynamics, characterised by the times $`T_e`$ and $`T_d`$ only slightly changes since both of these times depend on $`\sigma _{stick}`$ logarithmically. However, a short-time dynamics, previously described by decaying free oscillations with a period $`T_{osc}`$, changes radically and becomes chaotic (Fig. 8).

A typical time separation between consecutive velocity maxima ($`0.3`$ s for the first three maxima in the bottom part of Fig. 8) becomes somewhat smaller than free oscillation period $`T_{osc}0.7`$. Qualitatively, it happenes when during an oscillation the velocity increases so that shear stress drops below $`\sigma _{stick}`$, a sticking occurs that results in a rapid decrease of velocity, so the oscillation is “truncated” before it reaches its maximum. Another distinct feature of this regime is that the short-time chaotic modulation of the UP velocity continues through the whole explosion phase, while in the non-sticking regime, the oscilations decay after a few periods.

## VII Discussion

In the previous sections the long-term behavior and two possible scenarios of the short-time activity of the model based on viscoelastic hydrodynamics and stick-slip transition were described. Patterns in the time dependence of the UP velocity, presented in Figs. 7,8, qualitatively match the seismic and acoustic signals (Fig. 4) recorded at Karymsky volcano. Yet nothing has been said about a mechanism describing generation of seismic and acoustic waves themselves. The oscillating UP and related pressure oscillations in the LP can cause seismic waves either directly or act as a “pacemaker” for other processes. These processes can include a stimulated emission of gases either through creation and change of geometry of cracks in the UP or via degassing of magma under action of standing and running compression waves. Since these mechanisms have a strongly nonlinear character, their coupling to the UP oscillation can modify the oscillation period and decay patterns. We leave consideration of these phenomena for future study.

## VIII conclusion

As a result of the visual, seismic, and acoustic observations at Karymsky volcano, two scales of temporal periodicity were observed in the eruptive dynamics. First, the repetitive explosions which produced tha ash and gas clouds and volcanic bombs in 10-20 minutes intervals, and secondly, the periodic modulation of eruptive activity (chugging) with a typical timescale of order of 1 s. The existence of these two scales of temporal periodicity are not limited to the 1996-2000 eruption of Karymsky volcano, similar timescales were observed during previous Karymsky eruptions and on other andesite volcanos in the world.

A model of the motion of the magma column was suggested that accounts for both timescales in the eruption of Karymsky and other similar andesite Strombolian-type volcanos. Low- and high-frequency periodicity is caused by two distinct ways of accumulation of the elastic energy. The viscoelastic properties of the UP causes the high frequency tremor, while the compressibility of magma in the LP explains the low frequency cycles. The shear-stress-dependent stick-slip phase transition in the UP introduces the hysteresis into the dynamics of the magma column motion and explains the unharmonicity of the oscillations and long dormant periods in the eruption dynamics Physical simplicity of the model allows analytic estimates for the long and short periods of the explosive activity to be obtained.

## IX Acknowledgments

The authors are grateful to Martin Grant for stimulating discussions on the polymer extrusion problem, to Julie M. Kneller for carefully reading the manuscript, and especially to Emily Brodsky for providing important references, numerical values for magma parameters, and carefully reading the manuscript. This work was partially supported by the U.S. National Science Foundation.

## X Bibliography

* Bagdassarov N. S., Dingwell D. B. and Webb S. L., 1994. Viscoelasticity of crystal- and bubble-bearing rhyolite melts, Phys. Earth Planet. Inter., 83: 83-99.
* Bagdassarov N. and A. Dorfman A., 1998. Viscoelastic behavior of partially molten granites, Tectonophysics, 290: 27-45.
* Benoit J. P., McNutt S. R., 1997. New constraints on source of volcanic tremor at Arsenal Volcano, Costa Rica, using broadband seismic data. Geoph. Res. Lett., 24: 449-452.
* Bird R. B., Armstrong R. C., and Hassager O., 1987. Dynamics of Polymeric Liquids. Vol.1: Fluid mechanics, Wiley, New York.
* Borgia A. and Linneman S. R., 1990. On the Mechanism of Lava Flow Emplacement and Volano Growth: Arenal, Costa Rica. Lava Flows and Domes, ed. J.H. Fink, Springer-Verlag.
* Delinger R. P. and Hoblitt R. P., 1999. Cyclic eruptive behavior of silicic volcanoes. Geology, 27: 459-462.
* Farberov A. I., Strocheus A. B., and Pribulov E. S., 1983. Studies of weak seismicity of Karymsky volcano, August 1978. Vulkanologija i Seismologija, 3: 78-89 (in Russian).
* Griffths R. W. and J. H. Fink, 1993. Effects of surface cooling on the spreading of lava flows and domes. J. Fluid Mech., 252: 667-702.
* Hess P. S., 1989. Origins of Igneous Rocks, Harvard Univ. Press., Cambridge, USA.
* Jaupart C. and Vergniolle S., 1989. The generation and collapse of a foam layer at the roof of a basaltic magma chamber. J. Fluid Mech., 203: 347-380.
* Johnson J. B., 2000. Interpretation of infrasound generated by erupting volcanoes and seismological energy partitioning during Strombolian Explosions. PhD Thesis, Univ. of Washington, Seattle.
* Johnson J. B. and Lees J. M., 2000. Plugs and Chugs - Strombolian activity at Karymsky, Russia, and Sangay, Ecuador. J. Volc. Geotherm. Res., 101: 67-82.
* Johnson J. B., Lees J. M., and Gordeev E. I., 1998. Degassing explosions at Karymsky volcano, Kamchatka. Geophys. Res. Lett., 25: 3999-4000.
* Julian B. R., 1994. Volcanic tremor: Nonlinear excitation by fluid flow. J. Geoph. Res., 99: 11859-11877.
* Lees J. M. and Bolton E. M., 1998. Pressure cookers as volcano analogues. EOS, Trans. Am. Geoph. Un., 79 (45), Fall Meeting Supp., F620.
* Magus’kin M. A., Enman V. B., Seleznev B. V., and Shkred V. I., 1982. Peculiarities of the Earth core displacement on Karymsky volcano in 1970-1981 from geodesical and photogrammetric data. Vulkanologija i Seismologija, 4: 49-64 (in Russian).
* Manga M. 1996. Waves of bubbles in basaltic magmas and lavas. J. Geoph. Res., 101: 17457-17465.
* Murase T. and McBirney A. R., 1973. Properties of some common igneous rocks and their melts at high temperatures, Geol. Soc. Am. Bull., 84: 3563-3592.
* Persson B. N. J. , 1999. Sliding Friction. Springer-Verlag, Ch. 8.
* Shirokov V. A., Ivanov V. V., Stepanov V. V., 1988. On the deep structure of Karymksy volcano and peculiarities of its seismicity: local seismic network data. Vulkanologija i Seismologija, 3: 71-80 (in Russian).
* Shore J. D. , Ronis D., Piche L., and Grant M., 1997. Theory of melt instabilities in the capillary flow of polymer melts. Phys. Rev. E, 55: 2976-2992.
* Tokarev P. I. and Firstov P. P., 1967. Seismological studies of Karymsky volcano. Bulluten’ Vulkanologicheskih Stanciy, 43: 9-22 (in Russian).
* Webb S., 1997. Silicate melts: Relaxation, Rheology, and the Glass transition. Rev. of Geophys., 35: 191-218.
* Zubin M. I, Ivanov B. V., Shteinberg G. S., 1971. Structure of Karymsky volcano, Kamchatka, and some aspects of caldera genesis. Geologiya i Geofizika, 1: 73-81 (in Russian).


format: text

text:
# How Mira Variables Change Visual Light by a Thousand-fold

## 1 Introduction

Mira variables have long fascinated both amateur and professional astronomers. In the late 1500s, European and Asian astronomers noted the “disappearance” of a star, most likely $`o`$ Ceti, over a period of months. While this star dimmed dramatically, it did not really disappear and, in the early 1600s, Fabricus and Holwarda recognized that $`o`$ Ceti was a variable star (Hoffleit, 1997). Because of the extreme changes of its brightness, $`o`$ Ceti was called Mira (or “marvelous one”). Indeed, Mira variables are marvelous; some regularly change visual luminosity by $`8`$ magnitudes, or factors of more than 1000, over periods of hundreds of days (Merrill, 1940).

Mira variables are evolved red giants with extraordinarily large sizes; a typical Mira has a radius of $`2`$ AU (e.g., Pease 1931; Labeyrie et al. 1977; Lopez et al. 1997; Tej et al. 1999). Mira variables are thought to pulsate, owing to opacity from the variable ionization of hydrogen and helium below the stellar surface (e.g., Keeley 1970). These pulsations emerge as shock waves and propagate upwards into the lower photosphere near maximum visible light (e.g., Slutz 1976). A Mira variable can have a very complex, dynamic atmosphere, and much theoretical and computational work has been expended attempting to understand the physical processes involved (e.g., Wood 1979; Bertschinger & Chavalier 1985; Bowen 1988; Bessel et al. 1989; Plez, Brett & Nordlund 1992; Höfner et al. 1998; Woitke et al. 1999; Helling, Winters & Sedlmayr 2000). Because Miras are so large, their surface temperatures are low and many molecules form in their atmospheres. At distances of a few stellar radii, temperatures fall even lower, which allows the formation of silicate or graphite “dust.” Once dust grains form, radiation pressure rapidly disperses them into interstellar space, enriching the environment with heavy elements. The ability of red giant stars to form and disperse circumstellar material is crucial to the life-cycle of stars and for life as we know it. Thus, the study of the workings of the cool atmospheres of such stars has great significance to astrophysics.

Nearly 70 year ago, Pettit & Nicholson (1933) showed that Mira variables have temperature variations of $`30`$% and total luminosity variations of factors of $`2`$ over their roughly yearly cycle. The luminosity, $`L`$, the effective temperature, $`T_{eff}`$, and the radius, $`R`$, of a star are related by $`L=\sigma T_{eff}^4\pi R^2`$, where $`\sigma `$ is the Stephan-Boltzmann constant. Using this relation, measurements of $`L`$ and $`T_{eff}`$ give an estimate of $`R`$, and radii of Miras determined in this manner indicate size changes of $`40`$%.

The observed visual (Mattei, 2001) and infrared light (Lockwood & Wing, 1971) curves of $`\chi `$ Cygni, a bright oxygen–rich Mira variable, are presented in Fig. 1 as green and red crosses, respectively. Plotted with solid lines are the expected flux density as a function of time for a time–variable black-body model based on changes in temperature and size as determined by Pettit & Nicholson (1933). The infrared model falls only about 1 magnitude dimmer than the data and has approximately the correct amplitude. Thus, by adjusting the stellar size, this model could fit the infrared light curve well. However, the visible curve (green) starts at maximum light about 2 magnitudes too bright, and the amplitude is about 4 magnitudes too small. This is a general problem for Mira variables: the measured temperature and radius changes, determined mostly from infrared data, predict only about 3 to 4 magnitudes of visual light variation, whereas up to 8 magnitudes are sometimes observed.

Good fits to the observed light curves cannot be obtained by simply changing the variation of the temperature and radius of a stellar black-body. If the black-body temperatures are allowed to decrease considerably below those indicated by infrared data, one can generate a very large visual amplitude and reproduce the visual magnitude at minimum light. For example, in Fig. 1 we show light curves (dashed lines) for a stellar black-body whose temperature varies from 2300 K at maximum light to 1300 K at minimum light. While the visual light curve can be successfully modeled in this manner, the model for the infrared light curve degrades considerably. The problem is that time varying black-body models cannot simultaneously fit the optical and infrared light curves of a star like $`\chi `$ Cygni.

Perhaps one could explain Mira light curves through the effects of a strong opacity source that varies considerably with wave-band. Astronomers have long suspected that molecular absorption in the atmospheres of Mira variables contributed to the large change in visual luminosity (e.g., Pettit & Nicholson 1933; Smak 1964; Celis 1978). However, to our knowledge, this problem has never been directly addressed. In this paper, we present an explanation of how the formation of metallic oxides (notably TiO), as the star periodically cools, can solve the problem of Mira light curves.

## 2 Basic Physical Explanation

The fact that Mira variables change infrared (and total) luminosity by only about 1 magnitude suggests that the star is not undergoing exotic changes. Where then should one look to explain the dramatic optical variations? Consider a star with an average temperature of 2200 K and a radius of 2 AU. At a distance of 100 parsecs, this star would have visual magnitude, $`m_v`$, of $`+4.5`$ mag. Varying the radius of such a star by $`\pm 20`$% (from 2.4 to 1.6 AU) changes $`m_v`$ by only $`\pm 0.4`$ mag. Clearly, size changes are insufficient to explain the large observed swings in visual light. On the other hand, varying the temperature of the star by only $`\pm 15`$% (from 2500 K to 1900 K) strongly affects visible light — changing it by 3.4 mag. While this is inadequate to explain the large changes in the visual light of many Miras, it does suggest that the explanation comes predominantly through the effects of temperature variations.

The emission from a Mira variable with an 8 magnitude amplitude could be explained if, at minimum light, the radiation comes from material at about 1400 K. How far above the stellar surface111Throughout this paper we use the terms stellar surface, stellar size, and stellar radius to mean those quantities measured at infrared wavelengths in line-free regions. does this temperature occur? If, for example, the temperature decreases with radius as $`r^{1/2}`$, then a star with a minimum surface temperature of 1900 K will experience temperatures near 1400 K at about 1.8 stellar radii ($`R_{}`$). Thus, the problem of explaining the large amplitude of Mira variables may be reduced to finding a source of visual opacity sufficient to make the star appear 80% larger than the stellar size at minimum light, while remaining close the stellar size at maximum light. In order to accomplish this, the opacity in the stellar atmosphere should have the following characteristics: 1) The opacity must increase significantly as the temperature decreases from about 1900 to 1400 K. Were, for example, the opacity to remain high throughout the stellar cycle, the amplitude of the variations would only approach about 5 magnitudes. 2) The visual opacity must be great enough to achieve optical depths of $`1`$ for the low densities expected in the stellar atmosphere at $`r2R_{}`$.

It has long been known that oxygen-rich (M– and S–type) Miras have strong metallic oxide absorption bands across the visible spectrum, with titanium oxide (TiO) typically providing the dominant absorption (Merrill, 1940). In the extended photosphere of an oxygen rich Mira variable, the chemical reaction of Ti + O $``$ TiO proceeds from predominantly atomic to molecular phases for temperatures between about 2000 and 1600 K for total densities of $`10^{10}`$ and $`10^{14}`$ gm cm-3, respectively. So as the star cools when approaching minimum light, titanium can transition from mostly atomic to molecular form in the extended atmosphere of a Mira variable. Once TiO forms, it is extremely efficient at absorbing optical photons. This is why some sun screen products use titanium (or other metallic) oxides.

The implications for a star that naturally forms sun screen in its atmosphere can be dramatic. The optical depth to molecular absorption can be written as follows: $`\tau =n_m\alpha 𝑑\mathrm{},`$ where $`n_m`$ is the molecular density, $`\alpha `$ is the absorption cross-section, and $`d\mathrm{}`$ is path length along our line-of-sight. Can one achieve $`\tau 1`$ in the extended photosphere owing to TiO absorption? For an approximate estimate of opacity, assume a characteristic path length of $`\mathrm{}10^{13}`$ cm and a TiO cross-section, averaged over the visual band, of $`\alpha 10^{17}\mathrm{cm}^2`$ (Jorgensen, 1994). Then unit optical depth requires a characteristic TiO number density of $`n_m=10^4\mathrm{cm}^3`$. For solar abundances, this corresponds to a total mass density of $`10^{13}`$ gm cm-3. Such a mass density is reasonable in the extended stellar photosphere of a Mira variable. Thus, near minimum light TiO (and other metallic oxides) can absorb most of the optical light and produce a “visual surface” far above the true stellar surface.

The change in the visual appearance of a Mira between maximum and minimum visible light is schematically represented in Fig. 2. The left-hand panel of the figure represents the star near maximum visual light ($`\varphi =0`$); one can see through the extended atmosphere, in which metals such as Ti are predominantly in atomic form, almost down to the stellar surface. (Actually, even though the fraction of Ti in molecular form is very low at temperatures above 2000 K, a reasonable amount of TiO is found near the stellar surface, owing to the extremely high densities encountered there.) The right-hand panel of Fig. 2 represents the star near minimum light ($`\varphi 0.5`$); here the temperature of the star has declined and TiO is abundant in the extended photosphere. At this point, one can only “see” down to a radius of $`1.8R_{}`$. When this happens the visual intensity drops dramatically, which can explain why Mira variables change visual light by factors exceeding 1000.

## 3 Model Light Curves

It is exceedingly difficult, and probably beyond current capabilities, to model precisely all of the physical and chemical processes in the atmosphere of a Mira variable. However, in order to validate the basic physical explanation for Mira light curves, we need only approximate or representative light curves. This can be accomplished provided the temperature and density as a function radius and optical phase is available. We found four papers in the recent literature that contain time-dependent models of oxygen-rich Mira variables and for which temperature and density profiles are plotted for maximum and minimum visual light (Bertschinger & Chavalier 1985, Fig. 13; Bowen 1988, Fig. 6; Bessel et al. 1989, Figs. 2 & 4; Höfner et al. 1998, Fig. 9). These profiles were digitized and used as input physical models for calculations of apparent visual magnitudes. We note that there are very large differences in the temperature and density profiles among these models.

We used the following procedure to calculate apparent visual magnitudes. At each radius in the extended photosphere, we evaluated the chemical reaction for forming molecular TiO, assuming solar abundance ratios (\[Ti\]/\[H\] = $`9.8\times 10^8`$ and \[O\]/\[H\] = $`8.5\times 10^4`$), a TiO binding energy of $`6.9`$ eV, and chemical equilibrium. Since most of the oxygen is probably bound in CO and other molecules, we assumed that only about 1% of oxygen is available to bond to titanium. The absorption cross-section of TiO as a function of wavelength was estimated, based on calculations provided by R. Kurucz, for a range of temperatures and pressures appropriate for the conditions in the outer atmosphere of a Mira variable. The dependence of the cross-sections on pressure is insignificant, but the dependence on temperature is noticeable, and we interpolated the cross-sections in temperature between values calculated at 2200 and 1400 K.

We calculated numerically the transfer of radiation through the atmosphere using a method similar to that described by Reid & Menten (1997). At a given wavelength, we integrated the intensity contribution at each point along our line-of-sight toward the star, diminished by absorption along the path to the observer. This integration was repeated for different stellar impact parameters, and a monochromatic flux density was determined by summing the integrated intensity times the solid angle of an annulus appropriate for each impact parameter. Next we repeated the integrations over a range of wavelengths to sample the entire optical spectrum. Finally, the total flux density was obtained by summing over all wavelengths, accounting for the transmission characteristics of a visual filter (Johnson, 1952). In this manner, we generated an apparent visual magnitude, $`m_v`$, for a given stellar phase.

The calculated apparent visual magnitudes at maximum light, using the published temperature and density profiles from the four papers cited above, spanned a range of about 10 magnitudes. The model of Bowen (1988) returned $`m_v=6`$, which closely matched the $`\chi `$ Cygni data at maximum light (assuming a distance of 136 pc; see discussion in next section). Other models typically yielded visual magnitudes a few magnitudes brighter. At minimum light, the model of Bowen resulted in $`m_v=9`$, about 4.5 magnitudes too bright. The differences in visual light between maximum and minimum for most of the models were only about 3 magnitudes, and for two of the four models the light at minimum exceeded that at maximum. Thus, these published temperature and density profiles cannot explain the visual light curve of a large amplitude Mira variable like $`\chi `$ Cygni. This should not, perhaps, come as a surprise. The visual light is only a small fraction of the total light from the star, and thus it is relatively unimportant in determining the physical and dynamical conditions in the atmospheric models. Since the published models have not been optimized to fit visual light, it is not surprising that current models fail in this regard. However, the visual light is a sensitive probe of physical conditions in the atmosphere of a Mira variable, and we hope that future models would be optimized to fit visual light curves.

Since published models seem inadequate for our purposes, we constructed a simple hydrostatic model for the stellar atmosphere and evaluated the transfer of radiation numerically as described above. This model is computationally manageable and its parameters are easily adjusted to attempt to fit light curves. The model assumes a star whose surface temperature and radius vary out of phase as follows:

$`T_{}=T_0+\mathrm{\Delta }T\mathrm{cos}(\varphi ),`$ (1a)
$`R_{}=R_0+\mathrm{\Delta }R\mathrm{cos}(\varphi +\pi ),`$ (1b)
where $`\varphi `$ is the visual stellar phase. The phase shift of $`\pi `$ in Eq. (1b), forces the radius to be a minimum when the temperature is maximum (Pettit & Nicholson, 1933). This star is “embedded” in an extended photosphere. Following Reid & Menten (1997), we assume that temperature as a function of radius, $`T(r)`$, in the extended photosphere decreases as follows:

$$T(r)^4=T_{}^4\left(1\sqrt{1(R_{}/r)^2}\right).$$
(2)
This temperature profile decreases approximately as $`r^{1/2}`$ (for $`r1.2R_{}`$) and, for hydrostatic equilibrium, results in a density profile, $`n(r)`$, decreasing with radius given by

$$n(r)n_{}e^{C(1(R_{}/r)^{1/2})}.$$
(3)
In Eq. (3), $`C=2GM_{eff}m_H/kT_{}R_{}`$, $`M_{eff}`$ is the effective mass of the star, $`m_H`$ is the hydrogen mass, $`k`$ is Boltzmann’s constant, and $`G`$ is the gravitational constant. Examples of density and temperature profiles (for parameters given in §4) are shown in Fig. 3.

## 4 Observed and Modeled Light Curves for $`\chi `$ Cygni

We chose to model the star $`\chi `$ Cygni because 1) it displays a very large, consistent, change in visual light, and 2) excellent visual (Mattei, 2001) and infrared (Lockwood & Wing, 1971) light curves are available. We adopted a distance to $`\chi `$ Cygni of 136 pc, based on the period–luminosity relation of Feast et al. (1989) and the photometry of Haniff, Scholz & Tuthill (1995), and an effective stellar mass of $`1M_{}`$. We modeled the star as described in §3, adjusting the parameters of the model to fit the optical and infrared light curves. A good fit to the light curves of $`\chi `$ Cygni was obtained for a mean temperature and radius of $`T_0=2200`$ K and $`R_0=2.1`$ AU, respectively, for the central star. Changes in temperature, $`\mathrm{\Delta }T`$, and radius, $`\mathrm{\Delta }R`$, over the stellar cycle were set to $`\pm 20\%,`$ and a visual extinction of $`0.6`$ mag was used. The total mass in the extended atmosphere was a held constant over the stellar cycle by adjusting the density parameter $`n_{}`$.

We present the results of the radiative transfer calculations in Fig. 4. In addition to the visual light curve, we also required the model to reproduce the infrared light curve at a wavelength of 1.04 $`\mu `$m. At this wavelength, molecular opacity is very low (Lockwood & Wing, 1971), and the infrared radiation comes from the stellar surface. As one can see, this model successfully reproduces both the small variation in the infrared (where TiO opacity is low) and the extremely large changes in visual light. This demonstrates that, for reasonable physical parameters, one can explain the “marvelous” behavior of a Mira variable as outlined in §2.

Reid & Menten (1997) measured the radio continuum emission from a sample of Mira variables and found that the radio emission comes from $`2R_{}`$. This is the same region of the extended photosphere that is optically thick at visual wavelengths, owing to molecular opacity, near minimum light. Since the radio opacity is from an entirely different source (H- free–free interactions) than the optical opacity, one can check how well our simple model predicts photospheric emission at radio wavelengths. Thus, we followed the procedures described by Reid & Menten and calculated the emission at a wavelength of 3.6 cm that would result from the model for the optical and infrared emission described above. The predicted radio emission from the model of $`\chi `$ Cygni is also shown in Fig. 4. While a complete radio “light curve” is not yet available for $`\chi `$ Cygni, there are two measures of its 3.6 cm wavelength emission made with the NRAO222NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Very Large Array: $`0.27\pm 0.04`$ mJy at $`\varphi =0.21`$ (Reid & Menten) and $`0.19\pm 0.04`$ mJy at $`\varphi =0.43`$ (Reid, Menten & Goldston, 2001). We converted these flux densities to magnitude units and plot them in Fig. 4. The model, adjusted only to optical and infrared data, predicts the strength of the radio emission reasonably well. Also, the model predicts only a small variation of the radio flux density over the stellar cycle, which is consistent with that observed for other Mira variables ($`o`$ Ceti, R Leo, and W Hya) where better phase coverage exists.

We have shown that a simple model can successfully fit the optical, infrared, and radio light curves of $`\chi `$ Cygni. However, the temperature and density profiles of Mira variables are likely far more complex. So, why is our simple model so successful? One might be concerned that the model fits the data well only because it has many adjustable parameters. Indeed, the model has 9 parameters. However, all of the parameters are physical ones, required by the problem. The stellar distance comes from Feast et al. (1989) and the stellar mass was arbitrarily set at 1 M; neither were adjusted. The two stellar temperature and radius parameters are constrained mostly by the infrared light curves. The visual extinction only shifts the entire light curve up or down; the value adopted is small and not important. This leaves only two parameters that can be used to fit the visual light curve: the oxygen depletion and a single atmospheric density parameter ($`n_{}`$). So, the good fit of the model visual light curve to the observed data cannot be a result of over-parameterization.

How sensitive is the model visual light curve to the precise nature of the temperature and density profiles? In order to investigate this question, we changed the density profile from an exponential form, appropriate for a hydrostatic atmosphere, to a power-law form, appropriate for a smoothly expanding atmosphere. Using temperature and density profiles proportional to $`r^{1/2}`$ and $`r^5`$, respectively (see Fig. 3), we were able to reproduce the large amplitude of the visual light curve of $`\chi `$ Cygni nearly as well as the exponential model, while also obtaining satisfactory agreement with the infrared and radio data. This suggests that the model light curves are somewhat robust to changes in the forms of the temperature and density profiles. Details of the variation of temperature and density with radius do not strongly affect the amplitude of the visual light curve, provided that sufficient opacity occurs near minimum light at radii where the temperature falls to a low value as discussed in §2.

## 5 Comments and Caveats

In some previous attempts to explain the phenomenon of Mira variables, the effects of TiO absorption have been considered as “extinction”—that is simply absorbing the stellar light. However, this explanation is incomplete, because the medium containing the absorbing material also re-emits light. Consider, for example, the Sun. The solar photosphere has a strong opacity source from the reaction $`\mathrm{H}^{}\mathrm{H}+\mathrm{e}^{}`$ (Wildt, 1939); indeed, essentially all visible Sun light comes from the formation of H-. However, since density falls steeply with radius in the Sun, the region of high opacity is confined to a very small range of radii. Across this region, the temperature changes only slightly, and the visible light from the Sun departs only marginally from a blackbody at a single temperature. Thus, strong absorption of light in a stellar photosphere is not alone sufficient to change the visual intensity markedly. The absorption must happen over a large range in radius, where the temperature changes significantly, if it is to dramatically change the visual intensity. This occurs in Mira variables, owing to molecular (e.g, TiO) absorption, but it does not occur in the Sun, even though there is equally strong H- absorption.

Our simple model cannot explain all details of the wide variety of light curves observed. Some Miras exhibit asymmetric light curves, which decay slowly and rise rapidly toward maximum light; others display secondary “humps” in their curves. Also, Mira variables are not perfectly spherical—some are elongated by 10% or more (Karovska et al., 1991). Photospheric shocks and mass outflow generally tend to soften the density profile (i.e., increase the density at larger radii) above that obtained in pure hydrostatic equilibrium. On the other hand, there is clear evidence that Mira atmospheres are not smoothly expanding. Regions of large-scale infall are indicated by red-shifted absorption lines (Reid, 1976; Reid & Dickinson, 1976; Ferlet & Gillet, 1984) and they have been directly observed in the proper motions of some SiO masers (Diamond & Kemball, 1999). Thus, the combination of outflowing and infalling regions, perhaps better described as convective motions, when averaged over the entire atmosphere, might result in an average density profile that deviates less from a hydrostatic profile than expected from pure outflow. This needs to be studied. Still, explaining many details of Mira emissions requires more sophisticated modeling. The simple hydrostatic model used here was chosen for computational ease, since this paper is concerned with explaining the basic physical reasons for the large amplitude of the visual light curve and not with constructing a detailed atmospheric model. Ultimately, full dynamical models should be applied to this problem; these models should be constrained by observed light curves, in addition to observed spectra and measured sizes.

Our simplifying assumption, that the visual opacity is dominated by TiO, is only a first approximation. Other sources contribute to the overall opacity. For example, stars of spectral class S exhibit strong bands from ZrO. Indeed, $`\chi `$ Cygni has been variously classified as an M-type (strong TiO bands) and an S-type (strong ZrO bands in addition to TiO) by different observers at different epochs. Including sources of opacity beyond TiO would decrease somewhat the required total density to achieve similar total opacity. This would further enhance the effects leading to the large amplitude of the visual light curves and make our explanation for these light curves more robust. Some carbon-rich Mira variables of spectral class N and R are also known to have large amplitude light curves (Merrill, 1940). For these stars, essentially all oxygen is in CO and metallic oxides, such as TiO, cannot form abundantly. Thus, the explanation for the behavior of the light curves of carbon-rich Mira variables probably involves strong opacity from carbon-bearing molecules such as C2 and CN (Reid & Goldston, 2002). Finally, dust can form at temperatures below $`15001200`$ K in carbon or oxygen rich Miras. Indeed, dust shells are directly observed in some Miras, usually at radii greater than about $`3R_{}`$ (Danchi et al., 1994), and these can also contribute to the extinction of optical light (e.g., Winters et al. 1994).

In conclusion, we mention that strong evidence for our basic picture comes from interferometric observations at visible and radio wavelengths. Interferometric observations in the visible show large changes in the apparent size of the star for wavelengths in and out of TiO absorption bands (Labeyrie et al., 1977; Tuthill, Haniff & Baldwin, 1999; Young et al., 2000). For example, imaging observations made with the Palomar 5-m telescope clearly show a “shell” of high opacity at $`2R_{}`$ in a strong TiO absorption band (Haniff et al., 1992). Independent observations at radio wavelengths (Reid & Menten, 1997) directly measure the (electron) temperature at $`2R_{}`$ to be $`1500`$ K. This is within the range needed to cause the large changes in visual light as outlined in this paper. Thus, based on direct observational evidence, the simplicity of the mechanism described above, and the ability of our models to reproduce robustly the essential characteristics of Mira light curves, we feel that the physical picture outlined here may resolve the long-standing problem of the large amplitude of the light curves of Mira variables.

We thank R. Kurucz for providing the TiO cross-section data and the AAVSO for the visual data on $`\chi `$ Cygni from their International Database, which comprises observations by variable star observers worldwide.


text:
# 1 : 1a, 1b, 1c show the specklegram of HR4689, the autocorrelation of HR4689 and 𝜎_{𝑛𝑜𝑖𝑠𝑒} vs. WFP plot respectively. The axes of the figures 1(a) and 1(b) are the pixel value; each pixel value is 0.015 arc-seconds. Since the inherent property of the autocorrelation method is to produce double images with 180∘ ambiguities of a binary source, one of the two contours on either side of the central one (fig 1b) is the secondary component; the central contours represents the primary component. The contours at the corners are the artifacts.

Auto-correlation of Binary Stars

S. K. Saha111e-mail: sks@iiap.ernet.in and D. Maitra 222Summer project student from Indian Institute of Technology Kanpur 208016, India

Indian Institute of Astrophysics

Bangalore 560034

India

Abstract: Speckle interferometric technique is used to record a series of short exposure images of several close binary stars with sub-arcsecond separation through a narrow band filter centred at H$`\alpha `$ at the Cassegrain focus of the 2.34 meter (m) Vainu Bappu telescope (VBT), situated at Vainu Bappu Observatory (VBO), Kavalur . The auto-correlation method is developed under Image Reduction Analysis Facility (IRAF). Wiener filter is included in the programme to eliminate spurious high frequency contributions; a few sets of data provide the optimised results. The auto-correlated image of these stars gives the separation of the binary components.

Keywords: Speckle Imaging, Auto-correlation, Binary stars, Wiener filter.

PACS Nos: 07.60.Ly, 42.30.Wb, 95.55Br

1. Introduction

The speckle interferometric technique  is widely used to decode the deleterious effect of the atmospheric turbulence that limits the sensitivity of ground-based telescopes. Recent reviews by Saha  describe in depth the utility of this method. One of the important subjects in stellar astrophysics is the studies of close binary stars, that play a fundamental role in measuring stellar masses, providing a benchmark for stellar evolution calculations; a long-term benefit of interferometric imaging is a better calibration of the main-sequence mass-luminosity relationship . The procedure in obtaining masses of stars involves combining the spectroscopic orbit with the astrometric orbit, which is the projection of true orbit on the sky, from the interferometric data . In order to derive accurate orbital elements and masses, luminosities and distances, continuous observations of spectroscopic binaries are essential. Major contributions in resolving close binary systems, among others, came from the Center for High Angular Resolution Astronomy (CHARA) at Georgia State University, USA \[6-9\]. Survey of visual and interferometric binary stars have also been reported by the other groups \[10-19\].

Most of these results that appeared in various journals till date are from the usage of the auto-correlation technique though schemes, like Knox-Thomson algorithm , shift and add , triple correlation , and Blind Iterative deconvolution (BID) technique  restore the degraded images. All these schemes, barring BID, depend on the statistical treatment of a large number of images. Often, it may not be possible to record a large number of images within the time interval over which the statistics of the atmospheric turbulence remains stationary. In such cases, where only a few images are available, we describe an auto-correlation algorithm that is developed under IRAF; this algorithm has been applied on degraded images of several binaries, obtained at the 2.34 m VBT at Kavalur.

2. Outline of the theory of speckle interferometry

The instantaneous image intensity distribution $`I(𝐱)`$ recorded by the detector (CCD/photographic plate/eye etc.) is actually a convolution of the instantaneous ‘Point Spread Function (PSF)’ $`S(𝐱)`$, with the actual intensity distribution of the object $`O(𝐱)`$. Here $`𝐱`$ is the two dimensional position vector. The PSF describes how light from a point source (a two dimensional $`\delta `$ function on the sky) gets distributed across the detector. The variability of the corrugated wave-front yields ’speckle boiling’ and is the source of speckle noise that arises from difference in registration between evolving speckle pattern and the boundary of the PSF area in the focal plane. These specklegrams have additive noise contamination, $`N(𝐱)`$, which includes all additive measurement of uncertainties. This may be in the form of (i) photon statistics noise, and (ii) all distortions from the idealized isoplanatic model represented by the convolution of $`O(𝐱)`$ with $`S(𝐱)`$ that includes non-linear geometrical distortions.

If instead of dealing in the real space we choose to work in Fourier space, our life becomes simpler. This is because the operation of convolution in real space becomes simple multiplication in the Fourier space. Mathematically speaking:

$$I(𝐱)=O(𝐱)S(𝐱)+N(𝐱),$$
(1)
where ‘$``$’ signifies convolution operation.

In Fourier domain, equation (1) can be written as:

$$I(𝐮)=O(𝐮)S(𝐮)+N(𝐮),$$
(2)
where $`𝐮`$ is the two dimensional spatial frequency vector and ‘$``$’ implies the ordinary operation of multiplication. The quantity $`O(𝐮)`$ is called the Object Spectrum, and $`S(𝐮)`$ the Transfer function.

If the exposure time is short enough then the atmosphere can be taken to be static during the exposure. The common term that describes this is that ‘the atmosphere has been freezed during the exposure’. From experiments it has been found that for exposures below 20 ms the atmosphere can be considered as fairly static for reasonably stable nights.

A large number of such short exposure images (specklegrams) of the object of interest is taken. Also images of a reference star which is known to be a unresolved star (i.e., a point source, a two dimensional $`\delta `$ function on the sky) and which lies in the same isoplanatic patch as that of the program star, is taken at the same time. It is assumed that the atmospheric distortion affects the reference star in the same way as it does to the program star.

The ensemble average of the power spectrum is given by

$$=O(𝐮)^2<S(𝐮)^2>+<N(𝐮)^2>.$$
(3)
The reference star being a single star, the Fourier transform of its image is unity, i.e., $`O_{ref}(𝐮)`$ = 1.

The power spectrum of the program star is therefore given by

$$O_{pro}(𝐮)^2=\frac{<I_{pro}(𝐮)^2>}{<I_{ref}(𝐮)^2>},$$
(4)
and by Wiener-Kinchin theorem, the autocorrelation of the object is given by

$$A[O(𝐱)]=T^1[O(𝐮)^2],$$
(5)
where $`A`$ denotes Autocorrelation and $`T^1`$ denotes the operation of inverse Fourier transform.

3. Data Processing and Noise Reduction

The program for finding the power spectrum is written in FORTRAN. It uses IRAF commands and hence the macros in IRAF are called and used. The program is written in such a way that a large number of data files can be handled and large number of images can be averaged to increase the signal-to-noise ratio. The set of input specklegrams of object and the reference stars are read and their power spectrum is calculated. The object power spectrum and the image power spectrum are averaged in the Fourier space. Finally the autocorrelated image is calculated using another programme.

The disadvantage with a division as in equation (4), is that the zeros in the denominator will corrupt the ratio and spurious high frequency components will be created in the reconstructed image. Moreover, a certain amount of noise is inherent in any kind of observation. In our present case noise is primarily due to thermal electrons in the CCD interfering with the signal. Most of this noise is in the high spatial frequency regime.

In order to get rid of this high frequency noise as much as possible, a Wiener filter was used. Applying the Wiener filter is essentially the process of convolving the noise degraded image $`I(𝐱)`$ with the Wiener filter. The Wiener filter is applied in the frequency domain. Then the original image spectrum $`I(𝐮)`$ is estimated from the degraded image spectrum $`I^{}(𝐮)`$ by simply multiplying the image spectrum with the Wiener filter $`W(𝐮)`$. However, this will reduce the resolution in the reconstructed image. The advantage is the spurious higher frequency contribution are eliminated.

$$I(𝐮)=I^{}(𝐮)W(𝐮).$$
(6)
The Wiener filter, in the frequency domain, has the following functional form:

$$W(𝐮)=\frac{H^{}(𝐮)}{H(𝐮)^2+\frac{P_n(𝐮)}{P_s(𝐮)}}.$$
(7)
where, $`H(𝐮)`$ is the Fourier transform of the point-spread-function, $`P_s(𝐮)`$ the power spectrum of signal process, and $`P_n(𝐮)`$ the power spectrum of noise process.

The term $`P_n(𝐮)`$/$`P_s(𝐮)`$ can be interpreted as reciprocal of signal to noise ratio. In our case, the noise is due to the CCD. We developed an IRAF-based algorithm, where a Wiener parameter, $`w`$, is added to the PSF power spectrum in order to avoid zeros in the PSF power spectrum that helps in reconstructions with a few frames. The classic Weiner filter that came out of the electronic information theory where diffraction-limits do not mean much, is meant to deal with signal dependent ’coloured’ noise. In practice, this term is usually just a constant, a ’noise control parameter’ whose scale is estimated from the noise power spectrum. In this case, it assumes that the noise is white and that one can estimate its scale in regions of the power spectrum where the signal is zero (outside the diffraction-limit for an imaging system). The expression for the Wiener filter simplifies to:

$$I(𝐮)=\frac{P_{ref}(𝐮)}{P_{ref}(𝐮)+w},$$
(8)
where, $`w`$ is the noise-variance and termed as Wiener filter parameter in the program.

To get an optimally autocorrelated image, a judicious choice of the Wiener filter parameter is made according to the procedure described below:

For a very wide range of Wiener filter parameter values, the autocorrelated images are constructed. A small portion (16 X 16 pixels) of each image, far from the centre, is sampled to find-out the standard deviation in the intensity values of the pixels, $`\sigma _{noise}`$. The plot of standard deviation thus obtained against the Wiener filter parameter (Figures 1(c)) shows a minimum. The abscissa corresponding to this minimum gives the optimum Wiener filter parameter value.

The nature of $`\sigma _{noise}`$ vs. WFP plot is understood as follows:
The noise in our data is primarily in the high frequency region whereas the signal is at a comparatively low frequency. With zero WFP, there is no filtering and hence there is ample noise. As the WFP value is gradually increased from zero, more and more high frequency noise is cut off and the $`\sigma _{noise}`$ goes down. It attains a minimum when WFP value is just enough to retain the signal and discard the higher frequency noise. However, at higher WFP region, we are over compensating for the signal, and therefore blurring of the image starts to occur. This leads to the sharp increase in $`\sigma _{noise}`$. The ‘ringing’ effect due to sharp cutoff comes into play.

Computer simulations were carried out by convolving ideal star images having Gaussian profile with a random PSF to generate speckle pattern. The plot of $`\sigma _{noise}`$ vs. WFP also gives similar results though some of the major sources of noise, e.g., thermal noise due to electron motion in the CCD, effect due to cosmic rays on the frame, etc. were not incorporated in the simulation.

4. Instrumentation

The 2.34 m VBT has two accessible foci for back end instrumentation $``$ a prime focus (f/3.25 beam) and a Cassegrain focus (f/13 beam). The latter was used for the observations described in this paper. The Cassegrain focus has an image scale of 6.7 arcseconds per mm. This was further magnified using a microscope objective . This enlarged image was recorded through a 5 nm filter centred at H$`\alpha `$ using an EEV intensified CCD camera which provides a standard CCIR video output of the recorded scene. The interface between the intensifier screen and the CCD chip is a fibre-optic bundle which reduces the image size by a factor of 1.69. A DT-2851 frame grabber card digitises the video signal. This digitiser re samples the pixels of each row (385 CCD columns to 512 digitized samples) and introduces a net reduction in the row direction by a factor of 1.27. Finally an image-scale of 0.015 arcsecond per pixel is achieved. The frame grabber can store upto two images on its onboard memory. These images are then written onto the hard disc of a computer.

5. Results

The separations of HR4689 and HR5858, observed on 28 February 1997 were estimated using the procedure described above. The average seeing was around 1.5 arcsec. Figures 1a, 1b, 1c show the specklegrams of HR4689, minimum noise autocorrelated image of HR4689, and variation of noise with Wiener filter parameter respectively. The separation of HR4689 was found to be 0.17 $`\pm `$ 0.03 arcseconds and that of HR5858 was 0.23 $`\pm `$ 0.03 arcseconds. The separations observed are compatible with the values published in the CHARA catalogue .

6. Discussion and conclusions

The program of speckle imaging of binary stars at VBT, Kavalur has been a successful one. Several experiments and computer simulations have been carried out for improving the speckle interferometer at VBT . Algorithms are in stage of development for Fourier phase recovery too. Noise in the data is reduced to an optimum level with the judicious use of Wiener filter. Recent developments in CCD technology and associated photon counting devices of high quantum efficiency promise better signal to noise ratio, and hence better accuracy in determining the orbital parameters of binaries. However, in this algorithm, it gives a better signal to noise ratio if the reference star is brighter than the program star, both being in the same isoplanatic patch. This is very difficult to achieve in reality and when it is not fulfilled, we have tried to make a compensation by increasing the number of frames for the reference.

Acknowledgments: We would like to thank Mr. K. Sankarasubramanian for his constant help throughout the work, primarily in developing the algorithms and implementing them with success. The services rendered by Messrs V. Chinnappan, P. Anbazhagan, and V. Murthy are also greatfully acknowledged.

References

 A. Labeyrie, 1970, Astron. Astrophys., 6, 85.
 S. K. Saha, 1999a, Bull. Astron. Soc. Ind., 27, 443.
 S. K. Saha, 1999b, Ind. J. Phys., 73B, 553.
 H. A., McAlister, 1988, Proc. ESO-NOAO conf. ‘High Resolution Imaging Interferometry’, ed. F. Merkle, Garching bei München, FRG, 3.
 G. Torres, R. P. Stefanik, and D. W. Latham, 1997, Astrophys. J, 485, 167.
 W. I. Hartkopf, B. D. Mason, H. A. McAlister, N. H. Turner, D. J. Barry, G. Franz, and C. M. Prieto, 1996, Astron. J, 111, 936.
 W. I. Hartkopf, H. A. McAlister, and B. D. Mason, 1997, CHARA Contrib. No. 4, ‘Third Catalog of Interferometric Measurements of Binary Stars’, W.I.
 H. A. McAlister, W. I. Hartkopf, B. D. Mason, L. C. Roberts (Jr.), and M. M. Shara, 1996, Astron. J, 112, 1169.
 B. D. Mason, D. R. Gies, W. I. Hartkopf, W. G. Bagnuolo (Jr.), T. ten Brummelaar, and H. A. McAlister, 1998, Astron. J, 115, 821.
 D. Bonneau, Y. Balega, A. Blazit, R. Foy, F. Vakili, and J. L. Vidal, 1986, Astron. Astrophys. Suppl., 65, 27.
 A. Blazit, D. Bonneau, and R. Foy, 1987, Astron. Astrophys. Suppl., 71, 57.
 N. Miura, N. Baba, S. Isobe, M. Noguchi, and Y. Norimoto, 1992, J. Modern Opt., 39, 1137.
 I. I. Balega, Y. Y. Balega, I. N. Belkin, A. E. Maximov, V. G. Orlov, E. A. Pluzhnik, Z. U. Shkhagosheva, and V. A. Vasyuk, 1994, Astron. Astrophys. Suppl., 105, 503.
 E. P. Horch, Z. Ninkov, W. F. van Altena, R. D. Meyer, T. M. Girard, and J. G. Timothy, 1999, Astron. J, 117, 548.
 G. G. Douglass, R. B. Hindsley, and C. E. Worley, 1997, Astrophys. J Suppl., 111, 289.
 C. Leinert, A. Richichi, and M. Haas, 1997, Astron. Astrophys., 318, 472.
 S. K. Saha, and P. Venkatakrishnan, 1997, Bull. Astron. Soc. Ind., 25, 329.
 M. E. Germain, G. G. Douglass, and C. E. Worley, 1999, Astron. J, 117, 1905.
 S. K. Saha, R. Sridharan, and K. Sankarasubramanian, 1999a, ‘Speckle image reconstruction of binary stars’, Presented at the Astron. Soc. Ind. meeting, held at Bangalore, 1999.
 Knox K.T. & Thompson B.J., 1974, Ap.J. Lett., 193, L45.
 Lynds C.R., Worden S.P., Harvey J.W., 1976, Ap.J., 207, 174.
 Lohmann A.W., Weigelt G., Wirnitzer B., 1983, Appl. Opt., 22, 4028.
 Ayers G.R. & Dainty J.C., 1988, Optics Letters, 13, 547.
 S. K. Saha, G. Sudheendra, A. Umesh Chandra, and V. Chinnappan, 1999b, Experimental Astron., 9, 39.
 McAlister H. A. & Hartkopf W. I., 1988, CHARA contribution no. 2.


text:
# 1 Five parameters of the HLS at 𝜇=Λ and 𝑚_𝜌 for Λ_QCD=400 MeV. The unit of 𝐹_𝜋 is GeV.

## Acknowledgment

I would like to thank Professor Koichi Yamawaki for collaboration in Refs. References, References and References on which this talk is based. I would be grateful to the organizers for giving me an opportunity to present this talk. This work is supported in part by Grant-in-Aid for Scientific Research (A)#12740144.


text:
# Optimal orientation of striped states in the quantum Hall system against external modulations

## Abstract

We study striped states in the quantum Hall system around half-filled high Landau levels and obtain the optimal orientation of the striped state in the presence of an external unidirectional periodic potential. It is shown that the optimal orientation is orthogonal to the external modulation in the Coulomb dominant regime (the orthogonal phase) and is parallel in the external modulation dominant regime (the parallel phase). The phase boundary of these two phases is determined numerically in the parameter space of the strength and wave number of the external modulation at the half-filled third Landau level.

The modern semiconductor technology yields extremely pure two-dimensional (2D) electron systems in heterostructures. In the presence of a strong perpendicular magnetic field, various kinds of new phenomena in addition to the ordinary quantum Hall effect are found in the systems. The effect of a crystal structure has been believed to be ignored, because the magnetic length is much larger than the lattice constant of the host crystal. Hence, the system is supposed to have an orientational symmetry, that is, physics in the $`x`$ direction and $`y`$ direction is equivalent. In the end of the last century, however, highly anisotropic states, which have an enormous anisotropy in magnetoresistances, were observed around the half-filled high Landau level (LL). This observation agrees with the striped state which was predicted in a mean field theory. The charge density of the striped state is uniform in the one direction and periodic in the orthogonal direction. The Hartree-Fock (HF) theory and numerical calculations in a small system for striped states were studied recently at the filling factor where the anisotropy is observed and the results seem to support the striped state. The experiments show that the stripe direction is parallel to the specific crystallographic direction. The origin of the orientation of the striped state is still puzzling.

It is naturally considered that the origin of the orientation is related to yet unknown and very weak periodic structure in the sample. Therefore, we suppose that the origin to determine the orientation can be modeled by an external modulation in the 2D electron system. We compute the energy of striped states in the quantum Hall system with an external modulation and find the optimal orientation. It is naively expected that the external modulation breaks the orientational symmetry and makes the orientation of the striped state parallel to the modulation. Unexpectedly, we find the counter-intuitive phase in which the optimal orientation of the striped state becomes orthogonal to the external modulation, and the other phase in which the optimal orientation is parallel to the external modulation depending on the strength and wave number of the modulation. Our results are consistent with the recent experiments.

Let us consider a 2D electron system in a perpendicular uniform magnetic field $`B`$ and a unidirectional periodic potential. The total Hamiltonian $`H`$ of the system is written as $`H=H_0+H_1+H_2`$,

$`H_0`$ $`=`$ $`{\displaystyle \psi ^{}(𝐫)\frac{(\widehat{𝐩}+e𝐀)^2}{2m}\psi (𝐫)d^2r},`$ (1)
$`H_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \rho (𝐫)V(𝐫𝐫^{})\rho (𝐫^{})d^2rd^2r^{}},`$ (2)
$`H_2`$ $`=`$ $`g{\displaystyle \rho (𝐫)\mathrm{cos}(𝐊𝐫)d^2r},`$ (3)
where $`\widehat{p}_\alpha =i\mathrm{}_\alpha `$, $`_xA_y_yA_x=B`$, $`V=q^2/r`$, $`q^2=e^2/4\pi ϵ`$ ($`ϵ`$ is a dielectric constant). $`\psi (𝐫)`$ is the electron field, and $`\rho (𝐫)=\psi ^{}(𝐫)\psi (𝐫)`$. $`H_0`$ is the free Hamiltonian, which is quenched in the LL. $`H_1`$ is the Coulomb interaction term and $`H_2`$ is the external modulation term. We ignore the spin degree of freedom.

The electron field is expanded by the momentum state $`|f_l\beta _𝐩`$ in von Neumann lattice (vNL) formalism as

$$\psi (𝐫)=\underset{l=0}{\overset{\mathrm{}}{}}_{\mathrm{BZ}}\frac{d^2p}{(2\pi )^2}b_l(𝐩)𝐫|f_l\beta _𝐩,$$
(4)
where $`b_l(𝐩)`$ is the anticommuting annihilation operator with a LL index $`l`$ and a 2D lattice momentum $`𝐩`$ defined in the Brillouin zone $`|p_\alpha |\pi `$. The annihilation operator $`b_l(𝐩)`$ obeys a twisted periodic boundary condition $`b_l(𝐩2\pi 𝐍)=e^{i\pi (N_x+N_y)+iN_yp_x}b_l(𝐩)`$, where $`N_x`$, $`N_y`$ are integers. The momentum state is Fourier transform of the Wannier basis of vNL which localized at $`𝐫=a(r_sm,n/r_s)`$, where $`n`$, $`m`$ are integers. Here $`a=\sqrt{2\pi \mathrm{}/eB}`$ and $`r_s`$ is an asymmetric parameter. We consider only the $`l`$ th LL state and ignore the LL mixing. Hence, $`H_0`$ turns out to be constant and we omit the free Hamiltonian. Fourier transformed density operator $`\stackrel{~}{\rho }`$ in the $`l`$ th LL is written in vNL formalism as

$`\stackrel{~}{\rho }(𝐤)={\displaystyle _{\mathrm{BZ}}}{\displaystyle \frac{d^2p}{(2\pi )^2}}b_l^{}(𝐩)b_l(𝐩a\widehat{𝐤})f_l(k)`$ (5)
$`\times \mathrm{exp}[i{\displaystyle \frac{a}{4\pi }}\widehat{k}_x(2p_ya\widehat{k}_y)],`$ (6)
where $`\widehat{𝐤}=(r_sk_x,k_y/r_s)`$ and $`f_l(k)=L_l(\frac{a^2k^2}{4\pi })e^{\frac{a^2k^2}{8\pi }}`$, here $`L_l`$ is the Laguerre polynomial. We substitute Eq. (5) into $`H_1`$ and $`H_2`$ and find the ground state in two perturbative approaches. In the first approach, perturbative expansions with respect to $`g`$ in $`H_2`$, which describe the Coulomb dominant regime, are applied. In the second approach, perturbative expansions with respect to $`q^2/a`$ in $`H_1`$, which describe the external modulation dominant regime, are applied. The filling factor is fixed at $`l+1/2`$ in the following calculation and numerical calculations are performed at $`l=2`$.

(I) $`H_2`$ as a perturbation : We obtain the ground state of $`H_1`$ in the HF approximation first. Using the HF ground state, we treat $`H_2`$ as a perturbation. This approximation is relevant in the Coulomb dominant regime, $`gq^2/a`$. We use the striped state $`|\mathrm{\Psi }_1`$ which is uniform in the $`y`$ direction as the unperturbed ground state of the $`H_1`$. In the HF approximation, this is given as

$$|\mathrm{\Psi }_1=N_1\underset{|p_x|\pi ,|p_y|\pi /2}{}b_l^{}(𝐩)|0,$$
(7)
where $`|0`$ is the vacuum state for $`b_l`$ and $`N_1`$ is a normalization factor. The Fermi surface is parallel to the $`p_x`$ axis. The density of this state $`\mathrm{\Psi }_1|\rho (𝐫)|\mathrm{\Psi }_1`$ is uniform in $`y`$ direction and periodic in $`x`$ direction with a period $`ar_s`$. The orthogonality of the Fermi surface in the momentum space and the density in the coordinate space plays important roles and is reminiscent of the Hall effect. The one-particle spectrum is given by $`ϵ_{\mathrm{HF}}=ϵ_\mathrm{H}+ϵ_\mathrm{F}`$,

$`ϵ_\mathrm{H}`$ $`=`$ $`{\displaystyle \frac{2r_sq^2}{a\pi }}{\displaystyle \underset{n=odd}{}}f_l(2\pi n/ar_s)^2{\displaystyle \frac{\mathrm{cos}(np_y)}{n^2}},`$ (8)
$`ϵ_\mathrm{F}`$ $`=`$ $`{\displaystyle \frac{r_sq^2}{a}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _{\frac{\pi }{2}+p_y+2\pi n}^{\frac{\pi }{2}+p_y+2\pi n}}𝑑k_y{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dk_x}{2\pi }}h_l(𝐤),`$ (9)
where $`h_l(𝐤)=f_l(\sqrt{k_x^2+(k_yr_s/a)^2})^2/\sqrt{k_x^2+(k_yr_s/a)^2}`$. $`ϵ_{\mathrm{HF}}`$ depends on only $`p_y`$, and the self-consistency condition for $`|\mathrm{\Psi }_1`$ is satisfied. The Fermi velocity is in the $`y`$ direction. The HF energy per particle is calculated as $`E_{\mathrm{HF}}(r_s)=\mathrm{\Psi }_1|H_1|\mathrm{\Psi }_1/N`$ where $`N`$ is a number of electrons. $`E_{\mathrm{HF}}`$ is a function of $`r_s`$ and calculated as $`E_{\mathrm{HF}}=_{\pi /2}^{\pi /2}\frac{dp_y}{2\pi }ϵ_{\mathrm{HF}}(p_y)`$ at the half-filled $`l`$ th LL. The optimal value $`r_s=r_s^{\mathrm{min}}`$ is determined so as to minimize $`E_{\mathrm{HF}}(r_s)`$. The numerical value of $`r_s^{\mathrm{min}}`$ is $`2.47`$ at the half-filled third LL.

The perturbation energy per particle in the first order $`\mathrm{\Delta }E^{(1)}=\mathrm{\Psi }_1|H_2|\mathrm{\Psi }_1/N`$ is written as

$$\mathrm{\Delta }E^{(1)}=\frac{g}{2N}\mathrm{\Psi }_1|(\stackrel{~}{\rho }(𝐊)+\stackrel{~}{\rho }(𝐊))|\mathrm{\Psi }_1.$$
(10)
The operator $`\stackrel{~}{\rho }(𝐊)`$ moves an electron in the Fermi sea by $`a\widehat{𝐊}`$ in the momentum space. Therefore, except for the case that $`a\widehat{K}_y`$ coincides with a multiple of $`2\pi `$, $`\mathrm{\Delta }E^{(1)}`$ vanishes. We consider only the range $`|a\widehat{K}_y|<\pi `$ which is sufficient to compare our results with experiments.

The perturbation energy per particle in the second order $`\mathrm{\Delta }E^{(2)}`$ is written as

$$\mathrm{\Delta }E^{(2)}(g,\theta ,K)=_{\frac{\pi }{2}a\widehat{K}_y}^{\frac{\pi }{2}}\frac{dp_y}{2\pi }\frac{g^2f_l(K)^2}{ϵ_{\mathrm{HF}}(p_y+a\widehat{K}_y)ϵ_{\mathrm{HF}}(p_y)},$$
(11)
where $`\widehat{K}_y=K\mathrm{sin}\theta /r_s^{\mathrm{min}}`$, $`\theta `$ is an angle between the stripe direction and external modulation. We obtain the total energy per particle in the Coulomb dominant regime as

$$E^{\mathrm{Coul}}(g,K,\theta )=E_{\mathrm{HF}}(r_s^{\mathrm{min}})+\mathrm{\Delta }E^{(2)}(g,K,\theta ).$$
(12)
The $`\theta `$ dependence of $`\mathrm{\Delta }E^{(2)}`$ is shown in Fig. 1 at the half-filled third LL. As seen in Fig. 1, the energy is always minimum at $`\theta =\pi /2`$, that is, the optimal orientation of the striped state is orthogonal to the external modulation. We call this phase the orthogonal phase. Note that $`\mathrm{\Delta }E^{(2)}`$ vanishes and $`\theta `$ dependence disappears when $`K`$ equals the zeros of $`f_l(K)`$. In this case, the external modulation loses control of the stripe direction.

(II) $`H_1`$ as a perturbation : We diagonalize $`H_2`$ first by choosing the $`y`$ axis of vNL to be parallel to the external modulation and $`r_s=2\pi /aK`$, that is, the period of the striped states $`ar_s`$ equals the wave length $`2\pi /K`$ of the external modulation. Using this vNL basis, we treat $`H_1`$ as a perturbation. This approximation is relevant in the external modulation dominant regime, $`gq^2/a`$. Then the state $`|\mathrm{\Psi }_1`$ is the ground state of $`H_2`$. The external modulation term reads

$$H_2=|gf_l(K)|_{\mathrm{BZ}}\frac{d^2p}{(2\pi )^2}b_l^{}(𝐩)b_l(𝐩)\mathrm{cos}p_y.$$
(13)
The ground state energy per particle of $`H_2`$ is obtained as $`E_2(g)=\frac{2}{\pi }|gf_l(K)|`$. The perturbation energy per particle is calculated as $`E_{\mathrm{HF}}(\mathrm{r}_\mathrm{s})=\mathrm{\Psi }_1|\mathrm{H}_1|\mathrm{\Psi }_1/\mathrm{N}`$ with $`r_s=2\pi /aK`$. Hence, the total energy per particle in the external modulation dominant regime is given by

$$E^{\mathrm{ext}}(g,K)=E_{\mathrm{HF}}(2\pi /aK)\frac{2}{\pi }\left|gf_l(K)\right|.$$
(14)
In this state, the stripe direction is parallel to the external modulation, which we call the parallel phase.

Finally we compare the total energies of striped states obtained in (I) and (II) at the half-filled third LL. As a typical example, $`E^{\mathrm{Coul}}(g,K,\pi /2)`$ and $`E^{\mathrm{ext}}(g,K)`$ for $`aK=2`$ is plotted in Fig. 2. The bold line represents the lower energy state. As seen in Fig. 2, the orthogonal phase has lower energy and is realized at small $`g`$. The parallel phase has lower energy and is realized at large $`g`$. The phase boundary is calculated by solving the equation in $`g`$, $`E^{\mathrm{ext}}(g,K)=E^{\mathrm{Coul}}(g,K,\pi /2)`$ for various value of $`aK`$. The phase diagram in $`g`$-$`K`$ plane is shown in Fig. 3, where orthogonal and parallel phases are indicated by I and II respectively. The dashed lines correspond to the zeros of $`f_2(K)`$, $`aK=2.714`$ and $`6.550`$, at which the stripe direction is undetermined. At $`aK=2\pi /r_{\mathrm{min}}=2.544`$, where the period of the external modulation coincides with the optimal period of the stripe, the phase boundary touches the $`K`$ axis.

The direct verification of our results is made by observing a transition between the two phases by tuning the extenal modulation. The necessary wave length of the modulation for the verification is on the order of $`2a`$, which is about 100nm at $`B=2\mathrm{T}`$. Recently, the unidirectional lateral superlattice with a period 92nm is achieved on top of the 2D electron system. The experiment shows that the magnetoresistance orthogonal to the external modulation has a shallow and broad dent between two peaks around $`\nu =9/2`$. The magnetoresistance parallel to the external modulation does not have the same structure around $`\nu =9/2`$. Anisotropy observed in this experiment is small due to the low mobility compared with experiments of striped states . The strength of the modulation is estimated as $`g=0.015\mathrm{meV}`$. The parameters $`(g,K)=(0.006q^2/a,3.097/a)`$ correspond to this experimental setting, and are shown as $`X`$ in Fig. 4. $`X`$ belongs to the orthogonal phase. In this phase, the one-particle dispersion has no energy gap in the orthogonal direction and has an energy gap in the parallel direction to the external modulation. Hence, the magnetoresistance orthogonal to the external modulation is strongly modified by the injected electric current compared with the parallel magnetoresistance. This is consistent with the experiment. With a slightly larger period 115nm, orthogonal magnetoresistance around $`\nu =9/2`$ is structureless. In this case, the corresponding parameters $`(g,K)=(0.014q^2/a,2.478/a)`$ are shown as $`Y`$ in Fig. 4 . $`Y`$ belongs to the parallel phase. Hence, there is an energy gap in orthogonal direction to the external modulation and the magnetoresistance in this direction is not modified strongly by the injected electric current. This is also consistent with the experiment. We hope that a similar experiment with a higher mobility sample will give more clear evidence for our results.

In the half-filled lowest LL, the anisotropic effect is observed under the external modulation . The semiclassical composite fermion theory is proposed for this experiment. In this theory, it is assumed that anisotropic transport is caused by the density modulation. On the other hand, we study the striped state under the external modulation in the half-filled higher LL. Note that the origin of the anisotropy in the present case is spontaneous stripe formation rather than the external modulation.

It seems difficult to understand how the orthogonal phase is realized contrary to the naive expectation that two striped structures tend to be parallel. To understand the reason, it is convenient to consider the striped state in the momentum space. Since the Fermi surface of the striped state is flat as seen in Eq. (7), a perturbation with a wave number vector perpendicular to the Fermi surface affects the total energy most strongly. Therefore the orthogonal phase could be realized in a small external modulation. The point of our theory is that the mean field theory has the flat Fermi surface. The fluctuation around the mean field has been studied but discussions seem unsettled yet. Comparisons between experiments with an in-plane magnetic field and HF calculations indicate that the mean field energy is good approximation for the striped state. Higher order corrections are expected to be small because the Fermi velocity diverges due to the Coulomb interaction. We estimate $`\mathrm{\Delta }E^{(2)}`$ in the RPA approximation to the density correlation function as $`\pi _{00}^{\mathrm{RPA}}(k)=\pi _{00}(k)/(1\stackrel{~}{V}(k)\pi _{00}(k))`$. The results are shown in Fig. 1 by dashed lines. As seen in this figure the correction is small actually.

In summary it is shown that a weak external modulation determines the orientation of the striped state and there are two phases in the 2D parameter space of the strength and wave number of the external modulation, that is, the orthogonal phase and parallel phase. In the former phase, the optimal orientation of the striped state is orthogonal to the external modulation. In the latter phase, the optimal orientation is parallel to the external modulation. The phase diagram is obtained numerically at the half-filled third LL. We believe that our findings shed a new light on the origin of an orientation of striped states in quantum Hall systems.

We thank A. Endo and Y. Iye for useful discussions. This work was partially supported by the special Grant-in-Aid for Promotion of Education and Science in Hokkaido University provided by the Ministry of Education, Science, Sports, and Culture, and by the Grant-in-Aid for Scientific Research on Priority area (Physics of CP violation) (Grant No. 12014201).