format: text

text:
# 1 Introduction

## 1 Introduction

One of the challenges put forward by the Seiberg-Witten proposal  for an exact expression of the prepotential in $`N=2`$ extended $`d=4`$ Super Yang-Mills theory is its verification through microscopic instanton calculations. Dorey, Khoze and Mattis made the first steps into that direction determining the 1 and 2-instanton contribution to the prepotential . (The 1-instanton contribution was already calculated before in references  and .) These authors also derived a representation of the integration measure of instanton moduli in general . After this work further progress was hampered by the fact that beyond instanton number 3 the ADHM (Atiyah, Drinfeld, Hitchin, Manin ) moduli space is not known explicitly. In ,  we derived a standard algebraic-geometric form, to be reproduced below, representing the expansion coefficients of the prepotential as integral of the exponential of an equivariantly exact mixed differential form. This can be reduced in the $`k`$-instanton sector to an integral of a (4k-3)-fold wedge product of a closed differential 2-form, which is the formal representative of the Euler class of the moduli space viewed as U(1) principal bundle. So we were led to an algebraic-geometric interpretation of the coefficients. But we failed in our main objective, which was to calculate via localization technique concrete numbers. The reason for our failure was clarified and cured in a recent work by Hollowood , where it is argued that one should, to start with, resolve the short distance singularities of the moduli space111 R.F. thanks Francesco Fucito for a discussion of this point..

An elegant method to achieve the resolution consists in deforming the instanton configurations into a non-commutative domain . The parameter measuring the non-commutativity, a new length scale, acts as an ultraviolet regulator, while the integrals in question are - as formally shown below- insensitive to the deformation. The integrals can now be localized at the surface of critical points of the abelean vector field which goes with the remaining torus symmetry after spontaneous breakdown of the non-abelean gauge symmetry of the N=2 Yang-Mills theory222Localization without regularization renders a vanishing residuum at the corresponding critical surface..

The critical surface happens to coincide with the the manifold of non-commutative U(1) instantons (imbedded as submanifold into the ADHM moduli space of the non-abelean gauge group). The remaining integral along the critical surface has according to the recipes of the equivariant localization technique ,  the appearance of a sum of characteristic classes of the unitary bundle connected with the ADHM construction. Their evaluation is for general instanton number still a highly cumbersome task.

A conceivable way to overcome the difficulties relies on the deformation invariance of the characteristic classes. So one may follow backwards a gradient flow along a Gl(k) orbit  (chapter 6), , which relates the relevant matrices of the ADHM construction to diagonal matrices. We will not follow this approach here since we were reminded meanwhile by Nekrasov  of a technically simpler method333We still believe that it is worthwhile to pursue attempts either to evaluate the integrals without regularization or to handle the U(1)-instanton manifolds. It may lead to new mathematical insights.. The idea, already used some time ago in a macroscopical context , , is to modify the vector field underlying the localization method by adding pieces corresponding to space rotations. The critical set of the modified vector field consists of discrete points and the evaluation of the integrals is reduced to a manageable task: to find the eigenvalues of the vector field action on the tangent space of the respective critical point. The purpose of this paper is to show that the program can be implemented on the microscopic level following the conceptual lines of our previous papers , . Similar results have been announced in . We collect in the following section 2 some material concerning the instanton calculus. Section 3 is devoted to a description of the algorithm allowing to determine (in principle) the contributions from arbitrary instanton sectors. We will work out for the purpose of illustration the 2 ,3 and 4-instanton coefficient of the pure N=2 vector theory (i.e. without matter multipletts) with SU(2) as underlying gauge group in section 4. We end this section presenting a formula for the determinant of the deformed vector field action on the neighborhood of an arbitrary critical point in the general case of $`SU(N)`$ as gauge group.

## 2 Instanton calculus

We quote the ADHM data for the construction of $`SU(2)`$ instantons in a form which is generalizable to $`SU(N)`$, $`N2`$. (The restriction to $`SU(2)`$ is a matter of convenience; all our results can easily be extended to $`SU(N)`$). These data consist of complex matrices, two $`k\times k`$ matrices $`B_1`$, $`B_2`$, a $`k\times 2`$ matrix $`I`$ and a $`2\times k`$ matrix $`J`$, fulfilling a certain rank condition (cf. , ) and obeying the relations

$`[B_1,B_2]+IJ=0;`$ (2.1)
$`[B_1,B_1^{}]+[B_2,B_2^{}]+II^{}J^{}J=0,`$ (2.2)
where $`\mathrm{"}^{}\mathrm{"}`$ denotes hermitean conjugation. The data are redundant in the sense that sets of matrices related by $`U(k)`$ transformations,

$`𝒜(B_1,B_2,I,J)𝒜^{},`$
$`B_i^{}=gB_ig^1,I^{}=gI,J^{}=Jg^1;`$
$`i=1,2;gU(k)`$ (2.3)
give rise to the same self dual Yang-Mills configuration of topological charge $`k`$. Factoring out the redundancy from the above data (e.g. by imposing gauge fixing conditions) one is lead to the instanton moduli space as the set of gauge equivalence classes of solutions of Eq.’s (2.1), (2.2). This space is not smooth. It contains boundary points where the previously mentioned rank condition is violated (which corresponds to Yang-Mills configurations with topological charge density concentrated in part at single points). A resolution of singularities is achieved  through a modification of Eq. (2.2):

$$[B_1,B_1^{}]+[B_2,B_2^{}]+II^{}J^{}J=\zeta \mathbf{\hspace{0.17em}1}_{k\times k},$$
(2.4)
where $`\zeta `$ denotes a positive real number. Eq. (2.4) is in fact the starting point for the construction of solutions of the Yang-Mills selfduality equation on a non-commutative space-time  with $`\zeta `$ setting the scale of the non-commutativity. The appearance of the non-commutative space will be of no relevance for our purposes. Going from Eq. (2.2) to Eq. (2.4) we will content ourselves with an argument, to be given below, that the numbers we want to calculate do not depend on $`\zeta `$.

Dealing with a supersymmetric theory in quasi-classical approximation one has to take into account fermionic degrees of freedom which are here the Weyl zero modes of positive chirality in the selfdual Yang-Mills background. It turns out, , , that a neat realization of the fermions is supplied on the level of moduli through their identification with a basis of the cotangent bundle of these moduli. To adapt the cotangent space in a $`U(k)`$-invariant way we introduce some notation. Let $`𝐋`$ be a selfadjoint operator acting on the space of anti-hermitean $`k\times k`$ matrices $`M`$ by

$`𝐋M=\{II^{}+J^{}J,M\}+{\displaystyle \underset{l=1,2}{}}\left([B_l,[B_l^{},M]]+[B_l^{},[B_l,M]]\right).`$ (2.5)
It can be shown that $`𝐋`$ is invertible (in fact positive) as long as the above mentioned rank condition is satisfied (under the supposition of Eq.’s (2.1), (2.2). The invertibility is guaranteed per se if one relies on Eq.’s (2.1), (2.4) ). Let $`𝐗`$ be the matrix valued differential one-form

$`𝐗={\displaystyle \underset{i=1,2}{}}[B_i^{},dB_i]+J^{}dJdII^{}h.c.`$ (2.6)
We introduce a $`U(k)`$-covariantized exterior derivative on the ADHM matrices $`𝒜(B_1,B_2,I,J)`$ setting

$`𝒟𝒜d𝒜+Y𝒜,`$ (2.7)
with $`Y=𝐋^1𝐗`$ and $`Y𝒜([Y,B_1],[Y,B_2],YI,JY)`$. $`𝒟𝒜`$ satisfies by construction (cf. , ), the fermionic part of the supersymmetrized ADHM conditions (2.1) and (2.2) (or (2.4) resp.) thus giving rise to Weyl zero modes. It is worth to mention the geometric meaning of the $`U(k)`$ connection field $`Y`$. The flat euclidean space of data $`𝒜`$ (without Eq.’s (2.1), (2.2) ((2.4)) can be supplied with a Kähler metric. The Kähler property is preserved by the restriction imposed by Eq.’s (2.1), (2.2) ((2.4)) provided one projects in the tangent bundle of the moduli manifold onto the lifted (with respect to the $`U(k)`$ action) horizontal subspace , . To come to terms with the latter notion let us assume that the moduli space is realized through the imposition of a gauge fixing condition, e.g. by demanding that, say, the hermitean part of $`B_1`$ is diagonal. The tangent space of the larger manifold (without gauge fixing) may be decomposed into the subspace spanned by the infinitesimal $`U(k)`$ generators, the so called ”vertical” subspace and its orthogonal444orthogonal with respect to the induced metric. complement, the ”horizontal” subspace, the latter being isomorphic to the tangent space of the gauge fixed moduli manifold $`_k`$. Given a tangent vector $`T`$ on $`_k`$ one finds a unique vector $`V`$ in the vertical subspace s.t. $`(V+T)`$, called ”the horizontal lift of $`T`$”, is contained in the horizontal space. All that extends naturally to the corresponding cotangent bundles. The covariant derivative introduced in Eq. (2.7) serves just this purpose: it lifts the ordinary exterior derivative on $`_k`$ to the horizontal subspace of the larger cotangent bundle. A natural metric $`\stackrel{~}{g}`$ on the horizontal subspace induced from the flat one is given by

$$\stackrel{~}{g}(d𝒜^{},d𝒜)=g(𝒟𝒜^{},𝒟𝒜),$$
(2.8)
$`g`$ being the flat metric

$$g(d𝒜^{},d𝒜)=tr\left(\underset{l=1,2}{}dB_l^{}dB_l+dI^{}dI+dJdJ^{}\right).$$
(2.9)
One of the results of ,  is that the coefficient $`_k`$ of the $`𝒩=2`$ prepotential

$$(\mathrm{\Psi })=\frac{i}{2\pi }\mathrm{\Psi }^2\mathrm{log}\frac{2\mathrm{\Psi }^2}{e^3\mathrm{\Lambda }^2}\frac{i}{\pi }\underset{k=1}{\overset{\mathrm{}}{}}_k\left(\frac{\mathrm{\Lambda }}{\mathrm{\Psi }}\right)^{4k}\mathrm{\Psi }^2,$$
(2.10)
quoted here as a function of the $`𝒩=2`$ chiral vector superfield $`\mathrm{\Psi }`$ can be represented as an integral over the reduced moduli space $`_k^{}`$ \- that is the space from which the coordinates of the instanton center and its fermionic counterparts have been eliminated,

$$_k__k^{}e^{d_x\omega },$$
(2.11)
where $`\omega `$ is here the differential one-form

$$\omega =\mathrm{}e\left(𝒟I\overline{v}I^{}+J^{}v𝒟J\right),$$
(2.12)
$`v`$ denotes the Higgs field (the scalar component of the above mentioned superfield $`\mathrm{\Psi }`$) vacuum expectation value breaking the gauge group $`SU(2)`$ down to $`U(1)`$. We will chose $`v`$ to be of the form

$$v=\left(\begin{array}{cc}ia& 0\\ 0& ia\end{array}\right),$$
(2.13)
where $`a`$ is real. $`d_x`$ in Eq. (2.11) stands for an equivariant differential

$$d_xd+i_x,$$
(2.14)
$`d`$ being the ordinary exterior derivative and $`i_x`$ meaning contraction with the $`U(1)`$ vector field (denoted below by $`x`$ ) going along with the infinitesimal transformation

$`\delta B_i0;\delta IIv;\delta JvJ.`$ (2.15)
It is worthwhile to note, that the one-form $`\omega `$ is dual to the vector field (2.15) with respect to the metric $`\stackrel{~}{g}`$:

$$\omega =\mathrm{}e\stackrel{~}{g}(x^{},d𝒜).$$
(2.16)
The coefficient $`_k`$ may be deformed into

$$_k(t)__k^{}e^{\frac{1}{t}d_x\omega }$$
(2.17)
and we obtain

$$\frac{d}{dt}_k(t)=\frac{1}{t^2}__k^{}d_x\left(\omega e^{\frac{1}{t}d_x\omega }\right)==\frac{1}{t^2}__k^{}d\left(\omega e^{\frac{1}{t}d_x\omega }\right).$$
(2.18)
For the equality (2.18) use has been made of the equivariance of the one-form $`\omega `$,

$$d_x^2\omega =\left(di_x+i_xd\right)\omega =_x\omega =0,$$
(2.19)
where $`_x`$ is the Lie-derivative of the above introduced $`U(1)`$ vector field. To understand the second equality in (2.18) one has to note that the the integral over $`_k^{}`$ has to be taken with the top form from $`\mathrm{exp}(\frac{1}{t}d_x\omega )`$ inserted, which can be reached (being a top form ) only by $`d`$, not by $`i_x`$. The $`t`$-dependence of $`_k(t)`$ hinges on whether the total derivative integral in (2.18) picks up boundary terms. Boundary terms are apparently present in the unregularized version of $`_k^{}`$ defined by Eq.’s (2.1), (2.2) (see , ), but are not present, as has been convincingly argued by Hollowood , after $`\zeta `$-regularization555One may worry that the non-compactness of $`_k^{}`$ may create problems. This is not the case . Estimates involving the dilute gas approximation for instantons show that the integrands decay fast enough into the non-compact directions to avoid trouble.. $`_k`$ does nevertheless not depend on $`\zeta `$. Indeed, the integral (2.18) over manifolds $`_k^{}`$ associated to different parameters $`\zeta `$ are related by a simple rescaling which can be absorbed into a change of the parameter $`t`$ and so does not, according to Hollowood, alter $`_k`$.

The result

$$\frac{d_k}{dt}|_{\zeta 0}=0$$

suggests a saddle point evaluation of the integral over $`_k^{}`$, which should render an exact result. The saddle points are determined by $`i_x\omega `$,

$`i_x\omega =\stackrel{~}{g}(x^{},x)=tr({\displaystyle \underset{l=1,2}{}}[𝐋^1\mathrm{\Lambda },B_l^{}][𝐋^1\mathrm{\Lambda },B_l]`$
$`+\left(vI^{}I^{}𝐋^1\mathrm{\Lambda }\right)\left(Iv+𝐋^1\mathrm{\Lambda }I\right)`$
$`+(J^{}v+𝐋^1\mathrm{\Lambda }J^{})(vJJ𝐋^1\mathrm{\Lambda })).`$ (2.20)
The r.h.s. of Eq. (2.20) is a sum of squares, each of which should vanish for itself at a saddle point. In this way we are lead to the saddle point equalities

$`[𝐋^1\mathrm{\Lambda },B_l]=0;`$
$`Iv+𝐋^1\mathrm{\Lambda }I=0;`$
$`vJJ𝐋^1\mathrm{\Lambda }=0.`$ (2.21)
One should note that the $`𝐋^1\mathrm{\Lambda }`$ appearing in these equations are emerging from the $`U(k)`$ connections of the $`U(1)`$ vector field due to the horizontal lift of the latter. The interpretation of Eq.’s (2.21) is obvious: The saddle points are identical with the places at which the lifted vector field vanishes and can be shown, , to coincide with the union of direct products of two $`U(1)`$ instanton configurations with topological charges $`k_1`$ and $`k_2`$ satisfying the condition $`k_1+k_2=k`$. Using the general localization formula (see , p. 216 ff.) we can reduce the $`8k4`$ dimensional integral (2.11) over the moduli space $`_k^{}`$ to a sum of $`4k_1+4k_24=4k4`$ dimensional integrals over the spaces $`_{k_1,k_2}^{}_{k_1}(U(1))\times _{k_2}(U(1))\backslash 𝐂^2`$ (the common center of two $`U(1)`$ instantons is eliminated)

$$__k^{}e^{d_x\omega }=\underset{\genfrac{}{}{0pt}{}{_{k_1,k_2}}{k_1+k_2=1}}{}_{_{k_1,k_2}}\frac{1}{det^{1/2}\left(_𝒩+_𝒩\right)},$$
(2.22)
where $`_𝒩`$ is the action of the vector field on the subspaces of horizontal spaces orthogonal to $`_{k_1,k_2}`$ (these subspaces constitute the so called normal bundle) and $`_𝒩`$, a two form along $`_{k_1,k_2}`$ and a linear operator acting on these orthogonal subspaces, is the curvature of the normal bundle. Eq. (2.22) shows that each SW coefficient $`_k`$ is equal to a certain polynomial of the characteristic classes of the normal bundle. The formula (2.22) providing a drastic simplification of the initial problem, still does not allow us to perform computations beyond $`2`$-instantons666The $`2`$-instanton calculation goes smoothly and gives the desired result..

## 3 The modified vector field

We have seen in section 2, that the main building block of the superinstanton action is the vector field $`x`$ (see (2.11), (2.12), (2.16)). The difficulty we have encountered was a too large zero set of this vector field. We describe a natural deformation of this vector field whose zero set is only a discrete finite set. Such a vectorfield has been for good use in mathematics (cf. , and references therein777We became aware of the relevance of Nakajima’s work through .). We will heavily rely on some of the results of Nakajima in the following. Consider two independent rotations of space-time, first on the $`x_1`$, $`x_2`$ plane with the rotation angle $`ϵ_1`$ and the second one on the plane $`x_3`$, $`x_4`$ with rotation angle $`ϵ_2`$. It is convenient to introduce complex coordinates $`z_1=x_1+ix_2`$, $`z_2=x_3+ix_4`$ in (euclidean) space-time. The group element specified by the parameters $`(t_1,t_2)`$ ($`t_l\mathrm{exp}iϵ_l`$, $`l=1,2`$) acts on $`(z_1,z_2)`$ as $`(z_1,z_2)(t_1z_1,t_2z_2)`$. The respective action on the ADHM data is given by:

$$B_lt_lB_l;II;Jt_1t_2J.$$
(3.1)
The unbroken $`U(1)`$ subgroup of the gauge group acts as:

$$B_lB_l;IIt_v;Jt_v^1J,$$
(3.2)
where $`t_v=\mathrm{exp}ia\sigma _3`$ (our vector field $`x`$ is just the generator of the transformations (3.2) ). It is evident that the transformations (3.2), (3.1) act properly also on the moduli space $`_k`$, because they commute with $`U(k)`$. Let us combine the transformations (3.2), (3.1)

$$B_lt_lB_l;IIt_v;Jt_1t_2t_v^1J$$
(3.3)
and consider its generating vector field - to be denoted $`\stackrel{~}{x}`$ \- in the construction of a modified superinstanton action as the main building block substituting the vector field $`x`$. The components of the vector field $`\stackrel{~}{x}`$ are determined through the infinitesimal form of the transformations (3.3)

$$\delta B_lϵ_lB_l;\delta IIv;\delta JϵvJ,$$
(3.4)
where $`ϵ=ϵ_1+ϵ_2`$. To find the deformed one-form $`\stackrel{~}{\omega }`$ we substitute in Eq. (2.16) $`x`$ by $`\stackrel{~}{x}`$. As a result we get

$$\stackrel{~}{\omega }=\mathrm{}etr\left(i\underset{l=1,2}{}ϵ_l𝒟B_lB_l^{}𝒟IvI^{}+J^{}\left(viϵ\right)𝒟J\right).$$
(3.5)
For the bosonic part of the deformed superinstanton action $`d_{\stackrel{~}{x}}\stackrel{~}{\omega }`$ we obtain

$`i_{\stackrel{~}{x}}\stackrel{~}{\omega }=\stackrel{~}{g}(\stackrel{~}{x}^{},\stackrel{~}{x})=tr({\displaystyle \underset{l=1,2}{}}(iϵ_lB_l^{}+[𝐋^1\stackrel{~}{\mathrm{\Lambda }},B_l^{}])(iϵ_lB_l+[𝐋^1\stackrel{~}{\mathrm{\Lambda }},B_l])`$
$`+\left(vI^{}I^{}𝐋^1\stackrel{~}{\mathrm{\Lambda }}\right)\left(Iv+𝐋^1\stackrel{~}{\mathrm{\Lambda }}I\right)`$
$`+(J^{}(viϵ)+𝐋^1\stackrel{~}{\mathrm{\Lambda }}J^{})((v+iϵ)JJ𝐋^1\stackrel{~}{\mathrm{\Lambda }})),`$ (3.6)
where

$$\stackrel{~}{\mathrm{\Lambda }}=i_{\stackrel{~}{x}}𝐗=\underset{l=1,2}{}iϵ_l[B_l^{},B_l]+J^{}\left(v+iϵ\right)JIvI^{}h.c..$$
(3.7)
The deformed super-instanton action has an explicit dependence on the instanton center. Hence we are not allowed to restrict the moduli space of instantons imposing the condition $`trB_1=trB_2=0`$ which means that the center of instantons has zero coordinates. We define the deformed version of the Eq. (2.11) as

$$Z_k(a,ϵ_1,ϵ_2)__ke^{d_{\stackrel{~}{x}}\stackrel{~}{\omega }},$$
(3.8)
where integration is over the entire moduli space $`_k`$. The evaluation of this integral by means of localization technic is much simpler, because on the right hand side of the localization formula (2.22) we will have a finite sum over the zero locus set of vector field $`\stackrel{~}{x}`$ and therefore the curvature term of the normal bundle $`_𝒩`$ is absent. A zero locus of the vector field $`\stackrel{~}{x}`$ (equivalently a fixed point of the combined action (3.3)) on the moduli space) is defined by the conditions:

$$t_lB_l=g^1B_lg;It_v=g^1I;t_1t_2t_v^1J=Jg.$$
(3.9)
Using the freedom of $`U(k)`$ transformations we will assume that the group element $`g`$ out of $`U(k)`$ is diagonal

$$g=\left(\begin{array}{ccc}e^{i\mathrm{\Phi }_1}& \mathrm{}& 0\\ & \mathrm{}& \\ 0& \mathrm{}& e^{i\mathrm{\Phi }_k}\end{array}\right).$$
(3.10)
The fixed point condition (3.9) implies:

$`\delta B_{l,ij}\left(\mathrm{\Phi }_{ij}+ϵ_l\right)B_{l,ij}=0;`$
$`\delta I_{i\lambda }\left(\mathrm{\Phi }_i+a_\lambda \right)I_{i\lambda }=0;`$
$`\delta J_{\lambda i}\left(\mathrm{\Phi }_ia_\lambda +ϵ\right)J_{\lambda i}=0;`$
$`i,j=1,2,\mathrm{},k;l,\lambda =1,2;\mathrm{\Phi }_{ij}\mathrm{\Phi }_i\mathrm{\Phi }_j.`$ (3.11)
The geometrical meaning of these equations (as before for $`x`$) is the criticality of the horizontally lifted vector field $`\stackrel{~}{x}^{}`$ whose action in the neighborhood of a critical point is given by equation

$`\delta _{\stackrel{~}{x}^{}}B_{l,ij}\left(\mathrm{\Phi }_{ij}+ϵ_l\right)B_{l,ij};`$
$`\delta _{\stackrel{~}{x}^{}}I_{i\lambda }\left(\mathrm{\Phi }_i+a_\lambda \right)I_{i\lambda };`$
$`\delta _{\stackrel{~}{x}^{}}J_{\lambda i}\left(\mathrm{\Phi }_ia_\lambda +ϵ\right)J_{\lambda i}.`$ (3.12)
The classification of fixed points of the torus action (3.1) in the case of $`U(1)`$ (non-commutative) instantons can be found in  chapter 5.2. Fortunately, to adapt this to our case of $`SU(2)`$ (or $`SU(N)`$ in general) instantons and the combined action (3.3) only some minor modifications are needed.

There is one-to-one correspondence between the fixed points of $`\stackrel{~}{x}^{}`$ and ordered pairs of Young diagrams $`(Y_1,Y_2)`$ with $`k_1`$ and $`k_2`$ boxes resp. , s.t. $`k_1+k_2=k`$. Denote the number of boxes in the rows of the Young diagram $`Y_l`$ by $`\nu _{l,1}\nu _{l,2}\mathrm{}0`$ and the number of boxes in columns888Our convention is to numerate rows beginning from the top to the bottom and columns from left to right. by $`\nu _{l_,1}^{}\nu _{l,2}^{}\mathrm{}0`$. Then we distribute the phases $`\mathrm{\Phi }_1,\mathrm{\Phi }_2,\mathrm{},\mathrm{\Phi }_{k_1}`$ in the Young diagrams $`Y_1`$ and $`\mathrm{\Phi }_{k_1+1},\mathrm{\Phi }_{k_1+2},\mathrm{},\mathrm{\Phi }_k`$ in $`Y_2`$ subsequently filling row after row beginning from the upper left corner boxes. To the phase $`\mathrm{\Phi }_m`$ distributed in a box of $`Y_l`$ in the $`i`$-th row and $`j`$-th column we assign the value

$$\mathrm{\Phi }_m=a_l(i1)ϵ_2(j1)ϵ_1,$$
(3.13)
where $`a_1=a_2=a`$. The matrix element $`B_{1,mn}`$ ($`B_{2,mn}`$) is nonzero if and only if $`\mathrm{\Phi }_m`$ and $`\mathrm{\Phi }_n`$ are neighbors in a row (column). The only nonzero matrix elements of $`I`$ are $`I_{1,1}`$ and $`I_{k_1+1,2}`$. $`J`$ vanishes identically. These are some special points belonging to the manifold $`_{k_1,k_2}`$ introduced earlier. It is easy to see that the conditions (3.11) are fulfilled. Of course, to determine the actual values of nonzero matrix elements one should impose the ADHM equations (2.1), (2.4) but this is a more delicate problem. Below we will solve this problem for instanton charges $`k=2,3,4`$.

The next step is the construction of tangent spaces passing through the fixed points and to calculate the determinants of $`\stackrel{~}{x}^{}`$-action on this spaces. To do this, one needs to find all solutions of the linearized ADHM equations (the so called fermionic ADHM equations)

$`[\delta B_1,B_2]+[B_1,\delta B_2]+\delta IJ+I\delta J=0;`$
$`{\displaystyle \underset{l=1,2}{}}[\delta B_l,B_l^{}]+\delta II^{}J^{}\delta J=0`$ (3.14)
around the fixed points. Denote the fixed point ADHM data corresponding to the Young diagrams $`(Y_1,Y_2)`$ as $`𝒜_{Y_1,Y_2}`$. The function $`Z_k(a,ϵ_1,ϵ_2)`$ can be completely determined in terms of above data:

$`Z_k(a,ϵ_1,ϵ_2)={\displaystyle \underset{Y_1,Y_2}{}}{\displaystyle \frac{1}{det_{\stackrel{~}{x}^{}}}}|_{𝒜_{Y_1,Y_2}}.`$ (3.15)
The reason of appearance of $`1/det`$ instead of $`1/\sqrt{det}`$ as in Eq. (2.22) is due to our convention to consider complex tangent spaces and linear operators acting on them instead of their real forms.

Though the deformed superinstanton action $`d_{\stackrel{~}{x}^{}}\stackrel{~}{\omega }`$ tends to the initial one as $`ϵ_1,ϵ_20`$, it is not true for the $`Z_k(a,ϵ_1,ϵ_2)`$ which is highly singular at this limit. In his recent paper N.Nekrasov  brings some field theoretical arguments to come to the remarkable conclusion that the generating function

$$Z(q,a,ϵ_1,ϵ_2)1+\underset{k=1}{\overset{\mathrm{}}{}}Z_k(a,ϵ_1,ϵ_2)q^k$$
(3.16)
can be represented as

$$Z(q,a,ϵ_1,ϵ_2)=\mathrm{exp}\frac{1}{ϵ_1ϵ_2}(q,a,ϵ_1,ϵ_2),$$
(3.17)
where $``$ is regular at $`ϵ_1,ϵ_2=0`$ and that the SW coefficients $`_k`$ are nothing else than the coefficients of the Taylor expansion of $`(q,a,ϵ_1,ϵ_2)`$:

$$(q,a,ϵ_1,ϵ_2)=\underset{k=1}{\overset{\mathrm{}}{}}_k(a,ϵ_1,ϵ_2)q^k$$
(3.18)
at $`ϵ_1,ϵ_2=0`$. More precisely999The factor $`2^{2k2}`$ is needed to have the normalization adopted in .

$$_ka^{24k}=2^{2k2}_k(a,ϵ_1,ϵ_2).$$
(3.19)
We hope to return to this point in a future publication and to present an intrinsic explanation of the singularity structure (3.17) using directly the definition of $`Z(q,a,ϵ_1,ϵ_2)`$.

## 4 Low charge instanton calculations

In this section we carry out explicit calculations for instanton charges up to four. For simplicity we rescale ADHM data and set $`\zeta =1`$.

* One-instantons

This case is almost trivial. We have two pairs of Young diagrams101010We always note only nonzero $`\nu `$’s.:
a) $`Y_1=\left\{\nu _{1,1}=1\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$;
b) $`Y_1=\left\{\mathrm{}\right\}`$, $`Y_2=\left\{\nu _{1,1}=1\right\}`$. In the case a) $`B_1=B_2=0`$, $`I=\left(\begin{array}{cc}1& 0\end{array}\right)`$, $`\mathrm{\Phi }_1=a_1`$. The tangent space vectors are given by

$`\delta B_l=\left(\delta B_{l,11}\right);\delta I=\left(\begin{array}{cc}0& \delta I_{12}\end{array}\right);\delta J=\left(\begin{array}{c}\delta J_{11}\\ 0\end{array}\right).`$ (4.4)
The dimension of the tangent space is $`2\times 4=8`$ as it should. Taking into account (3.12) one easily calculates the determinant

$$det_{\stackrel{~}{x}^{}}=ϵ_1ϵ_2a_{21}\left(a_{12}+ϵ\right),$$
(4.5)
where $`a_{\lambda \mu }a_\lambda a_\mu ;\lambda ,\mu =1,2`$. There is no need to carry out calculation for the case b) because for interchanged $`Y_1`$ and $`Y_2`$ one obtains the same determinant with interchanged $`a_1`$ and $`a_2`$. Note also that simultaneous transposition of both Young diagrams gives rise to a determinant with interchanged $`ϵ_1ϵ_2`$. Below we will explicitly mention only one pair of such symmetry related diagrams, but of course we will take all of them into account in final expressions of $`Z_k`$ Eq. (3.15). Thus for the $`1`$-instantons

$$Z_1(a,ϵ_1,ϵ_2)=\left(ϵ_1ϵ_2a_{21}\left(a_{12}+ϵ\right)\right)^1+\left(ϵ_1ϵ_2a_{12}\left(a_{21}+ϵ\right)\right)^1.$$
(4.6)
From Eq. (3.16)-(3.18)

$$_1(a,ϵ_1,ϵ_2)=ϵ_1ϵ_2Z_1=\frac{2}{a_{12}^2ϵ^2}.$$
(4.7)
Taking the limit $`ϵ_{1,2}0`$ one obtains

$$_1=1/2.$$

* 2-instantons

There are two basic cases.
a) $`Y_1=\{\nu _{1,1}=2)\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$.
The fixed point ADHM data:

$`B_1=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right);B_2=0;I=\left(\begin{array}{cc}\sqrt{2}& 0\\ 0& 0\end{array}\right).`$ (4.12)
The tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cc}\delta B_{1,11}& \delta B_{1,12}\\ 0& \delta B_{1,11}\end{array}\right);\delta B_2=\left(\begin{array}{cc}\delta B_{2,11}& 0\\ \delta B_{2,21}& \delta B_{2,11}\end{array}\right);`$ (4.17)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\end{array}\right);\delta J=\left(\begin{array}{cc}0& 0\\ \delta J_{21}& \delta J_{22}\end{array}\right).`$ (4.22)
The determinant:

$$det_{\stackrel{~}{x}^{}}=2ϵ_1^2ϵ_2ϵ_{21}a_{21}\left(a_{21}ϵ_1\right)\left(a_{12}+ϵ\right)\left(a_{12}+2ϵ_1+ϵ_2\right).$$
(4.23)
b) $`Y_1=\{\nu _{1,1}=1)\}`$, $`Y_2=\left\{\nu _{2,1}=1\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_2`$.
Fixed point:

$`B_1=B_2=0;I=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).`$ (4.26)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cc}\delta B_{1,11}& \delta B_{1,12}\\ \delta B_{1,21}& \delta B_{1,22}\end{array}\right);\delta B_2=\left(\begin{array}{cc}\delta B_{2,11}& \delta B_{2,12}\\ \delta B_{2,21}& \delta B_{2,22}\end{array}\right);`$ (4.31)
$`\delta I`$ $`=`$ $`0;\delta J=0.`$ (4.32)
Determinant:

$$det_{\stackrel{~}{x}^{}}=ϵ_1^2ϵ_2^2\left(a_{12}^2ϵ_1^2\right)\left(a_{12}^2ϵ_2^2\right).$$
(4.33)
Using above data one first calculates $`Z_2`$ and then

$$_2(a,ϵ_1,ϵ_2)=ϵ_1ϵ_2\left(Z_2\frac{1}{2}Z_1^2\right).$$

The final result reads:

$$_2(a,ϵ_1,ϵ_2)=\frac{20a^2+7ϵ_1^2+16ϵ_1ϵ_2+7ϵ_2^2}{\left(4a^2ϵ^2\right)^2\left(4a^2\left(2ϵ_1+ϵ_2\right)^2\right)\left(4a^2\left(ϵ_1+2ϵ_2\right)^2\right)}.$$
(4.34)
Now let us take the limit $`ϵ_{1,2}0`$:

$$_2=\frac{5}{2^4}.$$

* 3-instantons

There are three basic cases.
a) $`Y_1=\left\{\nu _{1,1}=3\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_12ϵ_1`$.
Fixed point:

$`B_1=\left(\begin{array}{ccc}0& 0& 0\\ \sqrt{2}& 0& 0\\ 0& 1& 0\end{array}\right);B_2=0;I=\left(\begin{array}{cc}\sqrt{3}& 0\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.41)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{ccc}\delta B_{1,11}& \delta B_{1,12}& \delta B_{1,13}\\ 0& \delta B_{1,11}& \frac{1}{\sqrt{2}}\delta B_{1,12}\\ 0& 0& \delta B_{1,11}\end{array}\right);\delta B_2=\left(\begin{array}{ccc}\delta B_{2,11}& 0& 0\\ \delta B_{2,21}& \delta B_{2,11}& 0\\ \delta B_{2,31}& \frac{1}{\sqrt{2}}\delta B_{2,21}& \delta B_{2,11}\end{array}\right);`$ (4.48)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\\ 0& \delta I_{32}\end{array}\right);\delta J=\left(\begin{array}{ccc}0& 0& 0\\ \delta J_{21}& \delta J_{22}& \delta J_{23}\end{array}\right).`$ (4.54)
The determinant:

$`det_{\stackrel{~}{x}^{}}=6ϵ_1^3ϵ_2ϵ_{21}\left(ϵ_22ϵ_1\right)a_{21}\left(a_{21}ϵ_1\right)\left(a_{21}2ϵ_1\right)`$
$`\times \left(a_{12}+ϵ\right)\left(a_{12}+2ϵ_1+ϵ_2\right)\left(a_{12}+3ϵ_1+ϵ_2\right).`$ (4.55)
b) $`Y_1=\left\{\nu _{1,1}=2,\nu _{1,2}=1\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_1ϵ_1ϵ_2`$.
Fixed point:

$`B_1=\left(\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right);B_2=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 1& 0& 0\end{array}\right);I=\left(\begin{array}{cc}\sqrt{3}& 0\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.65)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{ccc}\delta B_{1,11}& 0& 0\\ 0& \delta B_{1,22}& \delta B_{2,22}\delta B_{2,11}\\ 0& \delta B_{1,32}& 2\delta B_{1,11}\delta B_{1,22}\end{array}\right);`$ (4.69)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{ccc}\delta B_{2,11}& 0& 0\\ 0& \delta B_{2,22}& \delta B_{2,23}\\ 0& \delta B_{1,11}\delta B_{1,22}& 2\delta B_{211}\delta B_{2,22}\end{array}\right);`$ (4.73)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\\ 0& \delta I_{32}\end{array}\right);\delta J=\left(\begin{array}{ccc}0& 0& 0\\ \delta J_{21}& \delta J_{22}& \delta J_{23}\end{array}\right).`$ (4.79)
Determinant:

$`det_{\stackrel{~}{x}^{}}=ϵ_1^2ϵ_2^2\left(2ϵ_1ϵ_2\right)\left(2ϵ_2ϵ_1\right)a_{21}\left(a_{21}ϵ_1\right)\left(a_{21}ϵ_2\right)`$
$`\times \left(a_{12}+ϵ\right)\left(a_{12}+2ϵ_1+ϵ_2\right)\left(a_{12}+2ϵ_2+ϵ_1\right).`$ (4.80)
c) $`Y_1=\left\{\nu _{1,1}=1,\nu _{1,2}=1\right\}`$, $`Y_2=\left\{\nu _{2,1}=1\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_2`$, $`\mathrm{\Phi }_3=a_2`$.
Fixed point:

$`B_1=0;B_2=\left(\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right);I=\left(\begin{array}{cc}\sqrt{2}& 0\\ 0& 0\\ 0& 1\end{array}\right).`$ (4.87)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{ccc}\delta B_{1,11}& 0& 0\\ \delta B_{1,21}& \delta B_{1,11}& \delta B_{1,23}\\ \delta B_{1,31}& \delta B_{1,32}& \delta B_{1,33}\end{array}\right);`$ (4.91)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{ccc}\delta B_{2,11}& \delta B_{2,12}& \delta B_{2,13}\\ 0& \delta B_{2,11}& \delta B_{2,23}\\ 0& \delta B_{2,32}& \delta B_{2,33}\end{array}\right);`$ (4.95)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta B_{2,23}\\ 0& 0\\ 0& 0\end{array}\right);\delta J=\left(\begin{array}{ccc}0& 0& 0\\ \delta B_{1,32}& 0& 0\end{array}\right).`$ (4.101)
Determinant:

$`det_{\stackrel{~}{x}^{}}=2ϵ_1^2ϵ_2^3ϵ_{12}a_{21}\left(a_{21}+ϵ_1ϵ_2\right)\left(a_{12}+ϵ_1\right)\left(a_{21}+ϵ_2\right)`$
$`\times \left(a_{12}+ϵ\right)\left(a_{12}+2ϵ_2\right).`$ (4.102)
Using these data we have calculated $`Z_3`$ and then

$$_3(a,ϵ_1,ϵ_2)=ϵ_1ϵ_2\left(Z_3+\frac{1}{3}Z_1^3Z_1Z_2\right).$$

Here is the final result:

$`_3(a,ϵ_1,ϵ_2)=`$ (4.103)
$`{\displaystyle \frac{16\left(144a^4+29ϵ_1^4+154ϵ_1^3ϵ_2+258ϵ_1^2ϵ_2^2+154ϵ_1ϵ_2^3+29ϵ_2^4+8a^2\left(29ϵ_1^2+71ϵ_1ϵ_2+29ϵ_2^2\right)\right)}{3\left(4a^2ϵ^2\right)^3\left(4a^2\left(2ϵ_1+ϵ_2\right)^2\right)\left(4a^2\left(3ϵ_1+ϵ_2\right)^2\right)\left(4a^2\left(ϵ_1+2ϵ_2\right)^2\right)\left(4a^2\left(ϵ_1+3ϵ_2\right)^2\right)}}.`$
Now let us take the limit $`ϵ_{1,2}0`$:

$$_3=\frac{3}{4}.$$

* 4-instantons

Now we need to investigate 7 different cases.
a) $`Y_1=\left\{\nu _{1,1}=4\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_12ϵ_1`$, $`\mathrm{\Phi }_4=a_13ϵ_1`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ \sqrt{3}& 0& 0& 0\\ 0& \sqrt{2}& 0& 0\\ 0& 0& 1& 0\end{array}\right);B_2=0;I=\left(\begin{array}{cc}2& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.112)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& \delta B_{1,12}& \delta B_{1,13}& \delta B_{1,14}\\ 0& \delta B_{1,11}& \sqrt{\frac{2}{3}}\delta B_{1,12}& \frac{1}{\sqrt{3}}\delta B_{1,13}\\ 0& 0& \delta B_{1,11}& \frac{1}{\sqrt{3}}\delta B_{1,12}\\ 0& 0& 0& \delta B_{1,11}\end{array}\right);`$ (4.117)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& 0& 0\\ \delta B_{2,21}& \delta B_{211}& 0& 0\\ \delta B_{2,31}& \sqrt{\frac{2}{3}}\delta B_{2,21}& \delta B_{2,11}& 0\\ \delta B_{2,41}& \frac{1}{\sqrt{3}}\delta B_{2,31}& \frac{1}{\sqrt{3}}\delta B_{2,21}& \delta B_{2,11}\end{array}\right);`$ (4.122)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\\ 0& \delta I_{32}\\ 0& \delta I_{42}\end{array}\right);\delta J=\left(\begin{array}{cccc}0& 0& 0& 0\\ \delta J_{21}& \delta J_{22}& \delta J_{23}& \delta J_{24}\end{array}\right).`$ (4.129)
The determinant:

$`det_{\stackrel{~}{x}^{}}=24ϵ_1^4{\displaystyle \underset{j=0}{\overset{3}{}}}\left(ϵ_2jϵ_1\right)\left(a_{21}jϵ_1\right)\left(a_{12}+ϵ+jϵ_1\right).`$ (4.130)
b) $`Y_1=\left\{\nu _{1,1}=3,\nu _{1,2}=1\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_12ϵ_1`$, $`\mathrm{\Phi }_4=a_1ϵ_2`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ \sqrt{2}& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right);B_2=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\end{array}\right);I=\left(\begin{array}{cc}2& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.143)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& \delta B_{1,12}& 0& 0\\ 0& \delta B_{1,22}& \sqrt{2}\delta B_{1,12}& \sqrt{2}\left(\delta B_{2,22}\delta B_{2,11}\right)\\ 0& 0& \delta B_{1,22}& \delta B_{1,34}\\ 0& 0& \delta _{1,43}& 3\delta B_{1,11}2\delta B_{1,22}\end{array}\right);`$ (4.148)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& 0& 0\\ \delta B_{1,34}& \delta B_{2,22}& 0& 0\\ 0& \sqrt{2}\delta B_{1,34}& \delta B_{2,22}& \delta B_{2,34}\\ 0& \sqrt{2}\left(\delta B_{1,11}\delta B_{1,22}\right)& \delta B_{1,12}& 3\delta B_{2,11}2\delta B_{2,22}\end{array}\right);`$ (4.153)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\\ 0& \delta I_{32}\\ 0& \delta I_{42}\end{array}\right);\delta J=\left(\begin{array}{cccc}0& 0& 0& 0\\ \delta J_{21}& \delta J_{22}& \delta J_{23}& \delta J_{24}\end{array}\right).`$ (4.160)
The determinant:

$`det_{\stackrel{~}{x}^{}}=4ϵ_1^3ϵ_2^2ϵ_{21}^2\left(3ϵ_1ϵ_2\right)a_{21}\left(a_{21}ϵ_1\right)\left(a_{21}2ϵ_1\right)\left(a_{21}ϵ_2\right)\left(a_{12}+ϵ\right)`$
$`\times \left(a_{12}+2ϵ_1+ϵ_2\right)\left(a_{12}+3ϵ_1+ϵ_2\right)\left(a_{12}+ϵ_1+2ϵ_2\right).`$ (4.161)
c) $`Y_1=\left\{\nu _{1,1}=2,\nu _{1,2}=2\right\}`$, $`Y_2=\left\{\mathrm{}\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_1ϵ_2`$, $`\mathrm{\Phi }_4=a_1ϵ`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ \sqrt{\frac{3}{2}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& \frac{1}{\sqrt{2}}& 0\end{array}\right);B_2=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ \sqrt{\frac{3}{2}}& 0& 0& 0\\ 0& \frac{1}{\sqrt{2}}& 0& 0\end{array}\right);I=\left(\begin{array}{cc}2& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.174)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& \delta B_{1,12}& 0& 0\\ 0& \delta B_{1,11}& 0& 0\\ \delta B_{1,31}& \delta B_{1,32}& \delta B_{1,11}& \sqrt{3}\delta B_{1,12}\\ 0& \sqrt{3}\delta B_{1,31}& 0& \delta B_{1,11}\end{array}\right);`$ (4.179)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& \delta B_{2,13}& 0\\ \delta B_{2,21}& \delta B_{2,11}& \delta B_{2,23}& \sqrt{3}\delta B_{2,13}\\ 0& 0& \delta B_{2,11}& 0\\ 0& 0& \sqrt{3}\delta B_{2,21}& \delta B_{2,11}\end{array}\right);`$ (4.184)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \delta I_{12}\\ 0& \delta I_{22}\\ 0& \delta I_{32}\\ 0& \delta I_{42}\end{array}\right);\delta J=\left(\begin{array}{cccc}0& 0& 0& 0\\ \delta J_{21}& \delta J_{22}& \delta J_{23}& \delta J_{24}\end{array}\right).`$ (4.191)
The determinant:

$`det_{\stackrel{~}{x}^{}}=4ϵ_1^2ϵ_2^2ϵ_{12}^2\left(2ϵ_1ϵ_2\right)\left(2ϵ_2ϵ_1\right)a_{21}\left(a_{21}ϵ_1\right)\left(a_{21}ϵ_2\right)\left(a_{12}+ϵ\right)^2`$
$`\times \left(a_{12}+2ϵ_1+ϵ_2\right)\left(a_{12}+2ϵ_2+ϵ_1\right)\left(a_{12}+2ϵ\right).`$ (4.192)
d) $`Y_1=\left\{\nu _{1,1}=1\right\}`$, $`Y_2=\left\{\nu _{2,1}=3\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_2`$, $`\mathrm{\Phi }_3=a_2ϵ_1`$, $`\mathrm{\Phi }_4=a_22ϵ_1`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& \sqrt{2}& 0& 0\\ 0& 0& 1& 0\end{array}\right);B_2=0;I=\left(\begin{array}{cc}1& 0\\ 0& \sqrt{3}\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.201)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& 0& 0& \delta B_{1,14}\\ \delta B_{1,21}& \delta B_{1,22}& \delta B_{1,23}& \delta B_{1,24}\\ \delta B_{1,31}& 0& \delta B_{1,22}& \frac{1}{\sqrt{2}}\delta B_{1,23}\\ \delta B_{1,41}& 0& 0& \delta B_{1,22}\end{array}\right);`$ (4.206)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& \delta B_{2,12}& \delta B_{2,13}& \delta B_{2,14}\\ 0& \delta B_{2,22}& 0& 0\\ 0& \delta B_{2,32}& \delta B_{2,22}& 0\\ \delta B_{2,41}& \delta B_{2,42}& \frac{1}{\sqrt{2}}\delta B_{2,32}& \delta B_{2,22}\end{array}\right);`$ (4.211)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ \sqrt{2}\delta B_{1,31}& 0\\ \delta B_{1,41}& 0\\ 0& 0\end{array}\right);\delta J=\left(\begin{array}{cccc}0& \sqrt{2}\delta B_{2,13}& \delta B_{2,14}& 0\\ 0& 0& 0& 0\end{array}\right).`$ (4.218)
The determinant:

$`det_{\stackrel{~}{x}^{}}=6ϵ_1^4ϵ_2^2ϵ_{21}\left(ϵ_22ϵ_1\right)a_{12}\left(a_{21}+3ϵ_1\right)\left(a_{12}^2ϵ_1^2\right)\left(a_{21}+ϵ_2\right)`$
$`\times \left(a_{21}+ϵ\right)\left(a_{21}+2ϵ_1+ϵ_2\right)\left(a_{12}2ϵ_1+ϵ_2\right).`$ (4.219)
e) $`Y_1=\left\{\nu _{1,1}=1\right\}`$, $`Y_2=\left\{\nu _{2,1}=2,\nu _{2,2}=1\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_2`$, $`\mathrm{\Phi }_3=a_2ϵ_1`$, $`\mathrm{\Phi }_4=a_2ϵ_2`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right);B_2=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right);I=\left(\begin{array}{cc}1& 0\\ 0& \sqrt{3}\\ 0& 0\\ 0& 0\end{array}\right).`$ (4.232)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& 0& \delta B_{1,13}& \delta B_{1,14}\\ 0& \delta B_{1,22}& 0& 0\\ \delta B_{1,31}& 0& \delta B_{1,33}& \delta B_{2,33}\delta B_{2,22}\\ \delta B_{1,41}& 0& \delta B_{1,43}& 2\delta B_{1,22}\delta B_{1,33}\end{array}\right);`$ (4.237)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& \delta B_{2,13}& \delta B_{2,14}\\ 0& \delta B_{2,22}& 0& 0\\ \delta B_{2,31}& 0& \delta B_{2,33}& \delta B_{2,34}\\ \delta B_{2,41}& 0& \delta B_{1,22}\delta B_{1,33}& 2\delta B_{2,22}\delta B_{2,33}\end{array}\right);`$ (4.242)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ \delta B_{1,31}+\delta B_{2,41}& 0\\ 0& 0\\ 0& 0\end{array}\right);\delta J=\left(\begin{array}{cccc}0& \delta B_{2,13}\delta B_{1,14}& 0& 0\\ 0& 0& 0& 0\end{array}\right).`$ (4.249)
The determinant:

$`det_{\stackrel{~}{x}^{}}=ϵ_1^3ϵ_2^3\left(2ϵ_1ϵ_2\right)\left(2ϵ_2ϵ_1\right)a_{12}^2\left(a_{21}+2ϵ_1\right)\left(a_{21}+2ϵ_2\right)\left(a_{21}+ϵ\right)^2\left(a_{12}^2ϵ_{12}^2\right).`$ (4.250)
f) $`Y_1=\left\{\nu _{1,1}=2\right\}`$, $`Y_2=\left\{\nu _{2,1}=2\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_2`$, $`\mathrm{\Phi }_4=a_2ϵ_1`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);B_2=0;I=\left(\begin{array}{cc}\sqrt{2}& 0\\ 0& 0\\ 0& \sqrt{2}\\ 0& 0\end{array}\right).`$ (4.259)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& \delta B_{1,12}& \delta B_{1,13}& \delta B_{1,14}\\ 0& \delta B_{1,11}& 0& \delta B_{1,13}\\ \delta B_{1,31}& \delta B_{1,32}& \delta B_{1,33}& \delta B_{1,34}\\ 0& \delta B_{1,31}& 0& \delta B_{1,33}\end{array}\right);`$ (4.264)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& \delta B_{2,13}& 0\\ \delta B_{2,21}& \delta B_{2,11}& \delta B_{2,23}& \delta B_{2,13}\\ \delta B_{2,31}& 0& \delta B_{2,33}& 0\\ \delta B_{2,41}& \delta B_{2,31}& \delta B_{2,43}& \delta B_{2,33}\end{array}\right);`$ (4.269)
$`\delta I`$ $`=`$ $`0;\delta J=0.`$ (4.270)
The determinant:

$`det_{\stackrel{~}{x}^{}}=4ϵ_1^4ϵ_2^2ϵ_{12}^2\left(a_{12}^2ϵ_1^2\right)\left(a_{12}^24ϵ_1^2\right)\left(a_{12}^2ϵ_2^2\right)\left(a_{12}^2ϵ_{12}^2\right).`$ (4.271)
g) $`Y_1=\left\{\nu _{1,1}=2\right\}`$, $`Y_2=\left\{\nu _{2,1}=1,\nu _{2,2}=1\right\}`$; $`\mathrm{\Phi }_1=a_1`$, $`\mathrm{\Phi }_2=a_1ϵ_1`$, $`\mathrm{\Phi }_3=a_2`$, $`\mathrm{\Phi }_4=a_2ϵ_2`$.
Fixed point:

$`B_1=\left(\begin{array}{cccc}0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right);B_2=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);I=\left(\begin{array}{cc}\sqrt{2}& 0\\ 0& 0\\ 0& \sqrt{2}\\ 0& 0\end{array}\right).`$ (4.284)
Tangent space:

$`\delta B_1`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{1,11}& \delta B_{1,12}& \delta B_{1,13}& \delta B_{1,14}\\ 0& \delta B_{1,11}& \delta B_{1,23}& 0\\ 0& 0& \delta B_{1,33}& 0\\ 0& \delta B_{1,42}& \delta B_{1,43}& \delta B_{1,33}\end{array}\right);`$ (4.289)
$`\delta B_2`$ $`=`$ $`\left(\begin{array}{cccc}\delta B_{2,11}& 0& 0& 0\\ \delta B_{2,21}& \delta B_{2,11}& 0& \delta B_{2,24}\\ \delta B_{2,31}& \delta B_{2,32}& \delta B_{2,33}& \delta B_{2,34}\\ \delta B_{2,41}& 0& 0& \delta B_{2,33}\end{array}\right);`$ (4.294)
$`\delta I`$ $`=`$ $`\left(\begin{array}{cc}0& \frac{1}{\sqrt{2}}\delta B_{1,23}\\ 0& 0\\ \frac{1}{\sqrt{2}}\delta B_{2,41}& 0\\ 0& 0\end{array}\right);\delta J=\left(\begin{array}{cccc}0& 0& \frac{1}{\sqrt{2}}\delta B_{1,14}& 0\\ \frac{1}{\sqrt{2}}\delta B_{2,32}& 0& 0& 0\end{array}\right);`$ (4.301)
The determinant:

$`det_{\stackrel{~}{x}^{}}=4ϵ_1^3ϵ_2^3ϵ_{12}^2a_{12}^2\left(a_{12}^2ϵ^2\right)\left(a_{12}ϵ_1\right)\left(a_{12}+ϵ_2\right)\left(a_{12}+2ϵ_1ϵ_2\right)\left(a_{21}+2ϵ_2ϵ_1\right).`$ (4.302)
The expression for

$$_4(a,ϵ_1,ϵ_2)=ϵ_1ϵ_2\left(Z_4Z_1Z_3\frac{1}{2}Z_2^2+Z_1^2Z_2\frac{1}{4}Z_1^4\right)$$

is very lengthy to present here. We only note here that it is indeed regular at $`ϵ_{1,2}0`$. Here are the first nontrivial terms of its expansion:

$$_4(a,ϵ_1,ϵ_2)=\frac{1469}{2^{15}a^{14}}+\frac{18445}{2^{15}a^{16}}\left(ϵ_1^2+ϵ_2^2\right)+\frac{15151}{2^{14}a^{16}}ϵ_1ϵ_2+\mathrm{},$$
(4.303)
Thus

$$_4=\frac{1469}{2^9},$$

which is as it should be.

Finally we quote a general formula for the determinant of the vector field action on the tangent space of a generic critical point at arbitrary instanton number and gauge group $`SU(N)`$. To obtain this formula we closely follow the line of arguments presented in , Section 5.2, where the characters of the torus action around fixed points are calculated. For our purposes we need to calculate the character of the representation of the group

$$U(1)^{N1}\times U(1)^2$$
(4.304)
(the first factor is the Cartan subgroup of the group $`SU(N)`$ and the second is the 2-torus acting in space-time and on ADHM data) in the tangent space at the fixed point specified by the Young diagrams $`Y_1,\mathrm{},Y_N`$. The result reads (cf. with the formula in proposition 5.8 page 67 of Nakajima’s book ):

$`\chi (a_l,ϵ_1,ϵ_2)={\displaystyle \underset{\lambda ,\mu =1}{\overset{N}{}}}T_{a_\mu }T_{a_\lambda }^1\left({\displaystyle \underset{sY_\lambda }{}}T_1^{h_\lambda (s)}T_2^{1+v_\mu (s)}+{\displaystyle \underset{s^{}Y_\mu }{}}T_1^{1+h_\mu (s^{})}T_2^{v_\lambda (s^{})}\right),`$ (4.305)
where $`h_\lambda (s)=\nu _{\lambda ,i}j`$, $`v_\lambda (s)=\nu _{\lambda ,j}^{}i`$ if the box $`s`$ is located on the $`i`$-th row and the $`j`$-th column of a Young diagram. It is assumed that $`\nu _{\lambda ,i}`$, $`\nu _{\lambda ,j}^{}`$ are defined for arbitrary positive integers $`i`$, $`j`$. For $`i>\nu _{\lambda ,1}^{}`$ and $`j>\nu _{\lambda ,1}`$ by definition they are identically zero. In (4.305) $`T_{a_\lambda }\mathrm{exp}ia_\lambda `$, $`T_l\mathrm{exp}iϵ_l`$ are elements of respective $`U(1)`$ factors of the group (4.304) taken in the fundamental representations. A term of the form $`T_{a_\mu }T_{a_\lambda }^1T_1^mT_2^n`$ in Eq. (4.305) indicates that the tangent space includes a (complex) one dimensional invariant subspace of our deformed vector field $`\stackrel{~}{x}^{}`$ action with eigenvalue $`a_\mu a_\nu +mϵ_1+nϵ_2`$. Multiplying all these eigenvalues for the determinant of $`\stackrel{~}{x}^{}`$ action we find111111A more complicated formula for the (inverse ) determinant is presented also in  Eq. (3.20) which has a discrepancy as compared to our formula (4.306). We have checked that in all cases explicitly described in this section the formula (4.306) gives correct results.

$`det_{\stackrel{~}{x}^{}}|_{(Y_1,\mathrm{},Y_N)}={\displaystyle \underset{\lambda ,\mu =1}{\overset{N}{}}}({\displaystyle \underset{sY_\lambda }{}}(a_{\mu \lambda }ϵ_1h_\lambda (s)+ϵ_2(1+v_\mu (s)))`$
$`\times {\displaystyle \underset{s^{}Y_\mu }{}}(a_{\mu \lambda }+ϵ_1(1+h_\mu (s^{}))ϵ_2v_\lambda (s^{}))).`$ (4.306)
The repercussions of this formula will be discussed in a subsequent publication.

## 5 Acknowledgements

This work has been supported by the Volkswagen Foundation of Germany. R. P. gratefully acknowledges the partial financial support of INTAS 00-561 and the Swiss grant SCOPE.


text:
# Mixing Time Scales in a Supernova-Driven Interstellar Medium

## 1 Introduction

Measurements of the O/H ratio in the interstellar medium (ISM) along lines of sight reaching altitudes above the galactic disk of up to 1500 pc by Meyer et al. (1998) and Cartledge et al. (2001) suggest that the ISM is well-mixed up to these heights. This suggestion is also supported by the O/H ratios observed in 41 HII regions in M101 (Kennicutt & Garnett, 1996), which show only a 0.1–0.2 dex dispersion about the radial gradient. In NGC 1569 and 4214 the dispersion is extremely low (Kobulnicky & Skillman 1996, 1997). Based on these observations it seems that there exists strong evidence pointing to a fairly well-mixed ISM in galaxies including the Milky Way.

On the other hand, recent observations using the Interstellar Medium Absorption Profile Spectrograph (IMAPS) suggest that there are local variations in the D/H ratio of up to a factor of three among the studied lines of sight (Jenkins et al. 1999; Sonneborn et al. 2000; for an up to date discussion on deuterium abundances see York 2002). The detection of deuterium in the ISM indicates that it was not processed into stars and the local variations reported in the D/H ratio suggest incomplete mixing in the ISM. However, such an inhomogeneous picture contradicts the almost homogeneous distribution seen in O/H ratioes.

Deuterium was produced during big bang nucleosynthesis and is destroyed during the pre-main sequence phase of stars. Thus, its ratio to hydrogen is initially constant in primordial gas and decreases as it is processed in stars. Deuterium may also be added into galactic interstellar gas by accretion of primordial gas, as required by some chemical evolution models in order to explain its present abundance in the Galaxy (see Tosi et al. 1998 for details.). Mulan & Linsky (1999) proposed that deuterium could be formed in stellar flares from M dwarf stars. However, the production of deuterium by this mechanism should be accompanied by a significant production of lithium, which has not been observed.

Contrary to oxygen, deuterium has a primordial origin and limited destruction mechanisms, making this species a good tracer of the mixing in the ISM. Abundances of deuterium are probably not influenced by depletions onto grains, as is the case with most other species. Although there has been some suggestion that deuterium may be preferentially bonded onto the surfaces of dust grains, this effect is not sufficient to explain the large variations in D/H seen in the ISM. Furthermore, deuterium and oxygen are the only species that have an ionization fraction that closely follows that of H. Deuterium and hydrogen have almost the same ionization potentials, photoionization cross sections and recombination coefficients, making it difficult for the D/H variations to be a result of any unusual ionization effects. Other effects such as fractionation by molecular chemistry have been considered by different authors, but there is no strong evidence that such effects can be responsible for the D/H variations seen in the ISM.

We face a set of observations that point towards contradicting conclusions. O/H observations support the idea that complete mixing occurs in the ISM on timescales comparable to those for chemistry to reach equilibrium, while D/H measurements indicate a not so well mixed ISM. Solving this contradiction is an important step towards the understanding of the mixing and variation of these elements in the ISM, and therefore, of the the interpretation of local measurements as constraints of its cosmic abundance.

Roy & Kunth (1995) discussed the problem of mixing in the ISM and concluded that the ISM should appear more homogeneous than it does. These authors suggested that mixing in the ISM is driven by three main mechanisms acting on different length scales. On the Galactic scale, they argued that turbulent diffusion of interstellar clouds in the shear flow of galactic differential rotation could remove the inhomogeneities in the ISM in less than $`10^9`$ yr, in particular the fluctuations seen in O/H. At scales of $`100\mathrm{}1000`$ pc, cloud collisions and expanding supershells driven by supernovas (SNs) from evolving massive star associations, differential rotation and triggered star formation (by stellar winds and SN explosions) could re-distribute and mix inhomogeneities efficiently in about $`10^8`$ yr; and at smaller scales of the order of $`1\mathrm{}100`$ pc, turbulent diffusion may be the dominant mechanism in cold clouds, while Rayleigh-Taylor and Kelvin-Helmoltz instabilities could develop in regions of gas ionized by massive stars, leading to complete mixing in less than 2 Myr. (This analysis does not take into account that during the mixing process new inhomogeneities are added into the ISM.)

Bateman & Larson (1993) already pointed out that turbulent mixing should be important in the ISM. However, they based their analysis on a picture of cloud collisions that ignores the likely origin of clouds as density enhancements in regions of converging flow in the turbulent ISM (e.g. Vázquez-Semadeni et al. 1995).

Tenorio-Tagle (1996) noted that classical diffusion theory on the molecular level predicts mixing timescales in the ISM for gas with $`T<10^5`$ K several orders of magnitude longer than that proposed by Roy & Kunth (1995). As a consequence, Tenorio-Tagle (1996) concluded that turbulent mixing as well as hydrodynamical instabilities, as proposed by Roy & Kunth (1995), are important mixing mechanisms that allow the reduction of the long mixing time scales obtained from classical diffusion theory.

The observations showing that the ISM has an inhomogeneous distribution pose the question of how it is possible in a SN-driven ISM for different species not to be mixed well enough to suppress local variations of their ratio to hydrogen along and among lines of sight. The ISM is regulated by SN explosions, and therefore, by well-structured and explosive flows rather than just by diffuse turbulence that acts on the smaller scales. To understand the mixing process in such a medium we carry out direct numerical simulations of the evolution of the ISM polluted with inhomogeneities of different length scales. We aim to understand if classical mixing length theory can be used to predict the mixing time scales as a function of size scale.

In § 2 we review the existing theoretical ideas on dispersion, diffusion and turbulent mixing and apply them to the ISM. Section 3 describes the modified version of the three-dimensional model of de Avillez (2000), how the runs used in the current study were set up, and the global evolution of the ISM as seen in these simulations. In § 4, the evolution of the tracer fields and the mixing process in the ISM are described. Section 5 discusses the results and their implications. Finally, § 6 presents a summary of the main results and final remarks.

## 2 Dispersion, Diffusion and Turbulent Mixing in the ISM

The ISM is a turbulent medium with inhomogeneous makeup driven by SN explosions. The dominance of SN explosions was deduced by Cox & Smith (1974) and McKee & Ostriker (1977) from the detection of hot gas in UV absorption lines and X-ray emission. That these explosions interact to drive a turbulent flow was demonstrated numerically by Rosen & Bregman (1995). The nature of that flow in three dimensions has been explored with high-resolution hydrodynamical simulations by de Avillez (2000) and de Avillez & Berry (2001), while Korpi et al. (1999a, 1999b) included the effect of magnetic fields, at lower resolution. The SN explosions individually and collectively transport hot gas over large distances in more or less laminar flows. Turbulence driven by the interactions of these flows forms eddies that mix the gas on smaller scales. On the smallest scales molecular diffusion finally completes the mixing process.

### 2.1 Molecular Diffusion

Molecular diffusion is the consequence of frequent, stochastically distributed collisions between gas particles. The individual particles scatter in a random walk around their initial location. The total squared displacement of a particle

$$<\mathrm{\Delta }x>^2=2Dt$$
(1)
where $`D`$ is the diffusion coefficient given by

$$D=\lambda u_{\mathrm{rms}},$$
(2)
$`u_{\mathrm{rms}}`$ is the root mean square velocity of the particles, and $`\lambda `$ is the mean free path of the particles. The distribution of the particles over time can be described by a Gaussian distribution with a standard deviation $`\sigma =\sqrt{u_{\mathrm{rms}}\lambda t}`$ around the initial location of the particles.

If collisions happen in a homogeneous gas enclosed in a fixed volume, the relevant quantity to describe the diffusive process is the mean free path. It does not make sense to talk about a diffusion coefficient or an expected displacement because, viewed from the outside, the collisions do not change the properties of the gas, only the positions of the individual molecules.

The situation is different if there is a composition gradient. Then there are more particles of the species under study in one part of the volume than in the other. A random walk carries more particles out from the higher density region than in from the lower density region. Net transport results, reducing the gradient and eventually leading to a homogenized distribution. The streaming $`𝐒`$ of particles can be described as

$$𝐒=\stackrel{}{𝐃}N$$
(3)
with $`\stackrel{}{𝐃}`$ being the diffusion tensor for anisotropic diffusion and $`N`$ the number of particles per unit volume. The gradient is the driving force for the flow, so a larger gradient leads to a larger flow. The flow also depends on the mobility of the particles, described by the diffusion tensor. Since the diffusion tensor depends on the particle speed and on the mean free path, for a given gradient, the streaming $`𝐒`$, as well as the average displacement $`<\mathrm{\Delta }x>`$, are largest for fast particles undergoing few collisions (thus having a large mean free path) and smallest for slow particles undergoing many collisions.

If the diffusion is isotropic, the diffusion tensor reduces to a diffusion coefficient and the streaming becomes $`S=DN`$. From the conservation of mass one derives the diffusion equation

$$\frac{N}{t}=(DN)$$
(4)
If the diffusion coefficient is also independent of the spatial coordinate the equation can be reduced to

$$\frac{N}{t}=D^2N.$$
(5)
This equation shows that the characteristic diffusion time over a typical length $`L`$, is

$$\tau =\frac{L^2}{D}=\frac{L^2}{\lambda u_{\mathrm{rms}}}.$$
(6)
The time evolution of the particle density being transported, say along the $`x`$direction, is given by

$$N(x,t)=\frac{1}{2\sqrt{\pi Dt}}_{\mathrm{}}^+\mathrm{}N_o(x^{})e^{(xx^{})^2/4Dt}𝑑x^{}$$
(7)
where $`N_o(x)`$ is the initial particle distribution.

Let us now consider diffusion of an idealized checkerboard pattern with square size $`\mathrm{\Delta }L/2`$ subject to diffusion with an isotropic, uniform diffusion coefficient $`D`$. In one dimension, the pattern can be described as a square wave in composition $`C`$ along the $`x`$-axis with wavelength $`\mathrm{\Delta }L`$ (see Figure 1), so that at $`t=0`$,

$$C=C_o(x_k^{})=\{\begin{array}{ccc}0\hfill & if& L_{k1}x_k^{}<B_{k1}^+\hfill \\ 1\hfill & if& B_{k1}^+x_k^{}<B_k^{}\hfill \\ 0\hfill & if& B_k^{}x_k^{}<L_k,\hfill \end{array}$$
(8)
where $`B_{k1}^+=L_{k1}+\mathrm{\Delta }L/4`$, and $`B_k^{}=L_k\mathrm{\Delta }L/4`$, The square wave evolves following the solution given in equation (7)

$$C(x,t)=\frac{1}{2\sqrt{\pi Dt}}\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}I_k$$
(9)
where the solution must satisfy the periodic boundary condition $`C(L_{k1},t)=C(L_k,t)`$, and the integral $`I_k`$ has the final form

$$\begin{array}{c}I_k=_{L_{k1}}^{L_k}C_o(x_k^{}))\mathrm{exp}(\frac{(xx_k^{})^2}{4Dt})dx_k^{}=\hfill \\ \\ \sqrt{\pi Dt}\left[\text{erf}\left(A\left(xB_{k1}^+\right)\right)\text{erf}\left(A\left(xB_k^{}\right)\right)\right],\hfill \end{array}$$
(10)
where $`\mathrm{\Delta }L=L_kL_{k1}`$, $`\text{erf}(y)`$ is the error function,

$$\text{erf}(y)=\frac{2}{\sqrt{\pi }}_0^ye^{\eta ^2}𝑑\eta ,$$
(11)
and if we use a scale of parsecs for $`x`$ and years for $`t`$, the value of the constant $`A=8.69/\sqrt{D_{27}t_{\text{yr}}}`$, with $`D_{27}=D/10^{27}`$ cm2 s-1. Figure 1 shows the time evolution of the solution (9) for a medium having inhomogeneities with a scale length of $`\mathrm{\Delta }L/\mathrm{\Delta }L_o`$, where $`\mathrm{\Delta }L_o`$ is a scale length of reference. The solution is given in units of $`x/\sqrt{D_{27}}`$, with $`x`$ again taken to be in parsecs.

The mixing time scales vary as $`(\mathrm{\Delta }L/\mathrm{\Delta }L_o)^2t_o`$, where $`t_o`$ is the time of evolution of a reference square wave with length $`\mathrm{\Delta }L_o/2=100`$ pc. This result is a direct conclusion from the diffusion timescale. As an example, in a medium with $`D=10^{27}`$ cm2 s-1 subject to mixing according to a diffusive mixing length theory, inhomogeneities with length scale of 100 pc would be erased in a time of $`10^6`$ yr, while those with a length scale of 1 kpc would be erased in a time of $`10^8`$ yr.

### 2.2 Turbulent Mixing

Turbulent diffusion occurs when gas is stirred by more or less random forces, forming eddies within which local composition gradients increase, enhancing local molecular diffusion.

Studies of isotropic, homogeneous, incompressible turbulence led to the classical Kolmogorov (1941) cascade picture of turbulent flows. Eddies with a certain size $`\mathrm{}`$ are assumed to have some typical velocity $`v`$ associated with them. The corresponding Reynolds number $`Re=\mathrm{}v/\nu `$ will be larger for larger eddies, so viscosity $`\nu `$ is not very important for them. The energy cascades from these large eddies to smaller ones. The size of the smallest eddies is set by the condition $`Re=1`$. Thus, the minimum size of an eddy is

$$\mathrm{}_d\frac{\nu }{v_d},$$
(12)
and the energy in these eddies is lost to viscous dissipation. Energy is fed at some rate $`ϵ`$ per unit mass per unit time to the largest eddies of size $`\mathrm{}_o`$ and velocity $`v_o`$ with respect to neighboring eddies, for which the Reynolds number $`Re=\mathrm{}_ov_o/\nu 1`$. This energy then cascades to smaller and smaller eddies until it reaches eddies with the smallest size possible, $`\mathrm{}_d`$, which dissipate the energy at the same rate as they are fed in order to maintain equilibrium. Energy does not accumulate at any scale; the intermediate eddies merely transmit this energy to the smaller eddies in a cascade. These intermediate eddies are characterized by their size $`\mathrm{}`$ and velocity $`v`$. Energy passes through the cascade at a constant rate $`ϵv^3/\mathrm{}`$, independent of wavenumber (Kolmogorov 1941). The velocity $`v`$ of the eddies scales as

$$v=ϵ^{1/3}k^{1/3},$$
(13)
resulting in an eddy lifetime

$$\tau _{eddy}=\frac{1}{kv}=ϵ^{1/3}k^{2/3},$$
(14)
where $`k1/\mathrm{}`$ is the wavenumber. Thus, the time taken by an eddy to disappear varies as $`\mathrm{}^{2/3}`$. An eddy of scale length $`\mathrm{}_{100}`$, in units of 100 pc, exists for a time

$$\tau _{100}=(97.8\text{ Myr})\frac{\mathrm{}_{100}}{v}$$
(15)
where $`v`$ is written in units km s-1. We can estimate the time scale for mixing of hot gas using this classical eddy turnover time. Taking $`v`$ as the root mean square velocity of the hot ($`u_{\mathrm{rms}}100`$ km s-1) gas and assuming that the eddies have a length scale $`\mathrm{}100`$ pc, the time scale for complete mixing to occur is roughly 1 Myr. Therefore, if incompressible turbulent mixing dominates the mixing process, hot regions will mix very quickly.

### 2.3 Dispersal

Dispersal describes the transport of gas across large regions by laminar flows, such as those within SN remnants or superbubbles. On a large enough scale, encompassing many regions of laminar flow, dispersion can still be treated as a diffusive process that satisfies the diffusion equation (5) with $`N`$ replaced by the average density of the ensemble of transported particles $`<N(𝐫,t)>`$. In this case it does not make sense to speak about the different phases of the transported gas, because all of it is dispersed with the flow, and therefore, the mean free path $`\lambda `$, of the molecular diffusion, is replaced by a mixing length path $`L`$, which is defined as the distance over which the ensemble of particles is transported before it halts. Note that it no longer makes sense to define the rms velocity of each gas component in the transported flow separately, because the whole ensemble is transported roughly at the same speed.

According to this approximation, the diffusion timescale of a system with parsec-scale motion can be written as

$$\tau _{\mathrm{dif}}=(301\text{ yr})\frac{L_{\mathrm{pc}}^2}{D_{27}},$$
(16)
where $`L_{\mathrm{pc}}`$ is the length over which the ensemble of particles is transported, in units of parsecs; $`D_{27}=lu_{\mathrm{rms}}`$ is the diffusion coefficient, written in units of $`10^{27}`$ cm2 s-1; and $`u_{\mathrm{rms}}`$ is the root mean square velocity of the emsemble of particles being transported. If one assumes that the characteristic length scale over which the flow is transported is the same as its mixing length, then the diffusion time scales as $`\tau _{\mathrm{dif}}L/u_{\mathrm{rms}}`$

## 3 Simulations

### 3.1 Supernova-Driven ISM Model

To study the mixing timescales in the ISM we added tracer fields to a modified version of the three-dimensional, SN-driven, ISM model of Avillez (2000). The model is run with an adaptive mesh refinement (AMR) hydrodynamics code on a $`1\times 1\times 20`$ kpc region of the Galactic disk. It includes a fixed gravitational field provided by the stars in the disk, radiative cooling assuming optically thin gas in collisional ionization equilibrium, and uniform heating due to starlight. The radiative cooling function is a tabulated version of that shown in Figure 2 of Dalgarno & McCray (1972) with an ionization fraction of 0.1 at temperatures below $`10^4`$ K and a temperature cutoff at 10 K. Background heating due to starlight varies with $`z`$ as described in Wolfire et al. (1995) and is kept constant in the directions parallel to the plane; in the Galactic plane at $`z=0`$ it is chosen to initially balance radiative cooling at 8000 K. The presence of background heating leads to the creation of thermally stable phases in the ISM. Due to the presence of a stable phase at low temperatures, the amount of cold gas seen in these simulations is larger than that found in the simulations presented in Avillez (2000).

The interstellar gas is initially distributed in a smooth disk with the vertical distribution of the cool and warm neutral gas given by Lockman et al. (1986) and summarized in the Dickey & Lockman (1990) distribution. In addition, an exponential component representing the $`z`$distribution of the warm ionized gas with a scale-height of 1 kpc in the Galaxy is included, as described in Reynolds (1987). The model includes supernovae types Ib+c and II. Supernovae of type Ia are not included in these simulations because their role is not important in the mixing process that occurs near the Galactic plane, although they may be important in the disk-halo interaction due to their scale-height.

The rates of occurrence of SNe types Ib and Ic in the Galaxy are $`2\times 10^3`$ yr-1, while those of type II occur at $`1.2\times 10^2`$ yr-1 (Cappellaro et al. 1997). We take the total rate of these SNe in the Galaxy to be $`\sigma _{\mathrm{gal}}=1.4\times 10^2`$ yr-1, corresponding to a rate of one SN every 71 yr. These rates are normalized to the volume of the stellar disk in which a specific SN type is found: a galactic radius of 12 kpc is assumed for all stellar disks, but their half thicknesses are assumed to be twice the scale height of the corresponding SNe distributed in the field. In particular as SNe Ib+c and II are assumed to have the same scale height then their stellar disks have the volume $`1.6\times 10^{11}`$ pc3.

SNe of types Ib+c, and II are explicitly set up at the normalized SN rate, with 40% of them placed at random locations distributed in an exponential distribution with a scale height of 90 pc, and 60% in locations where previous SNe occurred, representing OB associations, in a layer with a scale height of 46 pc following the distribution of the molecular gas in the Galaxy.

The first SNe in associations are chosen to occur in locations where the current local density is greater than 1 cm-3. No density threshold is used to determine the location where isolated SNe should occur, because their progenitors drift away from the parental association and therefore, their site of explosion is not correlated with the local density. Similarly, later SNe in associations are no longer determined by gas density.

The SNe are set up at the beginning of their Sedov phases, with radii determined by their progenitor masses, which are injected into the location of the explosion. Type II SNe come from early B stars with masses $`7.7M_{}M15M_{}`$ while type Ib+c SNe have progenitors with masses $`M15M_{}`$ (Tammann et al. 1994). In this model the maximum mass allowed for an O star is 30 $`M_{}`$. Avillez (2000) describes in detail the algorithm used to set up the isolated and clustered SNe during the simulations.

### 3.2 Numerical Method and Boundary Conditions

These simulations use the piecewise-parabolic method of Colella & Woodward (1984), a third-order scheme based on a Godunov method implemented in a dimensionally-split manner (Strange 1968) that relies on solutions of the Riemann problem in each zone rather than on artificial viscosity to follow shocks. The Godunov approach is typical of standard numerical techniques in regions where the solution is smooth. However, in regions with discontinuities, such as strong shocks, the Godunov method approximates the solution well by analytically solving an associated Riemann problem. This is an idealized problem describing the evolution of a simple pressure jump into shocks and rarefactions, with a contact discontinuity in between. Monotonicity constraints ensure that these discontinuities remain sharp and accurate as they traverse the computational grid. The higher-order spatial interpolation in the PPM allows steeper representation of discontinuities.

During the simulation, the mesh is refined periodically in regions with sharp pressure variations using the AMR scheme. The local increase of the number of cells corresponds to an increase in linear resolution by a factor of two (that is, every refined cell is divided into eight new cells). At every new grid the procedure outlined above is carried out, followed by the correction of fluxes between the refined and coarse grid cells. The adaptive mesh refinement scheme is based on Berger & Colella (1989), but the grid generation procedure follows that described in Bell et al. (1994).

The computational domain has an area of 1 kpc2 and a vertical extension of 10 kpc on either side of the midplane. In the simulations discussed here, AMR is used in the layer $`\left|z\right|500`$ pc. In the highest resolution runs, three levels of refinement are used, yielding a finest resolution of 1.25 pc. For $`\left|z\right|>500`$ pc the resolution is 10 pc. Periodic boundary conditions are used on the side boundaries, while outflow boundary conditions are used on the top and bottom boundaries.

### 3.3 Tracer Field Implementation and Mixing Runs

To follow how composition differences mix, we use a tracer field. A tracer field acts as a drop of ink in a fluid. In time the drop will spread and eventually mix completely into the fluid. This occurs because thermal motion leads to collisions between ink and fluid molecules, distributing both species uniformly. The simulations do not include any physical diffusion term, however, so contrary to what happens in a real fluid or gas, the mixing of the tracer field occurs by numerical diffusion, which will generally be faster and larger scale than the physical diffusion in astrophysical problems. As a consequence the mixing time of a tracer field in our model provides a rather strong lower limit for the timescale of mixing resulting from physical diffusion.

The tracer field is an intensive scalar whose advection is computed with the same algorithm as the density. We initially set its value to be either one or zero to represent regions of different composition. In the absence of diffusion, either numerical or physical, these values would remain fixed as the gas is advected. In the presence of diffusion, the tracer field values will change. In our simulations, this is a result of the numerical diffusion present in the Godunov advection scheme.

This scheme assumes that the flow solution is represented by a series of piecewise constant states, and therefore, that the numerical representation approximates the true solution near discontinuities. The numerical solution is evolved by considering the nonlinear interaction between these piecewise constant states. Viewed in isolation, each pair of neighbouring states constitutes a Riemann problem. The algebraic solution in the overall grid results from a collection of Riemann problems for all the interfaces between two sucessive cells at every time step. The discontinuity between two states evolves into some combination of shock and rarefaction waves, separated by a contact discontinuity (Figure 3). The tracer field is advected in the same way as the density and therefore, it diffuses between any rarefaction wave and the contact discontinuity formed as a result of the Riemann problem solution. This process happens between every pair of neighboring cells, leading sooner or later to a change of the tracer field value. In the case at hand we use tracer fields to follow the evolution and mixing of composition differences in the ISM. As mixing occurs by numerical diffusion, regions with values intermediate between zero and unity grow, until ultimately the entire volume is filled with gas having value roughly 0.5.

The tracer field is set up initially with values of $`C=0`$ and 1 on alternating squares of a checkerboard in the plane of the Galaxy with square sides of $`l=25,50`$ or 500 pc (Figure 2). The 500 pc scale corresponds to a tracer field that fills a quarter of a 1 kpc length square board. The distribution is uniform in the vertical direction; that is, the squares are extended into rectangular solids of uniform composition vertically. This setup allows a direct comparison to the analytic solution given in § 2.1 for the mixing of inhomogeneities of different length scales in the ISM due to classical turbulent diffusion.

In the current work we report on simulations using seven SN rates $`\sigma /\sigma _{\mathrm{gal}}=1`$, 5, 10, 15, 20, 30 and 50, having finest grid resolutions of 1.25, 2.5 or 5 pc. In general, the simulations were stopped between 50 and 70 Myr after complete mixing occurred. The simulations ran for 400 Myr for $`\sigma /\sigma _{\mathrm{gal}}=1`$, 300 Myr for $`\sigma /\sigma _{\mathrm{gal}}=5`$, and 200 Myr for the remaining SN rates. A summary of these runs is presented in Table 1.

## 4 Results

### 4.1 Overview of Results

In the first 50 Myr the system evolves into a statistical steady state on the global scale. Once disrupted by SN explosions, the disk never returns to its initial state, provided SNe continue to explode. Instead, regardless of the initial vertical distribution of the disk gas, a thin disk of cold gas forms in the Galactic plane, and, above and below, a thick inhomogeneous gas disk forms. Gas flows between the thin and thick gas disks, with upward and downward flowing gas coming into dynamical equilibrium. The upper parts of the thick disk form the disk-halo interface, where a large scale fountain into the halo is driven by hot ionized gas escaping in a turbulent convective flow from the disk-halo interface and from superbubbles blowing out of outer layers of the thin disk.

Table 1 gives the time for complete mixing, measured in the supernova forming stellar disk, which we take to be the time needed for the average value of the tracer field to saturate at a value $`0.5`$ (see §4.2). Results are shown for checkerboards with initial square sizes of $`l=25`$, 50, and 500 pc, and SN rates $`\sigma `$/$`\sigma _{gal}=1`$, 5, 10, 15, 20, 30, and 50, at resolutions $`\mathrm{\Delta }x`$ of the finest grid levels used during adaptive mesh refinement. A quick look at the table shows that the time required for complete mixing of inhomogeneities with $`l=25`$ and $`50`$ pc differs by only 5–20 Myr, while that for inhomogeneities with $`l=500`$ pc is larger by, at most, another 20 Myr. These differences decrease with increasing SN rate. Runs with different finest resolution $`\mathrm{\Delta }x`$ for the same SN rate show similar mixing time scales, suggesting that the results depend quite weakly on the actual strength of the diffusion.

Figures 4 and 5 present snapshots of the evolution of tracer fields in checkerboards with initial length scale of $`l=50`$ pc and SN rates $`\sigma /\sigma _{\mathrm{gal}}=1`$, 5 and 15. The snapshots in Figure 4 were taken at 50 and 126.6 Myr, while those in Figure 5 were taken at 50 Myr of evolution. The resolution of these snapshots is 1.25 pc in Figure 4 and 2.5 pc in Figure 5. We see that for the Galactic SN rate, the mixing process at 50 Myr is just starting to occur, while for larger SN rates, it is well advanced, with a larger fraction of the inhomogeneities reduced to weak structures that take a long time to be erased. For example, in the model with $`\sigma /\sigma _{\mathrm{gal}}=15`$ at a time of 50 Myr, the tracer field has a value varying between 0.4 and 0.6. For the other rates shown in the figures at 50 Myr there are still unmixed regions where tracer field has values of zero or one.

These Figures show large-scale laminar flow in the hot gas bubbles that transports tracer field over large distances. The expansion and interactions of the bubble drive shocks into the warm and cold gas. Vorticity generated at shock intersections drives eddies mixing the gas up to 100 pc scales, as well as Rayleigh-Taylor fingers that degenerate into smaller eddies at the tips of their caps. The eddies cascade from intermediate to small scales, further mixing the tracer field down to the diffusion scale. The mixing process thus proceeds by different mechanisms at different scales, because of the inhomogeneous, compressible nature of the flow.

The mixing process is independent of the spatial distribution of isolated SNe that occur randomly in the field. The locations of type II and Ib+c SNe occurring away from their parental OB associations are determined from a file of random numbers chosen at the beginning of a simulation. We ran simulations with different sets of random numbers and finest grid resolutions of $`\mathrm{\Delta }x=1.25`$ and 2.5 pc, and found mixing time scales virtually the same as the ones reported here, so we do not further discuss those models.

### 4.2 Average, Variance, and Maximum

The simulations start with exactly the same number of cells having value 1 and 0. Therefore, the average value of the tracer field, $`<C>=0.5`$ at $`t=0`$. By definition complete mixing occurrs when the tracer field at any point in the grid has $`C_{max}=C_{min}=0.5`$. However, due to round-off and phase errors of the numerical scheme, the actual value for completely mixed material may differ from 0.5 by a few percent. Once the tracer field has this constant value on the entire grid, no further variation occurs.

Figure 6 shows the time variation of $`<C>`$ for three SN rates $`\sigma /\sigma _{\mathrm{gal}}=1,5`$ and 15, and for two finest grid resolutions of 2.5 and 5 pc. During the first 50 Myr, $`<C>`$ shows small variations around 0.5, approaching a constant value as mixing increases and the small scale structure dominates the mixing process. The initial variations of $`<C>`$ result from the increased error occurring as steep color gradients interact with steep velocity gradients in SN-driven shock waves. As the color is mixed, the color gradients decrease, and $`<C>`$ approaches a constant value close to 0.5. The biggest errors, still only 2.4%, occur for the highest SN rates, where strong shocks are more widely distributed.

The average value of $`C`$ gives an indication of the time at which complete mixing occurred, but it does not provide any information of the mixing process. Further information may be obtained from the variance of the tracer field, $`\sigma _C^2`$. This evolves from $`\sigma _C^2=0.25`$ at $`t=0`$ to $`\sigma _C^2=0`$ when complete mixing occurs. The variance has a steady decrease with time until it reaches a value of $`10^3`$. For $`\sigma _C^2<10^4`$ the decrease is no longer steady, but rather sporadic, suggesting that we have reached a regime where numerical noise is dominating. Figure 7 shows the variance of $`C`$ (for $`\sigma _C^210^4`$) over time for inhomogeneities with different length scales, along with exponential fits to the steady decrease of $`\sigma _C^2`$ until it reaches a value of $`10^3`$. The fits shown are for the shorter length scales of the inhomogeneities. We find that the variance of the tracer field for inhomogeneities with length scales $`l=25`$ and 50 pc decays to $`10^3`$ remarkably close to exponentially, with $`\sigma _C^2e^{t/\tau }`$. For $`\sigma _C^2<10^3`$ such decrease is seen for $`\sigma /\sigma _{Galaxy}=1`$ or 5 if $`l=25`$ pc and $`\sigma /\sigma _{Galaxy}10`$ if $`l=50`$ pc.

Figure 8 shows the variation of the time constant $`\tau `$ with SN rate for inhomogeneities with initial $`l=25`$ and 50 pc, as well as the best power-law fits to the data points. The fits to an inverse square root law are remarkably good, with correlation coefficients of 0.997 and 0.998 for $`l=25`$ and 50 pc, respectively. For length scale of 50 pc,

$$\tau =(37.9\text{ Myr})(\sigma /\sigma _{\mathrm{gal}})^{0.5}.$$
(17)
In future work we will extend our parameter study sufficiently to define the dependence on inhomogeneity length scale $`l`$ fully.

This exponential decay law indicates that large variations in the tracer field decay quickly while smaller variations last much longer. Such a decay law for the variance of $`C`$ breaks down for larger inhomogeneities, and in particular it does not apply for inhomogeneities with $`l=500`$ pc. The extended parameter study should show at which length scale the deviation from this behavior begins to occur. We will also examine whether this result extends to smaller scales using finer grids.

The time variation of the maximum value of the tracer field, $`C_{max}`$ is shown in Figure 9. The exponential decay behavior is again seen, with a fast initial decline followed by a slow approach to complete mixing. Structures with small variations in composition will thus be present for extended periods of time. The figure also shows that the time taken for complete mixing to occur decreases with increasing SN rate until it saturates for $`\sigma /\sigma _{\mathrm{gal}}>20`$. This result is also seen in the variation of the mixing time with SN rate shown in Figure 10, which shows that for SN rates $`\sigma /\sigma _{\mathrm{gal}}>20`$ the mixing time scales differ by only a few Myr.

The mixing process has a rather weak dependence on the resolution used in the calculations, and thus of the actual strength and scale of diffusion. Running the scheme with a better resolution decreases the numerical errors, and therefore also the magnitude of numerical diffusion. Therefore, the mixing times become longer (going from right to left column in Table 2), but this increase is far slower than the increase in resolution, and declines with increasing resolution, as shown in Table 2 and Figures 11 and 10. This demonstrates that numerical convergence is occurring, and further that unresolved details of the mixing process at small scales do not affect the behavior at resolved scales. This appears to be because the mixing process depends primarily on motions near the SN driving scale of $`l100`$ pc, with the presence of small scale diffusion more important than its exact properties.

The similarity in the mixing time scales for $`\sigma /\sigma _{\mathrm{gal}}>20`$ occurs because for SN rates higher than this the entire disk quickly heats. Figure 12 shows the temperature probability distribution functions (PDFs) of the disk gas for $`\sigma /\sigma _{\mathrm{gal}}=20`$. After the first 50 Myr of evolution a large fraction of the disk gas already has temperature $`T>10^5`$ K and by 100 Myr of evolution all the gas has reached at least $`10^5`$ K. As a consequence of the high SN rate the disk becomes warmer and cooling becomes inefficient. The high rms velocities in this hot gas lead to quick mixing, as in hot regions in less energetic models. The complete lack of colder regions, however, means that all the gas mixes on the same short time scale. Once there is no further cold gas, increasing the SN rate further only slowly increases the rms velocity, leading to the saturation behavior seen.

We emphasize that classical mixing length theory badly fails to predict the mixing time scales that we compute. That theory predicts that the mixing time scales should depend on the length scales. That is, let a checkerboard be composed of squares with length scale $`L`$, and the squares be further divided into $`n`$ smaller squares with length scale $`l`$. Then mixing length theory would predict the diffusion time scales of these squares to be related by $`\tau _L/\tau _l=n^2`$ for constant diffusion coefficient. Thus, the predicted ratios of timescales (for example for $`L`$ and $`l`$ in parsecs) are $`\tau _{500}/\tau _{50}=100`$ and $`\tau _{50}/\tau _{25}=4`$. However, Table 1 shows that the actual ratioes are never greater than 1.2, a value far smaller than predicted.

### 4.3 Profiles

Another approach to measuring the mixing timescale is to directly compute the diffusion coefficient required to fit the average profile across an initially sharp color gradient, using the analytic theory of § 2.1. Figure 13 presents profiles of the tracer field averaged along the $`x`$-axis for $`y500`$ pc for the checkerboard with squares of initial length of 500 pc (shown in the top panel of Figure 2), overlaid by the classical diffusion solution from equation (10) with $`n=1`$. Initially the solutions fit the average profile of the tracer field well, but at later times the fits become somewhat worse. The correlation factors of these fits reduces as time evolves from 0.998 at 25 Myr to 0.94 at 225 Myr.

From the fits, we can deduce the diffusion coefficient $`D_{27}`$ written in units of $`10^{27}`$ cm2 s-1. The time variation of D27 for two different SN rates is shown in Figure 14. The straight lines are the best fit to the data points. We find that the diffusion coefficient varies exponentially with time, with the general expression of

$$\mathrm{log}D_{27}=A+t/\tau _D,$$
(18)
where, if $`t`$ is given in Myr, $`A=2.91`$ and $`\tau _D=107`$ Myr for $`\sigma /\sigma _{\mathrm{gal}}=1`$, and $`A=2.45`$ and $`\tau _D=67.1`$ Myr for $`\sigma /\sigma _{\mathrm{gal}}=5`$. If the mixing were regulated by classical diffusion, then we would expect that the diffusion coefficient would be nearly constant or vary around a mean value with time, rather than exponentially increasing.

## 5 Discussion

In this paper we try to understand how long it takes to erase inhomogeneities in the ISM and if this process depends on the inhomogeneity length scale. That is, does it take longer to mix big regions than small ones? We modified the three-dimensional model of Avillez (2000) to incorporate tracer fields and ran simulations with different SN rates, distributions of SNe in the field, and numerical resolutions. (Clustered SNe from OB associations are explicitly set up, with the first SN in an association occurring in a location where the local density is greater than 1 cm-3.)

The three-dimensional adaptive mesh refinement code used in our simulations, unlike 1-zone chemical evolution models (Matteucci et al. 1999; Tosi et al. 1998) and chemodynamical models (Samland et al. 1997; Samland 1998; Argast et al. 2000, to name a few), enables us to resolve local inhomogeneities in the ISM with a spatial resolution as good as 1.25 pc. This allows us to test the usual assumption of Galactic chemical evolution models that assume instantaneous mixing of the ISM at all times. Mixing in the real world takes finite time, and ultimately occurs on the molecular scale, through physical diffusion. We set a lower limit on the mixing time scale by assuming that it is set at the far larger grid scale by faster numerical diffusion. We find that the actual value of the mixing scale does not strongly influence the mixing timescale.

Matter can enter and leave the section of the Galaxy under study through the top and bottom outflow boundary conditions, as well as moving to neighboring periodic boxes (eg passing through a periodic side boundary condition), although we do not yet include Galactic rotation in our model. Chemical evolution models do not normally allow radial transport of matter between different rings, but it occurs in our models.

Because the domain extends to 10 kpc on both sides of the midplane, gas can cycle between the disk and the halo in a Galactic fountain as a result of the dynamics of the flow in the gravitational field of the Galaxy (see Kahn 1981; Avillez 1999). The mixing time scales of the disk gas are regulated mainly by the SNe, because the disk gas that enters the disk-halo cycle takes some 200 Myr to return to the disk (Kahn 1981). Furthermore before reaching the thin gas disk, the descending gas hits the thick gas disk and decelerates due to drag (Kahn 1981; Benjamin & Danly 1997). Therefore, most of the descending gas will not contribute strongly to the mixing process of the ISM with $`\left|z\right|250`$ pc.

This work only studies the mixing process as a way of erasing inhomogeneities in the ISM. It does not take into account continuing pollution of the ISM by SNe, planetary nebulae, and stellar winds, which we will study in future work. The actual enrichment of the disk gas depends on the sources and the efficiency of mixing. If the ejected metals are quickly distributed over a large volume, spatially homogeneous enrichment takes place. If the mixing volume is small, the ISM in the vicinity of a SN will be highly enriched, while large parts of the disk remain metal-poor until the gas mixes out on time scales of hundreds of Myr or more. In this case the ISM may remain chemically inhomogeneous, with newly formed stars having different compositions depending on where they form. Over gigayear timescales, mixing does appear to be quite effective over scales of many kiloparsecs.

This work shows that inhomogeneities that are present in the ISM, assuming that no further inhomogeneities are introduced into the system, can take a long time to be erased. If inhomogeneities are introduced at a time larger than the time scale for mixing, the system will show temporarily an homogeneous state, which ends when inhomogeneities are introduced. On the other hand, if inhomogeneities are introduced on a time scale smaller than the mixing time, then the system will show an inhomogeneous distribution during its evolution. Therefore, from the simulations we conclude that for the ISM to show a homogeneous distribution, no inhomogeneities should be introduced over periods of roughly a hundred megayears. This is longer than typical star formation time scales, so new SNe and pre-main sequence stars can continue to maintain an inhomogeneous composition so long as star formation continues.

The rate at which inhomogeneities are erased increases with the SN rate to values about ten times greater than the Galactic rate. Above that rate, no cold medium remains, mixing occurs quite quickly, but above this threshold, further increasing the SN rate does not markedly further increase the mixing efficiency.

The time for complete mixing is already longer than any time scale for chemical equilibrium to occur and therefore the ISM will show an uneven distribution. Even if the major fraction of the ISM gas is completely mixed, lower level inhomogeneities will still be seen as abundances are measured more carefully. Furthermore, as the mixing time scales are longer than the time interval for SNe type II to occur and enrich the ISM with metals, low-level inhomogeneities are continuously maintained, leading to a poorly mixed ISM, except in the hot regions where the mixing process is fast (of order 50 Myr).

For a galactic SN rate, the ISM is likely fairly inhomogeneous, since the SNe are spatially well separated and erasing these inhomogeneities takes a long time (of the order of 350 Myr). A single SN event can thus substantially influence its surroundings. If we consider, for example, a 20 pc sphere with number density $`n=1`$ around a SN with 1 M of metals in its ejecta, the ejecta represent a fraction $`0.1(Z/Z_{})`$ of the metals in the region. Thus, individual SNe can drastically change the local metallicity in low metallicity gas. Stars formed before mixing has finished can have a wide range of metallicities; further work will be needed to make quantitative predictions of the expected dispersion for comparison with observations, however.

The smallest mixing time scale that we find in our study is $`120`$ Myr and occurs for $`\sigma /\sigma _{\mathrm{gal}}20`$. After some $`80\%`$ of the gas has been mixed, the simulations show that the ISM remains inhomogeneous for at least twice as long on scales smaller than a kpc. Even when the rate of SNe is increased to ten times the Galactic rate, this uneven abundance distribution structure takes tens of Myr to disappear.

At the finest grid scale we suppress the turbulent diffusivity, as numerical viscosity becomes the dominant diffusion mechanism. At these scales there would physically be a continuing turbulent mixing that we do not resolve. The resolution study shown in Figures 10 and 11 demonstrates, however, that the mixing time is practically independent of the actual value of the diffusion scale. Changing the linear diffusion scale by factors of as much as four produced less than 20% changes in the diffusion time. This means that we can effectively assume that the mixing is dominated by an inertial scale decoupled from the diffusion scale, so that the details of how the diffusion occurs do not affect the larger-scale behavior.

## 6 Summary and Final Remarks

The main results of this study are:

* The mixing of a SN-driven ISM is characterized by three regimes: laminar flows dominate on scales close to a kpc, turbulent mixing dominates below about a hundred parsecs, and diffusive processes finally dominate at the smallest scale.
* Mixing at the larger scales appears nearly independent of the exact value of the diffusive scale, so long as some diffusive process ultimately operates at small scales.
* The time scale to erase inhomogeneities in the ISM on length scales from parsecs to kiloparsecs is nearly independent of the length scale; classical mixing length theories badly fail to describe this behavior.
* The time scale for complete mixing, providing no further inhomogeneities are introduced, is some 350 Myr for the Galactic SN rate, decreasing with increasing SN rate up to a threshold value. Above this threshold value, further increases in the SN rate have little further effect.
* The speed of mixing decreases almost exponentially, with a time constant dependent on the inverse square root of the SN rate: strong inhomogeities decay quickly, while weaker ones last much longer.
* The diffusion coefficient derived from fits of a solution of the diffusion equation to the simulation results does not remain constant as assumed by classical mixing length theory, but rather increases exponentially with time, with a time constant dependent on the SN rate.

This work is far from complete. At least four issues clearly remain to be addressed. First, we need to determine the actual length scale at which the exponential decrease in tracer field variance ceases, somewhere between 50 and 500 pc. Second, if further inhomogeneities are introduced, such as the chemical elements fed into the ISM by SNe, planetary nebulae, and stellar winds, large parts of the ISM, or indeed the whole ISM, may never reach homogeneity. Even though stronger inhomogeneities decay quickly, the exponential decline in mixing rate for weaker inhomogeneities means that a wealth of weaker inhomogeneities will remain. Determining quantitative values of the dispersion of metallicities will help solve outstanding questions in studies of stellar populations and the ISM. Third, the simulations do not include the Galactic magnetic field. One can speculate that the presence of magnetic field would contribute to an increase of the mixing time scales in the Galaxy. It could reduce or suppress turbulent flows, and also it would reduce the expansion of laminar flows. Furthermore, the random component of the field could enhance the diffusion in the SN shells, thereby fragmenting the shells more rapidly and affecting the mixing process. We plan to run a new set of simulations to determine dependence of the mixing process on field intensity and distribution. Finally, a theoretical approach, supported by direct numerical simulations, must be developed in order to explain the mixing process in the ISM, because as this work shows, the classical approach to the problem fails.

The authors thank D. York, E. Jenkins, and R. Ferlet for emphasizing the importance of this problem and them and M. Saha for useful discussions. This work was supported by an NSF CAREER grant (AST 99-85392), and has made use of the NASA Astrophysical Data System.