text:
# Introduction

## Introduction

The Protein Data Bank (http://www.rcsb.org/pdb)  has about 19,000 protein structures solved by several methods. There are many algorithms that use this information to derive understanding about protein interactions. Our method is based on using these equilibrium configurations and a maximum entropy algorithm to derive information about physical energy functions that could be used to approximate protein interactions.

## Our Method

The following is based largely on an algorithm originally developed by reference . The method is based on the following assumptions. The system is assumed to be in thermodynamic equilibrium. For proteins, this was shown to be a good assumption by Anfinsen . We also assume the energy can be written as a sum of terms which are products of parameters and functions of the configuration. $`E(\mathrm{\Gamma },\stackrel{}{P})=_ip_ih_i(\mathrm{\Gamma })=\stackrel{}{P}\stackrel{}{H}`$ where $`\mathrm{\Gamma }`$ represents the configuration of the system(s). $`\stackrel{}{P}=\{p_i\}`$ represents the set of parameters to be derived.

The probability of a configuration, given parameters, is given by the Boltzmann distribution $`Prob(\mathrm{\Gamma }|\stackrel{}{P})=e^{E(\mathrm{\Gamma },\stackrel{}{P})/kT}/Z=e^{(E(\mathrm{\Gamma },\stackrel{}{P})+F(\stackrel{}{P}))/kT}`$ ,where $`Z(\stackrel{}{P})=_\mathrm{\Gamma }Exp(\beta E(\mathrm{\Gamma },\stackrel{}{P}))`$ and $`F(\stackrel{}{P})=kTln(Z(\stackrel{}{P}))`$.

If we are given the exact equilibrium conformation, $`\mathrm{\Gamma }^{}`$ , the maximum likelihood of parameter values are those values for which the probability, $`Prob(\mathrm{\Gamma }^{}|\stackrel{}{P})`$ is a maximum wrt $`\stackrel{}{P}`$. Maximizing an exponential corresponds to maximizing the argument (ignoring the multiplicative constant $`\beta `$), $`E(\mathrm{\Gamma }^{},\stackrel{}{P})+F(\stackrel{}{P})=Q(\stackrel{}{P})`$. This also corresponds to extremizing the entropy $`TS=EF`$.

Our method is basically the multi-dimensional form of Newton’s method for optimizing functions. Maximizing $`Q(\stackrel{}{P})`$, Newton’s Method is

$$\stackrel{}{P}^{k+1}=\stackrel{}{P}^kD^2(Q(\stackrel{}{P}^k))^1D(Q(\stackrel{}{P}^k))$$
(1)
where $`(D^2)^1`$ represents the inverse Hessian matrix and D represents the gradient. In practice this is modified slightly,

$$\stackrel{}{P}^{k+1}=\stackrel{}{P}^k+ϵ(\mathrm{\Delta }\stackrel{}{P})$$
(2)
where the use of $`ϵ<1`$ corresponds to the ”Damped Newton’s Method”.

Using statistical mechanical definitions

$$\frac{(E(\mathrm{\Gamma }^{}))}{p_i}=h_i^{}\frac{^2(E(\mathrm{\Gamma }^{}))}{p_ip_j}=0$$
(3)
$$\frac{F}{P_i}|_{P^k}=kT(1/Z)\underset{\mathrm{\Gamma }}{}\beta h_ie^{\beta E(\mathrm{\Gamma },P)}=<h_i>$$
(4)
$$\frac{^2F}{p_ip_j}=\beta (<h_i><h_j><h_ih_j>)$$
(5)
Maximizing $`Q=E+F`$ wrt $`\stackrel{}{P}`$ leads to the following.

$$D(Q)_i=h_i^{}+<h_i>$$
(6)
$$D^2(Q)_{i,j}=\beta (<h_i><h_j><h_ih_j>)=\beta Cov(h_i,h_j)$$
(7)
Resulting in the following iterative equation where $`VCM(\stackrel{}{H})`$ is the variance-covariance matrix of $`\stackrel{}{H}`$

$$\mathrm{\Delta }\stackrel{}{P}=kTVCM(\stackrel{}{H})^1(<\stackrel{}{H}>\stackrel{}{H^{}})$$
(8)
The method is easily generalized to a distribution of equilibrium configurations.

$$\mathrm{\Delta }\stackrel{}{P}=kTVCM(\stackrel{}{H})^1(<\stackrel{}{H}><\stackrel{}{H}>_{Prob(\mathrm{\Gamma })})$$
(9)
$`<\mathrm{}>`$ represents a Boltzmann average and $`<\mathrm{}>_{Prob(\mathrm{\Gamma })}`$ represents an average over the given distribution.

Assuming the least prior information, the iteration starts with all parameters set to zero. This would not allow any useful MC evolution at all. The energy would be zero for any Monte Carlo move, thus not preferring any particular move. The energy is modified for the first few iterations with the addition of a clamping term.

$$E_{clamp}=P_{clamp}(rr^{})^2$$
(10)
This term makes the given conformation a minimum of the energy. Since this conformation is an equilibrium conformation, this seems to be a good approximation. Once the parameter values are sufficiently away from zero, this term is set to zero. If a distribution of conformations is given, this distribution can be used as the clamping terms.

### Computation

The basic algorithm is as follows.

* Given $`P^k`$, Monte Carlo simulations and averaging are used to find $`<h_j>`$ and $`<h_ih_j>`$.
* This leads to a matrix equation. $`\mathrm{\Delta }\stackrel{}{P}=kTVCM(\stackrel{}{H})^1(<\stackrel{}{H}>\stackrel{}{H}^{})`$
* Solve this for $`\mathrm{\Delta }\stackrel{}{P}`$
* $`\stackrel{}{P}^{k+1}=\stackrel{}{P}^k+ϵ(\mathrm{\Delta }\stackrel{}{P})`$
* Repeat until $`\stackrel{}{P}`$ converges.

### Computation: Convergence time

Assuming a system of M proteins, each containing N particles with P parameters. The Monte Carlo time required for useful averaging scales at least as $`N`$. The pair-wise energy calculation in each Monte Carlo step scales as $`N^2`$. Although for some models a radius cutoff and neighbor list were used, so this effect in practice is less than indicated. Solving a P by P matrix equation is a $`P^3`$ operation. Computation time for this method scales as

$$\alpha MN^3+\beta P^3$$
(11)
## Simulation

A Monte Carlo evolution was done on several systems. Four energy terms were used.

$$\pm \frac{a}{|r_ir_j|}+\frac{b}{(r_ir_j)^{12}}+c|r_ir_{i+1}|+d(r_ir_{i+1})^2$$
(12)
Non-bonded interactions were included similar to an electrostatic attraction and van der Waals repulsion. The sequence of positive and negative charges was randomly chosen. A covalent-type bond between one particle and the next in the sequence was given a quadratic form. For the simulation described below, the following parameters were used $`a=2,b=4,c=4,d=1`$. The results of the algorithm applied to this system are shown in figure 1 and table 1.

### Protein Model

Several major simplifications were made to allow convergence of real protein parameters in a reasonable time. The united residue approximation is used in the following model. Each of the residues is treated as one particle. This approximation is commonly used . This greatly simplifies the system, but should contain enough complexity to describe the system adequately.

Proteins seem to have a rugged free energy surface. To minimize this effect we use relatively few MC time steps per iteration. This keeps the protein in the local minimum of free energy even if parameters are far from correct.

Residues were placed at the $`C_\alpha `$ location, covalently bond only with next and previous particles on the chain. The energy function used for the covalent bonds was a normal distribution using the mean and variance derived from the data. The noise in the covalent bond parameter seemed to cause large perturbations in the convergence for other parameters. These parameters need to be derived separately.

The first model with consistent convergence used a statistical grouping of residues developed by Cieplak etal . The groupings were derived using a simplification of a statistically derived interaction matrix, the Miyazawa Jernigan (MJ) matrix . This is a simple, consistent grouping which decreases the number of parameters.

* Hydrophobic I (ave hydrophobicity scale value 2.6) (LFI) Leucine, Isoleucine and Phenylalanine
* Hydrophobic II (ave HP scale value 1.8 with large variance) (MVWCY) Methionine, Valine, Tryptophan, Cysteine and Tyrosine
* Polar I (ave HP scale value 1.15) (HA) Histidine, Alanine
* Polar II (ave HP scale value 0.6) (TGPRQSNED) Threonine, Glycine, Proline, Arginine, Glutamine, Serine, Asparagine Glutamic acid and Aspartic acid
* Lysine (ave HP scale value 1.9) (K) Lysine

Only one energy term was used corresponding to van der Waals attraction

$$4ϵ[(\sigma /r)^{12}(\sigma /r)^6]$$
(13)
Essentially this is a contact energy function. $`\sigma `$ was determined by comparing typical volumes and treating the residues as spheres. This gives radii from 2.4 - 3.8 Angstrom. $`\sigma `$ in the equation corresponds to where the core repulsion occurs (about 2\*radius) so a value of 5 was arbitrarily assumed.

In summary, the covalent properties were approximated from the mean and variance of bond lengths. The energy function has one term, a 6-12 combined term. A statistical grouping was used to further reduce parameters. This grouping and energy function model has 15 parameters. Results shown below were derived from proteins ranging in size from 20 to 400 residues. For consistency, only X-ray data and only complete proteins containing no extraneous molecules were used. The training set contained 821 proteins. All protein structures were obtained from the Protein Data Bank . Only 20 MC steps per particle were used.

### Protein Models - Results

Energy function: $`4ϵ[(\sigma /r)^{12}(\sigma /r)^6]`$ with $`\sigma =5`$

The zeros were artificially created, as the algorithm can not handle these parameters going negative. $`(ϵ0)`$ This would cause the MC to diverge. Despite this limitation, convergence was achieved. These results imply the least hydrophobic (most polar) group essentially has no non-bonded interactions. This is very similar to the HP model of polymers where hydrophobic collapse is modeled as HH atraction and other interactions (HP and PP) are ignored. This took 14 days on a Dual PIII 450MHz.

### Protein Models - Model Evaluation

This energy function was applied to several decoy sets. Despite the simplicity of the energy function, results were mixed with several very encouraging successes. All decoy sets were obtained from http://dd.stanford.edu/

### 4 State Reduced Decoy Set

### Local Minima Decoy Set (lmds)

Similar analysis was done for the lmds decoy set from C Kesar and M Levitt, 1999. These results were not very encouraging. The correct proteins were the worst scores in almost all cases and by a large amount. Currently the main weakness is due to the difference in derivation for the covalent bonds. This decoy set was created using a minimization of a backbone torsional energy function, hence was very different from our function. Essentially, our successful decoy set predictions are based on only non-bonded interaction calculations.

### Fisa Casp3 Decoy Set

### Protein Models - Discussion

The energy function used is extremely simple. Despite this, ranking of the native protein for some of the decoys was very encouraging. The poor performers probably need more complex energy functions. Correlation between RMDS and energy was investigated, but no simple relationship was found. For the best performers, the native energy was typically isolated at the lowest energy with most decoys concentrated a distinct difference in energy away.

The algorithm has potential, but several problems must be overcome. Complexity has not been handled very well and may be required for applicability. Inclusion of Coulomb, torsional, angular and backbone potential terms are required for realistic models. All atom and explicit solvent are further steps that can be taken. Terms of smaller magnitude are dominated by effects of terms with larger magnitude and have to be separately derived. Similarly, terms with smaller frequency of occurrence are dominated by more frequent terms. This separate derivation of terms can be organized and iterated consistently.

### Acknowledgements

J Deutsch, Lik Wee Lee, Leif Poorman, Stefan Meyer, TJ Cox, B Allgood


format: text

text:
# Weak approximation for linear systems of quadrics

## Abstract.

We give local conditions at $`\mathrm{}`$ ensuring that the intersection of $`n`$ quadrics in $`^N`$, $`Nn`$, satisfies weak approximation.

The second named author was partially supported by NSF grant DMS-0100537.

Let $`K`$ be a number field, $`V`$ a finite-dimensional vector space over $`K`$, and $`W\mathrm{Sym}^2V^{}`$ a $`K`$-vector space of quadratic forms on $`V`$. We can regard $`W`$ as a linear system of quadrics in $`V`$, and we denote the intersection of all quadrics in the system $`X_W`$. If $`S`$ is a finite set of places of $`K`$, we can ask whether $`X_W`$ satisfies weak approximation with respect to $`S`$, i.e., whether the diagonal map

$$X_W(K)\underset{sS}{}X_W(K_s)$$

has dense image. Assuming that $`X_W(K_s)\mathrm{}`$ for all $`sS`$, this is a refinement of the basic question of whether $`X_W(K)`$ is non-empty.

For $`dimW=1`$, the Hasse-Minkowski theorem implies $`X_W(K)`$ is non-empty if $`dimV5`$ and a generator of $`W`$ is indefinite at every real place of $`K`$. The weak approximation theorem for quadric hypersurfaces \[PR\] §7.1 Cor. 1 then applies. For $`dimW=2`$, J-L. Colliot-Thélène, J-J. Sansuc, and P. Swinnerton-Dyer \[CSS2, 10.4\] gave an affirmative answer to a question of A. Pfister \[Pf\]: if $`dimV>4dimW`$ and $`K`$ is totally imaginary (Pfister originally raised the question in the case $`K=(i)`$), must $`X_W(K)`$ be non-empty? D. Leep \[Le, 2.7\] proved a qualitative version of Pfister’s conjecture: if

$$dimV2(dimW)^2+2dimW+1,$$

and $`K`$ is totally imaginary, then $`X_W(K)`$ is non-empty. Using the circle method, W. Schmidt \[Sc\] proved an analogous result over $``$, under the hypothesis that $`X_W()`$ has a non-singular point.

Weak approximation is more delicate. A simple relation between the number of variables and the number of equations is no longer possible, even in the case of totally imaginary fields. For instance, if $`\alpha K`$ does not have a square root in $`K`$, then the quadric hypersurface given by $`x_1^2\alpha x_2^2=0`$ fails to satisfy weak approximation no matter how large its dimension. The obstruction can be understood in terms of the rank of the defining quadratic form being too low. For $`W`$ generated by a single indefinite quadratic form of rank $`5`$, weak approximation of $`X_W`$ is well known (see, e.g., \[PR\] S7.1 Cor. 1.) Colliot-Thélène, Sansuc, and Swinnerton-Dyer have a weak approximation theorem (\[CSS1\], \[CSS2\]) for $`dimW=2`$, $`dimV9`$, assuming every form is indefinite and has rank $`5`$. See also \[SS\] and \[Sko\] for good bounds when $`dimW3`$.

The main result of this paper is a weak approximation theorem for $`X_W`$ for arbitrary $`dimW`$ and arbitrary number fields $`K`$. The hypotheses make no explicit demands concerning the singular locus of $`X_W`$; instead, we impose conditions on the rank and signature of non-zero forms in $`W`$ with respect to the archimedean places of $`K`$.

There seems to be a general philosophy that a complete intersection $`X`$ of $`n`$ hypersurfaces of degree $`d`$ should have many rational points as long as there are no local obstructions (especially no obstruction at real places) and $`dimX`$ is large compared to $`n`$, $`d`$, and the dimension of the singular locus of $`X`$. This paper exhibits one concrete instance of this philosophy for $`d=2`$. For general $`d`$, C. Skinner \[Ski\] has proven weak approximation for smooth complete intersections as long as a certain auxiliary variety is not too badly singular. In particular, he has proven weak approximation for non-singular cubic hypersurfaces of sufficiently high dimension. It is not unlikely that his work can provide an alternative route to our result.

We would like to express our appreciation of the help we received from Colliot-Thélène, whose pointed questions and suggestions led to considerable improvements in the exposition of this paper. The original admissibility bound in our theorem was not stated (and was in fact worse than exponential in $`dimW`$). Both Colliot-Thélène and the referee strongly urged us to make it explicit, and thereby encouraged us to improve it to its current form (which is quadratic). We are grateful to both of them.

We begin with a definition intended at the same time to control the singular locus and to guarantee that there is no real obstruction:

###### Definition 1.

If $`V`$ is a vector space over $`k=`$ or $`k=`$, $`W\mathrm{Sym}^2V^{}`$ a vector space of quadratic forms on $`V`$, and $`m`$ is a positive integer, we say $`W`$ is *$`m`$-admissible* if for every non-zero $`wW`$, there are at least $`m`$ hyperbolics contained in $`w`$. When $`k=`$, this just means that the rank of $`w`$ is at least $`2m`$; over $``$ it means $`V`$ contains $`m`$-dimensional subspaces $`V_w^+`$ and $`V_w^{}`$ on which $`w`$ is positive-definite and negative-definite respectively. If $`V`$ is a vector space over a number field $`K`$, $`W\mathrm{Sym}^2V^{}`$, and $`m`$ is a positive integer, we say $`W`$ is *$`m`$-admissible* if $`W1`$ is an $`m`$-admissible space of quadratic forms on $`V_KK_s`$ for every archimedean completion $`k=K_s`$ of $`K`$.

We can now state the main theorem.

###### Theorem 2.

Let $`f(n)=2n^2+2n1`$. For every positive integer $`n`$, number field $`K`$, $`K`$-vector space $`V`$, and $`n`$-dimensional vector space of quadratic forms $`W\mathrm{Sym}^2V^{}`$ which is $`f(n)`$-admissible, the variety $`X_W`$ always satisfies weak approximation over $`K`$. That is, for every finite set $`S`$ of places of $`K`$, $`X_W(K)`$ is dense in $`_{vS}X_W(K_s)`$.

The growth rate of $`f(n)`$ is quadratic like those of \[Le\] and \[Sc\], but the actual minimum $`dimV`$ to which the theorem can be applied is larger than the bounds of Leep and Schmidt. We have $`f(2)=11`$, $`f(3)=23`$, so our result certainly does not improve existing bounds in low dimensions. The virtues of our theorem are that it gives full weak approximation rather than mere existence of solutions, and it does so for an arbitrarily large system of quadratic forms and without explicit hypotheses on the geometry of $`X_W`$. Moreover, the proof is elementary. To find points in $`X_W(K)`$, we choose a quadratic form $`w_0W`$ and a complementary subspace $`W^0W`$. We then proceed to find five vectors $`v_1,\mathrm{},v_5`$ which are mutually orthogonal with respect to $`W`$ and isotropic with respect to $`W^0`$ but not with respect to $`w_0`$. Then we apply weak approximation to the locus $`w_0=0`$ on the vector space spanned by the $`v_i`$.

The following notation will be useful. Given vector spaces $`V`$, $`W\mathrm{Sym}^2V^{}`$, and $`V_0V`$, we have a natural surjective linear transformation $`\mathrm{Sym}^2V^{}\mathrm{Sym}^2V_0^{}`$. The image of $`W`$ under this transformation will be denoted $`\mathrm{Res}_{V_0}W`$. Note that if $`\mathrm{ker}V^{}V_0^{}`$ is spanned by linear forms $`L_1,\mathrm{},L_r`$ on $`V`$, a quadratic form on $`V`$ lies in $`\mathrm{ker}\mathrm{Sym}^2V^{}\mathrm{Sym}^2V_0^{}`$ if and only if it is of the form $`_{i=1}^rL_iM_i`$ for $`M_iV^{}`$. In particular, the kernel of restriction consists only of forms of rank $`2r`$.

###### Lemma 3.

Let $`V`$ be a finite-dimensional vector space over $`k=`$ or $`k=`$ and $`W\mathrm{Sym}^2V^{}`$ be a vector space of quadratic forms. Let $`V_0V`$ be a subspace of codimension $`d`$. If $`W`$ is $`m`$-admissible for some integer $`m>d`$, then $`\mathrm{Res}_{V_0}W`$ is $`md`$-admissible.

###### Proof.

Let $`\overline{w}\mathrm{Res}_{V_0}W`$ be non-zero, and let $`w`$ be an element of $`W`$ mapping to $`\overline{w}`$. If $`k=`$, there there exists a $`2m`$-dimensional subspace $`V_wV`$ in which $`w`$ is non-degenerate. Then $`V_wV_0V_0`$ is a subspace of dimension $`2md>2(md)`$ on which $`\overline{w}`$ is non-degenerate. For $`k=`$, let $`V_w^+`$ (resp. $`V_w^{}`$) denote an $`m`$-dimensional subspace of $`V`$ on which $`w`$ is positive-definite (resp. negative-definite). Then $`\overline{w}`$ is positive-definite (resp. negative-definite) on $`V_wV_0`$, which has dimension $`md`$. ∎

Given a real or complex vector space $`V`$ and a vector space $`W\mathrm{Sym}^2V^{}`$ of quadratic forms, we define the (non-linear) *evaluation map* $`E:VW^{}`$ by the identity

$$E(v)(w)=w(v).$$

###### Proposition 4.

Let $`g(n)=n^2n+1`$. Let $`n`$ be a positive integer, $`V`$ a finite-dimensional vector space over $`k=`$ or $`k=`$, and $`W\mathrm{Sym}^2V^{}`$ an $`n`$-dimensional space of quadratic forms. If $`W`$ is $`g(n)`$-admissible, then the evaluation map $`E:VW^{}`$ is surjective.

###### Proof.

The idea of the proof is that for fixed $`v_1,v_2V`$, the map

$$(a,b)E(av_1,bv_2)$$

is of the form

$$(a,b)a^2w_1^{}+b^2w_2^{}+abw_{1,2}^{}$$

for suitable vectors $`w_1^{}=E(v_1)`$, $`w_2^{}=E(v_2)`$, and $`w_{1,2}^{}W^{}`$. In general, the image is therefore a cone in a 3-dimensional subspace of $`W^{}`$. If, however, it happens that $`w_{1,2}^{}`$ is a linear combination of $`w_1^{}`$ and $`w_2^{}`$, then this formula can be used to show that a wide range of linear combinations of $`w_1^{}`$ and $`w_2^{}`$ belong to $`E(V^{})`$. For general $`v_1`$ and $`v_2`$, there is no reason to believe this condition will hold, but we can use an inductive argument to prove that there exists $`v_3V`$ such that the analogous conditions hold for $`v_1`$ and $`v_3`$ and for $`v_2`$ and $`v_3`$, and moreover the vector $`w_3^{}=E(v_3^{})`$ belongs to the space of $`w_1^{}`$ and $`w_2^{}`$. The details are somewhat simpler in the case $`k=`$, so we begin there.

Suppose $`k=`$. It suffices to prove that $`E(V)`$ is a vector space, since if $`w(E(V))`$ = 0 for some vector $`wW`$, then $`w(v)`$ = 0 for all $`vV`$, and this implies $`w=0`$. Let $`v_1,v_2V`$, and let $`w_i^{}=E(v_i)`$. We would like to show that every linear combination of $`w_1^{}`$ and $`w_2^{}`$ lies in $`E(V)`$. As $`E(av)=a^2E(v)`$, this is obvious unless $`w_1^{}`$ and $`w_2^{}`$ span a two-dimensional space.

Let $`W_1=(W^{}/w_1^{})^{}`$ denote the hyperplane in $`W`$ determined by $`w_1^{}`$. Let

$$V_{1,2}=\{vVw(v+v_1)=w(v)+w(v_1),w(v+v_2)=w(v)+w(v_2)wW_1\}$$

be the subspace of $`V`$ consisting of vectors orthogonal to $`v_1`$ and $`v_2`$ with respect to every quadratic form in $`W_1`$ (i.e., with respect to every symmetric bilinear form associated to a quadratic form in $`W_1`$). It is enough to check this orthogonality on a basis for $`W_1`$, so $`V_{1,2}`$ is of codimension $`2n2`$ in $`V`$. By Lemma 3, $`\mathrm{Res}_{V_{1,2}}W_1`$ is $`g(n)(2n2)=g(n1)`$-admissible. Note that as $`g(n)2n2`$, the restriction map of $`W_1`$ to $`V_{1,2}`$ is an isomorphism. By the induction hypothesis, there exists $`v_3V_{1,2}`$ whose image in

$$\mathrm{Res}_{V_{1,2}}W_1^{}=W_1^{}=W^{}/w_1^{}$$

equals the class of $`w_2^{}`$ in $`W^{}/w_1^{}`$. Let $`w_3^{}=E(v_3)`$. Thus $`\mathrm{Span}(w_1^{},w_3^{})=\mathrm{Span}(w_1^{},w_2^{})`$ and $`v_1`$ and $`v_3`$ are orthogonal with respect to every form in $`W_1`$. This last condition implies that there exists $`c`$ such that

$$w(av_1+bv_3)=a^2w(v_1)+b^2w(v_3)+abcw_1^{}(w)=((a^2+abc)w_1^{}+b^2w_2^{})(w).$$

For any $`c`$, the function $`(a,b)(a^2+abc,b^2)`$ is a surjective self-map of $`^2`$, so indeed $`E`$ is surjective.

Now suppose $`k=`$. The formula $`E(av)=a^2E(v)`$ implies $`E(V)`$ is a cone. If $`E(V)`$ is contained in a closed half-space $`H_w=\{w^{}w(w^{})0\}`$ then $`w`$ is semidefinite, contrary to hypothesis. If $`E(V)`$ is convex, then by \[Ro, 11.7.3\], it follows that $`E(V)=W^{}`$. To prove convexity, it suffices to prove that the line segment connecting any $`w_1^{}=E(v_1)`$ and $`w_2^{}=E(v_2)`$ lies in $`E(V)`$. This is clear if $`w_1^{}`$ and $`w_2^{}`$ are linearly dependent, so we may assume that $`\mathrm{Span}(w_1^{},w_2^{})`$ is a plane. As $`E(V)`$ is a cone, it suffices to prove that there is a path connecting $`w_1^{}`$ and $`w_2^{}`$ which lies entirely in a closed half-plane within that plane.

We construct $`V_{1,2}`$, $`v_3`$, and $`w_3^{}`$ as in the complex case. We define $`c,d`$ such that

$$w(av_1+bv_3)=a^2w(v_1)+b^2w(v_3)+abcw_1^{}(w)$$

$$w(av_2+bv_3)=a^2w(v_2)+b^2w(v_3)+abdw_1^{}(w)$$

Replacing $`v_1`$ and $`v_2`$ by $`v_1`$ and $`v_2`$ respectively if necessary, we may assume $`c,d0`$. By construction, $`w_1^{}`$, $`w_2^{}`$, and $`w_3^{}`$ lie in a common half-plane. It follows that

$$E(av_1+bv_3)=a^2w_1^{}+b^2w_3^{}+abcw_1^{}$$

$$E(av_2+bv_3)=a^2w_2^{}+b^2w_3^{}+abdw_1^{}$$

all lie in this half-plane. Specializing to $`(a,b)=(t,1t)`$ for $`t[0,1]`$, we conclude that there exist paths from $`v_3`$ to $`v_1`$ and from $`v_3`$ to $`v_2`$ contained in a common half plane and therefore that there exists a path from $`v_1`$ to $`v_2`$ in that half-plane. The proposition follows.

We can now prove the main theorem.

Proof. We use induction on $`n`$. The case $`n=1`$ is well-known, so we may assume $`n2`$. We fix a non-zero $`w_0W`$ and a codimension-1 subspace $`W^0W`$ such that $`w_0W^0`$. Let $`S`$ be a finite set of places of $`K`$ and $`S_{\mathrm{}}`$ be the set of archimedean places of $`K`$. For each $`sS`$, we fix a non-empty open set $`U_s`$ in $`(V_KK_s)`$. We apply the induction hypothesis to $`W^0`$ and $`_{sS}U_sX_{W^0}(K_s)`$ to obtain a vector $`v_1V`$ such that $`v_1`$ is isotropic for all quadratic forms in $`W^0`$ and for all $`sS`$, the line $`K_sv_1`$ defines a point in $`(V_KK_s)`$ which lies in $`U_s`$. If $`w_0(v_1)=0`$, we are done, so we assume this quantity is not zero. Next we construct vectors $`v_2,\mathrm{},v_5V`$ iteratively with the following properties:

1. For $`i>j1`$, the vectors $`v_i`$ and $`v_j`$ are mutually orthogonal with respect to every form $`wW`$.
2. For $`i2`$, $`v_i`$ is $`W^0`$-isotropic
3. For $`i2`$, $`v_i`$ is not $`w_0`$-isotropic.
4. For every real embedding $`s`$ and every $`i2`$, we have $`\mathrm{sgn}s(w_0(v_i))=\mathrm{sgn}s(w_0(v_1)).`$

To achieve these properties, define for $`2i5`$

$$V_i=\{vVw(v+v_j)=w(v)+w(v_j)j<i\}.$$

Clearly $`V_i`$ is of codimension $`(i1)n4n`$ in $`V`$. By Lemma 3, $`\mathrm{Res}_{V_i}W`$ is $`f(n)4n`$-admissible, and as $`f(n)4n=f(n1)n^2n+1`$, Proposition 4 applies. As the rank of every non-zero $`wW`$ is $`2g(n)>4n`$, the restriction map $`W\mathrm{Res}_{V_i}W`$ is an isomorphism. This means that for every archimedean completion $`K_s`$ of $`K`$, there exists $`v_{i,s}V_i_KK_s`$ which is isotropic for every quadratic form in $`W^0`$ and satisfies $`w_0(v_{i,s})0`$. If $`K_s=`$, we can additionally choose the sign of $`w_0(v_{i,s})`$ to be $`\mathrm{sgn}s(w_0(v_1))`$. For each archimedean place $`s`$ of $`K`$, we define a non-empty neighborhood $`U_{i,s}`$ of $`K_sv_{i,s}`$ in $`X_{\mathrm{Res}_{V_i}W^0}(K_s)`$ given by

$$U_{i,s}=\{\begin{array}{cc}\{K_svvV_i,w_0(v)0\}\hfill & \text{if }K_s=\hfill \\ \{K_svvV_i,s(w_0(v)w_0(v_1))<0\}\hfill & \text{if }K_s=.\hfill \end{array}$$

We now apply the induction hypothesis to $`\mathrm{Res}_{V_i}W^0`$ and $`_{sS_{\mathrm{}}}U_{i,s}`$ to find the vector $`v_i`$. Condition (1) is automatic since $`v_iV_i`$. Condition (2) holds because $`Kv_i`$ lies in $`X_{\mathrm{Res}_{V_i}W^0}(K)`$. Conditions (3) and (4) are true because $`Kv_i`$ lies in $`_{sS_{\mathrm{}}}U_{i,s}`$.

Now consider the restriction of $`w_0`$ to $`V_0:=\mathrm{Span}(v_1,v_2,v_3,v_4,v_5)`$. We know that $`w_0`$ is diagonal in terms of the basis $`v_i`$, of full rank, and indefinite at all real places of $`K`$ by conditions (1), (3), and (4) respectively. Moreover, the intersection

$$PV_0\underset{sS}{}U_s$$

is open and non-empty (because it contains the line $`Kv_1`$). By the weak approximation theorem for indefinite quadratic forms of rank 5, there exists $`vV_0`$ which is isotropic for $`w_0`$ and such that $`Kv_{sS}U_s`$. Conditions (1) and (2) above now imply that $`v`$ is also isotropic for $`W^0`$ and therefore for all of $`W`$. The theorem follows.

We conclude with several questions. Does the theorem remain true if we replace our notion of $`m`$-admissibility by the condition that every non-zero $`wW`$ is indefinite and non-degenerate on an $`m`$-dimensional subspace of $`V`$? Do we get weak approximation if we assume only that $`dimV0`$ and every non-zero $`wW`$ is indefinite and of rank $`5`$?

Our proof of Proposition 4 depends on a trick that shows that under suitable admissibility conditions, the real locus of a system of quadratic equations (in real projective space) is connected. Can this trick be used to analyze higher connectivity? For instance, if $`V`$ is a real vector space, $`k`$ is a positive integer, and $`W`$ a vector space of quadratic forms which is $`m`$-admissible for some $`m`$ depending on $`k`$ and $`dimW`$, can we conclude that $`X_W()`$ has a $`k`$-connected double cover?


text:
# Orbital moment in CoO and in NiO

## Abstract

The total, orbital and spin moment of the Co2+ ion in CoO has been calculated within the quasi-atomic approach with taking into account strong correlations, crystal-field interactions and the intra-atomic spin-orbit coupling. The orbital moment of 1.39 $`\mu _B`$ amounts at 0 K, in the magnetically-ordered state, to more than 34% of the total moment (4.01 $`\mu _B`$) and yields the L/S ratio of 1.05, close to the experimental value. The same calculations yield for NiO the orbital and total moment of 0.46 $`\mu _B`$ and 2.45 $`\mu _B`$, respectively.

PACS No: 71.70.E; 75.10.D

Keywords: 3d magnetism, crystal field, spin-orbit coupling, orbital moment, CoO, NiO

1. Introduction

CoO and NiO attract a large attention of the magnetic community by more than 50 years. Despite of its simplicity (two atoms, $`NaCl`$ structure, well-defined antiferromagnetism (AF) with TN of 290 K and 525 K, respectively) and enormous theoretical and experimental works the consistent description of its properties, reconciling its insulating state with the unfilled 3$`d`$ band is still not reached .

The aim of this paper is to report the calculations of the magnetic moment of CoO and NiO, its direction and the value as well as the spin and orbital parts. In our approach, called Quantum Atomistic Solid-State theory (QUASST), we attribute the moment of CoO/NiO to the Co2+/Ni2+ ions. We have calculated the moment of the Co2+/Ni2+ ion in the CoO6/NiO6 octahedral complex, its spin and orbital parts. The orbital moment at 0 K as large as 1.38 $`\mu _B`$ and 0.46 $`\mu _B`$ has been found in CoO and NiO. The approach used can be called the quasi-atomic approach  as the starting point for the description of a solid is the consideration of the atomic-like low-energy electronic structure of the constituting atoms/ions, in the present case of the Co2+/Ni2+ ions.

2. Theoretical outline

We have treated the 7 outer electrons of the Co2+ ion in Mott insulator Co2+O2- as forming the strongly-correlated electron system 3$`d^7`$. The correlations among electrons in the unfilled 3$`d`$ shell are taken by two Hund’s rules, that yield the ground-term quantum numbers $`S`$=3/2 and $`L`$=3 related to the ground term $`{}_{}{}^{4}F`$ . Such the localized highly-correlated electron system interacts in a solid with the charge and spin surroundings. The charge surroundings has predominantly the octahedral symmetry owing to the $`NaCl`$-type of structure of CoO and NiO. Distortions are important for the detailed formation of the AF structure and influence the spin and orbital moments but the most essential physical interaction is the intra-atomic spin-orbit coupling. Our Hamiltonian for CoO consists of two terms: the single-ion-like term $`H_d`$ of the 3$`d^7`$ system and the $`d`$-$`d`$ intersite spin-dependent term $`H_{dd}`$, that is indispensable for the formation of the magnetic state. The calculations follow those, that we have performed for the description of FeBr2 and NiO . For the quasi-atomic single-ion-like Hamiltonian of the 3$`d^7`$ system we take into account the crystal-field interactions of the octahedral symmetry (the octahedral CEF parameter $`B_4`$=-30 K), the spin-orbit interactions with the spin-orbit coupling $`\lambda _{so}`$=-260 K and the tetragonal distortion approximated by the term $`B_2^0`$= -5 K. The single-ion states under the octahedral crystal field and the spin-orbit coupling have been calculated by consideration of the Hamiltonian:

$$H_d=B_4(O_4^0+5O_4^4)+\lambda _{so}LS$$
(1)
The calculated single-ion states under the octahedral crystal field and the spin-orbit coupling (the CoO6 complex) are presented in Fig. 1(c). The intersite spin-dependent interactions cause the (antiferro-)magnetic ordering. They have been considered in the mean-field approximation with the molecular-field coefficient n acting between magnetic moments $`m_d`$=($`L`$+2$`S`$) $`\mu _B`$. Having the CEF interactions quantified, here $`B_4`$ mostly from optical studies, the value of $`n`$ in the Hamiltonian

$$H_{dd}=n\left(m_dm_d+\frac{1}{2}m_d^2\right)$$
(2)
has been adjusted in order to reproduce the experimentally-observed Neel temperature. The fitted value of $`n`$, for TN of 290 K, has been found to be -2.43 meV/ $`\mu _B^2`$ (=42 T/$`\mu _B`$). It means that the Co-ion moment in the magnetically-ordered state experiences the molecular field of 169 T (at 0 K).

As the result of calculations we get the low-energy electronic structure, i.e. the energy of (2$`L`$+1)$``$(2$`S`$+1) states and their eigenfunctions. Having the energy of states we can calculate the free energy. It enables calculations of thermodynamical properties. The ground-state eigenfunction enables calculations of the total, orbital and spin moments.

In NiO the 8 outer electrons of the Ni2+ ion, forming the strongly-correlated electron system 3$`d^8`$, are described by $`S`$=1 and $`L`$=3 related to the ground term $`{}_{}{}^{3}F`$ . For the calculations of the quasi-atomic single-ion-like Hamiltonian of the 3$`d^8`$ system we take into account the crystal-field interactions of the octahedral symmetry and the spin-orbit coupling (the octahedral CEF parameter $`B_4`$=+2 meV and the spin-orbit coupling constant $`\lambda _{so}`$=-41 meV , p. 399). The preliminary results for NiO have been reported in .

3. Results and discussion

The calculated value of the total magnetic moment at 0 K in CoO and NiO in the magnetically-ordered state amounts to 4.01 and 2.45 $`\mu _B`$, respectively. In CoO it is built up from the spin moment $`m_s`$ of 2.62 $`\mu _B`$ ($`S_z`$=1.31) and the orbital moment $`m_o`$ of 1.39 $`\mu _B`$. In NiO it is built up from the spin moment $`m_s`$ of 1.99 $`\mu _B`$ and the orbital moment $`m_o`$ of 0.46 $`\mu _B`$.

The formation of the magnetic state and of the magnetic moment is associated to the polarization of two conjugate states as is shown in Fig. 2. The spin-like gap at 0 K amounts to 40 meV.

The orbital moment in CoO is quite substantial being more than 34% of the total moment. Our theoretical outcome, revealing the substantial orbital moment is in nice agreement with the very recent experimental result of 3.98$`\pm `$0.06 $`\mu _B`$ for the Co moment . The magnetic x-ray experiment has revealed the $`L`$/$`S`$ ratio of 0.95 \- the calculated by us values lead to $`L`$/$`S`$ ratio, in fact their $`z`$ components, of 1.05. Also our theoretical result about the substantial orbital moment in NiO, is in nice agreement with the very recent experimental result of 2.2$`\pm `$0.3 $`\mu _B`$ for the Ni moment at 300 K .

The effect of the small tetragonal distortion one can note comparing values of the moment in the paramagnetic state shown in Fig. 2 and the moment of the ground state in Fig. 1.

The calculated temperature dependence of the total moment, together with the orbital and spin moments, is shown in Fig. 3. These moments disappear above TN \- in the paramagnetic region the derived moment is the effective moment, that bears the information about the operators $`J^2`$ or $`S^2`$.

We would like to point out that the evaluation of the orbital moment is possible provided the spin-orbit coupling is taken into account. It confirms the importance of the spin-orbit coupling for the description of the 3$`d`$-ion compounds. The magnetically-ordered state of CoO has lower energy than the paramagnetic one by 380 J/mol (= 3.9 meV/ion) at 0 K. Of course, the energies of magnetic and paramagnetic states become equal at $`T_N`$.

Finally, we would like to point out that our approach should not be considered as the treatment of an isolated ion only - we consider the Co2+ ion in the oxygen octahedron. The physical relevance of our discussion to CoO is obvious - the $`NaCl`$ structure is built up from the edge sharing Co2+ octahedra and the same electronic structure occurs at each Co site. The sign of the octupolar CEF parameter $`B_4`$ is directly related to the oxygen anion surroundings and can be calculated from the first principles. Its change of sign from the Co2+ ion to Ni2+ ion is related to the sign reversal of the the octupolar moment of the charge cloud of the 3$`d^7`$ and 3$`d^8`$ configurations, seen in Stevens factor $`\beta `$. The value of $`B_4`$ of -30 K in CoO is fully consistent with a value of +260 K found in LaCoO3  taking into account the change of the valency and the average Co-O distance from 192 pm to 213 pm in CoO.

4. Conclusions

The orbital and spin moment of the Co2+ ion in CoO has been calculated within the quasi-atomic approach taking into account strong correlations, crystal-field interactions and the intra-atomic spin-orbit interactions by the same approach as we have earlier used for NiO , for FeBr2  and for LaCoO3 . The orbital moment of 1.39 $`\mu _B`$ amounts at 0 K, in the magnetically-ordered state, to more than 34% of the total moment of 4.01 $`\mu _B`$. This value is much larger than 3 $`\mu _B`$ expected for 3 spin-holes in case of seven 3$`d`$ electrons in the Co2+ ion. The presented approach explains in a very natural way the insulating state of NiO - in fact, CoO and NiO are the best insulators. Our atomic-like approach provides the discrete energy states for 3$`d`$ electrons in CoO and calculates, apart of the spin moment, the orbital moment. Our studies indicate that it is the highest time in solid-state physics to ”unquench” the orbital moment in the description of 3$`d`$-atom containing compounds, the more that it becomes visible in recent advanced experiments. The consistent description of the low-energy electronic structure and the ordered magnetic moment of four different compounds CoO, NiO, FeBr2 and LaCoO3 within the same atomistic approach we take as very strong argument for high physical adequacy of our approach pointing out the existence of discrete atomic-like states in 3$`d`$-atom containing compounds.


text:
# The signature of substructure on gravitational lensing in the ΛCDM cosmological model

## 1 Introduction

Whereas the current Cold Dark Matter (CDM) paradigm for structure formation is widely accepted, two major problems for CDM remain. While simulations predict cuspy dark matter halos (e.g. Moore 1994), observed rotation curves of low surface brightness galaxies indicate that their dark matter halos have more shallow cores (Kravtsov et al. 1998; Swaters et al. 2000; van den Bosch & Swaters 2001; de Blok & Bosma 2002). The other is the apparent over-prediction of the small-scale power in CDM simulations. As was shown by Moore et al. (1999) and Klypin et al. (1999), the number of satellite halos seen in N-body simulations appears to far exceed the number of dwarf galaxies observed around the Milky Way. Particular discrepancies have been found for satellite masses $`10^9M_{}`$.

Gravitational lensing is at present the only tool to investigate CDM substructure in galaxies outside the local group. As first noted by Mao & Schneider (1998), mass-substructure other than stars on scales less than the image separation can substantially affect the observed flux ratios in strong gravitational lens systems. Chiba (2002), Dalal & Kochanek (2002), Mao & Schneider (1998), Metcalf & Madau (2001), Metcalf & Zhao (2002), Metcalf et al. (2003), Keeton (2001), and Bradač et al. (2002) have all argued that substructure can provide the explanation for the flux anomalies in various systems. Dalal & Kochanek (2002) further conclude that the amount of substructure needed to explain the flux ratios of quadruply-imaged systems broadly agrees with the CDM predictions. At least for some systems the flux mismatches are probably not just an artifact of oversimplified macromodels of the main lens galaxy (see e.g. Evans & Witt 2003; Metcalf & Zhao 2002). As discussed by Keeton (2003) and Chen et al. (2003), fluxes can be further affected by clumps of matter at a redshift different from that of the lens, along the line of sight between the observer and the source; however, this effect is not dominant. It is also possible that the small scale structure does not consist of compact CDM clumps, also tidal streams or offset disc components can affect the flux ratios (see Möller et al. 2003; Quadri et al. 2003).

Keeton (2001) and Gaudi & Petters (2002) recently focused on the magnification relations that should be satisfied by particular four-image geometries (so called “fold” and “cusp” configurations). These relations are model-independent predictions for the magnifications of highly magnified images (Blandford & Narayan 1986; Blandford 1990; Schneider & Weiss 1992; Mao 1992). Strictly speaking, however, they hold only for ideal “fold” or “cusp” configurations and it is therefore in some cases hard to disentangle the effects of the source being further away from the cusp from the effect of substructure, purely by employing these relations.

The influence of substructure can not only be seen on image flux ratios, but also in the structure of multiple-imaged jets. The lens system B1152+199 consists of a doubly-imaged jet, one of which appears bent, whereas the other is not (Metcalf 2002). Alternative explanation is that an intrinsic bend in the jet is simply magnified in one image, and produces only a small effect in the other.

Microlensing can change the flux ratios not only in the optical (e.g. Woźniak et al. 2000), but also at radio wavelengths (Koopmans & de Bruyn 2000). Flux ratio anomalies can also be introduced by propagation effects in the interstellar medium (ISM) in the lens galaxy, by galactic scintillation, and scatter broadening (Koopmans et al. 2003). Fortunately, these effects are frequency dependent and one can distinguish them using multi-frequency observations. In addition, these electromagnetic phenomena are similar for images of different parities.

For substructure, on the other hand, Schechter & Wambsganss (2002) found that magnification perturbations should show a dependence on image parity. Microlensing simulations showed that the probability distributions for magnifications of individual images are not symmetric around the unperturbed magnification. The distribution depends on image parity and becomes highly skewed. The probability for the brightest saddle point image to be demagnified is increased.111Images form at extrema of the arrival time surface (i.e. Fermat’s principle). At saddle points, images have negative and at minima/maxima, they have positive parity. In quadruply-imaged system one observes two saddle point images, and two minima. The fifth image is a maximum, too faint to be observed. Observed lens systems also seem to show this image parity dependence (Kochanek & Dalal 2003), and this indicates that the flux ratio anomalies arise mainly from gravitational lensing, rather than propagation effects.

All these effects on flux ratios have placed some doubt as to whether the existence of substructure can be rigorously tested with strong lensing and what is expected signal. Several groups have tested the effects of substructure in strong lensing systems using a semi-analytic prescription for substructure (Metcalf & Madau 2001; Keeton et al. 2003; Dalal & Kochanek 2002; Kochanek & Dalal 2003). Recently Mao et al. (2004) showed using high-resolution numerical simulations that the fraction of surface mass density in substructures is lower than required by lensing; however both predictions are still uncertain.

In this paper we use different projections of two different galaxies obtained in N-body+gasdynamics simulations. Whereas a semi-analytic prescription overcomes the problem of shot-noise, and the problem of modelling becomes simpler (one has an analytic model for the underlying macro-distribution), by using the direct output of an N-body simulated galaxy, one does not make any simplifying assumption about the mass profiles of the macro model, or the substructure. Down to the resolution scales of the simulation we therefore believe we have a better comparison with a realistic galaxy.

Whereas higher-resolution DM-only simulations are available, the absence of baryons significantly affects lens properties. Those type of simulations are of limited use for the purpose of testing CDM substructure effect on strong lens systems. More precise, strong gravitational lensing is probing the galaxy potential on the inner $`510\text{ kpc}`$ (the typical size of the Einstein radius). It was shown by Treu & Koopmans (2004) that the range of projected baryonic fraction within the Einstein radius is $`f_{\mathrm{bar}}=0.30.6`$. It is therefore crucial to include the gravitational pull of baryons in our simulations.

This paper is structured as follows. In Sect. 2 we first give the main properties of the N-body simulations that we use. We also introduce an improved smoothing scheme compared to Bradač et al. (2002) and describe how to extract lensing properties from N-body simulations. In Sect. 3 we focus on the cusp relation for simulated lens systems. Sect. 4 describes the modelling of synthetic images and the phenomenon of suppressed saddle points. We conclude and give an outlook in Sect. 5.

## 2 Strong lensing by a simulated galaxy

N-body simulations can provide a powerful benchmark for testing the effects of substructure on strong lensing. One can simulate conditions in which propagation effects due to the ISM can be ignored and thus examine only the signature of substructure. The drawback of this method at present lies in the resolution available for simulations that include dark matter, gas and star particles. This limits our analysis to mass clumps of $`10^8M_{}`$. However, since the mass resolution is improving rapidly, this will soon be less of a problem.

As in Bradač et al. (2002), we used the nummerical N-body simulations for several realisations of galaxies including gas-dynamics and star formation (Steinmetz & Navarro 2003). We investigate two different halos, each of them in three different projections. The simulations were performed using GRAPESPH, a code that combines the hardware N-body integrator GRAPE with the Smooth Particle Hydrodynamics (SPH) technique (Steinmetz 1996).

In Table 1 the properties of the halos are listed. In both cases the original simulated field contains approximately 300 000 particles. The simulation is contained within a sphere of diameter $`32\mathrm{Mpc}`$ which is split into a high-resolution sphere of diameter $`2.5\mathrm{Mpc}`$ centred around the galaxy and an outer low-resolution shell. Gasdynamics and star formation are restricted to the high-resolution sphere, while the dark matter particles of the low-resolution sphere sample the large scale matter distribution in order to appropriately reproduce the large scale tidal fields (see Navarro & Steinmetz 1997 and Steinmetz & Navarro 2000 for details on this simulation technique). From the original simulated field we use a cube of size $`200^3\mathrm{kpc}^3`$ centred on a single galaxy. This volume lies well within the high-resolution sphere and is void of any massive intruder particles from the low-resolution shell.

All simulations were performed in a $`\mathrm{\Lambda }`$CDM cosmology ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_\mathrm{b}=0.019/h^2`$, $`\sigma _8=0.9`$). They have a mass resolution of $`1.26\times 10^7M_{}`$ and a spatial resolution of $`0.5\mathrm{kpc}`$. A realistic resolution scale for an identified substructure is typically assumed to be $`40`$ particles which corresponds to $`5\times 10^8M_{}`$. The quoted mass resolution holds for gas/stars. The high-resolution dark matter particles are about a factor of 7 ($`=\mathrm{\Omega }_0/\mathrm{\Omega }_\mathrm{b}`$) more massive. A detailed analysis of the photometric and dynamical properties of the simulated halos was carried out in Meza et al. (2003) for the elliptical and Abadi et al. (2003a, b) for the disc galaxy.

Early N-body simulations suffered from the problem of overmerging, i.e. most satellites were dissolved due to the tidal field of the host galaxy. Current N-body simulations demonstrate that quite often this tidal disruption is inefficient resulting in surviving substructures. While this qualitatively different result is usually accounted to the increased particle numbers of modern simulations, it is in fact rather caused by a more prudent choice of the numerical softening compared with a more rigorous limit on the numerical time stepping. Convergence studies using modern N-body simulations demonstrate that with a sufficiently small numerical softening length and sufficiently rigorous numerical time stepping, mass functions can be accurately produced down to clumps with only a very few tens of particles (see e.g. Springel et al. 2001).

### 2.1 Delaunay triangulation smoothing technique

From the irregularly sampled particle distribution in the simulation box, we reconstruct the density field. We apply a smoothing procedure and then project the resulting particle distribution to obtain the surface mass density $`\kappa `$. In Bradač et al. (2002) it was shown that a more sophisticated smoothing method should be employed for the data analysis than simply smoothing with a Gaussian kernel. A method is needed that adapts the kernel size in order to increase the signal to noise of the reconstructed field. For this purpose, we make use of the Delaunay tesselation technique from Schaap & van de Weygaert (2000).

We perform a three-dimensional Delaunay triangulation using the QHULL algorithm (Barber et al. 1996). The density estimator from Schaap & van de Weygaert (2000) is then evaluated at each vertex, and we interpolate values of the density at each three-dimensional grid point to obtain the $`\kappa `$ map. The resulting density field is then projected onto a two-dimensional grid.

Since the N-body simulations contain three independent classes of particles (gas, stars, and dark matter particles, each having different masses), we repeat the procedure described above for each group separately and the final $`\kappa `$ map is obtained by adding the contributions from all three classes.

The Delaunay tesselation method performs very well in comparison to the standard Gaussian smoothing technique used in Bradač et al. (2002), or an adaptive Gaussian smoothing technique (see also Schaap & van de Weygaert 2000). Since it is non-parametric, it adjusts the scale of the smoothing kernel such that regions of low noise (i.e. where the particle density is highest) are effectively smoothed less than regions with high noise. Also the shape of the kernel is self-adaptive. Hence, this method is very useful for the analysis of galaxies with high dynamic range and significant structure in the mass distribution (e.g. mass clumps and spiral arms).

### 2.2 Estimating the noise properties

A drawback of using the Delaunay tesselation method is that the signal-to-noise evaluation for the final surface mass density map is non-trivial. For example, with a Gaussian kernel one can determine the noise by simply looking at the number of particles in a smoothing element (for a more detailed estimate, see Lombardi & Schneider 2002). When using the tesselation technique, such an approach is not viable.

One approach for estimating the error is to use bootstrapping (see e.g. Heyl et al. 1994). Normally we calculate physical properties by using all $`n`$ particles in the simulation. To create a bootstrap image one has to randomly select $`n`$ particles out of this simulation with replacement; i.e. some of the particles from the original simulation will be included more than once and some not at all. In other words, we randomly generate $`n`$ integers from 1 to $`n`$, representing the bootstrapped set of particles. If a particle is included $`k`$ times in a bootstrapped map, we put a particle at the same position with $`k`$ times its original mass. One can then make an error estimation using the ensemble of such images and calculating the desired physical quantity for each of them.

Whereas the tesselation itself is done very quickly, the interpolation of density on a grid is a process that takes a few CPU days on a regular PC. We therefore limit ourselves to 10 bootstrapped maps and the elliptical halo only.

For each pixel $`i`$ we calculate the associated error using

$$\sigma ^2\left(\kappa _i\right)=\frac{1}{M1}\underset{m=1}{\overset{M}{}}\left(\kappa _i^m\kappa _i\right)^2$$
(1)
where $`X_i`$ is the average value of $`X`$ at pixel $`i`$ averaged over all $`M`$ bootstrapped maps; in our case $`M=10`$. This procedure gives us an estimate for the error on $`\kappa `$. We find that within the critical curves, i.e. $`\kappa 0.25`$ the average noise $`\sigma (\kappa _i)/\kappa `$ is of the order of $`5\%`$.

### 2.3 Strong lensing properties

Having obtained the $`\kappa `$-map, we then calculate the lens properties on the grid ($`2048\times \mathrm{\hspace{0.25em}2048}`$ pixels). The Poisson equation for the lens potential $`\psi `$

$$^2\psi (𝜽)=2\kappa (𝜽)$$
(2)
is solved on the grid in Fourier space with a DFT (Discrete Fourier Transformation) method. We used the publicly available C library FFTW (“Fastest Fourier Transform in the West”) written by Frigo & Johnson (1998). To reduce the boundary effects, padding was introduced. In particular, the DFT was performed on a $`4096\times \mathrm{\hspace{0.25em}4096}`$ grid, where one fourth of the grid contains the original $`\kappa `$ map and the rest is set to zero.

One can now proceed in two ways. Either, one can calculate the two components of the shear $`\gamma _{1,2}`$ and the deflection angle $`𝜶`$ by multiplying the potential in Fourier space by the appropriate kernel. However, since we calculate the Fourier transform of the potential on a finite grid, we filter out high spatial frequency modes. By multiplying the transform with different kernels for $`\gamma _{1,2}`$ and $`𝜶`$, these final maps do not correspond exactly to the same $`\kappa `$ map. The effect is small, but it shows up near the critical curves. Therefore, it is better to only calculate the lens potential $`\psi `$ using DFT and obtain the shear and deflection angle by finite differencing. The latter method is also less CPU-time and memory intensive.

The simulated galaxies are a field galaxy. However, most of the lenses in quadruple image systems are members of groups. To make our simulated galaxies more closely resemble a realistic system, we add an external shear to the Jacobi matrix (evaluated at each grid point). The shear components are the same for each projection and all halos. We use

$$\gamma _1^{\mathrm{ext}}=0.04,\gamma _2^{\mathrm{ext}}=0.16;$$

corresponding to the shear of the best-fit singular-isothermal ellipse model with external shear for the lens B1422$`+`$231 in Bradač et al. (2002). Figure 3 show the magnification maps of the elliptical and edge-on disk galaxy. The corresponding caustic curves (for a source at $`z=3`$) were obtained by projecting the points with vanishing determinant of the Jacobian matrix from the image plane to the source plane. The critical curves are plotted in white in Fig. 3a for the elliptical and Fig. 3b for the edge-on disk galaxy. The corresponding caustic curves are plotted in Fig. 3a and 5c respectively.

### 2.4 The importance of baryons

As mentioned above, the influence of baryons is very important in lens galaxies. Since we are interested in the effects of substructure, it would be desirable to use N-body simulations that have the highest possible dynamical range. At present this is achieved in high-resolution N-body simulations that only include dark matter. However, if baryons are not included, the central potential is more shallow than what we typically observe in lens galaxies. All quadruply-imaged systems for which the inner slope of the mass distribution has been measured, are well described by a total mass density profile $`r^\gamma `$ with $`\gamma 2`$ (Kochanek 1995; Cohn et al. 2001; Treu & Koopmans 2002a, b; Koopmans & Treu 2003; Treu & Koopmans 2004), consistent with the combined mass distribution of dark matter and baryons seen in the simulations. Hence, dark matter only simulations do not accurately represent the overall properties of lens galaxies, and instead we need to use hydro-dynamical simulations.

To test the importance of baryons we have simulated an elliptical halo using only dark matter particles and performed the same lensing analysis as described above. We project the density approximately along the long axis of the halo, thus maximising the central density. In Fig. 2 we plot the corresponding caustic curves for the halo simulations; in panel (a) the simulation including baryons (DM+B), for panel (b) we use the DM-only simulation (DM). The radial critical curve in the DM+B halo is located close to the galaxy centre and is therefore not well resolved (but also irrelevant for our purposes) - the corresponding caustic is therefore not plotted.

The two caustic curves are very different, indicating very different overall strong lensing properties of these two simulations. Whereas the DM+B simulation has a steep inner mass profile (very close to singular isothermal), the DM simulation has a caustic configuration typical for a lens with a shallow density profile (see e.g. Wallington & Narayan 1993). The radial critical curve has become smaller compared to the DM+B halo and the corresponding caustic curve is almost entirely enclosed by the asteroid caustic. The prominent naked cusp region is a three-image region. This configuration is extremely rare among the observed lensed systems. For a source located in such a region one would observe three highly magnified images. There is only one possible example of a triply imaged quasar out of $`50`$ doublets and triplets (Evans & Hunter 2002), namely APM 08275+5255 (Ibata et al. 1999). Further, the vast majority of systems similar to B1422+231 can not form in such a potential. We therefore conclude that the effects of baryons have to be included, and we will use only DM+B simulations from now on, discussing their limits where necessary.

## 3 The cusp relation

We generate different four-image systems using each simulated galaxy. For each projection, regions in the source plane where five images form are determined. The image plane is projected back to the source plane using the magnification and deflection angle maps. We used the grid search method from Blandford & Kochanek (1987) to find the pixels enclosed by the asteroid caustic and approximate image positions. Then the MNEWT routine from Press et al. (1992) is applied and the lens equation is solved for the image positions. For this step, we interpolate the deflection angle between the grid points. We use bilinear and bicubic spline interpolation, and both methods give comparable results. Once we have the image positions, their magnifications are calculated and the four brightest images are chosen. These represent the “observable” images; the fifth image is usually too faint (except in the regions where more than five images are formed) and therefore likely to escape observation.

We can now investigate the basic properties of the synthetic four-image systems. There are three basic configurations: the fold, cusp, and cross. They correspond to a source located inside the asteroid caustic, close to a fold, a cusp or near the centre, respectively (see e.g. Keeton 2001). All configurations have been observed, and even though one would naively think that fold and cusp images are rare among observed lenses, they are in fact frequently observed due to the large magnification bias. In this section we will mainly concentrate on cusp image configurations.

The behaviour of gravitational lens mapping near a cusp was first studied by Blandford & Narayan (1986), Blandford (1990), Schneider & Weiss (1992), and Mao (1992), who investigated the magnification properties of the cusp images and concluded that the sum of the signed magnification factors of three merging images approaches zero as the source moves towards the cusp. In other words,

$$R_{\mathrm{cusp}}=\frac{\left|\mu _\mathrm{A}+\mu _\mathrm{B}+\mu _\mathrm{C}\right|}{\left|\mu _\mathrm{A}\right|+\left|\mu _\mathrm{B}\right|+\left|\mu _\mathrm{C}\right|}0,\mathrm{for}\mu _{\mathrm{tot}}\mathrm{},$$
(3)
where $`\mu _{\mathrm{tot}}`$ is the unsigned sum of magnifications of all four images, and A, B and C is the triplet of images forming the smallest opening angle (Keeton 2001). The opening angle is measured from the galaxy centre and is spanned by the two images of equal parity. The third image lies inside this opening angle.

### 3.1 The cusp relation of an N-body simulated elliptical galaxy

The cusp relation (3) is an asymptotic relation and holds when the source approaches the cusp from inside the asteroid. One can derive the properties of lens mapping close to critical curves using a Taylor expansion of the Fermat potential around a critical point (see e.g. Schneider & Weiss 1992). Such calculations are very cumbersome and therefore it is difficult (if not impossible) to explore the influence of arbitrary substructure analytically. In practice, we can calculate $`R_{\mathrm{cusp}}`$ for the N-body simulated systems. Smoothing the original $`\kappa `$-map on different scales then gives an indication of the influence of substructure on $`R_{\mathrm{cusp}}`$.

The three cusp images \[designated as A, B and C in (3)\] are chosen according to the image geometry. Since we know the lens position, this procedure is straightforward and foolproof. We have identified the triplet of images belonging to the smallest opening angle (described above). Since we know the image parities and magnifications, one is tempted to identify the three brightest images as the cusp images and assign different parity to the brightest one than to the other two (e.g. as in Möller et al. 2003). However, due to the presence of shear and substructure this could lead to misidentifications.

Figure 3a shows the caustic curve in the source plane for the simulated elliptical at a redshift of $`z_\mathrm{l}=0.81`$. The source is at a redshift of $`z_\mathrm{s}=3`$. Approximately 30 000 lens systems are generated with source positions inside the asteroid caustic. $`R_{\mathrm{cusp}}`$ is plotted in gray-scale. The apparent discontinuities originate from different image identification. In the very centre of the caustic the meaning of “cusp image” is ill defined. As the source moves in the direction of the minor or major axes we chose different subsets of three cusp images and therefore the discontinuity arises.

The remaining panels of Fig. 3 show the effect of smoothing the small-scale in the surface mass density $`\kappa `$ map with Gaussian kernels characterised by standard deviation $`\sigma _\mathrm{G}`$. The values for $`\sigma _\mathrm{G}`$ were chosen to be $`\sigma _\mathrm{G}1,2,5\mathrm{kpc}`$ for panels (b), (c) and (d) respectively. Note that we do not smooth the $`\kappa `$-map directly. First we obtain the smooth model $`\kappa _{\mathrm{smooth}}`$ for the $`\kappa `$-map by fitting elliptical contours to the original map using IRAF.STSDAS package ellipse. We subtract $`\kappa _{\mathrm{smooth}}`$ from $`\kappa `$ and smooth the difference using different Gaussian kernels. We then add the resulting map back to the $`\kappa _{\mathrm{smooth}}`$. In this way the overall radial profile of the mass distribution is not affected.

The effect of smoothing on the cusp relation is clearly visible. In Fig. 3d one sees that the substructure is completely washed out when smoothing on scales of $`\sigma _\mathrm{G}5\mathrm{kpc}`$ is applied. As we go to smaller smoothing scales, the effects of substructure become clearly visible. In the extensions of swallowtails there is a region where the cusp relation is strongly violated (with $`R_{\mathrm{cusp}}0.6`$, where the smooth model predicts $`R_{\mathrm{cusp}}0.1`$). However, further out, a swallowtail can cause the cusp relation to change the trend and go to zero (due to high-magnification systems being formed in such region).

Finally, the cusp relation behaves differently for the source on the major or the minor axis (see especially Fig. 3d). This is a generic feature for smooth elliptical models and can easily be calculated for e.g. an elliptical isopotential model (see e.g. Schneider et al. 1992). We use this model since it is analytically tractable for source positions along the major and minor axis. In Fig. 4 we plot the cusp relation for the source moving along the major (minor) axis as a thick (thin) solid line for the elliptical isopotential model with $`ϵ=0.15`$. As the source approaches the cusp, $`R_{\mathrm{cusp}}0`$ for both source positions, however the slope is different. We also plot the total magnification factor of the three cusp images, i.e. $`\mu _{\mathrm{A}+\mathrm{B}+\mathrm{C}}=\left|\mu _\mathrm{A}\right|+\left|\mu _\mathrm{B}\right|+\left|\mu _\mathrm{C}\right|`$ as a thick(thin) dashed line.

### 3.2 The cusp relation in an N-body simulated disk galaxy

The procedure described above was applied to three different projections of the elliptical halo and very similar conclusions can be drawn. However, the question is, how much these conclusions depend on the specific morphological type of the galaxy. To investigate this we follow a similar procedure for a simulated disk galaxy.

In this case, however, we did not look at the effects of additional smoothing. It is difficult to subtract a smooth model, since the galaxy consists of a bulge, warped disk and extended halo, which can not simply be fitted by ellipses. If we were to smooth the edge-on projection with a Gaussian kernel, we would also wash out the (in our case warped) disc. We only include the analysis of the cusp relation for this halo in order to show the effect of the disc on $`R_{\mathrm{cusp}}`$.

Figure 5 shows the cusp relation of an N-body simulated disk galaxy in three projections. The disc clearly plays a role for the edge-on projection (see Fig. 5b,c), whereas the face-on projection is similar to the elliptical galaxy. Especially in Fig. 5c, where the disc extent in projection is smaller than the typical image separation, the cusp relation in the upper-right and lower-left cusp is strongly violated. This direction corresponds to the orientation of the disc.

### 3.3 Observed cusp relation

In this section we compare predictions of $`R_{\mathrm{cusp}}`$ from our simulated galaxies with the values from observed lens system. Unfortunately the number of observed systems is not large; seventeen four-image systems have been published (see Keeton et al. 2003 for a summary), and four of them are believed to show a typical cusp geometry (see Table 2).

Note that relation (3) is model independent and can only be broken in the presence of substructure on scales smaller than the image separation. Hence, if the smooth model provides an adequate description of the lens galaxy, one would expect $`R_{\mathrm{cusp}}0`$ for these lenses. This is clearly not the case and for this reason it is difficult to explain their flux ratios using simple, smooth models.

If we make a comparison with simulations, one can see that the pattern in Fig. 3d clearly does not explain the observed $`R_{\mathrm{cusp}}`$ of these four lenses. On the other hand, substructure on few $`\mathrm{kpc}`$ scales and below provides enough perturbations to $`R_{\mathrm{cusp}}`$ to explain the observed values.

The question arises, however, whether from the value of $`R_{\mathrm{cusp}}`$ we can conclude that we indeed see the effects of CDM substructure. Keeton et al. (2003) indeed argue that the cusp relation alone does not reveal anomalous flux ratios in B1422$`+`$231. Still, detailed modelling for this system shows that the flux ratios are anomalous. The difficulties in modelling B1422$`+`$231 are not only a consequence of a violation in cusp relation, but also that image D is a fainter than predicted from smooth models. On the other hand, the simulated disk galaxy shows violation of $`R_{\mathrm{cusp}}`$ even though there are no clear mass clumps present in the region where images are formed. Hence one has to be cautious when making conclusions about the presence of CDM substructure based on the value of $`R_{\mathrm{cusp}}`$ alone.

### 3.4 The influence of noise on the cusp relation

The simulated cusp-relation can be reliably compared with observations only if we know the noise properties of our simulations. The effects of noise and physical substructures need to be disentangled through a detailed analysis.

From the bootstrapping procedure (Sect. 2.2) we also get an estimate of the error in the cusp relation. We estimated the noise on $`R_{\mathrm{cusp}}`$ using a similar technique as for $`\kappa `$ in (1). Moreover, we do not perform the analysis directly in the source plane by subtracting the maps pixel-by-pixel. The problem is that bootstrapping somewhat changes the shape of the caustic curve (see also Fig. 7). Since we never observe the source plane directly, in reality we can not distinguish between two shifted, but identical caustics. We therefore have to match different bootstrap maps in the image plane. We compare the image positions generated by each source in the original frame with those generated by bootstrapped lenses. Thus for each source position in the original map (see Fig. 3a), we search for the source position in the bootstrapped map such that the four image positions from both maps differ as little as possible.

In principle, one can redo the ray-tracing and calculate image fluxes for the position in the original map using the bootstrapped-lens properties. However this is not necessary, since our grid in the source plane is fine enough and we only need an approximate error estimate. Figure 6 shows the estimated absolute error $`\sigma (R_{\mathrm{cusp},i})`$ for the elliptical halo from Fig. 3a. As described above, each value of $`\sigma (R_{\mathrm{cusp},i})`$, plotted in grey-scale, does not correspond to the error for a source position, but refers to the error of the system with image positions similar to the original ones.

The absolute error becomes slightly larger in the regions of the swallow tails. This is, however, not the effect of substructure vanishing in individual bootstrapped realisations. It is rather the effect of slight position changes of individual clumps. If one looks at the individual caustic curves of bootstrapped halos (see Fig. 7) the swallow-tails are present in all realisations, although they change their positions slightly. Since this hardly affects the image positions we cannot perfectly match the source positions $`i^{}`$ with the original source position $`i`$ in these regions; thus we are overestimating the true error.

We conclude that the values of $`R_{\mathrm{cusp}}`$ in the close proximity of the cusp can be as high as $`R_{\mathrm{cusp}}=0.6`$, with the error of $`\pm 0.1`$ as determined from the bootstrapping. Smoothing the substructure on scales as large as $`1\mathrm{kpc}`$ does not remove this effect, but reduces it. This is expected, since smoothing changes the profile of substructures.

Finally, we investigate how well we can sample the smooth mass distribution given the number of particles in the original simulation. We take an ellipsoidal power-law density profile following $`\rho (r)r^{2.9}`$. The power-law index was chosen, to closely reproduce the surface mass density $`\kappa `$ of the simulation which follows $`\kappa (\theta )\theta ^{1.9}`$ in the vicinity of the critical curve (see Sect. 4.1 for more details). The number of particles we use is the same as in the original N-body simulation. To each particle we assign the average mass of the original sample, leaving the total mass of the lens galaxy unchanged. We sample the density profile using a rejection method (e.g. Press et al. 1992). The resulting particle distribution was again adaptively smoothed using Delaunay tesselation and we perform the lensing analysis as in all previous cases.

The resulting cusp relation is given in Fig. 8a. We have chosen the profile such that the particle densities around the outer critical curves are similar in both cases. Only in that case can we reliably compare the noise properties of $`R_{\mathrm{cusp}}`$. The fact that the caustic is larger than in the case of the N-body simulated elliptical halo (compare to Fig. 3a) is here of lesser importance; it arises due the difference in the central profile, far from the critical curve.

The absence of strong violations of $`R_{\mathrm{cusp}}`$ close to the caustic in Fig. 8 as compared to Fig. 3a confirms that deviation of $`R_{\mathrm{cusp}}`$ from zero is not dominated by the shot noise of the particles, but is due to genuine substructure in the N-body simulation. For a more quantitative comparison we show the probability distribution of $`R_{\mathrm{cusp}}`$ for systems with $`\mu _{\mathrm{tot}}>20`$ in Fig. 8b as a solid line for the sampled smooth halo and as a dashed line for the original. The much tighter distribution for the sampled halo confirms the absence of strong violations of $`R_{\mathrm{cusp}}`$ in the sampled smooth halo compared to the original simulation.

This analysis also shows the advantage in smoothing with Delaunay tesselation. If we only use a Gaussian kernel (of the same size as in Fig. 3b), the deviations in the cusp-relation for the sampled particle distribution are much larger.

## 4 Lens models for synthetic image systems

To explore the saddle point demagnification phenomena we fit the synthetic image systems using a singular isothermal ellipse model with external shear (SIE+SH) (Kormann et al. 1994). We do not include the flux ratios in the fit, since they are affected by the substructure. Furthermore, we keep the lens position fixed. Using only image positions we thus have 7 parameters and 8 constraints (the parameters that we used are lens strength, two components of the ellipticity of the lens, two external shear components and the source position. The 4 image positions provide 8 constraints).

### 4.1 Surface mass density profile

We find that the average unsigned magnifications predicted by the SIE+SH model are higher than those from N-body simulations. This is true for all systems generated with the same halo; a consequence of the mass profile being steeper than isothermal around the typical position where the images are formed.

We have calculated the $`\kappa `$ profile for the N-body simulated elliptical by fitting the $`\kappa `$ with elliptical contours using the IRAF.STSDAS package ellipse. For the inner $`500\mathrm{pc}`$ the profile is very close to isothermal with slope $`1.0\pm 0.2`$, the profile becomes steeper than isothermal with slope $`1.9\pm 0.2`$. The break radius, where the profile becomes steeper than isothermal, is smaller than the radius where images form.

We therefore over-predict the magnifications using an isothermal profile. In order to deal with this problem, one would preferably use a power-law profile with the index mentioned above for the lens modelling, instead of SIE. However, these profiles require the evaluations of hypergeometric functions and/or numerical integrations (Grogin & Narayan 1996; Barkana 1998). Modelling is therefore computationally too expensive to apply to several 10 000 sources. Moreover, this is not critical, since we are not searching for the best fit macro-model, but rather pretending that we observe the systems and try to fit them with the model used for most observed lenses. Since we only observe the flux ratios and not directly the magnifications, it is impossible to compare magnification factors in practice. Thus one cannot spot the difference in profiles using only this consideration when dealing with real lens systems.

### 4.2 Suppressed saddle points

It was first noted by Witt et al. (1995) that the expected flux changes due to stellar-mass perturbers differ between saddle points and minima images. This was further investigated by Schechter & Wambsganss (2002) who conclude that for fold images a stellar-mass component added to a smooth mass distribution can cause a substantial probability for the saddle point image to be demagnified w.r.t. the minimum.

Recently, Kochanek & Dalal (2003) and Keeton (2003) investigated the effect of singular isothermal sphere (SIS) mass clumps on the flux perturbations for images of different types. Their conclusions are similar; SIS perturbers also cause the brightest saddle point image to behave statistically differently compared to those at minima.

Knowing whether the flux anomalies depend upon the image parity for observed lenses would be a major step forward in identifying flux anomalies either with substructure, or with propagation effects. Namely, if the observed flux anomalies depend upon the image parity and its magnification we can set limits on the influence of the ISM (see Kochanek & Dalal 2003).

To test whether saddle point images are also suppressed in our simulation we proceed as follows. Following Kochanek & Dalal (2003) we define $`C\left(<\mathrm{ln}(\mu _{\mathrm{obs}}/\mu _{\mathrm{mod}})\right)`$ to be the cumulative probability distribution of flux residuals, where $`\mu _{\mathrm{obs}}`$ is the “observed” magnification of a simulated image and $`\mu _{\mathrm{mod}}`$ is the magnification predicted by the best-fitting smooth SIE+SH model (as described above). In Fig 9a we plot $`C\left(<\mathrm{ln}(\mu _{\mathrm{obs}}/\mu _{\mathrm{mod}})\right)`$ for the systems with “observed” total unsigned magnifications $`\mu _{\mathrm{tot}}>20`$ for the simulated elliptical halo (see also Fig. 3a). We chose systems that have highly magnified images, because they are affected by substructure at most. We repeat the whole procedure with the same halo, but smoothed on a scale of $`\sigma _\mathrm{G}5\mathrm{kpc}`$; in this case essentially no difference is seen among images of different parity, as expected.

For the original halo the cumulative distribution is much broader for the brighter minimum and saddle point. This is in accordance with the conclusions from Mao & Schneider (1998); the higher the magnification, the more is the image flux affected by substructure. Among the two most magnified images, the saddle point is on average more demagnified compared to the brighter minimum. We have examined two other orthogonal projections of the mass density of the same halo and we find that qualitatively the results do not change with projection.

The effect, however, is not as pronounced as in Kochanek & Dalal (2003). The reason is two-fold. First, our simulations have a resolution of $`10^8M_{}`$ and structure at this scale and below is suppressed when using Delaunay tesselation. Kochanek & Dalal (2003), however, used SIS clumps with masses of $`10^6M_{}`$. Since these are more numerous, they can enhance the effect. Further, fitting SIE to the global mass profile is not fully justified. The mass profile is known only for a handful of observed lenses. Whereas the lens galaxies in MG1654$`+`$134, MG2016$`+`$112, 0047$``$281, and B1933$`+`$503 have a nearly isothermal profile, the one of PG1115$`+`$080 seems to be steeper (see Kochanek 1995; Cohn et al. 2001; Treu & Koopmans 2002a, b; Koopmans & Treu 2003). Besides the absence of substructure on scales $`10^7M_{}`$, our synthetic systems and their modelled quantities closely resemble the properties of realistic lenses and the way in which these are modelled.

In the analysis of Kochanek & Dalal (2003) the synthetic images were generated using an SIE macromodel with SIS substructure. This simplifies the model fitting and explains why they get a transition of cumulative distribution exactly at $`\mathrm{ln}(\mu _{\mathrm{obs}}/\mu _{\mathrm{mod}})0`$.

In addition, we have looked at the cumulative flux mismatch distribution in the bootstrapped maps, to investigate the significance of our results. In the bootstrapped images we confirm the broader distribution for bright minimum and saddle point images. In Fig 9b we plot the difference of the cumulative distribution of flux residuals $`\mathrm{\Delta }C\left(<\mathrm{ln}(\mu _{\mathrm{obs}}/\mu _{\mathrm{mod}})\right)`$ between the brightest saddle point and the brightest minimum images for the original halo (solid line with dots), bootstrapped images (solid lines) and halo when additionally smoothed on a scale of $`\sigma _\mathrm{G}5\mathrm{kpc}`$ (dashed line). Positive values of $`\mathrm{\Delta }C`$ thus denote the saddle point demagnification. In all bootstrapped images $`\mathrm{\Delta }C`$ is positive, except for few points (corresponding to $`C\left(<\mathrm{ln}(\mu _{\mathrm{obs}}/\mu _{\mathrm{mod}})\right)1`$). This confirms that the effect of the saddle point demagnification is present and comparable to the original halo, whereas in the smoothed halo this effect is not seen. We conclude that substructure on mass scales $`10^8M_{}`$ significantly contributes to the saddle point demagnification; possibly going a long way in explaining the observed saddle point demagnification.

## 5 Conclusion

In this work we have studied strong gravitational lensing properties of N-body simulated galaxies. In particular, we concentrated on the influence of substructure on flux ratios on highly magnified images. Such an analysis is crucial in order to fully understand the lensing signal that we observe in realistic lenses and to disentangle the influence on lensed-image fluxes due to propagation effects in the ISM and mass substructure.

We have examined two strong lensing signatures of substructure, i.e. the broken cusp relation observed in images that show a typical cusp configuration and the saddle point suppression. The saddle point suppression has been previously studied using semi-analytic prescriptions of substructure (e.g. Schechter & Wambsganss 2002; Kochanek & Dalal 2003). The effect of substructure on the cusp relation, however, has up to now not been studied in detail.

In order to determine the magnitude of both effects we use N-body simulated galaxies. The difference compared to the works of Schechter & Wambsganss (2002) and Kochanek & Dalal (2003) is that we are using a representation of substructure that is as realistic as possible and do not make any assumptions on the mass function and abundance of sub-halos. The drawback, however, is the resolution of the simulations. We are therefore not able to extrapolate the analysis to masses $`<10^8M_{}`$. Still, the signatures of both effects are clearly present, and in the future we plan to use higher-resolution N-body plus gasdynamical simulations to explore their effects in greater detail.

The main question when dealing with N-body simulations is how much are the magnification factors, that we use for synthetic image systems, affected by noise which can mimic substructure of $`<10^8M_{}`$. We show that the average relative noise in the surface mass density $`\sigma (\kappa )/\kappa `$ lies below the $`5\%`$ level for $`\kappa 0.25`$. Second, in Sect. 3.4 we show that the results for $`R_{\mathrm{cusp}}`$ are not significantly affected by the noise, and are dominated by physical substructure. The signal is dominated by several resolved mass clumps, which in projection lie close to the Einstein radius. Similarly, in the case of the phenomenon of suppressed saddle points, the bootstrap analysis shows that the signal also here is not dominated by noise.

The behaviour of $`R_{\mathrm{cusp}}`$ for sources close to a cusp is a very promising tool to detect substructure. Its main advantage is that it makes definite, model-independent predictions for the image magnifications. These predictions can only be broken in the presence of structure in the potential on scales smaller than the image separation. In Fig. 3, where we used a simulated elliptical galaxy to calculate $`R_{\mathrm{cusp}}`$, we clearly see these effects. When smoothing the substructure on larger scales we witness the transition to the pattern that is common for generic smooth SIE lens model.

However, the disc in the disk galaxy can also help to destroy the cusp relation for sources in the vicinity of a cusp. We have calculated the cusp relation pattern for the simulated disk galaxy, and even in the absence of obvious substructure in the form of clumps we can see strong violations of the cusp relation. This is expected, since the edge-on disk gives $`\kappa `$-variations on the scales smaller than the image separation and , similar to small-scale substructure. One cannot conclude from a broken cusp relation alone that we observe the signatures of mass substructure in the form of clumps. However, the observations show that most observed lenses are elliptical and therefore one can concentrate on this morphology. Still, detailed modelling is required in most cases (e.g. B1422$`+`$231) to clearly see the effect.

The phenomenon of suppressed saddle points is a very strong prediction that rules out a significant influence of the ISM on flux anomalies. If flux anomalies depend on parity and magnification, they must clearly be caused by lensing, if significant substructure is present. Observations so far, show a clear parity dependence, which is more obvious for highly magnified images.

Finally, our analysis shows that the two brighter images are more affected by substructure than the two fainter ones. In addition, we confirm that the brightest saddle point image in N-body simulated systems has a higher probability to be demagnified, in accordance with predictions from microlensing and from semi-analytic work by (Kochanek & Dalal 2003). It is therefore necessary that all mass scales are properly accounted for, in order to compare observations with theoretical predictions in detail.

For future work, we plan to look for jet curvature as seen from N-body simulations. At present, there is only one case of a curved jet observed that is likely the cause of gravitational lensing (Metcalf 2002). It will be interesting, however, to investigate the probability of more of these occurrences. We plan to investigate the signal one expects on average for multiple-imaged jets; this signal is also less affected by noise in the simulation and low-mass substructure.

In summary, gravitational lensing remains a very powerful tool for testing the existence of CDM substructure. N-body simulated galaxies do seem to produce the same effects as seen in observed lens systems. In addition, systematics on flux anomalies (scatter broadening, scintillation, microlensing) can be efficiently ruled out by multi-frequency and higher-frequency observations of lenses. Furthermore, the statistical analysis of large samples of lenses can directly probe the properties of CDM substructure in galaxies to a redshift of $`z1`$. This provides a unique tool to measure the evolution of these structures with cosmic time, as predicted in the hierarchical structure formation scenario.

###### Acknowledgements.

We would like to thank to our referee for his constructive comments, Oleg Gnedin, Chris Kochanek and Oliver Czoske for many useful discussions that helped improve the paper, Vincent Eke for the initial conditions of the simulations, and Richard Porcas for careful reading of the manuscript. This work was supported by the International Max Planck Research School for Radio and Infrared Astronomy at the University of Bonn, by the Bonn International Graduate School, by the Deutsche Forschungsgemeinschaft under the project SCHN 342/3–1, and by the NASA ATP program under grant NAG 5-10827.