text:
# Quantification of depth of anesthesia by nonlinear time series analysis of brain electrical activity

## I Introduction

Several recent publications have used analysis methods derived from the theory of nonlinear dynamics to describe aspects of neuronal activity. Measures like the estimated largest Lyapunov exponent (LLE) or the correlation dimension ($`D_2`$) were repeatedly applied  to quantify brain electrical activity, for example under the influence of epileptic processes, during different mental states, or in different sleep stages. Attempts to interpret the absolute numbers obtained in terms of “system complexity”, “attractor dimension”, “chaoticity”, or the like have met harsh and justified criticism . Any such interpretation would require some narrow assumptions (stationarity, predominant low dimensional determinism, etc.) which are expected to be violated when biological time series are analyzed. On the other hand, commonly applied linear or power spectral measures are based on different, but equally narrow assumptions, including (weak) stationarity and Gaussianity. The success of these latter methods is undisputed because of the practical usefulness of the results. In this paper, we will try to take such a pragmatic attitude and compare a number of indicators of human electroencephalogram (EEG) time series with respect to their power in monitoring depth of anesthesia.

Keeping patients at a well defined level of anesthesia is still a difficult problem in clinical practice. If anesthesia is too deep, a decompensation of the cardiovascular system is threatening. When anesthesia is too flat, the patient may wake up. Depth of anesthesia is expected to be reflected in the EEG. Some of the underlying mechanisms are known from animal studies. The EEG is generated by electrical discharges (i.e. postsynaptic potentials) of neurons located near to the surface of the brain (predominantly neurons from the neocortex) . During the awake state, complex tasks can be processed by these cells. Neuronal activity changes during anesthesia as well as during “synchronized” sleep (including all sleep phases except for the rapid-eye-movement-sleep, in which most dreams are experienced). In the middle of the brain (i.e. in the thalamus), neuronal populations start to produce an oscillating activity that increasingly synchronizes neo-cortical neurons. Simultaneously, the excitability of neo-cortical neurons decreases so that the amount of interactions between them decreases as well. Finally, nearly all neo-cortical neurons fall silent. This state will not be reached during normal sleep, it is only observed during coma or deep anesthesia. A small number of thalamic oscillators can even be active during this phase. Before a continuous zero line occurs in the EEG, a so-called burst-suppression phase can appear, where volleys from these thalamic neurons account for a cyclic excitation of the silenced neocortex .

The described features, except for the zero line, can be seen in Fig. 1. Note the marked hysteresis: The burst suppression disappears only about 1 min after the plasma concentration of the anesthetic is reduced. The scenario may be summarized by the hypothesis that with increasing depth of anesthesia, the highly complex activities of the awake state are increasingly superimposed and replaced by more simple oscillators. Of course, such a picture amounts to a gross simplification. The full mechanism will certainly not be well described by a single parameter like “depth of anesthesia”.

## II Validation of measures for depth of anesthesia

In current clinical practice, one or a few channels of EEG are routinely displayed during difficult anesthesias. Since the attending personnel has to monitor several critical parameters (blood pressure, heart rate, etc.), the vast amount of information contained in the EEG must be severely condensed in order to be useful. Only a few numbers may be monitored at a typical intervention time scale. Most pragmatically, a single number should be produced that indicates the instantaneous depth of anesthesia of the patient.

In order to derive and to validate such a measure of depth of anesthesia, the “true” level of anesthetic depth has to be known. One possible definition is to take the concentration of the anesthetic drug in the brain tissue as an indicator. This is reasonable if a monotonic relationship between the concentration of a single drug and the depth of anesthesia can be assumed. However, when multiple drugs are used, the drug interactions may become unpredictable. Unfortunately, even if a single drug is given, its brain concentrations cannot be measured directly in human patients. In order to circumvent this difficulty, in this study the gaseous (volatile) drug Sevoflurane was used which can be employed as single anesthetic. In that case the concentration of the drug in the expired air can be measured, which is proportional to the blood-plasma concentration of Sevoflurane. Then, by considering the brain as a hypothetical effect compartment, the equilibration of plasma concentration $`C_p`$ and the concentration in the brain $`C_e`$ can be modelled in a standard way by

$$\frac{dC_e}{dt}=k(C_pC_e)$$
(1)
where $`k`$ is the constant of equilibration that will vary from patient to patient. Since $`k`$ is unknown it has to be obtained by a fit of this pharmacokinetic model to the data. Equation (1) expresses the hysteretic response noted above to a change in the amount of Sevoflurane administered.

As a further complication, the pharmacodynamic relation between the concentration in the brain as the effect compartment $`C_e`$ and the neuronal effect $`E`$ is not known a priori. A reasonable assumption is that there exists a one-to-one, monotonic relationship. In that case, and if $`C_e(t)`$ and $`E(t)`$ are sampled at discrete times sufficiently far apart to assume statistical independence, this can be expressed by Spearman’s rank correlation. Independence is not a reasonable assumption for consecutive values of a slowly varying quantity. Therefore we cannot use the absolute value of Spearman’s rank correlation for a formal statistical test. However, since the conditions are identical for all the methods we want to compare, significant differences in the correlation may be evaluated.

The aim of our analysis here is to relate the measured $`C_p`$ and some estimated value $`\widehat{E}`$ of the neuronal effect $`E`$, assuming the kinetic model, Eq.(1). Let $`R_i`$ be the rank of $`C_e(i)`$ among all available $`C_e`$’s and $`S_i`$ the rank of $`\widehat{E}(i)`$ among all the corresponding $`\widehat{E}`$’s. An estimate for $`k`$ can then be found by minimizing the sum of the squared differences of ranks:

$$\underset{i=1}{\overset{N}{}}(R_iS_i)^2.$$
(2)
The corresponding rank correlation coefficient

$$\rho =\frac{_i(R_i\overline{R})(S_i\overline{S})}{\sqrt{_i(R_i\overline{R})^2}\sqrt{_i(S_i\overline{S})^2}}$$
(3)
quantifies the success of the fit and can thus be used as a performance index for a method to estimate the neuronal effect.

For our study, 17 patients undergoing an elective surgery were anesthetized with Sevoflurane. As described above, being unable to directly measure the brain concentration $`C_e`$, we measured the concentration in the breathed air at the end of expiration, which is proportional to the blood-plasma concentration $`C_p`$. The neuronal effect $`E`$ is estimated for each patient by four different methods based on bipolar EEG recordings from frontal scalp regions. All EEG recordings used in this study were sampled at 256 Hz with 12 bit resolution. For the time resolved evaluation, a window length of 30 s was used with windows taken every 5 s, that is, with 83% overlap. The total observation period for each patient varied between 13 min and 55 min, median 34 min.

## III EEG based measures of anesthesia depth

In this study, we have compared four methods to derive an indicator for depth of anesthesia from EEG data. Two of these measures are based on the power spectrum, one additionally takes the bispectrum into account, combined with a special feature detection scheme. The fourth is a nonlinear correlation index based on phase space reconstruction and correlation sums. Like in any other problem of EEG characterisation, none of the used methods is derived from first principles nor can any one be fully justified on physiological grounds. Rather, they aim at a useful reduction of the vast amount of information in an EEG trace to one or a few numbers that may be monitored with justifiable effort.

### A Power spectral measures

Apart from visual inspection, the most common method to analyse and quantify time series data in clinical medicine — often the only viable one with the given data quality — is the frequency power spectrum, usually computed using the Fourier transform. Spectral estimation in itself provides only a moderate reduction of information, except in the case that dominant oscillators can be identified by sharp spectral lines and overtones. The art of spectral analysis is to define key parameters that turn a power spectrum into a number, or at most a few essential quantities. For line spectra, the positions and widths of the resonances are obvious candidates. For broad band spectra like the EEG, one may consider the spectral content in certain frequency bands, characteristic frequencies, etc.

Here we use two “classical” spectral parameters  that have been previously applied for this purpose. The spectral edge frequency 95%  indicates the maximal relevant frequencies while the median frequency  gives a rough indicator over the overall typical frequency of the recording. The physiological motivation lies in the expectation that at deeper levels of anesthesia, high frequency components of the surface EEG should be suppressed due to the entrainment by slow oscillations in more central regions of the brain. Therefore, both indices are expected to decrease with increasing depth of anesthesia.

Both measures are based on the time windowed Fourier transform (using 5% cosine tapering of window edges) $`H(f)=_{t=0}^{N1}h(t)\mathrm{exp}(2\pi itf/N)`$ of the digitally sampled EEG time series $`h(t)`$. Inspecting the spectral band between $`f_{\mathrm{low}}=`$ 0.25 Hz and $`f_{\mathrm{high}}=`$ 30 Hz, the accumulated periodogram

$$P(f)=\underset{i=f_{low}}{\overset{f}{}}|H(i)|^2$$
(4)
is used to determine the two spectral parameters, the spectral edge frequency 95%, here denoted by $`f_{\mathrm{se95}}`$, and the median frequency $`f_{\mathrm{median}}`$. The parameter $`f_{\mathrm{se95}}`$ is defined as the frequency at which $`P(f_{\mathrm{se95}})=0.95P(f_{\mathrm{high}})`$ and the median frequency $`f_{\mathrm{median}}`$ as the frequency at which $`P(f_{\mathrm{median}})=0.5P(f_{\mathrm{high}})`$.

### B Bispectral index

A popular parameter in EEG instrumentation for anesthesia is represented by the bispectral index (bis. This composite parameter uses multiple measures of the EEG power spectrum and bispectrum as well as burst-suppression analysis and has been optimized to predict the degree of sedation. Unlike all other EEG parameters, it has been validated in large studies as a measure of sedation in anesthesia . Unfortunately, the precise algorithm is proprietary and has not yet been published. The bispectrum is one possible natural generalisation of the Fourier power spectrum to nonlinear signals. While the ordinary spectrum may be defined as the Fourier transform of the two-point auto-covariance function $`x(t)x(t\tau )`$ of a signal $`x(t)`$, the bispectrum is obtained by a two-dimensional Fourier transform of the three point auto-covariance function $`x(t)x(t\tau _1)x(t\tau _2)`$. The way the bispectral index condenses this two-variate function into a single number seems to have been established heuristically. At least no physiological motivation is given in the literature. In the absence of a published definition, we use the values of the bispectral index as computed by the Aspect 1000 EEG monitor and record it every 5 s.

The reported relative success  of the bispectral index suggests that properties of the EEG beyond the power spectrum are important for distinguishing of different states of anesthesia. We do not want to enter into a discussion here of what nature the extra structure is. It is very difficult to disentangle the effects of nonlinearity and non-stationarity, and, if nonlinearity is found, to what extent it is of dynamical nature or if it merely reflects the complicated transduction properties of the intervening tissue. From a physiological point of view, the situation is rather involved. The individual neurons clearly show nonlinear behaviour in their action potentials, but the law of large numbers suggests that this structure might be averaged out in surface EEG traces. Non-stationarity is very prominent in anesthesia, in particular in the burst suppression phase, where episodes of low amplitude activity are interrupted by short outbreaks of large amplitude waves. Spectral indicators are known to give spurious values in the presence of burst suppression. Therefore, the bispectral index algorithm contains a pattern recognition scheme to detect burst suppression, which is then treated separately.

### C Nonlinear correlation index $`D^{}`$

Early claims of low dimensional strange attractors in brain signals have been identified as spurious due to the lack of certain precautions when using nonlinear analysis tools. (See  for further discussion and references.) This does not exclude that tools and methodologies derived in the context of dynamical systems, or chaos theory, may be still be useful when employed with care. In particular, interpretations in terms of fractality, chaoticity, complexity, and the like, have to be attempted with utmost discretion.

In that spirit, we have adapted a prescription for an overall index of (nonlinear) coherence that has been found powerful for anticipating epileptic seizures from implanted electrode recordings  as well as in epilepsy models on a cellular level . This index, which we will call $`D^{}`$ contains many ingredients familiar from the Grassberger–Procaccia algorithm for the correlation dimension , whence the $`D`$. The  is there to remind us that there are important differences that preclude an interpretation in terms of fractal dimensions, number of degrees of freedom, etc.

Starting form a scalar time series $`s_n,n=1,\mathrm{},N`$, we first form delay embedding vectors as usual: $`𝐬_n=(s_{n(m1)\tau },\mathrm{},s_{n\tau },s_n)`$. The correlation sum $`C_m(r)`$ is then defined as usual by

$$C_m(r)=\alpha \underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{i\mathrm{\Delta }n}{}}\mathrm{\Theta }\left(r𝐬_i𝐬_j\right)$$
(5)
where $`\mathrm{\Theta }`$ is the step function and $`\alpha =2/(N\mathrm{\Delta }n)(N\mathrm{\Delta }n1)`$ is a normalisation constant. For the estimation of the correlation dimension, one would take $`\mathrm{\Delta }n`$ large to exclude temporal correlation and study the slope of a double logarithmic plot of $`C_m(r)`$ versus $`r`$ and look for a scaling region, that is, a range of values of $`r`$ where

$$C_m^{}(r)=\frac{d\mathrm{log}(C_m(r))}{d\mathrm{log}(r)}$$
(6)
is constant. For high enough embedding dimension $`m`$, that constant would be an estimator of the correlation dimension $`D_2`$. There is ample literature on using the correlation integral for dimension estimation .

Except for specific pathologies, we do not expect this procedure to arrive at a reasonable, finite value for $`D_2`$ when applied to the EEG. Since we cannot hope to estimate a proper dimension in the first place, there is no particular reason to concentrate on geometrical correlations only, or to aim at a proper phase space embedding. In principle, we could instead optimize parameters for the specific purpose at hand. However, a total of 17 patients is too small a population to be split into a training set and a test set and we are not interested in in-sample results that may be due to over-fitting. Therefore, we copy the parameter settings and operational procedures from previous studies on epilepsy prediction  without further optimization. Using unit delay, a range of embedding dimensions $`m=10,\mathrm{},25`$ was studied and only immediate temporal neighbours with $`ij<\mathrm{\Delta }n=10`$ were excluded. Each segment of 30 s was digitally band-pass filtered (0.25 Hz–30 Hz) before computing $`C_m(r)`$. The operational procedure of determining a “plateau” value $`D^{}`$ of $`C_m^{}(r)`$ is also described in Ref. . It is worth noting that this automatic procedure is designed to exclude both, large amplitude artefacts (large $`r`$) and small amplitude fluctuations (small $`r`$). Some conjecture on the reasons for the success of this procedure will be offered in the discussion section below.

## IV Clinical Results

In all 17 patients the following protocol was carried out. The concentration of Sevoflurane was at first slowly increased. As soon as a burst suppression pattern or even a zero line occurred in the EEG, Sevoflurane concentration was decreased again. Two typical examples of corresponding EEG time series and estimated measures of anesthesia depth are shown in Fig. 2. In case B, immediately a zero line occurs in the EEG rather than a burst suppression pattern. It is important to remember that the actual neuronal anesthetic effect is not simply proportional to the plasma concentration $`C_p`$ shown in Fig. 2 but follows hysteretically through the relaxation phenomenon we modelled by Eq. (1). Generally speaking, all four measures studied are expected to show a negative correlation with the depth of anesthesia.

In both cases, a slight decline of $`D^{}`$, but a small increase of bis precedes the burst suppression (resp. zero line) phase. The reaction of bis on these abrupt changes of EEG patterns is delayed. Depending on frequency characteristics of the background noise, spectral measures either decrease (A) or increase (B) during phases with zero-line EEG. For a quantitative evaluation, we need to estimate the effective brain concentration $`C_e`$ of Sevoflurane. For this purpose, we assume the pharmacokinetic model given by Eq. (1) and determine an equilibration-constant $`k`$ individually for each case by minimising the sum of the squared differences of ranks between $`C_e`$ and $`E`$, c.f. Eq. (2). The values for $`k`$ that were found  are summarized in Table I.

Compared to the variation from patient to patient, there were no dramatic differences among the results for the different measures of anesthesia. It is important to notice that the fit for $`k`$ is only important in order to establish a monotonic relationship between an EEG measure and assumed depth of anesthesia. The parameter $`k`$ describes the time course of the equilibration of anesthetic concentration between blood plasma (assumed to be equal to the concentration in the air at the end of expiration) and a hypothetical effect compartment (i.e. the brain or parts of it). It will be influenced by factors like brain perfusion, patients’ age, amount of body-fat, etc. However, since the value of $`k`$ is no longer needed once a monotonic relationship is assured, these influencing factors need not be under close control for this approach.

Based on the optimal values for $`k`$, for each case Spearman’s rank correlation coefficient $`\rho `$ was computed (cf. Eq. 3). Table II shows that the highest correlation coefficients were achieved for $`D^{}`$. In order to assess the statistical significance of this finding, we have performed a nonparametric test for the hypothesis that $`D^{}`$ performs better than bis. Since the populations of values where obtained on the same patients, we evaluated the Wilcoxon matched-pairs signed-rank test (see Fig. 3) and found that the correlation is higher for $`D^{}`$ with $`p<0.05`$.

For clinical applicability, it is not sufficient to have a good correlation for each patient individually. Rather, a monotonic relationship has to hold across patients. In order to illustrate the relevance of the relation between the brain concentration of Sevoflurane $`C_e`$ and its estimated EEG effect $`E`$ for the whole group of patients, for each investigated EEG measure a scatterplot of $`C_e`$ against $`E`$, including results from all patients, is shown in Figs. 4 and 5. The corresponding rank correlation coefficients $`\rho `$ are -0.06 for $`f_{\mathrm{median}}`$, -0.33 for $`f_{\mathrm{se95}}`$, -0.63 for bis, and -0.72 for $`D^{}`$. Due to the serial correlations within each trial, we don’t attempt a formal test of the hypothesis that a (negative) correlation exists. Visual inspection as well as the value for $`\rho `$ suggests that $`f_{\mathrm{median}}`$ is essentially useless as an indicator of depth of anesthesia in the investigated group of patients Also for the other spectral measure, $`f_{\mathrm{se95}}`$, the correlation is rather unreliable. It is conceivable that in a study that excludes burst-suppression EEG or with additional preprocessing (e.g., by introducing a burst-suppression detector) the situation may be more favourable. Both bis and $`D^{}`$ show a certain trend to decrease with increasing depth of anesthesia. Of all the studied measures, we find that $`D^{}`$ comes closest to an overall monotonic relation to the depth of anesthesia.

## V Discussion

In a relatively small but focused clinical study, we have compared four different methods to quantify depth of anesthesia by numerical analysis of EEG tracings. All four methods were defined outside the study whence the results can be considered to be out-of-sample. We have used the available data to answer the question which of the four measures comes closest to a monotonic relationship to the level of sedation. While for both of the two nonlinear measures, bis and $`D^{}`$, such a relationship seems to exist, the correlation is strongest for $`D^{}`$, despite the fact that bis had been specifically designed for this purpose in previous work.

With the limited data base available, any interpretation of the findings remains speculation. The success of $`D^{}`$ in previous studies on the anticipation of epileptic seizures has been discussed in Ref. , where it was pointed out that the connection with the epileptic alterations is lost if the temporal correlations are fully excluded from the correlation sum, as one would have to do for a dimension estimation . It can be suggested that geometrical correlations are not essential but the ability to select structures – temporal or geometric – by amplitude, or length scale in phase space, makes the correlation sum superior to other nonlinear statistics.

In earlier studies, $`D^{}`$ has been found helpful for the anticipation of epileptic seizures  from intra-cranial EEG recordings. It has also been shown  that even on a cellular level $`D^{}`$ could quantify the complexity of the synaptic input of a single neuron when inter-ictal epileptiform activity was induced in brain slices. The degree of autonomous activity, and hence the number of degrees of freedom of a neuronal network, is thought to decrease during the development of epileptiform activity. Indeed, a reduction of $`D^{}`$ extracted from intracellular recordings of single neurons was observed prior to visually detectable alterations .

In the common scenario of anesthesia, the number of degrees of freedom of human brain activity is lowered — as compared to the awake state — when nearly all neo-cortical neurons fall silent and only a central (thalamic) oscillator is able to evoke a cyclic excitation of the silenced neocortex . However, the number of degrees of freedom may still be huge and is probably not accessible from an EEG time series . It is even the more remarkable that $`D^{}`$ correctly ranks the qualitatively different EEG pattern of burst suppression as corresponding to the deepest level of anesthesia. At the onset of burst suppression, the EEG signature changes discontinuously with the estimated drug concentration, whence the relatively sharp transition in $`D^{}`$ between the preceding deep sleep phases and burst suppression in Fig. 5. Still, from a biological point of view it remains unclear which electrophysiological properties contribute to the monotonic decline of $`D^{}`$ under anesthesia — or epileptic activity.

In summary, extraction of $`D^{}`$ seems to be able to improve the quantification of depth of anesthesia from brain electrical activity, at least when Sevoflurane is used as anesthetic drug. For clinical applicability, other anesthetics have to be investigated. The size of this study is too small to admit any optimisation of the algorithm. It is conceivable, however, that the clinical reliability can be enhanced, and the computational burden relieved, by a more thorough understanding of which structures are really picked up by the algorithm.


format: text

text:
# References

In the last few years there has been a noticeable interest in the study of the thermodynamics and statistical mechanics of anomalous systems exhibiting nonextensivity . A system is extensive if its energy and entropy, as functions of intensive internal parameters (e.g. temperature), grow linearly with the size of the system (i.e. with $`N`$, the number of its microscopic components). Nonextensivity can be brought into scene by long-range interactions. From the static point of view, nonextensivity means that thermodynamic quantities like the internal or the free energies per particle are not constant with $`N`$, but rather diverge in the thermodynamic limit $`N\mathrm{}`$. From the dynamical point of view, one might think that long-range interactions induce persistent spatial and temporal correlations leading to the breaking of standard mixing and ergodicity properties, hence to a possible violation of the usual Boltzmann-Gibbs (BG) statistics. Classical Hamiltonian systems with many degrees of freedom are paradigmatic in the discussion of the foundations of equilibrium statistical mechanics. The maximal Lyapunov exponent (MLE) is a well established indicator of chaos. If it is positive, the system will generically be strongly chaotic and will satisfy the standard ergodic hypothesis attained through exponentially quick mixing. If it is instead zero non ergodicity could emerge, giving origin to less-than-exponential (typically power-law) mixing, hence to anomalous thermostatistical behavior. It is then interesting to study maximal Lyapunov exponents in systems with long-range interactions, where nonextensive thermodynamical and dynamical anomalies are expected. In a previous work  the maximal Lyapunov exponent of a one-dimensional system of $`N`$ planar rotators coupled with interactions decaying as the inverse power $`\alpha `$ of their distances (see model (1) below) has been studied as a function of $`N`$. For total energy above some threshold, the maximal Lyapunov exponent is proportional to $`1/N^\kappa `$, where $`\kappa (\alpha )`$ is a function which goes from $`\kappa (0)>0`$ to zero while $`\alpha `$ increases from zero to 1 (long range forces) and remains zero for $`\alpha >1`$ (short range forces). In this letter we want to investigate if the exponent $`\kappa `$, describing the weakening of the mixing properties of the dynamics in the long-range-interacting regime, is universal. In particular, we check the hypothesis that $`\kappa (\alpha ,d)`$ is a universal function of $`\alpha /d`$, where $`d`$ is the dimensionality of space.

Let us consider a simple model of $`N`$ planar rotators (XY spins) placed at the sites $`i=1\mathrm{}N`$ of a square and a cubic lattice (i.e., dimension $`d`$ = 2 and 3). The rotators have unit moment of inertia, angular momentum $`L_i`$ and position specified by the angle $`\theta _i[0,2\pi ]`$. Rotations occur inside an arbitrary reference plane. The model is described by the following classical hamiltonian in the conjugate canonical pairs $`\{\theta _i,L_i\}`$:

$$H=K+V=\frac{1}{2}\underset{i}{}L_i^2+\frac{1}{2}\underset{ij,ij}{}\frac{1\mathrm{cos}(\theta _i\theta _j)}{r_{ij}^\alpha };$$
(1)
where $`r_{ij}`$ is the distance between lattice site $`i`$ and $`j`$ with periodic boundary conditions in the form of nearest image convention; $`\alpha 0`$ is a variable parameter specifying the range of the interaction, which is defined long-range if $`\alpha d`$ and short-range if $`\alpha >d`$.

This system is nonextensive in the case of long-range interactions . To verify that, let us consider as an example a quantity related to the internal energy per rotator at $`T=0`$ (i. e., when $`K=0`$ and all rotators are equally oriented), namely the integral $`_1^{\mathrm{}}𝑑rr^{d1}r^\alpha `$, which converges for $`\alpha >d`$ and diverges for $`0\alpha d`$. In the latter case, the thermodynamics of the model is not well defined in the usual sense. However, it is possible to define energy quantities that are bounded in the thermodynamic limit if they are properly scaled (see ). It is interesting to remark that, in system (1) for $`0\alpha d`$, not only equilibrium (static) quantities, but also a quantity related to the local instability of the dynamics such as the maximal Lyapunov exponent might diverge with $`N`$ , and needs to be properly rescaled by a function of $`N`$ if one wants to evaluate the degree of chaoticity as a function of the size of the system.

An alternative procedure is to force the system to be extensive via a largely used mean-field-like rescaling  of the potential energy in (1) which then becomes:

$$H^{}=K^{}+V^{}=\frac{1}{2}\underset{i}{}L_i^2+\frac{1}{2\stackrel{~}{N}}\underset{ij,ij}{}\frac{1\mathrm{cos}(\theta _i\theta _j)}{r_{ij}^\alpha }$$
(2)
The rescaling factor $`\stackrel{~}{N}`$ in front of the potential energy $`V`$ of hamiltonian (2), is a function of $`\alpha ,d,N`$ and of the geometry of the lattice and has to be chosen to guarantee the existence of a bounded energy density $`U=<H^{}/N>`$ in the thermodynamic limit.

The model (1) has been introduced in  on a one-dimensional lattice and its thermodynamics has been studied by means of simulations in . When $`\alpha =0`$ and $`\stackrel{~}{N}=N`$ the model is also called the Hamiltonian Mean-Field XY Model (HMF) , whose dynamics and thermodynamics have been extensively studied . In this case the underlying lattice has no meaning and any dependence from its dimension $`d`$ disappears. It is found that the HMF model has a phase transition from a ferromagnetic state to a paramagnetic one at the critical temperature $`T_c=0.5`$ corresponding to a critical energy density $`U_c=0.75`$ ($`k_B=1`$). In this letter we use

$`\stackrel{~}{N}1+d{\displaystyle _1^{N^{1/d}}}𝑑rr^{d1\alpha }=\{\begin{array}{cc}\hfill \frac{N^{1\alpha /d}\alpha /d}{1\alpha /d}& \text{if }\alpha d\hfill \\ \hfill 1+\mathrm{ln}N& \text{if }\alpha =d\hfill \end{array}`$ (5)
whose properties are shown in figure 1. (Unity has been added to the integral in the definition of $`\stackrel{~}{N}`$ in order to have the convenient limit $`\stackrel{~}{N}1`$ when $`\alpha /d\mathrm{}`$.) Using a similar rescaling factor it was shown in  that model (2) and HMF model have the same thermodynamics, and in particular the same critical energy density.

Chaoticity properties and maximal Lyapunov exponent (MLE) of HMF model have been studied by Latora, Rapisarda and Ruffo . They calculated the curve MLE versus $`U`$ numerically for various $`N`$ and analytically for $`N\mathrm{}`$; MLE goes to 0 both when $`U`$ tends to $`0`$ and when $`U`$ tends to $`\mathrm{}`$. In these two limits the system becomes integrable, being represented by a set of harmonic oscillators around the ferromagnetic fundamental state of aligned rotators in the former case and by free rotators in the latter. Finite-size effects are present for $`U>U_c`$ where MLE vanishes in the thermodynamic limit as $`N^{1/3}`$. For fixed $`UU_c`$ the curve is smooth and presents a positive maximum at the critical energy $`U_c`$, where therefore in the thermodynamic limit exists a finite discontinuity.

The behavior of MLE in the one-dimensional model (2) has been studied, as said above, in , as a function of $`N`$ and $`\alpha `$ at a fixed high energy (above $`U_c`$). It has been shown that it goes to 0 in the thermodynamic limit when $`\alpha d`$, like in the HMF model, but it remains constant when $`\alpha >d`$. The way the exponent vanishes is still of the form $`N^\kappa `$ where now $`\kappa `$ is a function of $`\alpha `$. The curve $`\kappa (\alpha )`$ starts from the value $`1/3`$ at $`\alpha =0`$ and goes to 0 when $`\alpha 1^{}`$.

In the present work we extend this study to cases $`d>1`$. We still expect that when $`\alpha <d`$ and $`U>U_c`$ MLE vanishes in the thermodynamic limit with a power law of the type $`N^\kappa `$ with $`\kappa =\kappa (\alpha ,d)`$ satisfying $`\kappa (\alpha =0,d)=1/3`$  and vanishing when $`\alpha d^{}`$. Two questions seem quite relevant. Is this curve $`\kappa (\alpha ,d)`$ universal in the sense that it only depends on $`\alpha /d`$, or it also depends on the dimensionality $`d`$ of the system ? Moreover, is this curve independent from the energy density $`U`$ at which it is computed? In what follows we will positively answer to both questions.

We have simulated the constant energy dynamics of the nonextensive hamiltonian (1) and computed the maximal Lyapunov exponents.

Hamiltonians $`H`$ (eq. (1)) and $`H^{}`$ (eq. (2)) are strictly connected. If we multiply $`H^{}`$ by $`\stackrel{~}{N}`$, the potential energy $`V^{}`$ becomes $`V`$. Then, if we rescale time as $`t^{}=t\sqrt{\stackrel{~}{N}}`$, and using the fact that $`L_i`$’s involve first order time derivatives, we obtain that also the kinetic energy $`K^{}`$ becomes $`K`$. Finally we get $`H=\stackrel{~}{N}H^{}`$.

Integrating the equations of motion of $`H`$ is identical to integrating those of $`\stackrel{~}{N}H^{}`$ but for times longer by a factor $`\sqrt{\stackrel{~}{N}}`$. All thermodynamic quantities that can be obtained as temporal averages of time dependent observables are not influenced by the rescaling of time. In particular, once one has $`T,V,U`$ from $`H`$ simulations, one can get the corresponding values for $`H^{}`$ simply dividing by $`\stackrel{~}{N}`$. The evaluation of maximal Lyapunov exponents $`\lambda _{max}^{}`$ of hamiltonian $`H^{}`$ can also be obtained through rescaling of the one measured for $`H`$, noted $`\lambda _{max}`$ hereafter.

In order to calculate maximal Lyapunov exponents one must consider the limit,

$$\lambda _{max}=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{d(t)}{d(0)}=\underset{t\mathrm{}}{lim}\lambda (t)$$
(6)
with $`d(t)=\sqrt{_i(\delta \theta _i)^2+(\delta L_i)^2}`$ being the metric distance calculated from infinitesimal displacements at time $`t`$ obtained in turn through the double integration of both the normal equations of motion

$`\dot{\theta _i}`$ $`=`$ $`L_i`$ (7)
$`\dot{L_i}`$ $`=`$ $`{\displaystyle \underset{j,ji}{}}{\displaystyle \frac{\mathrm{sin}(\theta _j\theta _i)}{r_{ij}^\alpha }}`$ (8)
and the linearized ones

$`\dot{\delta \theta _i}`$ $`=`$ $`\delta L_i`$ (9)
$`\dot{\delta L_i}`$ $`=`$ $`{\displaystyle \underset{ji}{}}{\displaystyle \frac{\mathrm{cos}(\theta _j\theta _i)}{r_{ij}^\alpha }}(\delta \theta _j\delta \theta _i)`$ (10)
(we report the ones derived from $`H`$, the ones for $`H^{}`$ are exactly the same with the substitution $`r_{ij}^\alpha r_{ij}^\alpha \stackrel{~}{N}`$). Because of formula (6) which involves $`1/t`$, time rescaling yields $`\lambda _{max}^{}=\lambda _{max}/\sqrt{\stackrel{~}{N}}`$.

We computed $`\lambda _{max}`$ for typical values of $`N30`$ to $`4000`$ and a fixed energy density $`U=5.0`$ above the critical one $`U_c=0.75`$. We made simulations on a square and simple cubic lattices (i.e., $`d=2,3`$) with unitary lattice step, for various values of the $`\alpha /d`$ ratio.

Simulations started from an high temperature initial state: $`\theta _i`$’s are randomly extracted from a uniform distribution in $`[0,2\pi ]`$, momenta $`L_i`$’s are extracted from a uniform distribution in $`[0.5,0.5]`$ and then translated and rescaled to have zero total momentum and the desired total energy. We used the velocity-Verlet algorithm  with a time-step chosen to have a relative energy conservation of $`10^4`$ or better. Length of simulations were chosen looking at the asymptotic behavior of quantity $`\lambda (t)`$, eq. (6), where $`d(0)`$ is randomly chosen (see ). We report in figure 2, the curves $`\lambda _{max}^{}(N)`$ for various $`\alpha `$ in dimensions $`d=2`$ and $`d=3`$.

These curves extend those of fig. 3 in . In both cases we observe that, for growing $`N`$, if $`\alpha >d`$ then $`\lambda _{max}^{}`$ is positive and constant, whereas if $`\alpha <d`$ the maximal Lyapunov exponent tends to zero.

We fitted the data with the following functional form: $`\lambda _{max}^{}N^\kappa `$; $`\kappa `$ being the slope in the log-log plots of fig. 2.

In figure 3 we collect the slopes $`\kappa `$ as a function of the ratio $`\alpha /d`$, for $`d=1`$ (data from ), $`d=2`$ and $`d=3`$. Remarkably, through the simple scaling $`\alpha /d`$ we show that the $`N`$ dependence of the rescaled maximal Lyapunov exponent of model (1) is universal. All curves start at $`\alpha =0`$ from a value close to $`1/3`$, analytically and numerically found for the HMF model ; when $`\alpha /d`$ increases from zero to unity all our data collapse into a single curve, and then remain zero for $`\alpha /d`$ greater than unity. All of our data are roughly fitted by the heuristic expression $`\kappa =[1(\alpha /d)^2]/[3+(\alpha /d)^2/2]`$.

Extrapolating for $`N`$ going to infinity the results in fig. 2, it is interesting to note that the constant $`\kappa `$ behavior for $`\alpha >d`$ is consistent with the results found for the $`\alpha =\mathrm{}`$ (only first-neighboring coupling) model .

The universal function $`\kappa `$ shown in fig. 3 does not appear to depend on the energy density, provided this is greater than the critical one. For the case $`\alpha =0`$ this has been proved analytically in  and numerically in . We show in fig. 4 that this is also true for $`\alpha =0.8`$ and $`d=2`$, the asymptotic slope of $`\lambda _{max}^{}`$ being sensibly the same for the two energy densities $`U=5`$, $`U=8`$.

Summarizing our results, we can say that the classical hamiltonian $`H`$ is non extensive for $`0\alpha d`$ and extensive for $`\alpha /d>1`$. Through a $`\stackrel{~}{N}`$ rescaling of its potential energy it can artificially be made extensive, thus merging with a long standing tradition in mean-field-like approaches. Above some critical value for the (rescaled) total energy per particle, the weakening of the mixing properties is universal; more precisely, the maximal Lyapunov exponents scale as $`N^\kappa `$, $`\kappa `$ being a function only of the ratio $`\alpha /d`$ (and not $`d`$ nor $`U`$). In fact, it seems reasonable to expect that this $`\kappa (\alpha /d)`$ function is model-independent. More precisely, for total energies high enough, it is expected to have one and the same functional form $`\kappa (\alpha /d)`$ for all classical models whose potential is nowhere singular and whose attractive tail decays as $`1/r^\alpha `$. Consistently with this speculation, some specific gravitational and fluid models  do exhibit $`\kappa (0)=1/3`$, in agreement with the present results for the $`XY`$ models under analysis. An universal $`\kappa (\alpha /d)`$ fits in fact very well the generic expectation for such models within nonextensive statistical mechanics. Indeed, the available studies of hamiltonian systems including long range interactions provide a variety of indications that a crossover occurs, as a function of time, in the thermostatistical properties of the system. For times below a crossover time $`\tau (N)`$, several properties are anomalous (diffusion, velocity distribution and others), whereas they become normal for times much larger than $`\tau `$. The fact that $`\tau `$ seems to diverge with $`N`$ makes the anomalous regime to become the only observable one in the thermodynamic limit. It is precisely this regime which is focused on within that generalized formalism. For the specific case of the mixing properties, the present results suggest that a less-than-exponential mixing takes place until a crossover time after which the mixing becomes exponential, i.e., the usual one. This crossover time might well scale like $`1/\lambda _{max}^{}`$ (or a positive power of it). Therefore, since $`\lambda _{max}^{}`$ vanishes in the limit $`N\mathrm{}`$, we expect the crossover time to diverge, thus emerging the anomalous regime we were just discussing.


text:
# THE PARKER INSTABILITY IN A THICK GASEOUS DISK II: NUMERICAL SIMULATIONS IN 2D

## 1 INTRODUCTION

The interstellar disk in our Galaxy is thicker and more magnetized than previously assumed, and the transition between the gaseous disk and the halo has a complex magnetic field structure and extended gas components, including the ionized Reynolds layer (Hoyle & Ellis 1963; Reynolds 1989). This type of extended disk also seems to be common in other spiral galaxies, and thick diffuse gas layers, with worms and dusty filaments, have been observed, for instance, in NGC 891, NGC 5775, and NGC 4302 (Dettmar 1990; Rand, Kulkarni & Hester 1990; Rand 1995; Howk & Savage 1997). Similarly, the observed magnetic fields in edge–on galaxies can be traced well above the main gaseous disk (Hummel & Beck 1995).

Boulares & Cox (1990; hereafter BC) have incorporated the observed extended gas layers, the gradient of the magnetic field, and the cosmic ray pressure in a model for the thick disk of the Galaxy (see also Kalberla & Kerp 1998). The existence of magnetized gas layers with large scale heights has important consequences for the structure of the gaseous disk, and its response to large scale perturbations (see Cox 1990; Martos 1993; Franco et al. 1995). For instance, the total interstellar pressure in a disk with extended layers is substantially larger than that for thin disk models. Also, the interaction between the disk and the halo is more complex: spiral density waves can produce flows, shocks, and even star formation at high galactic latitudes (Martos & Cox 1998; Martos et al. 1999), whereas the penetration of high–velocity clouds into the disk is severely reduced by magnetic pressure and tension (e. g. Santillán et al. 1999). The stability of the system is also modified and, depending on the temperature and $`𝐁`$-field distributions, the equilibrium configuration of a thick disk can be Parker unstable. For the observed conditions at the solar neighborhood (BC), the linear stability analysis of a multi-component isothermal disk indicates that the fastest growing undular mode has a wavelength of the order of 3 kpc (Kim et al. 2000; hereafter Paper I). This is a factor of 8 larger than the one derived for the thin gaseous disk, and illustrates the variations that can be expected when considering the extended gas layers.

In Paper I we derived both the dispersion relations and the final equilibrium configuration in two-dimensions (2D) of the undular mode of the Parker instability. This mode generates high-density condensations at the magnetic valleys, which are eventually stabilized by magnetic tension and the system settles down into a new final equilibrium. The linear analysis provides the initial growth rates and inter–distance scales of the condensations, and the final equilibrium configuration (obtained by solving the magneto-hydrostatic equations with the flux freezing condition) provides their final steady-state structure. To complete the study, here we present 2D numerical simulations for the evolution of the instability. This enables us to verify the results of the linear analysis, and to follow the details of the non-linear regime until the final equilibrium is reached. As a check on the reliability of the numerical results, the numerical experiments have been performed with two different MHD numerical codes, ZEUS-3D and TVD. The results obtained with both codes are very similar, and confirm the growth rates and wavelengths of the fastest growing modes derived in Paper I. The plan of the present paper is as follows. The details of the numerical work are described in the next section. The results are described in Section 3, and a brief summary and discussion are given in Section 4.

## 2 NUMERICAL SETUP AND CODES

Since the first one-dimensional simulation of the problem (Baierlein 1983), there has been a number of numerical studies of the evolution of the Parker instability under different conditions. For instance, Matsumoto et al. (1988, 1990) and Matsumoto & Shibata (1992) have explored the nature of the undular and mixed modes for accretion disk environments with 2D and 3D simulations. Also, using the same numerical code, Shibata et al. (1989a, 1989b), Kaisig et al. (1990) and Nozawa et al. (1992) have used the Parker instability as the driving mechanism for emergent flux tubes in the Sun, and explained several features of the solar activity. In a different type of study, addressing the question of giant molecular cloud formation, Basu et al. (1997) and Kim et al. (1998) have performed 2D and 3D simulations of the instability in a thin gaseous disk under the influence of a simplified uniform gravity. Here we follow the evolution of the undular instability in a multi-component, isothermal and magnetized gaseous disk, with the properties of the solar neighborhood.

### 2.1 A Conservative Form of the Isothermal MHD Equations

We use two-dimensional Cartesian coordinates $`(y,z)`$, whose directions are azimuthal and perpendicular to the Galactic disk, respectively. When the vertical gravitational acceleration has only $`z`$-component, $`𝐠=g(z)\widehat{z}`$, the conservative form of the isothermal MHD equations is (see Paper I)

$$\frac{𝐪}{t}+\frac{𝐅_𝐲}{y}+\frac{𝐅_𝐳}{z}=𝐒,$$
(1)
$$𝐪=\left(\begin{array}{c}\rho \\ \rho v_y\\ \rho v_z\\ B_y\\ B_z\end{array}\right),$$
(2)
$$𝐅_𝐲=\left(\begin{array}{c}\rho v_y\\ \rho v_y^2+a^2\rho +\frac{1}{8\pi }(B_z^2B_y^2)\\ \rho v_yv_z\frac{1}{4\pi }B_yB_z\\ 0\\ v_yB_zv_zB_y\end{array}\right),𝐅_𝐳=\left(\begin{array}{c}\rho v_z\\ \rho v_yv_z\frac{1}{4\pi }B_yB_z\\ \rho v_z^2+a^2\rho +\frac{1}{8\pi }(B_y^2B_z^2)\\ v_zB_yv_yB_z\\ 0\end{array}\right),$$
(3)
$$𝐒=\left(\begin{array}{c}0\\ 0\\ gv_z\\ 0\\ 0\end{array}\right),$$
(4)
where $`a`$ is the isothermal speed for all interstellar components (i. e. the velocity dispersion of all gas components is set to the same value), and the rest of the symbols have their usual meaning.

### 2.2 Initial Condition: The Thick Galactic Disk

As the initial condition for our numerical experiments, we introduce the Galactic disk model that was originally proposed by Martos (1993). The model is the magnetohydrostatic equilibrium of a system composed of five interstellar gas layers, with a horizontal magnetic field and a realistic gravitational acceleration. Since the detailed building procedure of the model is given in Santillán et al. (1999) and Paper I, here we give only a brief description of the model.

The two building blocks of the model are the density distribution of the multi-component gaseous disk (BC),

$`n(z)`$ $`=`$ $`0.6\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(70\text{pc})^2}}\right]+0.3\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(135\text{pc})^2}}\right]+0.07\mathrm{exp}\left[{\displaystyle \frac{z^2}{2(135\text{pc})^2}}\right]`$ (5)
$`+`$ $`0.1\mathrm{exp}\left[{\displaystyle \frac{|z|}{400\text{pc}}}\right]+0.03\mathrm{exp}\left[{\displaystyle \frac{|z|}{900\text{pc}}}\right]\mathrm{cm}^3,`$
and a simple fitting formula (Martos 1993) to the gravitational acceleration at the solar neighborhood (Bienaymé, Robin & Crézé 1987),

$$g(z)=8\times 10^9\left[10.52\mathrm{exp}\left(\frac{|z|}{325\text{pc}}\right)0.48\mathrm{exp}\left(\frac{|z|}{900\text{pc}}\right)\right]\mathrm{cm}\mathrm{s}^2.$$
(6)
The mass density stratification, including 10% He by number, is $`\rho (z)=1.27m_\mathrm{H}n(z)`$, where $`m_\mathrm{H}`$ is the mass of a hydrogen atom. The total (gas plus magnetic) pressure, $`P_0(z)`$, is obtained from the integral of $`\rho (z)`$ and $`g(z)`$ with the boundary condition, $`P_0(z=10\mathrm{kpc})=0`$. Assuming that the gas is isothermal and setting the midplane value of the magnetic field, $`B_0(0)5\mu `$G (BC; Heiles 1996), we deduce the initial field stratification, $`B(z)`$, from the distribution of the total pressure $`P_0(z)`$. In contrast to thin disk models, where the magnetic-to-thermal pressure ratio $`\alpha `$ is usually held constant, our isothermal thick disk has a $`z`$-dependent $`\alpha `$ that increases with height (see Paper I). Thus, the high-latitude gas is largely supported by the magnetic field, and is unstable to the undular mode.

The isothermal thick disk model is specified by two parameters: the effective scale height, $`H_{\mathrm{eff}}=166`$ pc, and the isothermal sound speed $`a=8.4`$ km s<sup>-1</sup> (see Paper I). This scale height is simply defined as half of the total disk column density divided by the density at the midplane, and the sound speed is determined by the midplane values $`P_0(0)`$ and $`B_0(0)`$. These two parameters, $`H_{\mathrm{eff}}`$ and $`a`$, define our units of length and velocity. The time unit is defined as $`H_{\mathrm{eff}}/a=1.9\times 10^7`$ yr and, unless otherwise stated, all physical quantities are normalized to these three units.

### 2.3 Perturbations

The linear stability analysis in Paper I showed that the above equilibrium state is unstable against the undular instability when the perturbation wavelength is longer than about 1.5 to 1.8 kpc, depending on the location of the nodal point (i. e. the location where $`v_z=0`$; see below). These critical wavelength values correspond to nodal points located between 5 and 1.5 kpc above midplane, respectively. Thus, in order to initiate the instability in the numerical experiments, the equilibrium state is perturbed with appropriate wavelengths.

We have used three kinds of perturbations. The first type is the even parity midplane symmetric (MS) perturbation, which forms structures distributed symmetrically with respect to the midplane,

$$v_z=C\mathrm{sin}(\frac{2\pi y}{\lambda _y})\mathrm{sin}(\frac{2\pi z}{\lambda _z}),$$
(7)
where $`C`$ is the amplitude of the vertical velocity perturbation, and $`\lambda _y`$ and $`\lambda _z`$ are the horizontal and vertical wavelengths of the perturbations. The second type is the odd parity midplane antisymmetric (MA) perturbation,

$$v_z=C\mathrm{sin}(\frac{2\pi y}{\lambda _y})\mathrm{cos}(\frac{2\pi z}{\lambda _z}).$$
(8)
The initial amplitudes for both MS and MA perturbations are usually set to $`C=10^3a`$.

The last type is a random velocity perturbation, without any preferred wavelength or parity. The standard deviation of each velocity component is here set to 10$`{}_{}{}^{4}a`$.

### 2.4 Boundary Conditions

In deriving the dispersion relations in Paper I, perturbations were allowed to vary periodically along the horizontal axis. The boundary conditions at midplane determined the parity of the perturbation (i. e. $`v_z=0`$ at $`z=0`$ for the MS perturbation, and $`dv_z/dz=0`$ at $`z=0`$ for the MA case). However, regardless of parity, the vertical velocity perturbation at the upper (or lower) boundary was always set equal to zero. Therefore, the upper boundary becomes a nodal point.

In the numerical experiments, for consistency with the boundary conditions used in the linear analysis, we use a periodic condition at the $`y`$-boundaries of the computational domain, and a reflecting condition at the $`z`$-boundaries. The reflecting condition imposes $`v_z=0`$ on the vertical boundaries, which means that the boundaries mimic the nodal points of the linear analysis. Given that no flow is allowed through the $`z`$-boundaries and any outflow at one horizontal boundary is always compensated by an identical inflow at the other boundary, the mass in the computational domain is conserved at all times. This allows us to easily track the evolution of the energy content of the system.

The even parity perturbation is symmetric with respect to the midplane and has a nodal point at $`z=0`$. Thus, $`\lambda _z`$ in equation (7) should be set equal 2$`z_{\mathrm{node}}`$, and it is sufficient to use only one (the upper) hemisphere as the computational domain. For the case of an odd parity perturbation, $`\lambda _z`$ in equation (8) should be set equal to 4$`z_{\mathrm{node}}`$, and the simulations have to be performed with both hemispheres in the computational domain (see Table 1).

For completeness, we have also performed simulations with either outflow or inflow conditions at the $`z`$-boundaries. The results for the fastest growing modes are similar to ones presented here but, given that the magnetic and velocity fields change when the gas is allowed through the $`z`$-boundaries, some differences appear for perturbations with longer wavelengths. These results will be presented in a future communication.

### 2.5 Numerical Codes

The evolution of the undular instability is followed with both the isothermal MHD-TVD code described by Kim et al. (1999), and ZEUS-3D which has been described by Stone & Norman (1992a, 1992b). The first code is based on a second-order accurate, conservative, explicit, total variation diminishing (TVD) method. The $`𝐁=0`$ condition is maintained during the simulations by a scheme similar to a constrained transport (CT) scheme (Evans & Hawley 1988), which is reported in Ryu et al. (1998). The ZEUS-3D code, on the other hand, solves the ideal MHD equations on an Eulerian staggered mesh using an operator-split formalism wherein accelerations and heating are computed using centered differences followed by a conservative step. Monotonic upwind interpolation of either second or third order (van Leer or PPM advection) is used. Shocks are smeared using artificial viscosity. The CT algorithm is also adopted for the evolution of the magnetic field.

The main differences between the two codes are: i) the TVD code uses an approximate Riemann solver whereas ZEUS-3D uses a finite difference algorithm, ii) the present version of the TVD code solves the isothermal MHD equations (without the energy term) whereas ZEUS-3D has the adiabatic MHD with $`\gamma =1.01`$, and iii) ZEUS-3D allows less restrictive boundary conditions. The use of two different numerical algorithms, together with the linear analysis in Paper I, allows for a verification of the details of the numerical results. As described in the next section, and except for secondary details, the results from both codes are nearly identical. The advantage we have found in the isothermal TVD code is that it uses less cpu time per run. This is an expected advantage, however, because the code was optimized for isothermal problems and does not solve the energy equation. The important result of this exercise, then, is that both codes handle the problem in an adequate manner and with similar accuracies.

### 2.6 The Set of Numerical Experiments

In Table 1 we summarize the parameters of the numerical experiments presented in this paper. The simulations can be grouped into three categories. The first group is a set of $`(1+\frac{1}{2})`$D “equilibrium tests” (see below), performed with different resolutions and labeled as cases 1, 2, and 3. They were designed to check whether the codes were able to maintain the initial magnetohydrostatic equilibrium over long times. Within the accuracy of the algorithm for solving the equality in eqn. (1), a small force unbalance is allowed at every time step and its accumulative effect may affect the long-term evolution of the instability.

The second group of experiments is designed to compare the evolution, growth rates, and final equilibrium states for perturbations with different parities. This group is labeled as cases 4, 5, and 6. For this purpose, we have perturbed the initial state with MS and MA perturbations that have selected values of $`\lambda _y`$ and $`\lambda _z`$. The cases shown here have $`\lambda _y`$ equal to the most unstable wavelength and, as stated above, $`\lambda _z`$ is equal to an integer number of times $`z_{\mathrm{node}}`$.

The third group is illustrated with only one simulation, case 7, and is devised to pick up the value of the most unstable wavelength and the preferred parity (MS or MA) dominating the system. For this purpose, we generate random velocity perturbations in an extended computational domain.

## 3 RESULTS

As stated before, the coordinates ($`y`$, $`z`$) represent distances along and perpendicular to the midplane, respectively. The $`y`$-axis is anchored at the solar galactocentric radius and is quasi-azimuthal, defined by the tangent of the local field orientation in our Galaxy (see Heiles 1996 and Valleé 1997)

### 3.1 Equilibrium Tests: Thick Disk Oscillations

By $`(1+\frac{1}{2})`$D we mean that: i) the state vector $`𝐪`$ in eqn. (1) is a function of only $`z`$ and $`t`$, so $`𝐅_𝐲/y=0`$, but ii) the $`y`$-components of $`𝐯`$ and $`𝐁`$ survive and vary with $`z`$. In the isothermal TVD code, the source term of the gravitational acceleration is separated from the rest of the MHD equations by use of a splitting procedure (a detailed description of this procedure is given in Kim et al. 1999), but in ZEUS-3D the term is incorporated into the MHD equations (Stone & Norman 1992a,b). Here we verify if these different methods are able to maintain the magnetohydrostatic equilibrium over long periods of time.

After setting the initial condition described in §2.2, we trace the maximum vertical velocity as a function of time (Fig. 1). Numerical noise is produced by slight force imbalances at every cell, with wavelengths comparable to the cell-size, and propagates along the computational domain as compressive magnetosonic waves. The waves propagate through the stratified medium, and reach their peak velocities near the $`z`$-boundaries. The response of the gaseous disk to the continuous perturbation generated by this numerical noise is to oscillate with the natural frequency of the system. Obviously, the velocities involved in the normal mode of oscillation vary with height, from the ”shallow water wave” behavior near the plane to the more complex response at higher values of $`z`$ (see Martos 1993 and Martos & Cox 1998). The harmonic motion period in the linear gravity regime (for $`z`$ within 170 pc from midplane) is proportional to

$$\left(\frac{g(z)}{z}\right)^{1/2}|_{z=0}1.4\times 10^7\mathrm{yr},$$
(9)
and such approximation breaks down at high latitudes. Thus, the system acts as a complex set of coupled oscillators, and a modulated pattern should appear in our thick disk. This is indeed seen in the equilibrium tests performed with both codes.

In case 1, the simulation with the coarsest resolution, the peak velocity value shows modulated oscillations with a very small amplitude during the whole evolution. The resulting amplitude decreases to even lower values as the resolution is increased in cases 2 and 3. The top panel of Figure 1 shows the results for case 1 performed with ZEUS-3D. The oscillation amplitude is always less than 1% of the isothermal sound speed, and the plot shows the modulated oscillatory pattern. The oscillation is complex but three components with periods of about $`2.8\times 10^7`$, $`6\times 10^7`$ and $`11\times 10^7`$ year are revealed.

The lower panel of Figure 1 shows the same case performed with the isothermal TVD code. Since the TVD code is based on the Riemann solution, the gravity term should be treated separately. The code first updates the fluid variables (density, magnetic field, and velocity) without taking into account the gravitational acceleration. The force unbalance generated by this first step is then compensated by adding the vertical momentum due to gravity. There are residuals after this second step, because the code calculates the difference between two large numbers. Thus, due to a transient generated by the initial residuals, the first few velocity peaks are slightly larger than those obtained with ZEUS-3D, but the oscillation amplitude decreases below 1 % of the sound speed for the rest of the evolution. Also, the plot shows a cleaner oscillation mode, with a period of about $`6\times 10^7`$ year, because a fixed velocity pattern was enforced at the inner parts of the disk (i. e. this acts like a filter; the central parts of the disk are not allowed to oscillate at their natural frequencies, and one can isolate the oscillations of the extended layers).

The main conclusion of these tests, then, is that both ZEUS and TVD preserve the equilibrium state within an accuracy of a maximum spurious velocity of less than 1% after an elapsed time of about $`10^9`$ yr.

### 3.2 Sinusoidal MS and MA Perturbations

As stated above, cases 4, 5 and 6 have been constructed to follow the evolution of the undular instability triggered by MS and MA perturbations. From Paper I, the wavelength with the maximum growth rate in our thick disk model is about 18$`H_{\mathrm{eff}}`$, and this value is not sensitive to the parity of the perturbation. The three cases shown here have been done with this wavelength value but varying, aside from the perturbation parity, the domain size and grid resolution (see Table 1).

The four panels of Figure 2 show the rms value of the horizontal (dotted lines) and vertical (dashed lines) velocities as a function of time for both MS (left panels) and MA (right panels) perturbations. The solid line has a slope equal to the growth rate for the fastest growing mode with $`z_{\mathrm{node}}=9`$ derived in Paper I. The results obtained with the TVD code are displayed in the two upper panels, and the ones from ZEUS-3D are given in the lower panels. As mentioned before, the rms velocities are given in units of the isothermal sound speed (the panels have a natural logarithmic scale), and the units of time are $`H_{\mathrm{eff}}/a`$. The plots show a good consistency in the results attained with both codes.

The shapes of the velocity curves suggest that the evolution of the undular instability can be divided into three well defined stages: a linear growth, a nonlinear saturated regime, and a damping oscillatory stage. During the linear part (up to $`t=35`$ for the MS mode, and up to $`t=27`$ for the MA mode), the rms velocities increase at a constant growth rate. At early times, before $`t=20`$ for the MS mode and $`t=15`$ for the MA mode, the vertical rms velocity simply oscillates. The horizontal rms velocity, however, grows nearly steadily from zero (after a few oscillations) to its maximum value. A close look at the slopes for the horizontal rms velocities indicates that the slope of the MA case is slightly steeper than that of the MS case. This is a manifestation that the MA mode is slightly more unstable than the MS mode (Horiuchi et al. 1988). We did not derive this in the linear analysis of Paper I, but here show numerically that it is indeed the case in the next subsection. The end time for this phase can be defined as the moment when, due to the action of magnetic tension, the growth rate decreases abruptly to zero ($`t=35`$ for MS, and $`t=27`$ for MA). In a more strict sense, the term “linear” should be applied to times when the velocity has only a modest increase, smaller than the sound speed. Nonetheless, the rms values for the horizontal velocity at the end of the linear stage are nearly equal to the sound speed. This implies that the “linear approximation” in this case goes well beyond its conventional limits.

The end of the linear phase marks the beginning of the nonlinear saturated stage, in which the rms velocities reach their peak values. These peak values are maintained for a brief period of time, and then the velocities drop down and begin to oscillate about a nearly constant final value. It is a matter of taste to define the end epoch of this nonlinear stage. For simplicity, we take this moment as the time when the first oscillation begins, $`t=55`$ for the MS case and $`t=45`$ for the MA case. This corresponds to, approximately, the moment when the gas accumulated in the “magnetic valleys”, which is compressed by the ram pressure of the gas falling along the distorted field lines, begins to re-expand. After a large fraction of the halo gas has fallen into the valleys, the ram pressure decreases rapidly and the condensations readjust. The system also enters into its final equilibrium stage: as the rate of mass falling into the valley stalls, the weight of the condensations is stabilized and the tension of the distorted field lines is able to support them. Thus, the gas structures begin to stabilize in the damping oscillatory phase. As time goes by, the rms velocity gradually decreases due to numerical dissipation, and this is the reason why we use the term “damping” (see Kim et al. 1999).

Figure 3 shows snapshots of the density structures, velocity fields, and magnetic field topologies at three selected times. The snapshots obtained with both codes are nearly identical and, for simplicity, we show only the ones obtained with the TVD simulations. The left panels show case 4 with the MS perturbation, and the right panels show case 5 with the MA mode. The velocity fields are represented with arrows, and the standard velocity vector (with size $`2a`$) is shown at the top of the left panels. The $`B`$-field lines are chosen in such a way that the flux between consecutive lines is constant. The densities are color-coded from red to purple, as the value decreases. The initial sinusoidal perturbations propagate vertically in the form of magnetosonic waves and, as stated above, reach their peak values near the $`z`$-boundaries of the computational domain (see also Fig. 1).

The two upper panels, corresponding to the end time of the linear stage ($`t=35`$ for the MS case, and $`t=27`$ for the MA case), show that the instability is already fully developed. The midplane is slightly bended in the MA case. Later on, during the nonlinear stage, gas with supersonic velocities is continuously accumulated into the magnetic valleys (some weak shocks are generated during this stage), forming spur-like structures (see middle panels). The ram pressure of the falling gas keeps the layer thin and compressed, creating relatively large central densities. At the end of the nonlinear stage, the ram pressure drops and the over-compressed structures become thicker. The gas velocity near the edge of the spurs reverses its direction, and the system enters into the damping oscillatory state. The structures eventually move into a quasistatic state (lower panels), which is similar to the final equilibrium state obtained in Paper I.

The runs of the different energies in the system are plotted in Figure 4. The left side corresponds to the MS case, and the right side to the MA case. The internal, kinetic, magnetic, gravitational and heat-exchange energies are defined by

$$E_i=\frac{3}{2}p𝑑y𝑑z,E_k=\frac{1}{2}\rho (v_y^2+v_z^2)𝑑y𝑑z,$$
(10)
$$E_m=\frac{1}{8\pi }(B_y^2+B_z^2)𝑑y𝑑z,E_g=\rho \varphi 𝑑y𝑑z,E_h=p\mathrm{ln}pdydz,$$
(11)
where $`\varphi =_0^zg(z)`$ is the gravitational potential, and the double integrals are performed over the whole computational domain. Due to the isothermal condition, the internal energy is just proportional to the total system mass. Since there is no mass exchange between our system and the outside, the total mass is conserved in our simulations (see solid lines, labeled $`E_i`$). We use the internal energy as the normalization unit for all energies.

Following Mouschovias (1974), we define the energy integral (or total energy) as

$$E_t=E_k+E_m+E_g+E_h.$$
(12)
When shocks appear, as in our simulations, this energy is not conserved during the evolution (see also the simulations of Matsumoto et al. 1990). The curve labeled $`E_t`$ in Figure 4 is constant up to nearly the end of the linear stage, because the perturbation speeds are below $`a`$ at those stages. When shocks begin to form during the nonlinear regime, the total energy begins to decrease but it levels up later, when the gas speeds fall below the sound speed and the system enters into the damping oscillatory stage. This is particularly clear in the MS case. For the MA case, there are very weak shocks even at the damping oscillatory stage and the total energy curve is still decreasing at $`t=80`$ (see lower right panel in Fig. 3).

The kinetic energy, $`E_k`$, in both MS and MA cases shows small variations, even though shocks are generated during the nonlinear stage. This is due to the fact that these shocks occur in low density regions. The heat-exchange term, $`E_h`$, corresponds to the amount of work done by or on the gas (when $`\mathrm{\Delta }E_h`$ is positive, work is done by the gas; see Mouschovias 1974). The gravitational energy, $`E_g`$, decreases as gas from high latitudes falls down to the magnetic valleys. Comparing the runs for $`E_g`$ in both panels, one finds that a larger amount of gravitational energy is released in the MA case. This is because the MA oscillations are able to bend the midplane, and a larger amount of gas lowers its position to be collected in the high density condensations (see Paper I). The gravitational energy is first transformed into kinetic energy, and then it is either dissipated in shocks or goes into magnetic field compression and tension. In general, the instability proceeds at the expense of lowering the gravitational energy (Mouschovias 1974).

To compare the final equilibrium state derived in Paper I with the late times results of the numerical simulations, we obtain the column densities as a function of $`y`$. For a given horizontal location, the column density at any time $`t`$ is

$$N(y,t)=_{z(y,A[z=0,t=0])}^{z_{\mathrm{max}}}\rho (y,z,t)𝑑z,$$
(13)
where the upper limit of the integral corresponds to the upper boundary of the computational domain, and the lower one is the $`z`$-coordinate of the magnetic field line that was initially located at midplane (labeled $`A[z=0,t=0]`$; see Paper I). During the linear stage, the run of the column densities has only small variations. The important changes occur during the nonlinear stage, when the contrast between condensations and voids becomes maximum. Later on, during the oscillatory stage, the column densities maintain the values achieved during the end of the nonlinear regime and the values at $`t=80`$, normalized to the initial column density value $`N(t=0)`$, are plotted in Figure 5 for both MS (dotted line) and MA (dashed line) perturbations. The MA case has the larger density contrast between condensations and voids. As stated before, this is because the MA perturbations provide more room for the high density gas by bending the midplane. At this last time, the undular instability has almost reached the new final equilibrium in both MS and MA cases, and the final column density configurations are very similar to those shown in Figure 5 of Paper I. The differences in the numerical and analytical results for the MS case are of only about 10%. The differences in the MA case are in general larger, particularly in the voids, but they are also of the order of 10 % at the center of the magnetic valleys.

In Paper I, we found that the growth rate increases, up to a certain maximum value, with the location of the nodal point. This is because gravity, the driving force of the instability, increases up to a certain maximum value. To verify this dependence in the numerical experiments, we made a series of cases with different values for the $`z`$-boundaries. Here we show one simulation, case 6, with an extended vertical domain, $`12z12`$, and a sinusoidal MA perturbation. Figure 6 shows the resulting rms velocities (dotted and dashed lines as in Fig. 2) as a function of time. The solid line now has a slope equal to the growth rate corresponding to the nodal points for this case. As is clear from the figure, the growth of the linear part follows the rate predicted by the linear analysis for this case.

### 3.3 Random Perturbations

As a final issue of this part, we now verify the wavelength value and parity of the fastest growing mode in the thick disk system. For this purpose we have performed simulations with random velocity perturbations, without any preferred parity or wavelength. The horizontal size of the computational domain has been extended, and the size of the simulation in case 7 is three times longer than in the previous models. We also have increased accordingly (by a factor of three) the number of cells in the $`y`$-axis.

Figure 7 shows the run of the rms velocities for this case and, as before, the slope of the solid line corresponds to the growth rate. The horizontal component (dotted line) has an initial adjustment and enters into the linear growth at $`t8`$. The velocity grows linearly, up to $`t45`$, with the predicted maximum growth rate. As in the previous models, the rms of the vertical component has some initial fluctuations and then, after $`t32`$, it also grows with the predicted rate. The rest of the evolution is similar to the one already described in case 5. The only difference is that the evolutionary stages are now significantly delayed. The resulting structures are displayed in two snapshots, at $`t=45`$ and $`t=63`$ (the end of the linear and nonlinear stages, respectively), in Figure 8. The series of condensations are formed at slightly different times, but they have the inter-distance scale of the most unstable wavelength and the MA symmetry. This shows that, even when there was a parity degeneracy in the linear analysis, the odd parity mode is dominant.

## 4 DISCUSSION

In this paper, we have confirmed numerically the results of the linear analysis obtained in Paper I for the undular mode of the Parker instability in a thick, warm disk model of the Galaxy. The dependence of the growth rate on the position of the upper boundary of the numerical grid has been confirmed, and an explanation in terms of the action of the gravitational field is provided. The main properties of the final equilibrium configuration were also confirmed. The differences between the final states for the symmetric and antisymmetric modes were corroborated by the numerical study. The actual values of the column density distributions are slightly different, but the shapes and the final-to-initial column density ratios were confirmed. In addition, the present study also provides a demonstration that the antisymmetric mode is dominant and its signature should be present – that of an odd parity with respect to the central galactic plane – in structures formed by the instability. In addition to the above, we are confident of the new evolutionary details of the instability in a multi-component thick disk, since we have got the same results with two different numerical codes (TVD and ZEUS-3D).

These new details may have important implications in our understanding of the galactic disk. The extended vertical structure observed in our Galaxy and others, as the diffuse ionized warm layer, forces higher scale heights (for mass, pressure and $`B`$-fields) than assumed in previous calculations of the Parker instability. The resulting thick disk alters considerably the time and length scales involved in the process and, as long as the upper layers are not dominated by a hot coronal plasma phase, the disk is unstable to undular modes. Evidence for the pervasiveness of this hot phase is still under debate (see Cox 1995, and the conference proceedings edited by Arthur, Brickhouse & Franco 2000) but, regardless of the actual details of the gas temperature distribution, one can explore the consequences of the instability with the present isothermal thick disk model.

In considering some of the possible applications of this work, one has to keep in mind the fact that our results refer only to a 2–D analysis. The final equilibrium states, with a well defined periodicity of troughs and crests, are unstable to 3–D perturbations (Asseo et al. 1978). Nonetheless, tridimensional simulations in a thin disk (Kim et al. 1998) show that some of the 2–D structuring is still present in this case. The instability in this case is a combination of two modes, interchange and undular. The interchange mode corresponds to the case in which the field lines maintain their initial direction but some flux tubes are displaced upwards when they find themselves less dense than the surroundings. The undular mode (which is the one addressed in this study, and is simply referred to as the Parker instability) corresponds to the case in which the field lines are wavy above the midplane and gas aggregation occurs at the magnetic valleys. Three dimensional motions are more easily destabilized than those in two dimensions, for density fluctuations (which instigate the interchange mode) are related to the total (gas plus magnetic) pressure, whereas the undular mode depends on the gas pressure alone (Hughes & Cattaneo 1987).

The most unstable growth rate of the interchange mode has a vanishing wavelength along the third dimension, which is perpendicular to both gravity and initial field directions. This fact has been known from one of Parker’s (1967) pioneering works, and can limit the size and mass of the resulting structures. Our isothermal disk is unstable to the 3–D interchange mode (Martos 1993), and its fastest growth rate is that of an infinitesimally small wavelength perpendicular to the 2–D plane of motion considered in our calculations. As a result, and depending on the coupling of the field lines, small scale structuring or turbulence will be also made. Additional physical effects, such as rotation, self–gravity, or variable gravity in the $`z`$ direction, cannot stabilize the buckling of the field lines (Asseo et al. 1978; Elmegreen 1991). Thus, the instability will tend to proceed faster in 3–D than in 2–D, creating complex small scale structuring, but some of the 2–D signatures will be also present at the large, kiloparsec size scales (i. e. the large scale undular mode, with its odd parity with respect to the central galactic plane, will persist). The relative importance of the undular large scale structures depends on whether one is considering arm or inter-arm regions. The inter-arm regions are sheared, and the effects of rotation in the interchange mode lead to a preferential small scale structuring (see Kim, Ryu & Jones 2000). Thus, little or no undular signatures are found, and alternating dense and rarefied sheets are formed by the three-dimensional perturbations (Kim et al. 1998, 2000). The density structure of these sheets are far from those expected for giant molecular clouds, and the enhancement factors for the vertical column density are small. Thus, it is difficult to consider the Parker instability alone as the formation mechanism of giant clouds in the general interstellar medium (the resulting 3–D structures may be better more reminiscent of the worms observed in atomic hydrogen; see Koo, Heiles & Reach 1992).

In contrast, inside the arms, with a stronger gravitational field and no shear, the large scale undular structures are clearly present along the arms (Franco et al. 2000). The spacing between the resulting condensations is similar to the most unstable mode derived here, and the antisymettric mode creates a corrugation pattern along the spiral arm. An interesting test of this prediction is suggested by the structuring detected in the Carina arm, with a corrugation scale of 2.4 kpc (Alfaro, Cabrera-Caño & Delgado 1992). Early 3–D studies already showed that the spiral density wave triggers the instability in the thick gaseous disk (Martos & Cox 1994), and current work in 3–D (Franco et al. 2000) indicates that the undular mode indeed gives rise to a corrugation, or bending, along the arm. There is good agreement with the observed corrugation length when appropriate scaling of the gravitational field towards the inner Galaxy is accounted for. This feature, along with the hydraulic jumps generated by the action of spiral density waves (Martos & Cox 1998), can create strong spatial and velocity distortions in the vicinity of spiral arms. These issues certainly require further observational studies of corrugations and the velocity fields across spiral arms (Alfaro et al. 2000).

Physics not included in this study will alter some details of the picture discussed above. For instance, the Galactic field is mostly tangled and the random component must have a stabilizing influence (Zweibel & Kulsrud 1975). Also, self–gravity (Elmegreen 1982) and twisting of the field lines (Shibata & Matsumoto 1991; Kim, Ryu & Jones 2000) will contribute to the formation of coherent, finite (but small) size structures, setting more realistic lower limits to the theoretical zero wavelength predicted by theory for the fastest growing interchange mode. In addition, magnetic reconnection should have played an important role at the late evolutionary stages. The twisted random field components may lead to a rapid energy release (e. g. Lazarian & Vishniac 1999), and can create additional flows and modify the final energy budget of the process. We notice, however, that the scales of the problem addressed here are much larger than the cell size of the random component, 50 to 100 pc. The effect of galactic differential rotation, as stated above, will modify the outcome of the instability (e. g. Foglizzo & Tagger 1994, 1995). Finally, the timescales for our 2–D simulations are much larger than the cosmic ray residence time of 2.$`\times 10^7`$ years, making a causal association between cosmic ray escape and the Parker instability unlikely. But again, this, as most previous issues, demands precise three–dimensional studies, such as those presented by Kim et al. (1998) and Kim, Ryu & Jones (2000), and should be addressed in the near future.

We warmly thank Gene Parker and Don Cox for useful comments and suggestions during the development of this project, and Jane Arthur for a careful reading of an early version of this paper. The comments and suggestions made by the referee are gratefully acknowledged. We also thank M. Norman, M. MacLow and R. Fielder for continued advice on ZEUS-3D. This work has been partially supported by a bilateral CONACYT-Mexico and KOSEF-Korea agreement. JF, MM and AS acknowledge partial support from DGAPA–UNAM grant IN130698, and by a R&D grant from Cray Research Inc. The work by JK was supported in part by the Office of the Prime Minister through Korea Astronomy Observatory grant 99-1-200-00. The work by SSH was supported in part by grant BSRI-98-5411. The work by DR was supported in part by the grant KRF-99-015-DI0113. The numerical calculations were performed using UNAM’s ORIGIN-2000 supercomputer and SUN Enterprise 3500 at the Korea Astronomy Observatory.

FIGURE CAPTIONS

Fig. 1.- The panels show the maximum velocity of case 1 as a function of time. The upper plot shows the results obtained with the ZEUS code, and the lower plot shows the same case performed with the TVD code. The velocities and times are in units of $`a`$ and $`H_{\mathrm{eff}}/a`$, respectively.

Fig. 2.- The frames show the run of the rms velocities for the undular instability. The panels display the rms value of the horizontal (dotted lines) and vertical (dashed lines) velocities, as a function of time, for both MS (left panels) and MA (right panels) perturbations. The upper frames were obtained with the TVD code, and the lower ones with ZEUS-3D. The solid lines have a slope equal to the growth rate for the fastest growing mode with $`z_{\mathrm{node}}=9`$ derived in Paper I.

Fig. 3.- The Parker instability in 2-D for a magnetized multi–component gaseous disk, performed with the TVD code. The sequence shows the density (color logarithmic scale), velocity field (arrows) and magnetic fields (lines), at three selected times: $`t`$ = 35, 45, 80 for the MS mode (case4), and $`t`$ = 27, 36, 80 for the MA mode (case5).

Fig. 4.- The time evolution of different energies (internal $`E_i`$; kinetic $`E_k`$; magnetic $`E_m`$; gravitational $`E_g`$; heat-exchange $`E_h`$, and total $`E_t`$) for the simulations presented in Fig. 3. All energies are normalized to their initial values. The left panel corresponds to the MS case, and the right side to the MA case.

Fig. 5.- Final column density values as a function of horizontal location, for the MS (dotted line) and MA (dashed line) cases of Fig. 3. The column densities are normalized to the initial value, $`N(t=0)`$.

Fig. 6.- Same as in Fig. 2, but for case 6 performed with the TVD code. The slope of the solid line represents the growth rate corresponding to the fastest growing mode with nodal point $`z_{\mathrm{node}}=12`$ obtained in Paper I.

Fig. 7.- Same as in Fig. 2, but for case 7 performed with the TVD code. The slope of the solid line is the same as in Fig. 2.

Fig. 8.- Parker instability in a magnetized multi–component gaseous disk for case 7. The sequence shows the density (color logarithmic scale), velocity field (arrows) and magnetic fields (lines), at two selected times: $`t`$ = 45, 63.


text:
# Strong Discontinuities in the Complex Photonic Band Structure of Transmission Metallic Gratings

## Abstract

Complex photonic band structures (CPBS) of transmission metallic gratings with rectangular slits are shown to exhibit strong discontinuities that are not evidenced in the usual energetic band structures. These discontinuities are located on Wood’s anomalies and reveal unambiguously two different types of resonances, which are identified as horizontal and vertical surface-plasmon resonances. Spectral position and width of peaks in the transmission spectrum can be directly extracted from CPBS for both kinds of resonances.

Recently, metallic films with sub-wavelength apertures became the subject of increasing interest. Experiments show that 2D arrays of sub-wavelength holes in metallic films can lead to extraordinary transmission of light . These properties have been attributed to the excitation of coupled surface-plasmons on the upper and lower surfaces of the grating. The ability of these structures to control light has been shown , and applications have already been proposed in order to exploit these properties in different fields, such as electromagnetic filters or photolithography.

Electromagnetic calculations have been carried out by Porto *et al.*  in order to study the mechanisms that enhance the transmission of light through metallic gratings with very narrow slits. They distinguished two different mechanisms, that is the excitation of coupled surface plasmon polaritons on both surfaces of the metallic grating, and the coupling of incident plane waves with waveguide resonances located in the slits. They correspond respectively to surface-plasmon bands and flat bands in the energetic band structure. Similar resonances located in grooves have also been observed experimentally and numerically in the reflected light of metallic gratings .

Complex photonic band structures (CPBS) were previously calculated by Kuzmiak *et al.* for 2D periodic systems with metallic components , in order to obtain the attenuation and the lifetime of each mode. CPBS allowed a splitting of the lifetime of degenerate modes at Brillouin-zone boundaries to be observed.

In this Letter, we will show that, in the case of rectangular metallic gratings, complex dispersion curves demonstrate the existence of a new kind of discontinuities. Much stronger than for 2D periodic systems, they are located on Wood’s anomalies. They correspond to the transition between two different types of resonances, identified as horizontal and vertical surface-plasmon resonances. The first one is a periodic structure resonance, the second one is a Fabry-Perot like resonance. The discontinuity can be up to four orders of magnitude. Moreover, we will see on zero-order transmission spectra that the calculated complex frequencies are in excellent accordance not only with the spectral position, but also with the width of resonances’ peaks.

Fig. 1 shows the structure studied in this Letter, which is similar to the one of Porto *et al.* . In the following, the period of the grating $`d=3.5\mu \mathrm{m}`$ and the width of the slits $`a=0.5\mu \mathrm{m}`$ will be kept constant, whereas the height $`h`$ of the grating will vary. It is a symmetric structure, surrounded with air. The metal is assumed to be gold, whose dielectric function is taken from the reference used by Porto .

Our calculations were carried out for the TM polarization (the incident magnetic field is parallel to the direction $`y`$ of the grating, see Fig. 1). They are based on the exact modal expansion in region II  and $`S`$-matrices formalism. In the upper and lower regions (I and III), the electromagnetic field is expressed as superposition of plane waves $`e^{[i(k_x^{(n)}x\pm k_z^{(n)}z)i\omega t]}`$ (Rayleigh’s expansion). Their pure real wave vector component $`k_x^{(n)}`$, which represents the in-plane momentum of photons is given by Eq. (1):

$$k_x^{(n)}=k_x^{(0)}+\frac{2\pi n}{d}.$$
(1)
The two other terms $`k_z^{(n)}`$ and $`\omega `$ are complex quantities, according to the relation $`\left(k_x^{(n)}\right)^2+\left(k_z^{(n)}\right)^2=\left(\omega /c\right)^2`$.

For the calculation of CPBS, we use a $`S`$-matrix formalism, in which the amplitude of the reflected and transmitted diffracted waves are given by a matrix $`S`$ applied to a vector containing the amplitude of the upper and lower incident waves. An electromagnetic mode is defined by the condition of the existence of a wave in the structure without any incident wave, which is equivalent to say that the inverse ($`S^1`$) of the $`S`$-matrix has an eigenvalue that is null, with an eigenvector $`(R_n,T_n)`$, corresponding to an electromagnetic field in regions I and III:

$$\mathrm{\Psi }^I=\underset{n}{}R_ne^{[i(k_x^{(n)}x+k_z^{(n)}z)i\omega _0t]},\mathrm{\Psi }^{III}=\underset{n}{}T_ne^{[i(k_x^{(n)}xk_z^{(n)}z)i\omega _0t]},$$
(2)
with $`\mathrm{}(ik_z^{(n)})0`$.

For a given $`k_x=k_x^{(0)}`$, the electromagnetic modes are defined by complex frequencies $`\omega _0`$. The inverse of the imaginary part of the frequency provides the amplitude lifetime $`2\tau `$ of the resonance , which can be associated to the dimensionless quality factor $`Q=\mathrm{}(\omega _0)\tau `$.

The complex dispersion curves are lines in the three-dimensional space defined by the complex frequency $`\omega _0`$ of resonance modes and the parallel momentum of photons $`k_x`$. For a more comprehensive representation we are plotting two projections of $`\omega _0`$, the energy $`E=\mathrm{}\mathrm{}(\omega _0)`$ and the dimensionless quality factor $`Q=\mathrm{}(\omega _0)/\left(2\mathrm{}(\omega _0)\right)`$.

Fig. 2 shows the complex dispersion curves calculated for a height $`h=1.4\mu \mathrm{m}`$. The lines on the energetic band structure (Fig. 2(a)) represents the frontier of the Wood-Rayleigh’s anomalies. They correspond to the emergence (or the vanishing) of a diffracted wave in the region I or III, which wave is then at the grazing angle. These anomalies are observed on transmission spectra when the $`n^{\mathrm{th}}`$ parallel wave vector $`k_x^{(n)}`$ (given by Eq. (1)) is equal to the real incident wave vector $`k_0`$.

The energetic band structure in Fig. 2(a) shows two different types of modes, whose dispersion relations have already been studied by Porto. The bands located close to Wood’s anomalies, with slightly smaller energies, characterize horizontal surface-plasmon resonances. They are excited by the first non-propagating (or evanescent) diffracted wave, and their quality factors are greater than 10 and up to more than $`10^4`$ (see Fig. 2(b)). On the other hand, flat bands represent electromagnetic modes which are almost independent of the incident angle $`\theta `$, and with lower quality factors (below 20). They were previously reported as waveguide resonances . We will justify further in this Letter why these resonances are more precisely vertical surface-plasmon resonances.

The most remarkable result is given by the quality factor of the upper energetic band, which is broken into three parts. Considering the point A in Fig. 2(a), we see that the energetic band is continuous whereas the quality factor’s curve is discontinuous at the corresponding points A’ and A” in Fig. 2(b). On the left side of point A, the flat band corresponds to a vertical resonance with very low quality factors (of the order of 3), whereas the right side corresponds to an horizontal surface-plasmon resonance with high quality factors (up to more than $`10^4`$). The lifetimes of the resonances are $`\tau =2.9\mathrm{fs}`$ at A’ and $`\tau =33\mathrm{ps}`$ at A”.

In order to understand physically the behavior of the electromagnetic waves, it must be pointed out that the discontinuous transitions between the two types of resonance are located exactly on Wood’s anomalies (point A and its symmetrical point B), so that they are linked to the emergence or the vanishing of a diffracted wave. This phenomenon can be explained in considering the energetic band structure: below Wood’s anomaly, the first non-propagating (evanescent in $`z`$ direction) diffracted wave ($`n`$ = +1 at the point A and $`n`$ = -2 at the point B) is responsible for the excitation of the horizontal surface-plasmon resonance. When the resonance crosses Wood’s anomaly, the evanescent wave becomes a $`z`$-direction propagating diffracted wave (at grazing angle) which can not excite an horizontal resonance anymore, but is coupled to vertical surface-plasmons located on the vertical walls of slits. Then, we obtain a flat band which corresponds to a ”weaker” resonance in the Fabry-Perot like cavity delimited by the top and bottom surfaces of the grating. This vertical resonance is coupled with the last propagating diffracted wave, thus it ends at the crossing with Wood’s anomaly, where this wave vanishes.

The complex band structure in Fig. 2 has also a continuous transition between the horizontal and the vertical resonances in the lower band. The flat band at normal incidence corresponds to vertical resonances, and becomes a surface-plasmon band for greater incident angles. Its quality factor increases from 20 to about 300 during this smooth transition. The last flat band at about 0.35 eV corresponds to a vertical resonance for the whole range of $`k_x`$ in Fig. 2.

By varying the width $`a`$ of the slits, we have observed that the complex mode of the vertical resonance is nearly unchanged. In the calculation, it can be seen from the modal decomposition of the electromagnetic field, that this resonance corresponds mainly to the excitation of the first eigenmode of a $`z`$ infinitely extended region II. The wave is evanescent in the $`x`$ direction and propagates at the vertical interfaces air/metal ($`z`$ direction). This TM mode exists for any width of slit, and becomes equal to the wave vector $`k_{ps}`$ of the surface-plasmon of a semi-infinite air-metal interface in the limit of wide slits, so that it is more correct to call it a *vertical surface-plasmon mode* instead of a waveguide mode. The transitions between the regions II/I and II/III act as mirrors for the vertical surface-plasmon mode. Their reflection coefficients are low, so that the vertical resonance is a Fabry-Perot like resonance with low quality factors.

Complex dispersion curves shown above can be seen as convenient compact representations of the resonances, which can be associated to transmission spectra. For that purpose, we make the hypothesis that the complex frequency $`\omega _0^i`$ of a resonance $`i`$ is a single pole of the amplitude of the transmitted waves. Consequently, the zero-order transmitted energy can be approximated by a Lorentzian function near the resonance $`i`$. Assuming $`\omega `$ is the (real valued) frequency of the incident wave, we obtain:

$$T_0=\left|\frac{\alpha _i}{\omega \omega _0^i}\right|^2.$$
(3)
The complex band structure of the grating is shown in Fig. 3 for a height $`h=3\mu \mathrm{m}`$. Fig. 4 shows the zero-order transmission spectra calculated for the same height and two different angles of incidence. They correspond to $`\theta =0^{}`$ (a) and $`\theta =12^{}`$ (b) cuttings in Fig. 3(a). The Lorentzian curves corresponding to Eq. (3) are added to the spectra. The constant $`\alpha _i`$ is fitted to obtain the same maximum value as the spectra. At normal incidence (Fig. 4(a)), the broad resonances at 0.38 eV and 0.54 eV (not represented in Fig. 3) are not separated enough to allow the comparison to their respective Lorentzian functions to be made. In such a case, the amplitude of the transmitted wave is the sum of two complex Lorentzian functions, and phase considerations make it more difficult to fit to the spectra.

For the other peaks, Fig. 4 shows an excellent agreement of the Lorentzian curves with the spectra. Hence, the single-pole approximation is justified and Eq. (3) gives the width at half maximum:

$$\mathrm{\Delta }\omega _0=\mathrm{\hspace{0.25em}2}\mathrm{}(\omega _0^i).$$
(4)
It justifies also the above definition of the quality factor of the resonance, which corresponds to the more conventional $`Q=E/\mathrm{\Delta }E`$.

These transmission spectra allow some important characteristics of the resonances to be pointed out. At first, it must be noted that the flat band at about 0.35 eV has no solution at normal incidence, and the peak $`\gamma `$ disappears at $`\theta =0^{}`$. An energetic shift can be seen for the horizontal surface-plasmon resonances between $`\theta =0^{}`$ and $`\theta =12^{}`$, whereas the vertical resonance at about 0.17 eV does not move. Another point which comes out of the calculations is that vertical resonances in contiguous slits are nearly uncoupled. As a consequence, a periodic structure is not necessary for them to exist: the vertical surface-plasmon resonance could be excited in an unique slit, and fully coupled with a few ones. Hence, an important result is that the vertical resonance can be efficiently excited by a strongly focalized beam.

In conclusion, we have carried out the calculation of the complex dispersion curves of rectangular metallic gratings in TM polarization. CPBS are directly linked to experimental data, like transmission spectra. The quality factor and the lifetime of resonances, that they provide, are crucial characteristics in many applications already proposed for the extraordinary effects appearing in such transmission metallic gratings. CPBS show pronounced discontinuities located on Wood’s anomalies, which can not be evidenced in the usual energetic dispersion curves. They allow to distinguish unambiguously the two different natures of resonances: *horizontal and vertical surface-plasmon resonances*. Similar results can be obtained for reflection gratings. This work can also be extended in order to analyze resonances in more complex geometries, as Ebbesen’s structures  for example, for which a cutoff diameter exists for vertical modes. The characteristics of the resonances studied in this Letter open up new possibilities in photonic materials. Moreover, these structures allow huge concentrations of light in very small volumes, which could have applications in fundamental optics  and in optoelectronics.