text:
# Nonequilibrium orientational patterns in two-component Langmuir monolayers

## I Introduction

Spurred by experimental and technological developments in supramolecular chemistry and biophysics, there is a raising interest to study nonequilibrium structures related to self-organization phenomena in weakly condensed systems. In such systems, attractive potential interactions between constituent molecules are essential. Equilibrium structures in soft matter correspond to minima of free energy, and are caused by the competition between short-range attractive and long-range repulsive interactions . In contrast to this, nonequilibrium patterns require permanent supply of energy and correspond to persistent (stationary or time-dependent) kinetic states of a system. Typically, they result from an interplay between reactions, diffusion and phase transitions. Examples of nonequilibrium pattern formation in reactive soft matter include stationary Turing-like patterns in phase-separating binary polymer mixtures with chemical reactions  and in monomolecular adsorbates on metal surfaces . Theoretical analysis for two-component reactive adsorbates  and for three-component reactive polymer systems  has further shown that not only stationary patterns, but also traveling and standing waves are possible.

Langmuir films are monomolecular layers of amphiphilic molecules on an air-water interface. Such films are a classical example of soft matter and, generally, it should be expected that, in the presence of chemical reactions and energy flows, they would form nonequilibrium patterns. Though equilibrium properties of Langmuir monolayers are thoroughly investigated (see ), nonequilibrium phenomena in these systems still remain less explored. In an interesting series of experiments, Tabe and Yokoyama have studied illuminated Langmuir monolayers of amphiphilic derivatives of azobenzene by means of the Brewster-angle microscopy . In these experiments, transitions between trans and cis conformations of individual molecules were photo-induced by polarized light of a selected wavelength. Since physical properties of the two conformations are different, the trans and cis isomers essentially represented two different species. In addition to photo-induced periodic stationary patterns, these experiments have shown, for the first time, spontaneously emerging patterns of propagating waves of molecular reorientation under appropriate illumination conditions.

In the Letter , we proposed a model of reactive two-component Langmuir monolayers with orientational ordering. This model was taking into account phase separation in the two-component system, reaction interconverting both species, diffusion of reactants and processes of orientational ordering. Interactions between the components of the monolayer resulted not only from the positional, but also from the orientational order of the hydrophobic tails of constituent molecules, determined by their tilt. We have shown that this model already reproduces nonequilibrium traveling-wave patterns which arise as a consequence of a Hopf bifurcation with a finite wavenumber. To simplify the analysis, it was assumed that the azimuthal orientation of molecules remained fixed and uniform, so that only their tilts could vary. Moreover, some orientational-order contributions to the free energy of the monolayer were neglected. The inclusion of azimuthal variations is, however, important for a comparison with the experimental data yielded by the Brewster-angle microscopy that is sensitive to the local azimuthal orientational ordering.

The aim of the present paper is to formulate and to study a more general model that contains both orientational variables, and includes bend and splay distortion terms. Our analytical and numerical investigations show that this model has a significantly different phase diagram and new kinds of nonequilibrium patterns are possible here. After introduction of the model in Sec. II, we investigate in Sec. III the behaviour of the system in the equilibrium case, when the reaction is absent. The bifurcation analysis of the uniform steady state of the nonequilibrium system under illumination is presented in Sec. IV. Numerical simulations, revealing the formation of such nonequilibrium patterns as traveling droplets and stripes, stationary splay-defects, and complex azimuth reorientation kinetics, are reported in Sec. V. The paper ends with the conclusion and the discussion of the obtained results.

## II The Model

We consider a monolayer formed by two diffusive immiscible components, $`A`$ and $`B`$, that have strongly different shapes. Modeling the situation in the photo-isomerization experiments, molecules $`A`$ are supposed to have an elongated shape with a long tail (the trans isomer) and molecules $`B`$ to have a crumpled conformation (the cis isomer). According to this assumption, only molecules $`A`$ are subject to orientational order, whereas molecules $`B`$ play essentially the role of passive dilution with respect to such ordering. Furthermore, a photo-induced reaction interconverting $`A`$ and $`B`$ molecules is considered. The total concentration of components $`A`$ and $`B`$ in the monolayer is assumed constant. Therefore, the local composition of the monolayer is characterized by the concentration order parameter $`c`$ representing the local fraction of molecules $`A`$ (so that $`1c`$ gives the local fraction of molecules $`B`$). The local orientational order is described by the order parameter $`\stackrel{}{a}`$, that corresponds to the projection of the local mesoscopic average of the unit molecular director $`\stackrel{}{n}`$ of the elongated molecules $`A`$ onto the monolayer plane (see Fig. 1). The vector $`\stackrel{}{a}`$ is defined by its modulus $`a=|\stackrel{}{a}|=\mathrm{sin}\eta `$, where $`\eta `$ is the tilt angle, and its azimuth angle $`\phi `$.

After we have introduced the two order parameters $`c`$ and $`\stackrel{}{a}`$ needed to describe the system, the kinetic equations for their evolution should be formulated. Following the mesoscopic approach , the evolution equation for $`c`$ reads as

$$\frac{c}{t}=D^2c+\frac{D}{k_BT}\left[c(1c)\mu \right]+(c).$$
(1)
Here $`D`$ is the diffusion coefficient, $`T`$ is temperature, and $`\mu `$ is the chemical potential defined as $`\mu =\delta /\delta c`$, where $``$ is the free energy functional that will be specified below. In the presence of nonpolarized light, the reactive term $`(c)`$ is given by ,

$$(c)=\left[Ig(\lambda )+k_1^0\right](1c)\left[If(\lambda )+k_1^0\right]c,$$
(2)
where the first term corresponds to the conversion of molecules $`B`$ into molecules $`A`$ and the second term describes the reverse conversion process. In this equation, $`I`$ is the light intensity, $`\lambda `$ is the wavelength of light, and $`f(\lambda )`$ and $`g(\lambda )`$ are related to the surface molar extinction coefficients and quantum yields for the conformations of $`A`$ and $`B`$, respectively. For moderate and high illumination intensities, photo-induced conversion dominates over the thermal conversion processes, so that the thermal rate constants $`k_{\pm 1}^0`$ can be neglected. In this case, the ratio of conversion rates for forward and backward processes is independent of temperature or light intensity, and determined only by the light wavelength $`\lambda `$. For simplicity, we shall assume in our subsequent analysis that the wavelength $`\lambda `$ is chosen in such a way that $`f(\lambda )=g(\lambda )`$ . Under this condition, equation (2) takes the form

$$(c)=k(12c),$$
(3)
where $`k=If(\lambda )=Ig(\lambda )`$ is the reaction constant proportional to the intensity of the applied light.

The kinetic equation for the local orientation vector $`\stackrel{}{a}`$ is chosen as

$$\frac{\stackrel{}{a}}{t}=\mathrm{\Gamma }\frac{\delta }{\delta \stackrel{}{a}}kq(c)\stackrel{}{a}.$$
(4)
The first term on the right side corresponds to relaxation with a typical relaxation time $`\mathrm{\Gamma }^1`$. The second term takes into account that the reaction, interconverting one molecular conformation to another, also affects the local average orientation of molecules.

The choice of the function $`q(c)`$ should be based on the kinetic analysis with respect to the orientation of molecules involved in the reactive process. In our model, the reaction $`BA`$, transforming crumpled molecules into the elongated form, is assumed to be strongly energetically activated by light. This means that new molecules $`A`$ would be created with high initial energy and can adopt any orientation $`\stackrel{}{a}`$, independent of the orientation of surrounding molecules. Therefore, the orientation of newly created molecules $`A`$ is random and the statistical average of $`\stackrel{}{a}`$ is zero. Because any conversion from an elongated molecule $`A`$ to the crumpled molecule $`B`$ decreases the local order and the reverse conversion process does not increase it, the overall reaction effect is that it tends to destroy orientational ordering. Under such assumptions, local evolution for the orientation momentum $`c\stackrel{}{a}`$ is described by a simple equation $`_t(c\stackrel{}{a})=kc\stackrel{}{a}`$, that corresponds to the loss of $`c\stackrel{}{a}`$ when molecules $`A`$ transform into $`B`$, having no defined orientation. Splitting $`_t(c\stackrel{}{a})`$ and substituting the local variation of $`c`$ due to reaction, $`_tc=k(12c)`$, leads to $`_t\stackrel{}{a}=k\stackrel{}{a}(1c)/c`$. Therefore, we have $`q(c)=(1c)/c`$. This function $`q(c)`$ is used below.

The energy functional $``$ is constructed in terms of the order parameters $`c`$ and $`\stackrel{}{a}`$ of the system. We decompose the energy functional into two parts: one that accounts for the affinity between the two isomers in the monolayer, and the other due to the tail-tail orientational interactions. The first of these contributions, $`_c`$, depends exclusively on $`c`$, whereas the tail-tail interaction, $`_{or}`$, is sensitive to both $`c`$ and the orientation of the elongated molecules. Accordingly, $`=_c(c)+_{or}(c,\stackrel{}{a})`$. Note that entropic contributions are not considered, since they have been directly included as the diffusive term in the kinetic equation (1). The expression for $`_c`$ reads as ,

$$_c=𝑑x𝑑y\left[\stackrel{~}{\chi _0}c^2+\frac{1}{2}\stackrel{~}{\chi _2}\left(c\right)^2\right],$$
(5)
and describes phase separation due to lateral interactions between molecules. In the limit of short-range lateral interaction potentials, $`\stackrel{~}{\chi _0}`$ is determined by its strength, and $`\stackrel{~}{\chi _2}`$ can be estimated as $`\stackrel{~}{\chi _2}\frac{1}{2}\stackrel{~}{\chi _0}r_0^2`$, where $`r_0`$ is the characteristic radius of the interaction . Near the critical point of the equilibrium phase separation, the combination of Eqs. (1) and (5) is equivalent to the usual Cahn-Hilliard equation with the Landau free energy functional .

The part of the free energy functional associated with the distortion of the orientation of tilted molecules in the monolayer can be written as

$$_{or}=𝑑x𝑑y\left[\frac{1}{2}\stackrel{~}{p}(c)a^2+\frac{1}{4}\stackrel{~}{\beta }a^4+\frac{\stackrel{~}{K_s}}{2}\left(\stackrel{}{}\stackrel{}{a}\right)^2+\frac{\stackrel{~}{K_b}}{2}\left(\stackrel{}{}\times \stackrel{}{a}\right)^2+\stackrel{~}{\mathrm{\Lambda }}c\left(\stackrel{}{}\stackrel{}{a}\right)\right].$$
(6)
The first two terms correspond to the Landau expansion up to the quartic term for the modulus of the orientation vector . This expansion is only justified for sufficiently small $`\eta `$, since in this case $`a=\mathrm{sin}\eta \eta `$ is small as well. Generally, all coefficients in the Landau decomposition should depend on the local concentration $`c`$. We shall, however, consider only weakly nonuniform states, where local deviations of the concentration $`c`$ from the uniform stationary state $`\overline{c}=1/2`$ are small. Therefore, we neglect the dependence of the coefficient $`\stackrel{~}{\beta }`$ on the variable $`c`$. However, the coefficient $`\stackrel{~}{p}`$ of the quadratic term in the Landau free energy is already small near the instability and its dependence on $`c`$ must be retained. For condensed phases, lowering the lateral pressure of a Langmuir monolayer leads to an increase of its equilibrium tilt . Since we assumed that molecules $`B`$ play a role of passive dilution for the tilted molecules $`A`$, decreasing $`c`$ is roughly equivalent to decreasing the lateral pressure. Thus, we choose a linear dependence $`\stackrel{~}{p}(c)=\stackrel{~}{\pi }_0+2\stackrel{~}{\alpha }(0.5c)`$, where $`\stackrel{~}{\alpha }`$ is a positive coefficient, and $`\stackrel{~}{\pi }_0`$ is a decreasing function of the lateral pressure which determines the equilibrium tilt in the monolayer with $`c=\overline{c}`$.

The other contributions in Eq. (6) stand for the bend and splay distortion terms . More specifically, the third and fourth terms correspond to the classical Frank elasticity terms, that account for the splay and bend distortions, respectively. Normally, one takes the single Frank constant approximation ($`\stackrel{~}{K}=\stackrel{~}{K_s}=\stackrel{~}{K_b}`$). The fifth term is the lowest-order splay contribution, $`(\stackrel{}{}\stackrel{}{a})`$, that appears coupled to the composition order parameter $`c`$. Although some authors  prefer to couple the linear splay term to the tilt angle, we follow the suggestion by Selinger et al. for a two-component monolayer undergoing phase separation, as it is in our case. A similar approach was taken by Tabe et al.  who coupled the linear splay term with a certain order parameter related to the molecular density that varies across the monolayer. In general, for sufficiently strong coupling $`\stackrel{~}{\mathrm{\Lambda }}`$, the linear splay term destabilizes uniformly oriented phases, leading to equilibrium nonhomogeneous splayed states (see Sec. III). Finally, notice that a term $`(\stackrel{}{}\times \stackrel{}{a})`$ linear with respect to bend distortion is not considered, because it is not permitted by symmetry in nonchiral Langmuir monolayers.

Summarizing, the model presented here can be viewed as a Cahn-Hilliard equation for the composition variable $`c`$, coupled to a relaxational equation for the orientation order parameter $`\stackrel{}{a}`$. The contributions to the free energy have been derived considering that $`c`$ is close to its stationary uniform solution and for sufficiently small tilt angles $`\eta `$.

The analysis can be simplified by appropriately adimensionalizing energy, time and space. Energy is measured in units of $`k_BT`$, time in units of the relaxational time $`(\mathrm{\Gamma }k_BT)^1`$, and spatial coordinates are rescaled with the relaxational length $`\sqrt{D/(\mathrm{\Gamma }k_BT)}`$. The model is then characterized by the dimensionless parameters $`\chi _0=\stackrel{~}{\chi _0}/k_BT`$, $`\chi _2=\stackrel{~}{\chi _2}\mathrm{\Gamma }/D`$, $`\pi _0=\stackrel{~}{\pi }_0/k_BT`$, $`\alpha =\stackrel{~}{\alpha }/k_BT`$, $`\beta =\stackrel{~}{\beta }/k_BT`$, $`\kappa =k(\mathrm{\Gamma }k_BT)^1`$, $`K=\stackrel{~}{K}\mathrm{\Gamma }/D`$ and $`\mathrm{\Lambda }=\stackrel{~}{\mathrm{\Lambda }}\sqrt{\mathrm{\Gamma }(Dk_BT)^1}`$. With this choice, the final equations for the evolution of $`c`$ and the two components of vector $`\stackrel{}{a}`$ read as

$`{\displaystyle \frac{c}{t}}`$ $`=`$ $`^2c+\left[c(1c)\mu \right]+\kappa (12c),`$ (7)
$`{\displaystyle \frac{\stackrel{}{a}}{t}}`$ $`=`$ $`p(c)\stackrel{}{a}\beta a^2\stackrel{}{a}+K^2\stackrel{}{a}+\mathrm{\Lambda }\stackrel{}{}c\kappa {\displaystyle \frac{1c}{c}}\stackrel{}{a},`$ (8)
where

$`\mu =2\chi _0c\chi _2^2c+\alpha a^2+\mathrm{\Lambda }\stackrel{}{}\stackrel{}{a},`$

and $`p(c)=\stackrel{~}{p}(c)(k_BT)^1=\pi _0+2\alpha (0.5c)`$. This is the mathematical model which will be investigated below.

In order to obtain numerical results from the model, we numerically integrate Eqs. (8) on a $`100\times 100`$ square grid, using an explicit Euler scheme with a mesh size $`\mathrm{\Delta }x`$ and a time step $`\mathrm{\Delta }t`$ that assure a good numerical accuracy. Periodic boundary conditions are chosen to model the behaviour in a large system far from the boundaries. As initial conditions, small random perturbations around the stationary uniform states of the system $`\overline{c}`$ and $`\overline{a}`$, and a random distribution of azimuth angles $`\phi (0,2\pi )`$ are taken. To display simulation results, snapshots of the patterns after the transients are given in the figures. Each figure consists of two panels: the left panel shows in grey scale the spatial distribution of the composition variable $`c`$, with larger values corresponding to the darker color, and the right panel is used for visualization of the orientational field $`\stackrel{}{a}`$. The local directions of this field are visualized by using small arrows. Note that, for technical reasons, such arrows could not be used to indicate the states of all grid points in the simulations and therefore the visualization of the azimuthal orientation is rough. The grey color in the right panels is used to display the local tilt $`a`$, and again, darker regions correspond to the larger values of this variable. In some cases, videos of time-dependent patterns are additionally provided.

## III Equilibrium Patterns

Before addressing the nonequilibrium cases, we show some examples of equilibrium pattern formation in the considered system. The equilibrium conditions correspond to absence of illumination and are realized if all reactive terms in Eqs. (8) are omitted. One of the limitations of our previous simpler model  was that it did not exhibit any equilibrium pattern formation, despite the experimental evidence of spontaneous generation of striped patterns in non-illuminated monolayers . The linear splay term included in Eq. (6) can already lead to the formation of equilibrium orientational structures.

In the limit of $`\alpha 0`$ (i.e., in the absence of the quadratic Landau term for the tilt variation), the nonreactive version of the present model is similar to the description proposed by Selinger et al.  for nonchiral Langmuir monolayers. In their study, three different nonuniform phases (sinusoidal stripes, soliton stripes and square lattice of vortices) were found for smectic films undergoing chiral symmetry breaking under variation of a control parameter (corresponding to temperature), and similar phases for multicomponent Langmuir monolayers were predicted. As illustrated in Figs. 2, 3, 4 and 5, our model reproduces these regimes under appropriate choices of the parameter $`\chi _0`$ specifying the relative intensity of energetic lateral interactions.

When the characteristic energy of lateral interaction is much weaker than the orientation energy, a square array of alternating “inward” and “outward” splay defects (vortices) is formed (Fig. 2). The “inward” defects represent regions rich in elongated molecules $`A`$ which are oriented towards the center of the defect. In analogy to the “escape to the third dimension” found in defects in $`3D`$ nematics , the tilt in the center of these defects is almost zero, reducing the Frank elastic energy near that point. The “outward” splay defects are poor in elongated molecules $`A`$, but also exhibit vanishing tilt in their centers. Increasing the lateral interactions with respect to the splay coupling, the system organizes into a stripe phase with smooth variations of $`c`$ and $`a`$ (see Fig. 3). The profiles of variation of c, $`a`$ and $`\phi `$ across a stipe are displayed in Fig. 4. The equilibrium stripe patterns are similar to those observed in the experiments in absence of illumination . The modulation of the tilt (absent in the model by Selinger et al. because tilt variations were not included there) has a spatial frequency twice that of the azimuthal and concentration modulations. Moreover, in agreement with the experiments the amplitude of the tilt angle modulation depends on the stripe size, decreasing as the stripe widens (this is observed, for example, by decreasing the parameter $`\mathrm{\Lambda }`$).

Strong lateral attractive interactions (i.e., large interaction strengths $`\chi _0`$) lead to the formation of a striped phase with sharp wall domains (”soliton stripes”) in the modulation of composition, tilt and azimuth, which has also been reported by Selinger et al. . In Fig. 5 profiles of variation of c, $`a`$ and $`\phi `$ across such a stripe are plotted. Note again the double frequency modulation of the tilt with respect to the composition and azimuth variation.

We have checked that the inclusion of the quadratic Landau term for the tilt variable, neglected in the simulations described in this section, does not lead to significative differences in the properties of equilibrium patterns.

## IV Nonequilibrium Phenomena: Linear Stability Analysis

Full numerical exploration of the parameter space of the model is difficult because of the large number of relevant parameters. The linear stability analysis of the uniform stationary states can provide indications of what regions of the parameter space are worth being examined in search for a particular kind of a pattern. The stationary states of the system (8) are $`c=\overline{c}=1/2`$ and $`a=\overline{a}=\sqrt{\left(\pi _0\kappa \right)/\beta }`$, provided that $`\pi _0>\kappa `$. The azimuthal orientation is arbitrary in such a state. For convenience, we choose $`\overline{a_x}=\overline{a}`$ and $`\overline{a_y}=0`$. The linear stability of these uniform solutions is performed by adding small plane-wave perturbations $`\delta c`$, $`\delta a_x`$ and $`\delta a_y`$ proportional to $`\mathrm{exp}(iq_xx+iq_yy+\gamma (\stackrel{}{q})t)`$, and linearizing Eqs. (8). This leads to the linearization matrix

$$=\left(\begin{array}{ccc}\frac{q^2}{2}\left(\chi _0+\frac{q^2}{2}\chi _2+2\right)2\kappa & \frac{q^2}{2}\left(\alpha \overline{a}+i\frac{\mathrm{\Lambda }}{2}q_x\right)& i\frac{\mathrm{\Lambda }}{4}q_yq^2\\ 2\overline{a}(2\kappa \alpha )+iq_x\mathrm{\Lambda }& \pi _0\kappa 3\beta \overline{a}^2Kq^2& 0\\ iq_y\mathrm{\Lambda }& 0& Kq^2\end{array}\right),$$
(9)
where $`\stackrel{}{q}=(q_x,q_y)`$ is the wavevector of the considered mode. The first line in the matrix corresponds to the composition variable $`\delta c`$ and the next two lines stand for the orientation components $`\delta a_x`$ and $`\delta a_y`$.

The (complex) rates $`\gamma (\stackrel{}{q})`$ of various modes are determined by the eigenvalues of the linearization matrix $``$. The unstable modes are identified like those with Re$`(\gamma (\stackrel{}{q}))>0`$. If the imaginary part of $`\gamma (\stackrel{}{q})`$ is not zero for the first unstable mode, we have a wave instability (a Hopf bifurcation with a finite wavenumber), resulting in traveling or standing waves (cf. ). On the other hand, if Im$`(\gamma (\stackrel{}{q}))=0`$ for the first unstable mode, a Turing instability leading to nonequilibrium stationary periodic patterns is realized. The values of $`\gamma (\stackrel{}{q})`$ are yielded as the roots of the characteristic equation associated with the matrix $``$. This characteristic equation is, however, cubic and therefore its analytical solution is possible only in some special cases.

If splay coupling is absent ($`\mathrm{\Lambda }=0`$), the stability analysis is simplified and becomes equivalent to that of the previously studied reduced model . The phase diagram for the non-splay case in the plane ($`\pi _0`$,$`\kappa `$) is presented in Fig. 6a. We see five different regions, whose boundaries and marginal wave numbers can be obtained analytically . Region $`I`$ corresponds to the wave instability regime, where traveling or standing waves are expected. In region $`II`$ we have the Turing instability region, where stationary droplet-like structures with periodic variation of both local concentration and tilt are found. Region $`III`$ corresponds to the uniform tilted phase. When $`\pi _0<\kappa `$, only nontilted phases ($`\overline{a}=0`$) exist: a nonuniform phase due to a Turing instability in region $`IVa`$ and stable uniform phase in region $`IVb`$.

In order to understand the effects of the linear splay term, that was neglected above, we perform the stability analysis of the full model equations for a fixed wave plane $`q_y=0`$. In this case, the variable $`a_y`$ is decoupled from the other two degrees of freedom (actually, its dynamics is exclusively governed by the elastic damping term, see Eq. (9)), so that the stability analysis is reduced to a $`2\times 2`$ problem. The following results, summarized in Fig. 6b, are obtained: For nontilted phases ($`\pi _0<\kappa `$), the effect of the linear splay term is to move to larger reaction rates $`\kappa `$ the boundary between the Turing-instability region ($`IVa`$) and the uniform nontilted phase ($`IVb`$). For tilted phases ($`\pi _0>\kappa `$), the changes are more profound. All unstable modes in this area of the parameter space have now a nonvanishing imaginary part Im$`(\gamma (\stackrel{}{q}))`$. Therefore, the region $`II`$ with Turing instability and stationary tilted patterns completely disappears and becomes replaced by the region $`I`$ with traveling waves. Additionally, the region $`III`$ occupied by the tilted uniform phase is reduced.

Based on this (limited) stability analysis, one can expect that the traveling patterns would be found more often when the splay interactions are taken into account. Furthermore, it can be expected that the traveling patterns in the parameter region occupied by stationary Turing-like structures in absence of splay interactions (region $`II`$ in Fig. 6a), would be different for weak splay interactions from the traveling patterns in the parameter region where traveling waves are observed even without the splay effects (region $`I`$ in Fig. 6a). As $`\mathrm{\Lambda }0`$, traveling waves in the former region $`II`$ should slow down and become frozen at $`\mathrm{\Lambda }=0`$. Indeed, the analytical stability investigation shows that, in this region, the velocity of the most unstable mode is proportional to $`\kappa \mathrm{\Lambda }\overline{a}`$. On the other hand, the velocity of patterns in the former region $`I`$ remains finite in the limit $`\mathrm{\Lambda }0`$.

Numerical simulations of the model, which represent the main part of the reported study, agree with the predictions of the stability anaylsis.

## V Nonequilibrium Phenomena: Numerical Results

To facilitate the comparison with previous simulations of the model without azimuthal variation, we choose here the same numerical values of the common parameters $`\chi _2=0.0052`$, $`\chi _0=2.27`$, $`\alpha =1.5`$ and $`\beta =2`$. The above mentioned changes in the phase diagram due to the inclusion of the linear splay term indicate that this parameter region is worth examining with the model presented here in order to obtain spatio-temporal behaviors involving now composition, tilt and azimuth modulations. We examine numerically such parameter region, and the results are summarized in the three following subsections according to the value fixed for the splay coupling constant $`\mathrm{\Lambda }`$. All simulations are perfomed for a system with a linear size $`L=5`$.

### A No splay coupling

When splay coupling is absent ($`\mathrm{\Lambda }=0`$), azimuthal orientation of elongated molecules is not influenced by variations of the local composition $`c`$ or the local tilt $`a`$. Indeed, by looking at the second of equations (8) we notice that then all terms on the right hand side of this equation, except for the elastic term $`K^2\stackrel{}{a}`$, are proportional to $`\stackrel{}{a}`$ and therefore cannot change the direction of this vector. If the azimuthal orientation is initially uniform ($`\phi =const`$), this state is maintained at all times. For such a state, Eqs. (8) reduce to the model which was already investigated in Ref. and new numerical simulations are not needed. If the initial azimuthal orientation is not uniform ($`\phi const`$), the subsequent evolution of the orientation field is determined only by the elastic term.

In numerical simulations, we choose the parameters inside the region $`I`$ with traveling waves in the phase diagram shown in Fig. 6a, by fixing $`\pi _0=1.5`$, $`\kappa =1`$ and taking a small value $`K=0.001`$ of the elastic interaction constant. Random distribution of azimuthal orientations is chosen as the initial condition. The simulation results are presented in Fig. 7 (see also video Fig7.mpg). Because of elastic interactions, the molecules tend to have parallel orientation, and this leads, after some time, to a pattern characterized by a number of orientational defects that remains stationary henceforth. In the center of a defect, the molecules are non-tilted; the azimuthal direction changes by $`2\pi `$ after passing around a defect. The waves, similar to those described in Ref., travel on the background of the stationary orientational pattern, and their motion is not generally influenced by the azimuth angle distribution. The waves, however, break when they pass through a defect.

Inside the region $`II`$ in the phase diagram in Fig. 6a, stationary Turing-like patterns representing arrays of droplets are observed.

### B Weak splay coupling

The most important change in the phase diagram due to the inclusion of the linear splay term is that the region $`II`$ occupied by the Turing-like patterns in absence of splay coupling (Fig. 6a) is transformed into region $`I`$ where traveling patterns should be observed (Fig. 6b). However, as we shall see, the properties of traveling patterns in the parameter domain are different, for weak coupling, from those of the traveling waves in the old region $`I`$.

Figure 8 and video Fig8.mpg show the traveling pattern yielded by numerical simulation with the parameters $`\kappa =0.5`$, $`\pi _0=1.25`$, and $`\mathrm{\Lambda }=K=0.01`$ that correspond to the former region $`II`$. The pattern looks like a flow of droplets, with the local direction of flow determined by the azimuthal orientation of the elongated molecules in the monolayer. Such droplet-like structures of large $`c`$ and small $`a`$ emerge and move following the local molecular orientation path. Occasionally, rupture of the droplets, when they happen to approach divergent points of the molecular orientation field, is observed. Some of the droplets are pinned by the orientational defects and exhibit only rotation, but not a translational motion. Remarkably, the pattern of azimuthal orientation remains frozen under weak splay coupling, as it was the case in its absence.

In agreement with the linear stability analysis, the velocity of the emerging droplet structures is an increasing function of the splay coupling constant $`\mathrm{\Lambda }`$ and the reaction rate constant $`\kappa `$. Moreover, we have found that the size of such traveling structures is affected by the strength of elastic interactions: the larger $`K`$ the bigger are the droplets. For large $`K`$, however, the uniform state becomes stable, and no spatial organization is observed.

Figure 9 and video Fig9.mpg show traveling waves for a larger value of the reaction rate constant ($`\kappa =1`$), such that we are now inside the old region $`I`$ in Fig. 6a. The morphology of the pattern is now different and it resembles the pattern of traveling waves in absence of splay coupling (Fig. 7). The pattern is formed by stripes that move along the direction determined by local azimuthal orientation of the elongated molecules. Two orientational defects are seen in the upper right corner in Fig. 9,and they also correspond to defects of the traveling stripe pattern of the composition field. The waves are rotating around these orientational defects. Again, the pattern of azimuthal orientation becomes frozen after a rapid initial transient and it is not significantly affected by the traveling waves. Note from the simulation videos that the stripes move much faster than the droplets in Fig. 8.

Next, we examine more closely the profiles of the composition $`c`$ and the tilt $`a`$ in different traveling patterns. Figure 10a shows such profiles for a single traveling droplet from Fig. 8, also displayed in the inset in the left upper corner of this figure. The profiles in the cross sections which are parallel and perpendicular to the motion direction are presented here. The droplet corresponds to a local increase in the concentration of elongated molecules $`A`$ and a local decrease in the tilt of such molecules. It can be noticed that the tilt is also slightly increased along a ring surrounding the droplet, and that the droplet is not axially symmetric. Comparing the profiles for the perpendicular and parallel cross sections in Fig. 10a, we find that the droplet is slightly elongated in the direction parallel to the azimuth molecular orientation and, moreover, the tilt peak in the rear part of the moving droplet is higher than that at the front. This asymmetry determines the propagation direction of the droplet. Furthermore, we have checked that the asymmetry gets stronger when the splay coupling coefficient $`\mathrm{\Lambda }`$ is increased.

The profiles of a traveling stripe from Fig. 9 are displayed in Fig. 10b. We see that the profiles exhibit more smooth variation in this case, and the variation of $`c`$ and $`a`$ is closer to harmonical. Note that the spatial profile of the tilt $`a`$ is retarded with respect to that of the local composition $`c`$ and this again determines the motion direction.

In the simulations described above random initial conditions were chosen. We have also studied some special initial conditions for the orientation distribution that give rise to domain sinks or sources. In the simulation shown in Fig. 11 and video Fig11.mpg, the initial condition with two orientational domains ($`\phi =0`$ for $`x<L/2`$ and $`\phi =\pi `$ for $`x>L/2`$) has been chosen. The boundary between the domains plays here a role of source emitting waves. For other initial conditions, corresponding to outward vortices of molecular orientation, we have observed droplets that are generated in the center of a vortex and travel in the radial direction out of it. Generally, linear and point defects with positive (negative) splay $`(\stackrel{}{}\stackrel{}{a})`$ act as sources (sinks) for the traveling structures that involve tilt and composition variations.

### C Strong splay coupling

When the coefficient $`\mathrm{\Lambda }`$, specifying the intensity of splay coupling, is further increased, the azimuthal orientation of molecules becomes influenced by the traveling or stationary patterns. For strong splay coupling, stationary nonequilibrium patterns are usually observed. An example of such a pattern is shown in Fig. 12. Starting from random initial conditions, the system first develops a pattern of traveling stripes. Subsequently, the stripes undergo breakdown and a frozen array of orientation defects is produced. Inside each defect, the concentration of the elongated molecules is increased, and in general, these molecules are oriented towards the center of a defect. This stationary pattern is favoured in the case of strong coupling, since it minimizes the splay contribution towards the free energy of the system.

To avoid the formation of immobile splay defects, a region in the parameter space, characterized by the wave instability but lying closer to the boundary of a transition to the nontilted state, can be considered. Because the elongated molecules are only slightly tilted in this region, their reorientation is energetically easier than far from the boundary $`\overline{a}=0`$. If a relatively large value of the reaction rate constant $`\kappa `$ is additionally chosen, so that the system is driven further away from thermal equilibrium, traveling waves that are accompanied by azimuthal reorientation of molecules are formed.

An example of the developing wave pattern is shown in Fig. 13 and video Fig13.mpg. The temporal evolution of this pattern along the central vertical cross section is additionally displayed in Fig. 14. After a short transient starting from random initial conditions, the system soon developes a regular pattern of stripes traveling at a constant velocity.

The profiles of the composition $`c`$, the tilt $`a`$ and the azimuth angle $`\phi `$ along the line orthogonal to the propagation direction are shown in Fig. 15. An important difference with the profiles of equilibrium stationary stripes is that, for traveling stripes, the azimuth angle does not undergo $`2\pi `$ rotations from one stripe to another (seen in Fig. 4). Instead, a periodic angular modulation within the interval from zero to $`\pi `$ is observed here. Moreover, it can be noticed that the composition $`c`$ and the tilt $`a`$ vary completely in phase for such traveling stripes, i.e. the maximum of the concentration field coincides with a minimum of the tilt distribution. This differs from the behaviour of such variables in traveling stripes at weak splay coupling (cf. Fig. 10b). The propagation direction is now fixed by a small shift observed in the azimuth angle variation with respect to the profiles of variables $`c`$ and $`a`$.

## VI Conclusions

We have formulated a theoretical model to describe pattern formation in illuminated two-component Langmuir monolayers. In this model, nonequilibrium pattern formation results from an interplay between phase separation, optically induced transitions between the two immiscible conformations, and coupling of phase separation to orientational ordering of molecules in the monolayer. In comparison with the previous model , the current description is more complete because variations of the azimuthal variation of the elongated molecules in the monolayer are allowed. We have also incorporated into the model the splay coupling between azimuthal orientation and variations of local concentration. The situation corresponding to illumination with non-polarized light has been considered.

In contrast to the earlier simple model , the investigated system is characterized, in absence of illumination, by the formation of equilibrium stationary patterns representing arrays of orientational defects (vortices), stripes or soliton-like structures with sharp domain walls. All these patterns are caused by splay interactions and have a purely energetic origin, so that they can be also interpreted by considering minima of free energy (cf. ).

Linear stability investigation of the uniform state of the model under illumination conditions has been performed. It reveals that, generally, splay interaction favour the appearance of traveling patterns. Such patterns are observed even in the parameter regions occupied by stationary structures in the earlier model .

Numerical simulations of the proposed model have been undertaken. They show a rich spectrum of spatio-temporal structures. For weak splay coupling and relatively low reaction rates, a pattern of droplets slowly flowing along the local directions determined by azimuthal orientation of molecules has been observed. If the reaction is stronger, the droplets are replaced by a pattern of traveling stripes following the local azimuthal orientation. Interactions between droplets and stripes and the orientational defects have been seen. It is found that linear and point orientational defects play a role of sources or sinks of traveling structures in this system.

When splay coupling is weak, the spatial pattern of azimuthal orientation is not affected by local concentration variations. The temporal evolution of the orientation pattern in this case is governed by the elastic interactions between molecules. Starting from a random initial distribution, they produce after a short transient a stationary pattern of azimuthal orientation that remains frozen afterwards. In contrast to this, strong splay coupling with the concentration field leads to the appearance of stationary and spatiotemporal patterns where variations of azimuthal orientation are correlated with the the changes in the local composition of the monolayer. Often, stationary arrays of splay defects are observed in this case. However, when the system is close to the orientational phase transition (from the tilted to the non-tilted phase), patterns of traveling waves accompanied by strong variation of azimuthal orientation are found in the model.

This investigation, as well as the earlier Letter , are motivated by the experimental discovery of traveling structures in illuminated Langmuir monolayers by Tabe and Yokoyama . The extended model, which we have now explored, is able to reproduce some essential features of the experiments. In absence of illumination, equilibrium orientation patterns are yielded by the model. Under illumination for sufficiently strong splay coupling, traveling waves accompanied by changes in azimuthal orientation of molecules are also found. The experiments  were performed by using polarized light, and light polarization significantly influenced the properties of the observed patterns. In contrast to this, our model does not include the effects of light polarization and thus corresponds to illumination with non-polarized light. Though the inclusion of anisotropy effects due to light polarization is straightforward, the resulting model is more complicated and not yet discussed in the present paper.

Langmuir monolayers represent a classical example of soft matter and are closely related to biomembranes, playing a fundamental role in cell biology. Therefore, investigations of nonequilibrium pattern formation in such monolayers under reactive conditions can help in understanding of general mechanisms of nonequilibrium self-organization in soft matter. Our study provides evidence that traveling orientational wave patterns may represent a generic property of Langmuir monolayers that are subjected to chemical reaction and composition changes. We have found that various patterns, representing traveling droplets or stripes, wave sources and sinks, and orientational defects interacting with the traveling structures, are possible in such systems.

## Acknowledgments

This work was supported by the Dirección General de Investigación (Spain) through Project BXX2000-0638, and by the Comissionat per Universitats i Recerca (Generalitat de Catalunya) through Project 1999SGR00041. The study was also supported through the ESF REACTOR network programme.


format: text

text:
# On the turbulent 𝛼-disks and the intermittent activity in AGN

## 1 Introduction

The standard accretion disk (Shakura & Sunyaev 1973) is known to be subject to the thermal - viscous instability due to the partial hydrogen ionization (Meyer & Meyer-Hofmeister 1981; Smak 1984). As a result of this instability the disk cycles between the two states: a hot and mostly ionized state with a large local accretion rate and a cold, neutral state with a low accretion rate. This instability has originally been proposed to explain the large amplitude luminosity variations observed in cataclysmic variables (CV; Smak 1982; Meyer & Meyer-Hofmeister 1982). It is also believed that the same mechanism is responsible for the eruptions in soft X-ray transients (SXT; see e.g. Cannizzo, Ghosh & Wheeler 1982; Dubus, Hameury & Lasota 2001; Lasota 2001 for review). The ionization instability was also shown to operate in the disk around a supermassive black hole in active galactic nuclei (AGN; Lin & Shields 1986; Clarke 1988; Mineshige & Shields 1990; Siemiginowska, Czerny & Kostyunin 1996).

The characteristic timescales of a cycle activity scale roughly with the mass of a compact object (Hatziminaoglou et al. 2001). Therefore, the observed timescale of cycle of order of years in binaries translates into thousands to millions of years in galaxies, which harbor a supermassive black hole.

The variability amplitudes and timescales observed in Galactic binary systems require the viscosity to be smaller in the quiescent disk than during the outburst. Various ad-hoc viscosity scaling laws have been used in the disk evolution models, including the simplest one with a constant viscosity, however usually $`\alpha `$ is assumed to be 4-5 times larger during the outburst than in the quiescence (Cannizzo 1993). The value of $`\alpha _{\mathrm{hot}}`$ determines the timescale of an outburst, while $`\alpha _{\mathrm{cold}}`$ governs the separation between the subsequent outbursts. The difference between these two parameters determines the outburst amplitude.

The dependence of the viscosity parameter on the disk state is suggested to come out as a property of the magneto-rotational instability (MRI) mechanism (Balbus & Hawley 1991; 1998). Nonlinear development of this instability, which is assumed to be a primary source of the disk turbulence, is sensitive to the presence of resistive diffusion of magnetic field (Hawley, Gammie & Balbus 1996), measured by the magnetic Reynolds number $`Re_\mathrm{M}`$. Gammie & Menou (1998) calculated the Reynolds number for the two states of the disk in CVs and show that it is low in the quiescence. It means that the magneto-hydrodynamical (MHD) turbulence dies away and the matter accumulates in the outer disk.

Menou & Quataert (2001) argued that the Reynolds number will not be low enough to suppress the MHD turbulence in the quiescent disk of AGN. Therefore the efficiency of angular momentum transport should be comparable in the hot and cold states ($`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$). In such a case the outburst amplitude is dramatically reduced, and the ionization instability leads only to a small amplitude flickering (Siemiginowska et al. 1996).

The conclusion about the small amplitude of outbursts was obtained, however, for a geometrically thin, optically thick disk extending down to the marginally stable orbit. This assumption may not be correct. When the local accretion rate is low, the local radiative cooling may not be efficient (Rees et al. 1982; Begelman 1985). The innermost part of the disk may be replaced with an optically thin, hot and possibly two-temperature plasma. An example of such a solution, advection dominated accretion flow (ADAF) was introduced by Ichimaru (1977) and Narayan & Yi (1994). The development of the inner flow through the disk evaporation was recently discussed by several authors, starting with Meyer & Meyer-Hofmeister (1994).

In this paper we first perform a self-consistency check of the $`\alpha `$ description of the angular momentum transfer. Having justified $`\alpha _{\mathrm{cold}}`$ = $`\alpha _{\mathrm{hot}}`$ approach to AGN accretion disks, we consider the time evolution of the disk under the influence of the ionization instability. We also take into account the evaporation of the inner disk.

We define a parameter space where the instability zone (i.e. the partial hydrogen ionization zone) and the evaporation region overlap. The overlapping of these two regions strongly enhances the outburst amplitude due to the time evolution of the disk evaporation radius. We illustrate this behavior with exemplary lightcurves.

The structure of this article is as follows. In Section 2 we define the criterion for efficient development of the MHD turbulence in the accretion disk, taking into account its vertical structure. We confirm that indeed the bi-modal viscosity behavior is characteristic for Galactic binary systems. However, a single value of the viscosity parameter is appropriate in the case of AGN disks. In Section 3 we check if the zone of partial hydrogen ionization overlaps with the disk evaporation region and we calculate the evolution of the disk luminosity in case of the evaporated inner disk. Finally, in Section 4 we discuss our results and give conclusions.

## 2 Magneto-rotational turbulence

In this Section we estimate the strength of the MHD turbulence in the $`\alpha `$ accretion disks and investigate whether the viscosity prescription with a constant parameter $`\alpha `$ is appropriate for both Galactic systems (CVs and SXTs) and AGN. We only consider the instability due to partially ionized hydrogen and do not include regions affected by the radiation pressure instability (see Janiuk, Czerny & Siemiginowska 2002 for that model details). We assume that the viscous stress tensor is proportional to the gas pressure $`P_{\mathrm{gas}}`$ and therefore there is no radiation pressure dominated branch on the stability figures (so-called S-curves) described below.

We use the dimensionless accretion rate in Eddington units, $`\dot{m}=\dot{M}/\dot{M}_{\mathrm{Edd}}`$, assuming the efficiency of 1/12, as implied for the Schwarzschild black hole by the Newtonian potential:

$$\dot{M}_{\mathrm{Edd}}=3.52\frac{M}{10^8M_{}}[M_{}/yr].$$
(1)
The disk vertical structure (indicated by $`z`$) is calculated by solving the equations of viscous energy dissipation, hydrostatic equilibrium, and energy transfer:

$$\frac{dF}{dz}=\alpha P_{\mathrm{gas}}(\frac{d\mathrm{\Omega }}{dr})$$
(2)
$$\frac{1}{\rho }\frac{dP}{dz}=\mathrm{\Omega }^2z$$
(3)
$$\frac{dT}{dz}=\frac{3\kappa \rho }{4acT^3}F_\mathrm{l}$$
(4)
where we solve for temperature, density and pressure profiles, $`T(z)`$, $`\rho (z)`$ and $`P(z)`$. Here $`\mathrm{\Omega }`$ is the Keplerian angular velocity, $`a`$ and $`c`$ are physical constants. $`F_\mathrm{l}`$ is the energy flux transported locally in the direction perpendicular to the equatorial plane and carried either by radiation or by convection:

$$F_\mathrm{l}=F_{\mathrm{rad}}_{\mathrm{rad}}_{\mathrm{ad}},F_\mathrm{l}=F_{\mathrm{rad}}+F_{\mathrm{conv}}_{\mathrm{rad}}>_{\mathrm{ad}}$$
(5)
The frequency-averaged opacity $`\kappa `$ is taken to be the Rosseland mean and includes the electron scattering, free-free and bound-free transitions. The opacity tables are from Alexander, Johnson & Rypma (1983) and Seaton et al. (1994). The presence of dust and molecules is included in the opacity description. The details of the model were discussed in Pojmański (1986) and Różańska et al. (1999).

Figures 1 and 2 show, in $`\mathrm{\Sigma }`$ vs. $`T_{\mathrm{eff}}`$ (surface density vs. effective temperature) plane, the local solutions of the disk vertical structure calculated for a range of accretion rates for two extreme cases of the mass of the central object: 1M and 3 $`\times `$ 109M. The solutions located on both lower (below point A) and upper (above point B) branches of the S-curves are stable. The slope of the middle branch is, at least in some parts, negative which means that at that range of accretion rates the disk is thermally and viscously unstable.

Note, that in case of CVs (Figure 1) the location of the starting point of the instability is very sensitive to the adopted value of $`\alpha `$. For a small viscosity parameter, $`\alpha =0.02`$, the point $`A`$ is at much lower temperature than in case of large $`\alpha =0.1`$. The latter ’smoothens’ the S-curve in its lower part, so that only a single ’wiggle’ remains. Therefore the point $`A`$, defined as the first critical point above which the S-curve slope becomes negative, is shifted upwards in temperature in comparison with that characteristic for small viscosities. This fact should be taken into account when modeling the limit cycle in dwarf novae by means of the ’combined’ S-curves with $`\alpha _{cold}<\alpha _{hot}`$.

This behavior is however characteristic only for accretion disks around 1$`M_{}`$, while vanishes already for SXTs. Also for supermassive black holes there are always two ’wiggles’ in the S-curve, regardless of the viscosity parameter.

The ionization instability is characteristic for all accretion disks but the range of accretion rates (i.e. the location of the S-curve on the $`T\mathrm{\Sigma }`$ plane) depends on the chosen disk radius. We can invert this problem and say that for a fixed external accretion rates the instability will appear only for a certain range of radii. At larger radii the disk is on the lower stable branch while at lower radii it is on the higher, also stable branch. Simple analytical formulae for the unstable zone were provided by Siemiginowska et al. (1996) and applied by Menou & Quataert (2001) in their analysis.

In Figure 3 we show the radial extension of the unstable zone as a function of the accretion rate (in Eddington units) calculated numerically from our disk model. The position of the instability zone depends on the mass of the central object so we choose values representative for all objects from CVs to AGN with a large black hole mass. The plot shows the radial extension as a function of the accretion rate (in Eddington units). The horizontal slices of the shaded region correspond to the accretion rates, for which the disk is unstable at a particular radius. For $`\dot{m}0.01`$ (an accretion rate typical for many objects) the instability zone is located at $`10^5R_{\mathrm{Schw}}`$ for a CV, at $`2\times 10^3R_{\mathrm{Schw}}`$ for $`M10^6M_{}`$ and at $`500R_{\mathrm{Schw}}`$ for an extremely massive black hole of $`M=3\times 10^9M_{}`$.

Time evolution of the unstable part of the disk proceeds roughly in the form of oscillations between the upper and the lower stable branches. As argued by Gammie & Menou (1998), high temperature upper branch conditions are always favorable for the development of the efficient MRI and high values of viscosity are appropriate there. Lower branch conditions are different and the MRI mechanism may not be efficient. Therefore, we perfom our self-consistency check paying particular attention to the lower branch solutions at the vicinity of the turning point $`A`$.

### 2.1 Coupling of the magnetic field with the gas

The behavior of the magnetic field is governed by the fluid conductivity $`\sigma `$. The time-dependent magnetic field in the disk is described by the equation (Parker 1979):

$$\frac{B}{t}=\times (\stackrel{}{v}\times \stackrel{}{B})+\eta ^2\stackrel{}{B}$$
(6)
where $`\eta c^2/4\pi \sigma `$ denotes resistivity. The characteristic diffusion time in which the initial configuration of the magnetic field will decay is equal to $`\tau =L^2/\eta `$, where $`L`$ indicates the characteristic spatial scale. For timescales short in comparison with $`\tau `$ the second term in Equation 6 can be neglected and the magnetic field lines are frozen into the gas. Magnetic Reynolds number, defined as:

$$Re_\mathrm{M}=\frac{v\tau }{L}$$
(7)
can be used to distinguish between two cases: (1) the diffusion of field lines within the disk and (2) the lines are frozen in and carried along with the matter. Here we identify the velocity $`v`$ with the sound speed in the disk and the length $`L`$ of the magnetic field spatial variations with the disk thickness, as commonly used in the simulations (see e.g. Hawley et al. 1996; Gammie & Menou 1998).

Apart from the ohmic diffusion also the ambipolar diffusion may be important. The ambipolar Reynolds number is defined as:

$$Re_\mathrm{A}=\frac{\nu _{\mathrm{ni}}}{\mathrm{\Omega }}$$
(8)
where $`\nu _{\mathrm{ni}}`$ is the frequency of collisions between ions and neutral particles.

We can calculate the both Reynolds numbers locally in the disk. In order to determine the resistivity, $`\eta =(c^2m_\mathrm{e}/4\pi e^2)\times (\nu _{\mathrm{en}}/n_\mathrm{e})`$, we have to estimate the number density of electrons $`n_\mathrm{e}`$ and neutral particles $`n_\mathrm{n}`$. The frequency of collisions between electron and neutral particles is defined as $`\nu _{\mathrm{en}}=8.3\times 10^{10}T^{1/2}n_\mathrm{n}`$ in s-1 (Draine, Roberge & Dalgarno 1983). The ionization fraction $`x_\mathrm{e}=n_\mathrm{e}/n_\mathrm{n}`$ is calculated from the Saha equation for any given value of temperature and density in the disk. In this temperature range the free electrons originate mostly from the ionization of hydrogen, sodium, potassium and calcium.

The magneto-rotational instability acts over the entire disk thickness. Therefore, we need to consider the entire vertical disk structure when we estimate the coupling of the magnetic field to the gas. This is why instead of calculating $`Re_\mathrm{A}`$ and $`Re_\mathrm{M}`$ at the equatorial plane, we introduce the density weighted values of these parameters, integrated over the disk height:

$$\overline{R}e_\mathrm{A}=\frac{_0^HRe_\mathrm{A}\rho (z)𝑑z}{_0^H\rho 𝑑z}$$
(9)
and

$$\overline{R}e_\mathrm{M}=\frac{_0^HRe_\mathrm{M}\rho (z)𝑑z}{_0^H\rho 𝑑z}$$
(10)
where $`\rho (z)`$ is the density profile determined from the vertical structure model and $`Hz_{\mathrm{max}}`$ is the disk thickness.

Numerical simulations (Hawley et al. 1996; Hawley & Stone 1998) show that below the critical Reynolds numbers (defined as in Equations 7 and 8):

$$Re_\mathrm{M}^{\mathrm{crit}}=10^4$$
(11)
and

$$Re_\mathrm{A}^{\mathrm{crit}}=100$$
(12)
the action of MHD turbulence is inefficient and cannot be important in the angular momentum transport in the disk. As a consequence, the viscosity parameter ($`\alpha \alpha _{\mathrm{mag}}`$) is very small.

If the disk Reynolds numbers drop below critical values in the quiescence (e.g. on the lower branch of the S-curve), while rise above the critical values in the outburst (e.g. on the upper branch of the S-curve) then the assumption that $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$ (see Figure 1) is justified. Otherwise, the viscosity parameter should be the same on both cold and hot branches of the S-curve (see Figure 2). Below we calculate the Reynolds numbers in the quiescent disk and test them against the critical values in different types of accreting systems.

### 2.2 Results

#### 2.2.1 Accretion disks in binary systems

We consider an accretion disk in CVs, e.g. around a white dwarf with a mass M=1M and calculate the disk vertical structure at two exemplary radii: 104RSchw and 105RSchw. Figure 4 (a) shows the magnetic and ambipolar Reynolds numbers for the case of small viscosity parameter $`\alpha =0.02`$. It is plotted as a function of the accretion rate corresponding to the lower branch of the S-curve (e.g. in Figure 1 below the point $`A`$). Point $`A`$ indicates the highest accretion rate on the lower, stable branch of the S-curve. Both Reynolds numbers are below the critical values for most part of the lower branch. Magnetic Reynolds number is somewhat above the critical value close to the point $`A`$. The ambipolar Reynolds number is well below the critical value at the smallest radius and approaches it at the point A only for the largest radius. In general, the magnetic number at the lower branch is almost independent from the radius. The ambipolar number is increasing with the radius. However, the value of 105RSchw is rather an upper limit for the disk size in CV systems so any values of the Reynolds number larger than those in Figure 4 are not expected.

In Figure 4 (b) we plot the magnetic and ambipolar Reynolds numbers up to the point $`A`$ for large value of $`\alpha =0.1`$. In this case the Reynolds numbers are mostly small, but as the point $`A`$ reaches this time higher temperatures, also their values can exceed critical ones in the end of the extended stable branch.

Therefore, whenever the viscosity drops at the lower branch it will have a tendency to remain low. This result supports the view that the angular momentum transport due to MHD turbulence cannot be efficient when the disk is in the quiescence and the assumption that $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$ is justified.

We obtain a similar result for black hole X-ray transient systems with MBH $`10M_{}`$ as shown in Figure 5. The calculated magnetic and ambipolar Reynolds numbers fall generally below the critical values for the quiescent disk in these systems, regarding also the fact that this time there is no difference in the location of the point $`A`$ caused by the change in $`\alpha `$.

#### 2.2.2 AGN accretion disks

Now we consider accretion disks in AGN, e.g. systems harboring a supermassive black hole. Note, that for an accretion disk around a supermassive black hole the temperatures are lower than in the stellar mass black hole case. Therefore the ionization instability zone is located much closer to the central black hole and in our calculations we have to consider radii of order of a few hundreds $`R_{\mathrm{Schw}}`$ instead of hundreds of thousands.

We calculate the Reynolds numbers for a range of masses between $`10^6M_{}`$ and $`3\times 10^9M_{}`$. We found that the Reynolds numbers are orders of magnitude above the critical values, for the most part of the lower branch, and the greater the mass of a black hole, the larger $`Re`$. In Figure 6 we plot magnetic and ambipolar Reynolds numbers as a function of the accretion rate corresponding to the lower, cold branch of the S-curve for $`3\times 10^9M_{}`$.

Both Reynolds numbers clearly exceed the critical values: 102 and 104 for $`Re_\mathrm{A}`$ and $`Re_\mathrm{M}`$ respectively. Therefore only for very low accretion rates, $`\dot{m}<10^4`$, MHD turbulence in the disk would not develop and practically the entire cold branch of the S-curve is MHD turbulent. This result only weakly depends on the adopted viscosity parameter value. For $`\alpha =0.02`$ both Reynolds numbers had similar values to these presented in Figure 6 and were always greater than critical ($`\overline{R}e_\mathrm{M}>10^6`$ and $`\overline{R}e_\mathrm{A}>10^5`$ at the disk radius of $`\mathrm{log}R=16.5`$ \[cm\]) .

In order to understand whether the above result is valid for the entire disk we calculate the Reynolds numbers at different locations in the following way. At each radius we compute the entire S-curve, such as that in Figure 2, determine the position of the turning point $`A`$ and calculate the Reynolds numbers $`\overline{R}e_\mathrm{A}(A)`$ and $`\overline{R}e_\mathrm{M}(A)`$ at this point. Their values measured along the lower branch of the S-curve are lower (equal) to $`\overline{R}e_\mathrm{A}(A)`$ and $`\overline{R}e_\mathrm{M}(A)`$. We summarize this result in Figure 7, plotting the Reynolds numbers against the radius. Both $`\overline{R}e_\mathrm{A}`$ and $`\overline{R}e_\mathrm{M}`$ increase with radius implying that MRI is stronger further out in the disk. They decrease towards smaller radii but do not drop below the critical values even at 10 $`R_{\mathrm{Schw}}`$.

Our results support the view presented in Menou & Quataert (2001) and imply that in the thermally unstable AGN disks the viscosity should not change between the cold and hot states of the disk (i.e. $`\alpha _{\mathrm{cold}}=\alpha _{\mathrm{hot}}`$).

## 3 Evaporation of the disk

We concluded above that the viscosity parameter in AGN accretion disks does not vary between the cold and hot states of the disk. As a result the ionization instability will cause only small amplitude luminosity fluctuations not really detectable in the observations (Siemiginowska, Czerny & Kostyunin 1996; Menou & Quataert 2001). However, variations in the local accretion rate induced by the instability may result in the disk evaporation and affect the disk luminosity (Siemiginowska 1998). Here, we investigate the effects of the disk evaporation on the overall evolution of thermally unstable disk. We assume that $`\alpha _{\mathrm{cold}}=\alpha _{\mathrm{hot}}=0.1`$.

### 3.1 Evaporation radius

The accretion flow in low luminosity systems can occur either via a standard cold disk or via an advection dominated flow (ADAF). In the latter the bulk of gravitational energy is advected into the central object instead of being radiated (Ichimaru 1977; Narayan & Yi 1994). The ADAF solution forms below a certain transition radius that depends on the accretion rate. In classical ADAF model (Abramowicz et al. 1995; Honma 1996; Kato & Nakamura 1998) this relation is:

$$R_{\mathrm{evap}}=1.9\dot{m}^2\alpha _{0.1}^4R_{\mathrm{Schw}}.$$
(13)
The transition radius between the outer disk and inner ADAF was discussed in several papers (e.g. Liu et al. 1999; Różańska & Czerny 2000; Menou et al. 2000; Liu et al. 2002; Meyer & Meyer-Hofmeister 2002), under various assumptions. Here we adopt the canonical approach.

The ionization instability can develop only when the standard thin disk reaches the temperatures sufficient for the hydrogen ionization. If accretion rates are low the disk can evaporate without being ionized (e.g. without reaching point $`A`$ on the stability S-curve) and the instability does not develop. Therefore the evaporation radius cannot be very large for the ionization to occur:

$$R_{\mathrm{evap}}<R_{\mathrm{instab}}.$$
(14)
In Figure 3 we showed the radial extension of the unstable ionization zone as a function of the accretion rate. The instability cycle operates between the points $`A`$ and $`B`$ of the S-curve (cf. Figure 2), so we consider this entire range of accretion rates to be unstable. In Figure 3 we also indicate a location of $`R_{\mathrm{evap}}`$ for the ADAF prescription (note, that the evaporation radius does not depend on mass of a black hole). If the disk evaporates, only the unstable region above the solid line may exist and contribute to the total disk luminosity. The extension of the unstable region does not depend on $`\alpha `$ very much, while the formula defining the evaporation radius strongly depends on it. Therefore, we should emphasize that our result may be very sensitive to the adopted value of the viscosity parameter. Possibly, depending on $`\alpha `$, the entire unstable region of the disk evaporates, or, on the other hand, the evaporation radius is very small and does not affect the instability.

Because the location of the evaporation radius depends on viscosity, if it can be determined observationally it may provide a new method of estimating the value of the $`\alpha `$ parameter in accreting objects. For example, recently Nayakshin & Sunyaev (2003) proposed a new model explaining the X-ray flares in the Galactic Center by the star-disk interactions. This model would favor the case of the disk extending quite close to the central black hole. However, the presence of the disk in the Galactic Center is still unclear (see Quataert 2003 for a review).

### 3.2 Luminosity of the disk with inner ADAF

Below we study the evolution of the disk affected by the ionization instability assuming that the viscosity parameter is constant, as favored for accreting supermassive black holes. We used the method described in Siemiginowska, Czerny & Kostynuin (1996) and consider the two cases: (i) the accretion disk extends down to the marginally stable orbit and (ii) the central part of the disk is evaporated and forms an ADAF, while the standard thick disk remains at the outer radii. Since the accretion efficiency of an ADAF is low, we assume that this part of the flow does not contribute to the luminosity.

In order to determine the disk evaporation radius during the evolution we compare the local accretion rate with the critical (evaporation) accretion rate. We use the classical ADAF prescription, expressed by Equation 13. In Figure 8 we show radial profiles of the local accretion rate in several snapshots during one cycle of the disk instability. The straight solid line marks the critical accretion rate. Whenever the local accretion rate in the disk drops below this critical value, the disk evaporates and ADAF is formed. Hereafter, in our time dependent calculations we assume that the ADAF transition radius is the outermost radius at which the disk evaporates. Therefore even though the local accretion rate may be equal to the critical value at several locations in the disk, only the outermost one determines the transition radius and the standard disk cannot rebuilt itself below this radius.

In Figure 9 we present lightcurves for the two models: (i) the standard disk and (ii) the disk with inner ADAF. When the disk extends all the way to the marginally stable orbit, as considered in the previous time-dependent calculations, the luminosity variations are small and only a small flickering ($`\mathrm{\Delta }\mathrm{log}L_{\mathrm{disk}}0.2`$) can be seen in the resulting lightcurve. This is indicated by the upper solid line in Figure 9.

However, the situation changes dramatically, when the evaporation of the inner disk is taken into account. In this case the ionization instability leads to strong, large-amplitude outbursts. These outbursts are due to the large variations of the inner radius of the accretion disk (e.g. the size of the ADAF), as determined by the critical value of the local accretion rate. The outburst peak corresponds to the moment, when the inner disk radius is small and the entire accretion disk contributes to the total luminosity. The quiescent low luminosity state, on the other hand, is achieved when the local accretion rate drops below the critical value at the large distance from the center and most of the inner disk is evaporated.

## 4 Discussion and conclusions

We studied the effects of MHD turbulence and the efficiency of angular momentum transport in the quiescent accretion disks. We analyzed the full disk vertical structure, which enabled us to determine accurately the coupling of the accreting gas to the magnetic field over the entire disk height.

Confirming the previous results of Menou & Quataert (2001), we find that only in the case of accretion disks in binary systems, e.g. around white dwarfs (CVs) or stellar mass black holes (SXTs), the timescale of the magnetic field decay is too short for the field to be frozen into the gas. Therefore, MHD turbulence does not operate efficiently and the large amplitude outbursts observed in these sources can be accounted for by the thermal ionization instability cycle in which $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$. On the other hand, in the case of accreting supermassive black holes in AGN the gas in the disk is well coupled to the magnetic field.

This effect ($`Re_\mathrm{M}1`$ and $`Re_\mathrm{A}1`$) can be only slightly reduced by encounters of free electrons with charged dust particles. In principle, the presence of dust in the disk can affect the number density of free electrons. However, the grains evaporate if the temperature exceeds $`18002000`$ K, which is the case near the disk equatorial plane. As a result, the role of the dust in suppressing MHD turbulence is limited only to the uppermost disk layers and cannot have a global effect. We conclude, that the AGN accretion disks in quiescence have comparable angular momentum transport efficiency to the disks in outburst.

Therefore, we studied the ionization instability assuming a constant value of the viscosity parameter $`\alpha _{\mathrm{cold}}=\alpha _{\mathrm{hot}}`$. We took into account the evaporation of the innermost parts of the disk during the quiescence in AGN accretion disks. The transition to an ADAF below the evaporation radius during the evolution of the instability results in large luminosity variations ($`\mathrm{\Delta }\mathrm{log}L_{\mathrm{disk}}2`$). Because only the outer parts of the disk ($`R>R_{\mathrm{evap}}10^3R_{\mathrm{Schw}}`$) contribute to the total disk emission the minimum luminosity obtained within the instability cycle is lower than in the standard disk with no evaporation. We note, that also during the outburst the evaporation is important and ADAF remains in the innermost regions of the disk ($`R_{\mathrm{evap}}25R_{\mathrm{Schw}}`$), even at quite a high global accretion rate, $`\dot{M}=0.1\dot{M}_{\mathrm{Edd}}`$. The evaporation radius, however, strongly depends on the adopted value of the viscosity parameter $`\alpha `$ (cf. Equation 13).

The $`\alpha `$-disks in AGN are self-gravitating beyond a certain radius (Paczyński 1978; Hure 1998). A typical distance of the self-gravitating zone is about a thousand Schwarzschild radii (Hure 2000). However, the greater the black hole mass and lower the parameter $`\alpha `$, the closer to the center this zone is located. We can include a self-gravitation term equal to $`4\pi G\mathrm{\Sigma }`$ in Equation 3 of the disk vertical structure. This additional term modifies the S-curve, so there is a single ’wiggle’ on the middle branch of the curve, instead of two. The middle branch is also broader and reaches lower surface densities. The instability happens for a narrower range of effective temperatures, while the more strongly curved unstable middle branch covers broader range of surface density values. Therefore, the amplitude of the outburst is not expected to be significantly modified in our model by the presence of self-gravity. However, whether the model actually can be extended into the self-gravity region, remains an open question. In particular, our description of convection, which influences the shape of the S-curve, may not be applicable in case of a self-gravitating disk. Also, MRI was not extensively studied in the self gravity region (but see Fromang et al. 2003) and the use of viscosity parameter $`\alpha `$ may not be appropriate here. Instead, self-gravity may either provide a mechanism of angular momentum transfer (e.g. Paczyński 1978), or lead to the disk fragmentation and star formation (Collin & Zahn 1999). The physics of those processes is poorly known. Therefore the use of simple models like $`\alpha `$-disks and comparison of their predictions with observed behavior of AGN seems to be justified.

The variability amplitudes predicted by our model are about two orders of magnitude and support the observational evidence of intermittent activity of AGN. This kind of activity has been proposed to explain the formation of giant radio galaxies (Subrahmanyan et al. 1996) and the extended radio structures in the Giga-Hertz Peaked Spectrum (GPS) radio sources (Baum et al. 1990; Schoenmakers et al. 1998). The timescales of episodic activity in quasars were recently estimated by Martini & Schneider (2003). Recent observations of Compact Steep Spectrum (CSS) sources indicate that the evolutionary path of radio loud AGN is consistent with an on-off activity with duration governed by the mass of the central black hole (Marecki et al. 2003). Also the black hole masses implied by the HST observations of normal galaxies and AGN suggest that every galaxy must undergo an active phase in its life (Ferrarese et al. 2001).

The estimated timescales of the activity cycle are within $`10^5`$ to 107 years (e.g. Reynolds & Begelman 1997 for GPS sources). Di Matteo et al. (2003) proposed a feedback mechanism to explain the cyclic activity in M87, which occurs on a timescale of 107 years. In our time-dependent calculations, we adopted a simplified approach to the disk evaporation process, and therefore underestimated the repetition timescale. This is because we do not consider the time of the disk reconstruction after it has been evaporated, which increases the viscous timescale by a factor of $`2`$. Assuming that the outburst starts at the outer edge of the instability zone we obtain the viscous timescale of $`t_{\mathrm{visc}}=1/(\alpha \mathrm{\Omega })\times (r/H)^20.9`$ Myrs for $`M_{\mathrm{BH}}=10^8M_{}`$. Thus the outbursts will be less frequent than these presented in Figure 9. However, because during the outburst there is no ADAF and the entire disk is present, the decay timescale is calculated properly and the duration of an outburst is correct. On the other hand, the outburst profiles will be less steep than the ones calculated here, since the disk density rises gradually at the beginning of the instability cycle. Although these two effects are important, the extension of the unstable region is not very large in comparison with the disk size and the uncertainty involved in our calculations is not crucial. The detailed model of the time evolution of the accretion disk with transition to ADAF is clearly worth investigating and we address it to the future work.

Acknowledgments AJ is grateful to F. Meyer and E. Meyer-Hofmeister for interesting discussion. This work was supported in part by grants 2P03D00322 and 2P03D00124 of the Polish State Committee for Scientific Research. Partial support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number GO-01164X issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073.

## References

Abramowicz M.A., Chen X., Kato S., Lasota J.-P., Regev O., 1995, ApJ, 438, L37

Alexander D.R., Johnson H.R., Rypma R.L., 1983, ApJ, 272, 773

Balbus S.A., Hawley J.F., 1991, ApJ, 376, 214

Balbus S.A., Hawley J.F., 1998, Accretion Processes in Astrophysical Systems: Some Like it Hot!, eds. Holt S.S. & Kallman T.R., AIP Conf. Proceedings, 431, p. 79

Baum S.A., O’Dea C.P., de Bruyn A.G., Murphy D.W., 1990, A&A, 232, 19

Begelman M.C., 1985, Astrophysics of Active Galaxies and Quasi-Stellar Objects, p. 411 (Mill Valley: University Science Books)

Cannizzo J.K., Ghosh P., Wheeler J.C., 1982, ApJ, 260, L83

Cannizzo J.K., 1993, ApJ, 419, 318

Clarke C.J., 1988, MNRAS, 235, 881

Collin S., Zahn J.-P., 1999, A&A, 344, 433

Di Matteo T., Allen S.W., Fabian A.C., Wilson A.S., Young A.J., 2003, ApJ, 582, 133

Dubus G., Hameury J.-M., Lasota J.-P., 2001, A&A, 373, 251

Ferrarese L., Pogge R.W., Peterson B.M., Merritt D., Wandel A., Joseph C.L., 2001, ApJ, 555, L79

Fromang S., de Villiers J.P., Balbus S.A., 2003 (astro-ph/0309816)

Gammie C.F., Menou K., 1998, ApJ, 492, L75

Hatziminaoglou E., Siemiginowska A., Elvis M., 2001, ApJ, 547, 90

Hawley J.F., Gammie C.F., Balbus S.A., 1996, ApJ, 464, 690

Hawley J.F., Stone J.M., 1998, ApJ, 501, 758

Honma F., 1996, PASJ, 48, 77

Hure J.-M., 1998, A&A, 337, 625

Hure J.-M., 2000, A&A, 358, 378

Ichimaru S., 1977, ApJ, 214, 840

Janiuk, A., Czerny, B., & Siemiginowska, A., 2002, ApJ, 576, 908

Kato S., Nakamura K.E., 1998, PASJ, 50, 559

Lasota J.-P., 2001, NewAR, 45, L449

Lin D.N.C., Shields G.A., 1986, ApJ, 305, 28

Liu B.F., Yuan W., Meyer F., Meyer-Hofmeister E., Xie G.Z., 1999, ApJ, 527, L17

Liu B.F., Mineshige S., Meyer F., Meyer-Hofmeister E., Kawaguchi T., 2002, ApJ, 575, 117

Marecki A., Spencer R.E., Kunert M., 2003, PASA, 20, 46

Martini P., Schneider P., 2003 (astro-ph/0309650)

Menou K., Hameury J.-M., Lasota J.-P., Narayan R., 2000, MNRAS, 314, 498

Menou K., Quataert E., 2001, ApJ, 552, 204

Meyer F., Meyer-Hofmeister E., 1981, A&A, 104, L10

Meyer F., Meyer-Hofmeister E., 1982, A&A, 106, 34

Meyer F., Meyer-Hofmeister E., 1994, A&A, 288, 175

Meyer F., Meyer-Hofmeister E., 2002, A&A, 392, L5

Mineshige S., Shields G.A., 1990, ApJ, 351, 47

Narayan R., Yi I., 1994, ApJ, 428, L13

Nayakshin S., Sunyaev R., 2003, MNRAS, 343, L15

Parker E.N., 1979, Cosmical Magnetic Fields, (Oxford:Clarendon Press)

Paczyński B., 1978, Acta Astron., 28, 91

Pojmański G., 1986, Acta Astron. 36, 69

Quataert, 2003 in Astron. Nachr., Vol. 324, No. S1 (astro-ph/0304099)

Rees M.J., Begelman M.C., Blandford R.D., Phinney E.S., 1982, Nature, 295, 17

Reynolds C.S., Begelman M.C., 1997, ApJ, 487, L135

Różańska A., Czerny B., Życki P.T., Pojmański G., 1999, MNRAS, 305, 481

Różańska A., Czerny B., 2000, A&A, 360, 1170

Schoenmakers A.P., de Bruyn A.G., Rottgering H.J.A., van der Laan H., 1999, A&A, 341, 44

Seaton M.J., Yan Y., Mihalas D., Pradhan A.K., 1994, MNRAS, 266, 805

Shakura N.I., Sunyaev R.A., 1973, A.&A., 24, 337

Siemiginowska A., Czerny B., Kostyunin V., 1996, ApJ, 458, 491

Siemiginowska, A. 1998, Accretion Processes in Astrophysical Systems: Some Like it Hot!, eds. Holt S.S. & Kallman T.R., AIP Conference Proceedings 431, p. 211

Smak J. 1982, Acta Astron. 32, 199

Smak J. 1984, Acta Astron. 34, 161

Subrahmanyan R., Saripalli L., Hunstead R.W., 1996, MNRAS, 279, 257


text:
# The internal magnetic field in superconducting ferromagnets

\[

## Abstract

We have measured the nonlinear response to the ac magnetic field in the superconducting weak ferromagnet Ru-1222, at different regimes of sample cooling which provides unambiguous evidence of the interplay of the domain structure and the vorticity in the superconducting state. This is direct proof of coexistence of ferromagnetic and superconductive order parameters in high-$`T_c`$ ruthenocuprates.

\]

The problem of coexistence of superconductivity (SC) and ferromagnetism (FM) has been studied for almost 50 years starting from the theoretical work by Ginzburg  (see also ). Coexistence of weak-ferromagnetism (W-FM) and SC was discovered some time ago in RuSr2R2-xCexCu2O10 (R=Eu and Gd, Ru-1222) layered cuprate systems, and more recently  in RuSr2GdCu2O8 (Ru-1212). The SC charge carriers originate from the CuO2 planes and the W-FM is related to the Ru layers. In both systems, the magnetic order does not vanish when SC sets in at $`T_c`$, and remains unchanged and coexists with the SC state. The Ru-1222 materials (for R=Eu and Gd) display a magnetic transition at $`T_N=125180`$ K and bulk SC below $`T_c`$ = 25-50 K ($`T_N>T_c`$) depending on the oxygen concentration and sample preparation. This discovery has launched a new wave of investigations in this field . The problem is of general interest for condensed matter physics and is relevant for many materials, in particular unconventional superconductors (including some heavy fermions) with triplet pairing .

Despite a lot of work done in the past and recently, debates concerning whether such coexistence is genuine are still continuing. Evidence in favor of coexistence is mostly indirect and refers to some peculiarities of the magnetization curve. One of the most pronounced manifestations of SC-FM coexistence is the spontaneous vortex phase (superconducting vortices induced by the internal magnetic field from the FM magnetization). It explains well the magnetization curve of these materials (see  and references therein). However, this phase has not yet been observed experimentally (visualized as the more common mixed state of type II superconductors).

In the past the evidence of the FM-SC coexistence referred mostly to the magnetic properties of the materials affected by the presence of superconductivity. In this letter we present the first experimental evidence of the effect of the ferromagnetic order parameter, on the superconducting order parameter. The ferromagnetic order parameter, namely the spontaneous magnetization, is a source of an internal magnetic field inside a sample even without an external magnetic field $`H`$. On the other hand, the superconducting properties of type-II superconductors depend strongly on whether the sample was cooled to the SC state in zero magnetic field (ZFC) or in a finite magnetic field (FC). Here a “field” is supposed to be an external magnetic field. We show here that these properties depend also on the internal magnetic field during the cooling process. We exploited the procedure, which we shall call the internal-field cooling (IFC): The sample was cooled down to $`T_{IFC}`$ under an external magnetic field $`H_{IFC}`$, ($`T_{IFC}<T_N`$). At $`T_{IFC}`$, $`H_{IFC}`$ was turned off and further cool-down to $`T=5`$ K was done at $`H=0`$. It appears that, by using the IFC procedure, the properties of the SC state were different from those measured after the regular ZFC process from temperatures above $`T_N`$. Thus, in the SC state, the sample senses the internal magnetic field evolved from the remanent magnetization, which was formed in the normal ferromagnetic phase and then frozen at further cooling.

We measured the nonlinear response to the ac magnetic field, which is a sensitive probe of superconducting vorticity, as demonstrated by numerous investigations in the past . Ceramic sample of Gd1.5Ce0.5Ru2Sr2Cu2O10 (Ru-1222) with dimensions $`8\times 2\times 2`$ mm3 was prepared by a solid-state reaction as described in Ref. 2. In a nonlinear medium, magnetization oscillations, induced by an ac magnetic field $`h(t)=h_0\mathrm{sin}\omega t`$, may be expanded in a Fourier series:

$$M(t)=h_0\underset{n>0}{}\chi _n^{^{}}\mathrm{sin}(n\omega t)\chi _n^{^{\prime \prime }}\mathrm{cos}(n\omega t)$$
(1)
where $`\chi _n^{^{}}`$ and $`\chi _n^{^{\prime \prime }}`$ ($`n=1,2,3\mathrm{}`$) are the in-phase and out-of-phase components of the harmonic susceptibility. In all experiments described here we measured the voltage drop induced in a pickup coil, which is proportional to the time derivative of $`M(t)`$. Our home made experimental setup was adapted to a commercial MPMS SQUID magnetometer. An ac field $`h(t)`$ at a frequency of $`\omega /2\pi =1.5`$ kHz and an amplitude up to the $`h_0=3`$ Oe was generated by a copper solenoid existing inside the SQUID magnetometer. The temperature, dc magnetic field, and amplitude dependencies of the fundamental and third harmonic signals presented here have been measured by the two coils method  . In the present letter the results for the first and third harmonics will be discussed.

Figure 1 shows the temperature dependencies of the in-phase susceptibility $`\chi _1^{^{}}`$ and of the amplitude of the third harmonic $`A_{3\omega }h_0|\chi _3^{^{}}i\chi _3^{^{\prime \prime }}|`$, measured after the ZFC process at $`H=0`$. The temperature dependence of $`\chi _1^{^{}}`$ is typical for superconducting ferromagnets . This plot reveals three transitions: (i) the paramagnetic-antiferromagnetic transition at $`T_N125`$ K, (ii) the most pronounced transition, which corresponds to the peak at $`T_m78`$ K, and (iii) the transition into the SC state at $`T_c28`$ K. The nature of the second transition, which is evident both in the linear and the nonlinear response, is not yet completely clear and is discussed elsewhere . Ambiguity is connected with the magnetic phase between $`T_m`$ and $`T_N`$, which is characterized by low coercivity. On the other hand, the $`T_c<T<T_m`$ temperature region definitely corresponds to the weak ferromagnetic phase .

The third harmonic behavior is different for $`T>T_c`$ and $`T<T_c`$. For $`T>T_c`$ the behavior is typical for ferromagnetic materials and was known already from Rayleigh’s investigation on iron . The third harmonic response demonstrates a quadratic dependence on $`h_0`$ (inset in Fig. 2), which directly derived from the oscillatory motion of the domain walls . This signal should decrease at low temperatures and it becomes unobservable under our experimental conditions. For $`T<T_c`$ the third harmonic grows very fast with temperature decreasing (Fig. 1), and its dependence on the ac field amplitude (Fig. 2) is different from that at $`T>T_c`$, as evident from the saturation for the nonlinear response at high amplitude of excitation, instead of a quadratic growth. The growth of the nonlinear response in the superconducting materials was revealed in numerous previous experimental investigations . Various mechanisms were suggested for this nonlinear response based on the critical state model  and the presence of weak links . In particular, the response shown in Fig. 2 is well described by the Josephson-media model. We do not have to discuss these models, since all of them relate the response to the penetration of the magnetic flux (vortices) into the sample, and only this fact is essential for the present investigation. Thus it seems reasonable that the $`A_{3\omega }`$ at $`T<T_c`$ is an effective probe of the superconducting vorticity.

Figure 3 demonstrates ZFC dependence of $`A_{3\omega }`$ on the external magnetic field. One can see that $`A_{3\omega }`$ decreases with the magnetic field. At high magnetic fields $`A_{3\omega }`$ is a power function of the $`H`$: $`A_{3\omega }H^q`$, with $`q0.8`$. Suppression of the $`A_{3\omega }`$ by the magnetic field applied after ZFC was observed in the previous nonlinear studies and agrees with all suggested models of the nonlinear response. The nonlinearity under discussion is connected with a nonhomogeneous distribution of the magnetic flux, which penetrated into the sample, and the magnetic flux distribution becomes more and more uniform, when the vortex density increases. On the other hand, in the Meissner state the nonlinear response must be quite weak, and the magnetic field dependence of $`A_{3\omega }`$ should have a peak at some $`H`$, as was observed in some materials . But in ceramics with numerous weak links, such as our material, this field can be extremely small, and the peak is not observable. Moreover, we deal with the superconducting ferromagnets, where the spontaneous vortex phase can replace the Meissner state at $`H=0`$. Altogether this explains why we observe the maximum value of $`A_{3\omega }`$ at $`H=0`$.

Now let us consider the experimental results in the IFC process. After turning off the magnetic field $`H_{IFC}`$ at temperature $`T_{IFC}`$, the sample was cooled in $`H=0`$ down to $`T=5`$ K and the signal of the third harmonic at $`T=5`$ K was measured. Figure 4 shows $`A_{3\omega }(H_{IFC})`$ dependence for $`T_{IFC}=40`$ K and $`70`$ K. It is evident that the field $`H_{IFC}`$ suppresses the $`A_{3\omega }`$ signal similarly to the external field after ZFC Fig. 3 even though $`H_{IFC}`$ was turned off before the onset of superconductivity. Turning off $`H_{IFC}`$ at $`T=40`$ K affects $`A_{3\omega }`$ more strongly than for $`T=70`$ K due to larger remanent magnetization at $`T=40`$ K. This behavior is typical for the FM materials .

Figure 5 presents the signal of the third harmonic $`A_{3\omega }(T=5K)`$ as a function of $`T_{IFC}`$ after cooling in $`H_{IFC}=30`$ Oe. The signal of the $`A_{3\omega }(T=5K)`$ decreases for $`T_{IFC}<T_m`$. This demonstrates that the suppression of the third harmonic response by the internal magnetic field takes place only if the field cooling continues down to the weakly ferromagnetic phase with essential coercivity. It is known  that in idealized single-domain superconducting ferromagnets the internal magnetic field from the spontaneous magnetization $`4\pi \stackrel{}{M}`$ has the same effect on the phase diagram, i.e., on the magnetic flux penetrating into the sample, as the external field. This can be generalized in the more realistic case of a multi-domain sample with nonzero average internal field $`4\pi \stackrel{}{M}`$. On the basis of this argument we can use plot of $`A_{3\omega }(H)`$ (Fig. 3) as a calibration curve to estimate the magnitude of the frozen internal magnetic field ($`H_I`$). Namely, we take the value of $`A_{3\omega }`$ from the plot in Fig. 4, find the value of $`H`$, which corresponds to this value of $`A_{3\omega }`$ in Fig 3, and assume that this value of $`H`$ gives a reasonable estimation of $`H_I`$. Figure 6 presents the dependence of $`H_I`$ on $`H_{IFC}`$.

The internal magnetic field arises from the frozen remanent magnetization $`4\pi \stackrel{}{M}`$ after field cooling down to $`T_{IFC}`$. We have compared obtained in Fig. 6 with direct dc remanent magnetization measured in our previous studies . It appears that there is a reasonable agreement (with an accuracy of $`\pm 20\%`$) between the two methods, and confirms our scenario.

The phenomenon revealed in our experiment is possible if the domain structure formed in the ferromagnetic phase can be frozen down to the superconducting state. On the other hand, as was noted in the pioneering paper by Ginzburg  and confirmed by the detailed analysis in Ref. , superconductivity should strongly affect the equilibrium domain structure: Its period should grow, and in equilibrium any sample in the Meissner state is a single domain. But in our case we deal with a non-equilibrium domain structure, which is a metastable state due to coercivity.

The presence of the frozen internal field in the superconducting phase clearly demonstrates that the sample is in the mixed state with many vortices. One cannot call this state the spontaneous vortex phase because the latter refers to the equilibrium state, but we deal with a metastable state. We have analyzed here the nonlinear response, which is sensitive to the average internal field $`4\pi \stackrel{}{M}`$. The absolute value of the average magnetization $`\stackrel{}{M}`$ is less than the saturation magnetization $`M`$, which can determine the vortex density in a single-domain sample . However, the saturation magnetization may create vortices inside domains. Since $`\stackrel{}{M}`$ changes its direction from domain to domain, we obtain the vortex tangle, which does not contribute to the average internal field $`4\pi \stackrel{}{M}`$, studied here. This vortex tangle is expected to exist even after the ZFC process and contributes to the initial value of the third harmonic, which was detected without external or internal magnetic field. These arguments illustrate that the vorticity (magnetic flux) distribution in a real (especially ceramic) superconducting ferromagnet can be very complicated. Genuinely zero field cooling is practically impossible: if one cools a sample in zero external field, one cannot avoid internal magnetic fields from the spontaneous magnetization, even if these fields vanish on average but still remain inside the domains. A more detailed analysis of the magnetic-flux distribution would become possible if further investigations provided more information on the structure of the material: sizes of grains and domains, data on crystal anisotropy etc.

In summary, our measurements of the nonlinear response unambiguously demonstrate the coexistence of the superconducting and ferromagnetic order parameter in Ru-1222 samples below the superconducting critical temperature. Coexistence is manifested by the clear effect on the domain structure frozen from normal FM phase on superconducting properties. We tend to believe that the effect revealed in Ru-1222 is a general and can be observed in other materials with FM-SC coexistence.

The work was supported by the Klatchky Foundation and by a grant from the Israel Academy of Sciences and Humanities.