text:
# Numerical study of multilayer adsorption on fractal surfaces

## Abstract

We report a numerical study of van der Waals adsoprtion and capillary condensation effects on self-similar fractal surfaces. An assembly of uncoupled spherical pores with a power-law distributin of radii is used to model fractal surfaces with adjustable dimensions. We find that the commonly used fractal Frankel-Halsey-Hill equation systematically fails to give the correct dimension due to crossover effects, consistent with the findings of recent experiments. The effects of pore coupling and curvature dependent surface tension were also studied.

The adsorption of inert molecules on self-similar fractal surfaces has been a subject of much interest and debate for more than a decade . It is generally believed that the adsoption isotherm beyond the monolayer regime follows a scaling form known as the fractal Frankel-Halsey-Hill (FHH) equation

$$N=N_o[\mathrm{ln}(P_o/P)]^\zeta ,$$
(1)
where $`N`$ is the number of molecules adsorbed, $`P`$ is the vapor pressure above the fractal surface and $`P_o`$ is the saturation vapor pressure at the temperature of interest. The prefactor $`N_o`$ depends on system specific parameters such as the linear size of the surface $`L`$, the density of the film $`\rho `$, etc. The central issue has to do with how the exponent $`\zeta `$ depends on the fractal dimension $`D`$ of the surface. Some authors have suggested that

$$\zeta =(3D)/3$$
(2)
while others have aruged that

$$\zeta =3D.$$
(3)
The difference has to do with what one believes to be the correct length scale $`\mathrm{}`$ for measuring the volume of the adsorbed film. If the van der Waals adsorption potential is dominant, the average film thickness $`t`$ on a planar surface is expected to be given by

$$\alpha a^3/t^3=k_BT\mathrm{ln}(P_o/P),$$
(4)
where $`a^3`$ represents the volume occupied by a liquid molecule and $`\alpha `$ is a dimensionless constant. If the surface tension $`\gamma `$ between the film and the vapor is dominant, the well known Kelvin equation sets a minimum radius curvature $`r`$ for the liquid-vapor interface, given by

$$2\gamma a^3/r=k_BT\mathrm{ln}(P_o/P).$$
(5)
In both cases, $`k_BT\mathrm{ln}(P_o/P)`$ represents the chemical potential difference $`(\mu _o\mu )`$ between the bulk liquid and vapor phases at temperature $`T`$. Since the volume of the wetting film on self-similar surfaces is expected to scale as $`\mathrm{}^{3D}`$ , using either $`t`$ or $`r`$ as the length scale $`\mathrm{}`$ leads to the two different predictions for $`\zeta `$ above. However, Kardar and Indekeu pointed out that self-affine surface roughness would yet lead to a different prediction and a correct theory should include both the substrate potential and the surface tension . In this paper, we report a computer simulation that is aimed at clarifying the applicability of the theory to real experiments.

Before describing our numerical study, it is necessary to review the status of the experiments briefly. We note that Eqs. (1)-(3) have been widely used to analyze adsorption isotherm data for a long time without a great deal of success. In many cases, Eq. (1) fit the data over less than one decade of length scales . In some cases, using Eq. (2) to interpret the exponent gave surface fractal dimension $`D<2`$, which is unphysical . While there have been attempts to compare different methods of analyzing the isotherm data , verification of the results using other techniques such as small-angle scattering were rare . Thus the applicability of Eq. (1) to real systems has remained unclear. To address this issue, we recently carried out a nitrogen adsorption study of three shale samples with fractal pore surfaces . $`D`$ for each sample was independently determined by small-angle neutron scattering (SANS) . We found that using Eqs. (1)-(3) to analyze the adsoprtion data consistently gave lower values for $`D`$ and there were also noticeable systematic errors in the fits. Fig. 1 replots the data from one sample to illustrate the problem. The adsorbed amount is expressed in terms of the the ratio of the film volume $`V_{ad}`$ to the total pore volume $`V_{full}`$. The length scale is expressed in terms of the Kelvin radius $`r`$ measured in units of $`a`$ (0.4 nm for N2 molecules) and it is calculated using Eq. (5) for a range of pressure $`P`$ relative to the saturation value $`P_o`$ for liquid nitrogen at 65 K, the temperature of the experiment. The fit was made over a two-decade range (0.4 - 50 nm) that coincides with the SANS data. Fig. 1a shows that using a linear fit on log-log scales and interpreting the slope ($`=0.46`$) according to Eq. (3) give $`D=2.54`$, a result that is significantly less than the value $`2.83`$ obtained by SANS. Using Eq. (2) for the slope would give an unphysical value of $`1.62`$. The horizontal dotted line at $`V_{ad}/V_{full}0.1`$ represents the monolayer capacity as determined by the standard BET analysis  and the vertical dotted line at $`r/a3`$ represents the onset of hystereses between adsorption and desorption (not shown). The latter is a signature of capaillary condensation in small pores due to surface tension . There is clearly a shoulder above the fitted line near this point where the adsorbed volume is still less than two monolayers. At higher pressure or long length scales where capillary condensation is expected to be more dominant, the data shows a slight upward curvature. Fig. 1b shows that fitting only the data above the shoulder to Eq. (1) with an additive background gives $`D=2.37`$, if Eq. (3) is used, which makes the disagreement with SANS worse. The same systematic deviations were seen in all three samples with different fractal dimensions . The most likely explanation is that Eqs. (1) - (3) represent asymtotic behavior that may not be observable in real systems. The crossover between van der Waals adsorption and capillary condensation may not permit such a simple analysis. The best way to test this hypothesis is to carry out a computer simulation in a simple model with a known fractal dimension.

The model we used to represent a fractal surface is an assembly of independent spherical pores with a power law distribution in the radii $`R`$ . If the number density follows $`g(R)R^{(D+1)}`$, then the surface area measured on a length scale $`\mathrm{}`$ is $`S_{\mathrm{}}_{\mathrm{}}R^2g(R)`$d$`R`$ $`\mathrm{}^{2D}`$, as expected of a fractal surface. Hence $`D`$ can be easily adjusted without any complication in the calculation. For each pore radius $`R`$, how the film thickness $`t`$ grows with increasing pressure $`P`$ up to the point of capillary condensation can be calculated numerically by combining Eqs. (4) and (5) to account for both the substrate potential and the surface tension. The resulting equation in dimensionless form is

$$\mathrm{ln}(P_o/P)=\frac{\mu _o\mu }{k_BT}=B\left[\frac{a^3R^3}{t^3(2Rt)^3}+\frac{Ca}{Rt}\right]$$
(6)
where $`B=8\alpha /k_BT`$ and $`C=\gamma a^2/4\alpha `$. The first term on the right comes from the van der Waals potential inside a spherical pore of radius $`R`$ at a normal distance $`t`$ from the surface and $`B`$ is a measure of the adsorption energy relative to the thermal energy. The second term is due to the surface tension of a meniscus of radius $`Rt`$ and $`C`$ is a measure of the relative strength of surface tension to the substrate potential. In real systems, $`B`$ and $`C`$ should be of order unity or greater when adsorption and capillary condensation occur.

The solution of Eq. (6) can be best understood using the graphical approach of Cohen et al . Fig.2 is a plot of the equation with $`B=C=1`$ and $`R/a=10`$. The solid line shows how the right-hand side of Eq. (6) varies with the film thickness $`t/a.`$ The contributions from the two separate terms are indicated by the dashed lines. We note that while the van der Waals term falls off with increasing $`t`$, the surface tension term rises with $`t`$. Hence their sum has a minimum at the point $`(t_c,P_c)`$. For any nonzero pressure $`P`$, a horizontal (dotted) line can be drawn through the point $`\mathrm{ln}(P_o/P)`$ on the vertical axis. For low pressures, the line intersects the solid curve at two points that are the solutions of Eq. (6), but the physical film thickness is given only by the smaller of the two $`t`$ values. With increasing pressure, the horizontal line is shifted downward and $`t`$ increases. At the critical pressure $`P_c`$, the line is tangent to the solid curve and it gives the maximum film thickness $`t_c`$, which is the threshold of condensation. At pressure above $`P_c`$, Eq. (6) has no real roots for $`t`$, i.e., the pore is filled. We can see from Fig. 2 and Eq. (6) that increasing $`B`$ would shift both of the dashed lines upward, thereby lowering $`P_c`$ without affecting $`t_c`$. Increasing $`C`$ would shift only the dashed line on the right, thereby lowering both $`P_c`$ and $`t_c`$. In our study, we varied both $`B`$ and $`C`$ over the range of $`0.110`$. Each choice of $`(B,C)`$ gives a different isotherm. We chose three different $`D`$ values ($`2.5`$, $`2.7`$ and $`2.85`$) to roughly match our earlier experiment .

To simulate an adsorption isotherm, we used $`200`$ pore radii distributed uniformly on a logarithmic scale over three decades: $`2R/a`$ $`2000`$. For each $`R`$, Eq. (6) can be solved to find the condensation pressure $`P_c`$. In practice, we found for each pressure $`P`$ the critical radius $`R_c`$. All the pores with $`R<R_c`$ are filled, and for those with $`R>`$ $`R_c`$, the film thickness $`t`$ was calculated. The total adsorbed amount was obtained by summing the contributons over the power-law distribution $`g(R).`$ Fifty pressure steps were used to obtain the entire isotherm and the pressure range was chosen such that the Kelvin radius $`r`$ in Eq. (5) spanned four decades ($`0.1a`$ to $`1000a`$), much wider than any real experiment. In Fig. 3a, the open circles depict the simulated isotherm with $`B=C=1`$ for a sample with $`D=2.85`$, a value that matches closely to the real sample in Fig.1. The simulation data exhibit a shoulder just above $`r/a1`$ in a manner similar to the real data in Fig. 1, but here we can see explicitly that it caused by the onset of capillary condensation in the smallest pore: above this threshold, the increase of pressure (or $`r`$) causes a rapid rise in condensed volumes (solid triangles) and a corresponding decrease of film volume (crosses). That the adsorption is separated into these two components explains why using either Eq. (2) or Eq.(3) alone to interpret the data over the entire range would result in systematic errors. At best, one should expect Eq. (2) to apply below the condensation threshold and Eq. (3) above it.

Figure 3b shows how the local slope of the simulated isotherm varies in a log-log plot. In the limit of small $`r`$ (Region I) where there is no condensation, the slope of $`0.351/3`$ corresponds to $`D2`$ according to Eq. (2). This proved to be a robust result for all the cases we examined. It implies that the film thickness $`t`$ given by Eq. (4) is not a length scale that could be used to probe fractal features, because $`t`$ can be used only when the system is free of capillary condensation, but that guarantees it to be below the lower cut-off of the fractal features (the smallest pore size). As a result, the only possiblilty is to use the Kelvin radius $`r`$ and analyze the data in the limit of large $`r`$, because the volume of the adsorbed film is a smaller percentage of the pore volume in the larger pores. In our simulation, we found that this approach works only for cases with large values of $`C`$ and only near the upper limit of the range of $`r`$ we used. In other cases, we still found fractal dimensions clearly below the true value. For example, the data in Fig. 3b fit a slope of $`0.25`$ in the large $`r`$ limit (Region III). It corresponds to $`D=2.75`$ according to Eq. (3), which is appreciably less than the input value of $`2.85`$ used for the simulation. In the middle region (Region II) where $`r`$ is comparable to the real data in Fig.1, the slope is higher ($`0.36`$) and it results in an even lower observed dimension ($`D=2.64`$). This trend of disagreement is consistent with what we found between our SANS and adsorption experiments . Similar behavior was observed in our simulated samples with $`D=2.7`$ and $`D=2.5`$. There is always an anomaly at the onset of capillary condensation, albeit the appearance varies somewhat. The reason is that the slope in Region I is always close to $`1/3`$ while that in Region III varies with $`D`$. So how the data join together in the middle is not universal. For the $`D=2.5`$ case, the slope in Region III tends to $`0.5`$, which is much higher than that in Region I. Hence the shoulder structure is masked by the intrinsic increase of the slope from $`1/3`$ to $`1/2`$, making the anomaly appear more like a kink. For $`D=2.7`$, although the slope in Region III tends to $`0.3`$, a value nearly identical to that in Region I, the shoulder appearance remains pronounced.

Our model was extended in an ad hoc manner to include some effects that occur in real systems. For example, the density of a liquid film is typically higher near the substrate and liquid/vapor surface tension is expected to be stronger for concave surfaces with small radius of curvature. Both of these effects would increase the second term in Eq. (6) relative to the first term. We investigated these effects qualitatively by letting $`C`$ be a function $`C(R)=C_1+C_2\mathrm{exp}(R/\lambda )`$. From the data shown in Ref. , we estimated that $`\lambda 5.4a`$. In addition, for any real fractal surface, features of different sizes are not completely independent. Condensation in small pores can trigger the same to occur in neighboring large pores. We modelled such pore coupling effects by letting pores with radius $`R>R_c`$to condense with a probability $`K(R_c/R)^q`$, i.e., less likely for larger pores. Fig. 3c shows that with $`C_1=C_2=1,`$ $`K=0.3`$ and $`q=4`$, the isotherm in Fig. 3b acquires a more pronounced shoulder. The curvature in Region II also changes from concave downward to concave upward. A power law fit gives a slightly lower dimension ($`D=2.60`$). These results are in qualitative agreement with the real data in Fig. 1. They can be explained by the fact that the increased surface tension makes condensation in small pores to occur at a lower pressure and pore coupling makes some of the larger pores to follow. The net result is to enhance the anomaly at the onset of capillary condensation.

Despite the fact that our model is far from being a realistic representation of any real system, it clarifies a long-standing ambiguity associated with the interpretation of adsorption isotherm on fractal surfaces. By varying the parameters $`B`$, $`C`$ and $`D`$, we are able to see explicitly how the crossover between substrate driven adsorption and surface tension driven condensation prevents one from using the FHH equation to obtain the correct fractal dimension under common experimental conditions. We are also able to produce simulated isotherms that resemble the real data and understand the cause of their essential features. We hope that these understandings, though only qualitative, will serve as useful guides for interpreting adsorption isotherms with the FHH equation.

We thank J. Machta and R. Guyer for helpful discussions. This work was supported by The Petroleum Research Fund administered by the American Chemical Society under Grant Nos. 32191-AC2-SF98 and 33549-B9, and NSF Grant No. CTS-9803387.


format: text

text:
# RPA quasi-elastic responses in infinite and finite nuclear systems

## I INTRODUCTION

The quasi-elastic excitation of the nucleus is particularly interesting for nuclear structure studies because of the interplay between single particle and many-body effects. The excitation energies characterizing the quasi-elastic peak are well above the nucleon emission threshold, therefore, one or more nucleons are ejected from the nucleus. Before being ejected, any nucleon can interact with the others nucleons composing the nucleus. Considering this effect in a finite nucleus formalism is a difficult task because the nuclear final states are described in terms of their total angular momentum and the presence of the continuum implies the sum on a large amount of possible configurations.

In a nuclear matter formalism one takes advantage of the translational invariance to simplify the description of the final state of the hadronic system. The use of the infinite system formalism is rather appropriate for the quasi-elastic electron excitation, since excitation energies and transferred momenta are such that the excitation process is well localized within the nucleus and collective surface excitations are negligible. In any case additional approximations (say, a variable Fermi momentum or the local density approximation) have been considered in the literature to simulate finite size effects.

We can classify the necessary ingredients to describe the process in three main issues: the single particle basis, the inclusion of initial and/or final state interactions and the operator describing the action of the external probe. Regarding the first point we have already said that some work treats the nucleus as a finite system -, while other ones use a nuclear matter approach -.

Referring to the second point, one of the simplest approaches to include initial and final state interactions is the Random Phase Approximation (RPA), with or without exchange terms (the last approach named ring approximation). In the RPA -, - one particle-one hole ($`1p1h`$) excitations are summed up to infinite order. One step further in complexity is the so-called Second RPA (SRPA) ,  which, in addition to the $`1p1h`$ excitations, considers also those generated by $`2p2h`$. The Extended RPA (ERPA) , , ; contains ground state correlations beyond RPA. The Green function approach of Ref.  is based on a philosophy similar to that of the SRPA. In this approach the relationship between forward virtual Compton scattering and inclusive electron scattering is used to construct a one-body approximation to quasi-elastic electron scattering. The role of the short-range correlations has been investigated in the framework of the correlated basis function theory ,.

Concerning the third point, the external operator is usually represented by electromagnetic one-body operators, but in the transverse channel it is important to include two-body terms as the Meson Exchange Currents (MEC) , , , , and the excitation of the virtual or real $`\mathrm{\Delta }(1232)`$ resonance -.

In the present paper we investigate the sources of some inconsistencies between finite and infinite nucleus calculations. For example it seems that, once the residual interactions has been fixed, the effects of the ring approximation are larger in nuclear matter than in finite nuclei. To simulate the finite size effects in the quasi-elastic peak, we have done infinite systems calculations with a diffused Fermi surface. Then, we have studied the effects on the responses of the various characteristics of the interaction, such as the different channels, the range and the density dependence. The infinite systems results have been compared with the calculations done for the 12C and 40Ca nuclei, two doubly magic nuclei with the same number of protons and neutrons. The agreement we have finally obtained between the two kinds of calculations is satisfactory, even for the nucleus 12C, which is supposed to be relatively light to be well represented as an infinite system of nucleons.

Encouraged by this result, we have evaluated the response functions within the nuclear matter RPA framework of Ref.  to investigate the effects of the exchange terms neglected in the ring approximation.

The paper is organized as follows. In Sec. II we briefly present the effective theories used to perform the calculations. In Sec. III we compare the results obtained in finite nuclei and nuclear matter both for the free and the ring responses, and we present the nuclear matter RPA responses. Finally, in Sec. IV we draw our conclusions.

## II FORMALISM

The response function for inclusive quasi-elastic electron scattering is given by,

$$R(𝒒,\omega )=\frac{1}{\pi }Im<0|𝒪^{}(𝒒)G(\omega )𝒪(𝒒)|0>$$
(1)
where $`\omega `$ represents the excitation energy and $`𝒒`$ the three-momentum transferred by the electron. The nuclear ground state is denoted as $`|0>`$, while $`𝒪(𝒒)`$ is the excitation operator and $`G(\omega )`$ the polarization propagator,

$$G(\omega )=\frac{1}{\omega H+i\eta }\frac{1}{\omega +H+i\eta }$$
(2)
where $`H`$ is the nuclear Hamiltonian. We introduce the operator $`P`$, which projects onto $`1p1h`$ configurations. The way of building these configurations depends on the particular choice of the single particle basis. In nuclear matter the single particle wave functions are plane waves, while in a finite nucleus they are eigenfunctions of the one-body Schrödinger equation for a mean field potential. In our case we use a real, spherical, Woods-Saxon potential.

The response function generated by $`1p1h`$ excitations can be expressed as:

$$R_{PP}(𝒒,\omega )=\frac{1}{\pi }Im<0|𝒪^{}(𝒒)PG(\omega )P𝒪(𝒒)|0>$$
(3)
The evaluation of $`R_{PP}`$ is not straightforward since the nuclear Hamiltonian is in general not diagonal in the $`1p1h`$ basis. The solution of the problem is the RPA response. In the present work, we shall compare finite nucleus results obtained within the ring approximation of Ref. , with the corresponding nuclear matter responses.

In most calculations done in infinite systems the Fermi surface which separates, in momentum or energy space, the hole from the particle states, is a sharp step function. In the next section, we shall show that this approximation reproduces rather well the position and the shape of the free responses calculated in finite nuclei. However, the ring responses are appreciably different when evaluated in nuclear matter and in finite nuclei. The original motivation of the present work was to explain this discrepancy. The first attempt to improve our nuclear matter model, is to replace the step function representing the momentum distribution of particles and holes, by a more realistic one. For the holes with momentum $`𝒉`$ we make the substitution,

$$\theta (k_F|𝒉|)n(h)$$
(4)
where $`k_F`$ is the Fermi momentum, $`h|𝒉|`$, and

$$n(h)=\frac{1}{1+e^{(hk_F)/a}}$$
(5)
being $`a`$ a constant to be adjusted. An analogous expression can be obtained for particles.

In the next section we shall compare results obtained with this smoothed Fermi surface in ring approximation with those obtained in finite systems calculations. With the same smoothed momentum distribution we have done RPA calculations following the computational scheme developed in Ref.  which we briefly describe. In nuclear matter, direct RPA terms can be summed up to infinite order (ring series), but, in general, it is not possible to find a closed form to sum all the exchange terms. Normally these last terms are perturbatively considered and, for numerical reasons, this is done up to the second order. On the other hand, it possible to make the full summation of the exchange terms when a contact interaction is used. Full sums are also possible for separable interactions. To exploit this feature we rewrite the residual interaction $`V`$, as:

$$V=V_{contact}+\stackrel{~}{V}$$
(6)
The $`V_{contact}`$ term is a contact interaction conveniently chosen to make $`\stackrel{~}{V}`$ small. For $`V_{contact}`$ both direct and exchange terms are summed up to infinite order while $`\stackrel{~}{V}`$ is perturbatively considered up to the second order which we found to account reasonably well for the whole sum. Interference terms between $`V_{contact}`$ and $`\stackrel{~}{V}`$ are included up to infinite order in $`V_{contact}`$ and up to second order in $`\stackrel{~}{V}`$. A more detailed description of the method can be found in Ref. .

The finite nucleus calculations have been done within the Fourier-Bessel computational scheme adopted in Ref. . The single particle basis is generated by using a Woods-Saxon mean field whose parameters have been fixed to reproduce experimental rms radii and single particle energies of the bound states close to the Fermi surface. By neglecting the exchange diagrams, the RPA equations are rewritten in terms of local density functions which are expanded on a Fourier-Bessel basis. In this manner the problem to be solved, for every value of the excitation energy, is the diagonalization of a matrix whose dimensions are four times the number of the Fourier-Bessel expansion coefficients. More details about the method can be found in Ref. .

## III RESULTS

In this section we compare finite nucleus results with those of nuclear matter. This is done for two nuclei: 12C and 40Ca. The 12C nucleus is perhaps too light to be appropriately described in terms of nuclear matter. On the other hand we wanted to test our model also in extreme situations and, last but not least, the finite nucleus calculations are much less involved than in the 40Ca case.

The transferred momenta analyzed are $`q=400`$ MeV/c and $`q=500`$ MeV/c. These values are sufficiently large to eliminate the presence of collective surface vibrations, and at the same time, sufficiently small to require a limited number of partial waves. In the finite nucleus calculations we sum multipole excitations up to angular momentum J=12 .

The first step of our calculations consists in fixing the values of $`k_F`$ and $`a`$ in eq. (5) to reproduce the finite nucleus free responses. For 12C we have obtained the values $`k_F`$=0.85 fm-1 and $`a`$=0.20 fm-1 , and for 40Ca, $`k_F`$=1.0 fm-1 and $`a`$=0.17 fm-1. The comparison with the finite nuclei responses is shown in Figs. 1 and 2, where we have also added the results obtained with a step function Fermi surface represented by the dashed lines. In this last case we have used the procedure of Ref.  to fix the value of the Fermi momentum obtaining the values $`k_F`$=1.09 fm-1 for 12C and $`k_F`$=1.19 fm-1 for 40Ca. These values are noticeably different from those fixed by the smoothed momentum distribution. As expected, $`k_F`$ increases with increasing mass number. We succeeded in obtaining a satisfactory agreement with the finite nucleus responses, especially for 40Ca, in the case of a diffused Fermi surface. As expected the high energy tail of the finite nucleus responses can be reproduced only by the calculations with a diffused Fermi surface.

The comparison between the various responses calculated in ring approximation is shown in Figs. 3 and 4 for the 12C and 40Ca nuclei respectively. In these figures the meaning of the symbols is analogous to that of the previous figures. The interaction used in these calculations is the finite range polarization potential utilized also in Ref. . The nuclear matter calculations have been done with the values of $`k_F`$ and $`a`$ previously fixed.

Two observations should be done about these results. A first one is about the fact that, in general, the full curves reproduce better the finite nucleus results than the dashed ones. In a second place we observe that the longitudinal responses are better reproduced than the transverse ones.

The first observation induces to conclude that the sharp Fermi surface, even with an effective value of the Fermi momentum, is unable to reproduce the finite nucleus results. This feature does not depend from the residual interaction, as we show in Fig. 5 where the various 12C responses have been calculated, always in ring approximation, with a zero range Migdal interaction. We have used the following values of the parameters of this force: $`f_0=386`$, $`f_0^{}=289.5`$, $`g_0=106.2`$ and $`g_0^{}=135.1`$, expressed in MeV fm3 units. These values have been chosen to magnify some of the effects we want to discuss. Specifically, the big value of $`f_0`$, ten times larger than the one normally used , enhances the difference between the nuclear matter calculations in the longitudinal response. Here it becomes more evident the poor quality of the sharp Fermi surface calculations. Our diffused Fermi surface calculation is able to reproduce rather well the finite nucleus results even in these extreme conditions.

Concerning the observation that, the longitudinal responses are always better reproduced than the transverse ones, we have verified that this fact is due to the differences between infinite and finite systems. The isospin channel of the force does not contribute in ring approximation calculations of nuclear matter transverse responses, while it does in finite nuclei calculations. The effects of this difference become evident by comparing the transverse responses of Fig. 5 with Fig. 6. In this last figure the 12C transverse responses have been calculated with the contact interaction above described but without the isospin channel of the force, i.e. by setting $`f_0^{}=0`$. The agreement between the finite nucleus responses and those obtained with the diffused Fermi surface is comparable with that obtained in the longitudinal case.

We should remark the fact that the polarization potential has a density dependence in the scalar and isospin channels. In the nuclear matter calculations we have used the force parameters defined for the nuclear interior. We have checked the sensitivity of the results by switching off the density dependence in the finite nuclei calculations and comparing with the responses evaluated with the full interactions. The differences found between these two calculations are of the same order of the differences between nuclear matter and finite nuclei longitudinal responses. We conclude that for the calculations of quasi-elastic responses the density dependence of the force is not important.

In Figs. 7 and 8 we compare the RPA nuclear matter responses (full lines) with those evaluated in ring approximation (dashed lines). Both calculations have been done by using the polarization potential and the diffused Fermi surface. This comparison shows the effects of the exchange diagrams evaluated in RPA and neglected in ring approximation. For the particular interaction used these effects are noticeable, especially in the longitudinal responses and for low values of the momentum transfer.

In the same figures we also present the free responses (dotted lines) and the 12C and 40Ca experimental points. Our results show that the major source of disagreement between free responses and data is produced by correlations beyond the RPA. In Ref.  the role of the final state interactions and of the effective mass was pointed out. The inclusion of these two effects within a simplified model produces a good agreement with the 40Ca data . The same model is however unable to explain the 12C transverse response data.

## IV CONCLUSIONS

In the present work we have compared nuclear matter and finite nuclei quasi-elastic responses induced by electron scattering. We have proposed a computational scheme which is able to reproduce the main features of the finite nucleus and, at the same time, it has the computational advantages of nuclear matter. This has been achieved by using a diffused momentum distribution of particles and holes in nuclear matter calculations. We have described the momentum distribution with a simple Fermi function depending from two parameters whose values have been adjusted to reproduce the finite nuclei free responses.

With this simple model we have calculated the 12C and 40Ca responses for different values of the momentum transfer in ring approximation. The agreement with the finite nucleus calculations is excellent for the longitudinal responses. We traced the source of the small differences found in the transverse responses to the isospin part of the residual interaction which is not active in nuclear matter calculations.

Using the same residual interaction, the polarization potential, we have done RPA calculations, i. e. we have also considered the exchange diagrams. The aim of this calculation was to test its feasibility. Finite nuclei RPA calculations have been done -, but they require a large computational effort, while our method is simpler and more suitable to be used for those extensions beyond RPA which are necessary to describe the experimental points.

With respect to this last point some words of caution are necessary to avoid double counting. A first one is about the residual interaction which in RPA calculations cannot be the polarization potential constructed in Ref.  to consider in average manner the exchange diagrams. A second warning is about the diffuseness of the Fermi surface which is partially produced by correlations effects beyond RPA.

Our approach should be compared with the most commonly used method to account for finite nucleus effects in nuclear matter calculations: the local density approximation (LDA). Unfortunately this comparison is not straightforward. In LDA the responses are calculated for several values of the Fermi momentum and then they are appropriately averaged to obtain the final result. In the averaging procedure the way how the responses for each $`k_F`$ are weighted, differs from nucleus to nucleus. However, for each $`k_F`$ a sharp Fermi surface is employed and this implies that the second term in the right hand side of Eq. (2) does not contribute. This terms, however, gives a contribution when a smooth Fermi surface is considered. This contribution is normally small but it becomes appreciable when $`\omega `$ 20 MeV . This is the formal difference between our method and the LDA. From the pragmatical point of view the poor quality of LDA results in reproducing finite nuclei quasi-elastic responses has already been pointed out in Ref. , where the effective momentum approximation was found to be better. Here we have shown the superiority of the diffused Fermi surface scheme even with respect to the effective Fermi momentum approximation. One should also consider that also from the computational point of view our scheme is superior to the LDA. In the last case one has to calculate nuclear matter responses for different values of $`k_F`$, while in our case a single calculation is necessary. For the fast calculations done in ring approximation, this difficulty is irrelevant, but it becomes an handicap for RPA calculations, or for more elaborated ones like SRPA or ERPA.

###### Acknowledgements.

We thank A.M. Lallena for useful discussions.


text:
# Charmonium in finite temperature lattice QCD

## 1 Introduction

It is widely believed that the Quantum Choromodynamics (QCD) exhibits a phase transition at some temperature $`T_c`$, and quarks and gluons confined in the low temperature phase are liberated to form “quark gluon plasma”. At the beginning of 2000, CERN reported that the QGP state had been created in the heavy ion collision experiment. In this experiment, the $`J/\psi `$ suppression  is regarded as a key signal of QGP formation. More reliable inspection will be performed in RHIC project at BNL.

On the theoretical side, in spite of various approaches, we are still far from definite understanding of charmonium properties near the transition and the fate of the hadronic states in the plasma phase. We investigate this problem using lattice QCD, which enables us to incorporate the nonperturbative effect of QCD.

## 2 Our approach

The difficulties in studying the heavy quarkonium state at finite temperature can be categorized into two classes: (i) difficulties to obtain detailed information on the heavy quarkonium correlators at finite temperature, and (ii) analyses of the obtained correlators. Although it is difficult to overcome the latter, the former is possible to be solved by employing the techniques developed in recent years.

At high temperature, sufficient information on the temporal correlator requires high resolution in the temporal direction. To achieve large temporal cutoff with limited computational resources, we employ the anisotropic lattice, on which the temporal lattice spacing is finer than the spatial one . Calculations with a heavy quark need further implementation. The standard lattice action contains large systematic uncertainty for the quark mass larger than or comparable to the lattice cutoff. This problem was reexamined with the light of effective theoretical interpretation of the Wilson (and clover) quark action by FNAL group . Their approach is well-matched with to the anisotropic lattice, since it naturally introduces the anisotropy in the spatial and the temporal hopping parameters. We employ this approach to describe the heavy quark on present lattices, which have not sufficiently large spatial cutoff compared with the charm quark mass.

The latter problems, (ii), are much more difficult to be overcome. This is because one is inevitably enforced to extract the spectral properties from the data at the short temporal distance, where high frequency component of the dynamics is significant. At present, there is no well-established way of lattice simulations to attack the spectroscopy at $`T>0`$. For this reason, lattice studies of hadrons at finite temperature have been argued on the spatial correlation (screening mass) . Recently, QCD-TARO Collaboration (including three of us) analyzed the space-time structure of the mesonic correlators in the temporal direction to examine the temperature effect on them . They caught a sign of the bound state even in the deconfined phase, as well as an indication of the chiral symmetry restoration. Their strategy seems one of the best approaches also to the present subject, and is expected to give us novel information on the fate of charminium near the phase transition. Therefore we proceed the analysis with the following steps.

* At $`T=0`$ we construct good meson operators which has large overlap with the state of interest.
* Temperature effects are observed in the meson correlator between the operators constructed in the previous step.
* We define the “wave function” in the Coulomb gauge at finite temperature, and observe its $`t`$-dependence below and above $`T_c`$ to discuss the possibility of the bound state at $`T>T_c`$.

## 3 Fermilab action on anisotropic lattice

The quark action we use is in the following form.

$$S_F=\underset{x,y}{}\overline{q}(x)K(x,y)q(y),$$
(1)
$`K(x,y)=\delta _{x,y}`$
$`\kappa _\tau \{(1\gamma _4)T_{+4}+(1+\gamma _4)T_4\}`$
$`\kappa _\sigma _i\{(r\gamma _i)T_{+i}+(r+\gamma _i)T_i\}`$
$`(\kappa _\sigma c_Eg\sigma _{4i}F_{4i}(x)+r\kappa _\sigma c_B{\displaystyle \frac{1}{2}}g\sigma _{ij}F_{ij}(x))\delta _{x,y},`$ (2)
where $`T_{\pm \mu }=U_{\pm \mu }(x)\delta _{x\pm \widehat{\mu },y}`$. The bare anisotropy $`\gamma _F`$ of the fermion field appears in the ratio of the spatial and the temporal hopping parameters, $`\kappa _\tau =\gamma _F\kappa _\sigma `$. The Wilson parameter in this form is $`r=1/\xi `$, where $`\xi `$ is the renormalized anisotropy: $`\xi =a_\sigma /a_\tau `$. The hopping parameters are related to the bare quark mass as $`\kappa _\sigma =1/2(m_0+\gamma _F+3r)`$. We also define $`\kappa `$ with $`1/\kappa =1/\kappa _\sigma 2(\gamma _F+3r4)`$ as the parameter which plays the same role as $`\kappa `$ on the isotropic lattice. As the field strength $`F_{\mu \nu }`$, we use the standard clover-leaf definition. At the tree level, the coefficients of the clover terms, $`c_E`$ and $`c_B`$, are unity.

This action is obtained by generalizing the form in  to the anisotropic lattice, and by setting $`r=1/\xi `$. This choice of the Wilson term is adopted so that the action restores the axial interchange symmetry in the physical unit as $`m_00`$. Therefore we can use the same form in the light quark region as the $`O(a)`$ improved Wilson quark action . This is the reason that our form of the action is different from the form used in .

On the anisotropic lattice, we need additional tuning of the parameters called as “calibration”, to ensure the renormalized anisotropy is same for the gauge and the fermion fields . We define the fermionic anisotropy $`\xi _F`$ using the dispersion relation of the pseudoscalar and the vector mesons and tune the value of $`\gamma _F`$ so that $`\xi _F`$ coincides with the gauge field anisotropy $`\xi `$.

We incorporate the mean-field improvement : $`U_iU_i/u_{0\sigma }`$ ($`i=`$1,2,3) and $`U_4U_4/u_{0\tau }`$ with the mean-field values of the spatial and the temporal link variables, $`u_{0\sigma }`$ and $`u_{0\tau }`$ defined in the Landau gauge.

## 4 Result at zero temperature

Numerical calculation is done on the lattice of sizes $`12^3\times N_t`$ where $`N_t=72,20,16`$ and $`12`$, with the standard Wilson gauge action at the coupling $`\beta =5.68`$ and the bare anisotropy $`\gamma _G=4.0`$, in the quenched approximation . On $`N_t=72`$ lattice, the renormalized anisotropy and the spatial lattice spacing from the string tension are determined as $`\xi =5.3(1)`$ and $`a_\sigma ^1=0.85(3)`$ GeV. $`N_t=20`$, $`16`$ and $`12`$ lattices correspond to the temperatures $`T`$ 0.93 $`T_c`$, 1.15 $`T_c`$ and 1.5 $`T_c`$, respectively. We calculate the meson correlators at $`T=0`$ with the three sets of the parameters ($`\kappa `$, $`\gamma _F`$) listed in Table 1.

Each value of $`\gamma _F`$ is determined by the calibration procedure mentioned above. Observing the resultant meson masses summarized in Table 1, ($`\kappa `$, $`\gamma _F`$) $`=`$ (0.1013, 3.703) roughly corresponds to the charm quark. In the successive analysis, we use this set of parameters.

To optimize the meson correlator at $`T=0`$, we employ the variational analysis. This is implemented by diagonalizing the correlator matrix

$$C_{ij}(t)=_\stackrel{}{x}O_i(\stackrel{}{x},t)O_j^{}(0).$$
(3)
The operator $`O_i(\stackrel{}{x},t)`$ is defined with the smearing function $`\phi `$ in the Coulomb gauge as

$$O_i(\stackrel{}{x},t)=_\stackrel{}{y}\overline{q}(\stackrel{}{x}+\stackrel{}{y},t)\phi _i(\stackrel{}{y})\mathrm{\Gamma }q(\stackrel{}{x},t),$$
(4)
where the $`4\times 4`$ matrix $`\mathrm{\Gamma }`$ specifies the quantum number of the meson. As the smearing function $`\phi _i`$, we use the eigenfunction obtained by solving the Schrödinger equation for the $`S`$ state,

$`\left[{\displaystyle \frac{1}{2m_R}}{\displaystyle \frac{d^2}{dr^2}}+V(r)\right]y(r)`$ $`=`$ $`Ey(r),`$
$`y(r)`$ $`=`$ $`r\phi (r),`$ (5)
where $`V(r)`$ is the static quark potential obtained on the same lattice, and $`m_R=1.5/2`$ \[GeV\] is used.

Figure 1 shows the effective masses of the diagonalized correlators with two largest eigenvalues (smallest masses). Both correlators show plateaus beyond $`t10`$, though the statistical fluctuation is large for the correlator corresponding to the first excited state. We fit these data at $`t=16`$$`36`$ and find the meson masses in 2$`S`$ states as $`m_{Ps(2S)}=0.7677(77)`$ and $`m_{V(2S)}=0.7794(88)`$ in the lattice unit. Although these give 1$`S`$-2$`S`$ splitting smaller than the experimental value, considering the large systematic uncertainty on this coarse lattice, they are not inconsistent results. We regard the operator corresponds to the ground state in this analysis as the optimized meson operator at $`T=0`$.

## 5 Result at finite temperature

We start the argument at $`T>0`$ with the temperature dependence of the correlator between the optimized operators at $`T=0`$ difined in the previous section. Figure 2 shows the effective mass plot at finite temperatures in the vector channel. Below $`T_c`$ ($`N_t=20`$), the plateaus of the effective mass seem to appear. Although the effective mass of the ground state is slightly larger than that at $`T=0`$, it is difficult to identify the plateau precisely and determine the mass quantitatively with present statistics. More detailed analysis with higher statistics may open a stage to discuss the potential mass shift of charmonium near to the $`T_c`$ .

Above $`T_c`$ ($`N_t=16`$ and $`12`$), the correlators do not show any clear plateau in the effective mass plot. They decrease more rapidly than expected from the behavior at $`N_t=20`$ in both of the correlators correspond to the ground and the first excited states at $`T=0`$. The Ps channel shows similar feature, although the decrease of the effective mass is slightly milder. These behavior at least signal significant change in the nature of the correlators when the system crosses $`T_c`$. In , the (effective) masses increase as $`T`$ in the pseudoscalar and the vector channels. The observed behavior in present work, however, shows qualitatively different nature of the correlators.

To discuss the spatial correlation between quark and antiquark at $`T>0`$, we observe the “wave function” normalized at the spatial origin,

$$\varphi (\stackrel{}{r},t)=w_\mathrm{\Gamma }(\stackrel{}{r},t)/w_\mathrm{\Gamma }(\stackrel{}{0},t),$$
(6)
$$w_\mathrm{\Gamma }(\stackrel{}{r},t)=_\stackrel{}{x}\overline{q}(\stackrel{}{x}+\stackrel{}{r},t)\mathrm{\Gamma }q(\stackrel{}{x},t)O^{}(0).$$
(7)
If there is no bound state, like as with the free quark propagators, this wave function $`\varphi (\stackrel{}{r},t)`$ broadens with $`t`$ for any source smearing function. In contrast, the existence of the stable shape of $`\varphi (\stackrel{}{r},t)`$ give us a hint on the existence of the bound state.

Figure 3 shows the $`t`$-dependence of the wave function normalized at the spatial origin in the vector channel at $`N_t=`$20 and 12. The source function of these correlators are taken to be wider than the optimized one, to show clearly that $`\varphi (\stackrel{}{r},t)`$ keeps narrow shape. At the both $`N_t`$, the wave function gradually decrease as $`t`$, in contrast to the case with the free quark poropagators shown together. This indicates that the quark and antiquark tend to stay together even in the deconfined phase, up to 1.5 $`T_c`$.

These two observations are interesting, and at the same time puzzling. Obviously, the phase transition causes significant change in the charmonium correlators, especially in the vector channel. On the other hand, the result of the wave function suggests the persistence of the bound states. One natural picture is that the mesonic spectral function still have a peak with rather large width above $`T_c`$. This explains both of our observations. Direct extraction of the spectral function from the lattice data  may give us further information. For the definite understanding of the fate of chamonium at the phase transition, we need more systematic studies and development of procedures.

The simulation has been done on Intel Paragon XP/S and NEC HSP at INSAM, Hiroshima University. This work is supported by the Grant-in-Aide for Scientific Research by Monbusho, Japan (No.10640272,No.11440080)


text:
# Domain State Model for Exchange Bias

## Abstract

Monte Carlo simulations of a system consisting of a ferromagnetic layer exchange coupled to a diluted antiferromagnetic layer described by a classical spin model show a strong dependence of the exchange bias on the degree of dilution in agreement with recent experimental observations on Co/CoO bilayers. These simulations reveal that diluting the antiferromagnet leads to the formation of domains in the volume of the antiferromagnet carrying a remanent surplus magnetization which causes and controls exchange bias. To further support this domain state model for exchange bias we study in the present paper the dependence of the bias field on the thickness of the antiferromagnetic layer. It is shown that the bias field strongly increases with increasing film thickness and eventually goes over a maximum before it levels out for large thicknesses. These findings are in full agreement with experiments.

For a ferromagnet (FM) in contact with an antiferromagnet (AFM) a shift of the hysteresis loop along the magnetic field axis can occur which is called exchange bias (EB). Usually, EB is observed after cooling the system with the FM magnetized in saturation below the Néel temperature $`T_\mathrm{N}`$ of the AFM or after cooling the entire system in an external magnetic field. Although this effect is well known since many years its microscopic origin is still discussed controversially. For a review see a recent article by Nogués and Schuller .

In a previous Letter  we have reported on EB observed experimentally in Co/CoO bilayers as a function of volume defects in the antiferromagnet. Of particular importance in this study was the observation that it is possible to strongly influence EB in Co/CoO bilayers by diluting the antiferromagnetic CoO layer, i. e. by inserting nonmagnetic substitutions (Co1-xMgxO) or defects (Co1-yO) not at the FM/AFM interface, but rather throughout the volume part of the AFM. While the undiluted samples show only a very small EB, dilution increases EB dramatically. Since for all samples investigated a 0.4 nm thick CoO layer with minimum defect concentration was placed at the interface the observed EB is primarily not due to disorder or defects at the interface. Rather, the full antiferromagnetic layer must be involved and we have argued that in our systems EB has its origin in a domain state in the volume part of the AFM which triggers the spin arrangement and thus the FM/AFM exchange interaction at the interface. Indeed, in diluted antiferromagnets when cooled in external fields metastable domains occur carrying a surplus magnetization and having a very slow dynamics (for reviews see ). These domains are frozen to a large extend during hysteresis cycles and their frozen magnetization is the origin of EB. This domain state model for EB was supported by large scale Monte Carlo simulations performed at finite temperatures.

To gain additional evidence for the domain state model in the present paper we will concentrate on the dependence of EB on the thickness of the AFM. The system consists of a FM monolayer exchange coupled to a diluted AFM layer with $`t`$ monolayers. The FM is described by a classical Heisenberg model with vector spins $`\underset{¯}{S}_i`$ and exchange constant $`J_{\mathrm{FM}}`$. The dipolar interaction is approximated by an additional anisotropy term (anisotropy constant $`d_x=0.1J_{\mathrm{FM}}`$) which mimics the shape anisotropy leading to a magnetization which is preferentially in the $`yz`$-plane. Also, we introduce an easy axis in the FM ($`z`$-axis, anisotropy constant $`d_z=0.1J_{\mathrm{FM}}`$) in order to obtain well defined hysteresis loops. $`d_z`$ sets the Stoner-Wohlfarth limit of the coercive field, i. e. the zero temperature limit for magnetization reversal by coherent rotation ($`B_{\mathrm{SW}}=2d_z`$, in our units, for a field parallel to the easy axis). In view of the rather strong anisotropy in CoO we assume an Ising Hamiltonian for the AFM. Thus the Hamiltonian of our system is given by

$`=`$ $``$ $`J_{\mathrm{FM}}{\displaystyle \underset{i,j}{}}\underset{¯}{S}_i\underset{¯}{S}_j{\displaystyle \underset{i}{}}\left(d_zS_{iz}^2+d_xS_{ix}^2+\underset{¯}{S}_i\underset{¯}{B}\right)`$ (1)
$``$ $`J_{\mathrm{AFM}}{\displaystyle \underset{i,j}{}}ϵ_iϵ_j\sigma _i\sigma _j{\displaystyle \underset{i}{}}Bϵ_i\sigma _i`$ (2)
$``$ $`J_{\mathrm{INT}}{\displaystyle \underset{i,j}{}}ϵ_i\sigma _iS_{jz}.`$ (3)
with the antiferromagnetic nearest-neighbor exchange constant $`J_{\mathrm{AFM}}<0`$ and the effective in-plane magnetic field $`\underset{¯}{B}=B\underset{¯}{\overset{^}{z}}+B_y\underset{¯}{\overset{^}{y}}`$. The values of the magnetic moments are incorporated in $`B`$ and $`d`$, respectively, so that the quantities $`\underset{¯}{S}_i`$ denote unit vectors and $`\sigma _i=\pm 1`$ Ising spin variables. A fraction $`p`$ of the sites of the lattice is left without a spin (quenched disorder: $`ϵ_i=0,1`$). For the exchange constant of the AFM which mainly determines its Néel temperature (also depending on the dilution, of course) we set $`J_{\mathrm{AFM}}=J_{\mathrm{FM}}/2`$. There seems to be some evidence that the exchange coupling between Co and CoO is ferromagnetic  but its strength is not known experimentally. Therefore, we assume in our simulations a ferromagnetic coupling with ($`J_{\mathrm{INT}}=J_{\mathrm{AFM}}`$).

We use Monte Carlo methods with a heat-bath algorithm and single-spin flip methods for the simulation of the model explained above. The trial step of the spin update is a small variation around the initial spin for the Heisenberg spins and — as usual – a spin flip for the Ising spins . We perform typically 40000 Monte Carlo steps per spin (MCS) for a complete hysteresis loop. To observe the domain structure of the AFM we have to guarantee that typical length scales of the domain structure fit into our system. For the parameter values used in this simulation this is achieved for systems of lateral extension $`128\times 128`$.

In the simulations the system is cooled from above to below the ordering temperature of the AFM in an applied external cooling field $`\underset{¯}{B}=B_c\underset{¯}{\overset{^}{z}}`$ with $`B_c=0.25J_{\mathrm{FM}}`$ . The FM is then long-range ordered and its magnetization is practically homogeneous resulting in a nearly constant exchange field for the AFM monolayer at the interface. When the desired final temperature is reached a magnetic field $`\underset{¯}{B}=B\underset{¯}{\overset{^}{z}}+B_y\underset{¯}{\overset{^}{y}}`$ is applied which also has a small, constant perpendicular field component $`B_y=0.05J_{\mathrm{FM}}`$ in order to define a certain path for the rotation of the magnetization during field reversal and to avoid the system to be trapped in a metastable state. The $`z`$ component of the field $`\underset{¯}{B}`$ is then reduced in steps of $`\mathrm{\Delta }B=0.004J_{\mathrm{FM}}`$ from $`B=0.2J_{\mathrm{FM}}`$ down to $`B`$ and afterwards raised again up to the initial value.

Typical hysteresis loops are depicted in Fig. 1. Shown are results for the magnetization of the FM (upper figure) as well as that of the AFM interface monolayer (lower figure) for different thicknesses $`t`$ of the AFM. The hysteresis loops of the FM clearly show EB depending on the thickness of the AFM. The magnetization curves of the interface layer of the AFM are shifted upwards due to the fact that after field cooling the AFM is in a domain state with a surplus magnetization. This layer experienced during cooling the external field in addition to the exchange field of the FM both having the same direction. This shifted interface magnetization of the AFM acts as an additional effective field on the FM, resulting in EB. The hysteresis curve of the AFM interface layer follows that of the FM layer with a much lower saturation magnetization, however. With increasing thickness of the AFM layer the area of the hysteresis loop of the interface layer which is proportional to the energy losses in the AFM decreases indicating that the spin structure in the AFM is stabilized.

The domain structure in the AFM interface layer is shown in Fig. 2 for an AFM consisting of one monolayer (upper figure) and for ten monolayers (lower figure), respectively. The fractal structure of these domains is obvious. It has been observed previously in bulk systems and was analyzed in detail . The domain structures shown in Fig. 2 are to a large extend frozen. But during field cycles small spin arrangements in the domain boundaries can take place even at low temperatures resulting in the hysteresis loops shown in Fig. 1.

The structure of the domains depends on the thickness of the AFM. For an AFM monolayer the effective field acting on all AFM spins is the superposition of the strong exchange field and the external field. But it is well known that the size of the domains depends on the strength of the effective magnetic field. Large fields imply small domains and vice versa . The small domains seen in the upper part of Fig. 2 are thus due to the strong fields acting on the AFM monolayer. On the other hand if the AFM consists of ten layers nine of them are only exposed to the weak external field with the tendency to form larger domain sizes. The coupling of these layers to the AFM interface layer then results in a coarsening of the domains at the interface as seen in the lower part of Fig. 2. Note that the distribution of vacancies in the interface layer is exactly the same in both parts of Fig. 2. A further obvious consequence of this explanation is that the domain size becomes layer dependent and increases with increasing distance from the AFM interface. But after a certain distance from the interface the domain structure should become independent of the interface layer which means that also the bias field should become independent of the thickness of the AFM for large $`t`$. This behavior indeed is observed in our simulations. In Fig. 3 we show the dependence of the bias field on $`t`$ for different dilutions of the AFM. The bias field is determined as $`B_{\mathrm{EB}}=(B^++B^{})/2`$ where $`B^+`$ and $`B^{}`$ are those fields of the hysteresis loop branches for increasing and decreasing field where the easy axis component of the magnetization of the FM becomes zero. The absolute value of the bias field increases rapidly with film thickness $`t`$, goes eventually over a maximum and then levels out. This is in agreement with experiments . Note that for the system with the smallest dilution the absolute value of the bias field decreases for $`t>1`$ much stronger with increasing thickness than for the more diluted films. The reason is that the less diluted systems have a stronger tendency to order antiferromagnetically thus reducing the net magnetization at the AFM interface. Again, agreement with experiments is obtained.

In conclusion, we have found by Monte Carlo simulations further support for our domain state model for exchange bias. So far all our numerical results have been obtained for systems where the AFM has a very strong anisotropy, i.e. behaves Ising-like. It is of great interest to relax this condition, i.e. to consider vector spin models for the AFM. Work in this respect is in progress.

Acknowledgments: This work has been supported by the Deutsche Forschungsgemeinschaft through SFB 491.

Figure 1: Nowak et al

Figure 2: Nowak et al

Figure 3: Nowak et al


text:
# Can we look inside a dynamo?

## 1 Introduction

At the very beginning, the question of the origin of cosmic magnetic fields was formulated as an inverse problem. ”How could a rotating body such as the sun become a magnet?” was the question Larmor put in 1919 (Larmor 1919), i.e. he asked for the unknown cause of an evident effect. In the last decades, dynamo theory saw an enormous progress. Testing a large number of kinematic models in the sense of a forward problem, the inverse character of the original problem was partly pushed into the background. It seems to be only a few percent of papers on the dynamo topic which are dealing with an inverse problem in the strict sense.

In geophysics, there is a long tradition to reconstruct the tangential components of the velocity field at the core-mantle boundary from geomagnetic secular variations (see, e.g, Bloxham (1989), and references therein). But the value of this type of inverse frozen-flux modeling was called into question recently (Love 1999). The main argument of the criticism was that the removal of the diffusion term from the induction equation changes its differential order resulting possibly in a drastic change of the nature of the solutions.

Another approach to the inverse dynamo problem was formulated by Love and Gubbins (1996). In this paper, the eigenvalue equation was complemented by some regularizing functionals in order to optimize dynamo models with respect to energetic demands and smoothness properties of the arising magnetic field. However, geomagnetic observations were taken into account only in the sense that some smoothness properties of the model magnetic field were compared with those of the observed field and not in the sense of fitting the model to the actual data.

With the advent of the first successful hydromagnetic dynamo experiments (Gailitis et al. 2000; Müller and Stieglitz 2000) quite new perspectives for inverse dynamo theory come into play. Compared to cosmic bodies, much more measurement techniques and measurable quantities can be taken into account. One additional measurable quantity is the electric potential at the boundary of the dynamo module which can provide complementary information in addition to the magnetic field on the outside. Further on, in contrast to the case of cosmic bodies it is possible to apply various magnetic fields from outside and to determine the response of the dynamo on these primary fields. Those date can be collected in the regime below the critical magnetic Reynolds number ($`Rm`$) in order to predict the dynamo behaviour at larger $`Rm`$.

In connection with this, but up to now restricted to small $`Rm`$, there are some new developments concerning the determination of the velocity field from externally measured induced magnetic fields and electric potentials when primary magnetic fields are applied (Stefani and Gerbeth 1999, 2000a, 2000b). For a static homogeneous primary field and spherical geometry it was shown (Stefani and Gerbeth 2000a) that the defining scalars of the velocity can be determined except for an ambiguity of the radial dependence of their spherical harmonics expansion coefficients. Although one can try to get a reasonable guess of the unknown radial dependence by some regularization techniques, in the strict sense this dependence remains unobservable. A generalization of this method to higher $`Rm`$ would amount to a combination of an eigenvalue solver (similar to the philosophy of Love and Gubbins) with a data-fitting procedure. However this approach might look like in detail, we guess that the radial dependence will remain undetectable. Roughly speaking, two two-dimensional measured quantities (radial component of magnetic field plus electric potential) allow only to determine two two-dimensional desired quantities. Actually, this radial ambiguity was the first reason for the following investigation. At least for the small $`Rm`$ regime it seems clear that the necessary third dimension to determine the radial dependence can only be the time (or frequency).

Whereas a general inverse dynamo theory is still missing, it is worth to have a sidelong glance at the quite similar inverse problem in quantum mechanics. There is a huge literature on inverse scattering theory and inverse spectral theory (for an overview, see Chadan and Sabatier (1977) and Chadan et al. (1997)), and evidently many methods from quantum mechanics (Gelfand-Levitan equation, WKB methods, operator methods, etc.) remain to be adapted and utilized in inverse dynamo theory.

In the present paper, we will pursue a less ambitious program restricting ourselves to a pragmatic treatment of a simple inverse spectral dynamo problem. It is a simple problem as it concerns spherically symmetric $`\alpha ^2`$-dynamos in which, as a slight modification of the original model of Krause and Steenbeck (1967), the radial dependence of $`\alpha `$ is assumed to be unknown. The treatment will be pragmatic in the sense that we suppose eigenvalues of only a few magnetic field modes to be ”measured” and that we will try to infer the radial dependence of $`\alpha `$ from this restricted set of information.

Basically, it is our aim to get some acquaintance with the peculiarities if one tries to get information on the radial distribution of the source of dynamo action. Reminding quantum theory might be useful in the following respect: The existence of isospectral potentials is a well established fact in quantum theory. For example, it is possible to deform continuously the harmonic oscillator potential in such a way that all the eigenvalues are the same as for the original quadratic potential (De Lange and Raab 1991). Therefore, we should be aware of the possibility that any inversion procedure of any set (even of an infinite number) of measured eigenvalues might lead to a variety of radial $`\alpha `$-distributions which cannot be discriminated using only their eigenvalues. For the corresponding quantum inverse problem it is known that a potential can be uniquely determined if the spectra for two different boundary conditions were known. A similar behaviour should also be expected for the inverse dynamo problem, but this seems to be irrelevant for any realistic measurement strategy.

Our method might have application for cosmic dynamos as well as for laboratory dynamos. Of course, a spherically symmetric $`\alpha ^2`$-dynamo model is surely too simple to explain cosmic dynamos but the general method is easily extendable to more realistic models.

The paper is structured as follows: in the second section, the used eigenvalue solver is presented and five paradigmatic types of radial dependence of $`\alpha `$ are treated as a forward problem. In the third section, the evolutionary strategy to invert the ”measured” growth rates for a number of magnetic field modes into a radial dependence of $`\alpha `$ is presented. Then, we try to reconstruct five models from section 3. Numerically, we find quite different $`\alpha `$-profiles giving nearly the same eigenvalues. Additionally, we treat two other problems which can be, at least in principle, of interest for future dynamo experiments. At first we try to find such a function $`\alpha (r)`$ which yields zero growth rate for the lowest four field modes. Such a profile for $`\alpha `$ might be a candidate for a (hypothetical) dynamo experiment showing such interesting effects like mode-switching. Further on, we use our code for an ”energetic” optimization of $`\alpha ^2`$-dynamos. The paper closes with some remarks on possible future developments.

## 2 The forward problem

Shortly speaking, our approach to solve the inverse problem will be the most simple one: to solve the forward problem many times and to seek, in some appropriate manner, for those dynamo models which give eigenvalues as close as possible to the measured ones. Therefore, the present chapter on the forward problem describes also the kernel of the inverse problem solving code which will be presented in detail in the next section.

We start with the induction equation for a mean-field dynamo model restricted to a spherically symmetric $`\alpha `$-coefficient. For the sake of simplicity, we ignore any large-scale velocity field as well as any $`\beta `$-effect or any anisotropic $`\alpha `$-coefficients. The electrical conductivity $`\sigma `$ is assumed to be constant inside a sphere of radius $`R`$ and is assumed to be zero in the outer part. Then, inside the sphere, the magnetic field has to satisfy the equations (Krause and Rädler 1980)

$`{\displaystyle \frac{B}{t}}=\times (\alpha B)+{\displaystyle \frac{1}{\mu _0\sigma }}\mathrm{\Delta }B,B=0.`$ (1)
At the boundary with $`r=R`$, the magnetic field has to match continuously to a potential field.

In the usual manner, $`𝐁`$ is represented as a sum of poloidal and toroidal components,

$`B=\times (r\times S)r\times T`$ (2)
with the defining scalars $`S`$ and $`T`$ expanded in spherical harmonics according to

$`S(r,\theta ,\varphi )`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}RS_l^m(r)Y_l^m(\theta ,\varphi )\mathrm{exp}\lambda _lt`$ (3)
$`T(r,\theta ,\varphi )`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}T_l^m(r)Y_l^m(\theta ,\varphi )\mathrm{exp}\lambda _lt.`$ (4)
In order to simplify the notation, throughout the rest of the paper we will measure the length in units of $`R`$, the time in units of $`\mu _0\sigma R^2`$, and the coefficient $`\alpha `$ in units of $`(\mu _0\sigma R)^1`$.

Using (2), (3) and (4), the induction equation can be transformed into the following form (Rädler 1986):

$`\lambda _lS_l`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{d^2}{dr^2}}(rS_l){\displaystyle \frac{l(l+1)}{r^2}}S_l+\alpha (r)T_l`$ (5)
$`\lambda _lT_l`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{d}{dr}}(rT_l)\alpha (r){\displaystyle \frac{d}{dr}}(rS_l)\right){\displaystyle \frac{l(l+1)}{r^2}}(T_l\alpha (r)S_l).`$ (6)
Note that there is no coupling between field modes differing in the order $`l`$ of the spherical harmonics. Furthermore, we have skipped the index $`m`$ for the coefficients of the defining scalars as it does not show up in the equations. Note in addition that the eigenvalues $`\lambda `$ may be complex.

At the surface $`r=1`$, the boundary conditions

$`{\displaystyle \frac{dS_l}{dr}}+(l+1)S_l=T_l=0`$ (7)
must be fulfilled.

In order to have a comparable global measure of the $`\alpha `$-effect we introduce the following definition:

$`C=3{\displaystyle _0^1}\alpha (r)r^2𝑑r.`$ (8)
Of course, other definitions, e.g. one without the $`r^2`$-term in the integral, are as well conceivable. The definition used here might be sensible for the purpose of ”energetic” optimization as it reflects the volume averaged intensity of $`\alpha `$.

To solve the coupled system of equations (5) and (6) together with the boundary condition (7) we have applied a standard shooting procedure using Newton’s method. This code was validated for the case of constant $`\alpha `$ as well as for some known results with a radial variation of $`\alpha `$ (Rädler 1986).

In the following, we will present the results of the forward problem for five paradigmatic models. Table 1 shows the corresponding functions explicitely. The somewhat strange numerical factors arise due to the definition (8) The examples a, b, c, d, and e represent the following types of $`\alpha `$-profiles: type a1 gives a constant $`\alpha `$ whereas $`a2`$ is only a slightly modified version with a small oscillation around the constant value. Curves of type b and c are monotonically increasing and decreasing, respectively. Curves of type d and e have one maximum and one minimum between 0 and 1, respectively.

In Fig. 1 the dependence of the growth rates $`\lambda `$ on $`C`$ for the modes with $`1l6`$ is shown. For the sake of simplicity, we take into account only the lowest radial wavenumber $`n=1`$. Admittedly, it is not completely consistent to include, e.g., the mode with $`l=3`$ and $`n=1`$ having a growth rate of $`\lambda (C=0)=33.1`$ and to skip at the same time the mode with $`l=1`$ and $`n=2`$ having a growth rate of $`\lambda (C=0)=20.2`$ (Krause and Rädler 1980). In principle, there is no obstacle to take into account modes with higher radial wavenumbers. The restriction to $`n=1`$ might suffice for the simple toy model of inverse dynamo theory as it is considered here.

Fig. 1 gives clear evidence of a dependence of the curves for different values of $`l`$ and their relations to each other on the function $`\alpha (r)`$. For example, a monotonic increase of $`\alpha `$ (as in b1, b2) and a maximum between $`0`$ and $`1`$ (d1, d2) seems to make the curves for different $`l`$ ”converging” to each other. On the other hand, monotonic decrease (c1, c2) or a minimum (e1, e2) seems to make the curves more parallel or even more moving away. However this might be (and we will see later that this behaviour is not universal) we only note here that there are clear differences in the relations of the eigenvalues of different field modes which will be used in an inversion scheme to determine $`\alpha (r)`$. Apart from example a2, we used only curves with at most one maximum or minimum. In the next section we will face the problem that some oscillations (with at least one maximum and one minimum) around the treated paradigmatic model give eigenvalues which are (at least) very close to the original one.

Note that the restriction to values $`C<6`$ was not voluntary as we have found (for some of the models) regions of $`C`$ were no $`\mathrm{𝑟𝑒𝑎𝑙}`$ eigenvalues exist. This gives evidence for the existence of oscillatory dynamos. This possibility has been discussed in a number of papers (Rädler 1986, Rädler and Bräuer 1987) and is of some interest in its own right. Only for the sake of numerical simplicity we skip this problem here. There are no principle obstacles to make the eigenvalue code fit for complex computations, too. As for the inverse problem, which will be treated in the next section, we only remark that in case that some intermediate trial function of $`\alpha (r)`$ are found to have no real eigenvalues, these configurations are skipped and are replaced by the next trial function.

## 3 The inverse problem

### 3.1 Evolutionary Strategies

In this section, we want to present and apply an inversion scheme using ”measured” growth rates for a finite number of magnetic field modes to get the radial dependence of $`\alpha `$. From the very beginning, we will take into account the possibility that there could exist various solutions of this kind of inverse problem, i.e. various $`\alpha `$-profiles which can be distinguished considering, e.g., their mean quadratic curvature.

For that reason, we decided to use a regularization of the inverse problem. One convenient way of regularizing an inverse problem is to add to the usual functional of the mean squared residual deviation an additional weighted functional of some appropriate norm of quantities one is looking for and to search for minima of the arising total functional. Using this so-called Tikhonov regularization approach (for an overview see, e.g., Hansen 1992), one can scale through the regularization parameter. This way one usually gets Tikhonov’s L-curve. At the point of highest bending (the ”knee”) one can find a reasonable compromise between an optimal fitting of the model data to the measured data and a minimal norm of the desired quantity.

Here, we will slightly modify this regularization method in the following way. To begin with, we restrict the space of functions to polynomials of fourth order, excluding the term proportional to $`r`$ ($`\alpha (r)r`$ would correspond to a cusp of $`\alpha `$ at the origin which seems to be not very realistic). These polynomials are parameterized by values $`\stackrel{~}{\alpha }_i`$ at the fitting points $`r_i=0.25,0.5,0.75,1.0`$ where the polynomial expansion of $`\alpha (r)`$ has to fit optimally the values $`\stackrel{~}{\alpha }_i`$. In this first fitting procedure, however, we use an additional weighted functional which keeps the mean squared curvature more or less small. For a high regularization parameter any set of (four) values at the fitting points leads to relatively smooth curves in the space of polynomial expansion coefficients. When the regularization parameter goes to zero the fitting curve tends exactly to the parameterizing values at the fitting points. Thus, we are able to regularize the curves even before they are put into the kernel of the inversion scheme.

To put this procedure in formulae, let us assume that we have a trial set of values $`\stackrel{~}{\alpha }_i`$ at the fitting points $`r_1=0.25`$, $`r_2=0.5`$, $`r_3=0.75`$, and $`r_4=1.0`$. Fitting a function $`\alpha (r)=a+br^2+cr^3+dr^4`$ will be done by minimizing the functional

$`F={\displaystyle \underset{i=1}{\overset{4}{}}}(\alpha (r_i)\stackrel{~}{\alpha }_i)^2+p{\displaystyle _0^1}(2b+6cr+12dr^2)^2𝑑r`$ (9)
where $`p`$ denotes the regularization parameter.

After having explained the chosen parameterization, we will describe in the following the chosen minimum search procedure. We decided to use an evolutionary strategy (ES) for that purpose (for an overview see, e.g., Schöneburg, Heinzmann and Feddersen 1994). An ES is known to be an appropriate method in case of existence of several local minima. It is capable to leave local minima in order to find the global minimum of a given problem. We will sketch the philosophy of the used ES only shortly. For our particular code, we took over the scheme described by Küchler, Lehnert and Tschornack (1995).

An ES tries to utilize the principle of biological evolution for the numerical solution of (non-linear) optimization problems. Let us start with an ”population” of 30 ”individuals” which are vectors of parameterizing values at the fitting points. To each of those individuals a quality function $`Q`$ is ascribed. In order to compute this quality function $`Q`$, we have to solve, for every considered $`l`$, the eigenvalue equation for the trial function $`\alpha (r)`$ which results from the vectors of parameterizing values the way described above. Having computed the eigenvalues for the given individual we compare these eigenvalues with the measured ones and define the quality function as the mean squared residual deviation. Here, we assume the a-priori errors for all $`l`$-mode measurements to be equal. Of course, this might be corrected in any practical application.

The population evolves in the following manner: In every new generation only one child is created by one of three different kinds of reproduction which are mutation, cross-breeding, and crossing-over. For every generation the actual reproduction kind is used by random. For the relative frequency of mutation, cross-breeding, and crossing-over we have chosen the values 0.3, 0.3, and 0.4, respectively.

Mutation means that we select one ”parent” by random and create from this a ”child” by mutation of the parameters, i.e.

$`\stackrel{~}{\alpha }_i^C=\stackrel{~}{\alpha }_i^P+\mathrm{\Delta }_i\gamma s.`$ (10)
Here $`\stackrel{~}{\alpha }_i^C`$ is the i’th parameter of the child, $`\stackrel{~}{\alpha }_i^P`$ is the i’th parameter of the parent, $`\mathrm{\Delta }_i`$ is the spread of the corresponding parameter, $`\gamma `$ is a mutation width and $`s`$ is a normal distributed random number. It should be noted that $`\mathrm{\Delta }_i`$ is large at the beginning of the evolution when a large area of the parameter space is to be covered, but that it decreases when the whole population is running into the global minimum.

Cross-breeding means that the parameters of two parents $`P1`$ and $`P2`$ which are selected by random are averaged in the child:

$`\stackrel{~}{\alpha }_i^C=0.5(\stackrel{~}{\alpha }_i^{P1}+\stackrel{~}{\alpha }_i^{P2}).`$ (11)
Two parents are also selected in the crossing-over, but their parameters will be mixed in the sense that for each parameter $`\stackrel{~}{\alpha }_i^C`$ it will be selected by random whether it comes from parent $`P1`$ and $`P2`$.

After having created a child, it must be evaluated. At first, it has to be checked whether the parameters of the child still fulfill some reasonable given constraints. If this is not the case, a new child has to be created. In the positive case, however, it is checked whether the quality function of the child fulfills the condition

$`Q_c<Q_{worst}+0.4(Q_{worst}Q_{best})`$ (12)
where $`Q_{worst}`$ and $`Q_{best}`$ are the quality functions of the worst and the best individual of the population, respectively. The value 0.4 is to some extend arbitrary. Evidently this factor is responsible for the fact that the ES is able to leave local minima. The chosen factor 0.4 is relatively high. The effect is a certain slow-down of the convergence but a higher reliability to find the global minimum and not to be locked in a local one. If the condition (13) is fulfilled, after cross-breeding and crossing-over the qualitatively worse parent is replaced by the child. In case of mutation, the parent is replaced, except when the parent is the best individual of the population and the quality of the child is not better than the quality of the parent. In this special case, the worst individual of the population is replaced by the child.

Evidently, the number of individuals in every generation stays the same. The evolution can be stopped when all individuals of the population have gathered in one (hopefully: the global) minimum, e.g., when the spread of all parameters has become very small. Of course, if several minima of the quality function with nearly the same depth exist, the evolution might be locked in quite different of them, depending on the random numbers used in the evolution. In order to control this erratic behaviour, we have used the regularization scheme described above. Thus, scaling the regularization parameter $`p`$ from larger to smaller values, we expect the smoothest solution to show up first and solutions with higher curvature to appear only for smaller values of $`p`$. Than, however, the mentioned erratic behaviour might be observed.

### 3.2 Results for the paradigmatic examples

Fig. 2 shows the results of the used ES inversion scheme for the examples a1, b1, c1, d1, and e1 from section 2. Every plot on the right hand side shows nine curves $`\alpha (r)`$, the first representing the exact input function which was used in the forward task for the determination of the growth rates for different $`l`$-modes. The remaining eight curves give the output of the ES inversion scheme, differing in the value of the used regularization parameter $`p`$ in equation (9). As input for the inverse task, we have used the growth rates for the first 6 $`l`$-modes at $`C=4`$. Reminding that we parameterize the desired functions $`\alpha (r)`$ by means of 4 values, 6 growth rates give evidently some redundant information.

The plots on the right hand side of Fig. 2 are similar to what is called in regularization theory ”Tichonov’s L-curve”. On the abscissa the logarithm of the squared curvature of the output curves is given. The ordinate axis shows the logarithm of the rms of the residual deviation (i.e., the difference of the ”measured” growth rates and the growth rates resulting from the output curves). For further explanation, let us consider example b1 having a quadratic $`\alpha (r)`$-dependence as original function (example a1 is atypical, as the original curve has curvature zero). To begin with, consider in b1 the curves for $`p=10`$ and $`p=1`$. Evidently, the resulting functions $`\alpha (r)`$ are close to a constant, e.g., the large regularization parameter $`p`$ allows only for output curves of very small mean curvature. On the right hand side we see that the corresponding residual deviation for those curves are very large. Then, with decreasing $`p`$, we get three curves which are fitting reasonably to the original one. The corresponding residual deviations are decreasing drastically. But for the last three values of $`p`$ ($`p=10^4`$, $`p=10^5`$, and $`p=0`$) the resulting functions $`\alpha (r)`$ acquire some oscillatory behaviour. Looking at the right hand side plot, we see that the residual deviations for these functions are a bit higher compared to the minimum, but evidently the ES is not able to find the global minimum. A quite similar behaviour is found for all other examples (note that for example a1 ”Tikhonov’s L-curve” starts already at high values of $`p`$ with a small residual error, as the original curve has no curvature).

One might argue that another minimum search code would be better suited to find the absolute minimum. Note, however, that every real measurement is biased by measurement errors. Considering the value of the rms of the residual deviations on the right hand side, it seems clear that a realistic discrimination between the curves is hardly possible. Note, in addition, that our special parameterization (with a power series expansion up to fourth order) might prevent curves for smaller regularization parameter from fitting to ”measured” growth rates more exactly.

Connected with this problem one could also argue that the incorporation of more measured growth rates could improve the quality of the inversion. For this purpose, we carried out a numerical test taking into account different total numbers (denoted by $`NGRO`$) of ”measured” growth rates. For example a1 we have tested all possibilities with $`4NGRO12`$ at $`p=0`$. Fig. 3a shows the resulting functions. Figs. 3b, 3c, and 3d give the corresponding growth rates for the exact function and all functions resulting from the inversion (Fig. 3b shows the values for all modes with $`1l12`$, Fig. 3c shows the details for the first modes with $`1l3`$, Fig. 3d shows the details for the modes with $`10l12`$). The result is astonishing: all the quite different functions $`\alpha (r)`$ shown in Fig. 3a give nearly the same growth rates for the first modes with $`1l12`$, apart from very small differences. Of course, our numerical method and the used parameterization with a very limited number of parameters does not allow to construct exactly those different functions $`\alpha (r)`$ which might give exactly the same eigenvalues. At this point we would like to suggest analytical investigations in order to identify these isospectral functions $`\alpha (r)`$.

For any practical application where the accuracy of the measured values will be limited the consequences are clear: there are quite different functions $`\alpha (r)`$ giving nearly the same growth rates. Thus, a unique determination of the radial dependence of $`\alpha `$ from the eigenvalues of different $`l`$-modes seems to be impossible. However, using some kind of regularization, we can find those $`\alpha `$-functions giving the measured eigenvalues with minimal penalty function, e.g, with minimal quadratic curvature. Hence, it is possible to reconstruct at least $`\alpha `$-functions of the paradigmatic types as described in section 2.

### 3.3 Search for $`\alpha `$-profiles yielding level-crossing

In the following the described ES will be used to treat another problem. Whereas in the last subsection we had used growth rates resulting from known functions $`\alpha (r)`$ we will now search for a-priori unknown functions. Let us assume that we would know that a cosmic body has a steady magnetic field with some field modes, say, for $`1l4`$. If we would know (what is, of course, not very realistic) that the dynamo mechanism were only due to a spherically symmetric $`\alpha `$ we could ask for its radial dependence providing exactly these four marginal modes. This is a typical example our inversion scheme can be applied for.

Table 2 gives the functions resulting from the inversion, again for a number of decreasing regularization parameters $`p`$, together with the corresponding values for the first four growth rates. We actually give here the resulting power expansions in order to allow the reader to validate the results by means of a standard forward solver for spherically symmetric $`\alpha (r)`$.

These functions are depicted also in Fig. 4a. The dependence of the residual deviation on the curvature is again shown in Fig. 4b. Whereas for large $`p`$ the growth rates differ significantly from zero, the values for smaller $`p`$ are getting very close to zero. For the best fitting $`\alpha (r)`$ (at $`p=10^4`$) we have shown in Fig. 4c the dependence of the growth rates on $`C`$ for $`1l6`$. We see clearly the level-crossing at $`C8.5`$.

There might also be some relevance of the considered example for laboratory dynamos. Let us assume that we were able to build a spherical dynamo and that the dynamo action of its flow were due to a spherically symmetric $`\alpha `$. Then, a radial dependence of $`\alpha `$ corresponding to the given curves might lead to a very interesting behaviour when the back-reaction is taken into account. For example it might happen that the back-reaction changes the function $`\alpha (r)`$ in such a way that mode switching between the different $`l`$-modes occurs. We would like to suggest further investigation of such a model.

### 3.4 Optimization

Let us turn to a second problem which might also be important for laboratory dynamos. It concerns a kind of energetic optimization of the dynamo. Of course, the optimization criterion for a mean field dynamo is not as obvious as for large scale velocities. We have chosen to minimize the following functional:

$`C^{}=\left(3{\displaystyle _0^1}\alpha ^2(r)r^2𝑑r\right)^{1/2}`$ (13)
under the constraint that the mode with $`l=1`$ becomes marginal. Of course, the definition (13) is to some extend arbitrary. For comparison it is also interesting to consider the quantity

$`C^{\prime \prime }=3{\displaystyle _0^1}|\alpha (r)|r^2𝑑r`$ (14)
where we used the modulus of $`\alpha `$ instead of $`\alpha `$ itself in order not to get very small values of $`C^{\prime \prime }`$ only due to some cancellation effects along the $`r`$-axis.

To solve this problem, the described ES code had to be modified only slightly. The resulting functions, together with the corresponding values $`C^{}`$ and $`C^{\prime \prime }`$ are given in Table 3.

We do not claim to have found the function $`\alpha (r)`$ with the smallest value of $`C^{}`$ or $`C^{\prime \prime }`$. Remind that we have restricted the allowed functions to be a polynomial expansion up to fourth order. There might be quite different functions providing even lower values of $`C^{}`$ and/or $`C^{\prime \prime }`$.

Despite some arbitrariness of the definitions of $`C^{}`$ and $`C^{\prime \prime }`$, it is interesting to note that something like the volume averaged ”amount” of $`\alpha `$ can be decreased significantly compared to the well known critical value 4.49 for the Krause-Steenbeck dynamo model.

## 4 Conclusions

In this paper we have tried to reconstruct the radial distribution of $`\alpha `$ from growth rates of a few $`l`$-modes of the magnetic field.

It was shown that some paradigmatic functions $`\alpha (r)`$ can indeed be recovered by the used evolutionary strategy if a careful regularization of the trial functions is used. With decreasing smoothing effect of the regularization, the ES erratically jumps into various local minima and the resulting functions $`\alpha (r)`$ can differ drastically from the original one.

Although by our simple numerical means we are not able to test the isospectrality of different $`\alpha `$-profiles in a strict sense, the obvious existence of different minima with nearly the same quality functions seems to be more than an accident. A more precise identification of such isospectral $`\alpha `$-profiles is highly desirable and we would like to suggest further analytical and numerical investigation in this direction.

On the other side, for any practical application (e.g. with measurement errors and only a finite number of $`l`$-modes) it is quite clear that the function $`\alpha (r)`$ can be determined only roughly.

A unique determination of $`\alpha (r)`$ would probably require knowledge of the spectra for two different boundary conditions which is, even for laboratory dynamos, not very realistic.

An obvious future task is to generalize the presented method to the oscillatory case, possibly in connection with varying conductivity (for possible astrophysical applications of this case, see Schubert and Zhang 2000). In principle, the developed ES code can be combined with any forward solver for more complicated dynamo models in order to apply the method to real cosmic or laboratory dynamos. For other than spherical geometry it might also be useful to combine the method with an forward solver based on the integral equation approach (Stefani, Gerbeth, and Rädler 2000).