text:
# Random model for RNA interference yields scale free network

## I Introduction

Although biology on the whole is a “knowledge based” discipline, with a strong traditional bias towards a reverse engineering approach to evolution, and a “form follows function” approach, prominant workers in the field have stressed that natural selection could very well have operated on complex structures already present in the pre-biotic world. Eigen Eigen  pointed out that non-linear out of equilibrium systems were capable of amplifying random fluctuations participating in feed-back loops, while Kauffman Kauffman1 ; Kauffman2  introduced random Boolean networks as null models of the genetic regulatory mechanism, and showed that they could spontaneously give rise to a great degree of complexity. Meanwhile, within the last two decades, we have learned a great deal about Self-Organized Criticality Bak , namely the spontaneous emergence of scale free structures in open systems far from equilibrium, driven by conserved fluxes. Thus, within a statistical mechanics context, it is much more natural to consider an ensemble of different states among which complex structures arise spontaneously. Evolutionary processes can be regarded as inducing dynamics on the distribution of states in this high dimensional phase space.

In this paper we introduce a null model for gene interactions resulting in the regulation of gene expression. This model is based on sequence matching on a random bit-string representation of the chromosome and is therefore radically different from the random Boolean networks which have been considered before Kauffman1 ; Kauffman2 ; Weisbuch ; Derrida1 ; Derrida2 ; Flyvberg . It can be considered as an out of equilibrium system subject to a constant mutation rate, achieving a steady state invariant under further mutations.

We present simulation results which display many qualitative features of gene regulation networks found in nature Maslov1 ; Maslov2 . We find that the in- and out- degree distributions are qualitatively different from each other, the out-degree distribution exhibiting power law decay with $`n(k_{\mathrm{out}})k_{\mathrm{out}}^\gamma `$, with $`\gamma =1`$ and log- periodic oscillations for relatively small $`k`$. The in-degree distribution, on the other hand, is much more localized.

The network has smallest world characteristics, with the cluster diameter, or the average minimum path length being essentially independent of the cluster size as long as the network consists of one cluster.

The average clustering coefficients for the in and out bonds have been calculated as 0.648 and 0.034 respectively. It is interesting to note that the clustering coefficients for the in-bonds behave much like for classical random networks Erdos , while the out-bonds have a clustering coefficient typical of scale free networks. Dorogovtsev ; Barabasi

In section 2, we will first motivate and then define the model. In section 3, we will present our simulation results. In section 4, we consider a toy model which mimics some of the features we observe in our simulations, and indicate ways in which this toy model may be improved in order to provide insights into how our basic model works. Section 5 contains our conclusions and a discussion.

## II Modelling RNA interference

### II.1 Genomic regulatory networks

Protein networks, which are an important component of transcriptional gene regulation networks biobook1 ; Sole , display a scale free structure, with the out-degree distribution characterised by a power law $`n(k)k^\gamma `$, with the exponent $`\gamma 2.5`$Maslov1 ; Maslov2  Gene regulation networks actually operate at many different levels. Maniatis  Post-transcriptional gene regulation, or RNA interference Hannon , is a mechanism where RNA strips may go and directly bind upon complementary segments on messenger RNA destined to be translated into some protein, thereby suppressing the production of this protein. Although we are not aware of a scaling analysis of post-transriptional gene interaction, it would be a fair guess to assume that interaction complexes of various sizes may arise in this type of interaction as well.

In all protein-protein or intra-genomic interactions, as well as the transcription and translation mechanism itself, essential lock-and key mechanisms are in operation. For normal translation to take place, rRNA in the ribosomes must be able to recognize and match the different amino acids and the corresponding three-letter anti-codon on the mRNA. In this, the rRNA is aided by the intermediary tRNA, which assumes a very specific three dimensional structure depending on the amino acid to which it binds. In transcriptional gene regulation, certain proteins known as transcription factors (TF) must first be synthesized, and then go and bind onto specific “promoter” sites preceding the coding part of a gene, in order that the RNA polymerase might start transcribing the DNA code into the mRNA. Conversely, the binding of other proteins onto the same promoter sites may block the binding of the TF, and thereby block the production of the mRNA. Hannon ; biobook2  In both cases, such binding presupposes steric and chemical specificity. In RNA interference, the short interfering RNA (siRNA) strips bind onto complementary sequences on the mRNA, via Watson-Crick base pairing. Hannon ; biobook3

Although it seems, from the above, as if there is a great diversity in these lock-and-key mechanisms, it should be realised that they all eventually match linear codes, even though this matching may take a few intermediary steps. The three dimensional structures (so called secondary structures) which come into play either in tRNA or the TF, are actually determined by either the sequence of ribo-nucleic acids on the tRNA, or the sequence of amino acids (primary structure) of the protein constituting the TF. Clearly the simplest is direct Watson-Crick base pairing between complementary sequences, and it is this latter, as it appears in RNA interference, which we will take to be the paradigm for our model. In fact we will further simplify the matching condition to consist of the identity relation rather than complementarity, since both are one-to-one, we believe that this should not change our results.

### II.2 The Model

The model is defined as follows. We postulate a “chromosome” to consist of a sequence of fixed length L, of independently and identically distributed random numbers with the probability distribution

$$P(x)=p\delta (x2)+(1p)/2[\delta (x1)+\delta (x)]$$
(1)
We define a “gene” to consist of a sequence of 0’s and 1’s situated between the $`i`$th and $`i+1`$st occurance of the symbol “2,” and we will denote the $`i`$th gene by

$$G_i=\{x_{i,1},x_{i,2},\mathrm{}x_{i,\mathrm{}_i}\}i=1,\mathrm{},s$$
(2)
where $`x_{i,\mu }2,\mu =1,\mathrm{}\mathrm{}_i`$ and $`\mathrm{}_i`$ is the length of the $`i`$th gene. We have used periodic boundary conditions, but one could just as well agree to end the chromosome always with a “2” at the $`L+1`$st site. Then

$$\underset{i}{\overset{s}{}}\mathrm{}_i=Ls,$$
(3)
where $`s`$ is the number of genes (the number of times the symbol 2 appears) on the chromosome. Let $`n_{\mathrm{}}`$ be the number of genes of length $`\mathrm{}`$. It obeys the sum rule

$$\underset{\mathrm{}=0}{\overset{Ls}{}}n_{\mathrm{}}=s.$$
(4)
For a given number $`s`$ of genes, the number of possible realisations of a given set $`\{n_{\mathrm{}}\}`$ is

$$K[\{n_{\mathrm{}}\}]=\frac{s!}{_{\mathrm{}=0}^{Ls}n_{\mathrm{}}!},$$
(5)
and the most probable distribution $`\overline{n}(\mathrm{})`$ can easily be found by using Lagrange multipliers, to be $`\overline{n}(\mathrm{})=Lp^2(1p)^{\mathrm{}}`$, in the limit of large $`L`$.

With these definitions we obtain a sequence of genes, seperated from each other by the symbol 2. (See Fig. (1)). In case there are more than one consecutive 2’s, they will be considered to bracket null genes. Clearly, for large L, the number of non-null genes, $`N`$, will fluctuate around $`\overline{N}=LpLp^2`$.

Each of the non-null genes constitutes a node in our gene regulation network. The interactions do not depend on the proximity of the genes along the chromosome. We define the adjacency matrix $`w_{ij}`$ by the matching condition such that

$$w_{ij}=\{\begin{array}{cc}1\hfill & G_iG_j\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.$$
(6)
By $`G_iG_j`$ we mean $`x_{i,\mu }=x_{j,\mu +\nu }`$ for $`\mu =1,\mathrm{}\mathrm{}_i`$ for at least one integer $`\nu `$ such that $`0\nu \mathrm{}_j\mathrm{}_i`$. Note that $`w_{ij}=1`$ implies that $`\mathrm{}_i\mathrm{}_j`$; in the case of the equality, $`w_{ij}=1`$ if and only if the two sequences $`G_i`$ and $`G_j`$ are identical, i.e., congruent. Thus, two genes are said to interact if the sequence $`G_i`$ occurs at least once as an unbroken subsequence of $`G_j`$, i.e., if one can be embedded in the other at least once.

Clearly this adjacency (or connectivity) matrix is directed. Moreover connectivity is “transitive” in the sense that $`w_{ij}=w_{jk}=1`$ implies that $`w_{ik}=1`$. The latter condition gives rise to a preferential attachment of incoming bonds to large genes, while small genes have an enhanced distribution of out-bonds. However, we will see in the next section that the degree distribution is scale free for the out-bounds, but not for the in-bonds.

A simple argument tells us that our network is of the smallest world Dorogovtsev  type. If we take the out-bonds, we see that due to the transitivity, any two successive edges linking say vertex $`ij`$ and $`jk`$, necessarily imply the existence of another directed edge $`ik`$. Thus $`l_{\mathrm{min}}^{(o)}1`$, as long as the network consists of a single cluster. The in-bonds follow the same argument, giving $`l_{\mathrm{min}}^{(i)}1`$ for a network with one giant cluster.

The question of whether a giant cluster always exists or whether we can identify the analog of a “percolation threshold” we will address in the next section, where we will also report numerical results for the minimum path length for undirected bonds (i.e., ignoring the directionality of the edges.) We expect that the qualitative behaviour of the network should not depend on $`p`$ as long as $`p`$ is bounded away from the percolation threshold. Requiring the number of vertices, $`N`$, to be larger than unity, i.e., $`\overline{N}=LpLp^2=Lp(1p)>1`$ gives $`p(1p)>1/L`$ (however this lower limit turns out not to be tight enough, i.e., $`p_c>1/L`$). The average gene size, again in the limit of large $`L`$ is $`\mathrm{}=(1p)/p`$. The obvious requirement for a non-trivial network, that $`\mathrm{}>1`$ yields $`p<1/2`$. We therefore expect to find scaling behaviour, if any, for $`p_c<p<1/2`$.

The elements $`w_{ij}`$ of the connectivity matrix are equal to unity with probabilities

$$P_{ij}=P(\mathrm{}_i,\mathrm{}_j),$$
(7)
which depend on the lengths $`(\mathrm{}_i,\mathrm{}_j)`$ through $`\mathrm{}_i`$ and $`\nu =\mathrm{}_j\mathrm{}_i`$. It is trivial to see that

$$P(\mathrm{},\mathrm{})=\left(\frac{1}{2}\right)^{\mathrm{}}.$$
(8)
We may make a mean field theory type of approximation to $`P(\mathrm{},\mathrm{}+\nu )`$ by neglecting the correlations between overlapping subsequences of $`G_j`$, and obtain,

$$P(\mathrm{},\mathrm{}+\nu )\left(\frac{1}{2}\right)^{\mathrm{}}(1+\nu ).$$
(9)
In the Appendix we have computed $`P(\mathrm{},\mathrm{}+\nu )`$ explicitly for $`\nu =1,2`$. However, so far it has not been possible to extract the form of the degree distribution analytically.

The directed graph gives rise to different kinds of vertices, shown in Fig.(2), which allow different classes of clustering coefficients $`C`$. Let $`k_{\mathrm{out}}(i)`$, $`k_{\mathrm{in}}(i)`$, $`k(i)`$, be, respectively, the out- degree, in-degree and the total degree of the vertex $`i`$. The clustering coefficient at a given vertex $`i`$ is defined as

$$C_i=\frac{2E(i)}{k(i)[k(i)1]}$$
(10)
where $`E(i)`$ is the number of edges connecting the nearest neighbors of $`i`$. Thus it is the number of pairs of nearest neighbors directly connected to each other, normalized by the largest number of such connections possible. We may extend this concept to directed graphs and define,

$$C_{\mathrm{out}}(i)=\frac{2E_{\mathrm{out}}(i)}{k_{\mathrm{out}}(i)[k_{\mathrm{out}}(i)1]}$$
(11)
and similiarly for $`C_{\mathrm{in}}`$, where $`E_{\mathrm{out}}(i)`$ (respectively $`E_{\mathrm{in}}(i)`$) are the number of edges connecting out (in) nearest neighbors of $`i`$. Note that any pair of incoming and outgoing bonds passing through $`i`$ necessarily defines a triangle, due to transitivity, as shown in Fig.( 2). Thus, the clustering coefficient may be conveniently decomposed as

$$C_i=\frac{k_{\mathrm{in}}(i)k_{\mathrm{out}}(i)+E_{\mathrm{in}}(i)+E_{\mathrm{out}}(i)}{k(i)[k(i)1]/2}.$$
(12)
This is a null model in the sense that no assumptions have been made as to the fitness of any particular type of interaction; the resulting interactions depend only on the random sequences coded in the genes and on the distribution of gene lengths. This random network provides the “tabula rasa” on which we assume natural selection will subsequently act. As such, it is of great interest to determine the properties of the null network, which turn out to be highly non-trivial. This is the task to which we turn in the next section.

## III Simulation results

To characterize the network defined by our model, we generated random chromosomes as defined above. The statistical properties of the network obtained from the totally random chromosome were checked to be invariant under a constant mutational load, with a mutation probability of 0.01. A mutation is affected as follows. If the symbol ($`x`$) occupying the site to be mutated happens to be a 0 or 1, then it is flipped, i.e., we set $`x=\mathrm{mod}_2(x+1)`$. If $`x=2`$, then it exchanges places with either its right or left nearest neighbor, with equal probability. It should be noted that the first case corresponds to a substitutional mutation, whereas the second to a shifting of the position of the start sign which shifts the reading frame of the gene. The latter gives rise to far-reaching modifications, since the three-letter codes corresponding to the different amino acids will be completely modified along the whole gene if the reading frame is shifted.

The most remarkable property of the random network generated in this way is that it has a scale free out-degree distribution, and a qualitatively different in-degree distribution which is much less broad, with a pronounced narrow peak followed by a Gaussian peak. These qualitative behaviours match the results found on protein and genomic networks to a surprising extent. Maslov1 ; Maslov2

We display in Figs.(3,4) the out-degree distribution $`n(k_{\mathrm{out}})`$ for two sets of parameters, namely $`L=5000,p=0.1`$ and for $`L=15000p=0.05`$, with, the number of genes fluctuating around $`\overline{N}=450`$ and $`\overline{N}=712`$, respectively. The data have been averaged over 500 independent realizations. The distributions, which are strongly log periodic, with a power law envelope $`n_m(k_{\mathrm{out}})k_{\mathrm{out}}^\gamma ^{}`$, have the same characteristic behaviour for the two cases. The fits to the envelope of the peaks are shown in Figs.(5,6). We see a marked break in the log-log fit for larger $`p`$, with a crossover from $`\gamma ^{}1`$ to $`\gamma ^{}=0.45\pm 0.06`$ at about $`\mathrm{ln}(k_{\mathrm{out}})=2.5`$. This crossover behavior is just off-scale in Fig.( 6), as can be seen in Fig.( 7), with the incipient scaling with a power of $`1`$ again extending out to $`\mathrm{ln}(k_{\mathrm{out}})2`$.

It is also worthwhile mentioning that the integral over the peaks to the rightmost of the out-degree distribution simply gives the network size $`N`$, as it should, because these peaks come from the smallest genes of unit length which are embedded in all the other genes. Introducing a larger lower cutoff to the gene size (i.e., a cutoff bigger than unity) would somewhat restrict the range of $`k_{\mathrm{out}}`$, but also yield cleaner scaling behaviour for large values of $`k_{\mathrm{out}}`$ by effectively eliminating these peaks.

It should be noticed that for $`k_{\mathrm{out}}100`$, the dips in the degree distribution are about evenly spaced on the logarithmic scale, with a scale factor of $`1/2`$. We conjecture that each individual peak corresponds to a particular gene-length $`\mathrm{}`$. (This conjecture is borne out by comparing the integrals under the peaks with $`n_{\mathrm{}}`$ for $`p=0.1,L=5000`$ but is much more noisy for $`p=0.05,L=15000`$). Increasing $`\mathrm{}`$ by unity, exactly halves the leading contribution to the to $`k_{\mathrm{out}}`$, from genes of length $`\mathrm{}`$.

To get a better grip on the degree distribution, we numerically integrated the curves in Figs.(3,4). The result is shown for $`p=0.05`$ in Fig.(8). We find that

$$𝒩(k_{\mathrm{out}})^{k_{\mathrm{out}}}n(z)𝑑z\mathrm{ln}(k_{\mathrm{out}}),$$
(13)
or a very small power, $`k_{\mathrm{out}}^{\gamma +1}`$ where numerically $`\gamma +1<0.1`$, at least for $`k_{\mathrm{out}}`$ that are not too large. Thus we are led to believe that the overall out degree distribution scales like

$$n(k_{\mathrm{out}})k_{\mathrm{out}}^\gamma f(k_{\mathrm{out}}/\overline{N}),$$
(14)
with $`\gamma =1.0\pm 0.1`$. The scaling function $`f(x)\mathrm{const}.`$ for $`x1`$, and is log-periodic for intermediate values of its argument. Scaling breaks down for $`xO(1)`$. Note that $`\gamma `$ is smaller in absolute value than that reported for real gene regulation networks Sole ; Maslov1 ; Maslov2 , thus evolution seems to have narrowed the distribution somewhat, restricting the number of genes that may be affected by any one gene.

The in-degree distribution, which we display in Fig.(9) shows two peaks. The second can be fit to a Gaussian, as shown in Fig.(10). The first peak is more skewed than a Poissonian, but may be fitted reasonably well by a distribution of the form $`f(x)(xx_0)^\delta \mathrm{exp}[\xi (xx_0)]`$, where $`\delta =2.2`$ and $`\xi =0.15`$. (Fig.(11)).

The total degree distribution again displays a modulated structure, as can be seen in Fig.(12).

We have already discussed that the network will behave trivially either for very small $`N`$ (all very large genes with very small probabilities for interaction) or for large $`p`$ (i.e., $`p>1/2`$) where most of the chromosome will be occupied by non-coding partitions (the symbol 2), and many null genes. We have determined the threshold value, $`p_c`$, below which the the genes are too few and too long, such that the single giant cluster breaks up into more than one component. The results are reported in Table I, for $`L`$ ranging from 15000 to 1000, and undirected edges. Although we have only a few points, the dependence (see Fig. 13) fits a power law, with $`p_cL^\alpha `$, with $`\alpha 3/4`$, within the range explored. In the limit of $`L\mathrm{}`$, $`p_c0`$, but less fast than $`1/L`$.

On the other hand our expectations that there is a reasonably large interval of $`p`$’s ($`p_c<p<1/2`$) where the results do not depend quantitatively on the precise value of $`p`$ are borne out, as can be seen from our calculations for the minimum average path length for undirected edges.

We argued in section II that as long as there is a single giant cluster, the directed minimum path length is identically 1. For undirected paths this is no longer true and we must determine $`l_{\mathrm{min}}`$ numerically. Nevertheless, we can show that $`l_{\mathrm{min}}4`$. Note that for $`p>p_c`$, there will be an abundance of short genes of unit length, which will contain the symbol 0 or 1 with equal probability. These will have outgoing edges to all the typical genes which have an admixture of $`0`$s and $`1`$s, so that most such genes will be certainly linked by paths of length at most equal to 2. The “worst case” is that of atypical cases of all 1’s all or all 0’s which require at least one intermediate typical gene to link up ($`111111011100000000`$), which gives $`\mathrm{max}(l_{\mathrm{min}})=4`$. We have calculated $`l_{\mathrm{min}}`$ for fixed $`p`$ and different $`L`$. If we consider undirected edges and find that it at most depends very weakly on $`N`$, the number of vertices. We find, e.g., for $`L`$ ranging from 1650 to 15000, and $`p=0.05`$ that $`l_{\mathrm{min}}`$ ranges only from 1.673 to 1.699. (Table II)

Fixing the length $`L`$ and varying $`p`$ (Table III) shows again that as long as $`l_{\mathrm{min}}`$ is defined, i.e., just above the “percolation threshold” it is already very close to the value it will have for larger $`p`$.

The clustering coefficients are markedly different for the in- and out-degree distributions, as well as the total connectivity. We find, for $`p=0.05`$ and $`L=15000`$, $`<C_{\mathrm{out}}>=0.034`$, $`<C_{\mathrm{in}}>=0.648`$ and $`<C>=0.534`$. We have verified in the last case that this clustering coefficient is reproduced by computing the average total degree per node and deviding by the total number of nodes, i.e., $`<C><k>/N`$. Thus, the total connectivity behaves very much like in classical random graphs.

## IV Analytical results on toy model

In this section we would like to present a number of the results which may be obtained from a hierarchical model. In particular, we would like to exploit the transitivity property, to compute certain properties of a toy version of our original model.

Let us consider a tree network, with a branching ratio $`b`$. Take all connections emanating from nodes at the $`m`$th generation to the ones at the $`m+1`$st generation, to be directed. If one applies the transitivity rule, then this automatically generates further connections that link directly a node at the $`m`$th generation to all the nodes below it on the same branch of the tree. Then, the out degree (we will drop the index “out” from now on) of a node at the $`m`$th level will be

$$k(m)=\underset{m^{}=0}{\overset{Mm1}{}}b^m^{}=\frac{b^{Mm}1}{b1},$$
(15)
where $`M`$ is the total number of levels and $`k(m)b^{Mm}`$ for large $`Mm`$, which we may safely assume. The frequency of the nodes at the $`m`$th level is $`b^m`$. This gives, for the out-degree distribution,

$$n(k)\mathrm{const}.\left(\frac{1}{k}\right).$$
(16)
We now introduce redundancies within the same level. Let a fraction $`r`$ of the $`b`$ downstream neighbors emanating from any node be identical to each other, with $`r=a^m`$, $`b^1<a<1`$ being some constant. The nodes that are identical within any given level are connected to each other by two-way bonds within the sub-branch where they are located. I will call this subset of nodes the clones at any given level. This means that the out-degree $`k`$ of any one of these nodes is now $`k^{}=(a^mb)b^{Mm}`$, since each interconnected node inherits all the downstream neighbors of the $`a^mb`$ clones to which it is connected. The interconnections between nodes further downstream do not introduce any further change in $`k^{}(m)`$, since all the nodes that they can connect are already connected to the nodes at level $`m`$. The number of clone-nodes (not all identical to each other, but only in groups of $`a^mb`$) is clearly $`a^mb^m`$. Then, the out-degree distribution of the nodes belonging to the clone-sets becomes,

$$n(k^{})=b^{\gamma 1}\left(\frac{k^{}}{N}\right)^\gamma $$
(17)
where $`N=b^M`$ and with

$$\gamma =\frac{\mathrm{ln}a+\mathrm{ln}b}{\mathrm{ln}a\mathrm{ln}b}.$$
(18)
Note that under our assumptions above, $`1<\gamma <0`$. For large enough $`k`$ this then becomes the leading contribution to the out degree distribution, rather than $`\gamma =1`$.

## V Conclusions and Discussion

The model we have presented for the scaling behaviour of gene interaction, more specifically RNA interference, turns out to be a rich model which is worthwhile to investigate in its own right, as well as providing a highly structured starting point on which further evolutionary pressures could act.

We have checked that the statistical properties of the network are robust under random point mutations. Thus, the scale free random network may be seen as the attrator for a process starting from, say, a sequence of uniform genes, subjected to random point mutations over a very long time period. The scale free steady state can be considered as the outcome of a process of self-organization under random perturbations.

In conclusion it is worthwhile to note that in RNA interference, longer RNA chains are cut up by the “Dicer” enzyme into smaller segments, siRNA Hannon , which can be conveniently matched with complementary sequences in a larger number of mRNAs. This trick would change the power of the out-degree distribution from that found here.

Acknowledgements One of us (AE) would like to gratefully acknowledge partial support from the Turkish Academy of Sciences. It is a pleasure to thank B. Bollobas and J.-P. Eckmann for useful conversations.

APPENDIX We would like to compute the probabilities $`P(\mathrm{},\mathrm{}^{})`$ in Eq.(7). For successive $`\nu =\mathrm{}^{}\mathrm{}>0`$, it is convenient to define subsequences of $`G_j`$ of length $`\mathrm{}_i`$ with a shift $`\lambda `$, $`0\lambda \nu `$, such that $`g_{ij}^{(\lambda )}\{x_{j,1+\lambda },x_{j,2+\lambda },\mathrm{}x_{j,\mathrm{}_i+\lambda }\}`$. Then, we can obtain $`P(\mathrm{},\mathrm{}^{})`$ in terms of joint probabilities like $`𝒫(G_ig_{ij}^{(0)},G_i=g_{ij}^{(1)})`$, etc. Thus,

$`P(\mathrm{},\mathrm{})`$ $`=𝒫(G_i=g_{ij}^{(0)})`$
$`P(\mathrm{},\mathrm{}+1)`$ $`=𝒫(G_i=g_{ij}^{(0)})+𝒫(G_ig_{ij}^{(0)},G_i=g_{ij}^{(1)})`$
$`\mathrm{}`$
$`P(\mathrm{},\mathrm{}+\nu )`$ $`=𝒫(G_i=g_{ij}^{(0)})+\mathrm{}`$
$`+`$ $`𝒫(G_ig_{ij}^{(0)},G_ig_{ij}^{(1)},\mathrm{},G_i=G_{ij}^\lambda )`$
$`+`$ $`\mathrm{}`$
$`+`$ $`𝒫(G_ig_{ij}^{(0)},G_ig_{ij}^{(1)},\mathrm{},G_i=g_{ij}^{(\nu )})`$
Let us display $`𝒫(G_ig_{ij}^{(0)},G_i=g_{ij}^{(1)})`$ for greater clarity. We have

$`𝒫(G_ig_{ij}^{(0)},G_i=g_{ij}^{(1)})=`$
$`\left({\displaystyle \underset{\mu =1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mu 1}{}}}\delta _{x_{i,m}x_{j,m}}\right)\left(1\delta _{x_{i,\mu },x_{j,\mu }}\right){\displaystyle \underset{n=1}{\overset{\mathrm{}_i}{}}}\delta _{x_{i,n}x_{j,n+1}}`$
$`=2^\mathrm{}_i\left[12^\mathrm{}_i\right].`$
In particular we find

$`P(\mathrm{},\mathrm{}+1)`$ $`=`$ $`2^{\mathrm{}}\left[1+12^{\mathrm{}}\right]`$
$`=`$ $`2^{\mathrm{}}\left[22^{\mathrm{}}\right]`$
$`P(\mathrm{},\mathrm{}+2)`$ $`=`$ $`2^{\mathrm{}}\left[32^{\mathrm{}+1}{\displaystyle \frac{1}{3}}\left(12^2\mathrm{}\right)\right]`$
Note that the probabilities $`𝒫`$ do not depend upon the particular gene, except through the lengths of the sequences, as the genes do not overlap, and each of the symbols are independently and identically distributed, with equal probability, over 0 and 1.


format: text

text:
# Untitled Document

A Multiple Sum Involving the Möbius Function

Yoichi Motohashi

Abstract. We consider a multiple arithmetical sum involving the Möbius function which despite its elementary appearance is in fact of a highly intriguing nature. We establish an asymptotic formula for the quadruple case that raises the first genuinely non-trivial situation. This is a rework of an old unpublished note of ours. 2001 Mathematics Subject Classification: Primary 11A25; Secondary 11M06

1. Introduction. The aim of the present article is to discuss the asymptotics of the quantity

$$_k(z)=\underset{d_1z}{}\mathrm{}\underset{d_{2k}z}{}\frac{\mu (d_1)\mathrm{}\mu (d_{2k})}{[d_1,\mathrm{},d_{2k}]},k1,$$
$`(1.1)`$
as $`z`$ tends to infinity, where $`\mu `$ is the Möbius function and $`[d_1,\mathrm{},d_{2k}]`$ the least common multiple of positive integers $`d_1,\mathrm{},d_{2k}`$. We have the relation

$$_k(z)=\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{nN}{}\left(\underset{d|n,dz}{}\mu (d)\right)^{2k}.$$
$`(1.2)`$
Thus our problem is pertinent to the extremal behaviour of the truncated sum of the Möbius function over divisors, and somewhat remotely to the Selberg sieve (see the concluding remark). The case $`k=1`$ is treated in . The case $`k=2`$ is already quite involved and discussed in . The result and the outline of the argument there have been shown on a few occasions, first at the Problem Session of the Amalfi International Symposium on Analytic Number Theory, September 1989. There is, however, a minor error in , as is to be indicated below. We have reworked  because of its apparent relation with the recent important challenge  by Goldston and Yıldırım to the problem of finding small gaps between consecutive primes. Their argument depends on their previous work  which deals with higher correlations of short sums of

$$\mathrm{\Lambda }_z(n)=\underset{d|n,dz}{}\mu (d)\mathrm{log}(z/d),$$
$`(1.3)`$
an approximation to the von Mangold function. In their discussion it is needed, among other things, to study the asymptotic behaviour of the sum

$$\underset{nN}{}\mathrm{\Lambda }_z(n+j_1)\mathrm{\Lambda }_z(n+j_2)\mathrm{}\mathrm{\Lambda }_z(n+j_r),$$
$`(1.4)`$
where $`j_1,\mathrm{},j_r`$ are arbitrary non-negative integers. Expanding this via the definition $`(1.3)`$, we are led to an expression closely resembles to $`(1.1)`$. Goldston and Yıldırım applied to this expression an argument essentially the same as that of , apparently without being aware of our old unpublished work.

It appears to us, however, that their problem is less delicate than ours, as far as the handling of the relevant residue calculus is concerned. The factor $`\mathrm{log}(z/d)`$ makes their expression smoother than ours. Being translated into our situation, this is equivalent to having $`(s_1\mathrm{}s_{2k})^2`$ in place of the denominator $`s_1\mathrm{}s_{2k}`$ in $`(1.5)`$ below. Hence, both the convergence and the estimation issues are less troublesome with $`(1.4)`$, although the arithmetical issue can be highly involved when $`j_1,\mathrm{},j_r`$ are arbitrary. The argument of  starts with the following integral expression: For non-integral $`z`$

$$_k(z)=\frac{1}{(2\pi i)^{2k}}\mathrm{}M(s_1,\mathrm{},s_{2k})z^{s_1+\mathrm{}+s_{2k}}\frac{ds_1\mathrm{}ds_{2k}}{s_1\mathrm{}s_{2k}}$$
$`(1.5)`$
with

$$M(s_1,\mathrm{},s_{2k})=\underset{d_1=1}{\overset{\mathrm{}}{}}\mathrm{}\underset{d_{2k}=1}{\overset{\mathrm{}}{}}\frac{\mu (d_1)\mathrm{}\mu (d_{2k})}{[d_1,\mathrm{},d_{2k}]d_1^{s_1}\mathrm{}d_{2k}^{s_{2k}}},$$
$`(1.6)`$
where all integral are over vertical lines placed in the right half plane. Of course this is not a fully correct expression. We need to use, instead, a truncated version of Perron’s formula, and the vertical segments over which the integrations are performed should be placed in a well-poised way, as we shall show later.

We have, for $`\mathrm{Re}s_j>0`$ ($`j=1,\mathrm{},2k`$), the Euler product expansion

$$M(s_1,\mathrm{},s_{2k})=\underset{p}{}\left(1\frac{1}{p}+\frac{1}{p}\underset{j=1}{\overset{2k}{}}\left(1\frac{1}{p^{s_j}}\right)\right).$$
$`(1.7)`$
Thus

$$M(s_1,\mathrm{},s_{2k})=\frac{\zeta (1+s_{\lambda _1}+\mathrm{}+s_{\lambda _{2a}})}{\zeta (1+s_{\tau _1}+\mathrm{}+s_{\tau _{2b1}})}G(s_1,\mathrm{},s_{2k}),$$
$`(1.8)`$
where $`\zeta `$ is the Riemann zeta-function, and $`1\lambda _1<\mathrm{}<\lambda _{2a}2k`$, $`1\tau _1<\mathrm{}<\tau _{2b1}2k`$ with $`a,b1`$. The function $`G`$ is regular and bounded for $`\mathrm{Re}s_j>c(k)`$ ($`j=1,\mathrm{},2k`$) with a constant $`c(k)>0`$ which could be given explicitly. An appropriate shift of contours, along with I.M. Vinogradov’s zero-free region for $`\zeta (s)`$ (see \[6, Theorem 6.1\]), yields Theorem.As $`z`$ tends to infinity, we have

$$\begin{array}{ccc}\hfill _1(z)& =(1+o(1))\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}\frac{G_1(t)}{|\zeta (1+it)|^2}\frac{dt}{t^2},\hfill & (1.9)\hfill \\ \hfill _2(z)& =(1+o(1))\frac{3}{4\pi }(\mathrm{log}z)^2_{\mathrm{}}^{\mathrm{}}\frac{|\zeta (1+2it)|^2}{|\zeta (1+it)|^8}G_2(t)\frac{dt}{t^4},\hfill & (1.10)\hfill \end{array}$$

where $`G_k(t)=G(s_1,\mathrm{},s_{2k})`$ with $`s_1,\mathrm{},s_k=it`$, $`s_{k+1},\mathrm{},s_{2k}=it`$.  We note that $`G_k(t)>0`$. The formula $`(1.9)`$ is proved in  with an argument different from ours. The formula $`(1.10)`$ is a corrected version of the relevant claim made in ; there was an error in the computation of certain residues. The advantage of our argument over that of  is perhaps in that ours can give rise to $`(1.10)`$. Our argument should work, in principle, for any $`k`$. However, the mode of shifts of contours and the arrangement of residues become formidably complicated for $`k3`$. Thus the general case will probably require a new approach, though the case $`k=3`$ appears to be still manageable as a direct extension of the present work.

2. Proof of Theorem. We shall deal with the case $`k=2`$ only, for the case $`k=1`$ is analogous and in fact by far simpler. Also we shall assume that $`z`$ is half a large odd integer. Obviously this will make no difference. To begin with, let $`\alpha =(\mathrm{log}z)^1`$ and $`T=z^2`$. Then we have

$$\begin{array}{ccc}\hfill _2(z)=\frac{1}{(2\pi i)^4}& _{10\alpha 10Ti}^{10\alpha +10Ti}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}M(s_1,s_2,s_3,s_4)\hfill & \\ & \times z^{s_1+s_2+s_3+s_4}\frac{ds_1ds_2ds_3ds_4}{s_1s_2s_3s_4}+O\left(\frac{(\mathrm{log}z)^7}{z}\right),\hfill & (2.1)\hfill \end{array}$$

where the implied constant is absolute. To show this, we note first that for any positive integer $`d_1`$

$$\frac{1}{2\pi i}_{\alpha Ti}^{\alpha +Ti}\left(\frac{z}{d_1}\right)^{s_1}\frac{ds_1}{s_1}=\delta (d_1)+O\left(\frac{z}{Td_1^\alpha }\right),$$
$`(2.2)`$
where $`\delta (d)=1`$ if $`d<z`$ and $`0`$ otherwise. This is of course a crude consequence of Perron’s inversion formula. Multiply both sides by the factor $`(z/d_2)^{s_2}/s_2`$ with an integer $`d_2>0`$ and integrate with respect to $`s_2`$ as indicated by $`(2.1)`$. We have

$$\begin{array}{ccc}\hfill \frac{1}{(2\pi i)^2}& _{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{z^{s_1+s_2}}{d_1^{s_1}d_2^{s_2}}\frac{ds_1ds_2}{s_1s_2}\hfill & \\ & =\delta (d_1)\delta (d_2)+O\left(\frac{z\delta (d_1)}{Td_2^\alpha }\right)+O\left(\frac{z\mathrm{log}T}{T(d_1d_2)^\alpha }\right)\hfill & \\ & =\delta (d_1)\delta (d_2)+O\left(\frac{z\mathrm{log}T}{T(d_1d_2)^\alpha }\right).\hfill & (2.3)\hfill \end{array}$$

Repeating the same procedure, we get

$$\begin{array}{ccc}\hfill \frac{1}{(2\pi i)^4}& _{10\alpha 10Ti}^{10\alpha +10Ti}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{z^{s_1+s_2+s_3+s_4}}{d_1^{s_1}d_2^{s_2}d_3^{s_3}d_4^{s_4}}\frac{ds_1ds_2ds_3ds_4}{s_1s_2s_3s_4}\hfill & \\ & =\underset{j=1}{\overset{4}{}}\delta (d_j)+O\left(\frac{z(\mathrm{log}T)^3}{T(d_1d_2d_3d_4)^\alpha }\right).\hfill & (2.4)\hfill \end{array}$$

Then, we divide both sides by $`[d_1,d_2,d_3,d_4]`$ and sum the result. We find that the first term on the right of $`(2.1)`$ is equal to

$$\begin{array}{ccc}& _2(z)+O\left(\frac{z(\mathrm{log}T)^3}{T}\underset{d_1=1}{\overset{\mathrm{}}{}}\underset{d_2=1}{\overset{\mathrm{}}{}}\underset{d_3=1}{\overset{\mathrm{}}{}}\underset{d_4=1}{\overset{\mathrm{}}{}}\frac{|\mu (d_1)\mu (d_2)\mu (d_3)\mu (d_4)|}{[d_1,d_2,d_3,d_4](d_1d_2d_3d_4)^\alpha }\right)\hfill & \\ \hfill =& _2(z)+O\left(\frac{z(\mathrm{log}T)^3}{T}\underset{p}{}\left(1\frac{1}{p}+\frac{1}{p}\left(1+\frac{1}{p^\alpha }\right)^4\right)\right).\hfill & (2.5)\hfill \end{array}$$

Observing that this Euler product is $`O(\zeta ^4(1+\alpha ))`$, we end the proof of $`(2.1)`$. Now, let $`\beta =(\mathrm{log}z)^{3/4}`$. We shift the contour for the $`s_4`$-integral to the vertical segment $`[\beta 10Ti,\beta +10Ti]`$. In view of $`(1.8)`$, we encounter poles at $`s_4=s_1,s_2,s_3`$, and $`(s_1+s_2+s_3)`$. Computing respective residues, we have

$$_2(z)=\left\{_2^{(1)}+_2^{(2)}+_2^{(3)}+_2^{(4)}\right\}(z)+O\left(z^{\beta /2}\right).$$
$`(2.6)`$
by virtue of Vinogradov’s zero-free region for $`\zeta `$ together with the related bounds for $`\zeta `$, $`1/\zeta `$ (see \[6, Theorem 6.3; Lemma 12.3\]); the same combination will be implicitly invoked in what follows as well. Here

$$\begin{array}{ccc}\hfill _2^{(1)}(z)=& \frac{1}{(2\pi i)^3}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)\zeta (1+s_1+s_3)}{\zeta (1+s_1)\zeta (1+s_2)^2\zeta (1+s_3)^2}\hfill & \\ & \times \frac{\zeta (1+s_2+s_3)^2\zeta (1+s_2s_1)\zeta (1+s_3s_1)}{\zeta (1s_1)\zeta (1+s_1+s_2+s_3)\zeta (1+s_2+s_3s_1)}\hfill & \\ & \times G(s_1,s_2,s_3,s_1)z^{s_2+s_3}\frac{ds_1ds_2ds_3}{s_1^2s_2s_3},\hfill & (2.7)\hfill \\ \hfill _2^{(2)}(z)=& \frac{1}{(2\pi i)^3}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)\zeta (1+s_1+s_3)^2}{\zeta (1+s_1)^2\zeta (1+s_2)\zeta (1+s_3)^2},\hfill & \\ & \times \frac{\zeta (1+s_2+s_3)\zeta (1+s_1s_2)\zeta (1+s_3s_2)}{\zeta (1s_2)\zeta (1+s_1+s_2+s_3)\zeta (1+s_1+s_3s_2)}\hfill & \\ & \times G(s_1,s_2,s_3,s_2)z^{s_1+s_3}\frac{ds_1ds_2ds_3}{s_1s_2^2s_3},\hfill & (2.8)\hfill \\ \hfill _2^{(3)}(z)=& \frac{1}{(2\pi i)^3}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)^2\zeta (1+s_1+s_3)}{\zeta (1+s_1)^2\zeta (1+s_2)^2\zeta (1+s_3)}\hfill & \\ & \times \frac{\zeta (1+s_2+s_3)\zeta (1+s_1s_3)\zeta (1+s_2s_3)}{\zeta (1s_3)\zeta (1+s_1+s_2+s_3)\zeta (1+s_1+s_2s_3)}\hfill & \\ & \times G(s_1,s_2,s_3,s_3)z^{s_1+s_2}\frac{ds_1ds_2ds_3}{s_1s_2s_3^2},\hfill & (2.9)\hfill \\ \hfill _2^{(4)}(z)=& \frac{1}{(2\pi i)^3}_{4\alpha 4Ti}^{4\alpha +4Ti}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)\zeta (1s_1s_2)}{\zeta (1+s_1)\zeta (1s_1)\zeta (1+s_2)}\hfill & \\ & \times \frac{\zeta (1+s_1+s_3)\zeta (1s_1s_3)\zeta (1+s_2+s_3)\zeta (1s_2s_3)}{\zeta (1s_2)\zeta (1+s_3)\zeta (1s_3)\zeta (1+s_1+s_2+s_3)\zeta (1s_1s_2s_3)}\hfill & \\ & \times G(s_1,s_2,s_3,s_1s_2s_3)\frac{ds_1ds_2ds_3}{s_1s_2s_3(s_1+s_2+s_3)}.\hfill & (2.10)\hfill \end{array}$$

Let us first show that

$$_2^{(4)}(z)(\mathrm{log}z)^{3/2}.$$
$`(2.11)`$
We note that the bound $`_2^{(4)}(z)1`$ appears highly probable; in fact, this holds under the Riemann Hypothesis. To prove $`(2.11)`$ we observe first that

$$_2^{(4)}(z)=\frac{1}{(2\pi i)^3}_{\beta 4Ti}^{\beta +4Ti}_{\beta 2Ti}^{\beta +2Ti}_{\beta Ti}^{\beta +Ti}\{\mathrm{}\}\frac{ds_1ds_2ds_3}{s_1s_2s_3(s_1+s_2+s_3)}+o(1).$$
$`(2.12)`$
Then we shift the contour of the inner-most integral to

$$C=\{s_1:\frac{c}{(\mathrm{log}(2+|t|+|s_2|+|s_3|))^{3/4}}+it,TtT\},$$
$`(2.13)`$
where $`c>0`$ needs to be sufficiently small. We have

$$_2^{(4)}(z)=\frac{1}{(2\pi i)^3}_{\beta 4Ti}^{\beta +4Ti}_{\beta 2Ti}^{\beta +2Ti}_C\{\mathrm{}\}\frac{ds_1ds_2ds_3}{s_1s_2s_3(s_1+s_2+s_3)}+o(1).$$
$`(2.14)`$
This implies that

$$\begin{array}{ccc}\hfill _2^{(4)}(z)& (\mathrm{log}z)^{3/2}_{4T}^{4T}_{2T}^{2T}_T^T\mathrm{log}^{10}(2+|t_1|+|t_2|+|t_3|)\hfill & \\ & \times \frac{dt_1dt_2dt_3}{(1+|t_1|)(1+|t_2|)(1+|t_3|)(1+|t_1+t_2+t_3|)}+o(1).\hfill & (2.15)\hfill \end{array}$$

On the right side, the factor $`(\mathrm{log}z)^{3/2}`$ comes from the factor $`\zeta (1+s_2+s_3)\zeta (1s_2s_3)`$ in $`(2.14)`$. In the integrand, the denominator comes from that in $`(2.14)`$, and the logarithmic factor from those zeta-factors there, save for $`\zeta (1+s_2+s_3)\zeta (1s_2s_3)`$. One may see readily that

$$_T^T\frac{\mathrm{log}^{10}(2+|t_1|+|t_2|+|t_3|)}{(1+|t_1|)(1+|t_1+t_2+t_3|)}𝑑t_1\frac{\mathrm{log}^{11}(2+|t_2|+|t_3|)}{1+|t_2+t_3|}.$$
$`(2.16)`$
Thus we have

$$_2^{(4)}(z)(\mathrm{log}z)^{3/2}_{4T}^{4T}\frac{\mathrm{log}^{12}(2+|t_3|)}{(1+|t_3|)^2}𝑑t_3,$$
$`(2.17)`$
which proves our claim $`(2.11)`$. Let us treat $`_2^{(1)}`$. In $`(2.7)`$ we shift the $`s_3`$-contour to the segment $`[\beta 4Ti,\beta +4Ti]`$. We encounter poles at $`s_3=s_1,s_1,s_2`$. Computing the respective residues we get

$$_2^{(1)}(z)=\left\{_2^{(1,1)}+_2^{(1,2)}+_2^{(1,3)}\right\}(z)+O(z^{\beta /2}).$$
$`(2.18)`$
We have

$$_2^{(1,3)}(z)\mathrm{log}z.$$
$`(2.19)`$
In fact, $`_2^{(1,3)}(z)`$ is a linear polynomial in $`\mathrm{log}z`$, whose coefficients are bounded. More precisely, the leading coefficient is equal to

$$\begin{array}{ccc}\hfill \frac{1}{(2\pi i)^2}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}& \frac{\zeta (1+s_1+s_2)\zeta (1s_1s_2)\zeta (1+s_1s_2)\zeta (1s_1+s_2)}{(\zeta (1+s_1)\zeta (1s_1)\zeta (1+s_2)\zeta (1s_2))^2}\hfill & \\ & \times G(s_1,s_2,s_2,s_1)\frac{ds_1ds_2}{(s_1s_2)^2}.\hfill & (2.20)\hfill \end{array}$$

We shift the contour of the outer integral to

$$\{s_2:\frac{c}{(\mathrm{log}(2+|t|+|s_1|))^{3/4}}+it,2Tt2T\},$$
$`(2.21)`$
with a small $`c>0`$. We do not encounter any pole. The new double integral is bounded by a constant multiple of

$$_{2T}^{2T}_T^T\frac{\mathrm{log}^{12}(2+|t_1|+|t_2|)}{((1+|t_1|)(1+|t_2|))^2}𝑑t_1𝑑t_21,$$
$`(2.22)`$
as claimed. The constant term of the linear polynomial has more complicated expression than $`(2.20)`$, involving derivatives of the zeta-function. However, its treatment is analogous.

On the other hand we have

$$\begin{array}{ccc}\hfill _2^{(1,1)}(z)=& \frac{1}{(2\pi i)^2}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)^3\zeta (1+2s_1)}{\zeta (1+s_1)^3\zeta (1+s_2)^3}\hfill & \\ & \times \frac{\zeta (1+s_2s_1)}{\zeta (1s_1)\zeta (1+2s_1+s_2)}G(s_1,s_2,s_1,s_1)z^{s_1+s_2}\frac{ds_1ds_2}{s_1^3s_2},\hfill & (2.23)\hfill \\ \hfill _2^{(1,2)}(z)=& \frac{1}{(2\pi i)^2}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1s_1+s_2)^3\zeta (12s_1)}{\zeta (1s_1)^3\zeta (1+s_2)^3}\hfill & \\ & \times \frac{\zeta (1+s_2+s_1)}{\zeta (1+s_1)\zeta (12s_1+s_2)}G(s_1,s_2,s_1,s_1)z^{s_2s_1}\frac{ds_1ds_2}{s_1^3s_2}.\hfill & (2.24)\hfill \end{array}$$

In the latter we shift the $`s_1`$-contour to the segment $`[\alpha iT,\alpha +iT]`$. We do not encounter any pole. In the new $`s_1`$-integral we perform the change of variable $`s_1s_1`$. On noting $`G(s_1,s_2,s_1,s_1)=G(s_1,s_2,s_1,s_1)`$, we have

$$\left\{_2^{(1,1)}+_2^{(1,2)}\right\}(z)=2_2^{(1,1)}(z)+o(1).$$
$`(2.25)`$
We then shift the $`s_2`$-contour in $`(2.23)`$ to the segment $`[\beta 2iT,\beta +2iT]`$. We encounter poles at $`s_2=s_1,s_1`$, with the resulting double integral being $`O(z^{\beta /2})`$. The first pole contributes

$$\frac{1}{2\pi i}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+2s_1)^4G(s_1,s_1,s_1,s_1)z^{2s_1}}{\zeta (1+s_1)^6\zeta (1s_1)\zeta (1+3s_1)}\frac{ds_1}{s_1^4},$$
$`(2.26)`$
which is obviously $`O(z^{\beta /2})`$. Thus, computing the residue at $`s_2=s_1`$, we have

$$\begin{array}{ccc}& \left\{_2^{(1,1)}+_2^{(1,2)}\right\}(z)\hfill & \\ & =\frac{(\mathrm{log}z)^2}{2\pi i}_{\alpha iT}^{\alpha +iT}\frac{\zeta (1+2s_1)\zeta (12s_1)}{(\zeta (1+s_1)\zeta (1s_1))^4}G(s_1,s_1,s_1,s_1)\frac{ds_1}{s_1^4}+O(\mathrm{log}z).\hfill & (2.27)\hfill \end{array}$$

This error term is actually equal to a negligible term plus a linear polynomial of $`\mathrm{log}z`$, the coefficients of which are easily seen to be bounded.

From $`(2.18)`$, $`(2.19)`$ and $`(2.27)`$ we obtain

$$_2^{(1)}(z)=(1+o(1))\frac{1}{2\pi }(\mathrm{log}z)^2_{\mathrm{}}^{\mathrm{}}\frac{|\zeta (1+2it)|^2}{|\zeta (1+it)|^8}G(it,it,it,it)\frac{dt}{t^4}.$$
$`(2.28)`$
which ends our computation of $`_2^{(1)}(z)`$. Next, we shall consider $`_2^{(2)}`$; we may be brief. In $`(2.8)`$ we shift the $`s_3`$-contour to $`[\beta 4iT,\beta +4iT]`$. We encounter poles at $`s_3=s_2,s_2,s_1`$. Computing the respective residues, we have

$$_2^{(2)}(z)=\left\{_2^{(2,1)}+_2^{(2,2)}+_2^{(2,3)}\right\}(z)+O(z^{\beta /2}).$$
$`(2.29)`$
We have $`_2^{(2,3)}(z)\mathrm{log}z`$ similarly to $`(2.19)`$. We have

$$\begin{array}{ccc}\hfill _2^{(2,1)}(z)=& \frac{1}{(2\pi i)^2}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1+s_2)^3\zeta (1+2s_2)}{\zeta (1+s_1)^3\zeta (1+s_2)^3},\hfill & \\ & \times \frac{\zeta (1+s_1s_2)}{\zeta (1s_2)\zeta (1+s_1+2s_2)}G(s_1,s_2,s_2,s_2)z^{s_1+s_2}\frac{ds_1ds_2}{s_1s_2^3},\hfill & (2.30)\hfill \\ \hfill _2^{(2,2)}(z)=& \frac{1}{(2\pi i)^2}_{2\alpha 2Ti}^{2\alpha +2Ti}_{\alpha Ti}^{\alpha +Ti}\frac{\zeta (1+s_1s_2)^3\zeta (12s_2)}{\zeta (1+s_1)^3\zeta (1s_2)^3},\hfill & \\ & \times \frac{\zeta (1+s_1+s_2)}{\zeta (1+s_2)\zeta (1+s_12s_2)}G(s_1,s_2,s_2,s_2)z^{s_1s_2}\frac{ds_1ds_2}{s_1s_2^3}.\hfill & (2.31)\hfill \end{array}$$

In the latter we shift the $`s_2`$-contour to the segment $`[\beta 2iT,\beta +2iT]`$, and we get $`_2^{(2,2)}(z)z^{\beta /2}`$. On the other hand, in the former we shift the $`s_2`$-contour to $`[\beta 2iT,\beta +2iT]`$. We have

$$_2^{(2,1)}(z)=\frac{(\mathrm{log}z)^2}{4\pi i}_{\alpha iT}^{\alpha +iT}\frac{\zeta (1+2s_1)\zeta (12s_1)}{(\zeta (1+s_1)\zeta (1s_1))^4}G(s_1,s_1,s_1,s_1)\frac{ds_1}{s_1^4}+O(\mathrm{log}z).$$
$`(2.32)`$
Hence we obtain

$$_2^{(2)}(z)=(1+o(1))\frac{1}{4\pi }(\mathrm{log}z)^2_{\mathrm{}}^{\mathrm{}}\frac{|\zeta (1+2it)|^2}{|\zeta (1+it)|^8}G(it,it,it,it)\frac{dt}{t^4}.$$
$`(2.33)`$
It now remains for us to consider $`_2^{(3)}`$. This time we shift first the $`s_3`$-contour in $`(2.9)`$ to the segment $`[2\beta 4iT,2\beta +4iT]`$. We do not encounter any pole, and thus $`_2^{(3)}(z)`$ is equal to the new integral plus a negligible error. In the new integral we shift the $`s_2`$-contour to the segment $`[\beta 2iT,\beta +2iT]`$. We encounter only one pole at $`s_2=s_1`$. Computing the residue we get

$$_2^{(3)}(z)\mathrm{log}z.$$
$`(2.34)`$
Finally, collecting $`(2.6)`$, $`(2.11)`$, $`(2.28)`$, $`(2.33)`$ and $`(2.34)`$, we end our proof of $`(1.10)`$.

Remark. There is an old conjecture by P. Erdös about the size of the arithmetic function

$$\underset{z}{sup}\left|\underset{d|n,dz}{}\mu (d)\right|.$$
$`(2.35)`$
See  for details. Our problem is certainly related to the dual of Erdös’; that is, the supremum is taken in $`n`$ instead of $`z`$. As to the possible relation of our problem with the Selberg sieve, see , ,  and , in chronological order. In addition to these, see \[8, §1.3\] for an extension of $`(1.3)`$. It should be noted that ,  and  were developed in conjunction with the zero-density theory for the Riemann zeta- and Dirichlet $`L`$-functions. However, for this particular purpose those works turned out later to be redundant due to the observation \[8, (1.3.12)\].

References

 M.B. Barban and P.P. Vehov. An extremal problem. Trans. Moscow Math. Soc., 18 (1968), 91-99.

 F. Dress, H. Iwaniec and G. Tennenbaum. Sur une somme liée à la fonction Möbius. J. reine angew. Math., 340 (1983), 53-58.

 S. Graham. An asymptotic estimate related to Selberg’s sieve. J. Number Theory, 10 (1978), 83-94.

 D. Goldston and C.M. Yıldırım. Higher correlations of divisor sums related to primes. III. Preprint, September 2002.

 —. Small gaps between primes. Notes, March 2003.

 A. Ivić. The Riemann Zeta-Function. Theory and Applications. Dover Publ., Inc., Mineola, New York 2003.

 M. Jutila. On a problem of Barban and Vehov. Mathematika, 26 (1979), 62-71.

 Y. Motohashi. A problem in the theory of sieves. Kokyuroku RIMS Kyoto Univ., 222 (1974), 9-50.

 —. Sieve Methods and Prime Number Theory. Tata IFR, Lect. Math. Phys., 72, Springer Verlag, Berlin etc. 1983.

 —. Möbius function over divisors. Notes, March 1985.

Yoichi Motohashi Honkomagome 5-67-1 #901, Tokyo 113-0021, Japan

Email: am8y-mths@asahi-net.or.jp

HP: www.ne.jp/asahi/zeta/motohashi/


text:
# Two-surface wave decay: improved analytical theory and effects on electron acceleration

## I Introduction

The excitation of electron surface waves (ESWs) is a possible route to collisionless energy absorption in plasmas produced by the interaction of sub-picosecond, high intensity (typically $`10^{16}\text{W/cm}^2`$) laser pulses with solid targets. This regime is relevant for applications such as generation of ultrashort X-ray pulses, either as uncoherent thermal emission (gibbon  and references therein) or high laser harmonics vonderlinde , or production of energetic electrons and ions gibbon ; amiranoff ; pukhov ; umstadter .

Linear mode conversion of the laser wave into an ESW is not possible at a planar plasma-vacuum interface, because the phase matching between the ESW and the incident wave is not allowed; the process can take place for specially tailored density profiles, e.g. grating targets with a density profile modulated at the surface with a wavevector $`k_p`$. In this case, the condition for the excitation of a ESW with wavevector $`k_s`$ (parallel to the plasma surface) and the same frequency $`\omega _L`$ of the laser pulse is

$$k_L\mathrm{sin}\theta =k_s+k_p,$$
(1)
where $`k_L=\omega _L/c`$, $`\theta `$ is the angle of incidence of the laser pulse. Experimental investigations of laser absorption and X-ray emission in grating targets are described, e.g., in Refs.gauthier ; bastiani .

Nonlinear mode conversion, i.e. parametric excitation of two ESWs is possible also in a simple “step” profile. We name this process “two–surface wave decay” (TSWD) Macchi02 . Examples of similar TSWD processes were previously considered in regimes other than intense laser–plasma interactions and in the electrostatic limit only gradov ; stenflo . The general phase matching conditions for such a three–wave process are given by

$$k_0=k_++k_{},\omega _0=\omega _++\omega _{},$$
(2)
where $`k_0,\omega _0`$ and $`k_\pm ,\omega _\pm `$ are the wavevector and frequency of the “pump” wave and the two ESWs, respectively. If the electric field of the laser pulse acts as a “pump” for TSWD, then $`\omega _0=\omega _L`$, $`k_0=k_L\mathrm{sin}\theta `$, and equation (2) implies that two sub–harmonic ESW are generated with frequency around $`\omega /2`$ (“$`\omega \omega /2+\omega /2`$” TSWD). However, also the $`𝐯\times 𝐁`$ term of the Lorentz force may excite TSWD. This was observed in particle-in-cell (PIC) simulations Macchi01  for normal laser incidence; in this case, the $`𝐯\times 𝐁`$ force drives 1D oscillations at the plasma surface with frequency $`2\omega _L`$; after a few laser cycles the overlap of a standing wave with frequency $`\omega _L`$ was observed. This is a clear signature of a “$`2\omega \omega +\omega `$” process leading to two counterpropagating ESWs both having the frequency of the laser pulse.

According to the theoretical model Macchi02 ; MacchiLPB02 , the maximum growth rate of the $`2\omega \omega +\omega `$ case is found for normal incidence. As said above, TSWD does not need a structured target; however, it is worth noticing that in a grating target the wavevector $`k_p`$ of the surface modulation \[eq.(1)\] equals the wavevector of the two ESWs excited by the $`2\omega \omega +\omega `$ process at $`\theta =0`$. Hence, TSWD is enhanced in such a grating target. It is then interesting to notice that PIC simulations reported in Ref.bastiani , showing the generation of two countepropagating ESWs in grating targets at $`10^{16}\text{W/cm}^2`$, can now be interpreted as an evidence of TSWD seeded by the surface grating.

At very high intensities ($`10^{18}\text{W/cm}^2`$), simulations Macchi01  show that TSWD enters a strongly nonlinear regime leading to strong rippling of the plasma surface. This might be relevant to the surface instabilities which have been observed in experiments Tarasevitch  and play a detrimental role in high harmonic generation from solid targets norreys ; zepf ; tarasevitch2 ; such instabilities appear to grow even for pulses of few tens of fs tarasevitch2  and thus must be of electronic nature, since ions do not move on such a fast scale. Simulations also show that nonlinear TSWD affects fast electron generation. This effect is investigated in section III by means of test-particle simulations. Before that, in section II we present an improvement of the model of Ref.Macchi02  where temperature effects have been included. In both cases, we restrict ourselves to the case of the $`2\omega \omega +\omega `$ TSWD at normal laser incidence, and we neglect collisions for simplicity and because collective processes dominate absorption at intensities $`10^{16}\text{W/cm}^2`$.

## II Temperature and boundary effects on TSWD

In Ref.Macchi02  the TSWD growth rate was calculated for a step-like density profile \[$`n_i=n_0\mathrm{\Theta }(x)`$, being $`\mathrm{\Theta }(x)`$ the Heaviside step function\] using an Eulerian, fluid model with immobile ions and using the cold plasma approximation. We adopted the following expansion for all fields

$`f(x,y,t)`$ $`=`$ $`f_i(x,ty\mathrm{sin}\theta /c)+ϵf_0(x)e^{i\omega _0(ty\mathrm{sin}\theta /c)}+`$ (3)
$`ϵ^2[f_+(x)e^{ik_+yi\omega _+t}+f_{}(x)e^{ik_{}yi\omega _{}t}],`$
where $`ϵ`$ is a small expansion parameter, and $`f`$ stands for either the electron density or velocity or for the EM fields in the $`(x,y)`$ plane. The first term ($`f_i`$) of eq.(3) includes zero-order, unperturbed fields or oscillating fields that are non-resonant with the excited modes; the term $`f_0`$ represents the “pump” field at the frequency $`\omega _0`$; the last term is the sum of two counterpropagating surface modes. For the $`2\omega \omega +\omega `$ process at $`\theta =0`$, $`\omega _0=2\omega _L`$ and $`k_+=k_{}`$, and the zero– and first–order fields do not depend on $`y`$ (i.e. they are “1D” fields). The coupling between the pump and the surface modes (of order $`ϵ^3`$) originates from the nonlinear terms in the Euler equation ($`en_e𝐯\times 𝐁`$, $`m_en_e𝐯𝐯`$) and the current density $`𝐉=en_e𝐯`$.

The calculation for a warm, isothermal plasma proceeds very similarly battaglini  and thus only the differences from the “cold” case and their consequences are discussed below, while details of the calculation will be published elsewhere. The only difference in the starting Maxwell-Euler systems of equations comes from the pressure term in the Euler equation for electrons:

$$m_ed𝐯/dt=en_e(𝐄+𝐯\times 𝐁)k_BT_en_e,$$
(4)
where $`n_e`$, $`𝐯`$ and $`T_e`$ are the density, fluid velocity and temperature of electrons.

Assuming reflective conditions at the plasma boundary ($`x=0`$), the charge density $`\rho =e(n_in_e)0`$ only in a surface layer with a thickness on the order of the Debye length. At $`T_e=0`$ this corresponds to a surface charge layer \[i.e. $`\rho \delta (x)`$, where $`\delta (x)`$ is the Dirac delta function\], while the longitudinal velocity ($`v_x`$) is discontinous at the surface. This is relevant because in the $`T_e=0`$ case nonlinear coupling terms involving the product of fields that are singular at $`x=0`$ occur. The pressure term removes such singularities and enables to evaluate all such terms correctly. Therefore, by calculating the TSWD growth rate and then taking the $`T_e0`$ limit one achieves an improvement of the result obtained in the cold plasma case.

The final result of the calculation is the $`2\omega \omega +\omega `$ growth rate $`\mathrm{\Gamma }`$ shown in Fig.1 (the complete analytical expression is very lengthy battaglini  and is not reported here). With respect to the $`T_e=0`$ case, the two divergences of $`\gamma `$ are both quenched for growing $`T_e`$, but for different reasons. The divergence for $`\omega \omega _p/\sqrt{2}`$ is quenched since in this limit the ESWs have shorter wavelengths and are thus more affected by the thermal pressure that inhibits the formation of small-scale structures. The “pump” resonance at $`\omega \omega _p/2`$, due to the excitation of plasmons by the $`𝐯\times 𝐁`$ force, is quenched by energy transport out of the skin layer because for $`T_e0`$ the resonant plasmon propagates into the overdense plasma.

## III Electron acceleration

At high intensity of the laser pulse, it is well known that most of the absorbed energy goes into “fast” electrons injected into the overdense plasma region. At normal incidence, due to the leading action of the $`𝐯\times 𝐁`$ force, fast electron bunches are produced twice per laser cycle. The PIC simulations of Ref.Macchi01  showed that after the growth of nonlinear TSWD the generation of fast electrons was enhanced near spatial maxima of the standing surface oscillation (Fig.2, top). Thus, this effect may give an “imprint” for the formation of electron filaments, whose size and spacing would be close to (and scaling with) the laser wavelength as observed in other simulations lasinski ; sentoku , and affect energy transport by fast electrons into the plasma.

To investigate this effect further, we performed “test particle” simulations of electron motion in the overlapping “pump” and TSWD fields. In other words, we solve the equation of motion for electrons into a force field of the form given by eq.(3), i.e. the sum of an one-dimensional (1D) force of frequency $`2\omega `$ (whose analytical expression is obtained from the theoretical model) and the force field from a standing 2D surface wave of frequency $`\omega `$. A similar study, focused on the acceleration of electrons by a single ESW, was recently reported riconda .

The 1D and 2D fields vary in time as $`\mathrm{cos}2\omega t`$ and $`\mathrm{sin}\omega t`$, so that the temporal maxima and minima of the ESW field are always coincident in time with maxima of the 1D field, as it was found in theory and simulation Macchi02 ; Macchi01 . In what follows, $`x`$ is the direction normal to the plasma surface located at $`x=0`$ (the plasma occupies the $`x>0`$ region) and $`y`$ is the direction of propagation of ESW. The field amplitudes are chosen so that the system remains far from relativistic conditions and the ESW field can be considered as a perturbation. The ratio between the amplitudes of standing SW and the $`2\omega `$ pump was varied as a parameter. For the simulations reported in this paper, in terms of normalized amplitudes, $`a=eE/m_ec\omega `$ where $`E`$ is the amplitude of the electric field, the values are $`a_L=0.2`$ for the laser field and $`a_{ESW}=0.019`$ for the ESW field. The wavevectors of the laser wave and of the ESW are given by the known expressions for a cold plasma with density $`n_e/n_c=\omega _p^2/\omega ^2=5`$. To initialize the simulations, we gave the particles an initial position $`x(0)>0`$ such that the evanescent fields are negligible at that point, and an initial velocity $`v_x(0)<0`$ (on average $`v_x=0.1c`$) so that the particles move towards the surface region, where the ESW fields are localized. The particles are distributed in $`y`$ uniformly over a region of width $`\lambda _s`$, where $`\lambda _s`$ is the ESW wavelength.

The ($`y,p_x`$) projection of the test particle phase space (Fig.2, bottom) looks similar to the one from the fully self-consistent PIC simulation (Fig.2, top), showing that electron acceleration is indeed enhanced near spatial maxima of the standing wave. The overlap of the SW fields with the $`𝐯\times 𝐁`$-driven fields (homogeneous along the $`y`$-direction) thus leads to a modulation along $`y`$ of the longitudinal momenta.

The enhancement (or quenching) of the accelerating field occurs once per laser cycle (because the ESW frequency equals the laser frequency), but with opposite phase between contiguous maxima of the standing ESW. This is shown by Fig.3 where the complete ($`x,p_x`$) phase space projection, showing electron bunches generated at $`2\omega `$ rate, is compared with the same plot but including now only the electrons whose starting position lies within an interval of $`\lambda /4`$ width around the $`y=\lambda /4`$ spatial maximum of the ESW. It is thus evident how the most energetic electrons in the bunches penetrating in the $`x>0`$ region are generated near the ESW maximum and at $`\omega `$ rate, i.e. with the ESW frequency. Taking only the electrons around the $`y=3\lambda /4`$ spatial maximum gives an almost identical picture, except that the most energetic bunches are now out of phase by an angle of $`2\pi /\omega `$ with respect to those coming from $`y=\lambda /4`$.

Although the amplitude of ESW fields is 0.1 times the pump field amplitude, the relative modulation of the longitudinal momentum $`p_x`$ is about 30%. This can be qualitatively explained at follows. The electrons acquire energy from an evanescent, oscillating field if their transit time across the region where the field has non-vanishing amplitude is shorter than an oscillation period; this is the condition for non-adiabaticity of electron motion with respect to the field. The ratio between the transit time and the oscillation period is approximately given by the parameter $`\eta =L/v_0T`$ where $`L`$ is the evanescence length, $`v_0`$ is the average velocity of the electron and $`T`$ is the oscillation period. For the ESW, $`T=2\pi /\omega `$ and $`l_{ESW}=(c/\omega )\sqrt{\alpha 2}/(\alpha 1)`$ (where $`\alpha =\omega _p^2/\omega ^2`$) while for the 1D field at $`T=\pi /\omega `$ and $`l_{2\omega }=(c/2\omega )(1/\sqrt{\alpha 1})`$ and Macchi02 . Thus, $`\eta _{ESW}/\eta _{2\omega }=\sqrt{(\alpha 2)/(\alpha 1)}<1`$ which means that the electron motion (for a given $`v_0`$) is more non-adiabatic with respect to the ESW field rather than to the $`2\omega `$ oscillation. Thus, the contribution of the ESW in accelerating (or decelerating) electrons is enhanced by the relatively low frequency and short evanescence length.

The “imprint” effect of the standing ESW is also noticeable in the contourplot of the fluid velocity of electrons, shown in Fig.4. The velocity field has been computed by averaging the velocity over a spatial grid, as in PIC codes, and over time. As a consequence the electric current in the overdense plasma region has a spatial transverse modulation with the same wavelength of ESWs. The simulations of Ref.bastiani  for “grating” targets also show a modulation of the fast electron current correlated with same period of the grating, that may be due to the local enhancement of the longitudinal field by TSWD or to the geometrical “funnel” effect of the surface deformation ruhl . Our test particle results indeed show that due to TSWD a modulated fast electron current can be produced even by a plane–wave pulse on a flat surface.

From figures 2 and 3 it is also found that only near ESW spatial maxima a few electrons are ejected into vacuum, i.e. in the $`x<0`$ region (in Fig.3, for these electrons the oscillation of $`p_x`$ vs. $`x`$ is due to the effect of the $`𝐯\times 𝐁`$ force in vacuum). Their origin is likely to be due to the longitudinal field component directed into vacuum that is associated to the ESWs. In Ref.riconda , the features of electrons accelerated in vacuum by the field of a single ESW are investigated. Note, however, that electrostatic back-holding fields are not self-consistently included in test particle simulations and thus the number of electrons escaping in vacuum is likely to be overestimated. Nevertheless, it is interesting to notice that in the PIC simulations at very high intensities Macchi01 ; MacchiLPB02  “plumes” of electrons extending into the vacuum region are generated near the maxima of the standing ESW. A “plume” structure in the vacuum region is also evident in the velocity field shown in Fig.4.

### Acknowledgments.

This work was partly supported by the Italian Ministery of University and Research (MIUR) through the project “Generation of fast electron and ion beams by superintense laser irradiation”. Discussions with Serena Bastiani-Ceccotti, Caterina Riconda and Francesco Pegoraro are gratefully acknowledged.


text:
# Magnetic versus crystal field linear dichroism in NiO thin films

## Abstract

We have detected strong dichroism in the Ni $`L_{2,3}`$ x-ray absorption spectra of monolayer NiO films. The dichroic signal appears to be very similar to the magnetic linear dichroism observed for thicker antiferromagnetic NiO films. A detailed experimental and theoretical analysis reveals, however, that the dichroism is caused by crystal field effects in the monolayer films, which is a non trivial effect because the high spin Ni $`3d^8`$ ground state is not split by low symmetry crystal fields. We present a practical experimental method for identifying the independent magnetic and crystal field contributions to the linear dichroic signal in spectra of NiO films with arbitrary thicknesses and lattice strains. Our findings are also directly relevant for high spin $`3d^5`$ and $`3d^3`$ systems such as LaFeO3, Fe2O3, VO, LaCrO3, Cr2O3, and Mn4+ manganate thin films.

Magnetic linear dichroism (MLD) in soft-x-ray absorption spectroscopy (XAS) has recently developed into one of the most powerful tools to study the magnetic properties of antiferromagnetic thin films Thole85 ; Sinkovic90 ; Kuiper93 ; Alders95 ; Alders98 . The contrast that one can obtain as a result of differences in the magnitude and orientation of local moments is essential to determine the spin anisotropy and important parameters like the Néel temperature $`(T_N)`$, as well as to map out spatially the different magnetic domains that are present in antiferromagnetic films Stohr98 ; Spanke98 ; Stohr99 ; Scholl00 ; Nolting00 ; Ohldag01a ; Ohldag01b ; Zhu01 ; Hille01 . Such information is extremely valuable for the research and application of magnetic devices that make use of exchange bias.

Much of the modern MLD work has been focussed on NiO and LaFeO3 thin films, and the observed dichroism has been attributed entirely to magnetic effects Stohr98 ; Spanke98 ; Stohr99 ; Scholl00 ; Nolting00 ; Ohldag01a ; Ohldag01b ; Zhu01 ; Hille01 . Other sources that could contribute to linear dichroism, however, such as crystal fields of lower than octahedral symmetry, have been neglected or not considered. Indeed, one would expect that such low symmetry crystal fields are negligible for bulk-like NiO and LaFeO3 films, and, more fundamentally, that such crystal fields will not split the high-spin Ni $`3d^8`$ or Fe $`3d^5`$ ground state. We have illustrated this insensitivity in Fig. 1 for the Ni2+ case, where the energy level diagram in an $`O_h`$ environment is compared to that in a $`D_{4h}`$ point group symmetry SOremark . In contrast, an exchange field will split the Ni2+ ground state into three levels with $`M_S`$=-1,0,1 with an energy separation given by the exchange coupling $`J`$, see Fig. 1. The basis for obtaining strong dichroism in the Ni $`L_{2,3}`$ ($`2p`$$``$$`3d`$) absorption spectra is that dipole selection rules dictate which of the quite different final states can be reach and with what probability for each of the initial states. The isotropic spectrum of each of these three states will be the same, but each state with a different $`|M_S|`$ value will have a different polarization dependence Thole85 ; Sinkovic90 ; Kuiper93 ; Alders95 ; Alders98 . A completely analogous argumentation can be given for the orbitally highly symmetric high spin $`3d^5`$ and $`3d^3`$ cases, e.g. Mn2+, Fe3+, V2+, Cr3+, Mn4+.

In this paper we report on XAS measuments on single monolayer (ML) NiO films which are grown on a Ag(100) substrate and capped by a 10 ML MgO(100) film. We have observed strong linear dichroism in the Ni $`L_{2,3}`$ spectra, very similar to that measured for thicker NiO films. From a detailed theoretical and experimental analysis, however, we discovered that the dichroism can not be attributed to the presence of some form of magnetic order, but entirely to crystal field effects. The analysis provides us also with a practical guide of how to disentangle quantitatively the individual contributions to the linear dichroic signal, i.e. the contribution from magnetic interactions versus that from low symmetry crystal fields. This is important for a reliable determination of, for instance, the spin moment orientation in NiO as well as LaFeO3, Fe2O3, VO, LaCrO3, Cr2O3, and Mn4+ manganate ultra thin films, surfaces and strained films, where the low symmetry crystal field splittings may not be negligible as compared to the exchange field energies.

The polarization dependent XAS measurements were performed at the Dragon beamline of the National Synchrotron Radiation Research Center in Taiwan. The spectra were recorded using the total electron yield method in an XAS chamber with a base pressure of 3x10-10 mbar. The photon energy resolution at the Ni $`L_{2,3}`$ edges ($`h\nu 850880`$ eV) was set at 0.3 eV, and the degree of linear polarization was $`98\%`$. A NiO single crystal is measured simultaneously in a separate chamber upstream of the XAS chamber in order to obtain a relative energy reference with an accuracy of better than 0.02 eV. The 1 ML NiO film on Ag(100) was prepared in Groningen, by using NO2 assisted molecular beam epitaxy. Immediately after the NiO growth, the sample was capped in-situ with an epitaxial 10 ML MgO(100) film. Reflection high energy electron diffraction (RHEED) intensity oscillations recorded during growth of thicker films demonstrated the layer-by-layer growth mode and provided an accurate thickness calibration Altierithesis ; Altieri99 .

Fig. 2 shows the polarization dependent Ni $`L_{2,3}`$ XAS spectra of the 1 ML NiO film, taken at room temperature. The angle between the light polarization vector and the (001) surface normal is given by $`\theta `$ ($`\theta =90^{}`$ means normal light incidence). The general lineshape of the spectra is very similar to that of thicker NiO films and bulk NiO Alders98 .

Fig. 3 presents a close-up of the $`L_2`$ edge, the region most often used to measure the magnitude of the magnetic linear dichroic effect in antiferromagnetic NiO films Stohr98 ; Spanke98 ; Stohr99 ; Ohldag01a ; Ohldag01b ; Zhu01 ; Hille01 . The spectra of the 1 ML NiO film show a very clear polarization dependence. This linear dichroic effect is as strong as that for a 20 ML NiO film grown on MgO(100) (taken from Alders et al. Alders98 ), albeit with an opposite sign, as can be seen from Fig. 3. Relying on the analysis by Alders et al. Alders98  for the antiferromagnetic 20 ML film, one may be tempted to conclude directly that the spin orientation in the 1 ML film is quite different to that of the 20 ML film, i.e. that the spins for the 1 ML would be lying more parallel to the interface while those of the thicker films are pointing more along the interface normal. However, Alders et al. Alders98  have also shown that the magnetic ordering temperature of NiO films decreases strongly if the film is made thinner. In fact, for a 5 ML NiO film on MgO(100), it was found that $`T_N`$ is around or below room temperature, i.e. that no linear dichroism can be observed at room temperature. A simple extrapolation will therefore suggest that 1 ML NiO will not be magnetically ordered at room temperature. This is in fact supported by the 80 K data of the 1 ML NiO on Ag(100) as shown in Figs. 2 and 3: the spectra and the dichroism therein are identical to those at 300 K, indicating that $`T_N`$ must be at least lower than 80 K.

In order to resolve the origin of the linear dichroism in the 1 ML NiO system, we now resort to the Ni $`L_3`$ part of the spectrum. A close-up of this region is given in Fig. 4. We can easily observe that the strong polarization dependence of the spectra is accompanied by an energy shift $`\mathrm{\Delta }E`$ of 0.35 eV in the main peak of the $`L_3`$ white line. This shift seems small compared to the 852 eV photon energy being used, but it is very reproducible and well detectable since the photon energy calibration is done with an accuracy of better than 0.02 eV thanks to the simultaneous measurement of a NiO single crystal reference. We now take this energy shift as an indicator for the presence and strength of local crystal fields with a symmetry lower than $`O_h`$, i.e. crystal fields that do not split the ground state but do alter the energies of the XAS final states, and, via second order processes, also causes spectral weight to be transferred between the various peaks as we will show below.

To understand the Ni $`L_{2,3}`$ spectra quantitatively, we perform calculations for the atomic $`2p^63d^82p^53d^9`$ transitions using the same method as described earlier by Alders et al. Alders98 , but now in a $`D_{4h}`$ point group symmetry. The method uses the full atomic multiplet theory and includes the effects of the solid. It accounts for the intra-atomic $`3d`$-$`3d`$ and $`2p`$-$`3d`$ Coulomb and exchange interactions, the atomic $`2p`$ and $`3d`$ spin-orbit couplings, the O $`2p`$ \- Ni $`3d`$ hybridization with $`pd\sigma `$ = -1.29 eV, and an $`O_h`$ crystal field splitting of $`10Dq`$ = 0.85 eV. The local symmetry for the Ni ion sandwiched between the Ag(100) substrate and the MgO(100) film is in principle $`C_{4v}`$, but for $`d`$ electrons one can ignore the odd part of the crystal field, so that effectively one can use the tetragonal $`D_{4h}`$ point group symmetry. As we will explain below, the $`D_{4h}`$ parameters $`Ds`$ and $`Dt`$ Ballhausen62  are set to 0.12 and 0.00 eV, respectively, and the exchange field (the molecular field acting on the spins) to zero. The calculations have been carried out using the XTLS 8.0 programmTanaka94 .

The right panel of Fig. 4 shows the calculated $`L_3`$ spectrum for the light polarization vector perpendicular and parallel to the $`C_4`$ axis ($`\theta =90^{}`$ and $`\theta =0^{}`$, respectively). One can clearly see that the major experimental features are well reproduced, including the 0.35 eV energy shift between the two polarizations. This shift can be understood in a single electron picture. The ground state has the $`\underset{¯}{3d}_{x^2y^2}\underset{¯}{3d}_{z^2}`$ configuration, where the underline denotes a hole. The final state has a $`2p^5\underset{¯}{3d}_{x^2y^2}`$ or $`2p^5\underset{¯}{3d}_{z^2}`$ configuration. For $`z`$ polarized light the $`3d_{z^2}`$ state can be reached, but the $`3d_{x^2y^2}`$ can not, and the final state will be of the form $`2p^5\underset{¯}{3d}_{x^2y^2}`$. For $`x`$ polarized light the final state will be of the form $`2p^5\underset{¯}{3d}_{z^2y^2}`$$`=`$$`\sqrt{3/4}(2p^5\underset{¯}{3d}_{z^2})`$$`+`$$`\sqrt{1/4}(2p^5\underset{¯}{3d}_{x^2y^2})`$. In a pure ionic picture, the $`2p^5\underset{¯}{3d}_{x^2y^2}`$ state will be 4$`Ds`$+5$`Dt`$ lower in energy than the $`2p^5\underset{¯}{3d}_{z^2}`$ state. In the presence of the O $`2p`$ \- Ni $`3d`$ hybridization, we find that $`Ds`$=0.12 and $`Dt`$=0.00 eV reproduce the observed 0.35 eV shift.

Going back to the $`L_2`$ edge, we can see in Fig. 3 that the calculations can also reproduce very well the observed linear dichroism in the 1 ML NiO spectra. In fact, one now could also see the same 0.35 eV shift at this edge, although it is not as clear as in the $`L_3`$ edge. We would like to stress here that the good agreement has been achieved without the inclusion of an exchange splitting, i.e. the dichroism is solely due to the low symmetry crystal field splitting. It is a final state effect and the change in the ratio between the two peaks of the $`L_2`$ edge as a function of polarization can be understood as follows. In $`O_h`$ symmetry the first peak is due to two final states, one of $`T_2^{}`$ and one of $`E_1^{}`$ symmetry. The second peak is due to a final state of $`T_1^{}`$ symmetry. All three states have the $`2p^5\underset{¯}{3d}_{e_g}`$ configuration. If one reduces the crystal field to $`D_{4h}`$ symmetry the peaks will split. The $`T_2^{}`$ state will split into two states of $`B_2^{}`$ and $`E_1^{}`$ symmetry, the $`E_1^{}`$ into $`A_1^{}`$ and $`B_1^{}`$, and the $`T_1^{}`$ into $`A_2^{}`$ and $`E_1^{}`$. The energy splitting can be measured, but this is much easier done using the $`L_3`$ edge. We note that each of the two peaks in the L2 edge will have a state of $`E_1^{}`$ symmetry, so that these two will mix and transfer spectral weight. This can be seen with isotropic light, but will show up more pronounced as a linear dichroic effect if polarized light is used.

In contrast to the 1 ML NiO case, the good agreement between theory and experiment for the polarization dependent spectra of a 20 ML NiO film Alders98  have been achieved by assuming the presence of an antiferromagnetic order with an exchange field of about 0.16 eV in a pure local $`O_h`$ symmetry. It is surprising and also disturbing that a low symmetry crystal field could induce a spectral weight transfer between the two peaks of the $`L_2`$ white line such that the resulting linear dichroism appears to be very similar as a dichroism of magnetic origin. It is obvious that the ratio between the two peaks can not be taken as a direct measure of the spin orientation or magnitude of the exchange field in NiO films Alders98  if one has not first established what the crystal field contribution could be.

We now can identify two strategies for finding out which part of the linear dichroism is due to low symmetry crystal field effects. The first one is to study the temperature dependence as we have done above. Here we have made use of the fact that those crystal fields do not split the high spin ground state, so that there are no additional states to be occupied with different temperatures other than those already created by the presence of exchange fields. Thus there should not be any temperature dependence in the crystal field dichroism. The linear dichroism due to magnetism, however, is temperature dependent and scales with $`<`$$`M^2`$$`>`$. By going to temperatures high enough such that there is no longer any temperature dependence in the linear dichroism, i.e. when all magnetic ordering has been destroyed, one will find the pure crystal field induced dichroism.

The second strategy to determine the low symmetry crystal field contribution is to measure carefully the energy shift $`\mathrm{\Delta }E`$ in the main peak of the Ni $`L_3`$ white line for $`\theta `$=$`0^{}`$ vs. $`\theta `$=$`90^{}`$. We now calculate the ratio between the two peaks of the Ni $`L_2`$ edge as a function of $`\mathrm{\Delta }E`$, and the results are plotted in Fig. 5 for $`\theta =90^{}`$, $`\theta =0^{}`$, and $`\theta =54.7^{}`$ (isotropic spectrum). Since $`\mathrm{\Delta }E`$ is a function of $`Ds`$ and $`Dt`$ combined, we have carried out the calculations with $`Ds`$ as a running variable for several fixed values of $`Dt`$, and plotted the resulting $`L_2`$ ratios vs. $`\mathrm{\Delta }E`$. We now can use Fig. 5 as a road map to determine how much of the linear dichroism in the Ni $`L_2`$ edge is due to crystal field effects and how much due to magnetism. We can see directly that the 1 ML data lie on curves with the same $`Dt`$, meaning that the measured $`L_2`$ ratios are entirely due to crystal fields. The same can also be said for the 20 ML NiO on MgO at 528 K, which is not surprising since this temperature is above $`T_N`$. However, for the 20 ML NiO at 195 K, one can see that the data points do not lie on one of the $`Dt`$ curves (one may look for larger $`Dt`$ curves, but this results in lineshapes very different from experiment) or let alone on curves with the same $`Dt`$, indicating that one need magnetism to explain the $`L_2`$ ratios. In other words, knowing the $`L_2`$ ratio and $`\mathrm{\Delta }E`$ together allows us to determine the magnitude of the exchange interaction and the orientation of the spin moments. It is best to use the $`\theta =0^{}`$ spectra since here the $`L_2`$ ratio is determined almost by $`\mathrm{\Delta }E`$ alone and is not too sensitive to the individual values of $`Ds`$ and $`Dt`$.

To conclude, we have observed strong linear dichroism in the 1 ML NiO on Ag, very similar to the well known magnetic linear dichroism found for bulk like antiferromagnetic NiO films. The dichroism in the 1 ML, however, can not be attributed to the presence of some form of magnetic order, but entirely to crystal field effects. We provide a detailed analysis and a practical guide of how to disentangle quantitatively the magnetic from the crystal field contributions to the dichroic signal. This is important for a reliable determination of, for instance, the spin moment orientation in NiO as well as LaFeO3, Fe2O3, VO, LaCrO3, Cr2O3, and Mn4+ manganate ultra thin films, surfaces and strained films, where the low symmetry crystal field splittings may not be negligible as compared to the exchange field energies.

We acknowledge the NSRRC staff for providing us with an extremely stable beam. We would like to thank Lucie Hamdan for her skillful technical and organizational assistance in preparing the experiment. The research in Cologne is supported by the Deutsche Forschungsgemeinschaft through SFB 608.


text:
# X-Ray Sources Overdensity Around the 3C 295 Galaxy Cluster

## 1 Introduction

N-body and hydrodynamical simulations show that high redshift clusters of galaxies lie at the nexus of several filaments of galaxies (see e.g. Peacock 1999). Such filaments map out the “cosmic web” of voids and filaments of the large scale structure of the Universe. Thus, rich clusters represent good indicators of regions of sky where several filaments converge. The filaments themselves could be mapped out by Active Galactic Nuclei (AGNs), assuming that AGNs trace galaxies.

Cappi et al. (2001) studied the distribution of the X-ray sources around 3C 295 ($`z=0.46`$). The clusters was observed with the ACIS-S CCD array for a short exposure time ($`18`$ ks). They found a high source surface density in the 0.5-2 keV band which exceeds the ROSAT (Hasinger et al. 1998) and Chandra (Giacconi et al. 2002) logN-logS by a factor of $`2`$, with a significance of $`2\sigma `$ .

In this work an analysis of a deeper ($`92`$ ks) Chandra observation of 3C 295 is performed, to check whether the overdensity of such a field in the $`0.52`$ keV band is real or not, to define the structure of the overdensity, and to extend the above considerations to the $`210`$ keV band.

## 2 Observation and Data reduction

Chandra observed the $`16^{}\times 16^{}`$ field around the 3C 295 cluster with ACIS-I on May 18, 2001, for $`92`$ ks.

The data reduction was carried out using the Chandra Interactive Analysis of Observations software version 2.1.3.

All the data analysis has been performed separately in the $`0.52`$ keV, in the $`27`$ keV and in the $`0.57`$ keV bands.

An identical analysis was performed for the Chandra Deep Field South (CDFS) in the $`0.52`$ keV and $`27`$ keV bands for comparison, and to check our analysis methods.

The source detection was carried out using the ‘PWDetect’ algorithm (Damiani et al. 1997). We identified $`89`$ sources in the $`0.52`$ keV band, $`71`$ sources in the $`27`$ keV band and $`121`$ sources in the $`0.57`$ keV band (see fig. 1).

The counts in the $`0.52`$ keV, $`27`$ keV and $`0.57`$ keV bands were converted in $`0.52`$ keV, $`210`$ keV and $`0.510`$ keV fluxes using conversion factors appropriate for a $`\gamma =1.8`$ power law spectrum with a galactic absorption toward the 3C 295 field of $`N_H=1.33\times 10^{20}`$ cm<sup>-2</sup>. Such values take into account the quantum efficiency degradation of the CCD.

## 3 Analysis

The following analysis has been perfomed:

The sky coverage has been computed for 3C 295 and CDFS fields.

The LogN-LogS in the soft and hard bands for the whole filelds have been produced (figs. 2 and 3).

The LogN-LogS in the soft, hard and whole bands have been produced separately for each ACIS-I chip of the 3C 295 observation (see fig. 4 for the $`0.510`$ keV band).

3C 295 chip-to-chip and CDFS LogN-LogS have been fitted to evaluate slopes and normalizations (fig. 5 for the $`0.510`$ keV band).

A two dimensional Kolmogorov-Smirnov test has been applied to check whether the sources were uniformly distributed or not (tab. 1).

The angular correlation function has been computed in order to estimate at which scales the sources in the two fields were clustered; the errors associated to the function (fig. 6) have been calculated using both poisson and bootstrap statistics (Barrow, Bhavsar & Sonoda 1984).

## 4 Main results

The following results have been achieved:

3C 295 and CDFS LogN-LogS are in very good agreement both in the soft and hard band, and with the CDFS LogN-LogS by Rosati et al. 2002 (figs. 2 and 3).

The 3C 295 LogN-LogS in the soft, hard and broad bands computed separately for each ACIS-I chip show an overdensity of sources in the North-East (NE) chip (fig. 4 for the $`0.510`$ keV band) which reflects the clustering of sources clearly visible in fig. 1.

The discrepancy between the normalization of the LogN-LogS for the NE and SW chip is $`3.2\sigma `$, $`3.3\sigma `$ and $`4.0\sigma `$ in the soft, hard and broad band, respectively; fig. 5 shows the normalization vs. slope plot for the broad band.

The two dimensional Kolmogorov-Smirnov test shows that there is a considerable probability that CDFS sources are uniformly distributed ($`P15\%`$ in the soft and hard bands), while the probability that the 3C 295 sources are uniformly distributed is only a few per cent, and drops below 1 if we consider the $`0.57`$ keV band (tab. 1).

The angular correlation function of the 3C 295 sources features a strong signal on scales of a few arcmins, and also on lower scales in the $`0.57`$ keV band (fig. 6). Moreover, the function is above the value found by Vikhlinin & Forman (1995) for a large sample of ROSAT sources. On the other hand, no signs of a similar behavior is featured by the CDFS.

|  | 3C 295 | CDFS |
| --- | --- | --- |
| $`0.52`$ keV | $`\mathrm{3.09\hspace{0.33em}10}^2`$ | $`0.13`$ |
| $`27`$ keV | $`\mathrm{3.84\hspace{0.33em}10}^2`$ | $`0.17`$ |
| $`0.57`$ keV | $`\mathrm{6.32\hspace{0.33em}10}^3`$ | - |

Tab.1 KS test for CDFS and 3C 295 fields

## 5 Discussion

In this work we studied in great quantitative details the excess of sources clearly visible in the upper left corner of the Chandra observation of the 3C 295 galaxy cluster field (see fig. 1). Since N-body and hydrodynamical simulations show that clusters of galaxies lie at the nexus of several filament, this excess could represent a filament of the large scale structure of the Universe.

Moreover, if the redshift of our sources were the same of 3C 295 ($`z=0.46`$), we could compute the galaxy overdensity of the field assuming a spherical distribution. Under these assumption, we obtain a density contrast of $`70`$ which, despite large uncertainties, is intriguingly close to the expected galaxy overdensity of filaments $`10÷10^2`$, and much smaller than the overdensities of clusters of galaxies ($`10^3÷10^4`$).

Thus, in order to study this candidate filament, further Chandra observations and optical identifications are pursued to obtain the redshift of the sources, to map out the 3C 295 region, and delimiting the filament properties up to scales of $`24`$ arcmin (i.e., $`6`$ Mpc) from the 3C 295 cluster.

More details on the present work can be found in D’Elia et al., A&A, submitted.