text:
# Effect of point defects on heat capacity of yttria-stabilized zirconia

## Abstract

First-principles calculation and anharmonic dynamical theory were used in sequence to explain a large excess heat capacity observed in yttria-stabilized zirconia in comparison with the additive rule value. It is found that the excessive shape of heat capacity decays gradually with the Y2O3 doping when the number of environmentally different O sites falls to its zero value at 33 mol $`\%`$ Y2O3-ZrO2 due to the Y atoms adjacent. The model and results presented in this work provide a key insight into the complex behaviour and characterization of fast-ion conductors.

Point lattice defects both native and artificial, such as vacancies, self-interstitials, solute and substitute atoms, strongly affect the most fundamental properties of materials.Flynn  In yttria-stabilized zirconia (YSZ), the material with numerous commercial applications Subbarao , formed by the addition of Y2O3 to ZrO2, the trivalent dopant Y3+ substitute for some of the host cations and, in order to maintain charge neutrality, one O vacancy (V0) must be created for each pair of dopant cations. Y2O3 is used to stabilize the tetragonal ($`t`$) phase of (ZrO2)100-x(Y2O3)x over the composition range 2$`<x<`$9 mol $`\%`$ and the cubic ($`c`$) fluorite phase with 4$`<x<`$40 mol $`\%`$. The presence of relaxed V0 and Y substitute atoms make the local atomic environments of YSZ rather different from the stoichiometric high-temperature $`t`$\- and $`c`$-polymorphs of pure ZrO2 whose 8-fold cations are distorted and perfectly coordinated, respectively. In YSZ, the average cation coordination number, ranged between 7 and 8, is reduced gradually with increasing Y2O3. If the V0 associates with Zr ions, as the X-ray absorption findings Li  suggest, it may support a coordination-driven ordering model of stabilization.Stefanovich  Placing Y in the next nearest neighbour (NNN) cation positions allows the coordination of Zr in the nearest neighbours (NN) sites to the V0 to be similar to the monoclinic ZrO2 arrangement, with Zr located in 7-fold coordination environments, while Y remains 8-fold coordinated.

The electron energy-loss near-edge structure (ELNES) of YSZ demonstrates the features of the experimental O $`K`$-edge shapes which depend on the crystal structures and the Y2O3 composition.OCMVAFP  In modelling the most important and widely used properties of YSZ the relaxation mechanism is essential. The effects of relaxation were treated within the plane-wave, pseudopotential based free energy molecular dymanics (FEMD) technique.Ali  Using the relaxed configurations (RC) of point defects, the ELNES calculations were carried out.OCMVAFP  Since the results of calculation show very good agreement with the experimental O $`K`$-edge signal, it seems that theory reflects the realities of relaxation. There is much yet to be learned about the fundamental nature of transport processes in YSZ. Thus, our motivation is to make a contribution to the microscopic modelling and tailoring the heat capacity of YSZ.

A large excess heat capacity (EHC) has been observed heat\_cap  in YSZ at low and room temperatures for 7.8, 9.7 and 11.4 mol $`\%`$ Y2O3 compared to the heat capacity calculated from those of pure ZrO2 and Y2O3 by the additivity rule: $`\mathrm{\Delta }C_p`$ = $`C_p^{YSZ}`$ \- \[(1-$`x`$) $`C_p^{ZrO_2}`$ \+ $`x`$ $`C_p^{Y_2O_3}`$\]. The shape of $`\mathrm{\Delta }C_p`$($`T`$) is very similar to that of the Schottky disorder. The same behaviour of heat capacity may yield the two-level system (TLS), quite a simple model, used to simulate the low-temperature properties of amorphous insulators: $`C_p^{TLS}`$=$`k_B(\mathrm{\Delta }/T)^2`$ $`e^{\mathrm{\Delta }/T}`$/($`e^{\mathrm{\Delta }/T}`$+1)2, where $`\mathrm{\Delta }`$ is the energy difference between the levels. $`C_p^{TLS}`$ shows a maximum at temperature $`T_m`$ 0.42 $`\mathrm{\Delta }`$ when the filling-rate variation for the occupied TLS states reaches its maximum value. At $`TT_m`$, all TLS states are occupied equally and therefore heat capacity decreases as $`C_p^{TLS}1/T^2`$.

It appears that for an anharmonic oscillator (AO), which dynamics is defined by the double-well potential $`U`$($`x`$)=$`x^2`$+$`ax^3`$+$`bx^4`$, ($`a,b>`$0) with a maximum at $`x_0`$=0 and two minima at $`x_1`$ and $`x_2`$, separated by the energy $`\mathrm{\Delta }`$=$`U(x_2)U(x_1)>`$0, the temperature dependence of heat capacity is similar to that of TLS. The partition function $`Z_{AO}`$=$`\frac{1}{2\pi h}`$ $`𝑑x𝑑vexp[\frac{mv^2/2+U(x)}{k_BT}]`$ is calculated assuming $`U_0U_i`$. If $`U`$ is approximated near its minima by the second-degree polynomials, $`U_1(xx_1)^2\omega _1^2/2`$ and $`U_2(xx_1)^2\omega _2^2/2+\mathrm{\Delta }`$ where $`\omega _i`$ are the frequencies of small vibrations at $`x_i`$ then, using the general thermodynamic relations, one can obtain $`C_p^{AO}`$ = $`k_B`$$`\frac{\mathrm{\Delta }^2}{T^2}`$$`\frac{exp(\mathrm{\Delta }/T)}{exp(\mathrm{\Delta }/T+\omega _1/\omega _2)^2}`$ \+ $`k_B`$. In the latter, the temperature dependence of the first term, $`\mathrm{\Delta }C`$($`T`$), which defines the difference in heat capacity between anharmonic and harmonic oscillators, is close to $`C_p^{TLM}`$ and finally, at $`\omega _1/\omega _2`$=1, $`\mathrm{\Delta }CC_p^{TLM}`$. Hence, the EHC effect may appear as the result of strong anharmonic vibration modes either localized or delocalized ones.

The fact that at $`x>`$8 mol $`\%`$ Y2O3 $`\mathrm{\Delta }C_p`$ decreases heat\_cap  with increasing the dopant concentration may suggest that the delocalized anharmonic vibrations cause the EHC effect. It is known that the main contribution to the structural transformations comes from particular vibrational modes. In pure ZrO2, a zone-boundary soft phonon, $`X_2^{}`$, which breaks the $`c`$-symmetry of the O sublattice, displacing the O atoms toward their positions in the $`t`$-phase, may be responsible for the $`ct`$ phase transformation. In YSZ, the experimental data Exp-neutron  are not fully clarified for this soft mode because of the static disorder in the O sublattice. Using the RC of YSZ, the $`X_2^{}`$-like phonon was calculated OS2002  by means of the frozen-phonon method for each composition of (ZrO2)100-x(Y2O3)x between 3 and 10 mol $`\%`$ Y2O3. For 10 mol $`\%`$ Y2O3, the effective potential has a single minimum and, with decreasing Y2O3 content, the potential develops two minima. The temperature dependence of the phonon frequency, calculated within the modified pseudoharmonic approximation, quantifies accurately the transition temperature above which the $`c`$-phase is stabilized.OS2002  In the 3 mol $`\%`$ Y2O3 case, it may start around room temperature. The upper boundary of the single $`t`$-phase field in the experimental phase diagram Miller81  shows a similar rapid drop of temperature with Y2O3 concentration. However, the $`X_2^{}`$-like soft phonon of YSZ, in another words, the delocalized vibrations can not be responsible for the EHC effect because of the potential symmetry.

Regarding anharmonic effects, the atomic vibrations localized near the V0 should be considered as well. YSZ is a solid electrolyte in which ionic transport takes place by anions moving among their positions by the vacancy diffusion mechanism. If any two V0 sites, which an O atom occupies before and after a diffusive jump, are non-equivalent then theory enables the double-well asymmetric potential profile. In such a study of diffusion, the molecular dynamics technique can fruitfully be used. The principal output of the dynamic computation is the set of atomic trajectories from which all other properties are calculated. Molecular dynamics has been used in the past by many authors to calculate the potential energy dependant on the position of diffusion atom.Flynn  Even if one could calculate the induced trajectories through the saddle-point the results were strongly affected by the reaction coordinates of all other atoms. The benefits and limitations of this molecular dynamics approach have been discussed and it is generally accepted now that using the large unit cells this scheme is very complicated for precise calculation whereas the results, obtained in such a way, are not always satisfactory.Bennett

The features of the computed $`U`$ profile were demonstrated using the single-$`V_0`$ 11-atom Zr2Y2O7 cell, which contains two formula units of ZrO2 and one formula unit of Y2O3, that corresponds to a composition of 33 mol $`\%`$ Y2O3.OCMc+2000  The RC, shown in the (a) panel of Fig. 1, was used in a further static optimization using the displacements of O-2 along the $`z`$-direction while all other atoms are fixed at their relaxed positions. In this case, labeled in the panel (b) of Fig. 1 as a “solid pass” trajectory, the O-2 position in the $`z`$=0 plane corresponds to the saddle point with the potential barrier $`E_b`$ of $``$0.2 eV/cell. If Zr in the $`z`$=0 plane were allowed to relax when the O-2 is moved, as shown in Fig. 1 as the “soft pass” trajectory, then the equlibrium position of O-2 is much closer to the saddle point while $`E_b`$ becomes significantly less.

In Zr2Y2O7, where all O sites are equivalent to each other, the Zr/Y environment enables the symmetric form of $`U(x)`$ and therefore can not result in the EHC effect. The saddle-point neighorhood plays a key role on the instability-barrier formation, which is very sensitive to the model parameters. Thus, the dynamics theory of ion transport is still very much a semiempirical science and further simplification, based on the realities for the activation energy, should be made.

Atomistic RC of the (96-$`y`$)-atom supercells ($`y`$=1,2,…) allow to model (ZrO2)100-x(Y2O3)x. In the 95-atom cell, which corresponds to $``$3 mol $`\%`$ Y2O3, the RC with both Y atoms in the NNN shell was obtained.OCMVAFP  Fig. 2 shows the RC energies, expressed as an equivalent temperature per 12-atom cell, to give an idea of the relative stability and temperature required. In the 94-atom cell, the RC with two V0s along the $``$111$``$ fluorite direction results in a 6-fold coordinated and six 7-fold coordinated NN Zr atoms. This result supports the interpretation of the neutron data by Goff et al. Goff99 .

Electrostatic considerations suggest that V0s should repel. From this point of view, the tri-V0 RC, with all V0 separated by $`\sqrt{5}`$/2 from each other (in units of the $`c`$-ZrO2 lattice constant $`a`$), is rather reasonable. The repulsion tendencies between vacancies consistent with both electrostatics and the Zr-coordination lowering model. Since the shortest distance of $`a`$/2 between the V0 sites is not energetically favorable in YSZ at low Y2O3 concentrations we believe that each V0 has six NN O atoms each of which can jump into the empty site.

The intrinsic process of such chaotic O moves should occur in YSZ at finite temperatures via the vacancy mechanism. In particular, when an apical O in 3 mol $`\%`$ Y2O3 moves into V0 the RC with the Y NN and NNN atoms is created. The difference $`\mathrm{\Delta }`$=210 K between these two RC was used to construct the potential profile, as shown in the lower diagram of Fig. 2. If V0 moves consequently into the next O site of this anion tunnel along the $`z`$ axis, the Y NN-NNNN configuration is realized. We take $`E_b=E_a`$1 eV using the activation energy $`E_a`$ of YSZ, obtained from the Arrhenius relation for the ionic dc conductivity Arrhenius  $`\sigma \omega _0d^2exp(E_a/k_BT)`$, where $`\omega _0`$ is the attempt frequency. To form the effective potential, characterized by $`E_b`$, $`\mathrm{\Delta }`$ and the distance moved during a single jump $`d`$=$`x_1x_2`$, the three-parameter polynomial $`U(x)`$= $`ax^2+bx^3+cx^4`$ was used. Since $`E_b\mathrm{\Delta }`$, the O vibrations near the $`U`$ minima may be estimated as $`\omega \sqrt{\frac{32E_b}{md^2}}`$, where $`m`$ is the O mass. As a result, $`\omega `$ 5$`\times `$1014Hz is close to that found experimentally Ohta+2002  showing very reasonable model approximation.

In Fig. 3a the calculated EHC shape $`\mathrm{\Delta }C_p`$($`T`$) = $`𝒩^{TLS}C_p^{TLS}`$ is plotted, with its dependence on temperature and dopant concentration. Here $`𝒩^{TLS}`$ was obtained in such a way that each maximum coincides with the experimental EHC amount heat\_cap  at 7.76, 9.7 and 11.3 mol $`\%`$ Y2O3. The EHC observations of 7.76 mol $`\%`$ Y2O3 are shown in Fig. 3a to illustrate the similarity. All theoretical $`\mathrm{\Delta }C_p`$ shapes, with a maximum at 88 K, have been aligned with the experimental peak at $`T_m`$=75 K using the $`T/\mathrm{\Delta }`$ dimensionless scale.

The (b) panel of Fig. 3 shows that the $`𝒩^{TLS}`$/$`𝒩_y`$ ratio ($`𝒩_y`$ is the number of $`V_0`$) falls quickly with increase in Y2O3 content and, finally, $`𝒩^{TLS}`$/$`𝒩_y`$0 at $`x`$=33 mol $`\%`$ when the $`U`$ symmetry for diffusion atom should appear. At low concentrations of Y2O3, not all possible diffusion jumps lead to the asymmetric $`U`$. For example, in 3 mol $`\%`$ Y2O3, $`𝒩_y`$=1 whereas $`𝒩^{TLS}`$ falls from six to four since two O NN sites to V0 are energetically equivalent due to the Y atoms adjacent. The number of such environmentally equivalent O sites, which not contribute to the EHC effect, increases with increasing the dopant concentration.

The high-probability O hopping via $`V_0`$ between the low-energy equivalent sites, localized close up to each other, can be considered as some sort of the vacancy block, which results in the anisotropic absorption peak to the internal frictions measured in 9.5 mol $`\%`$ Y2O3.Ohta+2002  This observation is in agreement with our findings confirming earlier suggestions that some vacancies exist in bound states as di-$`V_0`$ along the  fluorite direction. As the number of $`V_0`$-blocks increases with Y2O3 doping, the peak magnitude for the localized relaxation may not scale simply with the dopant concentration, reaching a maximum at some critical concentration $`x_m`$ correlated with that of the EHC maximum. At $`x>x_m`$, the number of sites on the O sublattice with similar local environment increases significantly. This may lead to the long-range transport of O ions via vacancies reducing therefore the local hopping process in $`V_0`$-blocks. Regarding the concentration profile of localized relaxation explained, such observations at $`x>x_m`$ may be masked by overlap between the isotropic and anisotropic absorption peaks due to diffusion and $`V_0`$-block relaxation, respectively.

In summary, our study of YSZ being performed using the ab initio relaxed configurations and dynamical model of O hopping allows to explain the EHC anomaly observed with the Y2O3 doping. In YSZ, as it has been demonstrated, the anisotropic absorption peak of internal frictions and EHC effect have a common vacancy-block nature showing great potential for theory to deal with the realistic model of doped ceramics. Because of the general physics involved we believe that similar excessive behaviour of heat capacity might be observed in a wide class of ionic conductors where anion hopping goes through non-equivalent positions.

E.S. thanks Russian Basic Research Foundation (Grants No. 00-02-17426, 01-02-96462) for financial support.


format: text

text:
# STATISTICS OF COMPLEX DISSIPATIVE SYSTEMS ”Space Time Chaos: Characterization, Control and Synchronization”; June 19, 2000 - June 23, 2000; Department of Physics and Applied Mathematics and Institute of Physics, University of Navarra, Pamplona, Navarra, Spain

## 1 Summary

A complex quantum system with energy dissipation is considered. The quantum Hamiltonians $`H`$ belong the complex Ginibre ensemble. The complex-valued eigenenergies $`Z_i`$ are random variables. The second differences $`\mathrm{\Delta }^1Z_i`$ are also complex-valued random variables. The second differences have their real and imaginary parts and also radii (moduli) and main arguments (angles). For $`N`$=3 dimensional Ginibre ensemble the distributions of above random variables are provided whereas for generic $`N`$\- dimensional Ginibre ensemble second difference’s, radius’s and angle’s distributions are analytically calculated. The law of homogenization of eigenergies is formulated. The analogy of Wigner and Dyson of Coulomb gas of electric charges is studied. The stabilisation of system of electric charges is dealt with.

## 2 Introduction

We study generic quantum statistical systems with energy dissipation. The quantum Hamiltonian operator $`H`$ is in given basis of Hilbert’s space a matrix with random elements $`H_{ij}`$ . The Hamiltonian $`H`$ is not Hermitean operator, thus its eigenenergies $`Z_i`$ are complex-valued random variables. We assume that distribution of $`H_{ij}`$ is governed by Ginibre ensemble . $`H`$ belongs to general linear Lie group GL($`N`$, C), where $`N`$ is dimension and C is complex numbers field. Since $`H`$ is not Hermitean, therefore quantum system is dissipative system. Ginibre ensemble of random matrices is one of many Gaussian Random Matrix ensembles GRME. The above approach is an example of Random Matrix theory RMT . The other RMT ensembles are for example Gaussian orthogonal ensemble GOE, unitary GUE, symplectic GSE, as well as circular ensembles: orthogonal COE, unitary CUE, and symplectic CSE. The distributions of the eigenenergies $`Z_1,\mathrm{},Z_N`$ for $`N\times N`$ Hamiltonian matrices is given by Jean Ginibre’s formula :

$`P(z_1,\mathrm{},z_N)=`$ (1)
$`={\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{\pi j!}}{\displaystyle \underset{i<j}{\overset{N}{}}}|z_iz_j|^2\mathrm{exp}({\displaystyle \underset{j=1}{\overset{N}{}}}|z_j|^2),`$
where $`z_i`$ are complex-valued sample points ($`z_i𝐂`$). For Ginibre ensemble we define complex-valued spacings $`\mathrm{\Delta }^1Z_i`$ and second differences $`\mathrm{\Delta }^2Z_i`$:

$$\mathrm{\Delta }^1Z_i=Z_{i+1}Z_i,i=1,\mathrm{},(N1),$$
(2)
$$\mathrm{\Delta }^2Z_i=Z_{i+2}2Z_{i+1}+Z_i,i=1,\mathrm{},(N2).$$
(3)
The $`\mathrm{\Delta }^2Z_i`$ are extensions of real-valued second differences

$$\mathrm{\Delta }^2E_i=E_{i+2}2E_{i+1}+E_i,i=1,\mathrm{},(N2),$$
(4)
of adjacent ordered increasingly real-valued energies $`E_i`$ defined for GOE, GUE, GSE, and Poisson ensemble PE (where Poisson ensemble is composed of uncorrelated randomly distributed eigenenergies) .

There is an analogy of Coulomb gas of unit electric charges pointed out by Eugene Wigner and Freeman Dyson. A Coulomb gas of $`N`$ unit charges moving on complex plane (Gauss’s plane) C is considered. The vectors of positions of charges are $`z_i`$ and potential energy of the system is:

$$U(z_1,\mathrm{},z_N)=\underset{i<j}{}\mathrm{ln}|z_iz_j|+\frac{1}{2}\underset{i}{}|z_i^2|.$$
(5)
If gas is in thermodynamical equilibrium at temperature $`T=\frac{1}{2k_B}`$ ($`\beta =\frac{1}{k_BT}=2`$, $`k_B`$ is Boltzmann’s constant), then probability density function of vectors of positions is $`P(z_1,\mathrm{},z_N)`$ Eq. (1). Complex eigenenergies $`Z_i`$ of quantum system are analogous to vectors of positions of charges of Coulomb gas. Moreover, complex-valued spacings $`\mathrm{\Delta }^1Z_i`$ are analogous to vectors of relative positions of electric charges. Finally, complex-valued second differences $`\mathrm{\Delta }^2Z_i`$ are analogous to vectors of relative positions of vectors of relative positions of electric charges.

The $`\mathrm{\Delta }^2Z_i`$ have their real parts $`\mathrm{Re}\mathrm{\Delta }^2Z_i`$, and imaginary parts $`\mathrm{Im}\mathrm{\Delta }^2Z_i`$, as well as radii (moduli) $`|\mathrm{\Delta }^2Z_i|`$, and main arguments (angles) $`\mathrm{Arg}\mathrm{\Delta }^2Z_i`$.

## 3 Second Difference Distributions

We define following random variables for $`N`$=3 dimensional Ginibre ensemble:

$$Y_1=\mathrm{\Delta }^2Z_1,A_1=\mathrm{Re}Y_1,B_1=\mathrm{Im}Y_1,$$
(6)
$$R_1=|Y_1|,\mathrm{\Phi }_1=\mathrm{Arg}Y_1,$$
(7)
and for the generic $`N`$-dimensional Ginibre ensemble :

$$W_1=\mathrm{\Delta }^2Z_1,P_1=|W_1|,\mathrm{\Psi }_1=\mathrm{Arg}W_1.$$
(8)
Their distributions for are given by following formulae respectively :

$`f_{Y_1}(y_1)=f_{(A_1,B_1)}(a_1,b_1)=`$ (9)
$`={\displaystyle \frac{1}{576\pi }}[(a_1^2+b_1^2)^2+24]\mathrm{exp}({\displaystyle \frac{1}{6}}(a_1^2+a_2^2)),`$
is second difference distribution for 3-dimensional Ginibre ensemble.

$$f_{A_1}(a_1)=\frac{\sqrt{6}}{576\sqrt{\pi }}(a_1^4+6a_1^2+51)\mathrm{exp}(\frac{1}{6}a_1^2),$$
(10)
$$f_{B_1}(b_1)=\frac{\sqrt{6}}{576\sqrt{\pi }}(b_1^4+6b_1^2+51)\mathrm{exp}(\frac{1}{6}b_1^2),$$
(11)
$`f_{R_1}(r_1)=`$ (12)
$`\mathrm{\Theta }(r_1){\displaystyle \frac{1}{288}}r_1(r_1^4+24)\mathrm{exp}({\displaystyle \frac{1}{6}}r_1^2),`$
$`f_{\mathrm{\Phi }_1}(\varphi _1)={\displaystyle \frac{1}{2\pi }},\varphi _1[0,2\pi ],`$
are real part’s, imaginary part’s, modulus’s and angle’s probability distributions for 3-dimensional Ginibre ensemble, where $`\mathrm{\Theta }(r_1)`$ is Heaviside (step) function. Notice, that the angle $`\mathrm{\Phi }_1`$ has uniform distribution. Next, second difference’s distribution for $`N`$-dimensional Ginibre ensemble reads $`P_3(w_1)`$ :

$`P_3(w_1)=`$ (13)
$`=\pi ^3{\displaystyle \underset{j_1=0}{\overset{N1}{}}}{\displaystyle \underset{j_2=0}{\overset{N1}{}}}{\displaystyle \underset{j_3=0}{\overset{N1}{}}}{\displaystyle \frac{1}{j_1!j_2!j_3!}}I_{j_1j_2j_3}(w_1),`$
$`I_{j_1j_2j_3}(w_1)=`$ (14)
$`=2^{2j_2}{\displaystyle \frac{^{j_1+j_2+j_3}}{^{j_1}\lambda _1^{j_2}\lambda _2^{j_3}\lambda _3}}F(w_1,\lambda _1,\lambda _2,\lambda _3)|_{\lambda _i=0},`$
$`F(w_1,\lambda _1,\lambda _2,\lambda _3)=`$ (15)
$`=A(\lambda _1,\lambda _2,\lambda _3)\mathrm{exp}[B(\lambda _1,\lambda _2,\lambda _3)|w_1|^2],`$
$`A(\lambda _1,\lambda _2,\lambda _3)=`$ (16)
$`={\displaystyle \frac{(2\pi )^2}{(\lambda _1+\lambda _2\frac{5}{4})(\lambda _1+\lambda _3\frac{5}{4})(\lambda _11)^2}},`$
$`B(\lambda _1,\lambda _2,\lambda _3)=`$ (17)
$`=(\lambda _11){\displaystyle \frac{2\lambda _1\lambda _2\lambda _3+\frac{1}{2}}{2\lambda _1+\lambda _2+\lambda _3\frac{9}{2}}}.`$
Finally,

$`f_{P_1}(\varrho _1)=2\pi P_3(\varrho _1)=`$ (18)
$`=2\pi ^2{\displaystyle \underset{j_1=0}{\overset{N1}{}}}{\displaystyle \underset{j_2=0}{\overset{N1}{}}}{\displaystyle \underset{j_3=0}{\overset{N1}{}}}{\displaystyle \frac{1}{j_1!j_2!j_3!}}I_{j_1j_2j_3}(\varrho _1),`$
$`f_{\mathrm{\Psi }_1}(\psi _1)={\displaystyle \frac{1}{2\pi }},\psi _1[0,2\pi ],`$
are modulus’s and angle’s probability distributions for $`N`$-dimensional Ginibre ensemble. Notice, that again the angle $`\mathrm{\Psi }_1`$ has uniform distribution.

## 4 Conclusions

We compare second difference distributions for different ensembles by defining following dimensionless second differences:

$$C_\beta =\frac{\mathrm{\Delta }^2E_1}{<S_\beta >},$$
(19)
$$X_1=\frac{A_1}{<R_1>},$$
(20)
where $`<S_\beta >`$ are the mean values of spacings for GOE(3) ($`\beta =1`$), for GUE(3) ($`\beta =2`$), for GSE(3) ($`\beta =4`$), for PE ($`\beta =0`$) , and $`<R_1>`$ is mean value of radius $`R_1`$ for $`N`$=3 dimensional Ginibre ensemble .

On the basis of comparison of results for Gaussian ensembles, Poisson ensemble, and Ginibre ensemble we formulate :

Homogenization Law: Random eigenenergies of statistical systems governed by Hamiltonians belonging to Gaussian orthogonal ensemble, to Gaussian unitary ensemble, to Gaussian symplectic ensemble, to Poisson ensemble, and finally to Ginibre ensemble tend to be homogeneously distributed.

It can be restated mathematically as follows:

If $`H`$ GOE, GUE, GSE, PE, Ginibre ensemble, then Prob($`D`$)=max, where D is random event corresponding to vanishing of second difference’s probability distributions at the origin.

Both of above formulation follow from the fact that the second differences’ distributions assume global maxima at origin for above ensembles . For Coulomb gas’s analogy the vectors of relative positions of vectors of relative positions of charges statistically most probably vanish. It means that the vectors of relative positions tend to be equal to each other. Thus, the relative distances of electric charges are most probably equal. We call such situation stabilisation of structure of system of electric charges on complex plane.

## Acknowledgements

It is my pleasure to most deeply thank Professor Jakub Zakrzewski for formulating the problem. I also thank Professor Antoni Ostoja-Gajewski for creating optimal environment for scientific work.

## References


text:
# Dust emissivity in the Submm/Mm

## 1 Introduction

Despite being a fundamental parameter in Far Infrared, Submillimetre and Millimetre astronomy, few measurements of the dust emissivity111The emissivity proper is the emission \[absorption\] cross section normalized to the geometrical cross section of a dust grain. Other authors prefer the cross section normalized to the mass of the grain. Both quantities are the same when normalised to the analogous quantity for V-band extinction. are available (for a review, see Alton et al. 2000; James et al. 2002). The submm/mm dust emissivity is particularly important for star formation studies. Since molecules are known to deplete inside prestellar cores (Bergin et al. 2002), dust emission may represent the best tracer of the gas density distribution just prior to the onset of gravitational collapse. Thus measurements in the submm/mm define the initial conditions from which a core collapses to form a star. Alternatively, the density distribution can be mapped through extinction, by measuring the near-infrared colour excess towards giant stars in the background of a cloud (Lada et al. 1994). High resolution and S/N maps can be obtained for object in the foreground of dense stellar fields (Alves et al. 2001).

This is the case for the dark cloud Barnard 68, a starless globule seen in the foreground of the Galactic Bulge. Alves et al. (2001) have produced a high resolution extinction map of the cloud, measuring the H-K colour excess of nearly 4000 stars in its background. In this paper, we compare the extinction map with observations of submm/mm dust emission of similar resolution. Observations are described in Sect. 2. In Sect. 3 we derive the dust emissivity from the correlation between emission and extinction (in a way similar to Kramer et al. 1998, 2002). We will also adopt a temperature gradient within the cloud, which was derived from a model of dust heating. Finally, the derived emissivities are compared with other estimates from literature in Sect. 4.

## 2 Submm/mm Observations

Barnard 68 was observed with the Submillimetre Common User Bolometer Array (SCUBA) at the JCMT (Holland et al. 1999) and with the SEST Imaging Bolometer Array (SIMBA; Nyman et al. 2001). Both instruments have hexagonal arrays of 37 bolometers: SCUBA can observe in the submm at 850 $`\mu `$m with FWHM=14.5″ while SIMBA operates at 1.2 mm and has a beamsize of 24″.

SCUBA observed an area of 5′$`\times `$5′ around Barnard 68 in scan-map mode, chopping within the field of view. The resultant image is convolved with the chop function, which is removed by means of Fourier Transform analysis. The Emerson II technique is used to minimise the noise (Jenness et al. 1998): six scans of the field are made with different chop throws (20″, 30″ and 65″ in RA and Dec). In total, we coadded eight sets of scans, six observed in March/June 2002 and two in July 1998 (the latter retrieved from the SCUBA Archive). Standard data reduction was performed with the dedicated package SURF. After flat-fielding and masking of noisy bolometers, a baseline was removed from each bolometer. Maps were corrected for atmospheric extinction (the zenith optical depth was $`\tau _{850}`$=0.12,0.15,0.3 in the 1998 run and in the March/June 2002 runs, respectively) and for correlated sky noise. Calibration was obtained from scan-maps of Uranus. The final map has a sky noise of 20 mJy beam-1 (1-$`\sigma `$), equivalent to 3.5 MJy sr-1 for the SCUBA beamsize.

SCUBA scan-maps suffer from large scale undulations in the background of the source, which are due both to the poor sampling of small spatial frequencies and to uncertainties in the baseline removal (see, e.g. Visser et al. 2002). On our image, the undulations show up as a negative background around the source, of the order of the sky noise. We subtracted a mean background estimated on several sky apertures. The results of the next section are not significantly changed by the subtraction. The integrated flux at 850 $`\mu `$m inside the 2-$`\sigma `$ isophote is $`F_{850\mu \mathrm{m}}`$=4.0$`\pm `$1.0 Jy. SCUBA also took maps with the 450 $`\mu `$m array. The signal at this wavelength is poor, because of the background fluctuations and of the larger sky opacity ($`\tau _{450}0.5`$ for the best data). The peak emission is detected at 3-$`\sigma `$ level (1-$`\sigma `$ = 200 mJy beam-1; beamsize FWHM=10″) and it is not possible to conduct the surface brightness analysis described in Sect. 3. Nevertheless, we could estimate the integrated flux, which is $`F_{450\mu \mathrm{m}}`$= 13$`\pm `$4 Jy.

SIMBA observations were taken in fast scanning mode, without using a wobbling secondary mirror (Nyman et al. 2001). The bolometer array scans the observed area with a speed of 80″ s-1 in azimuth and a scan-to-scan separation of 8″ in elevation. Spatial frequencies are converted into temporal frequencies and the sky signal is filtered out with a low-cut filter. Each map of Barnard 68 covers $`600\mathrm{}\times 312\mathrm{}`$ in azimuth and elevation, and takes 7 minutes to complete. Barnard 68 was observed during June and October 2001. A total of 118 maps were coadded together, for a total integration time of 13.8 hours. The mean zenith optical depth is $`\tau _{1.2}`$=0.16.

Data were reduced with the software package MOPSI by R. Zylka. After deconvolution according to the frequency passband, subtracting a baseline to each scan, correcting for sky opacity and gain elevation the correlated sky noise fluctuations and artifacts due to electronics were removed. All maps of Barnard 68 were rebinned in a single image, with pixel size 8″. In an iterative way, a preliminary image was used as source model to improve the removal of correlated noise. Uranus maps were used to calibrate the data. Night-to-night fluctuations suggest a calibration uncertainty of 20% (1-$`\sigma `$). The final SIMBA image of Barnard 68 has a residual noise of 5.5 mJy beam-1 (1-$`\sigma `$), equivalent to 0.36 MJy sr-1 for the SIMBA resolution. The integrated flux is $`F_{1.2\mathrm{mm}}`$=0.7$`\pm `$0.2 Jy

The emission is compared to the extinction map of Alves et al. (2001) (the measured NIR colour excess has been converted to $`A_\mathrm{V}`$ using a standard reddening law; Rieke & Lebofsky 1985). All images have been smoothed to the SIMBA resolution, registered together (pointing accuracy is better than a few arcsec for the submm/mm observations) and resampled to a common pixel size of 12″. Fig. 1 shows the final SCUBA (left) and SIMBA (right) images, with superimposed $`A_\mathrm{V}`$ contours.

## 3 Analysis

The morphology of the extinction and emission maps in Fig. 1 is quite similar. Dust emission with $`S/N`$$`>`$3 traces regions with extinction $`A_\mathrm{V}`$$`>`$6 and $`A_\mathrm{V}`$$`>`$10, in the SCUBA and SIMBA images, respectively. In analogy with extinction, emission comes from a round region with a south east tail. At 1.2 mm, there are hints for a secondary emission peak. In Fig. 2 we show the pixel-to-pixel correlation of the emission maps with $`A_\mathrm{V}`$. As already seen in Fig. 1, there is a clear correlation between emission and extinction. For regions with $`A_\mathrm{V}`$$`>`$10 the scatter in the correlation with the submm data is of the order of the sky noise (1.7 MJy sr-1 for the smoothed image). For $`A_\mathrm{V}`$$`<`$10 the scatter increases because of the (possibly hotter) south-east tail. The pixels clearly belonging to the tail can be seen in Fig. 2 around $`A_\mathrm{V}`$=8 and $`I_\nu `$=10 MJy sr-1. The scatter at 1.2 mm is more uniform and slightly larger than the sky noise, because of the complex central morphology. A forthcoming paper will be devoted to the analysis of the central peaks and of the south-east tail.

For optically thin radiation (as for the submm/mm emission of Barnard 68) we can write

$$I_\nu =\frac{\kappa _\nu }{\kappa _\mathrm{V}}\times \frac{A_\mathrm{V}}{1.086}\times \frac{\rho B_\nu (T)𝑑l}{\rho 𝑑l},$$
(1)
where $`T`$ and $`\rho `$ are the temperature and density of dust grains, functions of the position within the cloud, $`B_\nu `$ is the Planck function and the integral extends along the line of sight through the dust. In the isothermal case, the last term on the right reduces to $`B_\nu (T)`$. Otherwise, a knowledge of the dust temperature and density distribution is needed.

In Fig. 3 we show the Spectral Energy Distribution (SED) of Barnard 68. Our integrated fluxes at 850 $`\mu `$m and 1.2 mm are similar to those presented by Ward-Thompson et al. (2002), from which we took the ISOPHOT fluxes at $`170\mu m`$ and $`200\mu m`$. We fitted a isothermal modified blackbody to the datapoints assuming a value for $`\beta `$, the emissivity spectral index ($`\kappa _\nu \nu ^\beta `$). For $`\beta `$=1.5-2 (Dunne & Eales 2001), we find acceptable fits for temperatures $`T_{\mathrm{iso}}=`$11-13 K. We adopt $`T_{\mathrm{iso}}=12\pm 2`$ K (1-$`\sigma `$) for the isothermal case, incorporating in the broad error the large uncertainties on $`\beta `$. The solid line in Fig. 3 is an indicative SED for $`\beta `$=1.7 and $`T_{\mathrm{iso}}`$=12 K. Similar temperature have been previously obtained for Barnard 68 (Ward-Thompson et al. 2002; Hotzel et al. 2002). Using $`T_{\mathrm{iso}}`$, we fitted Eq. 1 to the correlations of Fig. 2 and derived a single cloud-averaged value for the submm and mm emissivities. We obtained $`\kappa _{850\mu \mathrm{m}}/\kappa _\mathrm{V}=3.5\pm 1.010^5`$ and $`\kappa _{1.2\mathrm{mm}}/\kappa _\mathrm{V}=9.0\pm 3.010^6`$. Errors were estimated with a bootstrap technique and are dominated by the uncertainty on calibration and on temperature.

The fitted correlations are shown as solid lines in Fig. 2. While the fit is acceptable at 1.2 mm, a simple linear relation is unable to reproduce well the correlation observed at 850 $`\mu `$m for $`A_\mathrm{V}>20`$ since the $`I_\nu /A_\mathrm{V}`$ ratio decreases for increasing $`A_\mathrm{V}`$. This is expected if the temperature in the cloud core is lower than in the external part, as a result of the dust shielding of the external radiation field (Zucconi et al. 2001). Such a trend has been observed in other clouds (see, e.g. Kramer et al. 1998). A similar behaviour could also be explained with a isothermal dust distribution, if the dust emissivity is lower in denser regions. However, models and observations suggest that emissivities increase in dense cores (Ossenkopf & Henning 1994; Kramer et al. 2002).

We derived the temperature gradient inside the cloud from an improved version of the Zucconi et al. model . As in Alves et al. (2001), the adopted density distribution is that for a Bonnor-Ebert sphere, a pressure confined isothermal sphere in hydrostatic equilibrium. Dust is heated by a local Interstellar Radiation Field (Galli et al. 2002). The internal dust absorption is fixed by the measured extinction and the assumed absorption law $`\kappa _\nu /\kappa _\mathrm{V}`$, for which we used the tabulated values given by Ossenkopf & Henning (1994) for different models of dust coating/coagulation in dense clouds. While the absolute value of the temperature depends on the adopted dust model (for $`\lambda >30\mu `$m, $`\kappa _\nu /\kappa _\mathrm{V}`$ increases going from bare grains to grains with ice coating up to grains that undergo coagulation), the temperature gradient is found to be relatively independent of grain characteristics, with the temperature at the truncation radius $`T_{\mathrm{ext}}`$ ($`0.06`$ pc; Alves et al. 2001) about 1.5 times higher then the core temperature. This is because the absorption law does not change significantly with the dust model at shorter wavelengths, where most of the absorption occurs. In the following, we adopt the model radial gradient and derive $`T_{\mathrm{ext}}`$ from the data.

For the assumed density and temperature distributions, we simulated emission maps with 24″ resolution and 12″ pixels, in analogy with the data. The integrated SED is shown with a dashed line in Fig. 3. Again, a broad range of best fit temperatures is obtained, $`T_{\mathrm{ext}}`$=12-15 K. We use $`T_{\mathrm{ext}}=14\pm 2`$ K (1-$`\sigma `$). Finally, the simulated $`A_\mathrm{V}`$ vs $`I_\nu `$ correlation of Eqn. 1 is fitted to the observed data to derive $`\kappa _\nu /\kappa _\mathrm{V}`$. We obtained222Equivalent to $`\kappa _{850\mu \mathrm{m}}=1.5\pm 0.4\mathrm{cm}^2\mathrm{g}^1`$ and $`\kappa _{1.2\mathrm{mm}}=0.35\pm 0.1\mathrm{cm}^2\mathrm{g}^1`$ (assuming grain radius 0.1$`\mu `$m, grain density 3 g cm-3 and V-band extinction efficiency 1.5; Hildebrand 1983).

$`\kappa _{850\mu \mathrm{m}}/\kappa _\mathrm{V}`$ $`=`$ $`4.0\pm 1.010^5`$
$`\kappa _{1.2\mathrm{mm}}/\kappa _\mathrm{V}`$ $`=`$ $`9.0\pm 3.010^6.`$
These correlations are also shown in Fig. 2 with dashed lines. Although the fit to the submm data has improved, the emissivities are nearly the same as for the isothermal case.

## 4 Discussion

The derived submm/mm emissivities are shown in Fig. 4 together with estimates from literature. Emissivities derived from IRAS/COBE observations of diffuse Galactic dust (Boulanger et al. 1996; Bianchi et al. 1999) are within 10% of the values predicted by the popular Draine & Lee (1984) model, which is shown in Fig. 4 as a solid line. Dotted lines refer to the Ossenkopf & Henning (1994) models for bare grains, grains with a thin ice coating and grains with thin ice after coagulation has proceeded for 105 years in a gas with density 106 g cm-3 (from bottom to top, respectively).

James et al. (2002) derived the emissivity from a sample of local galaxies observed with SCUBA, estimating the dust mass from metals. Their value is compatible to that of Draine & Lee (1984), indicating similar properties for diffuse dust in external galaxies and in the Milky Way. Our 850 $`\mu `$m emissivity is higher than that of diffuse dust: a comparison with the models of Ossenkopf & Henning (1994) seems to suggest that grains in Barnard 68 possess molecular ice mantles or have coagulated into fluffy aggregates. However, $`\kappa _{850\mu \mathrm{m}}/\kappa _\mathrm{V}`$ is within 2$`\sigma `$ of that for bare grains. In a very recent paper, Kramer et al. (2002) measured $`\kappa _{850\mu \mathrm{m}}/\kappa _\mathrm{V}`$ for a molecular ridge and four embedded cores in the dark cloud IC5146. Their radially averaged values of the four cores, derived assuming $`\beta `$=2, range from $`2.110^5`$ to $`6.110^5`$ and are broadly compatible with our determination. At 1.2 mm, our emissivity is similar to that of diffuse dust but lower than the value previously estimated by Kramer et al. (1998) on a IC5146 core.

The ratio $`\kappa _{850\mu \mathrm{m}}/\kappa _{1.2\mathrm{mm}}`$ we derive implies a value for $`\beta `$ larger than what expected from measurements on Galactic diffuse dust and on laboratory cosmic dust analogues (Mennella et al. 1998). However, the error on the ratio is large and $`\beta =2`$ is still compatible (within 1-$`\sigma `$) with our result. Contamination of the submm and mm fluxes by the $`{}_{}{}^{12}\mathrm{CO}_{}^{}`$(3-2) and $`{}_{}{}^{12}\mathrm{CO}_{}^{}`$(2-1) lines (Avery et al. 1987) was found to be negligible.

Clearly, more observations are needed to narrow down the errors in the determination of the emissivity. The analysis presented here has to be repeated on a large sample of objects of relatively simple morphology like Barnard 68. Only with a statistical sample will it be possible to study the variation of dust properties with the environment.

###### Acknowledgements.

We are grateful to J. Alves, R. Cesaroni, A. Natta, C. Lada and C. Kramer (the referee) for useful comments and discussions.


text:
# Hard GRB spectra: thermal vs non–thermal emission

## 1. Introduction

Since their discovery, the time resolved spectral analysis of GRBs has been a testing ground–floor for the models proposed for the $`\gamma `$–ray prompt emission (Ford et al. 1995; Preece et al. 2000). Nonetheless, the mechanism(s) responsible for their emission is still an open problem (Ghisellini 2003). Different authors (Crider et al. 1997, 1999; Preece et al. 1998, 2002; Ghirlanda et al. 2002) have shown that the low energy photon spectrum (typically represented by a powerlaw $`N(E)E^\alpha `$) can be harder than the limit of the optically thin synchrotron model ($`N(E)E^{2/3}`$, Katz 1994; Tavani 1996) which is the most popular mechanism proposed for interpreting the burst emission. Modifications of the simplest optically thin synchrotron model (Lloyd & Petrosian 2000, 2002; Medvedev 2001) or variants based on Comptonization (Liang et al. 1997; Ghisellini & Celotti 1999; Lazzati et al. 2000) have been proposed to solve the inconsistency between theory and observations. All these models predict different limits for the low energy spectral hardness and can be directly tested by the comparison with the hardest bursts observed by BATSE.

We present the results of the time resolved spectral analysis of 5 bursts with a low energy spectral component harder than a flat photon spectrum, i.e. $`\alpha 0`$, for most of their duration. These results are shown to constrain the above emission models (Ghirlanda, Celotti & Ghisellini 2003).

Another interesting aspect of GRB spectra which we consider is their possible thermal character, which is predicted by the fireball model (Goodman 1986; Paczynski 1986). We find that the spectrum in the initial phase of these extremely hard bursts is successfully fitted with a thermal black body emission. The interpretations and implications of these results are discussed (Ghirlanda et al 2003) in the framework of a photospheric (e.g. Mészáros & Rees 2000, Daigne & Mockovitch 2002) and of the Compton drag model (Lazzati et al. 2000; Ghisellini et al. 2000).

## 2. Results

GRB 980306 represents a newly discovered case of hard burst. The spectral evolution of its first 6 s is reported in Fig.1 (left panel). The evolution of $`\alpha `$ (mid panel of Fig.1) indicates that, for most of the pulse, the low energy photon spectrum is harder than $`E^0`$. The maximum hardness ($`\alpha =1.1\pm 0.2`$) is reached at $`t1.5`$ s after the trigger.

In GRB 911118 (Fig.1, right panel) the low energy spectral index is $`\alpha 0`$ for the first $`6`$ s (32 spectra, phase A) The hardest photon spectrum (at $`t1.5`$ s after the trigger) has $`\alpha =0.74\pm 0.13`$.

In order to have homogeneous results we also re-analyzed three bursts (GRB 910807, GRB 910927, GRB 970111) already presented in the literature for their exceptional low energy spectral hardness (Crider et al. 1997, 2000; Frontera et al. 2000; Preece 2002).

Considering the relevance of these 5 GRBs for the comparison with the models, we also tested the independence of their low energy spectral hardness from possible sources of systematic error such as the detector response, the background calculation or the high energy spectral component. Moreover, we directly verified the low energy spectral hardness comparing the instrumental spectra with a simulated flat spectrum. These tests (Ghirlanda et al. 2003) indicate that most of the time resolved spectra of these bursts are harder than $`E^0`$ at more than 3$`\sigma `$.

As anticipated, we performed fits on the hardest spectra also with a black body model. In most of these bursts the initial 1.5 s of their emission is consistent with a black body with a temperature that decreases with time as roughly $`T_{BB}t^{1/4}`$ starting form $`kT100`$ keV (Ghirlanda et al. 2003).

## 3. Discussion

GRB 980306 and 911118, together with GRB 910807, 910927 and 970111 represent a challange for the proposed emission processes. Their low energy photon spectrum is harder than $`N(E)E^0`$ for a major part of the pulse/s reaching a maximum hardness of $`E^{1.5}`$. These results are inconsistent with the optically thin synchrotron model which can not account for spectra harder than $`E^{2/3}`$ (Katz 1994). Also its variants like synchrotron self absorption (Papathanassiou 1999; Granot, Piran & Sari 2000), small pitch angles of the emitting particles (Lloyd & Petrosian 2000) or Jitter radiation (Medvedev 2001) can produce a spectrum at most as hard as $`N(E)E^0`$, which is not compatible with the totality of these findings. Comptonization models (Liang et al. 1997; Ghisellini & Celotti 1999) also have difficulties, although they are consistent with very hard spectra in a limited range of energies.

The possible thermal character of the initial phase of all these 5 bursts can be interpreted as the emission of the fireball shells when they reach the transparency radius (Paczynski 1986; see also Daigne & Mockovitch 2002). For typical fireball parameters and using our observational findings this radius results $`R(\tau 1)5\times 10^{13}`$ cm and the corresponding bulk Lorentz factor is $`\mathrm{\Gamma }1000`$. An alternative scenario is the production of an observed thermal spectrum by the Compton drag of soft (seed) photons due to the fireball bulk motion (Lazzati et al. 2000). In this case the seed photons are required to have a temperature of $`5\times 10^4`$ K in order to produce the observed thermal spectrum. These results are consistent with the typical values assumed in the thoery of GRBs (e.g., Rees & Mészáros 1994).

At present, we are not able to discriminate between the photospheric and the Compton drag scenario. Consider also that emission at times greater than a few seconds can well be due to other processes, possibly linked to internal shocks starting to dominate at later phases. However, a key difference between the two scenarios is the fact that the seeds photons for the Compton drag process can be “used” only by the first shells, because the time needed to refill the scattering zone with new seeds exceeds the duration of the burst. Therefore observing black body emission for a long time, or during the rising phase of two time resolved peaks of the same GRB would be difficult to explain in terms of the Compton drag process.

### Acknowledgments.

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center On line Service, provided by the NASA/Goddard Space Flight Center. Giancarlo Ghirlanda and AC acknowledge the Italian MIUR for financial support.

## References

Crider A., Liang E.P., Smith I.A., 1997, ApJ, 479, L39

Crider A. et al., 1999, ApJ, 519, 206

Crider A. & Liang E.P., 2000, ApJ, 127, 283

Daigne F. & Mochkovitch R., 2002, MNRAS, 336, 1271

Ford L.A., Band D.L. & Matteson J.L., 1995, ApJ, 439, 307

Frontera F., Amati L., Costa E., et al., 2000, ApJS, 127, 59

Ghirlanda G., Celotti A. & Ghisellini G., 2002, A&A, 393, 409

Ghirlanda G., Celotti A. & Ghisellini G., 2003, subm. to A&A (astro-ph/0210693)

Ghisellini G. & Celotti A., 1999, ApJ, 511, L93

Ghisellini G., Lazzati D., Celotti A. & Rees M.J., 2000, MNRAS, 316, L45

Ghisellini G. 2003, astro-ph/0301256 see also these Proceedings

Goodman, 1986, ApJ, 308, L47

Granot J., Piran T. & Sari R., 2000, ApJ, 534, L163

Katz J.I., 1994, ApJ 432, L107

Lazzati D., Ghisellini G., Celotti A. & Rees M.J., 2000, ApJ, 529, L17

Lloyd N.M. & Petrosian V., 2000, ApJ, 543, 722

Lloyd N.M. & Petrosian V., 2002, ApJ, 565, 182

Liang E., Kusunose M., Smith I.A. & Crider A., 1997, ApJ, 479, L35

Medvedev M.V., 2001, astro-ph/0001314

Mészáros P. & Rees M.J., 2000, ApJ, 530, 292

Paczynski B., 1986, ApJ, 308, L43

Papathanassiou H., 1999, A&AS, 138, 525

Preece R.D. et al 1998, ApJ, 506, L23

Preece R.D. et al., 2000, ApJS, 126, 19

Preece R.D., Briggs M.S., Mallozzi R.S., Pendleton G.N., Paciesas W.S. & Band D.L., 2002, ApJ, 518, 1248

Rees M.J. & Mészáros P., 1994, ApJ, 430, L93

Tavani M., 1996, ApJ, 466, 768


text:
# Magnetic field evolution in galaxies interacting with the intracluster medium

## 1 Introduction

It is well known that galaxies in clusters have different physical and morphological properties than field galaxies. They are often Hi deficient, i.e. they have lost up to 90% of their atomic gas - mainly due to the interaction of the interstellar medium (ISM) with the intra-cluster medium (ICM) (Solanes et al. 2001). The cluster galaxies often possess Hi tails and strong compression regions (Cayatte et al. 1994, Bravo-Alfaro et al. 2000, Chyży et al. in prep.). These observed asymmetries of the gas distribution can be explained either by ram, pressure effects (Vollmer et al. 2001 and references therein) and/or by tidal interactions with nearby companions (e.g. NGC 4438, Combes et al. 1988). Recent radio continuum observations of Virgo Cluster galaxies show also magnetic field structures which are not present in normal galaxies (Soida et al. 1996, Vollmer et al. in prep.). NGC 4254, oriented face-on, shows strong radio polarization maxima distributed outside the optical spiral arms in the direction towards the Virgo Cluster center (Soida et al. 1996, Chyży et al., 2001). The highly inclined galaxy NGC 4522 also shows asymmetric radio polarization destinations (Vollmer et al. in prep.).

While the question of the gas reaction to the ICM ram pressure has a rather long history (Takeda, Nulsen, & Fabian 1984, Gaetz, Salpeter & Shaviv 1987; Balsara, Livio, & O’Dea 1994; SPH: Tosa 1994; Abadi, Moore, & Bower 1999, Schulz & Struck 2001; Vollmer et al. 2001), there were no investigations yet concerning the evolution of magnetic fields in cluster galaxies. Therefore, a comprehensive study of magnetic field behaviour under the influence of the gas stripping process is needed to provide a physical interpretation of the observed features in the radio continuum and polarization in stripped spiral galaxies.

The evolution of the magnetic field structure in galaxies is closely related to the amplification of the large scale magnetic field by a dynamo action. The main clue of the traditional MHD dynamo is the rate of small scale turbulent motions in the rotating medium producing large-scale poloidal magnetic field from the azimuthal one. Then the differential rotation reproduces the toroidal magnetic field with its strength amplification (e.g. Parker 1979 and references therein).

We solve the induction equation for the gas velocity fields obtained from the N–body simulations by Vollmer et al. (2001). With our simulations we analyze the problem of a total magnetic energy amplification in galaxies interacting with the ICM. In a first step no dynamo action is included.

The main goal of our project is to study a representative case of ICM–ISM interaction to investigate the characteristic features of the magnetic field evolution. This will be a first step to understand the observed magnetic fields in cluster galaxies.

The outline of this article is as follows: the model used to simulate the evolution of the magnetic field is introduced in Sect. 2. The resulting evolution of the magnetic field topology is shown in Sect. 3.1. In Sect. 3.2 the evolution of the total magnetic energy is investigated, followed by the construction of model maps of polarized radio emission (Sect 3.4). We discuss our results in Sect. 4 and summarize them in Sect. 5.

## 2 Model

We constructed a 3D MHD model of magnetic field evolution (see Otmianowska-Mazur et al. 1997 for more details) in a galaxy moving rapidly through the ICM. We apply the Zeus3D code (Stone & Norman 1992a and b) solving the induction equation:

$$B/t=\text{rot}(v\times B)\text{rot}(\eta \text{rot}B)$$
(1)
where $`B`$ is the magnetic induction, $`v`$ is the large-scale velocity of the gas, and $`\eta `$ is the coefficient of a turbulent diffusion. Realistic, time dependent gas velocity fields are provided by the 3D N-body sticky-particle simulations of Hi cloud complexes evolving in an analytically given gravitational potential of a disk, a bulge and a halo (Vollmer et al. 2001). The clouds collide inelastically and are affected by the assumed ICM ram pressure due to relative motion of the galaxy with respect to the ICM.

### 2.1 Model input parameters

We model a galaxy which is mildly Hi deficient (it has lost $``$30% of its gas) at the end of the simulation. It is on an eccentric orbit within the cluster and its disk plane is inclined by 20 with respect to the orbital plane. When the galaxy is approaching the cluster center the clouds are pushed and possibly stripped away in the direction opposite to its velocity respect to ICM. When the galaxy emerges from the cluster core, a process of cloud re-accretion starts due to the galactic gravitational potential. The sticky-particle code uses 10000 cloud complexes with different masses. Ram pressure is modeled as an additional acceleration on the clouds located at the surface of the gas distribution in the direction of the galaxy’s motion within the cluster. Thus, there is no explicit intra-cluster medium included in the model. The maximum ram pressure applied on the clouds corresponds to an intra-cluster medium density of $`2\times 10^3\mathrm{cm}^3`$ (see Vollmer et al. 2001 for more details) and a maximum velocity of the modeled galaxy of about 1500 km s-1 with respect to the cluster center. The time step of our calculations is $`10^7`$ yr. The time of closest approach to the cluster center is $`t`$=500 Myr after the start of the simulations.

The induction equation is analyzed in rectangular coordinates ($`XYZ`$). The number of grid points used is 171x171x71 along the $`X`$, $`Y`$ and $`Z`$ axis, respectively. This corresponds to the grid spacing of 200 pc in the galactic plane and of 280 pc in the $`Z`$ direction, resulting in a size of the modeled box of 34 kpc x 34 kpc x 20 kpc. Since the N-body code is discrete whereas the MHD code is using a grid, we have to interpolate the discrete velocities on the grid. This is done using a spline function with a density dependent smoothing length. It turned out that we had to use a large smoothing length to suppress the noise in the velocity field of the outer disk, which is due to a small, local particle density. In this way we avoid numerical artifacts of the magnetic field distribution at the outer disk. This has the consequence that the velocity field at the edge of the gas distribution is more extended than the gas distribution itself, which leads to an overstanding magnetic field at the edges of the gas distribution.

The large number of grid points gives rise to a computational difficulty: in our earlier papers (eg. Otmianowska-Mazur et al. 2002, Elstner et al. 2000) the cell size in $`Z`$ direction was 80 pc, whereas in the present simulations it is 280 pc, thus the value of the numerical diffusion in the $`Z`$ direction is certainly higher. In order to avoid such diffusion in the central part of our box during the first 20 time steps the velocity curve was smoothed from the distance of 5 kpc to the center with a Gaussian function with $`\sigma =2`$ kpc. This smoothing influences the magnetic field configuration in the later time steps only marginally, but allows to keep magnetic flux in the central disk constant during the first 200 Myr of its evolution. We assume the magnetic field to be partially coupled to the gas via the turbulent diffusion process (Otmianowska-Mazur et al. 2002). We made two basic experiments with two different behaviours of the turbulent diffusion coefficient: uniform (case A) and linearly growing with $`z`$ in the halo (case B) (see Otmianowska-Mazur et al. 2002 for a detailed explanation). Its value is $`5\times 10^{25}\mathrm{cm}^2\mathrm{s}^1`$ in the whole box for case A. For experiment B it has the same value in the disk plane. To check how numerical diffusion influences our results, two additional computations were made: no physical diffusion (only numerical diffusion) (case C), and uniformly distributed diffusion 5 times smaller than in case A (case D). The initial magnetic field is purely toroidal with a strength of 1 $`\mu `$G.

## 3 Results

### 3.1 Magnetic field geometry

All calculated models with different diffusion coefficient result in qualitatively similar magnetic field structures. Therefore, we discuss only case A (uniform diffusion) as the simplest one. The effects of the ICM ram pressure on the ISM of a galaxy moving rapidly through the ambient ICM can be divided into two main phases: the period of cloud pushing/stripping called also the compression phase and of cloud re-accretion (see Vollmer et al. 2001). Both of these processes influence significantly the magnetic field morphology. The stripping phase takes place roughly between 300 and 600 Myr. Re-accretion happens mainly between 700 and 1000 Myr.

In Figs.1-6 we present snapshots of 6 different timesteps: $`t`$=450, 550, 650, 750, 850, and 950 Myr. In each figure the magnetic field vectors in the midplane of the disk are shown while the total gas density is presented in a logarithmic scale. The galaxy is viewed face-on ($`X`$-$`Y`$ plane; left panel) or edge-on (right panel). In order to emphasize the most interesting features we show the $`X`$-$`Z`$ plane in Fig. 35 and the $`Y`$-$`Z`$ plane in Fig. 1, 2 and 6. The maximum lengths of the magnetic field vectors are $`B=5\mu `$G in the $`X`$-$`Y`$ plane and $`B=2.5\mu `$G in the $`X`$-$`Z`$ and $`Y`$-$`Z`$ planes respectively. The galaxy moves towards south-east, which is equivalent to the wind flowing from this direction. In the $`X`$-$`Z`$ and $`Y`$-$`Z`$ plane the wind also comes from the south-east.

#### 3.1.1 The compression phase

Shortly before the ram pressure reaches its maximum ($`t`$=450 Myr), the gas is strongly compressed forming a heavy spiral arm on the windward side of the disk (Fig. 1). The magnetic field is swept together with the gas and also forms a strong maximum in the direction of the wind (south-east). Since the compressed magnetic field rotates counter-clockwise together with the gas (frozen-in field), its maximum is located in the south-west of the galaxy. We observe a shift between the maximum of the magnetic field and that of the gas. Because there is no ISM located further out to the south–east than the gas maximum, the offset magnetic field can not be associated with the ISM. Three physical mechanisms can in principle create such an offset: (i) the pile–up of magnetic field associated with the ICM, (ii) turbulent diffusion, and (iii) shear amplification between the rotating galaxy and the intra-cluster medium. Since there is no explicit ICM and thus no magnetic field associated with the ICM included in the model and diffusion is much too slow to play a role, mechanisms (i) and (ii) are excluded. In order to estimate the growth of the magnetic field due to shear, we compare it to the growth due to gas compression:

$$\frac{B_\varphi }{t}=B_\varphi \frac{v_\mathrm{r}}{r}.$$
(2)
The growth of the magnetic field due to shear is:

$$\frac{B_\varphi }{t}=B_\mathrm{r}r\frac{\mathrm{\Omega }}{r},$$
(3)
where $`B_\varphi `$, $`B_\mathrm{r}`$ are the components of the magnetic field in azimuthal, radial direction, $`v_\mathrm{r}`$ is the velocity component in radial direction, $`\mathrm{\Omega }`$ is the angular velocity, and $`r`$ is the galactocentric distance. The derivative $`\frac{v_\mathrm{r}}{r}`$ can be estimated by $`\frac{\mathrm{\Delta }v_\mathrm{r}}{\mathrm{\Delta }r}`$ 40 km s-1/5 kpc. Since the wind accelerates the outer gas particles $`\frac{\mathrm{\Omega }}{r}`$ is reduced, thus’ $`r\frac{\mathrm{\Omega }}{r}<\mathrm{\Omega }<`$ 150 km s-1/5 kpc. In the south-east the magnetic field is mainly azimuthal $`B_\varphi B_\mathrm{r}`$. Thus, the amplification of the magnetic field due to gas compression exceeds that due to shear. We conclude that the shift between the gas and magnetic field maxima in the south-east is a pure numerical artifact. It is due to the use of a density dependent smoothing length together with a fixed minimum length. The smoothing length at the edge of the galaxy’s gas distribution is small, whereas slightly outside, it becomes much larger. The offset thus reflects the gradient of the smoothing length across the edge of the gas distribution.

In the north-western part of the disk the visible distribution of magnetic field structures is caused by stripping effects. In the western part of the outer disk the magnetic field forms small waves and changes its direction locally at an almost right angle. In the $`Z`$ direction (Fig. 1 right panel) an asymmetry in the distribution of magnetic field vectors can be observed in the compression region. Since the wind direction has a radial component and a component in the $`Z`$ direction, the magnetic field begins to be pushed radially inwards and towards positive $`Z`$.

At the timestep of 550 Myr (Fig. 2) the gas has formed a prominent accelerated arm in the west ($`X`$=10 kpc, $`Y`$=3 kpc). The magnetic field enhanced due to compression effects follows the rotation of the gas and moves counter-clockwise. In addition, being pushed by the wind it moves inwards. Its maximum strength is now located west ($`X`$=10 kpc, $`Y`$=0 kpc) of the galaxy closer to its center than for $`t`$=450 Myr. The magnetic maximum is offset from the gaseous one, again due to the gradient of the smoothing length across the edge of the gas distribution. The magnetic spiral arms in the north as well as the waving structure in the west are more pronounced than at $`t`$=450 Myr. In the $`Z`$ direction the compressed magnetic field of the south-eastern part of the galaxy is swept together with gas by few kpc above the disk plane (much higher than the initial magnetic field; Fig. 1).

At the time $`t`$=650 Myr (Fig. 3) the galaxy slows down on its orbit away from the cluster core. The accelerated arm has moved to the north of the galaxy and has a complex three dimensional structure. The gas maxima start to be redistributed and wind-compression features diminish their strength. However, at the south-east disk edge magnetic maxima are still present (Fig. 3 left). The northern part of the outer disk is filled with strong magnetic arms stretching away from the galaxy. At the south–west side of the disk ($`X`$=10 kpc, $`Y`$=-3 kpc) a counter-rotating gaseous arm forms. The visible maximum of the magnetic strength there is connected with the falling gas. In the $`X`$-$`Z`$ plane the magnetic field extends to even more positive $`Z`$ than for $`t`$=550 Myr, whereas the magnetic field at the location, where the counter-rotating arm forms, extends to negative $`Z`$.

#### 3.1.2 The phase of re-accretion

The period of re-accretion begins at $`t700`$ Myr. The galaxy moves slowly through outer regions of the cluster. The wind has ceased and a large number of gas clouds starts to fall back onto the galactic disk producing very violent gas fountains. The magnetic field is taken up and down together with the gas and forms poloidal loops extending high in the halo (even to 7 kpc, see Fig. 5 right panel).

At $`t`$=750 Myr (Fig. 4) the accelerated arm, which falls back onto the galactic disk in the north-east, produces a distribution of maximum gas density similar to that produced by the wind (Fig. 1). The counter-rotating arm, which has a very small mass, has no visible influence on the gas distribution, but generates magnetic field at negative $`Z`$ that are not visible in Fig. 4, because it is not located in the $`X`$$`Z`$ plane. Due to the strong shear motions at the trailing edge of the back-falling, accelerated arm, the magnetic field is enhanced there. It extends towards positive $`Z`$. A second maximum of the magnetic field is visible in the south-west where an expanding gas shell starts to form.

At $`t`$=850 Myr (Fig. 5) the ridge of the enhanced gas density has moved to the north-east, in a counter-rotating sense. The maximum of the magnetic field follows this ridge. Additionally, a strong magnetic arm is formed to the north, whose radial distance from the galactic center increases to the west. For the first time a weak magnetic arm is visible east of the galaxy. In the south-west ($`X`$=5 kpc, $`Y`$=-10 kpc) the expanding gas shell is now clearly visible. The magnetic field structure in the $`Z`$ direction became very complex. Several magnetic loops with heights up to several kpc above the disk plane can be observed in Fig. 5 (right panel).

At the last presented timestep ($`t`$=950 Myr, Fig. 6) we observe the beginning of a second expanding gas shell south of the galaxy center, which coincides with a strong ridge of magnetic field. Thus, the beginning of the expansion of a gas shell is accompanied by a strong magnetic arm. In the western part we observe strong open spiral arms. The magnetic field structure in the $`Z`$ direction shows now several loops extending up to 10 kpc in both direction of the disk plane.

### 3.2  Evolution of the total magnetic energy

The evolution of the total magnetic energy normalized to its initial value is presented in Fig. 7 for our 4 experiments: with a uniformly distributed magnetic diffusion coefficient (case A), with a diffusion coefficient growing into the halo (case B), without a physical diffusion (case C), and with a diffusion coefficient 5 times smaller than its basic value (case D). We have calculated the magnetic energy (i) in the whole box (34 kp$`\times `$34 kpc$`\times `$20 kpc) and (ii) in a thick disk centered on the galaxy center with a disk height of 1.8 kpc.

For cases A and B the total magnetic energy $`E`$ increases from its initial value $`E_0`$ and reaches its maximum $`E=2E_0`$ at $`t`$430 Myr. In the cases of the reduced physical diffusion (case C) and no physical diffusion (case D) the total magnetic energy rises to 3 times its initial value. The ram pressure compresses the gas at the windward side and gives rise to a radial component of the velocity field there. This enhances the magnetic field via the $`(v\times B)`$ term in the induction equation (Eq. 1). At the same time the wind expels the magnetic field out of the galaxy and out of the box of calculation. ¿From $`t`$=500 Myr to 700 Myr the sweeping of the magnetic field becomes more important than the generation of magnetic field due to compression. The magnetic energy in the disk and in the whole galaxy even decreases slightly below its initial value.

The process of re-accretion, which starts at about 650 Myr, causes a fast growth of the total magnetic energy in the whole galaxy from $`E4.5E_0`$ (case B, long dashed line) to $`E7.5E_0`$ (experiment A, dashed line). Within the thick disk the magnetic energy increases to $`E3E_0`$ and $`E3.5E_0`$, respectively. Moreover, the increase of the magnetic energy within the disk is slower than in the whole galaxy. In the re-accretion phase the magnetic energy within the whole box reaches a maximum for case A/B at $`t1050`$ Myr. The magnetic energy within the disk reaches its maximum $``$150 Myr later. We observe a clear decline of the magnetic energy within the disk and the whole galaxy after the maximum for case A and B.

The experiments with a reduced physical diffusion (case D) and without physical diffusion (case C) lead to a higher increase of the total magnetic energy than our basic two experiments A and B. The energy grows 10 times its initial value during the re-accretion phase. During the last $``$700 Myr the total magnetic energy stays almost constant for case D and decreases only very slowly for case C.

The influence of the physical diffusion can be directly observed in comparing case A with cases C/D (Fig. 7). Whereas the graphs for case C and D are very close to each other, those of case A (for the whole galaxy) are situated well below. Comparing cases A, C, and D for 1250 Myr$`<t<`$1350 Myr one can derive a roughly linear correlation between total magnetic energy and the diffusion coefficient. Thus, the physical diffusion clearly dominates over the numerical diffusion.

The linearly growing diffusion coefficient also has a non negligible influence on the evolution of the magnetic energy. During the re-accretion phase its maximum of $`E/E_0`$ is decreased by a factor $``$1.4 with respect to $`E/E_0`$ for the uniform diffusion coefficient. For later timesteps this factor decreases to a value of 2.3 at the end of the simulations. These factors are slightly lower for the magnetic energy computed only within the disk.

### 3.3 The influence of the cloud–cloud collision rate

In order to investigate the influence of the cloud–cloud collision rate on the evolution of the total magnetic energy, we have increased the radii of the clouds by a factor 6. This leads to an increase of the number of collisions per unit time of a factor 36. The evolution of the total magnetic energy for the simulations with the high and the low cloud–cloud collision rate is shown in Fig. 8.

In the case of a high cloud–cloud collision rate (solid line) the magnetic field energy grows faster and to higher values until $`t`$=450 Myr. This is due to formation of spiral structures within the disk, which is enhanced by a higher cloud–cloud collision rate. For $`t>450`$ Myr the magnetic energy evolves in the same way as that of the simulation with the low cloud–cloud collision rate, but with an approximately constant offset. At $`t=1000`$ Myr both curves reach the same maximum. For $`t>1000`$ Myr the total energy drops faster in the case of the high cloud–cloud collision rate. Within 500 Myr it decreases by less than half of its maximum value. This decrease is caused by the higher damping of vertical and shear motions due to frequent cloud–cloud collisions.

### 3.4 Polarization maps

In order to analyze the general properties of the magnetic field we use simulated maps of the polarized intensity (PI), constructed on the basis of our resulting magnetic fields. This process yields the information which can be directly compared with polarized radio continuum observations. We performed it in the way similar to that described in our previous works (e.g. Otmianowska-Mazur et al. 2000, Otmianowska-Mazur et al. 2002). The Stokes parameters I, Q and U are integrated along the line of sight and convolved with an assumed beam of 10 ′′. We then computed maps of polarized intensity and the angles of polarization B-vectors.

At selected time steps (see Fig. 9) of our model we construct maps of the radio polarization intensity (PI) with the inclination of $`0^{}`$ (face-on) to the sky plane. We use the density distribution of relativistic electrons that has the form

$$n_{\mathrm{rel}}\mathrm{exp}((R/R_0)^2)\mathrm{exp}((z/z_0)^2),$$
(4)
where $`R`$ is the radial distance to the galaxy center and $`z`$ is the vertical distance from the midplane of the disk. We set $`R_0`$=10 kpc and $`z_0`$=1 kpc. The distance of the galaxy is assumed to be $``$34 Mpc. Thus, one arcminute roughly corresponds to one kpc. The beamsize is 10′′. The PI maps of model A with high and low cloud–cloud collision rates are very similar. We show the polarized intensity together with the polarization B vectors and the gas surface density in Fig. 9 for model A with the uniformly distributed diffusion and the low cloud–cloud collision rate.

We will focus here only on the brightest PI emission features that arise from the ICM–ISM interaction. They can directly be compared with radio continuum observations of cluster galaxies. At $`t`$=450 Myr the gas is compressed in the south-east leading to an excess of density in the gas and a bright ridge of polarized emission. This ridge is symmetric and has the form of a half-ring. The offset in our model snapshot is again a numerical artifact due to the gradient of the smoothing length at the edge of the galaxy’s gas distribution (cf. Sect. 3.1.1).

Shortly after the galaxy’s nearest approach to the cluster center ($`t`$=550 Myr) the accelerated arm is formed in the west. It is again accompanied by a ridge of the strong polarized emission, which has its maximum in the west of the galaxy center. Again, the PI emission maximum is shifted outward in comparison with the maximum of the gas surface density.

At $`t`$=650 Myr the re-accretion phase begins and the gas falls back onto the disk in the north-east of the galaxy. In the east, where the gas hits the galactic disk, a strong polarized emission region forms. We observe also ridges of PI at the interface between the infalling and galactic gas in the north. In addition, one can still observe a polarized intensity ridge in the west, that represents the maximum PI emission. Both features form an incomplete ring structure.

This ring structure becomes more complete and symmetric at $`t`$=750 Myr with the maximum of PI located in the east of the disk being again slightly shifted outward with respect to the maximum of the gas surface density. This shift is real, because it is caused by backfalling gas that hits the disk in this region.

At $`t`$=850 Myr the gaseous clouds re-accrete form the north and hit the disk in the north-east. The ring structure in PI that is observed at $`t`$=750 Myr has now expanded. It becomes weaker where the expanding shell forms. The PI maximum is located in the north-east. A gas shell starts to expand to the south-west partly accompanied by the polarized intensity.

At the timestep of $`t`$=950 Myr a large and extended region of PI forms in the north-west where the gas of the expanding shell flows back onto the galactic disk. Most of the flux in polarized emission is located in the north of the galaxy. The most southern PI maximum is located where the expanding shell is anchored.

This north–western region of the extended polarized emission moves in the direction of rotation and appears in the north of the galaxy at $`t`$=1050 Myr. One can still observe a residual ridge of PI in the south-east and a new narrow feature of PI running north-south.

For the evolution of the PI, one can again distinguish two different evolutionary phases: (i) the compression phase and (ii) the re-accretion phase. While the ridges of maximum PI are mainly caused by the compression for $`t<`$600 Myr, they are amplified due to shear via re-accretion for $`t>`$600 Myr. Whereas the maxima of polarized emission should coincide with the maxima of the gas surface density in the compression phase, they can be located outside in the re-accretion phase.

## 4 Discussion

Our numerical experiment clearly indicates that the magnetic field is amplified in the whole disk and even outside the galaxy during the whole period of the evolution. There are several physical processes responsible for this fact. In the first phase, when the galaxy moves through the ICM, the mechanism of a strong compression takes place at the windward side of the galaxy. The azimuthal component of the magnetic field is easily amplified in this process, but also a radial component is created due to wavy gas motions (see Fig. 2 left panel). At later times of the galaxy evolution this radial field is transformed by the differential rotation into an azimuthal field (similarly to the MHD dynamo). This represents the most effective process of magnetic field growth. The re-accretion phase is much more complicated than the compression phase. The presence of very violent inflows and outflows, shears connected with them and a differential rotation (resulting in the reproduction of azimuthal field) gives again possibility of an efficient increase of the total magnetic energy. Inflowing gas moves in a spiral-like manner towards the disk forming strong magnetic field outside the disk and gaseous maxima at the interface with the disk gas.

In order to explain a complicated picture of relations of the velocity and magnetic field relations we present different velocity components as vectors projected onto the magnetic field strength shown as a grey-plot. Fig. 10 shows the galactic plane at the time step of $`t`$=450 Myr (the wind compression phase), while Fig. 11 presents two perpendicular cuts of the box at $`t`$=650 Myr (the re-accretion period). We choose the same time steps as in Fig. 1 and Fig. 3. Our modeled velocity is extrapolated into the whole box in order to avoid too strong velocity gradients in our calculations. This is the only way to avoid an artificial shear at the edges of the galaxy’s gas distribution.

The radial velocity (Fig. 10) serves as the best indicator of compression regions, where the strong gradients of this velocity component are observed. One of them is visible at the south-east part of the galaxy and is placed in the direction of the wind causing the formation of the magnetic field maximum in the central part of this quarter (near the gaseous arm in Fig. 1, left panel). Above this region in the the north-eastern part of the map one can see the radial velocity, which is related to a gas flow from the north. A part of the gas is swept from the galaxy to the west. The presented radial velocity component creates a radial magnetic field component. The differential rotation transforms this radial component into an azimuthal field and amplifies in this way the magnetic stripes occuring in this region. West of the galaxy center the radial velocity shows strong outflow of gas. The magnetic field there is amplified due to the same mechanism: strong compression causes the creation of a radial magnetic field component. The differential rotation transforms it into a toroidal component amplifying the already existing toroidal magnetic field. (cf. Fig. 1, left panel).

In order to present processes of magnetic field amplification connected with the re-accretion period, the left panel of Fig. 11 shows the radial velocity component vectors in the plane of the disk at the beginning of the re-accretion phase ($`t`$=650 Myr). The strong radial velocity connected with the wind is not longer present south-east of the galaxy. In the south-west the gas is moving outwards. This is the beginning of the shell expansion. The velocities of the back-falling counter-rotating arm are not visible, because this arm is not located within the disk plane. In the north the gas falls onto the galactic disk causing the formation of large number of magnetic maxima there.

The right panel of Fig. 11 presents velocity vectors in the $`X`$-$`Z`$ plane of the modeled box. At the east the magnetic field is amplified by the flow. In the lower right quarter of the image the velocity vectors are almost parallel to the disk, whereas those of the upper right quarter point towards the disk plane. Such a complicated velocity picture should be studied with a detailed analysis of the 3D flow, in order to investigate which mechanism (compression, shear, inflow, outflow) is responsible for the magnetic field amplification.

The most interesting and new result of our calculations is the growth of the total magnetic energy during the phase of re-accretion. We found that the total magnetic energy grows up to a few times its initial value (cf. Fig. 7) within the thick disk. Three processes are responsible for this growth: (i) compression of magnetic field by re-accreting gas, (ii) enhancement of magnetic field due to shearing motions, and (ii) a new kind of dynamo mechanism: the re-accreting material has important velocities in the radial and the vertical direction. Together with the differential rotation the galaxy’s velocity field becomes very complex. For $`t>700`$ Myr large gas loops are formed in the halo with heights up to several kpc. The magnetic field associated with them has an important poloidal component. The differential rotation in the disk creates the toroidal component of the magnetic field from the poloidal component (see Ruzmaikin et al. 1988). This may constitute a mechanism of the magnetic field generation. Thus, we find a new kind of the dynamo working in cluster galaxies, which have undergone a ram pressure event and are now in the phase of re-accretion.

The modeled polarization maps show a characteristic evolution of the distribution of the polarized intensity (PI). This evolution can be divided into two distinct phases:

* the compression phase and
* the re-accretion phase.

During the compression phase the maximum of the PI should coincide with the maximum of the gas surface density. This is not necessarily the case during the re-accretion phase. Thus, the polarized radio continuum emission could represent a tool to discriminate between different evolutionary stages of a cluster galaxy that undergoes or has undergone an ICM–ISM interaction. Together with the gas distribution and velocity field the PI distribution gives an important additional constraint on dynamical models. We plan to check this hypothesis in our future study. Observationally two Virgo galaxies that were observed in PI with the VLA at 2.8 cm and 6 cm (NGC4254: Soida et al. 1996 and NGC4522: Vollmer et al. in prep.). Both show asymmetric ridges of PI. Our future project is to compare the modeled PI maps to the high-frequency (Faraday rotation-free) polarized radio emission observations of Virgo cluster galaxies.

## 5 Conclusions

The evolution of the large-scale galactic magnetic field in a cluster galaxy undergoing an intracluster – interstellar interaction is investigated using 3D numerical simulations. In order to solve the induction equation (Eq. 1) we use velocity fields resulting from an N-body, sticky particle code including ram pressure effects (Vollmer et al. 2001).

We considered the general problem of magnetic field evolution in galaxies interacting with the ICM. Two experiments are performed with different values of magnetic diffusion. The classical $`\alpha `$-effect has not been involved in our computations. We perform the simulations of a galaxy, which is mildly HI deficient after the stripping event, as a representative model of the N-body calculations. This simulation presents the most characteristic features of magnetic field evolution during and after an ICM–ISM interaction.

The modeled magnetic fields are used to construct maps of the distribution of polarized radio continuum emission, which can be directly compared to observations.

We find that:

* interactions with the intracluster medium (i.e. the galaxy motion with respect to the ambient medium) lead to the formation of a heavy spiral arm along with a strong magnetic field concentration, on the disk windward side. On the other side of the disk the gaseous and magnetic arms stretch away from the disk center.
* Due to the gas re-accretion after the stripping event the vertical magnetic fields arise, forming fragments of a large-scale poloidal field, resembling those created by a dynamo process.
* Due to the ICM–ISM interaction, the total magnetic energy in cluster galaxies, which are affected by ram pressure, grows without any dynamo action. The interaction with the ambient gas may serve as a magnetic field amplification process in cluster spirals.
* The distribution of polarized radio continuum emission shows characteristic features during the galaxy evolution within the cluster. The comparison of maps of polarized radio continuum emission with our simulations will help to determine the evolutionary stage of a galaxy that undergoes an ICM–ISM interaction. It will also serve as an important additional constraint on dynamical modeling of ICM–ISM interactions.

###### Acknowledgements.

The authors express their attitude to Prof. Chantal Balkowski and Prof. Marek Urbanik for helpful discussions. K.O.-M. and B.V. thank Dr. Marian Soida for a detailed discussion of the polarization maps. This work was partly supported by a grant from the Polish Committee for Scientific Research (KBN), grant no. 0249/P03/2001/21.