text:
# Colloidal interactions in two dimensional nematics

## I Introduction

In the last few years, there has been considerable theoretical and experimental interest in colloidal systems with anisotropic host fluids .

Since the temperature changes the anisotropic properties of the host fluid, the dispersed particles exhibit a rich variety of collective behaviors, leading to unexpected phase diagrams . Depending on the details of the system, chain-like , crystal  or cellular  structures have been observed.

In the case of nematic hosts in three dimensions (3D), the colloidal interactions result from the competition between the anchoring properties of the spherical colloidal surfaces favoring e.g., radial (homeotropic) nematic orientation (anchoring) and the bulk elasticity that favors a uniform nematic. For sufficiently strong homeotropic anchoring $`W`$, the spherical colloids may induce topological defects or singularities in the nematic director field. If the size of the colloidal particle $`a`$ is large when compared to the nematic correlation length $`\zeta `$, a hyperbolic point defect appears at a certain distance from the particle (satellite dipolar configuration ). For smaller particles, however, a ring disclination around the equatorial plane of the sphere may have a lower energy (saturn-ring quadrupolar configuration ). Finally, if $`Wa`$ is smaller than the typical nematic elastic constant $`K`$, the director field varies smoothly and there are no defects (surface-ring quadrupolar configuration ). Depending on the symmetry of the distortion , the dominant long range interactions are dipolar  or quadrupolar . The short range interactions are more difficult to predict, as nonlinear terms in the elastic free energy come into play. Experimental evidence from a system where silicon oil droplets were dispersed in a nematic host  reveals that when the nematic configuration around each droplet has dipolar symmetry, the attractive long range interaction between aligned droplets becomes strongly repulsive at short distances. A nematic topological defect appears between the particles and prevents their collapse. By contrast, when the long range distortion exhibits quadrupolar symmetry no such defects appear and the oil droplets coalesce.

2D colloidal systems were also investigated recently and similar qualitative features were observed. In these systems the interactions between the dispersed particles are also mediated by the elastic distorsions of the anisotropic host. However, for 2D nematic hosts and homeotropic anchoring conditions, only quadrupolar configurations are stable . Dipolar configurations may be found in 2D smectic C films , a system that is similar to a 2D nematic, but where the director does not exhibit mirror symmetry, thus excluding configurations with half integer topological charges. In this case, Pettey et al.  obtained analytic expressions for the free energy of a system of multiple disk-defect pairs, valid in the linear regime where the particles are sufficiently far apart. We investigated the short range interactions in this system by using a numerical method capable of describing the nematic orientation profiles induced by a small number of colloids with arbitrary shape in 2D . Numerical investigations of this type are faced with a problem due to the very different length scales characterizing the colloids and the topological defects. We have solved it , by using finite elements with adaptive meshing to minimize the Frank free energy and obtained very accurate interaction energies between disks at arbitrary separations. We have found a strong repulsion at short range followed by the expected long range dipolar interaction with a pronounced minimum at an intermediate distance. The equilibrium director field configuration exhibits a topological defect at the mid-point between the disks.

In this article, we use similar numerical methods  to minimize the Landau-de Gennes free energy of two circular disks in a 2D nematic. The great advantage of this model is that no special treatment of the topological defects is required, by contrast with the Frank elastic free energy description, where the defects appear as singularities of the director field $`𝐧`$. In the next section, we describe the nematic free energy based on the tensor order parameter formalism. Then, we generalize the analytical results of Pettey et al.  to nematic configurations with half integer topological defects and obtain analytic expressions for the free energy of the quadrupolar nematic distortion for a single disk and the long range quadrupolar interaction between disks. Finally, we present the results of extensive numerical calculations for the interaction between disks at arbitrary separations.

## II The Landau-de Gennes free energy

In the tensor order parameter formalism the 2D nematic order is represented by the (tilt) angle of the director field $`𝐧=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$, that describes the mean orientation of the molecules, and the orientational order parameter $`Q`$, that measures the degree of alignment of the molecules with respect to the mean orientation. In the isotropic phase, where the molecules point in every direction, $`Q=0`$. By contrast, in the nematic phase the molecules are preferentially aligned and the orientational order parameter takes a well defined equilibrium value $`Q=Q_{eq}`$ that depends on the temperature.

We define a 2D symmetric traceless tensor field, from the tilt angle $`\theta (𝐫)`$ and the orientational order parameter $`Q(𝐫)`$,

$$Q_{ij}(𝐫)=Q(𝐫)\left(n_i(𝐫)n_j(𝐫)\delta _{ij}/2\right)$$
(1)
that is invariant under the nematic symmetry operations $`𝐧𝐧`$ and $`𝐫𝐫`$.

The Landau-de Gennes free energy may be understood as a Taylor expansion in terms of the tensor order parameter $`Q_{ij}`$ (the bulk free energy) and its derivatives $`_kQ_{ij}`$ (the elastic free energy) :

$$F=_\mathrm{\Omega }\left(\frac{A}{2}Q_{ij}^2+\frac{C}{4}(Q_{ij}^2)^2+\frac{L}{2}(_kQ_{ij})^2\right)d^2𝐫$$
(2)
where the summation over repeated indices is implied. The integration is over the 2D space occupied by the nematic $`\mathrm{\Omega }`$. Only terms with the nematic symmetry are allowed in the free energy expansion. In addition, stability requires that the total free energy is bounded from below. In the nematic $`A`$ and $`C`$ are positive constants, defining the equilibrium value of the orientational order parameter $`Q_{eq}=\sqrt{2A/C}`$. For simplicity, we adopted the one-constant approximation for the elastic free energy. In the general 2D case two elastic terms – corresponding to different invariant derivatives of the tensor order parameter – have to be included in the expansion. The elastic constant $`L`$ is related to the Frank elastic constant through $`K=4LA/C`$. The total free energy can be written as a functional of the tilt angle $`\theta `$ and the orientational order parameter $`Q`$:

$$F=_\mathrm{\Omega }\left(\frac{A}{4}Q^2+\frac{C}{16}Q^4+\frac{L}{4}(Q)^2+LQ^2(\theta )^2\right)d^2𝐫.$$
(3)
Let $`a`$ be a typical length of the system, e.g. the radius of the colloidal particles. If $`L/a^2A`$, the energy associated with the variation of the tilt angle $`\theta `$ is much smaller than the energy associated with the orientational order parameter. In this case we may set $`Q=Q_{eq}`$ everywhere in the nematic and the free energy simplifies to

$$F_{el}=\frac{K}{2}_\mathrm{\Omega }(\theta )^2d^2𝐫$$
(4)
which is the one-constant elastic Frank free energy.

For strong distortions, however, topological defects may appear. In this case the director field becomes singular and within a small region around the defect the orientational order parameter vanishes and the contribution of the first three terms in Eq.(3) has to be taken into account. This energy – the core energy of the defect – arises from a small region of size $`\xi `$, and may be estimated by a boundary layer analysis .

Let us consider a topological defect of charge $`q`$ placed at the center of the 2D space. In polar coordinates $`𝐫=(r,\varphi )`$, the defect may be described by $`\theta (\varphi )=q\varphi `$, asymptotically valid in the core region $`r<\xi `$, far from any other boundaries. The orientational order parameter, constant outside the core region, decreases rapidly to $`Q(r=0)=0`$ inside the core in order to regularize the nematic free energy. A simple description of this behavior of $`Q`$ is given by the ansatz

$$Q(r)=Q_{eq}(1e^{\frac{r}{\xi }}).$$
(5)
By substituting this expression into the free energy (3) we obtain the simple result

$$F(\xi )=F_c(\xi )+F_{el}(\xi )=f\frac{A^2}{C}\xi ^2+q^2\pi K\mathrm{ln}\frac{R_{\mathrm{}}}{\xi }$$
(6)
where $`f0.97`$ and $`R_{\mathrm{}}`$ is the size of the system. Differentiation with respect to $`\xi `$ yields $`\xi =|q|\zeta `$, where $`\zeta =\sqrt{(2\pi /f)(L/A)}`$ is the nematic correlation length. If $`\zeta `$ is much smaller than the other characteristic lengths of the system, we can use the Frank elastic free energy (4) to obtain the director configuration in systems with topological defects. The Frank elastic free energy, however, is defined outside the core region only. The total free energy can be computed by additing to it the core energy of each defect

$$F_c=q^2\frac{\pi K}{2}.$$
(7)
## III Isolated disk

Fukuda and Yokoyama  studied the distortion of a 2D nematic around an isolated disk, with homeotropic anchoring, by minimizing the Landau-de Gennes free energy (2) using numerical methods with adaptive grids. They have shown that a director configuration with a dipolar defect (topological charge $`q=1`$) is unstable, and changes into a configuration with a pair of defects with topological charge $`q=1/2`$ that are located symmetrically around the particle.

The energy of this quadrupolar configuration may be calculated analitically and compared with the energy of the dipolar case, obtained by Pettey et. al.  for 2D smectic-C films.

Consider the disk at the center of the reference frame $`𝐫=(x,y)`$, or $`𝐫=(r,\varphi )`$ in polar coordinates, where $`\varphi `$ is defined counterclockwise, starting from the $`x`$-axis, in the range $`\pi <\varphi <\pi `$ (see Fig.1). At the outer boundary, for large $`r`$, the nematic is uniform with $`\theta =\pi /2`$. Homeotropic boundary conditions apply at the disk’s boundary, i.e., $`\theta (\varphi )=\varphi `$ on the upper contour and $`\theta (\varphi )=\pi +\varphi `$ on the lower one. Minimization of the Frank elastic free energy (4) yields Laplace’s equation $`\mathrm{\Delta }\theta =0`$, everywhere (in the nematic) except within the defects.

An image solution with two half integer topological charges located symmetrically at a distance $`r_d`$ from the center of the disk, satisfying the homeotropic boundary conditions at the disk’s contour and a uniform director at infinity, is given by

$$\theta =\frac{\pi }{2}2\psi +\frac{1}{2}\psi _{r_d}+\frac{1}{2}\psi _{r_d}+\frac{1}{2}\psi _{\stackrel{~}{r}_d}+\frac{1}{2}\psi _{\stackrel{~}{r}_d}$$
(8)
where $`\stackrel{~}{r}_d=a^2/r_d`$ is the position of a virtual half integer topological defect inside the disk. The field

$$\psi _{x_i}(x,y)=\mathrm{arctan}\frac{xx_i}{y}$$
(9)
describes a topological defect of charge $`q=1`$ located at $`𝐫_i=(x_i,0)`$ and satisfies Laplace’s equation everywhere except at $`𝐫_i`$.

The elastic Frank free energy can be calculated using Green’s Theorem:

$$_\mathrm{\Omega }(\theta )^2d^2𝐫=_\mathrm{\Omega }\theta \frac{\theta }{𝐧}𝑑𝐥_\mathrm{\Omega }\theta \mathrm{\Delta }\theta d^2𝐫,$$
(10)
where $`\mathrm{\Omega }`$ is the boundary of the 2D nematic domain, $`d𝐥=(dx,dy)`$ is an oriented contour element, and $`\theta /𝐧=(\theta /y,\theta /x)`$. If the tilt angle $`\theta `$ satisfies Laplace’s equation, the second term on the right hand side of (10) vanishes, and the free energy becomes

$$F_{el}=\frac{K}{2}_\mathrm{\Omega }\theta \frac{\theta }{𝐧}𝑑𝐥.$$
(11)
Close to the outer boundary, for $`ra`$, the director is nearly uniform and the contribution of the integral over this contour may be neglected. The contribution of the integral over the upper disk’s contour is calculated by taking $`\gamma _1(t)=(x(t),y(t))=(a\mathrm{cos}t,a\mathrm{sin}t)`$ with $`\theta (t)=t`$, yielding

$`{\displaystyle _{\gamma _1}}\theta {\displaystyle \frac{\theta }{𝐧}}𝑑𝐥=\pi \left[\mathrm{ln}\left(1{\displaystyle \frac{a}{r_d}}\right)+\mathrm{ln}\left(1+{\displaystyle \frac{a}{r_d}}\right)\right]`$ (12)
By symmetry the integral over the lower disk’s contour is identical.

Now consider the integrals over the straight lines along the $`x`$-axis from the disk to the defects. On $`\gamma _2`$ and $`\gamma _6`$, the tilt angle is zero and there is no contribution to the free energy. On $`\gamma _4(t)=(t,0)`$ (and $`\gamma _8`$), $`\theta (t)=\pi `$ and a simple calculation yields

$`{\displaystyle _{\gamma _4}}\theta {\displaystyle \frac{\theta }{𝐧}}𝑑𝐥`$ $`2\pi (2\mathrm{ln}{\displaystyle \frac{r_d}{a}}{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{\xi }{r_da}}{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{2r_d}{r_d+a}}`$ (13)
$`{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{r_d+a}{a}}{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{r_d^2+a^2}{a(r_d+a)}}).`$ (14)
in the limit $`r_d>>\xi `$. If $`a>>\xi `$ we may neglect the contributions from the circular contours around the defects. Finally, the Frank elastic free energy is given by

$$F_{el}(r_d)\frac{\pi K}{2}\left[\mathrm{ln}\frac{a}{\xi }\mathrm{ln}2\left(\frac{a^3}{r_d^3}\frac{a^7}{r_d^7}\right)\right].$$
(15)
Minimizing with respect to $`r_d`$ yields $`(r_d)_{eq}=\sqrt[4]{7/3}a`$, as obtained by Fukuda and Yokoyama  using a simple force balance argument. The total free energy for this quadrupolar configuration is obtained by adding the core energy of the defects $`2F_c=\pi K/4`$, and using the appropriate core’s size $`\xi =\zeta /2`$.

The Frank elastic free energy for a dipolar configuration, with a topological defect of charge $`q=1`$, with the preferred position at $`𝐫_d=(0,r_d)`$ was calculated by Pettey et. al. :

$$F_{el}(r_d)\pi K\left[\mathrm{ln}\frac{a}{\xi }\mathrm{ln}\left(\frac{a^2}{r_d^2}\frac{a^4}{r_d^4}\right)\right]$$
(16)
yielding the equilibrium position $`(r_d)_{eq}=\sqrt{2}a`$. The core energy of the defect is now $`F_c=\pi K/2`$, and the core’s size is $`\xi =\zeta `$.

The difference between the free energy of these configurations, when the core energies of the defects are taken into account, is a function of the correlation length $`\zeta `$:

$`F_{dip}F_{quad}`$ $``$ $`{\displaystyle \frac{\pi K}{2}}\left(2.08+\mathrm{ln}{\displaystyle \frac{a}{\zeta }}\right).`$ (17)
For $`a>\zeta `$, the difference is positive, and thus the quadrupolar configuration has the lowest free energy.

This result contrasts with that for spherical colloidal particles in 3D nematic liquid crystals. In the latter case, the dipolar configuration has a punctual defect (with zero core energy). The quadrupolar configuration exhibits a disclination ring, and its elastic free energy depends logarithmically on the ratio between the correlation length and the particle’s radius. For small particles, the nematic distortion is of the quadrupolar type. As we increase the ratio $`a/\zeta `$, the energy of the quadrupolar configuration increases, and eventually exceeds that of the dipolar configuration. This crossover occurs at $`a/\zeta 10^3`$, that corresponds to particles with radii of the order of $`\mu `$m.

## IV Long range interaction

If the separation $`𝐑=(R\mathrm{cos}\alpha ,R\mathrm{sin}\alpha )`$ between particles is large, the nematic distortion will be approximately the sum of the isolated quadrupolar solutions, $`\theta \theta _1+\theta _2`$, with an interaction energy

$`F_{int}`$ $``$ $`K{\displaystyle _\mathrm{\Omega }}\theta _1\theta _2d^2𝐫`$ (18)
$`=`$ $`K{\displaystyle _{\mathrm{\Gamma }_1}}\theta _1{\displaystyle \frac{\theta _2}{𝐧}}𝑑𝐥+K{\displaystyle _{\mathrm{\Gamma }_2}}\theta _2{\displaystyle \frac{\theta _1}{𝐧}}𝑑𝐥`$ (19)
where the subscripts refer to the different disks (see Fig.2).

Consider the integral over the contour $`\mathrm{\Gamma }_2`$ around the second disk. We use the same notation as for isolated disks and neglect the contributions from the circular contours around the defects. To calculate the contribution from the upper circular contour around disk 2 we write $`\gamma _1(t)=(R\mathrm{cos}\alpha +a\mathrm{cos}t,R\mathrm{sin}\alpha +a\mathrm{sin}t)`$. On this contour $`\theta _2(t)=t`$, and $`\theta _1`$ can be expanded as a power series in $`a/R`$ yielding

$`{\displaystyle _{\gamma _1}}\theta _2{\displaystyle \frac{\theta _1}{𝐧}}𝑑𝐥`$ $`=`$ $`{\displaystyle _\pi ^0}t({\displaystyle \frac{\theta _1}{x}}|_𝐑a\mathrm{cos}t+{\displaystyle \frac{\theta _1}{y}}|_𝐑a\mathrm{sin}t`$ (20)
$`+`$ $`{\displaystyle \frac{^2\theta _1}{x^2}}|_𝐑a^2\mathrm{cos}^2t+{\displaystyle \frac{^2\theta _1}{y^2}}|_𝐑a^2\mathrm{sin}^2t`$ (21)
$`+`$ $`2{\displaystyle \frac{^2\theta _1}{xy}}|_𝐑a^2\mathrm{sin}t\mathrm{cos}t)dt+\mathrm{}`$ (22)
When the contribution from the lower circular contour around the disk is added, some of the terms in (22) cancel. After integration, the sum of these two contributions gives

$$_{\gamma _1+\gamma _5}\theta _2\frac{\theta _1}{𝐧}𝑑𝐥=\pi a^2\frac{^2\theta _1}{xy}|_𝐑+\mathrm{}$$
(23)
The integrals over the straight lines can be computed in a similar manner. On $`\gamma _2`$ and $`\gamma _6`$ the tilt angle $`\theta _2=0`$, and there is no contribution to the free energy. On $`\gamma _4(t)=(R\mathrm{cos}\alpha +t,R\mathrm{sin}\alpha )`$ (or $`\gamma _8`$), $`\theta _2(t)=\pi `$ and $`\theta _1`$ may be expanded as a power series in $`a/R`$. The result is

$$_{\gamma _{rl}}\theta _2\frac{\theta _1}{𝐧}𝑑𝐥=\pi (r_d^2a^2)\frac{^2\theta _1}{xy}|_𝐑+\mathrm{}$$
(24)
The derivative $`^2\theta _1/xy`$ is obtained from the quadrupolar expression

$$\theta _1=\frac{\pi }{2}(r_d^2+\stackrel{~}{r}_d^2)\frac{xy}{(x^2+y^2)^2}+\mathrm{}$$
(25)
valid for large $`r=\sqrt{x^2+y^2}`$. Finally, adding the symmetric contribution from the contour $`\mathrm{\Gamma }_1`$, the total interaction energy becomes

$$F_{int}6\pi K(r_d^4+a^4)\frac{12\mathrm{sin}^22\alpha }{R^4}.$$
(26)
The effective interaction between disks decays as $`R^4`$ and is strongly anisotropic: repulsive if the particles are aligned horizontally or vertically ($`\alpha =0`$ or $`\alpha =\pi /2`$, respectively), and attractive for intermediate oblique orientations (the preferred orientation is $`\alpha =\pi /4`$).

## V Short range interaction

At small disk separations, the nematic deformation is no longer the sum of the isolated quadrupolar solutions (26).

The interaction between two circular disks was studied numerically, using finite elements with adaptive meshes to minimize the total nematic free energy at a given separation. A first triangulation respecting the predefined physical geometry (disks separation and orientation) is constructed . The triangulation space corresponds to the domain of integration of the free energy. It is delimited by the disks’ boundaries and an outer boundary at infinity, in practice a large distance from the disks, of the order of $`20a`$. Although analytically simpler, the Frank free energy (4) is plagued with divergences of $`\theta `$ at the defects, and thus it is not the most adequate for numerical computations. For this reason, we used the Landau-de Gennes free energy (2) for a 2D nematic. The functions $`Q_{xx}`$ and $`Q_{xy}`$ are set at all vertices of the mesh and are linearly interpolated within each triangle. Using standard numerical procedures the elastic free energy is minimized under the constraints imposed by the boundary conditions, i.e., homeotropic strong anchoring at the disks perimeter and constant alignment at the outer boundary. Finally, to overcome the differences in length scales set by the disks and the defects, of up to 2 orders of magnitude (we chose $`\zeta ^2L/A=10^3`$), new adapted meshes are generated iteratively from the result of the previous minimization. The new local triangle sizes are related with the variations of the previous solution. If the variations are strong, finer meshes are required in order to guarantee an almost constant numerical weight of each minimization variable . The final meshes with $`2\times 10^4`$ variables and a minimal edge length of $`10^4a`$ provided excellent accuracy, yielding free energies with a relative error of the order of $`10^4`$.

To test the numerical accuracy of the procedure we calculated the nematic configuration for a single isolated disk, and found a solution with two half integer defects at $`r_d=1.23a\pm 0.01a`$, in very good agreement with the analytical result.

The effective interaction between two disks is plotted in Fig.3 as a function of the distance $`R`$ between the disks, at three different orientations ($`\alpha =0,\pi /4,\pi /2`$). A simple fit of these curves confirmed the quadrupolar long range decay $`R^4`$, at large disk separations. However, for smaller separations, the free energy changes dramatically. In particular, orientations that are repulsive at large separations ($`\alpha =0,\pi /2`$) become attractive at different threshold distances ($`R_{th}5a`$ and $`R_{th}3a`$, respectively). Also, the free energy corresponding to the parallel orientation ($`\alpha =0`$) becomes smaller than the free energy of the long range preferred oblique orientation ($`\alpha =\pi /4`$) at short distances.

To elucidate this behavior we plotted on the left of Fig.4 the elastic free energy as a function of the orientation $`\alpha `$, at three different separations ($`R=2.4a,3.0a,4.0a`$). The minimum of each curve defines a preferred orientation $`\alpha ^{}`$, that decreases as the disks approach each other. On the right of Fig.4 we see that at very large separations, the disks have an oblique orientation (with $`\alpha ^{}=\pi /4`$), that changes to a parallel orientation (with $`\alpha ^{}=0`$) as their separation decreases.

In the following figures we show the nematic order parameter for several separations, when the orientation of the disks is parallel ($`\alpha =0`$, Fig.5) or perpendicular ($`\alpha =\pi /2`$, Fig.6). The nematic order parameter varies between $`Q=0`$ (grey regions) and $`Q=Q_{eq}`$ (white). These regions around the disks correspond to the nematic defects.

In both disk alignments, if the separation is large enough, the half integer defects around each disk stay at the same positions as in the isolated disk configuration. However, at a certain separation $`R_{th}`$, exactly when the interaction becomes attractive for that orientation, the position of the defects starts to change, in order to minimize the total free energy.

Finally, we analyzed in detail the interaction between the disks for very small separations. For strong (fixed) homeotropic anchoring, we observed a repulsion when the disks are nearly at contact (for all orientations) at a separation $`R2.1a`$. This repulsion may not occur if we do not impose fixed homeotropic boundary conditions at the disks. If we consider a finite anchoring strength we have to add to the free energy (2) a surface term

$$F_a=\frac{W}{2}_{\mathrm{\Lambda }_1+\mathrm{\Lambda }_2}(Q_{ij}Q_{ij}^n)^2𝑑l$$
(27)
where $`Q_{ij}^n`$ are the values of the tensor order parameter for homeotropic alignement.

In Fig.7 we plot the total elastic free energy as a function of the distance $`R`$ between disks, for three different reduced anchoring strengths ($`W/Ka=1000;40;30`$). These anchoring strengths are large enough to induce defects in the nematic host. Clearly the repulsion appears at small separations for strong anchoring ($`W/Ka=1000`$) and vanishes at a critical anchoring strength that lies between $`30<W/Ka<40`$. In the weak anchoring regime the coalescence of the dropplets may occur.

## VI Conclusions

We compared the (analytic) elastic free energies of 2D nematic configurations with one (satellite) and two (ring) defects around isolated disks with homeotropic boundary conditions, and found that the ring configuration is the stable one for realistic values of the reduced correlation length $`\zeta /a`$.

Using the analytic ring solution we derived the long range interaction between two disks immersed in 2D nematic hosts. We found an anisotropic interaction, that is repulsive for orientations close to perpendicular or parallel orientations and attractive for oblique alignement, with a preferred orientation at $`\alpha =\pi /4`$. In this limit, the interaction decays as $`R^4`$.

Finally, the analytic results were compared with very accurate numerical calculations, valid at arbitrary separation. At small disk separations, the defects may change their positions, leading to an unexpected complex interaction with the long range repulsive orientations becoming attractive. As a result as we decrease the distance between the disks, their preferred orientation changes from oblique to parallel alignment. For sufficiently large anchoring strengths, at very short range, we found a new repulsion for all orientations, that may prevent coalescence. This repulsion vanishes for small anchoring strengths.

The equilibrium structure of inverted nematic emulsions is complex and depends sensitively on the radii of the disks, the nematic correlation length and anchoring strength. In particular short range nonlinear effects may change the interaction between disks dramatically. It is thus reasonable to assume that three body interactions may be relevant to describe the formation of the the structures observed experimentally. Finally, in order to understand the observed structures dynamical aspects may also have to be taken into account.

Acknowledgments

We acknowledge the support of the Fundação para a Ciência e Tecnologia (FCT) through a running grant (Programa Plurianual) and grants No. SFRH/BPD/1599/2000 (MT) and No. SFRH/BPD/5664/2001 (PP).


format: text

text:
# Products of non-stationary random matrices and Multiperiodic equations of several scaling factors

## 1. Introduction

Kingman’s subadditive ergodic theorem was originally proved in 1968 \[Ki1, Ki2\]. A more recent proof was given by Katznelson and Weiss in 1982 \[KW\]. It is one of the most important results in ergodic theory. In this paper we consider the following set-up which resembles a dynamical system without invariant measure and try to get results similar to Kingman’s theorem. Let $`\beta >1`$ be a positive real number. Let $`\{f_n\}`$ be a sequence of uniformly almost periodic functions (i.e. in the sense of Bohr, see Section 2.1) defined on the real line $``$. Suppose the following subadditivity condition is fulfilled

$$f_{n+m}(x)f_n(x)+f_m(\beta ^nx)\text{for a.e.}x\text{and all}n,m.$$

where $`a.e.`$ refers to the Lebesgue measure. We would like to study the almost everywhere convergence of $`n^1f_n(x)`$. The Kingman theorem applies in the special case where $`\beta >1`$ is an integer and the $`f_n`$’s are periodic. The typical case in our mind is

$$f_n(x)=\mathrm{log}M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x)$$
(1.1)
where $`M:\text{GL}(^d)`$ is a matrix-valued uniformly almost periodic function. We will prove that the limit $`lim_n\mathrm{}n^1f_n(x)`$ exists almost everywhere (a.e. for short) with respect to the Lebesgue measure under the condition that the $`n^1f_n(x)`$ have *joint periods* (see Theorem 2.5). As a consequence, an Oseledec type theorem is proved for the matrix products involved in (1.1) (see Theorem 2.9). It is proved that the condition on the existence of joint periods is satisfied when $`\beta `$ is a PV-number (see Section 3).

Our consideration is partially motivated by the study of multiperiodic functions, already investigated by Strichartz et al. \[JRS\], Fan and Lau \[FL\], and Fan \[F\]. By a Multiperiodic function of one real variable we mean any function $`F:`$ which is a solution of a functional equation of the following form

$$F(\xi )=f_1\left(\frac{\xi }{\rho _1}\right)F\left(\frac{\xi }{\rho _1}\right)+\mathrm{}+f_d\left(\frac{\xi }{\rho _d}\right)F\left(\frac{\xi }{\rho _d}\right)$$

where $`d1`$ is an integer; $`\rho _1>1,\mathrm{},\rho _d>1`$ are $`d`$ real numbers, called scaling factors; $`f_1,\mathrm{},f_d`$ are $`d`$ complex valued functions defined on the real line, called determining functions. The equation will be called a multiperiodic equation.

We will assume that the determining functions $`f_j`$ are periodic or almost periodic in the sense of Bohr, as is the case in most applications. We will also assume that the scaling factors $`\rho _j`$ are commensurable in the sense that $`\rho _j`$ are powers of some real number $`\beta >1`$. Without lost of generality, we assume that $`\rho _j=\beta ^j`$ for $`1jd`$. Then the multiperiodic equation becomes

$$F(\xi )=f_1\left(\frac{\xi }{\beta }\right)F\left(\frac{\xi }{\beta }\right)+\mathrm{}+f_d\left(\frac{\xi }{\beta ^d}\right)F\left(\frac{\xi }{\beta ^d}\right)$$
(1.2)
As far as we know, there is few work done for the non–commensurable case which is much more difficult.

In the literature, the case where $`d=1`$ and $`\beta =2`$ (or an arbitrary integer) has been studied, especially in the theory of wavelets \[D\]. In fact, the scaling function $`\phi `$ of a wavelet satisfies a scaling equation

$$\phi (x)=a_n\phi (2xn).$$

The Fourier transform of $`\phi `$ satisfies a multiperiodic equation of the form (1.2) with only one scaling factor $`\beta =2`$ and only one determining function $`f_1(x)=f(x)=\frac{1}{2}_na_ne^{inx}`$.

The scaling functions in wavelets constitute a class of functions sharing a kind of similarity. More generally, multiperiodic functions arise as Fourier transforms of self-similar objects such as Bernoulli convolution measures ($`d=1`$, $`\beta >1`$ being a real number and $`f`$ being a trigonometric polynomial), inhomogeneous Cantor measures ($`d`$ may be greater than $`1`$) or more general self-similar measures produced by iterated function systems (see \[S\]).

In the case of one scaling factor (i.e. $`d=1`$, then we write $`f_1(x)=f(x)`$), the existence of the solution of the multiperiodic equation (1.2) is simple and is assured by the consistency condition $`f(0)=1`$ and a regularity condition, say $`f`$ is Lipschitz continuous. Actually the solution can be written as an infinite product

$$F(x)=\underset{n=1}{\overset{\mathrm{}}{}}f\left(\frac{x}{\beta ^n}\right).$$

For the existence of the general equation (1.2), we have

Theorem A. Let $`d1`$. Suppose that the determining functions $`f_1,\mathrm{},f_d`$ are Lipschitz continuous, and satisfy the consistency condition

$$f_1(0)+\mathrm{}+f_d(0)=1.$$

Suppose furthermore that $`f_1(0),\mathrm{},f_d(0)[0,+\mathrm{})`$. Then equation (1.2) admits a unique continuous solution $`F`$ such that $`F(0)=1`$.

The proof of this theorem is postponed to Section 4.2.

Our study of equation (1.2) is converted to that of vector valued equations of the form

$$G(x)=M\left(\frac{x}{\beta }\right)G\left(\frac{x}{\beta }\right)$$
(1.3)
where $`M:_{d\times d}()`$ is a matrix valued determining function and $`G:^d`$ is a vector valued unknown function. Matrix products will be involved in the study of equation (1.3), which produces some difficulties. However, equation (1.3) is a simple recursive relation because it contains only one scaling factor. Equation (1.2) is equivalent to equation (1.3) with $`M(x)`$ and $`G(x)`$ equal respectively to

$$\left(\begin{array}{ccccc}f_1(x)& f_2(\frac{x}{\beta })& \mathrm{}& f_{d1}(\frac{x}{\beta ^{d2}})& f_d(\frac{x}{\beta ^{d1}})\\ 1& 0& \mathrm{}& 0& 0\\ 0& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& & 0& 0\\ 0& \mathrm{}& 0& 1& 0\end{array}\right)\text{and}\left(\begin{array}{c}F(x)\\ F(\frac{x}{\beta })\\ \mathrm{}\\ F(\frac{x}{\beta ^{d1}})\end{array}\right)$$
(1.4)
We would like to know the asymptotic behavior at infinity of the solution $`G`$. This is a natural question because $`G`$ often represents Fourier transform of a function (a measure or a distribution) and the asymptotic behavior at infinity describes quantitatively the regularity of the solution. Unfortunately, there is no closed form formula for $`G`$ in general and the behavior of $`G`$ is rather complicated, as is shown by the Fourier transform of the Cantor measure ($`d=1`$, $`\beta =3`$ and $`f(\xi )=\mathrm{cos}\xi `$).

Following \[JRS\], we will study the pointwise asymptotic behaviors of

$$h_n(x):=\frac{1}{n}\mathrm{log}|F(\beta ^nx)|$$

as $`n\mathrm{}`$.

We will prove that, under some conditions, the limit $`lim_n\mathrm{}h_n(x)`$ exists and is equal to a constant almost everywhere with respect to Lebesgue measure.

This answers partially a questions in \[JRS\]. More precisely, we have the following results, whose proofs are postponed to Section 4.3.

Theorem B. Suppose that the conditions in Theorem A are satisfied. Furthermore, suppose that the determining functions $`f_1,\mathrm{},f_d`$ are either identically zero or strictly positive and $`1`$-periodic, and that $`\beta >1`$ is a Pisot number. Let $`F`$ be the solution of the equation (1.2). Then there is a constant $``$ such that

$$\underset{n\mathrm{}}{lim}\frac{1}{n}F(\beta ^nx)=a.e.$$

The constant $``$ in the theorem is the leading Liapunov exponent of the matrix $`M(x)`$ above, defined by

$$(M)=\underset{n1}{inf}\frac{1}{n}𝕄\mathrm{log}M(\beta ^{n1}x)M(\beta ^{n2}x)\mathrm{}M(\beta x)M(x).$$

where $`𝕄f`$ denotes the Bohr mean of an almost periodic function $`f`$ (see Section 2.1 below).

Theorem C. Suppose that the conditions in Theorem A are satisfied. Furthermore, suppose that the determining functions $`f_1,\mathrm{},f_d`$ are $`1`$-periodic, Lipschitz, and that $`\beta >1`$ is a Pisot number with maximal conjugate of modulus $`\rho `$. Let $`F`$ be the solution of the equation (1.2). If

$$\underset{x}{sup}\frac{(1+|f_1(x)|+\mathrm{}+|f_{d1}(x/\beta ^{d2})|)(|f_1(x)|+\mathrm{}+|f_d(x/\beta ^{d1})|)}{|f_d(x/\beta ^{d1})|}<\rho ^1$$

then there is a constant $`\lambda `$ such that

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}\underset{j=0}{\overset{d1}{}}|F(\beta ^{nj}x)|=\lambda \text{a.e.}$$

In this case we do not know if the constant $`\lambda `$ equals the leading Lyapunov exponent of the matrix.

Theorem A will be proved in Section 4.2 as a special case of a more general result (Theorem 4.1), Theorm B and Theorem C in Section 4.3. Both Theorm B and Theorem C are consequences of our Kingman’s Theorem (Theorem 2.5) and Oseledec’s Theorem (Theorem 2.9) which are discussed in Section 2. In Section 3, we prove that the joint period condition required in both Kingman’s Theorem and Oseledec Theorem is satisfied when $`\beta >1`$ is a Pisot number.

## 2. Kingman theorem and Oseledec theorem

### 2.1. Total Bohr ergodicity and joint $`ϵ`$-period

Let us first recall the definition of *uniformly almost periodic* functions and some of their properties (see \[Bo\]). Next we will introduce the notions of total Bohr ergodicity and a joint $`ϵ`$-period.

Let $`f`$ be a real or complex valued function defined on the real line. A number $`\tau `$ is called a translation number of $`f`$ belonging to $`ϵ0`$ (or an $`ϵ`$-period) if

$$\underset{x}{sup}|f(x+\tau )f(x)|ϵ.$$

We say that $`f`$ is a uniformly almost periodic (u.a.p.) function if it is continuous and if for any $`ϵ>0`$ the set of its translation numbers belonging to $`ϵ`$ is relatively dense (i.e. there exists a number $`\mathrm{}>0`$ such that any interval of length $`\mathrm{}`$ contains at least one such translation number). H. Bohr proved that the space of all u.a.p. functions is a closed sub-algebra of the Banach algebra $`C_b()`$ of bounded continuous functions equipped with the uniform norm and that it is the closure of the space of all (generalized) trigonometric polynomials of the form

$$\underset{\mathrm{finite}}{}A_ne^{i\mathrm{\Lambda }_nx}(A_n,\mathrm{\Lambda }_n).$$

For any u.a.p. function $`f`$, as is proved by Bohr, the following limit exists

$$𝕄f=\underset{T\mathrm{}}{lim}\frac{1}{T}_0^Tf(x)𝑑x.$$

It is called the Bohr mean of $`f`$. For any locally integrable not necessarily u.a.p. $`f`$, we define $`𝕄f`$ as the limsup instead of the limit.

###### Definition 2.1.

A sequence of real numbers $`(u_n)_{n0}`$ is said to be totally Bohr ergodic if for any arithmetic subsequence $`(u_{am+b})_{m0}`$ ($`a1,b0`$ being fixed) and for any real $`p>0`$, the sequence $`(u_{am+b}x)_{m0}`$ is uniformly distributed (modulo $`p`$) for almost every $`x`$ with respect to the Lebesgue measure.

The following is the main property of totally Bohr ergodic sequences that we will use.

###### Lemma 2.2.

Suppose that $`(u_n)_{n0}`$ is a totally Bohr ergodic sequence. Then for any u.a.p. function $`f`$ and any integers $`a1,b0`$, we have

$$\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{n=0}{\overset{N}{}}f(u_{an+b}x)=𝕄fa.e.$$

###### Proof.

It is a consequence of the fact that $`u_{an+b}x`$ is uniformly distributed ($`modp`$) for almost all $`x`$, for any real $`p>0`$, the fact that $`f`$ can be uniformly approximated by trigonometric polynomials and the Weyl criterion. ∎

###### Remark 2.3.

Suppose the $`(u_n)_{n0}`$ is a sequence such that $`inf_{nm}|u_nu_m|>0`$, then the sequence $`(u_nx)`$ is uniformly distributed for almost every point $`x`$ \[KL\]. Consequently, the sequence $`(u_n)`$ is totally Bohr ergodic. A more special case is $`u_n=\beta ^n`$ with $`\beta >1`$. This is the most interesting case for us. On the other hand, no bounded sequence can be totally Bohr ergodic.

###### Definition 2.4.

Let $`(F_n)_{n0}`$ be a sequence of u.a.p. functions. Let $`ϵ>0`$ and $`N`$. A real number $`\tau `$ is called a joint $`ϵ`$-translation number for $`(F_n)_{nN}`$ if

$$\underset{nN}{sup}\underset{x}{sup}|F_n(x+\tau )F_n(x)|ϵ.$$

If for any $`ϵ>0`$ there exists $`N(ϵ)`$ such that such joint $`ϵ`$-translation numbers for $`(F_n)_{nN(ϵ)}`$ are relatively dense, we say that $`(F_n)_{n0}`$ has joint periods.

### 2.2. Kingman’s theorem

Following ideas of Katznelson and Weiss \[KW\] we prove the following version of Kingman’s Theorem. The difficulty in our case is that we have to deal with an infinite measure space. We are also dealing with non–stationary sequences.

###### Theorem 2.5.

Let $`(u_n)_{n0}`$ be a totally Bohr ergodic sequence of real numbers and $`(f_n)_{n0}`$ be a sequence of uniformly almost periodic functions. Suppose
(i) The sequence $`(n^1f_n)`$ has joint periods.
(ii) The following subadditivity is fulfilled

$$f_{n+m}(x)f_n(x)+f_m(u_nx)\text{for a.e. }x\text{ and all }n,m.$$

(iii) For any $`n1`$

$$\underset{m}{sup}\left(f_m(u_nx)f_m(x)\right)<\mathrm{}\text{for a.e. }x.$$
(2.1)
Then the following limit exists and is a constant

$$\underset{n\mathrm{}}{lim}\frac{1}{n}f_n(x)=\underset{n}{inf}\frac{1}{n}𝕄f_n\text{for a.e. }x.$$

###### Proof.

The proof is a modification of Katznelson-Weiss’ proof \[KW\]. Without loss of generality we assume that $`u_0=1`$.

Let $`\gamma =inf_n\frac{1}{n}𝕄f_n`$. Let us put

$$f^{}(x)=\underset{n\mathrm{}}{lim\; inf}\frac{1}{n}f_n(x),f^+(x)=\underset{n\mathrm{}}{lim\; sup}\frac{1}{n}f_n(x).$$

We remark that the subadditivity implies that $`f^\pm (x)f^\pm (u_nx)`$ for all $`n`$ and a.e. $`x`$. In a finite measure space this would imply the invariance a.e. In our case it is the boundedness (2.1) which implies that $`f^\pm (u_nx)f^\pm (x)`$ a.e., what makes the functions $`f^\pm `$ invariant in the sense that $`f^\pm (x)=f^\pm (u_nx)`$ a.e. ($`n`$).

The first part of the proof, i.e. $`f^+\gamma `$ a.e., is simple. We just exploit the fact that any (infinite) arithmetical subsequence of $`(u_nx)`$ is Bohr-uniform distributed. This provides us with a kind of ergodic theorem. In fact, fix an integer $`N`$. For any integer $`n`$ write $`n=mN+r`$ with $`0r<N`$. We have

$$f_n(x)\underset{k=0}{\overset{m1}{}}f_N(u_{kN}x)+f_r(u_{mN}x).$$

The $`N`$ functions $`f_r`$ ($`r=0,1,\mathrm{},N1`$) being bounded, by Lemma 2.2, this readily implies by the Bohr-uniform distribution of $`(u_{kN}x)_{k0}`$ that

$$f^+(x)\underset{m\mathrm{}}{lim}\frac{1}{mN+r}\underset{k=0}{\overset{m1}{}}f_N(u_{kN}(x))=\frac{1}{N}𝕄f_N$$

for a.e. $`x`$. Hence $`f^+(x)\gamma `$ for a.e. $`x`$.

Next, we want to prove that $`f^{}(x)\gamma `$ for a.e. $`x`$. For this we assume that $`\gamma >\mathrm{}`$, otherwise it is trivially true. Adding to each $`f_n`$ the constant value $`nf_1`$ creates a new subadditive sequence $`\stackrel{~}{f}_n:=f_nnf_1`$ with $`\stackrel{~}{f}_n0`$ ($`n1`$), $`f^{}=\stackrel{~}{f}^{}+f_1`$ and $`\gamma _f=\gamma _{\stackrel{~}{f}}+f_1`$. So, we may assume that $`f_n0`$ ($`n1`$). We can furthermore set $`f_1=0`$ (observe that this will not affect the subbadditivity condition since $`f_n0`$). Then for any $`\mathrm{\Delta }>0`$ we truncate the function $`f_n`$ in the way

$$f_{n,\mathrm{\Delta }}=\mathrm{max}(f_n,n\mathrm{\Delta }).$$

Note that the sequence $`f_{n,\mathrm{\Delta }}`$ fulfills the assumptions of the theorem. Note that in this case $`\gamma _\mathrm{\Delta }\mathrm{\Delta }`$ and also $`f_\mathrm{\Delta }^{}(x)\mathrm{\Delta }`$ for all $`x`$. It is clear that

$$f_\mathrm{\Delta }^\pm (x)=\mathrm{max}(f^\pm (x),\mathrm{\Delta }).$$

Assume that we proved the theorem for the sequence $`f_{n,\mathrm{\Delta }}`$ for any $`\mathrm{\Delta }`$. Then we claim that $`f_\mathrm{\Delta }^{}(x)f^{}(x)`$ for a.e. $`x`$ as $`\mathrm{\Delta }`$ goes to $`\mathrm{}`$. In fact, if $`f_\mathrm{\Delta }^{}(x)>\mathrm{\Delta }`$ for some $`\mathrm{\Delta }`$ then $`f_\mathrm{\Delta }^{}(x)=f^{}(x)`$ and $`\gamma _\mathrm{\Delta }=\gamma _f>\mathrm{}`$. On the other hand if $`f_\mathrm{\Delta }^{}(x)=\mathrm{\Delta }`$ for all $`\mathrm{\Delta }`$ then obviously $`f^{}(x)=\mathrm{}`$ and $`\gamma _\mathrm{\Delta }=\gamma _f=\mathrm{}`$. This proves the theorem for the sequence $`f_n`$.

¿From now on we assume that $`f_1=0`$, $`f_n0`$ and $`f_n`$ is truncated and we skip the subscript $`\mathrm{\Delta }`$. Let $`ϵ>0`$. By the hypothesis (i) on the joint periods, there is an integer $`N(ϵ)`$ such that the joint $`ϵ`$-translation numbers are relatively dense. For these numbers $`\tau `$ we have

$$\left|\frac{f_n(x+\tau )}{n}\frac{f_n(x)}{n}\right|ϵ(x,nN(ϵ))$$
(2.2)
Notice that there is no loss of generality to suppose that $`N(ϵ)`$ increases as $`ϵ`$ decreases to $`0`$. We define

$$n_ϵ(x)=\mathrm{min}\left\{nN(ϵ):\frac{1}{n}f_n(x)f^{}(x)+ϵ\right\}$$

Let $`A_K^ϵ=\{x:n_ϵ(x)>K\}`$. Notice that if $`ϵ^{}<ϵ^{\prime \prime }`$, we have $`N(ϵ^{})N(ϵ^{\prime \prime })`$ and $`n_ϵ^{}(x)n_{ϵ^{\prime \prime }}(x)`$, so that $`A_K^{ϵ^{\prime \prime }}A_K^ϵ^{}`$ for $`K>N(ϵ^{})`$.

We claim that

###### Lemma 2.6.

For any $`ϵ>0`$, we have

$$\underset{K\mathrm{}}{lim}\widehat{𝕄}(A_K^ϵ)=0$$

where $`\widehat{𝕄}A=𝕄1_A`$ denotes the Bohr mean of the characteristic function of the set $`A`$ (defined if necessary with the limsup).

In order to prove this Lemma 2.6 we need the following lemma which says that $`A_K^ϵ`$ is to some extent periodic.

###### Lemma 2.7.

For any joint $`ϵ`$-translation number $`\tau `$ of $`(f_n/n)_{nN(ϵ)}`$, we have $`A_K^{2ϵ}+\tau A_K^ϵ`$.

Let us first prove Lemma 2.7. Suppose $`xA_K^{2ϵ}+\tau `$, i.e. $`x\tau A_K^{2ϵ}`$, then

$$\frac{f_n(x\tau )}{n}>f^{}(x\tau )+2ϵ(N(2ϵ)nK).$$

This, together with the fact that $`\tau `$ is a joint $`ϵ`$-translation number for all $`f_n/n`$ with $`nN(ϵ)(N(2ϵ))`$ (see (2.2)), implies

$$\frac{f_n(x)}{n}\frac{f_n(x\tau )}{n}ϵ>f^{}(x\tau )+ϵ(N(ϵ)nK).$$

That means $`xA_K^ϵ`$. Thus we have finished the proof of Lemma 2.7.

Now let us prove Lemma 2.6. Since joint $`\frac{ϵ}{2}`$-translation numbers are relatively dense there exists $`L=L(\frac{ϵ}{2})>0`$ such that any interval of length $`L`$ contains such a joint $`\frac{ϵ}{2}`$-translation number. Since $`_KA_K^{\frac{ϵ}{2}}=\mathrm{}`$, for any $`\eta >0`$ there exists $`K_0>0`$ such that

$$\left|A_K^{\frac{ϵ}{2}}[L,L]\right|<L\eta (KK_0)$$
(2.3)
where $`||`$ denotes the Lebesgue measure (see the definition of $`n_ϵ(x)`$). We claim that $`\widehat{𝕄}(A_K^ϵ)\eta `$ ($`KK_0`$). Otherwise $`\widehat{𝕄}(A_K^ϵ)>\eta `$ for some $`KK_0`$. Then by the definition of $`\widehat{𝕄}(A_K^ϵ)`$ there exists $`x_0`$ such that

$$_{x_0}^{x_0+L}\chi _{A_K^ϵ}(x)𝑑xL\eta .$$

Take a joint $`\frac{ϵ}{2}`$-translation number $`\tau [x_0L,x_0]`$, i.e. $`Lx_0+\tau 0`$. Then by Lemma 2.7, we have

$`\left|A_K^{\frac{ϵ}{2}}{\displaystyle [L,L]}\right|`$ $``$ $`{\displaystyle _{x_0+\tau }^{x_0+\tau +L}}\chi _{A_K^{\frac{ϵ}{2}}}(x)𝑑x`$
$`=`$ $`{\displaystyle _{x_0}^{x_0+L}}\chi _{A_K^{\frac{ϵ}{2}}}(y+\tau )𝑑y`$
$``$ $`{\displaystyle _{x_0}^{x_0+L}}\chi _{A_K^ϵ}(y)𝑑y`$
For the first inequality we have used the fact that $`[x_0+\tau ,x_0+\tau +L][L,L]`$ and for the last inequality we have used Lemma 2.7. What we have deduced contradicts (2.3). Thus Lemma 2.6 is proved.

We continue our proof of Theorem 2.5. Let $`S:=f^{}_{\mathrm{}}<\mathrm{}`$. Let $`K`$ be such that $`\widehat{𝕄}(A_K^ϵ)ϵ/S`$. We define first

$$g(x)=\{\begin{array}{cc}f^{}(x)\hfill & \text{ if }xA_K^ϵ\hfill \\ 0\hfill & \text{ if }xA_K^ϵ\hfill \end{array}\text{and}m(x)=\{\begin{array}{cc}n^ϵ(x)\hfill & \text{ if }xA_K^ϵ\hfill \\ 1\hfill & \text{ if }xA_K^ϵ\hfill \end{array}$$

Lemma 2.6 implies that

$$𝕄g𝕄f^{}+ϵ\text{a.e.}.$$
(2.4)
Moreover by the invariance of $`f^{}`$ we have (remember that $`u_0=1`$)

$$g(x)g(u_kx)\text{for a.e. }x\text{ and all }0km(x)1.$$
(2.5)
Then we have

$$f_{m(x)}(x)(g(x)+ϵ)m(x)\underset{k=0}{\overset{m(x)1}{}}g(u_kx)+ϵm(x).$$
(2.6)
We define inductively $`m_0(x)=0`$ and

$$m_k(x)=m_{k1}(x)+m(u_{m_{k1}(x)}x).$$

Now choose $`R>K`$ and let $`k(x)`$ be the maximal $`k`$ for which $`m_k(x)R`$. Note that $`Rm_{k(x)}(x)<K`$. Now we get by the subadditivity and Equation (2.6)

$$\begin{array}{cc}\hfill f_R(x)& \underset{k=0}{\overset{k(x)1}{}}f_{m(u_{m_k(x)}x)}(u_{m_k(x)}x)+\underset{0}{\underset{}{f_{Rm_{k(x)}(x)}(u_{m_{k(x)}(x)}x)}}\hfill \\ & \underset{k=0}{\overset{k(x)1}{}}\underset{j=m_{k1}(x)}{\overset{m_k(x)}{}}g(u_jx)+(m_k(x)m_{k1}(x))ϵ\hfill \\ & \underset{j=0}{\overset{m_{k(x)}(x)}{}}g(u_jx)+m_{k(x)}(x)ϵ\hfill \\ & \underset{j=0}{\overset{R1}{}}g(u_jx)+Rϵ+KS.\hfill \end{array}$$

Taking the Bohr mean, using that $`𝕄g=𝕄(gu_j)`$ and dividing by $`R`$ gives

$$\begin{array}{cc}\hfill \frac{1}{R}𝕄f_R& 𝕄g+ϵ+\frac{KS}{R}\hfill \\ & 𝕄f^{}+2ϵ+\frac{KS}{R}\hfill \end{array}$$

by Equation (2.4). Now we let $`R\mathrm{}`$ and we get

$$\gamma 𝕄f^{}.$$

We claim that $`f^{}\gamma `$ implies $`f^{}=\gamma `$ for a.e. $`x`$. Suppose this was not the case, then one could find $`ϵ>0`$, $`\delta >0`$ and an interval $`J=(0,L)`$ of length $`|J|=L=L_ϵ`$ such that $`|AJ|=\delta >0`$, where

$$A=\{x:f^{}(x)<\gamma ϵ\}$$

By the invariance of $`f^{}`$ we have $`u_kA=A`$ for all $`k`$. Hence, for all $`k`$

$$\frac{1}{u_kL}_0^{u_kL}\chi _A𝑑x=\frac{1}{L}_0^L\chi _A𝑑x>\frac{\delta }{L}>0.$$

Since $`lim\; sup_ku_k=+\mathrm{}`$ (see remark 2.3), we have $`\widehat{𝕄}A>\frac{\delta }{L}`$, and thus

$$\underset{¯}{𝕄}f^{}<\gamma \left(1\frac{\delta }{L}\right)+(\gamma ϵ)\frac{\delta }{L}<\gamma .$$

###### Remark 2.8.

One can prove a similar theorem for more general sequences $`(u_n(x))`$. In this case it seems to be necessary to assume $`L^1`$ Bohr-uniform distribution.

### 2.3. Oseledec theorem

Kingman’s theorem implies the following Oseledec type theorem (see Ruelle \[Ru\]).

###### Theorem 2.9.

Let $`\beta >1`$ be a real number. Let $`M:GL_d()`$ be a uniformly almost periodic function. Write

$$M_x^n=M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x).$$

Suppose the $`q`$-exterior products $`\frac{1}{n}\mathrm{log}(M_x^n)^q`$ have joint periods, for $`q=1,\mathrm{},d`$. Then there is $`\mathrm{\Gamma }`$ with $`\beta \mathrm{\Gamma }\mathrm{\Gamma }`$ of full Lebesgue measure (in the sense that $`\mathrm{\Gamma }`$ has 0 measure) such that if $`x\mathrm{\Gamma }`$ then

a) $`lim_n\mathrm{}(M_x^nM_x^n)^{\frac{1}{2n}}=\mathrm{\Lambda }_x`$ exists.

b) Let $`\mathrm{exp}\lambda _x^{(1)}<\mathrm{}<\mathrm{exp}\lambda _x^{(s)}`$ be the eigenvalues of $`\mathrm{\Lambda }_x`$ (where $`s=s(x)`$ and the $`\lambda _x^{(r)}`$ are reals), and $`U_x^{(1)},\mathrm{},U_x^{(s)}`$ the corresponding eigenspaces. Let $`m_x^{(r)}=dimU_x^{(r)}`$. We have $`\lambda _{\beta x}^{(r)}=\lambda _x^{(r)}`$ and $`m_{\beta x}^{(r)}=m_x^{(r)}`$ and

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}M_x^nv=\lambda _x^{(r)}\text{when}vV_x^{(r)}V_x^{(r1)}$$

for $`r=0,\mathrm{},s`$ where $`V_x^{(0)}=\{0\}`$ and $`V_x^{(r)}=U_x^{(1)}+\mathrm{}+U_x^{(r)}`$.

c) Moreover $`V_x^{(r)}`$ depends measurably on $`x`$ and $`M_xV_x^{(r)}=V_{\beta x}^{(r)}`$.

d) In addition the functions $`\lambda _x^{(r)}`$ and $`m_x^{(r)}`$ are constant a.e.

###### Proof.

This theorem follows in the standard way from the a.e. convergence of

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}\left(M_x^n\right)^q$$

for $`1qd`$, which in the classical case is insured by Kingman’s theorem. So we only need to check that for the functions $`f_n^{(q)}(x)=\mathrm{log}\left(M_x^n\right)^q`$ and the sequence $`u_n=\beta ^n`$ the assumptions of Theorem 2.5 hold.

First we note that $`M^1`$ is uniformly almost periodic because $`M`$ is uniformly almost periodic and $`M^1(x)GL_d()`$. Second we note that $`M^q`$ and $`(M^q)^1`$ are again uniformly almost periodic, since each entry is a rational function of the entries of $`M`$ and $`M^1`$, respectively. Hence,

$$\underset{x}{sup}M(x)^q+\underset{x}{sup}(M^q)^1(x)=W_q<\mathrm{}.$$

Subadditivity is obviously fulfilled since $`(M_x^{n+m})^q=(M_{\beta ^nx}^m)^q(M_x^n)^q`$.

Condition (2.1) follows from

$$(M_{\beta x}^n)^q=M(\beta ^nx)^q(M_x^n)^q(M(x)^q)^1(M_x^n)^q+W_q.$$

Finally by Remark 2.3 the sequence $`\beta ^n`$ is totally Bohr ergodic, so Theorem 2.5 applies. Assertions (a), (b) and (c) follows from Proposition 1.3 (see also the proof of Theorem 1.6) in \[Ru\].

Now we prove d). By Kingman’s theorem (Theorem 2.5), we have for almost all $`x`$

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}M_x^{nk}=\underset{n1}{inf}\frac{1}{n}𝕄\mathrm{log}M_x^{nk}(1ks).$$
(2.1)
On the other hand, by the properties of exterior product, if we write $`k=_{i=1}^{j1}m_x^{(si)}+\mathrm{}`$ with $`0j<s`$ and $`0\mathrm{}m_x^{(sj)}`$ we have

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}M_x^{nk}=m_x^{(s)}\lambda _x^{(s)}+\mathrm{}+m_x^{(sj+1)}\lambda _x^{(sj+1)}+\mathrm{}\lambda _x^{(sj)}.$$
(2.2)
We can solve $`\lambda _x^{(r)}`$ ($`1rs`$) from the system (2.1)-(2.2). The solution is independent of $`x`$ since it depends only on the right hand side terms in (2.1). Consequently $`m_x^{(r)}`$ is also independent of $`x`$. ∎

## 3. When $`\beta `$ is a Pisot–Vijayaraghavan number

We restrict our attention to the special case where $`\beta >1`$ is a Pisot–Vijayaraghavan (PV) number and $`f_n`$ are defined by (1.1). We will prove that, under some extra condition, the sequence $`n^1f_n`$ has joint periods and the Kingman theorem and the Oseledec theorem apply. To do this, we need a distortion lemma and some properties of PV-numbers.

### 3.1. Distortion lemmas

###### Lemma 3.1.

Let $`M:GL_d()`$ such that

$$D:=\underset{x}{sup}M(x)M(x)^1<\mathrm{}$$
(3.3)
Let $`(x_k)`$ and $`(y_k)`$ be two sequences in $``$ and let $`\theta _k=M(x_k)M(y_k)`$. Then for all $`n`$ and $`0𝐯^d`$ we have

$$\frac{|M(x_1)\mathrm{}M(x_{n1})M(x_n)𝐯|}{|M(y_1)\mathrm{}M(y_{n1})M(y_n)𝐯|}\mathrm{exp}\left(C\underset{k=1}{\overset{n}{}}D^k\theta _k\right).$$

where $`C=(sup_xM(x))^1`$.

###### Proof.

Let

$`Q_1(x)`$ $`=`$ $`𝕀,Q_k(x)=M(x_1)M(x_2)\mathrm{}M(x_{k1})(1<kn)`$
$`Q^n(x)`$ $`=`$ $`𝕀,Q^k(x)=M(x_{k+1})M(x_{k+2})\mathrm{}M(x_n)(1k<n).`$
We can write

$$\frac{|Q_n(x)𝐯|}{|Q_n(y)𝐯|}=\underset{k=1}{\overset{n}{}}\frac{|Q_k(y)M(x_k)Q^k(x)𝐯|}{|Q_k(y)M(y_k)Q^k(x)𝐯|}.$$
(3.4)
Setting $`E_k=M(x_k)M^1(y_k)𝕀`$ we get the following estimate for the numerator of the general term in the above product

$`|Q_k(y)M(x_k)Q^k(x)𝐯|`$ $`=`$ $`|Q_k(x)(𝕀+E_k)M(y_k)Q^k(x)𝐯|`$
$`=`$ $`|Q_k(x)(𝕀+E_k)Q_k(x)^1Q_k(x)M(y_k)Q^k(x)𝐯|`$
$`=`$ $`|(𝕀+\stackrel{~}{E}_k)Q_k(x)M(y_k)Q^k(x)𝐯|`$
where $`\stackrel{~}{E}_k=Q_k(x)E_kQ_k(x)^1.`$ It follows that

$$\frac{|Q_k(y)M(x_k)Q^k(x)𝐯|}{|Q_k(y)M(y_k)Q^k(x)𝐯|}𝕀+\stackrel{~}{E}_k.$$
(3.5)
It is obvious that

$$\stackrel{~}{E}_kD^{k1}E_k.$$
(3.6)
On the other hand

$$E_k\underset{x}{sup}M(x)^1\theta _k.$$
(3.7)
By combining (3.4), (3.5), (3.6) and (3.7), we obtain

$$\frac{|M(x_1)\mathrm{}M(x_{n1})M(x_n)𝐯|}{|M(y_1)\mathrm{}M(y_{n1})M(y_n)𝐯|}\underset{k=0}{\overset{n1}{}}\left(1+CD^k\theta _k\right).$$

If $`M(x)`$ is non-negative, the next lemma shows that condition (3.3) is not needed for positive vectors

###### Lemma 3.2.

Let $`M:GL_d()`$ be such that the entries of $`M(x)`$ are either identically zero or bounded from below by a positive number $`\delta >0`$ (independent of entries). Then for any sequences $`(x_k)`$ and $`(y_k)`$ in $``$ and for any non-negative vector $`𝐯`$ we have

$$\frac{|M(x_1)\mathrm{}M(x_{n1})M(x_n)𝐯|}{|M(y_1)\mathrm{}M(y_{n1})M(y_n)𝐯|}\mathrm{exp}\left(\frac{1}{\delta }\underset{k=1}{\overset{n}{}}\theta _k\right),$$

where $`\theta _k=M(x_k)M(y_k)`$ and the norm $`|v|=_{i=1}^d|v_i|`$ on $`^d`$ is chosen.

###### Proof.

We may write

$`|M(x_1)\mathrm{}M(x_{n1})M(x_n)𝐯|`$
$`=`$ $`{\displaystyle \underset{i_0,i_1\mathrm{},i_n}{}}M(x_1)_{i_0,i_1}M(x_2)_{i_1,i_2}\mathrm{}M(x_n)_{i_{n1},i_n}v_{i_n}`$
We have a similar expression for $`|M(y_1)\mathrm{}M(y_{n1})M(y_n)𝐯|`$. Now compare the two expressions term by term. By the hypothesis, both quantities $`M(x_1)_{i_0,i_1}`$ and $`M(y_1)_{i_0,i_1}`$ are either zero or larger than $`\delta `$. So, using the trivial inequality $`x/ye^{x/y1}`$ we have

$$M(x_1)_{i_0,i_1}M(y_1)_{i_0,i_1}e^{\delta ^1\theta _k}.$$

The same estimates hold for other pairs $`M(x_k)_{i_{k1},i_k}`$ and $`M(y_k)_{i_{k1},i_k}`$. The desired inequality follows. ∎

### 3.2. Two properties of PV–numbers

Let $`\beta >1`$ be a PV–number of order $`r`$. We denote its conjugates by $`\beta _1^{},\mathrm{},\beta _{r1}^{}`$. Then for $`n1`$, denote

$$F_n=\beta ^n+\beta _{}^{}{}_{1}{}^{n}+\mathrm{}+\beta _{}^{}{}_{r1}{}^{n}.$$

###### Lemma 3.3.

The number $`F_n`$ is an integer and we have

$$|\beta ^nF_n|(r1)\rho ^n(n1)$$

where $`\rho =\mathrm{max}_{1jr1}|\beta _j^{}|<1.`$

Given any real number $`\beta >1`$ (not necessarily integral), we can expand each number $`x[0,1)`$ in a canonical way into its $`\beta `$-expansion \[Re\] (see also \[P, Bl\]):

$$x=\underset{n=1}{\overset{\mathrm{}}{}}\frac{ϵ_n(x)}{\beta ^n}$$

where $`(ϵ_n(x))_{n1}`$ is a uniquely determined sequence in $`\{0,1,\mathrm{},[\beta ]\}^{}`$. We may also call $`(ϵ_n(x))_{n1}`$ the $`\beta `$-expansion of $`x`$. We note that not all sequences in $`\{0,1,\mathrm{},[\beta ]\}^{}`$ are $`\beta `$-expansions. Let $`D_\beta `$ be the set of all possible $`\beta `$-expansions of numbers in $`[0,1)`$. A finite sequence $`ϵ=(ϵ_1,\mathrm{},ϵ_n)`$ (of length $`n`$) in $`\{0,1,\mathrm{},[\beta ]\}^n`$ is said to be admissible if it is the prefix of the $`\beta `$-expansion of some number $`x`$. For such an admissible sequence, we define

$$I(ϵ_1,\mathrm{},ϵ_n)=\{x[0,1):ϵ_1(x)=ϵ_1,\mathrm{},ϵ_n(x)=ϵ_n\}.$$

It is known that if $`D_\beta `$ is endowed with the lexicographical order, the map which associates $`x`$ to its $`\beta `$-expansion is strictly increasing. The set $`I(ϵ)=I(ϵ_1,\mathrm{},ϵ_n)`$ is an interval, called a $`\beta `$-interval of level $`n`$. Its length is denoted by $`|I(ϵ)|`$.

###### Lemma 3.4.

Suppose $`\beta >1`$ is a PV–number. There is a constant $`C>0`$ such that

$$C^1\beta ^n|I(ϵ_1,\mathrm{},ϵ_n)|C\beta ^n$$

for any integer $`n1`$ and any $`\beta `$-interval $`I(ϵ_1,\mathrm{},ϵ_n)`$.

See \[F\] for proofs of Lemma 3.3 and Lemma 3.4.

### 3.3. Existence of joint periods

###### Definition 3.5.

Let $`\beta >1`$ be a positive real number and let $`M:_{d\times d}()`$. If the entries of $`M`$ are functions of the form $`f(\beta ^nx)`$ ($`n`$) where $`f`$ is $`1`$–periodic continuous, we say that $`M`$ is $`\beta `$–adapted u.a.p.

###### Remark 3.6.

The matrix $`M(x)`$ defined by (1.4) associated to a multiperiodic function is $`\beta `$–adapted.

###### Proposition 3.7.

Let $`\beta >1`$ be a PV–number and let $`M:\text{GL}_d()`$ be $`\beta `$–adapted and $`\alpha `$–Hölder continuous. Suppose that

$$D\rho ^\alpha <1$$
(3.8)
where $`\rho `$ is the maximal modulus of the conjugates of $`\beta `$ and $`D`$ is the same as in the distortion lemma (Lemma 3.1). Then for any $`1qd`$ the sequence $`n^1f_n^{(q)}(x)`$ has joint periods, where

$$f_n^{(q)}(x)=\mathrm{log}(M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x)))^q.$$

###### Proof.

Since $`M`$ is $`\beta `$–adapted, the entries of $`M(\beta ^kx)`$ are all of the form $`h_{i,j}(\beta ^{\mathrm{}_{i,j}}x)`$ with 1–periodic function $`h_{i,j}(x)`$ and integer $`\mathrm{}_{i,j}0`$ for sufficiently large $`k`$. So, if necessary, we consider $`\mathrm{log}(M(\beta ^{n1}x)\mathrm{}M(\beta ^{k_0}x))^q`$ for some sufficiently large but fixed $`k_00`$.

Consider $`\tau =\beta ^m\eta _m+\mathrm{}+\beta \eta _1+\eta _0`$ where $`m1`$ and $`0\eta _i\beta `$ are integers. We are going to show that all such $`\tau `$ are joint $`ϵ`$-translation numbers for $`n^1f_n^{(n)}(x)`$ with $`nN(ϵ)`$, where $`N(ϵ)`$ depending on $`ϵ`$ is an integer to be determined.

By Lemma 3.3, we have

$$\underset{j}{inf}|\beta ^k\tau j|C^{}\rho ^k$$

for all $`k`$ and some constant $`C^{}`$ independent of $`k`$ and $`\tau `$. For $`k0`$ each entry of $`M^q(\beta ^{k+k_0}x)`$ is a degree $`q`$ polynomial in $`d^2`$ variables of the form $`h(\beta ^{\mathrm{}+k}x)`$, with $`h`$ $`\alpha `$-Hölder, 1–periodic, and $`\mathrm{}0`$. Notice that we have $`h(\beta ^{\mathrm{}+k}(x+\tau ))=h(\beta ^{\mathrm{}+k}x)+O(\rho ^{k\alpha })`$, hence

$$M(\beta ^{k+k_0}(x+\tau ))^qM(\beta ^{k+k_0}(x))^q=C_q\rho ^{k\alpha }$$

for some constant $`C_q`$. By the distortion Lemma 3.1 and the above estimate we have

$`|f_n^{(q)}(x+\tau )f_n^{(q)}(x)|`$
$`=`$ $`\left|\mathrm{log}{\displaystyle \frac{(M(\beta ^{n1}x+\beta ^{n1}\tau )\mathrm{}M(\beta x+\beta \tau )M(x+\tau ))^q}{(M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x))^q}}\right|`$
$`=`$ $`C_qC^{}{\displaystyle \underset{k=1}{\overset{n}{}}}D^k\rho ^{(kk_0)\alpha }{\displaystyle \frac{CC^{}D\rho ^{k_0}}{1D\rho ^\alpha }}=:𝒞`$
So, we may choose $`N(ϵ)=𝒞/ϵ`$. In order to finish the proof, it suffices to notice that Lemma 3.4 implies that all these $`\tau `$ form a subset with bounded gap in $``$. ∎

###### Proposition 3.8.

Let $`\beta >1`$ be a PV–number and let $`M:\text{GL}_d()`$ be $`\beta `$–adapted and $`\alpha `$–Hölder continuous. Suppose that the entries of $`M(x)`$ are either identically zero or larger than a constant $`\delta >0`$. Then $`n^1f_n(x)`$ has joint periods, where

$$f_n^{(q)}(x)=\mathrm{log}(M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x))^q).$$

###### Proof.

The proof is the same as the last proposition. But we use Lemma 3.2 instead of Lemma 3.1. ∎

## 4. Multiperiodic functions

As we pointed out in the introduction and as we will see in Section 4.2, our scalar equation (1.2) can be converted to the vector equation (1.3). So, we first study the vector equation (1.3).

### 4.1. Equation $`G(x)=M(x/\beta )G(x/\beta )`$

Let $`M:_{d\times d}()`$ be a matrix valued function. We consider the following vector valued equation

$$G(x)=M\left(\frac{x}{\beta }\right)G\left(\frac{x}{\beta }\right)$$

where the unknown $`G:^d`$ is a vector valued function.

###### Theorem 4.1.

Let $`\beta >1`$ be a real number and $`M(x)`$ be a complex matrix valued function. Suppose that $`M`$ is Lipschitzian and that $`M(0)`$ is non negative and has $`1`$ as a simple eigenvalue with a corresponding strictly positive eigenvector $`𝐯`$. Then there exists, up to a multiplicative constant, a unique continuous solution $`G(0)0`$ of the equation $`G(x)=M(x/\beta )G(x/\beta )`$. The solution can be defined by

$$G(x)=\underset{n\mathrm{}}{lim}M\left(\frac{x}{\beta }\right)M\left(\frac{x}{\beta ^2}\right)\mathrm{}M\left(\frac{x}{\beta ^n}\right)𝐯$$

where the convergence is uniform on every compact subset in $``$.

###### Proof.

Write $`𝐯=(v_1,\mathrm{},v_d)^t`$. We introduce the following norm for $`^d`$

$$z=\underset{1jd}{\mathrm{max}}\frac{|z_j|}{v_j}(z=(z_j)_{1jd}^d).$$

Then a matrix $`A=(a_{i,j})_{d\times d}()`$, considered as an operator on the normed space $`(^d,)`$, admits its operator norm

$$A=\underset{1id}{\mathrm{max}}\frac{1}{v_i}\underset{j=1}{\overset{d}{}}|a_{i,j}|v_j.$$

Notice that $`M(0)=1`$ because $`M(0)𝐯=𝐯`$.

Since the eigenvalue $`1`$ of $`M(0)`$ is simple (and isolated), and $`M(x)`$ is Lipschitz continuous, by the perturbation theory of matrices, there is a neighborhood of $`0`$, say $`[\delta ,\delta ]`$ ($`\delta >0`$), such that for any $`x[\delta ,\delta ]`$, $`M(x)`$ has a simple eigenvalue $`\lambda (x)`$ and a corresponding eigenvector $`v(x)`$ satisfying

$$|\lambda (x)1|C|x|,v(x)𝐯C|x|(x[\delta ,\delta ])$$
(4.1)
for some constant $`C>0`$. We claim that the limit

$$G(x)=\underset{n\mathrm{}}{lim}M\left(\frac{x}{\beta }\right)M\left(\frac{x}{\beta ^2}\right)\mathrm{}M\left(\frac{x}{\beta ^n}\right)𝐯$$
(4.2)
exists (uniformly on any compact set). It is clear that the limit function is a solution.

Denote

$$Q_n(x)=M\left(\frac{x}{\beta }\right)M\left(\frac{x}{\beta ^2}\right)\mathrm{}M\left(\frac{x}{\beta ^n}\right).$$

The proof of the existence of the limit in (4.2) is based on the following lemma.

###### Lemma 4.2.

For any $`\delta >0`$, there exists a constant $`D>0`$ such that for any $`n1`$ and any $`x[\delta ,\delta ]`$ we have

$$Q_n(x)DQ_n(x)𝐯𝐯D|x|.$$

To get the boundedness of $`Q_n(x)`$, it suffices to notice that

$$Q_n(x)\underset{j=1}{\overset{n}{}}f\left(\frac{x}{\beta ^j}\right)$$

where the scalar function $`f(x)=M(x)`$ is Lipschitzian and $`f(0)=1`$ (we have used our choice of the norm of $`^d`$), and that the products converge uniformly on $`[\delta ,\delta ]`$ to a continuous function \[FL\]. Now we prove that

$$Q_n(x)𝐯Q_{n1}(x)𝐯C^{}\frac{|x|}{\beta ^n}$$
(4.3)
where $`C^{}>0`$ is some constant. In fact, since $`M(x)v(x)=\lambda (x)v(x)`$, we have

$$M\left(\frac{x}{\beta ^n}\right)𝐯𝐯=M\left(\frac{x}{\beta ^n}\right)\left[𝐯v\left(\frac{x}{\beta ^n}\right)\right]+\left[\lambda \left(\frac{x}{\beta ^n}\right)v\left(\frac{x}{\beta ^n}\right)𝐯\right].$$

Multiplying both sides by $`Q_{n1}(x)`$, we get

$`Q_n(x)𝐯Q_{n1}(x)𝐯`$
$``$ $`Q_n(x)\left(𝐯v\left({\displaystyle \frac{x}{\beta ^n}}\right)\right)+Q_{n1}(x)\left(𝐯\lambda \left({\displaystyle \frac{x}{\beta ^n}}\right)v\left({\displaystyle \frac{x}{\beta ^n}}\right)\right)`$
Notice that

$$\lambda (x)v(x)𝐯\lambda (x)1v(x)+v(x)𝐯.$$

Using the last inequality, the estimates in (4.1) and that we have just proved $`Q_n(x)D`$, we obtain (4.3). Then for $`n>m`$

$$Q_n(x)𝐯Q_m(x)𝐯\underset{k=m+1}{\overset{n}{}}Q_k(x)𝐯Q_{k1}(x)𝐯\frac{C^{}|x|}{\beta ^{m1}(\beta 1)}.$$

That means $`Q_n(x)𝐯`$ is a Cauchy sequence in the space $`C([\delta ,\delta ])`$ of continuous functions equipped with uniform norm. Since for any fixed integer $`n_0`$, we have

$$\underset{n\mathrm{}}{lim}Q_n(x)𝐯=Q_{n_0}(x)\underset{n\mathrm{}}{lim}Q_n\left(\frac{x}{\beta ^{n_0}}\right)𝐯,$$

it follows that the uniform convergence of $`Q_n(x)𝐯`$ on $`[\delta ,\delta ]`$ implies its uniform convergence on any compact set.

The uniqueness of solution is easy. Let $`G0`$ be a solution. First notice that $`G(0)`$ is an eigenvector of $`M(0)`$ associated to the simple eigenvalue 1. Hence we may assume that $`G(0)=𝐯`$. By iterating the equation, we get

$`G(x)=Q_n(x)G\left({\displaystyle \frac{x}{\beta ^n}}\right)=Q_n(x)𝐯+Q_n(x)\left(G\left({\displaystyle \frac{x}{\beta ^n}}\right)𝐯\right)`$
The last term converges to zero (uniformly on any compact set) because of $`Q_n(x)D`$. So, $`G(x)`$ must be the limit of $`Q_n(x)𝐯`$. ∎

###### Remark 4.3.

In the theorem, neither the almost periodicity of $`M(x)`$ nor the positivity of $`M(x)`$ is required, but only the positivity of $`M(0)`$. That $`1`$ is an eigenvalue of $`M(0)`$ is necessary for the equation (1.3) to have a solution $`G(x)`$ such that $`G(0)0`$.

###### Remark 4.4.

The Lipschitz continuity is not really necessary. Hölder continuity or even Dini continuity is sufficient.

###### Remark 4.5.

If the entries of $`M`$ are (real) analytic, then the solution $`G`$ is also analytic. Because, for any $`x_0`$, there is a disk on the complex plane centered at $`x_0`$ on which $`Q_n(x)𝐯`$ (as functions of complex variable $`x`$) uniformly converges.

### 4.2. Existence of multiperiodic functions

Here we give a proof of Theorem A based on Theorem 4.1.

Let $`M(x)`$ be as in (1.4). It is easy to see that the characteristic polynomial of $`M(0)`$ takes the form

$$P(u)=u^df_1(0)u^{d1}f_2(0)u^{d2}\mathrm{}f_{d1}(0)uf_d(0)$$

The consistency condition implies that $`1`$ is an eigenvalue of $`M(0)`$. Notice that

$$P^{}(1)=f_1(0)+2f_2(0)\mathrm{}+(d1)f_{d1}(0)+df_d(0)>0.$$

So, the eigenvalue $`1`$ is simple. By Theorem 4.1, there is a unique solution of $`G(x)=M(x/\beta )G(x/\beta )`$. Let

$$G(x)=(G_1(x),G_2(x),\mathrm{},G_d(x))^t$$

Then $`G_1(x)`$ is a solution of (1.2). If $`F`$ is a solution of (1.2). Let

$$\stackrel{~}{G}_1(x)=F(x),\stackrel{~}{G}_2(x)=F(x/\beta ),\mathrm{},\stackrel{~}{G}_d(x)=F(x/\beta ^{d1}).$$

Then $`\stackrel{~}{G}=(\stackrel{~}{G}_1,\mathrm{},\stackrel{~}{G}_d)^t`$ is a solution of $`G(x)=M(x/\beta )G(x/\beta )`$. Thus the uniqueness of the solution of Equation (1.3) implies that of Equation (1.2).

### 4.3. Asymptotic behavior of multiperiodic functions

Let us consider the asymptotic behavior of a multiperiodic function, or more generally the asymptotic behavior of a solution $`G`$ of Equation (1.3) provided it exists (the existence may be guaranteed by Theorem 4.1.

###### Theorem 4.6.

Let $`\beta >1`$ be a PV–number whose maximal conjugate has modulus $`\rho `$. Let $`M:\mathrm{GL}_d()`$ be a $`\beta `$–adapted u.a.p. Hölder function of order $`\alpha >0`$. Suppose that $`G`$ is a solution of $`G(x)=M(x/\beta )G(x/\beta )`$. Suppose furthermore that one of the following conditions is satisfied

(i) $`D\rho ^\alpha <1`$ where $`D=sup_xM(x)M(x)^1`$ (NB. $`\beta `$ must be Pisot).

(ii) The entries of $`M(x)`$ are either identically zero or larger than a constant $`\delta >0`$.

Then for a.e. $`x`$ the limit

$$h(x)=\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}|G(\beta ^nx)|$$

exists and is independent of $`x`$.

###### Proof.

We first consider the case (i). By Proposition 3.7, Theorem 2.9 applies. Hence for a.e. $`x`$, if we denote by $`r(x)`$ the integer such that the vector $`G(x)V_x^{(r)}V_x^{(r1)}`$ we get

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}|G(\beta ^nx)|=\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}|M_x^nG(x)|=\lambda _x^{(r(x))}.$$

But $`G(\beta x)=M(x)G(x)`$ hence $`G(\beta x)V_{\beta x}^{(r)}V_{\beta x}^{(r1)}`$, from what follows $`r(\beta x)=r(x)`$, i.e. $`r`$ is invariant. Hence constant a.e. because of the total Bohr ergodicity of the sequence $`\beta ^n`$.

Case (ii). We use the notation of Proposition 3.8. Since $`n^1f_n`$ has joint periods, Theorem 2.5 applies (see the proof of Theorem 2.9 for details). Hence the following limit exists a.e.

$$\underset{n\mathrm{}}{lim}\frac{1}{n}f_n(x)=,$$

where $`=inf_n\frac{1}{n}𝕄(f_n)`$. In view of $`G(\beta ^nx)=M(\beta ^{n1}x)\mathrm{}M(x)G(x)`$ the positivity of $`M`$ and $`G`$ gives

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}|G(\beta ^nx)|=\text{a.e.}$$

Note that when $`G`$ is the solution of equation (1.3) with $`M`$ and $`G`$ given by (1.4) we have

$$|G(x)|=|F(x)|+|F(x/\beta )|+\mathrm{}+|F(x/\beta ^{d1})|,$$
(4.4)
thus the asymptotic behavior of $`\frac{1}{n}\mathrm{log}|G(\beta ^nx)|`$ and $`\frac{1}{n}\mathrm{log}_{j=0}^{d1}|F(\beta ^{nj}x)|`$ are the same. Thus Theorem C follows as an immediate corollary of Theorem 4.6. This partially answers a question in \[JRS\] (Conjecture 4.1., p. 263).

We prove now Theorem B. By the primitivity of $`M(0)`$ and the hypothesis, there exists an ingeter $`\tau 1`$ such that $`\stackrel{~}{M}(x):=M(x/\beta ^{\tau 1})\mathrm{}M(x/\beta )M(x)`$ has all its entries strictly positive (even larger than $`c\delta ^\tau `$ for some constant $`c>0`$). Consider the equation

$$G(x)=\stackrel{~}{M}(x/\beta )G(x/\beta ^\tau ).$$

We examine the first entries of both sides. We can find two constants $`0<c_1<c_2`$ such that we get

$$c_1F(x)F(x/\beta ^{\tau +1})+\mathrm{}+F(x/\beta ^{\tau +d})c_2F(x).$$

Thus Theorem B follows from Theorem 4.6.

###### Theorem 4.7.

Under the same conditions as Theorem 4.6, for any $`q^+`$, the following limit exists

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}_0^1M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x)^q𝑑x.$$

###### Proof.

Write

$$Z_n=_0^1P_n(x)^q𝑑x\mathrm{with}P_n(x)=M(\beta ^{n1}x)\mathrm{}M(\beta x)M(x).$$

It suffices to show that there is a constant $`C>0`$ such that

$$Z_{n+m}CZ_nZ_m(n1,m1).$$

We assume that $`q=1`$, just for simplicity. We will use the fact that there is a constant $`L>0`$ such that

$$M(x)M(y)_2L|xy|^\alpha (x,y).$$

We use the notation $`\mathrm{\Pi }_{i=0}^nM_i=M_nM_{n1}\mathrm{}M_1M_0`$ for the (noncommutative) product of the matrices $`M_0,\mathrm{},M_n`$. Write

$`Z_{n+m}`$ $`=`$ $`{\displaystyle \underset{ϵ}{}}{\displaystyle _{I(ϵ)}}{\displaystyle \underset{k=0}{\overset{m1}{}}}M(\beta ^{n+k}x){\displaystyle \underset{j=0}{\overset{n1}{}}}M(\beta ^jx)𝑑x`$
$``$ $`{\displaystyle \underset{ϵ}{}}{\displaystyle _{I(ϵ)}}{\displaystyle \underset{j=0}{\overset{n1}{}}}M(\beta ^jx){\displaystyle \underset{k=0}{\overset{m1}{}}}M(\beta ^{n+k}x)𝑑x`$
where the sum is taken over all $`\beta `$-intervals $`I(ϵ)`$ of level $`n`$ (see Lemma 3.4). Let $`a_ϵ`$ be the left endpoint of $`I(ϵ)`$. The integral in the last sum , after the change of variables $`\beta ^n(xa_ϵ)=y`$, becomes

$$\beta ^n_0^{\beta ^n|I(ϵ)|}\underset{j=0}{\overset{n1}{}}M(\beta ^ja_ϵ+\beta ^{n+j}y)\underset{k=0}{\overset{m1}{}}M(\beta ^ky+\beta ^{n+k}a_ϵ)𝑑y.$$

Notice that if $`ϵ=(ϵ_1,\mathrm{},ϵ_n)`$, then

$$\beta ^{n+k}a_ϵ=\beta ^{n+k}\left(\frac{ϵ_1}{\beta }+\mathrm{}+\frac{ϵ_n}{\beta ^n}\right)=\beta ^{n+k1}ϵ_1+\mathrm{}+\beta ^kϵ_n.$$

So, by Lemma 3.3, there is an integer $`n_ϵ`$ such that

$$|\beta ^{n+k}a_ϵn_ϵ|=O(\rho ^k+\rho ^{k+1}+\mathrm{}+\rho ^{n+k1})=O(\rho ^k)$$
(4.5)
By the distortion lemma (Lemma 3.1, Lemma 3.2), we have

$$\underset{j=0}{\overset{n1}{}}M(\beta ^ja_ϵ+\beta ^{n+j}y)C\underset{j=0}{\overset{n1}{}}M(\beta ^ja_ϵ)$$

$$\underset{k=0}{\overset{m1}{}}M(\beta ^ky+\beta ^{n+k}a_ϵ)C\underset{k=0}{\overset{m1}{}}M(\beta ^ky)$$

Therefore, we get

$$Z_{n+m}C\beta ^n\underset{ϵ}{}P_n(a_ϵ)C^{}Z_nZ_m.$$

###### Corollary 4.8.

Let $`F`$ be the multiperiodic function defined by (1.2). Suppose that $`\beta >1`$ is a PV–number and that $`f_1,\mathrm{},f_d`$ are either identically zero or larger than a constant $`\delta >0`$. Suppose further that $`M_x^{\mathrm{}}:=M(\beta ^\mathrm{}1x)\mathrm{}M(\beta x)M(x)`$ has strictly positive entries for some integer $`\mathrm{}>0`$. Then for any $`q^+`$, the following limit exists

$$\underset{T\mathrm{}}{lim}\frac{1}{\mathrm{log}T}_0^TF(x)^q𝑑x.$$

###### Proof.

Without loss of generality, we may only consider the subsequence $`T_n=\beta ^n`$. Since $`|G(x)|=_{j=0}^{d1}|F(x/\beta ^j)|`$ where $`G`$ is the solution of the associated vector equation (1.3), we have only to show the existence of the limit

$$\underset{n\mathrm{}}{lim}\frac{1}{n}_0^{\beta ^n}|G(x)|^q𝑑x.$$

Making the change of variables $`x=\beta ^ny`$, we are led to prove the existence of the limit

$$\underset{n\mathrm{}}{lim}\frac{1}{n}_0^1|G(\beta ^nx)|^q𝑑x.$$

Notice that $`G(\beta ^nx)=M_x^nG(x)`$. Notice also that $`G(x)`$ has strictly positive entries by the hypothesis on $`M_x^{\mathrm{}}`$. So, for the non negative matrix $`M_x^n`$ we have

$$C^1M_x^n|G(\beta ^nx)|CM_x^n(x[0,1])$$

for some constant $`C>0`$. By the proof of the last theorem, $`\mathrm{log}_0^1|G(\beta ^nx)|^q𝑑x`$ is subadditive. ∎

###### Remark 4.9.

Let $`M(x)`$ be the matrix defined by (1.4). Let $`\stackrel{~}{M}`$ be the numerical matrix obtained by replacing $`f_j(x)`$ in $`M(x)`$ by $`0`$ or $`1`$ according to $`f_j(x)0`$ or not. Then $`\stackrel{~}{M}^{\mathrm{}}>0`$ implies $`M_x^{\mathrm{}}>0`$. In particular, $`\stackrel{~}{M}^d>0`$ if $`f_j(x)`$ are all strictly positive.

###### Example 4.10.

Let $`\beta >1`$ be a PV–number. Let $`f_1(x)`$ and $`f_2(x)`$ be two strictly positive $`1`$–periodic Hölder continuous functions such that $`f_1(0)+f_2(0)=1`$. There is a unique multiperiodic function $`F`$ defined by

$$F(x)=f_1\left(\frac{x}{\beta }\right)F\left(\frac{x}{\beta }\right)+f_2\left(\frac{x}{\beta ^2}\right)F\left(\frac{x}{\beta ^2}\right)$$

For almost every $`x`$, $`n^1\mathrm{log}F(\beta ^nx)`$ has a limit as $`n\mathrm{}`$; for any $`q^+`$, $`(\mathrm{log}T)^1_0^TF(x)^q𝑑x`$ has a limit as $`T\mathrm{}`$.

###### Example 4.11.

Let $`\beta >1`$ and $`a,b`$. Consider the contractive transformations on $``$ defined by

$$S_1x=\frac{x+a}{\beta },S_1x=\frac{x+b}{\beta ^2}.$$

For any $`0<p<1`$, there exists a unique probability measure $`\mu `$ with compact support such that

$$\mu =p\mu S_1^1+(1p)\mu S_2^1.$$

It is a self-similar measure. Its Fourier transform satisfies the equation

$$\widehat{\mu }(x)=f_1(x/\beta )\widehat{\mu }(x/\beta )+f_2(x/\beta ^2)\widehat{\mu }(x/\beta ^2)$$

with $`f_1(x)=pe^{2\pi iax}`$ and $`f_2(x)=qe^{2\pi ibx}`$ with $`q=1p`$ . This is a special case of the equation (1.2). The corresponding matrix defined by (1.4) and its inverse are respectively equal to

$$M(x)=\left(\begin{array}{cc}pe^{2\pi iax}\hfill & qe^{2\pi ibx/\beta }\hfill \\ 1\hfill & 0\hfill \end{array}\right),M(x)^1=q^1e^{2\pi ibx/\beta }\left(\begin{array}{cc}0\hfill & qe^{2\pi ibx/\beta }\hfill \\ 1\hfill & pe^{2\pi iax}\hfill \end{array}\right)$$

If we take the norm $`|v|=\mathrm{max}(|v_1|,|v_2|)`$ on $`^2`$, the operator norms for $`M(x)`$ and $`M(x)^1`$ are respectively $`M(x)=1`$ and $`M(x)^1=\frac{1+p}{1p}`$. So, when $`\beta `$ is a PV–number, under the condition $`\frac{1+p}{1p}<\frac{1}{\rho }`$, for almost all $`x`$ the following limit exists and does not depend on $`x`$

$$\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}\left(|\widehat{\mu }(\beta ^nx)|+|\widehat{\mu }(\beta ^{n1}x)|\right).$$


text:
# Are pulsar glitches triggered by a superfluid two-stream instability?

## Abstract

Mature neutron stars are expected to have several superfluid components. Strong evidence for this is provided by the glitches that have been observed in dozens of pulsars. The underlying idea behind most glitch models is that, as the neutron star crust spins down due to the emission of electromagnetic radiation, the superfluid component lags behind until a critical point is reached and angular momentum is transferred from the superfluid to the crust, leading to the spin-up associated with the glitch. In this Letter we describe a superfluid analogue of the two-stream instability that is well known in plasma physics, and provide evidence that this instability is likely to be relevant for neutron stars. This is a new physical mechanism which may play a key role in explaining the glitch mechanism and which could also prove to be relevant in laboratory experiments on superfluid Helium.

Neutron stars are cosmic laboratories of extreme physics glen . With a mass of roughly one and a half times that of the Sun compressed inside a radius of ten kilometres, they have density several times that of nuclear matter, strong magnetic fields and may spin rapidly. By now nearly 1300 pulsars (rotating neutron stars) have been observed lori . In fact, estimates based on the local population and the supernova rate suggest that the total number of neutron stars in the galaxy may be as large as $`2\times 10^8`$. Despite decades of theoretical effort the physics of these compact objects is not well understood. For example, there are open issues concerning the pulsar radiation mechanism and we do not know the detailed supranuclear equation of state. The main reason for the many uncertainties is that the description of these objects requires physics well beyond our laboratory experience.

Certain aspects of neutron star physics are, however, well established. It is clear that the overwhelming majority of these objects must be very cold in the sense that their temperatures are significantly below the Fermi temperatures of neutrons and protons (typically about $`10^{12}`$ K). Nuclear physics calculations of the transition temperature to superfluidity consistently yield a value of the order of $`10^9`$ K for the neutrons and protons. A newly born neutron star is expected to cool below this temperature within a few weeks to months following the supernova explosion. Thus most neutron stars in our galaxy are expected to contain superfluid components sauls . The solid ion crust that forms the outer kilometer or so of the star will be permeated by superfluid neutrons, while the fluid in the core is expected to contain superfluid neutrons and superconducting protons. The deep core may contain exotic phases of matter like hyperons and/or deconfined quarks, which may also become superfluid taka ; alfo . This means that astrophysical neutron stars provide a manifestation of large scale quantum phenomena, and if we want to understand their dynamics we need to allow for the presence of, potentially weakly coupled, superfluid components.

Even though a superfluid is locally irrotational, it can mimic ordinary fluid rotation by forming an array of vortices (with typical separation smaller than $`10^{11}`$ m). In the core of each vortex the superfluidity is destroyed, the particles are in an ordinary fluid state and can thus carry non-zero vorticity. Superfluids have zero viscosity, but dissipative mechanisms can exist when vortices are present. Particularly important is the direct friction between the different fluid components (known as ”mutual friction”). Another consequence of the interaction between the fluids is the non-dissipative entrainment effect. It arises because the bare neutrons (or protons) are ”dressed” by a polarization cloud of nucleons comprised of both neutrons and protons. Since both types of nucleon contribute to the cloud the momentum of the neutrons, say, is modified so that it is a linear combination of the neutron and proton currents. When one of the nucleon fluids begins to flow it will, through entrainment, induce a momentum in the other. This has the effect that a portion of the protons (and electrons) will be pulled along with the superfluid neutrons that surround each vortex.

Strong observational support for the presence of superfluid components in a neutron star is provided by the so-called glitches lyne : sudden spin-up events where the observed rotation rate (presumably that of the stars crust) may increase by as much as one part in $`10^6`$. The glitches typically have a relaxation time of the order of days to months, which is much longer than can be explained in terms of normal viscosity. Very soon after the first glitch in the Vela pulsar rada ; reich  was observed in 1969 it was argued baym  that the long relaxation time indicates the presence of a neutron superfluid. The glitch events, and the associated relaxation, are usually explained in terms of superfluid vortex dynamics. This is natural since the superfluid can only change its rotation rate if the vortices move. The standard model describes a glitch as an event in which a significant number of vortices are suddenly unpinned from the crust nuclei itoh , angular momentum is transferred to the crust and the vortices are eventually repinned. In particular, the so-called ”vortex creep” model cheng  is known to provide an adequate description of the glitch relaxation. The actual mechanism that triggers the glitch remains unspecified in virtually all models, however. This may be one of the most important unresolved issues in this area of research.

In this Letter we propose that a superfluid two-stream instability (analogous to that known to be relevant for electrons streaming through ions in plasma physics and ) plays a key role in triggering pulsar glitches. The key parameter governing the onset of this dynamical instability is the rotational lag between a superfluid component and the rest of the star. The existence of such a lag is generally assumed in most glitch models and there is observational evidence that glitches occur once this lag builds up to a critical level. From a statistical analysis of 48 observed glitches in 18 pulsars, Lyne, Shemar and Graham Smith lyne  deduce that the glitch activity depends linearly on the spin-down rate (as one would expect if a glitch is due to the release of a built up rotational difference). They estimate that a typical value for the critical velocity is $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_\mathrm{p}5\times 10^4`$ (where $`\mathrm{\Omega }_\mathrm{p}`$ and $`\mathrm{\Omega }_\mathrm{n}`$ are the rotation rates of the crust and the superfluid neutrons, respectively, and $`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{n}\mathrm{\Omega }_\mathrm{p}`$ represents the extent to which the superfluid lags behind as the crust spins down).

In a recent paper acp1  we have studied a superfluid analogue of the two-stream instability. As far as we are aware, this mechanism has not previously been discussed in this context. Yet, it could turn out to be of extreme importance for all superfluid systems that exhibit relative motion. The two-stream instability is analogous to the familiar Kelvin-Helmholtz instability draz . Its key distinguishing feature is that the two fluids are interpenetrating rather than in contact across an interface. In principle, the instability may operate in any system with two interpenetrating components moving at different rates. Unstable waves are such that they can be associated with a negative energy (in the sense that the energy in the perturbed flow is smaller than the energy of the unperturbed system). For this to be possible the wave must be such that it moves forwards with respect to one of the fluids, but backwards according to an observer riding along with the other fluid. Here we consider a neutron star in which the two fluids rotate at constant rates. A pulsation mode of the star, which is proportional to $`e^{i(\omega t+m\phi )}`$, will carry negative energy if the angular velocity of the mode pattern (the “pattern speed”) $`\sigma _p=\omega /m`$ lies in between the two rotation rates, i.e. when $`\mathrm{\Omega }_\mathrm{p}<\sigma _p<\mathrm{\Omega }_\mathrm{n}`$, assuming that the superfluid neutrons lag behind as the charged component is being spun down through electromagnetic braking.

The results described in this Letter were borne out of an effort to understand the dynamics of neutron stars with one, or several, superfluid components. The simplest models involve two interpenetrating fluids, eg. the superfluid neutrons in the inner crust and the core and all charged constituents (the ions in the crust, core protons, and all electrons) and the entrainment effect that acts between them ac ; pca ; pr . We have recently demonstrated acp1  that the two-stream instability may operate in these systems provided that the two fluids are coupled, either ”chemically” or via the entrainment. In order to provide support for the idea that the superfluid two-stream instability is relevant for pulsar glitches we consider a model problem of two fluids, allowed to rotate at different rates, inside a thin shell. By assuming that the shell is infinitesimally thin we may ignore radial motion, which means that the system permits only toroidal velocity perturbations. In fact, all oscillation modes of this shell model are closely related to the inertial r-modes of rotating single fluid objects pap . A somewhat laborious calculation acp2  leads to the following dispersion relation for oscillation modes of this system

$`\left\{l(l+1)[1x_\mathrm{p}(1+\epsilon )](\kappa +\mathrm{\Lambda })2(1x_\mathrm{p})\mathrm{\Lambda }+x_\mathrm{p}\epsilon (l1)(l+2)(1\mathrm{\Lambda })\right\}`$
$`\times \{l(l+1)[1\epsilon ](\kappa +1)2\epsilon (l1)(l+2)(1\mathrm{\Lambda })\}[l(l+1)]^2x_\mathrm{p}\epsilon ^2(\kappa +\mathrm{\Lambda })(\kappa +1)=0.`$ (1)
Here we have assumed that the waves are expanded in standard spherical harmonics $`Y_l^m(\theta ,\phi )`$ with integer indices $`l`$ and $`|m|l`$ . Furthermore, $`x_\mathrm{p}0.15`$ is the proton fraction, $`\kappa =\omega /m\mathrm{\Omega }_\mathrm{p}`$ represents the frequency of oscillation as measured with respect to the rotation rate of the charged component, and we have introduced a dimensionless measure of the relative rotation $`\mathrm{\Lambda }=\mathrm{\Omega }_\mathrm{n}/\mathrm{\Omega }_\mathrm{p}`$. The entrainment effect is represented by the coefficient $`\epsilon `$, which is expected to lie in the range $`0.4\epsilon 0.7`$ (see Ref. pca  for further discussion). A detailed analysis of the roots of the above dispersion relation shows that complex frequency modes will exist (for any multipole $`l`$) provided that the relative rotation rate $`\mathrm{\Lambda }`$ is sufficiently large.

As a ”physically realistic” example, we consider the case when $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_\mathrm{p}=5\times 10^4`$ which corresponds to $`\mathrm{\Lambda }=1.0005`$. We then find that we must have $`l>65`$ in order for the region of instability to overlap with the anticipated range of the proton fraction and the entrainment. This indicates that the instability will operate on length scales shorter than $`\pi R/l500`$ m (if we take the shell radius to be $`R=10`$ km, which would be the size of a typical neutron star). The region of instability for this rotational lag and $`l=100`$ is illustrated in Figure 1. Since our results only depend on the azimuthal index $`m`$ through the scaling $`\mathrm{Im}\omega =m\mathrm{\Omega }_\mathrm{p}\mathrm{Im}\kappa `$, and the dispersion relation does not depend explicitly on $`m`$, the fastest growth time corresponds to the $`m=l`$ mode. By analyzing the large $`l`$ limit of the dispersion relation one finds that the fastest possible growth time for an unstable mode can be approximated by

$$t\frac{6.7\times 10^2}{l}\left(\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }_\mathrm{p}}\right)^1\left(\frac{P}{0.1\text{ s}}\right)\text{ s}$$
(2)
This simple formula provides a useful estimate of the fastest growth rate of the two-stream instability for different parameter values, and we find that for a star rotating at the rate of the Vela pulsar, $`P=89`$ ms, we would have $`t1.2`$ s for $`l=m=100`$. Modes corresponding to higher multipoles may grow even faster, but as we will discuss later they are also faster damped by viscosity. Interestingly, this predicted growth time is significantly shorter than the resolved rise time of a large Vela glitch dods : $`t<40`$ s.

Our results suggest that the two-stream instability is generic in dynamical superfluids, and could possibly limit the relative flow in any superfluid system. Furthermore, the results for the shell model problem demonstrate that entrainment provides a sufficiently strong coupling for the instability to operate at astrophysically relevant relative flows. We believe that these results indicate that the two-stream instability is relevant for neutron stars, and that it may be the agent that triggers pulsar glitches. However, much further research is required if we are to understand this new mechanism in detail. For example, it is crucial to establish that the unstable modes grow faster than the relevant dissipation timescales. In the case of a superfluid neutron star core the main dissipation mechanisms are likely to be mutual friction and shear viscosity due to electron-electron scattering. Using estimates due to Mendell (see Eq. 38 in Ref. mend ) one can argue that the shear viscosity will be the dominant dissipation mechanism for all but the lowest values of $`l`$. To estimate the shear viscosity damping we can use results obtained for the secular r-mode instability akrev . Combining Eq. 43 in Ref. ksteg  with the viscosity coefficient for electron-electron scattering akrev  we estimate the viscous damping time as

$$t\frac{1.2\times 10^4}{l^2}\left(\frac{1.4M_{}}{M}\right)\left(\frac{R}{10\text{ km}}\right)^5\left(\frac{T}{10^7\text{ K}}\right)^2\text{ s}$$
(3)
It is interesting to compare this damping timescale to the growth rate of the unstable modes in our shell model. We find that, in order to overcome the viscous damping, an unstable mode must correspond to

$$l<90\left(\frac{\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_\mathrm{p}}{5\times 10^4}\right)\left(\frac{0.1\text{ s}}{P}\right)\left(\frac{1.4M_{}}{M}\right)\left(\frac{R}{10\text{ km}}\right)^5\left(\frac{T}{10^7\text{ K}}\right)^2.$$
(4)
In the case of the Vela pulsar (the archetypal glitching neutron star) the core temperature can be estimated to be $`5\times 10^7`$ K (roughly two orders of magnitude higher than the observed surface temperature page ). At this temperature modes with $`l>2500`$ or so are likely to be stabilized by shear viscosity. Given that our results indicate that the two-stream instability is active for much smaller values of $`l`$, see for example Figure 1, we conclude that dissipation is unlikely to suppress the instability in sufficiently young neutron stars. Naturally, the situation changes as the star cools further. Our estimates show that shear viscosity will suppress all unstable modes if the core temperature is below a few times $`10^6`$ K. Hence, the two-stream instability may not be able to overcome viscosity in a sufficiently cold neutron star. This would be consistent with the absence of glitches in mature pulsars lyne .

If the superfluid two-stream instability is, indeed, relevant for pulsar glitches then what is its exact role? Much further work is required to answer this question, but it is interesting to speculate about the possibilities. Most standard models for glitches are based on the idea of catastrophic vortex unpinning in the inner crust itoh ; cheng . An interesting scenario is the thermally induced glitch model due to Link and Epstein link . They show that the deposit of $`10^{42}`$ erg of heat in the crust would be sufficient to induce a large Vela glitch. The mechanism that triggers the glitch, eg. by depositing heat in the crust, is typically not identified in these models. We propose that the two-stream instability may fill this gap in the theory. Of course, glitches need not originate in the inner crust. Jones pbj  has argued that the vortex pinning is too weak to explain the size and frequency of the Vela glitches. Indeed, the observational evidence for free precession stairs  in PSR B1828-11 indicates the absence of significant vortex pinning. If this argument is correct then the large pulsar glitches must be due to some mechanism operating in the core fluid. Since our model problem was based on equations relevant for the conditions that prevail in a neutron star core it is clear that the two-stream instability may equally well serve as a trigger mechanism for glitches originating there. One possibility would be that the unstable mode grows so large that the superfluid degeneracy is broken, which would lead to an immediate coupling to the other normal fluid components. The key requirement for the two-stream instability to operate is the presence of a rotational lag. Such a lag can build up both when there is a strong coupling between the two fluids (i.e. when the vortices are pinned to the crust) and when this coupling is weak lang . A key issue for future theoretical work on pulsar glitches concerns to what extent a rotational lag can build up in various regions of the star.

Before we conclude this Letter it is worth emphasizing that we expect the two-stream instability to be generic in superfluid systems. The two-fluid equations that we have employed are analogous to the standard Landau model for superfluids, which means that the instability should be relevant also for laboratory systems. It is, in fact, not inconceivable that the notion that this instability triggers pulsar glitches may be tested using experiments on superfluid 4He.

NA and RP acknowledge support from the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137). NA acknowledges support from the Leverhulme Trust in the form of a prize fellowship. GC acknowledges partial support from NSF grant PHY-0140138.