text:
# Pattern and wavenumber selection in ferrofluids

## I Introduction

When a ferrofluid layer is subjected to a vertically oriented and uniform magnetic field, above a critical value of the field strength a hexagonal pattern of peaks appears on the surface of the liquid. This Rosensweig or normal field instability was first observed by Cowley and Rosensweig in 1967 . Further increase of the magnetic field up to a second threshold gives rise to a transition from the hexagonal to a square planform .

The arrangement of peaks resulting from the Rosensweig instability into patterns of different geometry is just one particular example from an impressive variety of pattern formation in physical systems . It is well known that although the instability threshold itself can be obtained from a linearized version of the underlying equations, the pattern selection problem requires the inclusion of non-linear terms. A standard procedure to probe the non-linear regime perturbatively is by means of amplitude equations .

However, unlike many other examples of pattern forming physical systems discussed in the literature the surface profile of a ferrofluid in a static magnetic field is an equilibrium structure. Accordingly the relative stability of planforms and possible transitions between different patterns can be investigated theoretically by studying the appropriate thermodynamical potential. Still the problem is a complicated non-linear task since the local magnetic field determining the surface profile in turn depends on the surface deflection via boundary conditions. As a consequence the variational minimization of the thermodynamic potential in the surface profile cannot be accomplished exactly and one has to resort to approximate procedures applicable in the experimentally relevant situation of slightly overcritical magnetic fields.

In classical investigations along these lines Gailitis  and Kuznetsov and Spektor  analyzed the stability of the different patterns of the Rosensweig instability by means of an energy minimization principle. These as well as related investigations using methods of functional analysis  or a generalized Swift-Hohenberg equation  were confined to fluid layers of infinite depth.

On the other hand experimental investigations are usually done with rather thin fluid layers and the effects of finite thickness on the linear regime have been studied recently in great detail. For example the dispersion relation of surface waves on a ferrofluid layer of arbitrary thickness was determined in  and the influence of viscosity on the linear dynamics was elucidated in . Taking advantage of the dependence of the threshold of instability and the wavelength of the most unstable mode on the thickness of the fluid layer in a clever way, it is e.g. possible to measure both the normal and the anomalous dispersion branch of surface waves on ferrofluids .

In the present paper we complement these investigations by a thorough study of the weakly non-linear regime of the Rosensweig instability slightly above the critical magnetic field for a ferrofluid of arbitrary depth. Taking the limit of infinite layer thickness we also critically discuss the classical findings obtained in . Moreover we are able to quantify the restriction to sufficiently small susceptibilities $`\chi `$ of the fluids which was always used in previous studies.

Our method of investigation is a generalization of the variational minimization of an energy functional already used in . Near the instability this functional may be written as a power series in the amplitude of the surface deflection and the minimization can be performed explicitly. Moreover our approach allows the theoretical investigation of the wavenumber selection problem addressed also in recent experiments . Including the wavenumber $`k`$ into the set of variational parameters we determine the dependence of the wavenumber of the patterns on the magnetic field and investigate the influence of a varying wavenumber on the stability of hexagons and squares. Therefore by a slight extension of our method we are able to study questions which are not easily accessible to other non-linear methods, as e.g. amplitude equations.

The paper is organized as follows. In section II the basic equations are collected and transformed into a form suitable for the calculation of the energy. In section III a perturbation ansatz for the surface deflection is put forward and the important issue of its consistency and region of validity is discussed. Section IV is devoted to the comparison of our findings with the classical results of pattern selection in a fluid with infinite depth. Subsequently in section V we consider a ferrofluid layer with arbitrary thickness and analyze the effects of the finite depth on the stability of the patterns. In section VI we address the wavenumber selection problem. Finally section VII contains a summary and compares our findings with recent experimental results.

## II Basic equations

We consider the situation sketched in Fig. 1. A horizontally unbounded ferrofluid layer is subjected to an external magnetic field $`𝐇_0`$ which in the absence of any magnetic permeable material is of the form $`𝐇_0(x,y,z)=H_0𝐞_z`$. The incompressible magnetic fluid of density $`\rho `$, surface tension $`\sigma `$ and susceptibility $`\chi `$ is bounded from below at $`z=d`$ by an impermeable material and has a free surface described by $`z=\zeta (x,y)`$ with the magnetically impermeable air above. The gravitational acceleration $`𝐠=g𝐞_z`$ acts parallel to the z-axis. Our aim is to determine which static profile $`\zeta (x,y)`$ develops for a magnetic field $`𝐇_0`$ strong enough to destabilize the flat surface $`\zeta (x,y)=0`$.

Stable configurations of the ferrofluid surface with infinite horizontal extension are given by minima of the thermodynamic potential per unit area in the $`xy`$ plane

$$f[\zeta (x,y)]=\frac{\rho g}{2}\zeta ^2(x,y)\frac{\mu _0}{2}\underset{d}{\overset{\zeta (x,y)}{}}𝑑z𝐇_0\chi 𝐇(x,y,z)+\sigma \sqrt{1+(_x\zeta (x,y))^2+(_y\zeta (x,y))^2}.$$
(1)
Here $`\mu _0`$ is the permeability of free space, $`𝐇(x,y,z)`$ is the magnetic field in the presence of the ferrofluid, and the brackets denote the average over the $`xy`$ plane defined by

$$F(x,y):=\underset{L\mathrm{}}{lim}\frac{1}{4L^2}\underset{L}{\overset{L}{}}𝑑x\underset{L}{\overset{L}{}}𝑑yF(x,y).$$
(2)
The three terms on the r.h.s. of Eq. (1) describe the hydrostatic energy, the magnetic energy , and the surface energy respectively. Note that both the hydrostatic and the surface energy increase when the interface profile starts to deviate from the flat reference state, whereas the magnetic energy decreases. For sufficiently large $`𝐇_0`$ this gives rise to the normal field or Rosensweig instability .

The magnetic fields $`𝐁`$ and $`𝐇`$ are determined by the static Maxwell equations,

$$𝐁=0\mathrm{and}\times 𝐇=\mathrm{𝟎},$$
(3)
together with the boundary conditions,

$$\underset{z\pm \mathrm{}}{lim}𝐇(x,y,z)=H_0𝐞_z,$$
(4)
$$\left[(𝐁_\mathrm{a}𝐁_{\mathrm{ff}})𝐧\right]|_{z=\zeta }=0\mathrm{and}\left[(𝐇_\mathrm{a}𝐇_{\mathrm{ff}})\times 𝐧\right]|_{z=\zeta }=\mathrm{𝟎},$$
(5)
$$\left[(𝐁_{\mathrm{ff}}𝐁_\mathrm{b})𝐞_z\right]|_{z=d}=0\mathrm{and}\left[(𝐇_{\mathrm{ff}}𝐇_\mathrm{b})\times 𝐞_z\right]|_{z=d}=\mathrm{𝟎},$$
(6)
where $`𝐧`$ denotes the normal vector on the surface $`\zeta (x,y)`$ and the respective subscripts indicate the fields above, in and below the ferrofluid. The minimization problem for the thermodynamic potential becomes non-trivial since the magnetic field $`𝐇(x,y,z)`$ depends on the profile $`\zeta (x,y)`$ of the surface.

Throughout the paper we will assume the linear relation

$$𝐁=\mu _0(1+\chi )𝐇$$
(7)
between the magnetic induction and the magnetic field. For the experimentally relevant magnetic fluids and fields  this approximation deviates about 5% from the exact relation between $`𝐁`$ and $`𝐇`$.

It is convenient to introduce a scalar magnetic potential $`\psi (x,y,z)`$ defined by

$$𝐇=\psi $$
(8)
which by (3) has to satisfy the Laplace equation

$$\mathrm{\Delta }\psi =0.$$
(9)
The two characteristic scales of the problem are the critical wavenumber at the onset of the instability,

$$k_{\mathrm{c},\mathrm{}}=\sqrt{\frac{\rho g}{\sigma }},$$
(10)
for a ferrofluid with infinite depth $`d\mathrm{}`$ and the corresponding critical magnetic field,

$$H_{\mathrm{c},\mathrm{}}=\sqrt{\frac{(1+\chi )(2+\chi )2\sqrt{\rho g\sigma }}{\chi ^2\mu _0}},$$
(11)
which were first derived in  from a linear stability analysis. Henceforth we therefore measure all lengths in units of the capillary length $`k_{\mathrm{c},\mathrm{}}^1`$, all wavenumbers in units of the critical wavenumber $`k_{\mathrm{c},\mathrm{}}`$, the magnetic field $`𝐇_0`$ in units of the critical value $`H_{\mathrm{c},\mathrm{}}`$, the scalar magnetic potential $`\psi `$ in units of $`H_{\mathrm{c},\mathrm{}}/k_{\mathrm{c},\mathrm{}}`$, and energies per unit area in units of $`\sigma `$. Moreover it is convenient to introduce the following rescaled magnetic potentials in the space above, in and below the ferrofluid respectively:

$`\psi _\mathrm{a}`$ $`:=`$ $`\psi {\displaystyle \frac{(2+\chi )}{H_0\chi }},`$ (12)
$`\psi _{\mathrm{ff}}`$ $`:=`$ $`\psi {\displaystyle \frac{(1+\chi )(2+\chi )}{H_0\chi }},`$ (13)
$`\psi _\mathrm{b}`$ $`:=`$ $`\psi {\displaystyle \frac{(2+\chi )}{H_0\chi }}.`$ (14)
Each of these potentials has to fulfill the Laplace equation (9). Using the abbreviation

$$\eta :=\frac{\chi }{2+\chi },$$
(15)
the asymptotic form of the magnetic field for $`|z|\mathrm{}`$ as specified by Eq. (4) gives rise to the requirements

$$\underset{z+\mathrm{}}{lim}_z\psi _\mathrm{a}(x,y,z)=\frac{1}{\eta }=\underset{z\mathrm{}}{lim}_z\psi _\mathrm{b}(x,y,z)$$
(16)
for the magnetic potentials. Moreover the boundary conditions (5) and (6) for the magnetic fields translate into the following boundary conditions for the potentials:

$$\left[_x\psi _\mathrm{a}_x\psi _{\mathrm{ff}}\right]|_{z=\zeta }_x\zeta +\left[_y\psi _\mathrm{a}_y\psi _{\mathrm{ff}}\right]|_{z=\zeta }_y\zeta \left[_z\psi _\mathrm{a}_z\psi _{\mathrm{ff}}\right]|_{z=\zeta }=0,$$
(17)
$$\frac{1+\eta }{1\eta }\psi _\mathrm{a}|_{z=\zeta }\psi _{\mathrm{ff}}|_{z=\zeta }=0,$$
(18)
$$\left[_z\psi _{\mathrm{ff}}_z\psi _\mathrm{b}\right]|_{z=d}=0,$$
(19)
$$\psi _{\mathrm{ff}}|_{z=d}\frac{1+\eta }{1\eta }\psi _\mathrm{b}|_{z=d}=0.$$
(20)
Using Eqs. (8) and (13) and exploiting the fact that $`𝐇_0`$ is parallel to the z-axis, we finally get the energy $`f`$ as functional of the surface deflection $`\zeta =\zeta (x,y)`$ in the form

$$f[\zeta ]=\frac{\zeta ^2}{2}H_0^2\left(\psi _{\mathrm{ff}}|_{z=\zeta }\psi _{\mathrm{ff}}|_{z=d}\right)+\sqrt{1+(_x\zeta )^2+(_y\zeta )^2}.$$
(21)
As stated above, the main problem, rendering a straightforward minimization of $`f`$ in $`\zeta (x,y)`$ impossible, is the rather implicit dependence of the potential $`\psi _{\mathrm{ff}}(x,y,z)`$ on the surface deflection $`\zeta (x,y)`$ specified by the boundary conditions (17) and (18).

## III The perturbation ansatz

The variational problem for the energy functional posed in the last paragraph can in general not be solved exactly. In order to make analytic progress, we restrict ourselves to the vicinity of the critical magnetic field and assume that the amplitude of the surface deflection is still small. It is then possible to expand the energy in this amplitude and to retain only the first terms. In this paper we will consider an expansion of the energy up to fourth order in the amplitude. The applicability of this approach will be critically discussed below.

Generalizing the perturbation expansions used in  for the one-dimensional variant of the Rosensweig instability, we write the surface profile in the form

$$\zeta (x,y)=\underset{n=1}{\overset{4}{}}A_n\mathrm{cos}(𝐤_n𝐫)+\underset{n=5}{\overset{17}{}}B_n\mathrm{cos}(𝐤_n𝐫),$$
(22)
where $`𝐫=(x,y)`$ and $`𝐤=(k_x,k_y)`$ are two-dimensional vectors. The terms with $`n=1,\mathrm{},4`$ constitute the main modes. Their corresponding wavevectors $`𝐤_1`$ to $`𝐤_4`$ have all the same modulus $`k`$ whereas their mutual orientation is chosen such that the ansatz allows the description of ridges (e.g. $`A_1=A,A_2=A_3=A_4=0`$), squares (e.g. $`A_1=A_4=A,A_2=A_3=0`$), and hexagons (e.g. $`A_1=A_2=A_3=A,A_4=0`$) (see Fig. 2).

The terms with $`n=5,\mathrm{},17`$ are higher harmonics with wave vectors $`𝐤_5`$ to $`𝐤_{17}`$ being linear combinations of two wave vectors of the main modes and amplitudes $`B_n`$ of order $`O(A^2)`$. These terms are needed to satisfy the minimum conditions for the energy functional to the required order $`O(A^3)`$. The intuitive meaning of this fact is that reliable results on the relative stability of different planforms requires some information on the deviation of the nonlinear surface profile from the simple cosine shape describing the linear instability. On the other hand there is no need to include even higher harmonics in the ansatz (22) since the corresponding contributions would average to zero in the $`xy`$integrations in the definition (21) of the energy functional. In conclusion the chosen ansatz is the only consistent perturbation ansatz for the energy which includes terms up to fourth order in the amplitude of the surface deflection. The values of the amplitudes $`A_1,\mathrm{},A_4`$ and $`B_5,\mathrm{},B_{17}`$ as well as the wave vector modulus $`k`$ of the main modes are the free parameters which may be used to minimize the energy. In particular the possible minimization in $`k`$ is a special advantage of our approach since it allows an analysis of the non-linear wave vector selection problem not easily accessible to other approaches.

Because of the boundary conditions at the surface of the fluid the magnetic potential $`\psi `$ must have a similar dependence on $`x`$ and $`y`$ as $`\zeta (x,y)`$. The ansatzes

$`\psi _\mathrm{a}`$ $`=`$ $`{\displaystyle \frac{z}{\eta }}+{\displaystyle \underset{n=0}{\overset{17}{}}}u_n\mathrm{e}^{|𝐤_n|z}\mathrm{cos}(𝐤_n𝐫),`$ (23)
$`\psi _{\mathrm{ff}}`$ $`=`$ $`{\displaystyle \frac{z}{\eta }}+{\displaystyle \underset{n=1}{\overset{17}{}}}[v_n^+\mathrm{e}^{|𝐤_n|z}+v_n^{}\mathrm{e}^{|𝐤_n|z}]\mathrm{cos}(𝐤_n𝐫),`$ (24)
$`\psi _\mathrm{b}`$ $`=`$ $`{\displaystyle \frac{z}{\eta }}+{\displaystyle \frac{2d}{1+\eta }}+{\displaystyle \underset{n=0}{\overset{17}{}}}w_n\mathrm{e}^{|𝐤_n|z}\mathrm{cos}(𝐤_n𝐫),`$ (25)
fulfill the Laplace equation (9) and the asymptotic boundary conditions (16). Omitting the mode $`𝐤_0:=(0,0)`$ in the potential $`\psi _{\mathrm{ff}}`$ fixes the free constant of the potential $`\psi `$. Plugging Eq. (25) into the four remaining boundary conditions Eqs. (17)-(20) we can determine the parameters $`u_n,v_n^+,v_n^{}`$ and $`w_n`$ as functions of the surface deflection amplitudes by expanding the equations up to third order in the $`A_n`$. The resulting expressions specify the dependence of the magnetic field on the surface profile which is then used to calculate the energy $`f`$ using Eq. (21) up to fourth order in the amplitudes $`A_n`$. The resulting dependence of $`f`$ on the higher order amplitudes $`B_n`$ is simple and the minimization in the $`B_n`$ can be performed explicitly. Subtracting from the resulting expression the reference value of the energy for a flat interface we finally arrive at

$`f_d(\{A_n\},k)`$ $`:=`$ $`f(\{A_n\},k)f(\zeta =0)`$ (26)
$`=`$ $`{\displaystyle \frac{1}{2}}l(ϵ,k)[A_1^2+A_2^2+A_3^2+A_4^2]\gamma (ϵ,k)[A_1A_2A_3]+{\displaystyle \frac{1}{4}}g(ϵ,k)[A_1^4+A_2^4+A_3^4+A_4^4]+`$ (28)
$`+{\displaystyle \frac{1}{2}}g_h(ϵ,k)[A_1^2A_2^2+A_2^2A_3^2+A_3^2A_1^2]+{\displaystyle \frac{1}{2}}g_t(ϵ,k)[A_2^2A_4^2+A_3^2A_4^2]+{\displaystyle \frac{1}{2}}g_n(ϵ,k)[A_1^2A_4^2]+O(A^5).`$
This expression gives the energy of a ferrofluid layer of depth $`d`$ with surface profile $`\zeta (x,y)`$ as specified by (22) at arbitrary strength of the external field $`H_0`$ up to fourth order in the deflection amplitudes $`A_n`$. The dependence on $`H_0`$ is expressed via the supercriticality parameter

$$ϵ=\frac{H_0^2}{H_{\mathrm{c},d}^2}1.$$
(29)
The function $`l(ϵ,k)`$ is given by

$$l(ϵ,k)=\frac{1}{2}+(1+ϵ)H_{\mathrm{c},d}^2\frac{1e^{2kd}\eta }{1e^{2kd}\eta ^2}k\frac{1}{2}k^2.$$
(30)
It fixes the critical wavenumber $`k_{\mathrm{c},d}`$ for the onset of the instability and the corresponding critical magnetic field $`H_{\mathrm{c},d}`$ for a ferrofluid with finite depth $`d`$ via the conditions

$$l(ϵ=0)|_{k=k_{\mathrm{c},d}}=0$$
(32)
and

$$\frac{l(ϵ=0)}{k}|_{k=k_{\mathrm{c},d}}=0.$$
(33)
Determining the minima of $`f`$ with respect to the amplitudes $`A_n`$, one finds that $`ϵ0`$ with $`A_n0`$ as expected: the external field will be near to its critical value if the equilibrium values of the amplitudes $`A_n`$ are small. For ridges and squares the explicit scaling is $`ϵA^2`$ for $`A0`$, in the case of hexagons one finds $`ϵA`$ for small $`A`$. This in turn implies that to the desired order $`O(A^4)`$ of the expansion of $`f`$ in the amplitudes $`A_n`$ only the $`ϵ`$-dependence of $`l`$ and $`\gamma `$ has to be retained. The exact form of these dependencies is $`l(ϵ,k)=l(0,k)+ϵl(ϵ,k)/ϵ`$ and $`\gamma (ϵ,k)=(1+ϵ)\gamma (k)`$. In this way we arrive at our main result for the energy $`f_d(\{A_n\},k)`$ of the surface deflection which describes the equilibrium configurations consistently up to fourth order in the amplitudes $`A_n`$:

$`f_d(\{A_n\},k)`$ $`=`$ $`{\displaystyle \frac{1}{2}}l(ϵ,k)[A_1^2+A_2^2+A_3^2+A_4^2](1+ϵ)\gamma (k)[A_1A_2A_3]+{\displaystyle \frac{1}{4}}g(k)[A_1^4+A_2^4+A_3^4+A_4^4]+`$ (35)
$`+{\displaystyle \frac{1}{2}}g_h(k)[A_1^2A_2^2+A_2^2A_3^2+A_3^2A_1^2]+{\displaystyle \frac{1}{2}}g_t(k)[A_2^2A_4^2+A_3^2A_4^2]+{\displaystyle \frac{1}{2}}g_n(k)[A_1^2A_4^2]+O(A^5).`$
The explicit expressions for the various coefficients in (35) are rather long. We therefore display their form here only for the comparatively simple situation of a fluid of infinite depth, $`d\mathrm{}`$:

$$\gamma _{\mathrm{}}(k)=\frac{3}{4}\eta k^2,$$
(37)
$$g_{\mathrm{}}(k)=\frac{1}{16}k^3(83k)\frac{k^4\eta ^2}{14k+4k^2},$$
(38)
$$g_{h,\mathrm{}}(k)=\frac{1}{4}k^3(11\eta ^27\eta ^2\sqrt{3}+3\sqrt{3}\frac{3}{4}k3)\frac{3}{8}\frac{k^4\eta ^2(198\sqrt{3})}{12\sqrt{3}k+3k^2},$$
(39)
$`g_{t,\mathrm{}}(k)`$ $`=`$ $`k^3({\displaystyle \frac{3}{4}}\sqrt{3}\sqrt{2}1+4\eta ^2{\displaystyle \frac{7}{4}}\eta ^2\sqrt{3}\sqrt{2}{\displaystyle \frac{5}{16}}k)+`$ (41)
$`{\displaystyle \frac{1}{8}}{\displaystyle \frac{\eta ^2k^4(18+\sqrt{3}4\sqrt{2}4\sqrt{3}\sqrt{2})}{2k\sqrt{3}\sqrt{2}+k^2+k\sqrt{2}\sqrt{3}}}{\displaystyle \frac{1}{8}}{\displaystyle \frac{\eta ^2k^4(18\sqrt{3}+4\sqrt{2}4\sqrt{3}\sqrt{2})}{2k\sqrt{3}\sqrt{2}+k^2k\sqrt{2}+\sqrt{3}}},`$
$$g_{n,\mathrm{}}(k)=k^3(\frac{1}{8}k3\eta ^2\sqrt{2}1+4\eta ^2+\sqrt{2})\frac{k^4\eta ^2(94\sqrt{2})}{12k\sqrt{2}+2k^2}.$$
(42)
## IV Ferrofluid layer with infinite depth

In this section we consider a ferrofluid of infinite depth $`d\mathrm{}`$ and fix the wavenumber modulus $`k`$ at its critical value $`k_{\mathrm{c},\mathrm{}}=1`$. This enables us to compare our findings with the well-known results of previous theoretical investigations.

Possible surface patterns together with their amplitudes are given by stationary points of the energy function $`f_{\mathrm{}}(\{A_n\},k=1)`$. We will in particular study ridges (“rolls”) given by $`A_1=A_R`$ with $`A_2=A_3=A_4=0`$, squares represented by $`A_1=A_4=A_S`$ with $`A_2=A_3=0`$, and hexagons described by $`A_1=A_2=A_3=A_H`$ with $`A_4=0`$. The last planform is shown in Fig. 1. Using Eqs. (30), (35) and (III) we find for the corresponding amplitudes

$$A_R(ϵ)=\sqrt{\frac{ϵ}{g_{\mathrm{}}}},$$
(43)
$$A_S(ϵ)=\sqrt{\frac{ϵ}{g_{n,\mathrm{}}+g_{\mathrm{}}}},$$
(44)
and

$$A_H(ϵ)=\frac{1}{2}\frac{\gamma _{\mathrm{}}(1+ϵ)+\sqrt{\gamma _{\mathrm{}}^2(1+ϵ)^2+4ϵ(2g_{h,\mathrm{}}+g_{\mathrm{}})}}{2g_{h,\mathrm{}}+g_{\mathrm{}}}$$
(45)
respectively. Together with $`\chi `$ also $`\gamma `$ is always positive<sup>*</sup><sup>*</sup>*This holds true also in the case of a finite layer depth $`d`$.. Hence only up-hexagons with a peak instead of a dip in the center are stable.

Our results for the hexagons differ slightly from the classical ones reported in  in which the dependence of the cubic term in Eq. (35) on $`ϵ`$ is missing. In  this dependence was omitted because the discussion was restricted to fluids with $`\gamma 1`$ so that the approximate scaling $`ϵA^2`$ is very similar to the true one $`ϵ\gamma A+O(A^2)`$. Otherwise our results coincide with those of Gailitis, in particular using $`k=k_{\mathrm{c},\mathrm{}}=1`$ in (III) the respective expressions found in  are reproduced. The investigations of Kuznetzov and Spektor  were restricted to first order in $`\eta `$ and are contained in both our results and those of Gailitis.

Before addressing the relative stability of the different planforms introduced above it is important to characterize the domain of validity of our central expression (35) for the energy. Being a perturbative result the corresponding amplitudes $`A_n`$ must not be too large. Fig. 3 gives a quantitative impression of what that means by displaying lines in the $`\chi ϵ`$ plane corresponding to fixed values of the amplitudes of different patterns. It is in particular important to realize that for sufficiently large values of the susceptibility $`\chi `$ the fourth order coefficients in the energy functional may change sign. The corresponding patterns than appear through a backward bifurcation . At the same time higher order terms (i.e. sixth or higher order in the amplitudes $`A_n`$) in the energy are necessary to saturate the instability. The concrete values of $`\chi `$ for the planforms considered are $`\chi _R2.54`$ ($`g_{\mathrm{}}`$ becomes negative), $`\chi _S0.41`$ ($`g_{n,\mathrm{}}+g_{\mathrm{}}`$ becomes negative), and $`\chi _H1.05`$ ($`2g_{h,\mathrm{}}+g_{\mathrm{}}`$ becomes negative). Beyond these values of $`\chi `$ our ansatz (22) is unable to describe the arising pattern. The value of $`\chi _R`$ was first calculated in .

Keeping this limitation of our approach in mind we have investigated the pattern selection problem by studying the character of the extremum of the energy functional (35) at the ridge, square and hexagon solutions given in Eqs. (43)-(45). This is rather similar to a linear stability analysis of the fixpoint solutions of the corresponding amplitude equations . The results are as follows: Ridges are never stable, since $`g_{n,\mathrm{}}<g_{\mathrm{}}`$ always. We have always $`g_{\mathrm{}}+g_{n,\mathrm{}}<g_{h,\mathrm{}}+g_{t,\mathrm{}}`$ and $`g_{n,\mathrm{}}<g_{\mathrm{}}`$. Hence (for $`\chi <\chi _S`$) squares are stable if $`ϵ>ϵ_S={\displaystyle \frac{(g_{\mathrm{}}+g_{n,\mathrm{}}g_{h,\mathrm{}}g_{t,\mathrm{}}+\sqrt{(g_{\mathrm{}}+g_{n,\mathrm{}}g_{h,\mathrm{}}g_{t,\mathrm{}})^24\gamma _{\mathrm{}}^2(g_{\mathrm{}}+g_{n,\mathrm{}})})^2}{4\gamma _{\mathrm{}}^2(g_{\mathrm{}}+g_{n,\mathrm{}})}}`$. In order for hexagons to be stable (for $`\chi <\chi _H`$) the following three conditions have to be fulfilled.

* $`ϵ>ϵ_F={\displaystyle \frac{(\gamma _{\mathrm{}}^2+2g_{h,\mathrm{}}+g_{\mathrm{}}\sqrt{(\gamma _{\mathrm{}}^2+2g_{h,\mathrm{}}+g_{\mathrm{}})(2g_{h,\mathrm{}}+g_{\mathrm{}})})(2g_{h,\mathrm{}}+g_{\mathrm{}}\sqrt{(\gamma _{\mathrm{}}^2+2g_{h,\mathrm{}}+g_{\mathrm{}})(2g_{h,\mathrm{}}+g_{\mathrm{}})})}{\gamma _{\mathrm{}}^2\sqrt{(\gamma _{\mathrm{}}^2+2g_{h,\mathrm{}}+g_{\mathrm{}})(2g_{h,\mathrm{}}+g_{\mathrm{}})}}}`$,
* either $`ϵ<ϵ_H={\displaystyle \frac{\gamma _{\mathrm{}}^2(g_{h,\mathrm{}}+2g_{\mathrm{}})}{(\gamma _{\mathrm{}}^2+g_{h,\mathrm{}}g_{\mathrm{}})(g_{h,\mathrm{}}g_{\mathrm{}})}}`$ or $`g_{h,\mathrm{}}<g_{\mathrm{}}`$,
* either $`ϵ<ϵ_h={\displaystyle \frac{\gamma _{\mathrm{}}^2(g_{n,\mathrm{}}+2g_{t,\mathrm{}})}{(\gamma _{\mathrm{}}^2+2g_{h,\mathrm{}}+g_{\mathrm{}}2g_{t,\mathrm{}}g_{n,\mathrm{}})(2g_{h,\mathrm{}}+g_{\mathrm{}}2g_{t,\mathrm{}}g_{n,\mathrm{}})}}`$ or $`2g_{h,\mathrm{}}+g_{\mathrm{}}<2g_{t,\mathrm{}}+g_{n,\mathrm{}}`$.

These results are displayed in Fig. 4. For regions in which two different patterns are stable, i.e. in which different local minima of the energy exist, we have also determined the respective Maxwell values of $`ϵ`$ at which these minima have the same value. For $`ϵ=ϵ_{HS}`$ the energy of the hexagon pattern is equal to the energy of the square pattern, for $`ϵ=ϵ_{HF}`$ it is equal to the energy of the flat surface.

For all susceptibilities either the relation $`g_{h,\mathrm{}}<g_{\mathrm{}}`$ or the relation $`2g_{h,\mathrm{}}+g_{\mathrm{}}<2g_{t,\mathrm{}}+g_{n,\mathrm{}}`$ is violated and consequently hexagons always become unstable for sufficiently large $`ϵ`$.

Of particular importance is the line denoted by $`ϵ_h`$. For $`ϵ>ϵ_h`$ the hexagon pattern is unstable to squares which in turn for $`\chi >\chi _S`$ are not saturated by the fourth order term in the energy. In that case within our perturbation ansatz we are not able to predict which pattern will show up. For $`\chi >\chi _O0.56`$ this transition occurs right at onset and our analysis of hexagons is therefore restricted to $`\chi <0.56`$. This stability criterion for hexagons was to our knowledge not discussed in the literature before. It was not found by Silber and Knobloch , since with their approach the relative stability between hexagons and squares could not be determined. It was also missed by Gailitis  and by Kuznetsov and Spektor  since they considered only the case $`\chi 1`$.

## V Ferrofluid layer with arbitrary thickness

The aim of this section is to investigate how the finite depth $`d`$ of a ferrofluid layer affects the pattern selection. Only ferrofluids with relatively small susceptibilities $`\chi <\chi _S`$ will be considered in order not to leave the region of applicability of our perturbation approach (see the discussion above).

From the linear stability analysis  it is well known that a thin layer retards the onset of the instability and shifts the wavenumber of the unstable mode to larger values. The dependence of the critical magnetic field $`H_{\mathrm{c},d}`$ and of the corresponding critical wavenumber $`k_{\mathrm{c},d}`$ on the layer thickness $`d`$ can be determined from Eqs. (30). It is shown in Fig. 5. If the depth $`d`$ of the fluid is larger then the critical wavelength at infinite depth $`\lambda _{\mathrm{c},\mathrm{}}=2\pi /k_{\mathrm{c},\mathrm{}}`$, which is typically of the order of $`10`$ mm, the finite thickness of the layer can be ignored. Even if the layer is very thin, the parameters characterizing the linear instability of the flat surface are modified only within a few percent for fluids with small susceptibilities.

To start with the non-linear analysis we have first to determine the dependence on $`d`$ of the critical values $`\chi _R,\chi _S`$ and $`\chi _H`$ beyond which higher order terms in the energy functional are necessary. The corresponding results are shown in Fig. 6. All values of the critical susceptibilities get eventually smaller if the layer depth decreases. They do not cross each other such that $`\chi _S`$ is the smallest one for all values of $`d`$. We have included the result for ridges also although these are unstable in the setup of Fig. 1 for all $`d`$ and $`\chi `$. However, ridges can be stabilized by an additional magnetic field tangential to the undisturbed surface  and therefore our results give first informations on the behavior in this modified experimental setup.

For very shallow layers the amplitude of the arising pattern may become larger than the layer depth already at onset. As a consequence the ferrofluid layer disintegrates into unconnected regions and “dry spots” occur. This topological change is outside our theoretical analysis. The corresponding value of $`\chi `$ denoted by $`\chi _d`$ is also shown in Fig. 6. As expected it becomes the smallest of all critical $`\chi `$-values for $`d0`$.

In order to investigate the stability of the patterns arising on a layer with finite thickness we can perform basically the same analysis as in the previous section. We just have to fix the wavenumber $`k`$ to its critical value $`k_{\mathrm{c},d}`$ corresponding to the given thickness $`d`$ and to use the appropriate value for $`H_{\mathrm{c},d}`$ in the function $`l(ϵ,k=k_{\mathrm{c},d})`$ (see Eq. (30)). Moreover we have to take into account the dependence of the coefficients $`\gamma ,g,g_h,g_t`$ and $`g_n`$ on the depth $`d`$. In this way we can finally determine for any given depth $`d`$ the stability regions of the different patterns by means of the corresponding energy function $`f(\{A_n\},k=k_{\mathrm{c},d})`$ as given by Eq. (35).

The results for the various stability boundaries and corresponding Maxwell points for a fluid with susceptibility $`\chi =0.05`$ are shown in Fig. 7. All the critical values of $`ϵ`$ at which transitions between different planforms arise increase in absolute value with decreasing layer thickness $`d`$. Simultaneously the hysteretic behaviour becomes more pronounced. As in the linear analysis we find that the effects of finite depth can be neglected, as long as the fluid layer is thicker then the critical wavelength $`\lambda _{\mathrm{c},\mathrm{}}`$. However, for a thin layer the stability of the hexagons and squares is changed noticeably. For instance the value of $`ϵ_H`$ at which the transition from hexagons to squares takes place may increase by a factor of $`3`$ implying an increase of the magnetic field of about $`10\%`$. As discussed above for very thin layers our theory is no longer applicable since the dips of the hexagonal pattern at onset already reach the bottom of the fluid. The corresponding values of $`d`$ are indicated by the shading in Fig. 7.

Fig. 8 displays the analogous results as obtained for a ferrofluid with the larger susceptibility of $`\chi =0.30`$. Whereas the situation for $`ϵ<0`$ is rather similar to the one shown in Fig. 7, the behaviour for positive $`ϵ`$ is qualitatively different. This is due to the fact that although $`\chi _S(d\mathrm{})>\chi `$ with decreasing depth $`d`$ both $`\chi _S(d)`$ and $`\chi _O(d)`$ eventually become smaller than $`0.30`$ (cf. Fig. 6). This gives rise to the divergencies in the energy discussed above and changes the stability chart accordingly.

## VI Wavenumber selection

In many pattern forming systems the wavelength of the first unstable mode is often a good estimate for the typical length scale of the emerging pattern. If a whole band of modes is unstable the one with maximal growth rate is likely to dominate the arising structure . On the other hand the developed pattern is largely determined by the non-linearity in the problem and therefore its length scale may also substantially deviate from the relevant scales of the linearized theory.

A particularly gratifying feature of our variational approach to pattern formation in the normal field instability in ferrofluids is the possibility to include the wave number modulus $`k`$ into the set of parameters varied to minimize the energy functional. This allows to directly determine the dependence of the pattern periodicity on the external magnetic field and to investigate the influence of the variation in wavenumber on the stability of the different patterns.

We therefore investigate in the present section the wavenumber selection problem for the Rosensweig instability. For simplicity we will consider a ferrofluid with infinite depth only. As shown in the previous section this is a very good approximation as long as $`d>\lambda _{\mathrm{c},\mathrm{}}`$. We will compare our findings with recent experiments in which this condition is met.

The necessary calculations are similar to those reported in section IV with the only extension that we no longer fix the wavenumber modulus $`k`$ in advance but include it into the set of parameters with respect to which the energy is to be minimized.

A typical result is shown in Fig. 9 displaying the maximal surface deflection $`\zeta (0,0)`$ and the corresponding wavenumber $`k`$ of the patterns. As can be seen from the figure, the wavenumber of the patterns decreases with increasing supercriticality. This is in particular pronounced for the square pattern and implies that in this case the system may lower its energy very efficiently by increasing the wavelength of the pattern.

This effect introduces an additional threshold for $`ϵ`$ into the stability chart shown in in Fig. 10. If $`ϵ>ϵ_\lambda `$ we find $`k0`$ and the pattern disappears by increasing its wavelength to infinity. This process is accompanied by an unbounded increase of the amplitudes $`A_1,A_4`$ which carries us outside of the validity of our perturbation ansatz (22). The stability regions of the other pattern remain qualitatively the same as without minimization in $`k`$ as can be seen from a comparison between figures 10 and 4.

For supercritical magnetic fields a whole band of wavenumbers giving rise to stable patterns exists. For hexagons this stability region is shown in Fig. 11. At fixed $`ϵ`$ the hexagonal planform becomes unstable if the wavenumber is increased sufficiently. In  such an increase of the wavenumber was achieved experimentally by compressing the hexagonal pattern in a hopper. The observed transition from a hexagonal to a square pattern with smaller wavenumber occurs when the limit $`ϵ_h`$ of the stable $`k`$-band is exceeded and is hence in qualitative agreement with our theory.

## VII Discussion

In the present paper we have theoretically investigated the formation of patterns on the free surface of a ferrofluid resulting from the Rosensweig instability. For a ferrofluid layer of arbitrary depth $`d`$ in an external magnetic field $`𝐇_0`$ the equilibrium surface profile was determined by minimizing the appropriate thermodynamic potential. We have assumed a linear dependence of the magnetization on the magnetic field which is at most 5% off the correct values for the experimentally situations with which we compare our findings. The analysis is restricted to the vicinity of the critical magnetic field at which the flat interface becomes unstable since it uses a perturbation expansion for the energy up to fourth order in the amplitude of the surface deflection. Our findings generalize the classical results of Gailitis  and Kuznetsov and Spektor  to layers of finite depth $`d<\mathrm{}`$.

We have found that, at least for the slightly supercritical magnetic fields considered, the finite thickness of the ferrofluid layer can be neglected as long as it remains larger than the critical wavelength of the linear instability which in typical experiments is about 1 cm. For thinner layers nonlinear effects become increasingly important. This shows up quite generally in more pronounced hysteresis effects for the transition between different planforms. We also find smaller values of the susceptibilities for which squares and ridges may appear through backward bifurcations and correspondingly higher order terms in the energy are necessary to saturate the instabilities. The critical values of the susceptibility $`\chi `$ beyond which our theory breaks down are displayed in Fig. 6. Together with the findings shown in Fig. 3 they quantify the qualitative statement $`\chi 1`$ used in . An extension of our analysis including higher orders similar to what has been done in  for the one-dimensional situation seems feasible but rather tedious.

It is possible to understand the enhanced hysteretic behaviour for thin layers qualitatively. A thin layer inhibits surface deformations resulting from the Rosensweig instability (see Fig. 5) since the magnetic field energy is suppressed due to the finite depth of the layer. On the other hand for the developed pattern the magnetic flux is concentrated at the peaks of the surface so that the effective layer depth becomes larger. Correspondingly the pattern is additionally stabilized and decays only after a reduction of the field stronger than expected from the linear theory. A similar reasoning explains why backward bifurcations of ridges and squares are facilitated in thin layers.

In addition to determining the stable planforms by minimizing the energy functional in the amplitudes of the different modes present in our ansatz for the surface profile, we have also addressed the problem of wavenumber selection by including the wavenumber modulus $`k`$ into the set of variational parameters. We found quite generally a decrease of the wavenumber of the pattern with increasing magnetic field. This is in qualitative agreement with recent experiments by Abou et al.  in which the the field intensity was abruptly increased to overcritical values and the corresponding wavelength $`\lambda =2\pi /k`$ of the resulting hexagonal pattern was found to increase with increasing field. Likewise it was found that the wavenumber $`k`$ of the hexagonal pattern is already for slightly supercritical magnetic fields smaller then the critical wavenumber $`k_\mathrm{c}`$. The corresponding difference $`\mathrm{\Delta }k`$ is also found in the theory and shown in Fig. 9. Unfortunately a quantitative comparison with the experimental results is impossible since in the experiment a ferrofluid with a susceptibility of $`1.4>\chi _H`$ was used. For such a susceptibility the theoretical analysis would require higher order terms in the perturbation expansion of the energy. Nevertheless we found that $`\mathrm{\Delta }k`$ grows with increasing susceptibility such that the experimentally observed value of a few percent agrees well with our result for $`\mathrm{\Delta }k`$ at smaller susceptibilities $`\chi <\chi _H`$.

Moreover the theoretical analysis has shown that when including the wavenumber into the set of variational parameters the square pattern gets rather susceptible to an unbounded increase in the wavenumber. This gives rise to a reduced stability domain for squares and may be related to the fact that in the experiments using jumps in the field intensity always hexagonal arrays of peaks were found . Quite generally our results verify that the resulting planform may depend on the details of the experimentally chosen way to reach overcritical magnetic fields.

Let us finally stress that the variation of the wavenumber of the patterns with the external field as found here theoretically and observed in the experiments does not invalidate our ansatz with a fixed wavenumber used before. As reported in  the wavenumber of the arising hexagonal pattern coincides with $`k_\mathrm{c}`$ if the magnetic field is increased and decreased in a quasi-static way. This is probably due to the fact that the boundary conditions suppress the emergence or disappearance of peaks in a developed pattern. Therefore using the quasi-static process a metastable pattern with $`k_\mathrm{c}`$ is produced whereas a jump in field intensity gives rise to the most stable pattern corresponding to a smaller wavenumber.

###### Acknowledgements.

We would like to thank Bérengère Abou for explaining her experimental findings to us and Adrian Lange for interesting discussions. This work was supported by the Deutsche Forschungsgemeinschaft under the project En 278/2-1.


format: text

text:
# Acknowledgement

## Acknowledgement

I am very much indebt to Prof.P.S.Letelier, Prof.Harald Soleng,Prof.J.G.Pereira and Prof.I.D.Soares for helpful discussions on the subject of this paper.Financial support from CNPq. and UERJ is gratefully acknowledged.


text:
# Superluminal motion in a compact steep spectrum radio source 3C 138

## 1 Introduction

The radio source 3C 138 (= 4C 16.12 = J0521+1638) is identified with a quasar with *m*v=18.84 and *z*=0.759 (Hewitt & Burbidge hewitt (1989)). It is a compact, powerful quasar with a well-defined turn-over in the spectrum at about 100 MHz and a steep high frequency spectrum of 0.65 (S $``$ $`\nu ^\alpha `$), making it a prototype of the compact steep spectrum (CSS) radio sources (Fanti et al. fanti90 (1990)). Its spectrum remains straight up to 22 GHz, supporting that no dominant flat spectrum core exists (Kameno et al. kameno (1995)).

The radio structure of 3C 138 consists of a core with a bright jet and compact lobe to the east and a fainter, more diffuse lobe to the west from the VLA and MERLIN observations (Akujor et al. akujor (1993) and references therein; Lüdke et al. ludke (1998)). High resolution VLBI observations (Geldzahler et al. geldzahler (1984); Fanti et al. fanti89 (1989); Nan et al. nan (1992); Dallacasa et al. dallacasa (1995); Cotton et al. cotton (1997)) reveal several knots within the main jet emission extending about 400 mas in a position angle of 65 and a low brightness counter-jet about 250 mas away in the opposite direction. There are two knots dominating the central 10 mas region. Fanti et al. (fanti89 (1989)) claimed on the basis of apparent size that the eastern knot is the core. Cotton et al. (cotton (1997)) argued that the western one is the location of the nuclear core mainly because of its very weak polarization ($`<0.4\%`$) compared to a peak polarized intensity of $`3.5\%`$ for the eastern component at 5.0 GHz. Cotton et al. (cotton (1997)) also suggested a relative proper motion of 0.27 mas yr-1 between the two knots, corresponding to an apparent velocity of 9.7c (assuming H0 = 65 km s-1 Mpc-1 and q0 = 0.5, hereafter).

In general, the radio emission from CSS sources is thought to be relatively free of beaming effects, consistent with the fact that CSS sources are the least variable of compact extragalactic radio sources, and CSS sources are rarely found to be superluminal (O’Dea odea (1998)). When combined with the observed two-sided jet structure, such a superluminal motion in 3C 138, if confirmed, will strictly confine the jet velocity and its direction to the line of sight of the observer. We report on 5.0 GHz VLBI observations of 3C 138 as an attempt to confirm its proper motion and probe its compact radio core area.

## 2 Observations and data reduction

The VLBI observations were carried out at 5.0 GHz on November 8, 1997 with the enlarged European VLBI Network composed of Effelsberg (100-m), Hartebeesthoek (26-m), Jodrell MK2 (25-m), Medicina (32-m), Noto (32-m), Onsala (25-m), Shanghai (25-m), Torun (32-m), Urumqi (25-m) and Westerbork (25-m). Sources were observed in a snapshot mode, i.e., 5-12 thirteen-minute scans spreading across a wide range in hour angle to obtain a good u-v coverage. The left-circular polarized radio signals were recorded in Mark III mode C with a total bandwidth of 14 MHz (7 IFs and 2 MHz per IF band) at each station. The correlation of the data was performed on the MPIfR MK III Correlator in Bonn, Germany with an output averaging time of 4 s.

Then, a priori visibility amplitude calibrations for each antenna were applied using the antenna gain and the system temperatures measured during the observations. After the data were globally fringe-fitted in AIPS (Schwab & Cotton schwab (1983)), the subsequent data inspection, imaging and modelfitting were carried out with DIFMAP (Shepherd shepherd (1997)).

Because of the large separation ($``$400 mas) between the eastern compact lobe and central core region in 3C 138, special attention was paid to the smearing effects of time averaging and bandwidth in order to image an undistorted full field. The 2 MHz data in each IF were not averaged through 7 IFs, and the averaging time was only 12 s to minimize smearing effects.

Images were produced with the usual iterations of cleaning and phase self-calibration. Only a constant gain factor for each IF band of each antenna was implemented in the later stage of imaging. These gain corrections are consistent with those obtained for other sources (Shen et al. in preparation). To improve the sensitivity to the low surface brightness regions, the visibility data were naturally weighted and tapered with FWHM = 30M$`\lambda `$.

## 3 Results

### 3.1 Overall structure

The overall morphology of 3C 138, as shown in Fig. 1, resembles past VLBI images. There are two notable emission regions: the central core and the eastern lobe, which are separated by 400 mas at a position angle of 70. In between them are several discrete emission components.

The total cleaned flux density in the image (Fig. 1) is $``$0.85 Jy, which is the same as the correlated flux density on the shortest baseline, Westerbork-Jodrell MK2, but far less than the values of 3.24 Jy and 3.58 Jy by Dallacasa et al. (dallacasa (1995)) and Cotton et al. (cotton (1997)), respectively. The single-dish measurements from the University of Michigan Radio Astronomy Observatory (UMRAO) show a very stable total flux density of 3.80 Jy at 4.8 GHz with a variation less than 6% over past 10 years. The source is heavily resolved on baselines longer than 3 M$`\lambda `$ (see Fig. 2c in Fanti et al. fanti89 (1989)). So, the main reason for the shortage of cleaned flux density in our image is the lack of shorter baselines like MERLIN or VLA-Pie Town (cf. Dallacasa et al. dallacasa (1995)). As a result, the more extended, lower surface brightness regions seen in previous images (e.g. Fanti et al. fanti89 (1989)) are not well recovered in our image, which causes discontinuity in the main jet emission towards the eastern lobe. The connection of separated emissions (including the core region), however, exhibits a sign of oscillation seen from other observations. No counter-jet was significantly detected in our image though there are some very weak diffuse emissions at about 3 times noise level around the counter-jet region of the 1.7 GHz image by Cotton et al. (cotton (1997)). The RMS noise in the image is 1.2 mJy/beam, about twice the expected thermal noise for the total on-source time of 1.3 hr.

Two compact components (labeled F and G after Fanti et al. fanti89 (1989)) in the eastern lobe can be clearly seen. They are aligned nearly orthogonal to the jet axis, which may be interpreted as a working surface of the jet due to its interaction with the surrounding medium.

### 3.2 Central core region

For the purpose of identifying the center of activity, a full resolution image is displayed in Fig. 2. This image was made from the same calibrated data as Fig. 1, except that the baselines shorter than 30 M$`\lambda `$ were removed and the uniform weighting was adopted to ensure a high resolution. To yield a quantitative description, modelfitting to both amplitudes and phases in the calibrated visibility data was applied. Because of the complicated overall structure and the limited observational data, we fit a model to core region only while keeping the cleaned components for the remaining structure (cf. Piner et al. piner (2000)). In this way, we can minimize the effect mainly from the components F and G. We also tried a modelfitting to the longer baselines data corresponding to Fig. 2. We even made an image of the central region only with a 1-minute averaging time and then fit the visibility data. All these results are quite consistent. In Table 1 is presented the best fitting model consisting of three circular Gaussian components.

The sizes of components CE and CW are about the same, but CE is much brighter. The flux density of CE is stable over past 12 years (cf. Akujor et al. akujor (1993); Dallacasa et al. dallacasa (1995)). The two point spectral indices ($`\alpha _{5.0\mathrm{GHz}}^{8.6\mathrm{GHz}}`$) obtained by combining the high resolution VLBA results at 8.6 GHz (Fey & Charlot fey (1997)) are almost the same for both CE and CW (see Table 1). Component CE2 is less well constrained because it is much too resolved. Actually, CE2 emission has appeared in several previous 5 GHz images (e.g. Fig. 13 in Cotton et al. (cotton (1997))), but was first fitted at 8.6 GHz by Fey & Charlot (fey (1997)).

## 4 Discussion

### 4.1 Proper motion

There are three published 5 GHz VLBI measurements made at epochs 1985.5 (Fanti et al. fanti89 (1989)), 1989.72 (Dallacasa et al. dallacasa (1995)) and 1994.97 (Cotton et al. cotton (1997)). Fanti et al. (fanti89 (1989)) resolved the central region with two compact components, but did not present the separation at the time. Dallacasa et al. (dallacasa (1995)) made the first measurement of the apparent separation of 5.0 mas. Cotton et al. (cotton (1997)) estimated the separation of 3.9 and 6.5 mas for epochs 1985.5 and 1994.97 respectively, and obtained an averaged rate of separation of 0.27 mas yr-1 from all three-epoch observations (see Table 1 and Fig. 14 of Cotton et al. (cotton (1997))). A separation of 6.5 mas at epoch 1994.97 is larger than that of 5.76 mas at a later epoch 1997.86 in Table 1. This would be indicative of a back flow in the jet motion, which is very unlikely if not impossible. The 8.6 GHz VLBA observations (Fey & Charlot fey (1997)) measured a separation of only 5.4 mas at epoch 1995.78, which is also much less than that at epoch 1994.97. To further investigate this, we measured the peak-to-peak separation from images directly, and then compared them with those modelfitting results. As a result, we got an excellent agreement between the modelfitted separation and the peak-to-peak separation at epochs 1989.72 and 1997.86 at 5.0 GHz and at epoch 1995.78 at 8.6 GHz, but the discrepancy at epochs 1985.5 and 1994.97 was as large as 0.5 and 0.9 mas, respectively. This difference is about one-half and one-quarter of the corresponding beam size, and quite difficult to understand for us at present.

Therefore, we decided to estimate the proper motion by using two epoch (1989.72 and 1997.86) 5 GHz data (denoted by crosses in Fig. 3) only. As shown in Fig. 3, this produces a relative proper motion of 0.093 mas yr-1 between CE and CW, which corresponds to a superluminal speed of 3.3c. We used other two 5.0 GHz observations to evaluate its uncertainty. The peak-to-peak separation (denoted by open squares in Fig. 3) of 4.4 and 5.6 mas at the corresponding epochs 1985.5 and 1994.97 gave an uncertainty of 5%, while the published points (denoted by filled squares in Fig. 3) would introduce an uncertainty of 16%. This proper motion is much slower than the reported 0.27 mas yr-1 (Cotton et al. cotton (1997)). The ongoing monitoring observations confirm that the knot is definitely moving much slower than suggested by the earlier data (Bill Cotton, private communication 2000).

### 4.2 Synchrotron self-absorption

There is a relative position offset of 0.17 mas at epoch 1995.78 for CE between the position interpolated from 5 GHz data (fitting line in Fig. 3) and the measurement made at 8.6 GHz (Fey & Charlot fey (1997)) denoted by a open triangle in Fig. 3. Such a frequency dependence of the observed position can be explained by synchrotron self-absorption (Lobanov lobanov (1998)). To estimate this position offset using eq. (11) of Lobanov (lobanov (1998)), we have searched the literature for its radio spectral energy distribution. As a result, the total synchrotron luminosity of 1.5$`\times `$1045 erg s-1 is estimated. A superluminal speed of 3.3c implies a maximum angle $`\theta `$ to the line of sight, a minimum jet speed $`v(=\beta `$c) and a minimum Lorentz factor $`\gamma `$ of 34, 0.96c and 3.5, respectively. The projected opening angle $`\varphi _\mathrm{o}`$ of 5.8 is measured from CE component size over its distance to CW at 8.6 GHz. These values predict a maximum position offset of 0.22 mas between 5.0 and 8.6 GHz under the assumption of energy equipartition between radiating particles and magnetic field. The upper limit is mainly due to the uncertainty in the determination of the viewing angle $`\theta `$ and the resultant lower limits to both $`\beta `$ and $`\gamma `$. Thus, such an observed offset in the central region of 3C 138 can be well explained as a frequency dependent shift of the synchrotron self-absorbed core. Assuming a conical jet geometry in which the magnetic field decreases linearly with its distance to the core, we can further derive the magnetic field in the 5 GHz VLBI core region of 0.12, 0.02 and 0.01 G for components CW, CE and CE2, respectively (cf. Lobanov lobanov (1998)).

Recently, a significant position offset has been reported for another superluminal CSS source 3C 216 (see Fig. 3 of Paragi et al. paragi00 (2000)). These may suggest that the core region in CSS sources could also be described by ultracompact jet models, which have been successful in interpreting the central region of radio-loud active galactic nuclei (cf. Lobanov lobanov (1998) and references therein).

We are aware, however, that this offset may be due to the positional uncertainties of measurements and limited resolution and blending at lower frequency. Future high resolution, multi-frequency VLBI observations will be critical to clarify this.

### 4.3 Core identification

It is hard to determine the real core in 3C 138 based on the apparent source size and spectral index available now. From their high resolution polarization sensitive imaging, Cotton et al. (cotton (1997)) have ruled out the possibility that CE was the core and further suggested that CW harbors the central energy source. The non variation in the flux density of CE, discussed in Sect. 3.2, can also be viewed as an indication that CE probably is not related to the center of activity. Future observations with better resolution at lower frequencies (such as the 1.6 GHz space VLBI (VSOP) imaging observations (cf. Hirabayashi et al. hirabayashi98 (1998))) and better sensitivity at higher frequencies should be able to clearly separate and detect the two compact components CW and CE and, thus enable us to undoubtedly determine if the spectral index, as suggested by Cotton et al. (cotton (1997)), is flatter in component CW.

Our new measurements confirm the systematic increasing in the knot separation with time. But the rate is three times slower than the previous one. The detected superluminal speed, when combined with the observed two-sided jet structure in 3C 138, poses very strong constraints on the jet kinematics.

The expected ratio of the flux density between jet and counter-jet is given by $`\mathrm{R}=(\frac{1+\beta \mathrm{cos}\theta }{1\beta \mathrm{cos}\theta }`$)2+α. Here, intrinsically symmetrical two-sided jets with the same speed in two opposite directions are assumed.

The ratio of the flux density of the eastern to western lobe seen in the 1.7 GHz image (Cotton et al. cotton (1997)) is about 16. Using the spectral index between 1.7 and 5.0 GHz of 0.75 (Fanti et al. fanti89 (1989)), we obtained the product $`\beta \mathrm{cos}\theta =0.46`$, in good agreement with the value estimated at 5 GHz (Fanti et al. fanti89 (1989)). But this restricts to a maximum apparent speed of 1.6c, which is significantly smaller than 3.3c. Therefore, the superluminal motion of 3.3c in the central core region, though much less than earlier reported one, is still too high to allow the overall counter-jet structure to be seen in terms of the projection effect and Doppler boosting of an intrinsically identical jet.

One solution to circumvent this is to assume a non-uniform distribution of R for the whole jet since usually lobes are observed to move at a much lower speed. This can naturally explain the absence of parsec-scale counter-jet in the central region. Cotton et al. (cotton (1997)) obtained a lower limit to the ratio R of 20 within the central area. Our new observations slightly raised this lower limit to 30 at 5 times RMS noise level. These lower limits are consistent with the detected superluminal motion in the inner region being due to the relativistic boosting in the jet.

Such a change in R could be caused by the strong interaction of the jet with the surrounding medium in the immediate region of the core. This is supported by the existence of a large amount of ionized gas in the central region (Cotton et al. cotton (1997)). The bending in the jet or a helical jet model is a possible explanation, but unlikely to affect R significantly in case of 3C 138 (cf. Fanti et al. fanti89 (1989)). A superluminal speed of 3.3c implies a ratio R greater than 360 (using the spectral index of component CE in Table 1), 12 times the lower limit estimated above. On the other hand, the possibility that the observed overall asymmetry is intrinsic is still there, but difficult to test because of the limited data gathered. The nature of the core region is worthy of a more detailed study.

###### Acknowledgements.

We wish to thank Carla Fanti, the referee, for a careful reading of the manuscript and very useful comments, and Bill Cotton for providing us with the 1.7 GHz VLBI image and for communication of information prior to publication. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with the National Aeronautics and Space Administration. This research has made use of data from the University of Michigan Radio Astronomy Observatory which is supported by the National Science Foundation and by funds from the University of Michigan.


text:
# Fast Ground-State Reordering from Off-Diagonal Fluctuations of Interaction

\[

## Abstract

It is shown that off-diagonal fluctuations of interaction lead to the breakdown of perturbation theory in the immediate vicinity of the ground-state of a fermionic system at a rate inversely proportional to the number of considered particles. This rate is parametrically much faster than previously expected. It is shown that this follows from the coherent addition of many small second order contributions to the energy of low-lying levels.

\]

Random interactions naturally occur in various fields of physics. Assuming the symmetry of the probability distribution of the one-body Hamiltonian under unitary basis transformation in Hilbert space directly results in the randomness and a similar symmetry for the distribution of the Interaction Matrix Elements (IME). Typically, this symmetry is satisfied up to an energy scale $`E_c`$ (the Thouless energy) so that models with random interactions may be expected to describe interacting systems provided one restricts oneself to energy scales of an order of magnitude given by $`E_c`$ . In condensed-matter systems, this corresponds to one-body excitation energies inside a window of width $`E_c=g\mathrm{\Delta }`$ around the Fermi energy, hence containing a finite number of levels given by the conductance $`g`$ ($`\mathrm{\Delta }`$ is the mean level spacing).

Random Interaction Models (RIM) first appeared in nuclear physics and early investigations focused on the deviations from Random Matrix Theory introduced by considering interactions of finite order $`kn`$ ($`n`$ is the number of particles, for $`k=n`$ one recovers the Dyson ensembles of random matrices) . The complexity of the excited spectrum of rare-earth atoms motivated the introduction of RIM in atomic physics , and RIM reproduced statistical features obtained in nuclear shell model calculations . Last but not least, RIM were used in condensed matter physics to study the quasiparticle lifetime  and the conductance peak spacings and heights distributions  in quantum dots in the Coulomb blockade regime. More recently a lot of attention has been devoted to the properties of the ground-state of RIM  and several spectral features presented in these latter investigations indicate that RIM are realistic down to very low excitation energies. It is the purpose of the present letter to investigate in detail the ground-state structure of systems of few ($`n`$) fermions as the relative strength $`U/\mathrm{\Delta }`$ of fluctuations of interaction increases. The main result of these investigations is a breakdown of perturbation theory for the ground-state at a critical value $`U_c(n)\mathrm{\Delta }/n`$. In particular, $`U_c`$ is parametrically smaller by a factor $`1/n`$ than previously suggested . Once this threshold is reached, the ground-state starts to be composed of an exponentially large number of Slater determinants, an occurrence which goes far beyond a simple renumbering of one-particle orbitals. We will see that this is due to the coherent addition of large numbers of small second order contributions for perturbation theory in the immediate vicinity of the ground-state.

We investigate a model for spin-$`1/2`$ fermions

$`H`$ $`=`$ $`{\displaystyle ϵ_\alpha c_{\alpha ,s}^{}c_{\alpha ,s}}+{\displaystyle U_{\alpha ,\beta }^{\gamma ,\delta }c_{\alpha ,s}^{}c_{\beta ,s^{}}^{}c_{\delta ,s^{}}c_{\gamma ,s}}`$ (1)
The greek indices label $`m/2`$ different one-body energies which are distributed as $`ϵ_\alpha [0;m/2]`$ so that the mean spacing between spin-degenerate levels is $`\mathrm{\Delta }1`$ and $`s^{(^{})}=,`$ are spin indices. The Hamiltonian is Spin Rotational Symmetric (SRS): the interaction commutes with both the total spin $`|\stackrel{}{S}|`$ and its projection $`S_z`$. The corresponding eigenvalues $`\sigma `$ and $`\sigma _z`$ are good quantum numbers and the Hamiltonian matrix acquires a block structure where blocks are labelled by $`\sigma _z`$ and $`\sigma `$. From now on we consider only the $`\sigma =\sigma _z=0`$ (1/2) block for $`n`$ even (odd). The results and arguments presented below can however be generalized to higher spin blocks or models of spinless fermions.

The Hamiltonian (1) can be viewed as a model of interacting fermions expressed in the basis of Slater determinants constructed from the set of one-particle eigenfunctions $`\psi _\alpha `$ of $`H_0ϵ_\alpha c_{\alpha ,s}^{}c_{\alpha ,s}`$. In this basis, the IME are given by $`U_{\alpha ,\beta }^{\gamma ,\delta }=_{i,j}U(ij)\psi _\alpha ^{}(i)\psi _\beta ^{}(j)\psi _\gamma (i)\psi _\delta (j)`$, where $`U(ij)`$ is the interaction potential in real space; $`i`$ and $`j`$ label sites on a lattice. Assuming that the presence of disorder or chaotic boundary scattering confers a random character to the $`\psi _\alpha `$’s directly leads to fluctuations of the $`U_{\alpha ,\beta }^{\gamma ,\delta }`$ around their average value. Only diagonal matrix elements have a nonzero average leading to mean-field charge-charge, spin-spin and BCS interactions  which we will neglect as they have at most a marginal influence on the object of interest here. We therefore take the IME in (1) to be randomly distributed inside a zero-centred Gaussian distribution of width $`U`$, $`P(U_{\alpha ,\beta }^{\gamma ,\delta })e^{(U_{\alpha ,\beta }^{\gamma ,\delta })^2/2U^2}`$.

We start with a perturbative calculation of the ground-state energy up to the second order in the small parameter $`U/\mathrm{\Delta }1`$. The first order correction vanishes on average, and the second order correction can be written

$$\mathrm{\Delta }_0^{(2)}=\overline{}\frac{(U_{\alpha ,\beta }^{\gamma ,\delta })^2}{ϵ_\gamma +ϵ_\delta ϵ_\alpha ϵ_\beta }$$
(2)
$`\overline{}`$ indicates that the sum is taken over occupied levels $`\gamma `$ and $`\delta `$ and unoccupied ones $`\alpha `$ and $`\beta `$ (taking into account spin degeneracy), giving a total number of contributions $`Kn^2(mn)^2`$. For the ground-state, one has $`ϵ_\gamma +ϵ_\delta <ϵ_\alpha +ϵ_\beta `$, and each term in the sum (2) is negative. The contributions from each scattering process therefore add coherently resulting in a strong reduction of the ground-state energy. The latter can be estimated for an equidistant spectrum by replacing sums by integrals in (2). This gives an homogeneous polynomial of order three in $`n`$ and $`m`$, each term being multiplied by a logarithmic correction. In the dilute limit $`1nm`$ the $`m^3`$ and $`m^2n`$ terms drops out exactly so that the dominant contribution to the second order correction in the ground-state energy reads

$$\mathrm{\Delta }_0^{(2)}A(U^2/\mathrm{\Delta })n^2m\mathrm{log}(m)$$
(3)
The prefactor $`A`$ can be estimated $`A7`$; details of the perturbative treatment are presented elsewhere .

In Fig.1 we show the behavior of the ground-state energy for $`n=2,`$ 4, 6 and 8. It is clearly seen that the average ground-state energy exhibits the quadratic dependence in $`U`$ predicted by (2-3) at weak $`U`$, while at larger $`U`$, this dependence becomes linear. Most importantly, this crossover occurs at a weaker $`U`$ for larger $`n`$. Both these features can be understood as one realizes that for dominant fluctuations of interaction $`U/\mathrm{\Delta }1`$ the Hamiltonian depends linearly in $`U`$. In this regime, and for $`n1`$, the average many-body density of states is Gaussian with a width $`B`$ well approximated by $`BnmU`$ . This width thus gives a good parametric estimate for the average ground-state energy in the asymptotic regime $`_0^{(2)}BnmU`$  so that equating the asymptotic and perturbative expressions for the ground-state energy ($`\mathrm{\Delta }_0^{(2)}nmU`$) gives a critical value

$$U_c\mathrm{\Delta }/n$$
(4)
The crossover between the two regimes occurs at a rate given by $`U_c`$ and is inversely proportional to the number of particles. As can be seen in Fig.1, below this threshold, the ground-state energy is well approximated by the second order perturbation expression (3), while above it, its $`U`$-dependence becomes linear which cannot be captured by any finite order of perturbation theory, and we are forced to conclude that $`U_c`$ gives the radius of convergence for perturbation theory in $`U`$.

For $`U>U_c`$, the ground-state starts to be mixed with higher excited states and the number of Slater determinants $`|J`$ necessary to construct it becomes much larger than one. A measure of this number is given by the Participation Ratio (PR), which, for a given many-body eigenstate $`|\mathrm{\Psi }_\lambda (U)`$, is defined as $`\xi _\lambda =(_J|J|\mathrm{\Psi }_\lambda (U)|^2)^2/_J|J|\mathrm{\Psi }_\lambda (U)|^4`$. If $`|\mathrm{\Psi }_\lambda (U)`$ is made up of one single Slater determinant (as it is for $`U=0`$), one has $`\xi _\lambda =1`$ whereas when $`|\mathrm{\Psi }_\lambda (U)`$ is homogeneously spread over the $`N`$ Hilbert space basis states, $`\xi _\lambda N`$. To check the above threshold (4), we numerically calculated the PR for the ground-state wavefunction $`|\mathrm{\Psi }_0(U)`$. Results obtained from the exact diagonalization of small systems of $`n8`$ particles are shown in Fig.2. We first note a faster increase of $`\xi _0`$ with $`U`$ at larger $`n`$ which is not due to an increase of the Hilbert space size as it occurs before a significant portion of the Slater determinant basis starts to be occupied. Most importantly, this increase is exponential (note the vertical log-scale on Fig.2) : fixing e.g. $`U/\mathrm{\Delta }=0.2`$, one gets $`\xi (n)\mathrm{exp}(n/1.6)`$. Since each order in perturbation theory gives only an algebraic number of contributions, this exponential increase of the PR cannot be captured by a finite, $`n`$-independent number of orders, and therefore indicates the divergence of the perturbation series as $`n`$ increases. Another important feature of Fig.2 is the quite small value $`U\mathrm{\Delta }`$ at which $`\xi _0`$ starts to increase. This is quite surprising as $`\mathrm{\Delta }`$ is the only energy scale close to the ground-state. When the PR increases beyond $`\xi _02`$, second (and higher) order terms become more and more important, until eventually perturbation theory breaks down. One can thus extract parametrically the radius of convergence of perturbation theory from the PR by defining a critical interaction strength $`U_c`$ from the condition $`\xi _0(U_c)=2`$. (This choice is arbitrary as long as one chooses a threshold value $`\xi _0(U_c)`$ much smaller than the total size of the Hilbert space as this would introduce an artificial $`n`$ and $`m`$ dependence in $`U_c`$.) The inset to Fig.2 shows a $`n`$ (and absence of $`m`$) dependence of $`U_c`$ clearly confirming (4).

As yet another check of the threshold (4), we consider the onset of level repulsion for the spacing $`s_1=_1_0`$ between the ground-state energy $`_0`$ and the first excited level $`_1`$. At $`U=0`$, this spacing corresponds to a one-body excitation, and considering a randomly distributed one-body spectrum, it is Poisson distributed $`P_P(s_1)=\mathrm{exp}(s_1)`$. Switching on $`U`$ induces off-diagonal matrix elements connecting the two lowest levels and leads to repulsion between them. Making a two-level approximation incorrectly gives the onset of this repulsion at $`U=\mathrm{\Delta }`$. A correct treatment goes as follows. Most terms in the second order contribution (2) are common to both the energies of the ground- and first excited states. The difference corresponds to scattering from or onto the $`(n/2)^{\mathrm{th}}`$ and $`(n/2+1)^{\mathrm{th}}`$ orbitals, giving a contribution to the energy difference between the lowest two levels

$$\overline{}\frac{(U_{\alpha ,n/2}^{\gamma ,\delta })^2}{ϵ_\alpha +ϵ_{n/2}ϵ_\gamma ϵ_\delta }\frac{(U_{\alpha ,n/2+1}^{\gamma ,\delta +1})^2}{ϵ_\alpha +ϵ_{n/2+1}ϵ_\gamma ϵ_{\delta +1}}$$
(5)
(There is an additional similar contribution arising from the scattering from the $`n/2`$ and $`n/2+1`$ orbitals.) We first make the substitution $`\sigma ((U_{\alpha ,n/2}^{\gamma ,\delta })^2(U_{\alpha ,n/2+1}^{\gamma ,\delta +1})^2)=U^2`$ for IME fluctuations. Then, and considering an equidistant spacing, a simple counting shows that the number of contributions in the above expression corresponding to a given energy denominator $`p\mathrm{\Delta }`$ grows as $`p^2`$ up to $`p=n/2`$ and decreases afterward. Neglecting logarithmic corrections, one gets an estimate for the contributions to the relative fluctuations given by $`(U^2/\mathrm{\Delta })n^2`$. This estimate can also be obtained by replacing the sums in (5) with integrals. Once these fluctuations become of the order of the spacing between $`_0`$ and $`_1`$, $`(U^2/\mathrm{\Delta })n^2\mathrm{\Delta }`$, an avoided crossing occurs. This signals the onset of level repulsion which therefore occurs at the rate given by (4).

This reasoning is confirmed on Fig.3 (a,b,c) where we present the distribution $`P(s_1)`$ of the first level spacing for Hamiltonian (1) and different number of particles $`n`$ at fixed interaction strength $`U=\mathrm{\Delta }/n`$. The two main features that can be seen are : 1) $`P(s_1)`$ at fixed $`U=U_c`$ is $`n`$-independent, except that 2) $`P(s_1)`$ shows a clear even-odd dependence, in particular there is more level repulsion for odd $`n`$. This second fact is understood once one realizes that due to SRS, and within a given $`(\sigma _z,\sigma )`$-block, the first $`U=0`$ excitation is doubly degenerate for $`n`$ odd, while it is not degenerate for $`n`$ even. This results in a stronger effect of the interaction in the former case and the even-odd effect seen on Fig.3 (a,b,c). It is to be stressed that this is purely a spin effect, this even-odd dependence does not appear for spinless fermions. On Fig.3 (d) we also show a plot of the low-lying part of the spectrum which shows an avoided crossing at a value $`U/\mathrm{\Delta }0.05`$ much too small to be explained by a two-level approximation. We stress that in (d) the avoided crossing is not accompanied by an increase of the level spacing as would be the case in presence of non-fluctuating IME, and which would delay the breakdown of perturbation theory.

Equation (4) predicts a parametrically small radius of convergence for perturbation theory beyond which the ground-state is composed of exponentially many Slater determinants. It is somehow in opposition with earlier works  which, based on extrapolations from higher excitation energies, predicted a border $`\mathrm{\Delta }`$. The mechanism leading to the breakdown of perturbation theory at higher excitation energies is however due to mixing between two-body levels which are very close in energy, and thus give the dominant contribution in a second order approach. In the immediate vicinity of the ground-state these processes with small denominators do not exist as the lowest energy scale is $`\mathrm{\Delta }`$ \- the approaches used in  correspond then to a two-level approximation involving ground- and first excited states. The presented theory show that the correct procedure is to retain a large number of small second order contributions which have an enhanced effect close to the ground-state, as they acquire a well-defined sign and hence add coherently to give the small threshold (4). One has to conclude that the approach used in , while fully justified at a sufficient excitation energy, must break down somewhere close to the ground-state.

Random interaction models of the kind investigated here describe interacting systems restricted to energy differences of the order of the Thouless energy $`E_c`$ . In disordered diffusive systems one has $`E_c=g\mathrm{\Delta }`$ ($`g`$ is the conductance) and introducing a cut-off at $`E_c`$ leads to $`m`$, $`ng`$. From (4) one would expect the breakdown of perturbation theory to occur for $`U_c\mathrm{\Delta }/g`$, precisely at the amplitude of fluctuations of IME in these systems . The studied model neglects however two ingredients which must be considered when dealing with more realistic systems. First, correlations between the two terms in (5) significantly reduce the value of this latter expression below the uncorrelated estimate $`U^2(\mathrm{\Delta }/g)^2`$. Second, the mean-field charge-charge interaction (that we voluntarily neglected), when incorporated self-consistently into the one-body spectrum, leads to a strong reduction of the one-body density of states , hence an increase of the level spacing at the Fermi level against which (5) must be weighted, and an increase of the denominators for the infrared contributions to second order perturbation theory. It is expectable that these two facts ensure the convergence of perturbation theory and the stability of the ground-state against off-diagonal fluctuations in diffusive metallic systems.

In summary, based on second order perturbation theory we predicted a fast reordering of the ground-state for systems of randomly interacting fermions which we confirmed numerically. While the influence of this reordering on the quasiparticle excitations is at present unclear, we found that it is accompanied by a change of statistics of the first energy excitation above the ground-state. Similar investigations in lattice models where this onset of level repulsion may influence the system’s transport properties , and for which correlations between IME should play an important role (neglected in the presently studied model) are highly desirable.

Numerical computations were performed at the Swiss Center for Scientific Computing. Work supported by the Swiss National Science Foundation.