text:
# Role of interactions in the far-infrared spectrum of a lateral quantum dot molecule

## Abstract

We study the effects of electron-electron correlations and confinement potential on the far-infrared spectrum of a lateral two-electron quantum dot molecule by exact diagonalization. The calculated spectra directly reflect the lowered symmetry of the external confinement potential. Surprisingly, we find interactions to drive the spectrum towards that of a high-symmetry parabolic quantum dot. We conclude that far-infrared spectroscopy is suitable for probing effective confinement of the electrons in a quantum dot system, even if interaction effects cannot be resolved in a direct fashion.

Nanoscale semiconductor structures are very promising for future components of microelectronic devices. These systems are also scientifically very interesting as they exhibit novel and fundamental quantum effects. The most prominent difference between the two-dimensional artificial atoms, or quantum dots (QD), and their normal counterparts are the enhanced correlation and magnetic field effects. In addition to that, the new features can be controlled by the tunable system parameters, both experimentally and in theoretical models. Coupling together QDs, one can construct QD molecules (QDM). The electronic structure of these is also quite intriguing. For example, the two-electron QDM has a highly non-trivial spin-phase diagram and composite-particle structure of the wave function Ari . One experimental possibility to probe the new physics of nanostructures is use far-infrared (FIR) absorption spectroscopy HeitmannKotthaus ; JacakQD . For highly symmetric QDs, the FIR spectra reflect only the center-of-mass properties of the system. If the symmetry is lowered, the relative motion of the electrons starts to couple to the center-of-mass motion, enabling one to study the many-body correlation effects using FIR.

In this Letter we present FIR spectra for a two-electron QDM using exact diagonalization. The method used solves accurately the quantum-mechanical model for a QDM, revealing the important electron-electron correlation effects. Our calculated spectra show clear deviations from the spectra of highly symmetric QDs. To analyze these deviations, we vary our system parameters in order to separate the correlation effects from the single-particle effects of the low-symmetry confinement potential. This analysis shows, however, that these effects are ultimately entangled. Furthermore, we find, very surprisingly, that the correlation effects actually compensate the lowered symmetry of the QD confinement potential, resulting in spectra much closer to an electron in a high-symmetry QD than in a low-symmetry one.

The starting point for understanding FIR spectra is the highly symmetric, parabolic QD. This spectrum consists of two Kohn modes, whose dispersion does not depend on the number of confined electrons nor their interactions MaksymChakra . Several experiments on QDs have shown deviations of the FIR spectra from the Kohn modes HeitmannKotthaus ; Demel ; Dahl ; Dahl2 ; Hochgrafe . It is clear that these deviations require a non-parabolic QD, but the detailed cause of the deviations, and thus the interpretation of the measured spectra, is difficult to obtain. An especially interesting question is how the electron-electron interaction and correlation effects emerge in the FIR spectrum when the symmetry of the QD is lowered. In our model, we can tune both the deviation of the potential from being perfectly parabolic as well as the electron-electron interaction strength, which enables us to analyze the source of the non-trivial details seen in the FIR spectra. This type of analysis is extremely important for explaining the details of the spectra.

The calculated FIR spectra have shown that a broken $`xy`$-symmetry (e.g. elliptic or rectangular) is required to observe the zero-field splitting of the Kohn modes as well as additional modes, and a broken rotational symmetry (e.g. square) is required to observe anticrossings Veikko ; Pfannkuche ; Madhav ; Magnusdottir ; Ullrich ; Rodriguez . However, the role of interaction effects and lowered QD symmetry in FIR spectra is not fully analyzed in these studies. Our QDM model suits perfectly for this detailed analysis as the FIR spectra show all the above-mentioned features. One interesting finding is that the electron-electron interactions drive the FIR spectrum closer to the Kohn modes instead of making the deviations more pronounced. This means that the electrons feel an effective potential that is more parabolic than the the bare confinement one. Thus the interactions smoothen the non-parabolicity of the external potential.

We model the two-electron QDM by a 2D Hamiltonian

$$H=\underset{i=1}{\overset{2}{}}\left(\frac{(i\mathrm{}_i\frac{e}{c}𝐀)^2}{2m^{}}+V_\mathrm{c}(𝐫_i)\right)+C\frac{e^2}{ϵr_{12}},$$
(1)
where $`V_\mathrm{c}`$ is the external confining potential, for which we choose $`\frac{1}{2}m^{}\omega _0^2\mathrm{min}\{(xd/2)^2+y^2,(x+d/2)^2+y^2\}`$. This potential separates to two QDs at large inter-dot distances $`d`$, and with $`d=0`$ it simplifies to one parabolic QD. We use the GaAs material parameters $`m^{}/m_e=0.067`$ and $`ϵ=12.4`$, and the confinement strength $`\mathrm{}\omega _0=3.0`$ meV. $`𝐀`$ is the vector potential of the magnetic field (along the $`z`$ axis) taken in the symmetric gauge. We use $`C`$ to scale the Coulomb interaction strength. The wave functions and energies of the ground and excited states are found by the exact diagonalization of the Hamiltonian matrix as in Ref. \[Ari, \].

Using the Fermi golden rule within the electric-dipole approximation for the perturbing electromagnetic field, the transition probability from the ground state to the $`l`$’th excited state can be calculated as

$$𝒜_{l,\pm }\left|\mathrm{\Psi }_l\left|e^{\pm i\varphi }_{i=1}^2𝐫_i\right|\mathrm{\Psi }_0\right|^2\delta (E_lE_0\mathrm{}\omega ),$$
(2)
where we have used right- or left-handed circular polarization of the field, corresponding to the two possible signs marked with ’$`\pm `$’. We mainly present results for unpolarized light, averaging over the two circular polarizations. In addition to that, we decompose one selected spectrum to two circular as well as to two axial polarizations.

The calculated FIR spectra are shown in Fig. 1 for several different system parameters. We plot the transition energy as a function of the magnetic field using a line whose width is proportional to the transition probability. We have only included the transitions which have a probability of more than one percent of the maximum value. To make it easier to compare the relative strengths of the different resonances, we have plotted the probabilities also separately below each spectrum. In addition to that, we show the energies of the two Kohn modes of a parabolic QD using open circles. The upper row of figures present spectra for the total spin $`S=0`$ and the lower one for $`S=1`$. The spectra for both spin types are shown for two inter-dot distances $`d`$ and three Coulomb strengths $`C`$. The ground state of our system changes from $`S=0`$ to $`S=1`$ when the magnetic field is increased (see Ref. \[Ari, \] for more details), and this transition is shown by a vertical line in the four spectra of the fully interacting cases $`C=1`$. One should note that the transition probabilities between different spins are zero, and thus $`S`$ is conserved in the allowed transitions.

As a general feature of the calculated spectra shown in Fig. 1 one can see that each spectrum has as a major component two branches, where the higher one ($`\omega _+`$) has a positive dispersion and the lower ($`\omega _{}`$) a negative one. For the parabolic QD, these branches are the two Kohn modes (circles in the figures). Now the lowered symmetry of the QD and the correlation effects show up in the FIR spectra as a deviation from the Kohn modes. These deviations are the splitting of $`\omega _+`$ and $`\omega _{}`$ at $`B=0`$, the anticrossings in $`\omega _+`$, and an additional mode $`\omega _{+2}`$ above $`\omega _+`$ in the $`S=1`$ spectra. Next we will analyze these deviations in detail.

Figs. 1 (a) and (b) show the ($`C=1`$) spectra for the distance $`d=10`$ nm which corresponds to closely coupled QDs. One can see that there are only minor deviations from the Kohn modes. For $`S=0`$, one can see a small $`B=0`$ splitting of the two modes. There are also anticrossings in the upper modes, but these are hardly visible as the energy gap is much smaller than the line width. These anticrossings are, however, visible as oscillations in the transition probabilities shown below the spectra. These probabilities also show that the upper mode of the $`S=1`$ spectra actually consists of two levels. These levels are energetically nearly degenerate, but the probabilities vary as can be seen in lower part of Fig. 1 (b). The sum of the two upper modes adds up to a probability which is nearly constant and close to the probability of the lower mode.

If we increase the distance $`d`$ between QDs, the confinement potential becomes less parabolic. In the extreme limit, however, the two QDs decouple and one is left with two electrons localized in separated dots. Figs. 1 (c) and (d) show the spectra for $`d=30`$ nm which is in the interesting range of coupling between the two QDs. One can see that the features discussed for $`d=10`$ nm are again present, but at this time much clearer in the spectra. The $`B=0`$ splitting is rather large for $`S=0`$, and there is also a small splitting in the $`S=1`$ case. There are many more anticrossings than in the $`d=10`$ nm case, and those for weak $`B`$ now have a clear energy gap. The two nearly degenerate modes of the $`S=1`$ case are now split in energy, and the higher one ($`\omega _{+2}`$) has a smaller transition probability than $`\omega _+`$. The comparison of the spectra of Figs. 1 (a) to (c) and (b) to (d) shows that the first anticrossing point is rather insensitive to the change of $`d`$ even if the gap in energy opens up clearly. A further comparison shows that at $`B=0`$, the $`\omega _+`$ mode stays at the energy of 3 meV, but $`\omega _{}`$ is lowered. The splitting of the modes results from the fact that the symmetry in $`x`$\- and $`y`$-directions is broken in the confinement potential for a non-zero $`d`$. For our Hamiltonian, the confinement along $`y`$-direction is kept parabolic even for $`d>0`$, and the strength in this direction is 3 meV. For the $`x`$-direction, the potential is non-parabolic for $`d>0`$, and the effective strength of an approximating parabolic confinement is smaller than 3 meV for any finite $`d`$. In this way, the $`B=0`$ modes correspond to excitations to either the short or long axis of the system. This fact can be further demonstrated by considering the spectra of $`x`$\- and $`y`$-polarized radiation, presented in Fig. 2 (a).

One can see that at $`B=0`$ the transitions are given purely by either of the two linear polarizations. For a non-zero $`B`$, the two linear polarizations mix due to a term $`𝐀`$ in the Hamiltonian. Measurements using linearly polarized radiation lead to similar conclusions, see Refs. \[Dahl, ; Dahl2, ; Hochgrafe, \]. For a strong magnetic field, the magnetic field length ($`l_B=\sqrt{\mathrm{}/eB}`$) becomes smaller than the confinement length ($`l_0=\sqrt{\mathrm{}/\omega _0m^{}}`$) at $`B1.7`$ T. Thus at greater $`B`$, the harmonic Landau potential dominates over the QDM potential and the electrons start to localize in the minima of the confinement potential Ari . Therefore, the non-parabolicity is diminished. The same effect can also be seen in the absorption of circularly polarized radiation, shown in Fig. 2.

To separate the effects of the non-parabolic confinement potential from those of interactions on the FIR spectra, we have performed calculations with reduced Coulomb interaction strengths $`0C<1`$. We present results for the $`d=30`$ nm case as it shows highly non-trivial spectra. Figs. 1 (g) and (h) show a non-interacting spectra $`C=0`$, and Fig. 1 (e) shows a $`S=0`$ spectrum for $`C=0.2`$ and Fig. 1 (f) the $`S=1`$ one for $`C=0.5`$. Different intermediate values of $`C`$ are chosen for different $`S`$ as the ones shown are the most representative ones. Comparison of the $`C<1`$ spectra to those of Figs. 1 (c) and (d) with $`C=1`$ shows a very surprising feature: the spectra for the reduced interaction strengths differ much more from the Kohn modes than the fully interacting cases. The splitting of the modes at $`B=0`$ grows as the interactions are reduced, the anticrossing gaps are in general clearer, although a direct comparison is not very easy as there are many shifts in the anticrossing positions. In addition, the $`\omega _{+2}`$ mode of $`S=1`$ becomes more distinct in the spectra as $`C`$ is reduced.

We start the analysis of the $`C<1`$ spectra from the zero field splitting of the two modes. The potential along the $`y`$-axis is parabolic, and $`\omega _+`$ at $`B=0`$ is not affected by the interactions. On the other hand, the excitation energy along the non-parabolic axis is clearly influenced by scaling $`C`$. The Coulomb repulsion effectively steepens the confinement which leads to an increase in the excitation energy.

In experiments for elliptic QD lattices, the zero field resonance of $`\omega _{}`$ was found to exhibit an unexplained and clear dependence on the number of electrons in the QD Hochgrafe . The lowest energy was observed with the smallest number of electrons. As the non-interacting $`S=0`$ spectrum coincides with the spectrum of a one-electron QDM, one can study the particle number dependence of $`\omega _{}`$ for the smallest particle numbers. Interestingly, we find a very similar dependence in our data by comparing the corresponding spectra in Figs. 1 (g) and (c). The reduction of the $`B=0`$ excitation in Ref. Hochgrafe  occurs only in the long-axis direction. This implies a non-parabolic potential in that direction, and on the other hand, a parabolic one in the short-axis direction.

Another interesting question is how the anticrossings change when the interactions are scaled down with $`C<1`$. For all interaction strengths, the energy gap at the anticrossing point gets smaller as one moves to stronger $`B`$. This can be seen to result from the localizing effect of the magnetic field: as the electrons localize around the potential minima, they feel a more parabolic external confinement. Now when the interactions are made weaker, the positions of the anticrossings move, in general, to higher $`B`$. However, the energy gaps associated to fully and non-interacting cases are nearly the same; the change is typically less than 10%. Thus the reduced interaction is compensated by the stronger localizing effect of $`B`$ at the new anticrossing position of higher $`B`$. This can be understood by noting that interactions enhance localization: the higher kinetic energy of localization is compensated by the reduced interaction energy. In the extreme limit, a Wigner molecule of electrons is formed AriCL . It is now tempting to interpret this balance between the two localizing effects in the way that the interactions affect the anticrossings only indirectly via the effective potential. Furthermore, as one cannot evidently identify any clear detail of the spectra to result from the electron-electron interactions, the possibility of FIR spectroscopy to reveal correlation effects is somewhat questionable. Its ability, however, to probe for the effective confinement potential of electrons is clear. This is especially true if various polarizations of radiation are used.

To summarize, we have presented FIR spectra of a lateral two-electron quantum-dot molecule and studied the effects of confinement potential and correlations on the spectra. We have shown that lowering the symmetry of the confinement potential induces changes in the spectra. These changes are deviations from the two Kohn modes, and include the splitting of zero-field resonances, anticrossings in $`\omega _+`$, and an additional mode above $`\omega _+`$ in the spin-triplet spectrum. We have studied the effect of the electron-electron interactions on these deviations, and surprisingly, the interactions drive spectra towards those of a high-symmetry QD. All the deviations from Kohn modes listed above are reduced as the interactions are turned on. This leads to an interesting conclusion: the deviations from the Kohn modes result merely from the non-parabolic external confinement potential. The interactions change the effective potential of the electrons, but as our results show, the resulting potential is more parabolic than the bare external confinement. In this way, the resulting FIR spectra are closer to Kohn modes than the non-interacting ones. In addition, our results show that probing interaction effects by FIR spectroscopy is not straightforward, as clear signatures of interactions are not present even in the two-electron case where the interaction effects should be strongest. On the other hand, FIR spectroscopy is able to reveal the non-parabolicity of the external confinement potential very efficiently, particularly if polarized radiation is used.

###### Acknowledgements.

We acknowledge support by the Academy of Finland’s Centers of Excellence Program (2000-2005).


format: text

text:
# NUMERICAL EXAMINATION OF THE STABILITY OF AN EXACT TWO-DIMENSIONAL SOLUTION FOR FLUX PILE-UP MAGNETIC RECONNECTION

## 1 INTRODUCTION

The current consensus is that the explosive release of magnetic energy by virtue of magnetic reconnection in the corona is responsible for the observed solar flares (e.g. Tsuneta, 1996; Shibata, 1999; Priest & Forbes, 2000; Moore et al., 2001). Reconnecting current sheets in the solar corona are believed to form and evolve on the scale of days and weeks, allowing the free magnetic energy in the corona to be accumulated in the sheet. Eventually the sheet is disrupted by an instability, leading to the flare energy release on the scale of minutes and hours.

A well-known difficulty in the flare theory is the so-called trigger problem, which relates to the fact that a relatively slow evolution of the sheet, described by ideal magnetohydrodynamics (MHD), is followed by rapid destabilization when resistive effects in the sheet become important (e.g. Wang & Bhattacharjee, 1994; Hirose et al., 2001). The Kelvin–Helmholtz (KH) and tearing instabilities are typically invoked to interpret the transition to the latter stage of the evolution (e.g. Knoll & Chacon, 2002; Biskamp, 1994, for a review).

Given possible applications to the problem of current sheet evolution and disruption in the solar corona, it is interesting to determine whether the standard results on the KH and tearing instabilities are robust with respect to the influence of the reconnection flow in the sheet and the global geometry of the reconnection solution. This provides the motivation for the present work.

A weakness of most previous analyses of the current sheet instabilities is that they assumed an initial unperturbed state that was not a solution to the resistive steady MHD equations. Previous studies typically assumed a “quasi-equilibrium” one-dimensional $`\mathrm{tanh}`$-like profile of magnetic field either with zero velocity (Malara et al., 1991, 1992; Dahlburg et al., 1992) or with a one-dimensional velocity proportional to the magnetic field (Dahlburg & Einaudi, 2001). Effects of the plasma flow and the resistive diffusion were not exactly balanced for these initial conditions unless a nonuniform electric resistivity in the sheet was postulated (Schumacher & Seehafer, 2000). Left alone, such current sheet would resistively decay. One might also argue that heuristic one-dimensional models for the unperturbed magnetic and velocity field could only describe the immediate vicinity of the current sheet rather than a global MHD solution for magnetic reconnection. We are aware of one work in which the unperturbed state is an exact solution of MHD equations with a uniform and constant electric resistivity: Phan & Sonnerup (1991) investigated linear stability against tearing mode of a one-dimensional current sheet supported by a stagnation-point flow in the plane perpendicular to the magnetic field vector.

In this paper, we use MHD simulations to examine the stability of a reconnecting current sheet described by an exact two-dimensional solution for magnetic reconnection (Craig & Henton, 1995). The Craig–Henton solution provides analytical description of flux pile-up magnetic reconnection in incompressible plasma of arbitrary uniform electric resistivity. The resistive diffusion is exactly balanced by a sheared stagnation-point flow, so that steady-state MHD equations are satisfied globally. Notably, we limit consideration to two-dimensional instabilities of the Craig–Henton solution, although previous studies indicated that three-dimensional ideal instabilities of two-dimensional current sheets can grow much faster than any two-dimensional instabilities (e.g. Schumacher & Seehafer, 2000; Dahlburg et al., 1992).

Recall that flux pile-up current sheets have been repeatedly studied and generalized not only in two but also in three dimensions (e.g. Craig et al., 1995; Litvinenko & Craig, 2000). It is understood that both the Craig–Henton solution and those considered by Sonnerup & Priest (1975) and Phan & Sonnerup (1991) are particular cases of a general family of solutions for flux pile-up merging and reconnection. The analytical nature of the solutions made it possible to demonstrate explicitly many of the characteristics of magnetic reconnection, such as the appearance of small length scales defined by resistivity, magnetic sling shots, and Alfvénic outflows (e.g. Litvinenko & Craig, 1999, 2000). It is worth stressing that, although the pile-up reconnection solution is locally equivalent to the Sweet–Parker current sheet as far as reconnection scalings with resistivity are concerned, the two models are not identical globally. The Sweet–Parker model postulates a uniform advection region outside the current sheet, whereas the Craig–Henton solution explicitly describes the magnetic field amplification caused by the non-uniform velocity field outside the sheet.

The paper is organized as follows. In Section 2, we briefly review the steady-state flux pile-up reconnection solutions in two dimensions and identify two types of the solutions depending on the sign of a plasma velocity component at the boundary. In Section 3, we present the numerical methods for examining the stability of the solutions and describe the numerical results. We compare the results with predictions of the linear theory of the KH and tearing instabilities, discuss the limitations of the local approach and suggest possible solar applications in Section 4. We summarize in Section 5.

## 2 REVIEW OF THE EXACT FLUX PILE-UP RECONNECTION SOLUTION

We begin by reviewing the Craig–Henton solution that satisfies the steady-state equations of incompressible non-viscous resistive MHD in a planar geometry. Using the Poisson bracket $`[\psi ,\varphi ]\psi _x\varphi _y\psi _y\varphi _x`$, the non-dimensionalized basic equations can be written as

$`[^2\varphi ,\varphi ]=[^2\psi ,\psi ],`$ (1)
$`E+[\psi ,\varphi ]=\eta ^2\psi ,`$ (2)
where $`\psi `$ and $`\varphi `$ are, respectively, a flux function for the magnetic field $`𝑩=\times \psi \widehat{𝒛}`$ and a stream function for the velocity field $`𝒗=\times \varphi \widehat{𝒛}`$. The problem is specified by two dimensionless parameters $`E`$ and $`\eta `$. The uniform electric field $`E`$, normalized by $`v_\mathrm{A}B_0`$ (in electromagnetic units), characterizes the reconnection rate. The inverse Lundquist number $`\eta =\eta _0/(v_\mathrm{A}L)`$ is based on a uniform magnetic diffusivity $`\eta _0`$. Here and in what follows, $`B_0`$, $`\rho _0`$, and $`L`$ are the reference values of the magnetic field, plasma density, and spatial length scale, $`v_\mathrm{A}=B_0/\sqrt{4\pi \rho _0}`$ is the reference Alfvén speed.

An exact reconnection solution satisfying the equations above is given by

$`𝑩(x,y)=(\beta x,\beta y)+(0,{\displaystyle \frac{E}{\eta \mu }}\mathrm{daw}\left(\mu x\right)),`$ (3)
$`𝒗(x,y)=(\alpha x,\alpha y)+(0,{\displaystyle \frac{\beta }{\alpha }}{\displaystyle \frac{E}{\eta \mu }}\mathrm{daw}\left(\mu x\right)),`$ (4)
where $`\mathrm{daw}(x)_0^x\mathrm{exp}\left(t^2x^2\right)𝑑t`$ is the Dawson function, and the dimensionless parameters $`\alpha `$ and $`\beta `$ define, respectively, the amplitude and the shear of the velocity field (Craig & Henton, 1995). These parameters are used to define $`\mu `$ according to

$$\mu ^2=\frac{\beta ^2\alpha ^2}{2\alpha \eta }.$$
(5)
Clearly $`\alpha `$ and $`\beta `$ should satisfy $`(\beta ^2\alpha ^2)/\alpha >0`$ for the physically interesting case of a localized current sheet, in which case the sheet thickness scales as $`\mu ^1`$.

The distribution of plasma pressure,

$$p(x,y)=p_0\frac{1}{2}\left(\alpha ^2\left(x^2+y^2\right)+\left(\frac{E}{\eta \mu }\mathrm{daw}\left(\mu x\right)\right)^2\right)\beta y\frac{E}{\eta \mu }\mathrm{daw}\left(\mu x\right),$$
(6)
is derived from the original momentum equation

$$\left(𝒗\right)𝒗=^2\psi \psi p,$$
(7)
assuming a constant base pressure $`p_0`$.

Equations (3) and (4) demonstrate that both the magnetic field and the velocity field have a stagnation-flow component proportional to $`(x,y)`$, and a shear-flow component proportional to $`(0,\mathrm{daw}(\mu x))`$. The shear width scales as $`\mu ^1`$.

Assuming that $`\beta 0`$, the solutions are categorized into two groups according to the flow direction along the $`x`$-axis at the external boundary: the outflow solutions with $`\alpha >0`$ and $`\beta >\alpha `$ and the inflow solutions with $`\alpha <0`$ and $`|\alpha |>\beta `$. In the outflow solution the flux pile-up is accomplished mainly by the advection of the $`x`$-component of magnetic field by the flow parallel to the $`y`$-axis. In the inflow solution $`\beta =0`$ corresponds to the anti-parallel merging with stagnation flow (Sonnerup & Priest, 1975), and $`\beta >0`$ describes the fully two-dimensional solution for flux pile-up reconnection, which is also called “reconnective annihilation” (Priest & Forbes, 2000).

Viscous effects, ignored in the solution above and hence throughout this paper, are likely to be important in the region of high shear coinciding with the current sheet. This is potentially more important for the outflow solution that is characterized by a strong velocity shear. The Craig–Henton formalism has been generalized to include viscosity by Fabling & Craig (1996). The latter solution is more complex mathematically. This is why in this paper we concentrate on a simpler non-viscous case in which a simple analytical expression for the initial state is available. Although viscosity effects are neglected without much justification in the outflow case, we wish to stress from the outset that the neglect is justified much better in the inflow case when the flow is only weakly sheared. As we show in Section 4.4, it is the inflow case that is of primary interest to us in the solar flare context. The inflow solution can be tearing unstable, which may explain flare energy release by tearing-induced reconnection.

## 3 NUMERICAL EXAMINATION OF THE STABILITY OF THE EXACT SOLUTIONS

Typical solutions in the outflow case and the inflow case are shown in Figure 1. As noted in the previous section, a shear of width $`2\mu ^1`$ is present along the $`y`$-axis in both the magnetic and velocity fields. Therefore, we can expect that these solutions can be subject to the KH instability and the tearing instability. We investigated the stability of the solutions using MHD numerical simulations. Specifically we adopted the analytical exact solution as the initial condition, perturbed it slightly, and examined the time development of the system.

### 3.1 Numerical Method

The basic equations for the numerical simulations are compressible resistive MHD equations in the following dimensionless form:

$`{\displaystyle \frac{\rho }{t}}+\left(\rho 𝒗\right)=0,`$ (8)
$`{\displaystyle \frac{}{t}}\left(\rho 𝒗\right)+\left(\rho 𝒗𝒗𝑩𝑩+\left(p+{\displaystyle \frac{𝑩^2}{2}}\right)𝐈\right)=0,`$ (9)
$`{\displaystyle \frac{𝑩}{t}}+\times 𝑬=0,`$ (10)
$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{1}{2}}\rho 𝒗^2+{\displaystyle \frac{p}{\gamma 1}}+{\displaystyle \frac{𝑩^2}{2}}\right)+\left(\left({\displaystyle \frac{1}{2}}\rho 𝒗^2+{\displaystyle \frac{\gamma p}{\gamma 1}}\right)𝒗+𝑬\times 𝑩\right)=0,`$ (11)
$`𝑬=𝒗\times 𝑩+\eta \times 𝑩,`$ (12)
where $`\rho `$ is the plasma density and $`\gamma =5/3`$ is the ratio of specific heats. The above equations are normalized in the same way as those in Craig & Henton (1995) and in the previous section, using the reference values of length $`L`$, density $`\rho _0`$, and Alfvén speed $`v_\mathrm{A}`$. The time is normalized with the global Alfvén time $`t_\mathrm{A}L/v_\mathrm{A}`$.

To simulate the Craig–Henton solutions derived for incompressible flow, we set the constant base pressure $`p_0`$ in equation (6) large enough so that the sound speed is at least ten times greater than both the fluid speed and the Alfvén speed. This ensures that the compressible fluid behaves almost as the incompressible one.

The numerical scheme used in the simulations is based on an upwind scheme with Roe’s numerical fluxes (Powell et al., 1999). Second-order accuracy both in space and time is accomplished by variable extrapolation of Monotone Upstream-centered Schemes for Conservation Laws (MUSCL) and multi-step time integration, guaranteeing Total Variation Diminishing (TVD) with a non-linear limiter function (e.g. Hirsh, 1990). Constraint Transport (CT) scheme is combined with the above scheme to develop the magnetic field, ensuring $`𝑩=0`$ (Balsara & Spicer, 1999).

The simulation box is set to $`|x|<2.5,|y|<2.5`$, and $`500\times 500`$ constant grids are assigned. This resolution is determined so that at least ten grids are assigned in the current sheet thickness. The analytical solutions of Craig & Henton (1995) are adopted as the initial conditions for the velocity, magnetic field, and pressure. The initial dimensionless density is set to unity: $`\rho (x,y,0)=1`$. In addition, a small random velocity perturbation (1% of the typical Alfvén velocity $`v_\mathrm{A}`$) is added in the region of $`|x|<1,|y|<1`$. The inverse Lundquist number is uniform and constant: $`\eta (x,y,t)=\eta `$. The boundary values are fixed to the initial values for the whole simulation time. The base dimensionless pressure $`p_0`$ is fixed at $`510^4`$ in all cases. The simulations are performed with the fixed Courant number of 0.8, and the corresponding number of time steps required for a typical run is of order $`10^5`$.

As far as the accuracy of the numerical method is concerned, we verified that the density fluctuations in the simulation are typically of order 0.01% and never exceed 0.54%. The typical value of $`𝐯`$, normalized by the ratio of the sound speed and the grid size, is of order $`10^4`$, which is consistent with the magnitude of density fluctuations. Therefore the above value for the base pressure $`p_0`$ is indeed large enough to mimic the incompressible behavior. We also confirmed that the compressible code can maintain the unperturbed equilibrium for stable cases. Additionally we performed convergence tests for several cases and confirmed that the stability does not depend on the resolution.

### 3.2 Overview of the Numerical Results

Parameters of the solutions of which we examined the stability are summarized in Table 1. In Table 1 and hereafter, we adopt the Alfvén time $`T_\mathrm{A}L/V_\mathrm{A}`$ for the normalization of time; $`T_\mathrm{A}`$ is defined in each case based on the Alfvén speed $`V_\mathrm{A}`$ at the entrance to the current sheet. This normalization allows us to separate the instability effects under study from amplification of magnetic field carried by the reconnection flow. The relation between $`T_\mathrm{A}`$ and $`t_\mathrm{A}=L/v_\mathrm{A}`$ used in the normalization of the basic equations is

$$\frac{T_\mathrm{A}}{t_\mathrm{A}}=\frac{v_\mathrm{A}}{V_\mathrm{A}}\frac{\mu \eta }{E},$$
(13)
since $`v_\mathrm{A}`$ is the Alfvén speed at the location where the magnitude of the magnetic field in equation (3) is unity. Here we have defined $`V_\mathrm{A}`$ ignoring a factor $`0.5`$, which is the maximum value of the Dawson function, for simplicity.

The growth of the system in each case is summarized in Figure 2. We quantify the perturbation of the velocity distribution in the sheet by plotting the time evolution of

$$I(t^{})\mathrm{log}\left(_{|x|<\mu ^1,|y|<2}|𝒗(x,y,t^{})𝒗(x,y,0)|𝑑x𝑑y/_{|x|<\mu ^1,|y|<2}𝑑x𝑑y\right),$$
(14)
where $`t^{}`$ is the time normalized by $`T_\mathrm{A}`$. The duration of the simulation is $`0<t^{}<5`$ for the outflow solutions and $`0<t^{}<50`$ for the inflow solutions. We stopped the simulations before the evolution of the system could be affected by the fixed boundary conditions.

Figure 2 shows that a choice of parameters leads to instability of both the outflow and the inflow solutions. The nature of the detected instability is, however, entirely different for the two types of solutions. This can be clearly seen in Figure 3, in which the time evolution of some typical unstable solutions in both cases is shown. In the outflow solutions (Fig. 3 (a)(b)), the current sheet along the $`y`$-axis is warped, and then vortices with magnetic islands are formed, suggesting that the KH instability occurs. On the other hand, in the inflow solutions (Fig. 3 (c)(d)), the time development of the current distribution is symmetric with respect to the $`y`$-axis, and long and narrow magnetic islands are formed, which are characteristic features of the tearing instability.

The difference of the unstable behavior comes from the fact that the velocity shear is stronger than the magnetic shear in the outflow solutions, $`\beta >\alpha `$, and vice versa in the inflow solutions, $`\beta <|\alpha |`$ (see the terms multiplying the Dawson function in eqs.  and ). These features can be seen in Figures 4 and 5, which show the distributions of the ratio of the kinetic energy to the magnetic energy. In the outflow (inflow) solutions, the kinetic (magnetic) energy is dominant near the current sheet.

## 4 INTERPRETATION OF THE NUMERICAL RESULTS

In this section, we discuss the parameter dependencies of the instabilities on the basis of the standard linear theories of the KH instability for the outflow solutions and the tearing instability for the inflow solutions. We also point out that the standard results, derived for uniform initial conditions outside the current sheet, have limited applicability to the pile-up solutions, characterized by strongly nonuniform profiles of the magnetic and velocity fields outside the sheet.

### 4.1 KH Instability of the Outflow Solution

The standard linear growth rate $`\omega `$ for the KH instability is given by the following expression (Chandrasekhar, 1981; Michael, 1955):

$$\omega =k\sqrt{\left(\frac{\mathrm{\Delta }U}{2}\right)^2V_\mathrm{A}^2},$$
(15)
where $`k`$ is the wave number, $`V_\mathrm{A}`$ is the local Alfvén speed, and $`\mathrm{\Delta }U`$ is the difference in plasma velocities at the opposite edges of the sheet. Note that Chandrasekhar (1981) discussed in detail the case of a uniform magnetic field whereas Michael (1955) allowed for different magnetic field strengths on the two sides of a discontinuity. The expressions for $`\omega `$ are identical, however, if the magnetic field is symmetrical on either side of the current sheet, which is the case for the Craig–Henton solution.

When applying the formula for the instability growth rate to the Craig–Henton solution, it should be remembered that a continuous velocity profile has the effect of limiting the KH instability to large wave lengths, so that the most unstable mode has $`k`$ of order the inverse of the sheet thickness $`l`$ (Chandrasekhar, 1981). In this case, since $`l=L/\mu `$,

$$k_{\mathrm{max}}=\mu /L.$$
(16)
Recalling equations (4) and (13), it is straightforward to determine the maximum dimensionless linear growth rate:

$$\omega _{\mathrm{max}}t_\mathrm{A}=(k_{\mathrm{max}}L)\sqrt{\left(\frac{\mathrm{\Delta }U/2}{v_\mathrm{A}}\right)^2\left(\frac{V_\mathrm{A}}{v_\mathrm{A}}\right)^2}=\mu \sqrt{\left(\frac{\beta }{\alpha }\frac{E}{\eta \mu }\right)^2\left(\frac{E}{\eta \mu }\right)^2}=\frac{E}{\eta }\frac{\sqrt{\beta ^2\alpha ^2}}{\alpha }.$$
(17)
The corresponding growth rate normalized with $`T_\mathrm{A}`$ (see eq. ) is given by

$$\omega _{\mathrm{max}}T_\mathrm{A}=\omega _{\mathrm{max}}t_\mathrm{A}\frac{T_\mathrm{A}}{t_\mathrm{A}}=\frac{\beta ^2\alpha ^2}{\sqrt{2\alpha ^3\eta }}.$$
(18)
Linear theory predicts that, in order to overcome the stabilizing effect of the magnetic field, the shear has to be large enough, $`\beta >\alpha `$, which is the case for the outflow solution.

A major difference between the geometry of the standard KH instability and the outflow Craig–Henton solution is the strongly nonuniform profile of the reconnecting magnetic field and the associated presence of the current sheet at the surface $`x=0`$. The stabilization condition, however, should not depend on whether or not the magnetic field is uniform. In fact it is clear on physical grounds that in either case the solution is stable as long as the average Alfvén speed exceeds the fluid velocity difference (e.g. Biskamp, 2003). On the other hand, we will see below that the growth rate of unstable solutions can be strongly modified by the global nonuniformity of the field outside the current sheet.

Note that the resistivity in our simulations is small enough for the time scale of the instability $`(\omega _{\mathrm{max}})^1`$ to be shorter than the diffusion time scale $`\tau _\mathrm{d}`$ over the current sheet thickness $`l`$. The diffusion time scale is defined as $`\tau _\mathrm{d}=l^2/\eta _0`$ or $`\tau _\mathrm{d}/T_\mathrm{A}=E/(\mu ^3\eta ^2)`$ in the dimensionless form. The ratio of $`\tau _\mathrm{d}/T_\mathrm{A}`$ to $`(\omega _{\mathrm{max}}T_\mathrm{A})^1`$ is typically several tens, and its smallest value is $`5.77`$ in Cases C and D.

The dimensionless wave lengths of the fastest growing modes in the unstable cases can be explained by the linear theory, but the growth rates are smaller by an order of magnitude, compared with the predicted values in Table 1. For example, in Case B, the fastest growing mode has the wave length of $`0.5`$ (see Fig. 3(a): $`t^{}=3.0`$), which agrees with the value predicted by the linear theory ($`=0.51`$), while the growth rate ($`2.7`$) is much smaller than the predicted value ($`=21.21`$). One factor that may be responsible for this significant reduction in the growth rate is the magnetic field gradient outside the sheet, which is associated with the reverse electric current in the Craig–Henton solution. The magnetic field pile-up on either side of the current sheet corresponds to the field variation of order $`E/(\eta \mu )`$ over the length scale $`\mu ^1`$ (eq. ). Since the corresponding energy density $`(E/(\eta \mu ))^2`$ is comparable to the kinetic energy of the shearing motion (eq. ), the KH instability may be suppressed by the nonuniform magnetic field outside the sheet.

### 4.2 Tearing Instability of the Inflow Solution

The maximum growth rate of the tearing instability and the corresponding wave number are given by the following formulas, ignoring factors of order unity (Furth, Killeen, & Rosenbluth, 1963):

$`\omega _{\mathrm{max}}={\displaystyle \frac{1}{\sqrt{\tau _\mathrm{d}\tau _\mathrm{A}}}},`$ (19)
$`k_{\mathrm{max}}=\left({\displaystyle \frac{\tau _\mathrm{A}}{\tau _\mathrm{d}}}\right)^{\frac{1}{4}}{\displaystyle \frac{1}{l}},`$ (20)
where $`\tau _\mathrm{d}=l^2/\eta _0`$ and $`\tau _\mathrm{A}=l/V_\mathrm{A}`$ are the diffusion and Alfvén time scales in the current sheet of thickness $`l`$. Using equations $`\eta =\eta _0/(v_\mathrm{A}L),l/L=1/\mu `$, and (13), we have the following expressions for the dimensionless linear growth rate and wave number for the tearing instability:

$`\omega _{\mathrm{max}}t_\mathrm{A}=\sqrt{{\displaystyle \frac{V_\mathrm{A}}{v_\mathrm{A}}}{\displaystyle \frac{\eta _0}{v_\mathrm{A}L}}\left({\displaystyle \frac{L}{l}}\right)^3}=\sqrt{{\displaystyle \frac{E}{\eta \mu }}\eta \mu ^3}=\sqrt{{\displaystyle \frac{\left(\beta ^2\alpha ^2\right)E}{2\alpha \eta }}},`$ (21)
$`k_{\mathrm{max}}L=\left({\displaystyle \frac{\eta _0}{Lv_\mathrm{A}}}{\displaystyle \frac{L}{l}}{\displaystyle \frac{v_\mathrm{A}}{V_\mathrm{A}}}\right)^{\frac{1}{4}}{\displaystyle \frac{L}{l}}=\left(\eta \mu {\displaystyle \frac{\mu \eta }{E}}\right)^{\frac{1}{4}}\mu =\left({\displaystyle \frac{\left(\beta ^2\alpha ^2\right)^3}{8\alpha ^3\eta E}}\right)^{\frac{1}{4}},`$ (22)
where it should be remembered that $`\alpha <0`$ and $`\beta <|\alpha |`$. The corresponding growth rate normalized with $`T_\mathrm{A}`$ is

$$\omega _{\mathrm{max}}T_\mathrm{A}=\omega _{\mathrm{max}}t_\mathrm{A}\frac{T_\mathrm{A}}{t_\mathrm{A}}=\frac{\alpha ^2\beta ^2}{2|\alpha |\sqrt{E}}.$$
(23)
Thus stronger inflows (large $`\alpha `$) and weaker two-dimensionality (small $`\beta `$) are the factors that make the current sheet more unstable with respect to tearing. Theory demands that $`\lambda _{\mathrm{max}}/l`$ be greater than $`2\pi `$ for the instability, which is the case for the inflow solutions (see Table 1).

The analytical estimate for $`\omega _{\mathrm{max}}t_\mathrm{A}`$ ignores the effect of the reconnection flow on the instability. It is easy to see on physical grounds that a strong flow would carry the growing modes out of the reconnection region, thus suppressing the instability (see Biskamp, 1994, and references therein). A rough condition for such stabilization is simply $`|\alpha |>\omega _{\mathrm{max}}t_\mathrm{A}`$, which becomes

$$\frac{|\alpha |}{1\beta ^2/\alpha ^2}>\frac{E}{\eta }$$
(24)
for the flux pile-up reconnection solution under consideration. For typical simulation parameters, however, this condition is difficult to achieve unless $`\beta |\alpha |`$, indicating a broad current sheet extending to the external boundaries.

Chen & Morrison (1990) and Ofman et al. (1991) studied the influence of the velocity shear on the tearing instability and found that the instability is suppressed when the flow shear is larger than the magnetic field shear at the neutral plane. It is worth remarking that the flow geometry assumed in our study is different. Recall that the inflow Craig–Henton solution is defined by a weak shear, specifically $`\beta <|\alpha |`$. This condition is required in order to have a localized current sheet (see Section 2). Another difference with the problem investigated by Chen & Morrison (1990) is that the Craig–Henton solution possesses the reconnection outflow that also have a stabilizing effect. There is a formal similarity between stabilization by the velocity shear and the reconnection outflow, although the flow geometry is different. One also should remember that when the flow shear is large, the behavior of the reconnecting current sheet is more complex because it is no longer easy to distinguish the effects of the tearing and KH instabilities, and both effects are generally present (Ofman et al., 1991).

Both the wave lengths and the growth rates in the unstable cases can be basically explained by the linear theory. For example, in Case K, the fastest growing mode has the wave length of $`2.5`$ (see Fig. 3(c): $`t^{}=25.0`$ ) and the growth rate of $`0.3`$, which agree fairly well with the values predicted by the linear theory ($`\lambda _{\mathrm{max}}/L=3.51`$ and $`\omega _{\mathrm{max}}T_\mathrm{A}=0.33`$). Furthermore, comparing (c) and (d) in Figure 3, the fastest growing mode in Case N has a shorter wave length and a larger growth rate than that in Case K, which is also consistent with the linear theory (see Table 1).

### 4.3 Effects of Nonuniform Initial Conditions

Equations (18) and (23) predict that the KH instability occurs in the outflow solutions ($`\beta >\alpha `$) and the tearing instability occurs the inflow solutions ($`\beta <|\alpha |`$) for any values of the parameters, $`\alpha `$, $`\beta `$, $`E`$, and $`\eta `$. The instabilities in our simulations can be interpreted using the standard linear theories. We found, however, that the Craig–Henton solution is stable for some values of the parameters. This indicates that the standard linear results, based on uniform initial conditions, are of limited use when analyzing the stability of global nonuniform MHD solutions.

In the application of linear theories to the solutions discussed above, we used the values of the velocity and magnetic field at the entrance to the current sheet. The solutions at hand, however, are neither one-dimensional nor uniform. We expect that, if the length scale of one-dimensionality is not large enough compared with the wave length for the most unstable mode, the instabilities will be modified owing to the resulting interaction with the external fields. In Figure 4 (Figure 5), the dark (light) color region outside the current sheet is the region where the KH (tearing) instability will be influenced by the external fields, since the magnetic (kinetic) energy is dominant there. For example, in Case F, the kinetic-energy dominant region (light color region) is smaller than that in Case B, and thus the KH instability occurs in a smaller region compared with Case B (see Fig. 3 (a)(b)).

It may be argued, therefore, that the instabilities can be suppressed when the effects of two-dimensionality and nonuniformity, ignored in the standard treatment, become important in the initial configuration. Equations (3) and (4), describing the structure of the Craig–Henton solution, can be used to quantify the expected effect. The two-dimensional magnetic field in the outflow solution, for example, should start influencing the shearing velocity profile when the first term of equation (3) becomes large compared with the second term of equation (4), that is when $`\alpha `$ increases (Case H), $`E`$ decreases (Case C), or $`\eta `$ increases (Case D). This may explain why cases C (small $`E`$) and D (large $`\eta `$) in the outflow solution appear to be stable. A similar argument can be given for the inflow solution by comparing the first term of equation (4) and the second term of equation (3). The two-dimensionality should come into play when $`\alpha `$ increases (Case Q), $`\eta `$ increases (Case M), or $`E`$ decreases (Case L). Perhaps this is why we did not detect instability in cases L (small $`E`$), M (large $`\eta `$), and Q (large $`\alpha `$) in the inflow solution.

The lack of instability in cases L, M, and Q is also consistent with the qualitative equation (24). Interestingly, it appears that this order-of-magnitude stability condition can be made quantitative if its left-hand side is multiplied by a factor of order $`10^2`$. Recall that a similarly large factor for the aspect ratio of a dynamic Sweet–Parker current sheet, stable with respect to tearing, was suggested by Biskamp (1994).

### 4.4 Tearing in the Solar Corona

This work is partly motivated by the idea that tearing of a large-scale current sheet in the solar corona can lead to impulsive energy release in solar flares (e.g. Sturrock, 1994). The flux pile-up solution investigated in this paper is a member of the only exact analytical MHD family of solutions for magnetic reconnection that we are aware of, and analytical scalings based on these solutions may be of general interest. One important question is whether the growth rate of the tearing mode is sufficiently rapid if the realistic electric resistivity in the solar corona is assumed. Our numerical results are in quantitative agreement with the linear theory of tearing instability as long as the stabilizing effects of two-dimensional geometry are small enough. We now proceed to show that the instability scalings do imply the possibility of rapid conversion of magnetic energy in a large-scale current sheet in the corona.

The time-scale of destabilization, caused by tearing of a current sheet, is conveniently defined as $`\tau =\omega _{\mathrm{max}}^1`$. As discussed above, our simulations confirm that the growth rate for the flux pile-up solution $`\omega _{\mathrm{max}}`$ can be estimated using equation (21). While solar observations put some constraints on the dimensions of reconnection regions, reconnection rates, and the electric resistivity, the observations are not detailed enough to specify the geometry of reconnective flows, given by the parameters $`\alpha `$ and $`\beta `$. This is why, in order to apply the formula for $`\omega _{\mathrm{max}}`$ to unstable current sheets in solar active regions, we rewrite equation (21) as follows:

$$\frac{\tau }{t_\mathrm{A}}=\frac{S^{1/2}}{\stackrel{~}{B}_s^{1/2}}\left(\frac{l}{L}\right)^{3/2},$$
(25)
where $`S=\eta ^1`$ is the Lundquist number, and the dimensionless pile-up factor $`\stackrel{~}{B}_s=V_A/v_A`$ is simply the magnetic field at the entrance to the sheet relative to the reference value $`B_0`$. This equation for $`\tau `$ reduces to the standard result for tearing in the Sweet–Parker current sheet, provided the pile-up is absent, in which case $`\stackrel{~}{B}_s=1`$ and $`l/L=S^{1/2}`$.

Adopting $`L=10^{9.5}`$ cm, $`B_0=10^2`$ G, and $`v_A=10^9`$ cm s<sup>-1</sup> as the reference values in a solar active region, we can estimate $`S10^{14}`$ based on the classical value of electric resistivity. Models for the structure of flaring current sheets suggest that their thickness is unlikely to exceed $`l10^5`$ cm (e.g. LaRosa, 1992; Litvinenko & Craig, 2000), which leads to $`l/L10^{4.5}`$. Coronal magnetic fields are unlikely to exceed $`10^3`$ G, hence we should assume a modest pile-up factor $`\stackrel{~}{B}_s10`$. These numerical values result in the instability time-scale of order the Alfvén time $`t_A10^{0.5}`$ s. Although nonlinear effects eventually come into play, this estimate for $`\tau `$ confirms that pile-up current sheets can be unstable to tearing and can be destroyed quite rapidly in the corona. This conclusion should remain qualitatively valid if current-driven instabilities give rise to plasma turbulence and enhanced electric resistivity in the sheet. Although the sheet should thicken up due to increased resistivity, this effect will be partly balanced by the corresponding decrease in the effective Lundquist number.

A relevant question at this point is whether flux pile-up reconnection solutions adopted for our simulations are possible in the solar corona at all. It might be argued that, since the Craig–Henton solution is unstable for a wide range of parameters, it is not realized on the Sun. We wish to stress, however, that the simulations also identified a parameter regime in which the Craig–Henton solution is stable. In particular the stability with respect to tearing in cases L, M, and Q is consistent with the qualitative stabilization condition given by equation (24). The stability appears to be intimately related to the two-dimensionality of the solution. Equation (24) indicates that instability is easier to achieve when $`\beta 0`$, which corresponds to one-dimensional merging of magnetic fields. The shear parameter $`\beta `$ would be difficult to determine observationally. As in our estimate for $`\tau `$, however, we can eliminate $`\beta `$ in equation (24) by specifying the current sheet thickness $`l`$. In order to give an illustrative example, we assume $`\alpha 1`$ and $`l10^6`$ cm. Then the stabilization condition can be written as $`E<(l/L)^210^7`$, which corresponds to $`E<10^4`$ V cm<sup>-1</sup> in dimensional units. We suggest that stable two-dimensional solutions can describe relatively thick (large $`l`$), slowly reconnecting (small $`E`$), stable current sheets in the solar corona.

Finally, based on this sensitivity of our numerical results to the global geometry of magnetic reconnection, it would be tempting to suggest that changes in the geometry of the solution may be responsible for rapid changes in the rate of magnetic energy release by magnetic reconnection and thus may trigger solar flares. One caveat is that we only performed a parametric study of stability of the steady flux pile-up solutions. A rigorous approach to the trigger problem would require the solution of a much more complicated task: we would have to track the dynamical evolution of a current sheet with time-dependent parameters that reflect, for example, a slow change in boundary conditions at the photosphere due to the emergence of magnetic flux and photospheric flows.

## 5 DISCUSSION

Exact MHD solutions are now available, which describe genuine flux pile-up magnetic reconnection, as opposed to simple one-dimensional merging (Craig et al., 1995). Since these analytical solutions explicitly demonstrate many of the properties of fast reconnection, their analysis can be a powerful tool in the reconnection theory (e.g. Litvinenko & Craig, 2000; Heerikhuisen, Litvinenko, & Craig, 2002). Furthermore, the solutions are known quantitatively to describe the results of numerical simulations performed using a time-dependent code (Heerikhuisen, Craig, & Watson, 2000).

We have examined in this paper the stability of an exact two-dimensional solution for flux pile-up reconnection in incompressible resistive plasma (Craig & Henton, 1995). We used numerical simulations to show that the solutions of the outflow type can be Kelvin–Helmholtz (KH) unstable and the solutions of the inflow type can be tearing unstable. Notably, we used the steady Craig–Henton solution as an initial condition, thus allowing enough time for the instability to develop. Our approach should be contrasted with that of a time-dependent simulation (Heerikhuisen, Craig, & Watson, 2000) in which a localized current sheet was present for a finite (perhaps insufficient) time in a doubly periodic geometry. Although the latter simulation provided some evidence that the sheet could break up into magnetic islands, it was not clear whether the tearing instability was responsible for the phenomenon.

In the outflow solutions, where a current sheet is formed by strong shearing flows, the KH instability was observed to develop. The numerical results can be qualitatively explained by the standard linear theory of the KH instability (eq. ). The absolute values of the growth rates, however, are smaller than those predicted. This may be due to the strongly nonuniform profile of the magnetic field at the entrance to the current sheet. The inflow solutions are characterized by a fast, weakly sheared inflow that leads to a strong magnetic field pile-up and the formation of a localized current sheet. These are the conditions under which the tearing instability is likely to develop. The computed growth rates can be explained fairly well by the standard linear theory of the tearing instability (eq. ). The instability, however, is suppressed when the effects of two-dimensionality are large enough to modify the evolution of the flux pile-up current sheets.

Thus the numerical results often can be interpreted by applying standard results for the KH and tearing instabilities to the particular geometries of magnetic field and plasma flow, specified by the Craig–Henton solution. These instabilities cause turbulent and non-steady evolution of the current sheet, which is favorable for fast magnetic reconnection (e.g. Shibata & Tajima, 2001; Tajima & Shibata, 1997) or particle acceleration (e.g. Litvinenko, 2003; Aschwanden, 2002). In some cases, however, the global nature of the solution leads to the suppression of the instabilities. Of particular interest is the dependence of the growth rate and the most unstable wave length on the parameters of the reconnection solution: reconnection rate $`E`$, resistivity $`\eta `$, reconnection flow amplitude $`\alpha `$ and shear $`\beta `$. Finally, we note that the sensitivity of the results to the nonuniform profiles of the magnetic and velocity fields in the reconnection region, as well as the effects of two-dimensionality of the solution, may be responsible for rapid changes in the rate of magnetic energy release by magnetic reconnection, which in turn may help in solving the trigger problem in solar flare research.

Numerical computations were carried out on VPP5000 at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan (project ID: yst15a and hst21a), which is an inter-university research institute of astronomy operated by Ministry of Education, Science, Culture, and Sports, and on VPP5000 at the Super Computer Center, Nagoya University, by the joint research program of the Solar-Terrestrial Environment Laboratory, Nagoya University. This study was initiated as a part of the ACT-JST summer school for numerical simulation in Astrophysics and Space Plasma. One of the authors (Y.E.L.) acknowledges fruitful discussions with I. J. D. Craig and the support by NSF grants ATM-0136718 and IP-9910067 and NASA grant NAG5-10852.


text:
# Measuring the matter distribution within 𝑧=0.2 cluster lenses with XMM-Newton

## INTRODUCTION

Being the most massive gravitationally bound objects in the Universe, clusters of galaxies are prime targets for studies of structure formation and evolution.

Of particular interest are the most X-ray luminous clusters. Indeed, it has been convincingly demonstrated that high X-ray luminosity is a reliable indicator of cluster mass (e.g. \[Allen, 1998\]). Therefore, X-ray selected clusters samples (e.g. \[Ebeling et al., 1996\]) are the best choice for undertaking representative studies of the mass distribution in massive, X-ray luminous clusters. In particular, as they are easily detected out to very large redshift, we can hope to probe cluster evolution back from early epochs. However, X-ray luminosity is also sensitive to the presence of cooling flows, mergers and non-thermal effects, thus the Luminosity-Mass and Temperature-Mass relation needs to be clearly understood before any strong ascertion is made on cluster evolution (\[Smith et al., 2003\]). Proper mass estimate is likely to arise only by combining multiwavelength observations of clusters of galaxies: namely strong and weak lensing observations (\[Smith et al., 2001\]), velocity distribution and dynamics of the cluster galaxies (in a similar way as in \[Czoske et al., 2001\]) and through accurate, spatially resolved X-ray temperature measurements as we are concentrating on here.

Following this premise, we are conducting a survey (Table 1) of the most X-ray luminous clusters ($`L_{X_{[0.12.4keV]}}810^{44}ergs^1`$) in a narrow redshift slice at $`z0.2`$, selected from the XBACS (X-ray Brightest Abell-type Clusters Sample; \[Ebeling et al., 1996\]). As XBACS is restricted to Abell clusters (\[Abell, Corwin & Olowin, 1989\]), it is X-ray flux limited but not truly X-ray selected. However, a comparison with the X-ray selected ROSAT BCS (Brightest Cluster Sample, \[Ebeling et al., 2000\]) shows that $`75\%`$ of the BCS clusters in the redshift and X-ray luminosity range of our sample are in fact Abell clusters. Hence, our XBACS subsample is, in all practical aspects, indistinguishable from an X-ray selected sample.

## THE ANALYSIS PIPELINE

As to the version `5.2` of the official XMM-Newton data analysis software (XMMSAS), the extended sources could not be handled properly, regarding the background subtraction and vignetting corrections. The development of methods, where weights are applied to events according to their energies and positions as well as calibration data so as to correct the event list as if the instrument response were uniformely identical to the on-axis response, started as soon as 2001 (\[Arnaud et al., 2001, Majerowicz et al., 2002, Marty et al., 2002, Marty et al., 2003\]). We concentrated on our side to embedding these algorithms into an environment capable of pipeline processing a series of data from clusters of galaxies observations. Indeed, the manual operation of XMMSAS (using either external correction routines or integrated algorithms from versions `5.3` and later) is extremely time-consuming: several days for each dataset.

Thanks to our pipeline, we reduced the required time down to less than $`10h`$ per dataset, adding a working day for eye-inspection of the full sample results and another for compilation and final manual calculations.

This pipeline integrates the latest calibration data from the EPIC consortium as well as specific fitting routines written in IDL and XSPEC environments. It is fully configurable to allow eye-inspection and partial or total re-runs if needed. It is extensively detailed by \[Marty et al., 2003\], but useful information about “blank” and “dark” fields (for background estimation and subtraction) are also respectively given by \[Lumb et al., 2002, Marty et al., 2002\]. The general processing scheme is the following:

| $``$ | re-process each cluster raw data to account for latest calibration knowledge; |
| --- | --- |
| $``$ | re-project “blank” and/or “dark” data onto telescope aspect solution used for each cluster; |
| $``$ | extract valid events within good-time-intervals (to reject noise and proton flares); |
| $``$ | detect and mask out spurious components (field sources and remaining bad pixels); |
| $``$ | weight the events according to calibration data; |
| $``$ | find the best cluster centroid; |
| $``$ | process radial brightness profiles and fit a standard $`\beta `$-model convolved by the instrument PSF; |
| $``$ | process full maps in physical flux units and adaptively smooth according to PSF description; |
| $``$ | process spectra and call XSPEC for fitting (global for the mean temperature, annuli for the profile); |
| $``$ | compute hardness maps and derive another temperature profile. |

## THE PIPELINE RESULTS

### Surface brightness profiles

The core radius $`R_c`$, power index $`\beta `$, normalization $`f_0`$ and background constant $`BKG`$ are coming from the standard $`\beta `$-model fitting. Note that the surface brightness profiles have been extracted between $`0.3`$ and $`6.3keV`$, logarithmically rebinned to enhance their signal-to-noise ratios and that the inner $`40arcsec`$ of each cluster profile have been ignored for the fit so as to avoid fitting any emission excess, associated for example to a cooling flow. The background subtraction concerned only the instrumental component so that the $`BKG`$ constant is a measure of the sky X-ray emission background. The $`\beta `$-models are convolved by the instrument PSF before fitting.

The $`R_{1\sigma }`$ parameter is a measure of the radius inside which the cluster signal-to-noise ratio is greater than $`1`$. The $`F_X(data)`$ and $`F_X(model)`$ parameters then are the integrated fluxes within this detection limit radius, from the data and from the best fit $`\beta `$-model; in each case, these fluxes are cleaned from the sky background component; the difference between them thus shows the contribution of any excess emission.

### Surface brightness maps

Surface brightness maps that have been produced in the same energy band during the pipeline analysis are not presented herein, because of room restrictions and because they are not necessary for the following calculations. However, they may be found on the internet111Only compiled on ftp://www-station.ias.u-psud.fr/pub/epic/HiLx/2002Marty\_santiago\_poster.jpg at the moment, but individual maps will soon be available in the same ftp directory. if so desired.

### Spectral analysis

Spectra are extracted in the $`0.2`$ to $`12.0keV`$ band and rebinned according the Poisson statistics requirements (at least $`25`$ net counts per bin). They are background subtracted from both instrumental components (estimated using the “closed” data) and sky components (from a spectrum extracted in a peripheral annulus of each dataset, between the $`600`$ and $`800arcsec`$ radii) and then fitted against the $`TBabs\times MEKAL`$ photo-absorption and bremsstrahlung models.

The redshifts $`z`$ and hydrogen column densities $`N_H`$ have been respectively fixed to the same optical redshift references as used in the original XBACS database, and to the weighted means of the $`N_H`$ measurements from the \[Dickey & Lockman, 1990\] database included within $`1\mathrm{deg}`$ of the cluster center. The resulting temperatures $`kT`$, abundance relative to solar composition $`\mu `$ and normalization $`S_0`$ are then reported. These measurements come from spectra extracted inside a circle of $`400arcsec`$ radius, and hence integrate any central excess emission; that is why temperatures may appear a bit low.

### Direct temperature profiles

Spectra have also been extracted within $`40arcsec`$-wide concentric annuli following the same algorithm than for global spectra. Fitting them using the same XSPEC settings hence led to a direct measure of the temperature profile of each cluster. However, no PSF correction have been applied from one annulus to another thus probably leading to a general underestimation in the fitted temperatures, specially toward the center (\[Markevitch, 2002\]).

A last spectrum has been extracted and fitted, in the same region than that used for the $`\beta `$-model fitting, i.e. an annulus starting at $`40arcsec`$ and extending upto the arbitrary limit of $`400arcsec`$. This allows to compare between the global temperature with the one found outside the core region, which may host a cooling flow.

### Intrinsic parameters calculation

The luminosity distance $`D_L`$ and $`1arcsec`$ equivalent length $`R_{1\mathrm{"}}`$ at the cluster position are deduced from the cluster redshift $`z`$ and the default cosmology ($`\mathrm{\Omega }_0=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ and $`H_0=50h_{50}km/s/Mpc`$).

The virial radius $`R_V`$ has been identified with the distance $`r_{200}`$ from the cluster center where the mean enclosed overdensity equals $`200`$, as calibrated upon the measured X-ray temperature by \[Evrard et al., 1996\].

The X-ray luminosity $`L_X`$ was integrated between our XMM-Newton energy bandwith ($`0.36.3keV`$) within the virial radius defined above. We applied a numerical bolometric correction to deduce the absolute luminosity $`L_X(bol)`$.

The galaxies velocities dispersions $`\sigma _{wl_D}`$ and $`\sigma _{wl_B}`$ have been reprinted from two weak lensing analysis, respectively from NOT/ALFOSC and UH8k detectors (\[Dahle et al., 1990\]), and CFHT12k instrument (\[Bardeau & Kneib, 2003\]). The corresponding $`\beta `$ parameters were computed according to the relationship with the measured X-ray temperature in the isothermal hydrostatic equilibrium case (\[Cavaliere & Fusco-Femiano, 1976\]).

The electronic density $`n_e(0)`$ at the center of the cluster has been estimated using the standard isothermal $`\beta `$-model equations (\[Cavaliere & Fusco-Femiano, 1976\]) and the previously measured X-ray temperature and luminosity. By integration within the virial radius defined above, we also estimated the cluster gas mass. Finally, the total virial mass has been deduced from two different calibrated relations upon the measured X-ray temperature ($`M_B`$ and $`M_T`$ resp. from \[Cavaliere & Fusco-Femiano, 1976, Arnaud & Evrard, 1999\]) as well as from the virial theorem relation upon the galaxies velocities dispersion ($`M_V`$ from the $`\sigma _{wl_B}`$ defined above, except for A773 for which we used $`\sigma _{wl_D}`$). The gas fraction $`f_g`$ simply is the ratio of the gas mass to the total mass (we chose $`M_B`$ as the isothermal reference).

## DISCUSSION

All analysed clusters show an excess in central brightness that seems correlated with the presence of a very bright central galaxy in optical images. They hence are all candidates for cooling flow hosts, but A267 and A383 temperature profiles do not exhibit significant cooling cores. In addition, all clusters but A1835 tend to follow the universal temperature profile calibrated by \[Loken et al., 2002\] (Figure 3), radii being rescaled either to core radius or virial radius natural units, with a declining slope toward the edges. They are also very close to the \[DeGrandi & Molendi, 2002\] polytropic interpretation of their Beppo-SAX observed clusters data, which show slopes more horizontal at high radii for those hosting a cooling flow than for non-cooling flow specimens, although our sample does not present a steeply cooling core in average.

If the problem of the cooling core may be bypassed by using only outer regions for temperature and $`\beta `$-model fitting, these declining profile still questions the hypothesis of isothermality followed throughout this first analysis. Nevertheless, we expect from the general good agreement of the different mass (Figure 3) and $`\beta `$ parameters estimations that our optical (CFHT12k) and X-rays (XMM-Newton) measurements are still both quite consistent with the isothermal model. But we note also that all calculations rely on a virial radius that has been estimated using a relation based on the isothermal hypothesis; that illustrates the need for another independant estimator of this radius. Even the clue given by the temperature profiles (Figure 3), that seem to indicate that the virial radius is about $`10`$ times the core radius, is not so helpful if models fitted to either optical or X-ray data still rely upon isothermality.

We plotted a temperature-luminosity logarithmic diagram (Figure 3) for our subsample and compared it to the model proposed by \[Arnaud & Evrard, 1999\]. We observe that five clusters are indeed well co-aligned but following a shallower slope. A209 seems to depart from this alignment, as well as A1835 in a very significant manner. We also compared a temperature-velocity dispersion diagram with the model from \[Girardi et al., 1996\] and found it to be very close to the lower $`90\%`$ confidence level boundary, with again the exception of A1835.

Finally, the temperature-mass diagram is remarkably consistent with a mean power index of $`1.576`$, while the gas mass fraction seems to average about $`10\%`$ except again for A1835, as well as A209 and A773.

## CONCLUSIONS & PERSPECTIVES

While early XMM-Newton observations (\[Arnaud et al., 2001, Majerowicz et al., 2002\]) did not confirm declining temperature profiles as seen by late ASCA observations (\[Markevitch et al., 1998\]) and Beppo-SAX (\[DeGrandi & Molendi, 2002, Loken et al., 2002\]), we find in the present analysis of 7 XBACS clusters that 6 of them do follow this kind of polytropic model. The seventh, A1835, systematically departs from the average in all diagrams (temperature-mass, temperature-luminosity), which confirms the rather complicated nature of its internal dynamics (\[Majerowicz et al., 2002\]).

However, we still are working toward the comparison of these emission weighted temperature profiles, that have been obtained without correction for PSF and may thus have been underestimated (\[Markevitch, 2002\]), with profiles built on the basis of hardness ratio maps, that may be adaptively handled to account for the PSF and binned between concentric isophotes rather than concentric annulii to adapt to non-spherical morphologies. Also, more work is intended so as to derive accurate error bars on each diagram, and attempt to estimate the virial radii from independent optical datasets.

Finally, this work has been made possible through the use of a pipeline analysis environment, specifically aimed at dealing with extended sources data from XMM-Newton instruments. It has confirmed the importance of maintenance and upgrades, to keep up-to-date with the instruments calibration database, as well as the need for such automatic tools in the purpose of building catalogs from the latest X-ray telescopes, so that massive cluster investigation may be possible in the future: access to large statistical samples, cross-compare mass and morphological measurements with optical data, correlate X-ray and Sunyaev-Zel’Dovich data…The ESA’s Planck satellite is promising thousands of millimetric clusters detections by the end of the decade, which could lead to narrower constraints on cosmology, provided an extensive X-rays database, like BAX (\[Sadat et al., 2002\]), also exists.

## ACKNOWLEDGEMENTS

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the NASA. This works also benefits from observations conducted with the CFH12k camera at the Canada France Hawaii Telescope.

The authors wish to thank Bruno Altieri, Christian Erd, David Lumb, Richard Saxton, as well as the whole “epic-cal” team for their help. Our thanks also go to Pasquale Mazzotta for fruitful discussions and comparisons between results from Chandra and XMM-Newton observatories.

J.P.K. acknowledges support from an ATIP/CNRS grant and from PNC. I.R.S. acknowledges support from the Royal Society and Leverhulme Trust.

E-mail address of P.B. Marty: philippe.marty@ias.u-psud.fr