text:
# The specific heat and optical birefringence of 𝐹𝑒_0.25𝑍𝑛_0.75𝐹₂

## Abstract

The specific heat $`(C_m)`$ and optical birefringence $`(\mathrm{\Delta }n)`$ for the magnetic percolation threshold system $`Fe_{0.25}Zn_{0.75}F_2`$ are analyzed with the aid of Monte Carlo (MC) simulations. Both $`\mathrm{\Delta }n`$ and the magnetic energy $`(U_m)`$ are governed by a linear combination of near-neighbor spin-spin correlations, which we have determined for $`\mathrm{\Delta }n`$ using MC simulations modeled closely after the real system. Near a phase transition or when only one interaction dominates, the temperature derivative of the birefringence $`[d(\mathrm{\Delta }n)/dT]`$ is expected to be proportional $`C_m`$ since all relevant correlations necessarily have the same temperature dependence. Such a proportionality does not hold for $`Fe_{0.25}Zn_{0.75}F_2`$ at low temperatures, however, indicating that neither condition above holds. MC results for this percolation system demonstrate that the shape of the temperature derivative of correlations associated with the frustrating third-nearest-neighbor interaction differs from that of the dominant second-nearest-neighbor interaction, accurately explaining the experimentally observed behavior quantitatively.

Measuring the linear optical birefringence $`(\mathrm{\Delta }n)`$ in anisotropic, antiferromagnetic crystals undergoing magnetic phase transitions is a powerful way of determining the magnetic specific heat ($`C_m`$) critical behavior. It has been shown that the temperature derivative of the optical birefringence $`[d(\mathrm{\Delta }n)/dT]`$ is proportional to $`C_m`$. In many cases, the birefringence technique has provided the most precise experimental determinations of universal critical behavior parameters in pure and randomly mixed and dilute magnetic Ising systems. For the case of the three-dimensional $`(d=3)`$ random-field Ising model (RFIM), which applies for a dilute antiferromagnet with an applied field along the ordering direction, birefringence measurements yielded evidence of a phase transition. The advantages of the technique are threefold: the technique is typically much easier to employ than traditional heat pulse techniques; the effects of concentration gradients inevitably present in mixed and dilute crystals can be greatly reduced; and the phonon contributions to the specific heat are greatly suppressed in the birefringence. Since the transition typically varies with concentration in mixed and dilute systems, the critical behavior is often masked by the concentration gradients quenched into the system during growth. The laser beam used in the optical technique can be aligned perpendicular to the gradient, often reducing the gradient effects by an order of magnitude. This has been crucial in the study of random-exchange Ising model (REIM) and RFIM systems in $`d=2`$ and $`3`$. The virtual elimination of the phonon background has allowed detailed analysis of the specific heat in $`d=1`$ and $`2`$ systems. For $`d=2`$ this has allowed a detailed comparison of the magnetic specific heat of the pure system and the Onsager solution to the $`d=2`$ Ising model and, for dilute systems, a scaling analysis of the destruction of the phase transition by random fields in the $`d=2`$ Ising model. For $`d=1`$, the technique has been used successfully to determine the exchange constants in systems where the phonon background overwhelms the magnetic contributions to conventional specific heat data.

That $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are proportional in certain cases is a consequence of the fact that both $`\mathrm{\Delta }n`$ and the magnetic energy $`(U_m)`$ are governed by a linear combination of near-neighbor spin-spin correlations . This proportionality has been experimentally shown explicitly for the $`d=3`$ pure Ising system $`FeF_2`$ as well as its site-diluted counterparts $`Fe_{0.93}Zn_{0.07}F_2`$ and $`Fe_{0.46}Zn_{0.54}F_2`$ by comparing directly the pulsed heat and optical data. This proportionality has been demonstrated both in zero applied field ($`H=0`$) and in applied fields. The magnetically dilute crystals of $`Fe_{1x}Zn_xF_2`$ have been the most extensively studied realizations of the RFIM, an important class of systems with randomness and frustration. Such experimental verification has been crucial in this controversial field since some authors have expressed reservations about the proportionality for $`H>0`$, i.e. the RFIM case.

Despite the experimental verification of the proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$, it is important to verify the details of the mechanism for the proportionality near phase transitions in a specific case. One opportunity to analyze the detailed relationship between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ presents itself in the compound $`Fe_{0.25}Zn_{0.75}F_2`$, which, with respect to the dominant exchange interaction, $`J_2`$, is at the magnetic percolation threshold, the concentration below which no ordering is possible. Throughout the rest of this work we will refer to the percolation threshold concentration as the appropriate one if only $`J_2`$ were being considered. The presence of other interactions make the concept of percolation more complicated and they can play an important role in the ordering processes even when they are very small. Interestingly, $`d(\mathrm{\Delta }n)/dT`$ data for $`Fe_{0.25}Zn_{0.75}F_2`$ exhibit a change in sign near $`T=8`$ K. While the high temperature behavior of $`d(\mathrm{\Delta }n)/dT`$ is accurately proportional to $`C_m`$, the sign reversal in $`d(\mathrm{\Delta }n)/dT`$ at low temperature is not reflected, of course, by the behavior of $`C_m`$. We present data obtained from both the optical and pulsed heat techniques that show the nonproportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ at low $`T`$, as shown in Fig. 1. We also present Monte Carlo (MC) simulation results that can be used to explain why the proportionality breaks down between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ at low temperatures. $`U_m`$ and $`\mathrm{\Delta }n`$ each depend on the near-neighbor spin-spin correlations, but with different weights. We show that whereas $`U_m`$ is dominated by the second-nearest-neighbor correlation, $`\mathrm{\Delta }n`$ is predominantly governed by the third-nearest-neighbor correlation that shows a different temperature dependence below $`15`$ K. The results of this study are consistent with theory regarding the proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ near a phase transition, where the different spin-spin correlations must necessarily have the same temperature dependence. We also demonstrate that an applied field has a predictable effect on both $`U_m`$ and $`\mathrm{\Delta }n`$.

The compound $`Fe_{1x}Zn_xF_2`$ is ideal for RFIM studies. Pure $`FeF_2`$ is well modeled by the Hamiltonian

$$=\underset{l}{}\underset{i<j}{}J_l\stackrel{}{S_i}\stackrel{}{S_j}+D\underset{i}{}\stackrel{}{S_i}^2.$$
(1)
The $`\stackrel{}{S}=2`$ system has the interaction strengths $`J_1=0.069`$ K, $`J_2=5.278`$ K, and $`J_3=0.279`$ K, as determined from neutron inelastic scattering measurements. The large crystal-field anisotropy $`D=9.29`$ K persists while the exchange interactions remain largely constant as the magnetic spins are diluted in $`Fe_{1x}Zn_xF_2`$, making this an excellent Ising system for all magnetic concentrations $`x`$. Single crystals can be grown for all $`x`$ with very small concentration gradients  and with superb structural quality. The magnetic interactions are dominated by the antiferromagnetic second-nearest-neighbor super-exchange interaction, $`J_2`$, between the body-center and body-corner ions. All other interactions are negligible except near the percolation threshold concentration, where a small frustration due to the third-nearest-neighbor interaction between ions of the same sub-lattice along the direction perpendicular to the spin-ordering $`c`$-axis, $`J_3`$, becomes important. At low temperatures, the Heisenberg character of the spin-spin interaction is less important and we might expect $`U_m`$ to be fairly represented, for $`H=0`$, by the REIM Hamiltonian,

$$=\underset{l}{}\underset{i<j>}{}J_lϵ_iϵ_jS_iS_j,$$
(2)
where $`S_i=\pm 2`$, $`ϵ_i=1`$ if site $`i`$ is occupied and zero otherwise, and $`J_l`$ is the strength of the $`l^{\mathrm{th}}`$ nearest-neighbor interaction. $`C_m`$ is then given by a sum of temperature derivatives of spin-spin correlations,

$$C_m=\frac{d(U_m)}{dT}=\underset{l}{}\underset{i<j}{}J_lϵ_iϵ_j\frac{d(<S_iS_j>_l)}{dT},$$
(3)
where the $`J_l`$ are assumed independent of temperature. The temperature dependent birefringence $`(\mathrm{\Delta }n)`$, the difference between the indices of refraction along the spin ordering axis and perpendicular to it, depends only on a linear combination of the same correlation functions. The temperature derivative yields

$$\frac{d(\mathrm{\Delta }n)}{dT}=\underset{l}{}\underset{i<j}{}I_lϵ_iϵ_j\frac{d(<S_iS_j>_l)}{dT}.$$
(4)
In general the values of the coefficients $`I_l`$ associated with $`\mathrm{\Delta }n`$ are not related to the respective values $`J_l`$. Nevertheless, we will find a proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ in a temperature regime where the temperature dependence of the contributing correlations $`<S_iS_j>_l`$ are all proportional. This is necessarily the case in the critical region of a phase transition where the correlation length grows larger than any relevant interaction length and all spin-spin correlations necessarily have the same temperature dependence. We may also find an excellent proportionality between $`C_m`$ and $`\mathrm{\Delta }n`$ in the case where one magnetic interaction dominates in both $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$. Such is the case of $`MnF_2`$, where the proportionality holds accurately over a very wide temperature range, $`5<T<100`$ K. A lack of proportionality has been observed in pseudo-low-dimensional systems where neither of the above conditions holds. The percolation threshold presents another opportunity to observe a breakdown in proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$, since small frustrating interactions can affect spins with few or no neighbors connected by the dominant exchange interaction. We will show that this is indeed the case for the percolation threshold sample $`Fe_{0.25}Zn_{0.75}F_2`$. We will also show that the behavior is readily interpreted using the results of MC simulations of this Ising antiferromagnet.

The $`Fe_{0.25}Zn_{0.75}F_2`$ crystal was cut from a boule grown  at the University of California, Santa Barbara. For the adiabatic heat pulse technique, the $`2.24`$ g sample was mounted on a thin sapphire plate using GE7031 varnish. A small Stablohm 800 wire heater, connected to a four-wire constant-power supply, was wound and varnished onto the sapphire plate. A shielded carbon thermometer was attached with varnish to the sample and was connected using a four-wire technique to a current ratio transformer resistance bridge. The sample was suspended inside a sample chamber by $`0.0254`$ mm Be-Cu wires used to connect to the thermometer and heater. While providing good electrical connections, the alloy wires provide only a small heat leak from the sample to the copper sample chamber. In this way, the heat generated by the thermometer leaks away from the sample in a way controlled by a temperature difference, $`\delta T`$, between the sample chamber and the sample. In the absence of a heat pulse, the sample temperature can remain constant if $`\delta T`$ is properly controlled. This is done by controlling $`\delta T`$ with a bridge and controller using the sample thermometer and a thermometer located in a cold finger, commercially calibrated for $`H=0`$, upon which the sample chamber is mounted. A second thermometer in the cold finger is used with a second bridge to determine the absolute temperature. A preliminary calibration is performed to determine $`\delta T`$ versus $`T`$ such that the sample temperature does not drift in the absence of a heat pulse. The calibration is incorporated into computer control algorithms so that the balance is automatically preserved over the entire temperature range during a specific heat experiment. The computer applies a pulse and determines the resulting change in temperature, from which the specific heat is calculated. The thermometry sensitivity is approximately $`50`$ $`\mu `$K. The pulse energy ranged from approximately $`3`$ $`\mu `$J at the lowest temperatures to $`8`$ mJ at the highest temperatures where measurements were made. After collecting specific heat data for the sample and addenda, the sample is removed from the sample chamber and the thermometer is fixed to the sapphire plate. The specific heat of the thermometer, the sapphire plate, the varnish and the wires is then measured at $`H=0`$. This background specific heat is subtracted from all of the specific heat data before further analysis. The phonon contribution to the total specific heat is approximated by the Debye model, valid at temperatures small relative to the Debye temperature. The Debye temperature for pure non-magnetic $`ZnF_2`$ is $`250`$ K and that of $`FeF_2`$ is $`256`$ K. We weighted the two pure Debye temperatures by the respective concentrations to obtain a Debye temperature of $`252`$ K for $`Fe_{0.25}Zn_{0.75}F_2`$. The Debye specific heat is calculated with this Debye temperature and is subtracted from the data leaving the magnetic component of the specific heat ($`C_m`$) which can then be used in comparisons with $`\mathrm{\Delta }n`$ and MC results.

For the birefringence $`(\mathrm{\Delta }n)`$ measurements, two faces were polished parallel to the spin-ordering $`c`$-axis. A linearly polarized laser beam ($`\lambda =632.8`$ nm) impinges normally upon the sample with a polarization $`45^{}`$ to the $`c`$-axis. The beam traverses a distance of $`8.05`$ mm through the sample. A $`0.5`$ mm pin-hole in front of the sample minimizes the effects of concentration gradients and vibrations. The Sénarmont technique is used to measure $`\mathrm{\Delta }n`$ with a resolution of $`2\times 10^9`$. Calibrated carbon resistance thermometers, used for their low field dependence and high sensitivity, yield a temperature stability better than $`50`$ $`\mu `$K. The sample was zero-field-cooled (ZFC) to $`5`$ K before slowly raising the temperature in approximately $`0.1`$ K steps. About $`400`$ s were required to establish equilibrium and to measure $`\mathrm{\Delta }n`$ at each step. Reasonable variations in the rates of heating and cooling and stabilization times had no observable effect on the data.

The correlation functions for the first three nearest-neighbor pairs have been calculated using MC simulations. The magnetic lattice corresponding to the body-centered-tetragonal $`Fe_{1x}Zn_xF_2`$ lattice is described as two cubic sub-lattices of size $`L\times L\times L`$ each, delineated as one dimensional arrays bit coded to accommodate large lattice sizes. All of the results reported here were obtained with $`L=256`$, corresponding to more than $`3.3\times 10^7`$ sites magnetically occupied with probability $`x`$. Beginning at high temperature each magnetic site is randomly visited many times and flipped with a probability given by either the heat bath or metropolis algorithm in temperature steps of $`0.01`$ K. Both periodic boundary conditions and free boundary conditions were applied to the lattices as they were cooled and then warmed in magnetic fields of $`H=0`$ and $`2`$ T. The first three nearest-neighbor interactions were included in the Hamiltonian with values taken from spin-wave dispersion measurements for $`FeF_2`$, where $`J_1=0.069`$ K, $`J_2=5.278`$ K, and $`J_3=0.279`$ K. The correlation functions for these three types of neighbors were calculated at each temperature step and averaged over ten different configurations of magnetic spins. Increasing the number of MC steps per temperature step dramatically did not change the results once a minimum number of steps is performed so that the simulation stays in quasi-static equilibrium at each temperature. The proportionality constant between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ is found to be the same in $`Fe_{0.93}Zn_{0.07}F_2`$ and $`Fe_{0.46}Zn_{0.54}F_2`$ such that, $`AC_m`$ = $`d(\mathrm{\Delta }n)/dT`$, where $`A=9.17\times 10^6`$. Our sample $`Fe_{0.25}Zn_{0.75}F_2`$ is near the percolation threshold concentration and cannot obtain long-range order associated with $`J_2`$. However there is short-range order which produces a rounded peak in the $`C_m`$ data at about $`T=23`$ K. $`d(\mathrm{\Delta }n)/dT`$ vs. $`T`$ data agree with the $`C_m`$ data near this peak when the proportionality constant $`A=9.17\times 10^6`$ is used as shown in Fig. 1. However, at low temperature $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are no longer proportional.

MC simulations provide the opportunity to investigate the relative importance of the different nearest-neighbor interactions on the birefringence. The second-nearest-neighbor interaction dominates the energy Hamiltonian. Therefore the specific heat is dominated by the temperature derivative of the second-nearest-neighbor correlation at the percolation threshold concentration. However the simulations indicate that $`d(\mathrm{\Delta }n)/dT`$ is significantly governed by the temperature derivative of the third nearest-neighbor correlation. This is evident because the two different correlations have contrasting low temperature behaviors, with only the third-nearest-neighbor one changing the sign of its slope at low temperature as shown in Fig. 2.

In Fig. 3 we show the MC results for the Ising $`C_m`$ in the percolation threshold concentration simulation along with the $`d=1`$ Ising and Heisenberg exact models. It is reasonable to expect that the $`d=1`$ Ising model should correspond roughly to the behavior of the simulation. In both systems the spins have an average of two neighbors, there is a similar energy gap and each should have a vanishing specific heat at high temperature. The $`d=1`$ Heisenberg model, on the other hand, has a finite specific heat at zero temperature (the model is not strictly followed by any real system at low temperature). This comparison does help to explain the fact that the MC data, shown again in Fig. 4 as the small amplitude curve, do not correspond very well to the experimental data. Although the real system has a large anisotropy, it nevertheless has a Heisenberg interaction. While the system exhibits an energy gap and its specific heat falls to zero as the temperature is decreased, a significant amount of specific heat near $`T=10`$ K is attributable to this Heisenberg interaction. Since we are interested in the relative contributions to the biregringence from the frustrating third nearest interaction, we choose to multiply the pure Ising behavior of the MC data by a factor of $`1.5`$, thereby obtaining the good fit for $`T<20`$ K exhibited in Fig. 4.

As mentioned above, at magnetic concentrations $`x=0.93`$ and $`x=0.46`$ it has been shown that the proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ is approximately the same as for pure $`FeF_2`$. If we assume the same proportionality for the percolation threshold concentration in the temperature regime where $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are proportional and that the third nearest-neighbor correlation dominates the behavior of the birefringence, we obtain

$$AJ_2\underset{i<j}{}ϵ_iϵ_j\frac{d(<S_iS_j>_2)}{dT}=I_3\underset{i<j}{}ϵ_iϵ_j\frac{d(<S_iS_j>_3)}{dT}.$$
(5)
Furthermore, in this temperature regime we expect

$$\underset{i<j}{}ϵ_iϵ_j<S_iS_j>_2=2\underset{i<j}{}ϵ_iϵ_j<S_iS_j>_3,$$
(6)
where the factor of two arises because there are eight second-nearest neighbors and only four third-nearest neighbors. This allows us to calculate $`I_3`$ in the previous equation to be $`4.85\times 10^5`$. With this value of $`I_3`$ and the third nearest-neighbor correlation from MC simulations we obtain a good fit to $`d(\mathrm{\Delta }n)/dT`$ at low temperature using the same factor of 1.5 used for $`C_m`$, as shown in Fig. 5. We stress that we have calculated $`I_3`$ in a regime where $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are proportional and used this value to fit the simulation results to $`d(\mathrm{\Delta }n)/dT`$, so it is striking that the fit works well below $`T=8`$ K where $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are no longer proportional. This provides strong evidence that the third nearest-neighbor correlation dominates the birefringence at low temperature and in the region where $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ are proportional. However a small contribution from other correlations can not be ruled out. Detailed fits between the correlations calculated by MC simulation and the birefringence below $`8`$ K suggest such contributions are less than ten percent.

To understand why the small, frustrating third-nearest-neighbor interaction affects the third-nearest-neighbor correlation but has little effect on the second-nearest-neighbor correlation, we turn to the computer simulations. When the simulation is run with only the $`J_2`$ interaction included in the Hamiltonian, all three nearest-neighbor correlations have the same temperature dependence as expected. Close inspection of the MC simulation result show that 1.60% of the sites flip at low temperature as a result of the frustration. This is close to the probability of having a site that contains a spin having $`J_3`$ interactions but no $`J_2`$ interactions which is 1.71%. This suggests that it is the spins that are not dominated by the $`J_2`$ interactions that can order via the $`J_3`$ interaction. The ordering of these spins is thereby dominated by the $`J_3`$ interaction and the temperature derivative of the difference between third nearest-neighbor correlations calculated with and without the $`J_3`$ interaction present shows a typical Schottky peak around a temperature comparable to the interaction strength times average number of neighbors. Since the probability of finding such spins decreases dramatically above the percolation threshold concentration, no such effect should be observed at larger concentrations.

We finally turn to the effect of an applied field on the behavior of $`d(\mathrm{\Delta }n)/dT`$ and, hence, the third-nearest-neighbor correlations. Having determined the zero-field proportionality between the birefringence data and simulations of the $`J_3`$ correlations at low temperature, we can now test whether the same proportionality holds in an applied field even though the correlations themselves change significantly with an applied field at low temperatures. In Fig. 5 we show the measured behavior of $`d(\mathrm{\Delta }n)/dT`$ vs. $`T`$ in an applied field $`H=2`$ T along the $`c`$-axis. Also shown in the figure is the MC prediction for the behavior based on the value of $`I_3`$ obtained for zero field. The fit is clearly good. Hence, the value of $`I_3`$ has not changed significantly and the simulation correctly describes the temperature dependence of the birefrincence in the presence of an applied field. This result is important since the birefringence technique has been used extensively to study the RFIM which is realized by applying a field to dilute anisotropic antiferromagnets above the percolation threshold concentration. Near a phase transition the proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ holds in a field and the birefringence technique will yield the proper critical behavior. This is important because the birefringence technique invariably yields higher quality data.

We have demonstrated that a small frustrating interaction, $`J_3`$, affects the corresponding correlation, causing a reversal in the sign of the temperature derivative at temperatures comparable to the interaction strength $`J_3`$. In turn, this causes a reversal in the temperature derivative of the birefringence in the case of $`Fe_{0.25}Zn_{0.75}F_2`$. The behavior of the birefringence in this system is fortuitious because it allows us to determine that third-nearest-neighbor correlation is the dominant one responsible for the birefringence. In most other cases of $`d=3`$ systems, one cannot determine the source of the birefringence because all the correlations have similar temperature dependences. The fact that the proportionality, given by $`A`$, is nearly constant for $`x`$ between the percolation threshold and pure concentrations in $`Fe_xZn_{1x}F_2`$ strongly suggests that the birefringence is dominated by the third-nearest-neighbor correlations for all concentrations. We used the value of $`J_3`$ measured for pure $`FeF_2`$. The good agreement between the MC simulation and the data suggests that the strength of the frustrating interaction $`J_3`$ is close to the same value in the percolation threshold concentration sample. Spin-glass-like behavior near the percolation threshold has strongly hinted at frustration in $`Fe_{0.25}Zn_{0.75}F_2`$, but the strength of the frustration has not been previously determined. Far above the percolation threshold concentration, or at high temperatures near percolation, this frustrating interaction has little affect on the ordering processes or correlations, which all have the same temperature dependence.

By comparing experimental data for $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ and interpreting the results using MC simulation results, we have been able to explain why the proportionality between $`C_m`$ and $`d(\mathrm{\Delta }n)/dT`$ fails at low temperature near the percolation threshold concentration but holds at higher temperatures. From the modeling we have done, we have added to the understanding of the origins of the birefringence in the experimental system. We also have shown that the birefringence technique is well behaved even in a field and, as is well documented, the critical behavior of the RFIM is faithfully characterized with the birefringence technique.

This work was funded by Department of Energy Grant No. DE-FG03-87ER45324.


format: text

text:
# Mean Field Dynamo Saturation: Toward Understanding Conflicting Results

## 1. Field Growth and Constraining the Turbulent EMF

Unlike the turbulent amplification of small scale magnetic energy to near equipartition with the the kinetic energy spectrum, the mean field dynamo (MFD) field generation (c.f. Parker 1979, Kulsrud 1999) on scales $`>`$ turbulent input scale is controversial (c.f. Field et al. 1999). The MFD equation is $`_t\overline{𝐁}=\times 𝐯\times 𝐛+\lambda ^\mathrm{𝟐}\overline{𝐁}+\times (\overline{𝐕}\times \overline{𝐁})`$, with the turbulent EMF $`𝐯\times 𝐛=\alpha \overline{𝐁}\beta \times \overline{𝐁}`$, and pseudo-scalar $`\alpha `$ and scalar $`\beta `$. How well MFD growth applies when the when the dynamic magnetic backreaction is included depends on the survival of $`𝐯\times 𝐛`$. Blackman & Field (2000a) used Ohm’s law and mean field theory (e.g.$`𝐁=𝐛+\overline{𝐁}`$; $`𝐛`$ =0) to constrain the dynamic value of $`𝐯\times 𝐛`$ analytically. Deriving $`𝐯\times 𝐛\overline{𝐁}/𝐜=\eta 𝐣𝐛+𝐞𝐛`$ and then expanding the fluctuating electric field $`𝐞`$ into its potentials, gives $`𝐞𝐛=_𝐭𝐚𝐛/\mathrm{𝟐}𝐜+𝐚\times 𝐞\varphi 𝐛`$. Thus, for $`𝐯\times 𝐛`$ not to be resistively limited, there must be time variation of $`𝐚𝐛`$, or non-vanishing boundary terms. When such terms vanish, helical turbulence without mean field gradients gives $`\alpha (b/\overline{B})^2\alpha _0/R_m`$, where $`\alpha _0`$ is the kinematic value of the psuedoscalar coefficient $`\alpha `$, and $`R_m`$ is the magnetic Reynolds number. There is thus an ambiguity in interpreting all existing numerical experiments suggesting $`\alpha `$ quenching (e.g. Cattaneo & Hughes 1996); the quenching might not be dynamical, but may be due to boundary conditions.

## 2. Magnetic Helicity Escape, Dynamo Action, & Coronal Activity

The above result highlights the role of total magnetic helicity $`H^M=_V𝐀𝐁d^3x`$ (Elsässer 1956), where $`V`$ is a volume of integration, and $`𝐀`$ is the vector potential. That MFD growth involves a magnetic helicity inverse cascade was demonstrated by Pouquet et al. (1976). The $`\alpha `$ effect conserves $`H^M`$ by pumping a positive (negative) amount to scales $`>L`$ (the outer turbulent scale) and a negative (positive) amount to scales $`L`$. Brandenburg’s (2000) simulations confirm this inverse cascade and the role of $`H^M`$ conservation.

A large-scale field can be generated only as fast as $`H_M`$ can be removed or dissipated. Presently, simulations have invoked boundary conditions for which the growth of large scale field is resisitvely limited. Large $`R_m`$ systems must rely on open boundary conditions. To see this, note that $`H^M`$ satisfies $`_t(𝐀𝐁)+c(𝐄\times 𝐀+A_0𝐁)=2c𝐄𝐁`$ where $`𝐄=\frac{𝐕}{c}\times 𝐁`$. Consider two cases. Case (1): The mean scale = universal scale, or the integration is over periodic boundaries. Then boundary terms vanish, so $`_t𝐀𝐁=2c𝐄𝐁=2c\overline{𝐄}\overline{𝐁}2c𝐞𝐛=\eta 𝐉𝐁`$ and $`_t(\overline{𝐀}\overline{𝐁})=2c\overline{𝐄}\overline{𝐁};_t𝐚𝐛=2c𝐞𝐛`$. This is the case of section 1. Dynamo action is resistively limited. Case (2): The system (e.g. Galaxy or Sun) mean volume V $`<<`$ universal volume. Here we must use the gauge invariant relative helicity $`H_R^M`$ inside and outside of the spherical or disk rotator (Berger & Field 1984). The integral over the the universal volume then satisfies $`_t_U𝐀𝐁d^3x=_tH_{R,in}^M+_tH_{R,out}^M=2c_U𝐄𝐁d^3x0`$. The formulae for the $`H_R^M`$ of the mean and fluctuating quantities inside the rotator are $`_tH_{R,in}(\overline{𝐁})=2c_{in}\overline{𝐄}\overline{𝐁}d^3x+2c_{S_{in}}(\overline{𝐀}_p\times \overline{𝐄})𝑑𝐒`$ and $`_tH_{R,in}(𝐛)=2c_{in}𝐞𝐛d^3x+2c_{S_{in}}𝐚_p\times 𝐞𝑑𝐒.`$ In a steady state, $`_tH_{R,in}=0=_tH_{R,out}`$. This and $`𝐄𝐁0`$ imply that the above surface terms above must be equal and opposite. Moreover, the surface term balances the $`\overline{𝐄}\overline{𝐁}`$ term in the $`_tH_{R,in}^M`$ equation. The boundary term can thus allow for a significant turbulent EMF because the latter is contained in $`\overline{𝐄}\overline{𝐁}`$. Dynamo action unrestricted by resisitivity is possible only in case (2). This is consistent with Pouquet et al. (1976) and Brandenburg (2000).

If $`H_M`$ flows through the boundary, then so does magnetic magnetic energy. Blackman & Field (2000b) showed that a typical minimum power leaving the system when a MFD is operating is given by $`\dot{E}^M\frac{k_{\mathrm{min}}}{8\pi }|\dot{H}^M|=\frac{k_{\mathrm{min}}}{6}\left|\alpha \overline{B}^2\right|V`$, where $`V`$ is the system volume. Dynamos operating in the Sun, accretion disks, and the Galaxy would then lead to a net escape of magnetic energy and small and large scale magnetic helicity. Coronal activity from the emergence and dissipation of helical magnetic flux is thus a prediction of the MFD in all of these cases, and is observed directly in the Sun (c.f. Pevtsov et al. 1999). For the Galaxy, $`\dot{E}^M>(\pi R^2)\alpha \overline{B}^210^{40}(R/12\mathrm{k}\mathrm{p}\mathrm{c})^2(\alpha /10^5\mathrm{cm}/\mathrm{s})(\overline{B}/5\times 10^6\mathrm{G})^2\mathrm{erg}/\mathrm{s}`$, in each hemisphere. Blackman & Field (2000b) discuss how this relation may be consistent with coronal energy input rates required by Savage (1995) and Reynolds et al. (1999).

## 3. Open Questions

An MFD unlimited by resitivivity requires the helicity to flow through the boundary AND that there be some mechanism that enables this flow. Thus there are two separate issues. Even if the boundary conditions allow it, does it actually happen? One may have to include the dynamics of buoyancy or winds to fully demonstrate the non-resistive MFD. Note that turbulent diffusion of the mean magnetic field (not necessarily the actual field) across the boundary is required to maintain a quadrupole field in the Galaxy with a net flux inside the disk. Similarly, for the Sun, the solar cycle requires net diffusion through the boundary. The flow of helicity would appeal to the same dynamics needed by these constraints.

The analytic and numerical studies that we have seen which show catastrophic suppression of the dynamo coefficients, or resistively limited dynamo action, either (1) invoke periodic boundary conditions, and/or (2) are 2-D, or (3) do not distinguish between zeroth order isotropic components of the turbulence and the higher order anistropic perturbations for a weak mean field (Blackman & Field 1999). This means that there always seems to be an alternative explanation. and the observed suppression is then ambiguous as we have described. This does not mean that some of the physical concepts found in the strong suppression results are invalid, but just that they may be valid only for the restricted cases considered. For example, the observation that the Lagrangian chaos properties of the flow are changed in the presence of a weak mean field for turbulence in a periodic box (e.g. Cattaneo et al. 1996) needs to be understood in the relation to the imposed boundary conditions, and the shape of the magnetic energy spectrum (e.g. dominated at small or large scales?). Along these lines, note that the helicity constraint is global, but also becomes a constraint for any sub-volume of a periodic box once the system is fully mixed.

## 4. References

Blackman, E. G. & Field, G. B. 1999, ApJ 521 597.

Blackman, E. G. & Field, G. B. 2000a, ApJ 534 984.

Blackman, E. G. & Field, G. B. 2000b, MNRAS, in press, astro-ph/9912459.

Berger M.C. & Field G.B., 1984, JFM, 147 133.

Brandenburg A., 2000, submitted to ApJ, astro-ph/0006186.

Cattaneo, F., & Hughes, D.W. 1996, Phys. Rev. E. 54, 4532.

Cattaneo, F., & Hughes, D.W., Kim E.-J, 1996, PRL 76, 2057.

Elsässer, W. M. 1956, Rev. Mod. Phys. 23, 135.

Field, G. 1986, in Magnetospheric Phenomena in Astrophysics, R. Epstein & W. Feldman, eds., AIP Conf. Proc. 144 (Los Alamos: Los Alamos Sci. Lab.), 324.

Field, G. B., Blackman, E. G., & Chou, H. 1999, ApJ 513, 638.

Frisch, U., Pouquet, A., Léorat, J. & Mazure, A. 1975, JFM 68, 769.

Krause F. & Rädler K.-H., 1980, Mean-field magnetohydrodynamics and dynamo theory, (New York: Pergamon).

Kulsrud R., 1999, ARAA, 37 37.

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

Pevtsov, A., Canfield, R. C., & Brown, M. R. eds. 1999, Magnetic Helicity in Space and Laboratory Plasmas, (American Geophysysical Union).

Pouquet, A., Frisch, U., & Leorat, J. 1976, JFM 77, 321.

Reynolds R.J., Haffner L.M., Tufte S.L., 1999, ApJ 525 21.

Savage B.D., 1995, in The Physics of the Interstellar Medium and Intergalactic Medium, A. Ferrara, C.F. McKee, C. Heiles, & P.R. Shapiro eds. ASP conf ser vol 60. (San Francisco: PASP) p233.


text:
# Solitons in the one-dimensional forest fire model

\[

## Abstract

Fires in the one-dimensional Bak-Chen-Tang forest fire model propagate as solitons, resembling shocks in Burgers turbulence. The branching of solitons, creating new fires, is balanced by the pair-wise annihilation of oppositely moving solitons. Two distinct, diverging length scales appear in the limit where the growth rate of trees, $`p`$, vanishes. The width of the solitons, $`w`$, diverges as a power law, $`1/p`$, while the average distance between solitons diverges much faster as $`d\mathrm{exp}(\pi ^2/12p)`$.

PACS numbers: 05.65.+b, 05.45.Df

\]

The Bak-Chen-Tang (BCT) forest fire model is a simplified model of turbulent phenomena. Propagating fires dissipate or burn trees at an average rate determined by the tree growth probability, $`p`$, representing power fed into the system. The fires are self-sustaining, with no spontaneous ignition.

The fires exhibit non-trivial spatio-temporal correlations in the slow driving limit, where the growth rate of trees, $`p`$, vanishes . In two and three dimensions, a “scale-dependent” critical behavior emerges , which is not adequately understood at present. The fractal dimension of fires depends on the length scale of observation, up to the correlation length, $`\xi `$, where the fires become space filling. This correlation length diverges with the inverse growth rate $`1/p`$ to a power $`2/d`$ for $`d=2,3`$ . The behavior in the one dimensional ($`d=1`$) BCT model, where the fires are constrained by dimensionality to have a much simpler structure, has not previously been investigated.

Here we show that in one dimension, as $`p0`$, fires self-organize into a state characterized by propagating, branching and annihilating solitons. These are localized objects, compared to their scale of separation, moving with constant velocity. Outside of the solitons only solid forest exists. As an advancing fire moves into solid forest, each advancing fire is trailed by a cloud of smaller forests separated by holes, with possibly a few fires. Each fire advancing into solid forest together with its trail of finite forests and fires constitutes a soliton. New solitons are created by a process where forward moving solitons emit opposite moving ones.

The width of the solitons, $`w`$, diverges as $`1/p`$. Within the solitons, the finite forests, which are connected segments of trees, have a distribution of sizes, $`P(s)1/s^2`$ up to the cut-off scale $`1/p`$ imposed by the width of the solitons. The other statistical properties are completely determined by the branching rate for solitons, which can also be calculated analytically. This branching process creates new solitons, and must be balanced by the collision rate of oppositely moving solitons, which annihilates both solitons. We find that the density of fires remains finite but decreases extremely fast, $`n\mathrm{exp}(\pi ^2/12p)`$, in the limit $`p0`$. Thus, even in one dimension, the fires are self-sustaining with no spontaneous ignition in the slow driving limit. Although the solitons are diverging in width, their width vanishes compared to their separation, and they are well-defined in the limit $`p0`$.

All of this behavior is entirely different from the behavior of the Drossel-Schwabl (DS) forest fire model, which self-organizes into a critical state with a power law distribution of forest fires. (The BCT model does not have a power-law distribution of forest fires.) In the DS model, fires are not self-sustaining but are injected at a small rate, $`f`$. The DS model is critical in the limit $`f/p0`$ where a separation in time scales between burning of entire forests and the growth of trees occurs. Since the forests burn down instantly, and the fires are injected rather than self-sustaining, this limit gives a completely different physical picture than the propagating fires in the strict $`f=0`$ problem, corresponding to the original BCT model. Solitons, obviously, cannot appear in the DS model.

The BCT forest fire model is defined as follows. Each site on a $`d`$ dimensional lattice of linear extent $`L`$, can be in one of three states: occupied by a fire, occupied by a tree, or empty. During a time step, all sites which contain fire burn down, leaving an empty site, and ignite neighboring trees on the lattice. Then all empty sites are independently occupied with trees with probability $`p`$. This two step process is repeated indefinitely. The model can be studied with either open or periodic boundary conditions, and the results discussed below are found in either case. After a transient period, the system enters a statistically stationary state with a complex distribution of fires, and forests, which are connected clusters of trees. The fire is self-sustaining and fed by the tree growth process, for $`p>1/\mathrm{ln}L`$ (for $`d=1`$). Otherwise the fire dies out completely.

\[

\]

In order to get a clear picture, it is useful to view the space-time dynamics of the burning process. The numerical simulations were initiated with a random distribution of trees and fires. The spatio-temporal behavior in the steady-state is shown in Figure 1 for $`p=0.15`$. The fires propagate either left or right, moving into solid forests, with unit velocity. Each such fire leaves a trail of forests, empty sites, and possibly some fires in their wake. Each moving trail, led by a fire advancing with unit velocity into a fully connected forest, constitutes a soliton. After some time, of order $`1/p`$, the small forests trailing in the wake of the advancing fire have grown to form a completely connected forest again, with no trace of the passing fire. Thus the width of the soliton scales as $`1/p`$. For small $`p`$, such as that shown in Figure 1, the solitons are well separated from each other and collisions are rare.

When fires moving in opposite directions collide, they annihilate each other, leaving clusters of small forests. After a time interval of order $`1/p`$, the small forests have all healed and joined the large surrounding forest, and the remnants of both solitons have completely disappeared.

Now and then, solitons emit solitons propagating in the opposite direction. This is the fundamental process for the self-organizing dynamics. Consider a state with a density $`n`$ of random placed solitons, each with a randomly chosen direction. In the stationary state, the rate of soliton death due to collisions must balance the rate of new solitons emitted as back-fires. The expected life-time of a soliton is equal to the average distance between solitons, $`d=1/n`$. Assume that each soliton emits back-propagating solitons at a rate $`r`$. Stationarity requires $`r=1/d=n`$.

In order for a backward-moving fire to survive, and emerge as a soliton, the new fire must survive a number of time steps, by using the new, growing forests inside the “parent” soliton as stepping stones, before it can escape from its parent soliton and propagate in a fully-connected forest. Since the system is one-dimensional, one might suspect that this probability will vanish, for any $`p<1`$, since sooner or later the newly created fire would always meet an empty site inside its parent soliton. Actually, for any finite $`p`$ this is not the case.

At the first step, where the initial fire moves forward one unit into a fully connected forest, there is a probability $`p`$ that this fire will branch backwards. In this case, there will be two fires at the same time at neighboring sites. This means that a site that was on fire in the previous time step, $`t1`$, is regrown with a tree and set again on fire at time $`t`$. Thus $`p`$ is the basic branching probability for new fires to be created.

We will consider the process where the newly emitted fire only moves backward (see Fig. 1), opposite from the original fire, with no forward motion, or wandering. This is the dominant process for $`p0`$. What is the probability that such a newly created fire will survive one additional time step? It is $`(1(1p)^3)`$, since it is now three time steps since the forward-moving, parent fire passed the advancing position of the new, daughter fire. At the $`m`$’th time step for the new fire, the conditional probability of surviving one more time step is $`1(1p)^{2m+1}`$. Thus, the probability for surviving $`m`$ steps is

$$P_{surv}(m)=p\underset{m^{}=1}{\overset{m}{}}(1(1p)^{2m^{}+1}).$$
(1)
Taking the limit $`m\mathrm{}`$ and $`p0`$ gives

$$\underset{p0}{lim}\mathrm{ln}P_{surv}=\frac{1}{2p}_0^{\mathrm{}}\mathrm{ln}(1e^u)𝑑u=\frac{\pi ^2}{12p}.$$
(2)
This is the rate of emitting back-propagating solitons. From the previous argument $`P_{surv}=n`$. This result has been compared with numerical simulations as shown in Fig. 2. Numerical simulation results for $`0.1<p<0.3`$, and soliton distances ranging from $`20`$ to $`10000`$, yield an excellent fit to the form $`n=\mathrm{exp}(C/p)`$, with $`C=0.834`$, within 2% of the exact value $`\pi ^2/12`$.

Note that there are two length scales, each diverging with $`1/p`$, but in very different ways. One is the width of solitons, $`w1/p`$. The other is the much larger distance between solitons, $`d=\mathrm{exp}(\pi ^2/12p)`$. Because of the exponentially growing distance between solitons, it is numerically feasible to simulate the model only for $`p>0.1`$. However, the length scales are already extremely well separated for these values, establishing the soliton picture.

Also for any finite $`p`$ there is always a finite probability that a soliton will annihilate itself by encountering a rare hole in the nominally fully connected forest. However the rate of this annihilation process for a soliton is $`e^{apd}`$, where $`a`$ is a constant of order unity. This is derived by considering that such a rare hole left over was created by a soliton approximately $`d`$ time steps in the past. This rate vanishes, in the limit $`p0`$, compared to the rate $`1/d`$ for the pair-wise annihilation described above.

The distribution of the forest sizes within the comparatively narrow solitons has a structure determined entirely by a cascade process in which isolated clusters of trees form bigger and bigger forests . This is caused by tree growth filling in the empty sites between neighboring forests. Forests of size $`s`$ are both created and destroyed by this same process. We obtain analytically a power-law distribution of forests within solitons, $`P(s)1/s^2`$, up to a cut-off, which is the size of the soliton (of order $`1/p`$). The same cascade process occurs in the Nagel-Schreckenberg traffic model , where it leads to $`1/f`$ noise  and also in the Drossel-Schwabl forest fire model , which has spontaneous ignition of fires and instantaneous burning of forests  (and no solitons), unlike the model discussed here. It is very important to distinguish these small forests, forming the micro-structure of the solitons, from the exponentially larger forests separating solitons. In the limit $`p0`$, the solitons are point like, viewed at the large scale of their separation, and fractal-like, viewed at scales less than their width.

In one dimension, a soliton can propagate only if the forest density in front of it is unity, so the solitons do not feel the effects of previous solitons, and, except for branching, there is no interaction other than at contact. This behavior is somewhat reminiscent of one-dimensional Burgers turbulence, where the dynamics is also governed by colliding solitons, or shocks . Persistent, localized propagating structures that interact via collisions have also been observed in some cellular automata, which unlike the model discussed here, are deterministic . One distinct feature of the BCT forest fire model is that a single soliton can emit other solitons, unlike shocks in Burgers turbulence, or solitons in cellular automata.

In higher dimensions, the propagating fires interact with each other at a distance. This interaction is mediated by their effect on the forest density and tree-tree correlations. The resulting spatio-temporal structure of dissipation is much more complicated. Thus it is not clear if the soliton description will remain valid or useful for the higher dimensional BCT models, or other systems where the interactions between dissipating structures are long-range. Nevertheless, the simple results found here for the one-dimensional BCT forest fire model suggests that solitons, or shocks, may govern the dynamics of other self-organized critical models of turbulent, intermittent phenomena. In fact, we have observed similar behavior in an interface depinning model . In this case the sites on fire correspond to the sites which are moving, and the limit $`p0`$ corresponds to the limit where the driving force approaches the depinning threshold from the moving phase.

PB and MP thank the National University of Singapore for great hospitality, and working environment.


text:
# ARCS, The Arcminute Radio Cluster-lens Search I. Selection Criteria and Initial Results

## 1 Introduction

Observations of gravitational lensing provide a probe of the distribution of mass in the universe on scales ranging from micro-lensing by sub-solar mass objects in the Galactic halo to weak-lensing by large-scale structures of the universe. In order to survey matter both light and dark on different size scales, unbiased systematic searches for lensing events within a well-defined sample of background sources must be carried out. Since the discovery of giant arcs in the region of clusters A370 and Cl2244 (Lynds & Petrosian 1986; Soucail et al. 1987) more than 50 clusters have been found to contain arc-like images of background galaxies, including a clear example of multiple imaging by the cluster 0024+1654 (Colley, Tyson & Turner 1996). The lensed images provide diagnostics for the cluster potential, allowing the distribution of mass to be mapped independently of the distribution of light. However, the currently known examples cannot tell us about the general prevalence of mass structures on the scale of clusters of galaxies (CGs) because they are in no sense statistically complete, having largely been discovered through observations of previously known clusters. Only by carefully examining a large well-defined sample of background sources for cases of gravitational lensing can we hope to place constraints on the general distribution of cluster-sized masses.

$`N`$-body cold-dark-matter (CDM) simulations of the universe make predictions about the degree of clustering of dark matter haloes at different epochs (for example Jenkins et al. 1998). The statistics of strong gravitational lensing events within a sample of background sources should therefore allow us to distinguish between models with different rates of cluster evolution. Early simulations by Wambsganss et al (1995) suggested that there could be as many strongly-lensed sources with image separations $`>`$5″, associated with groups and clusters of galaxies, as there are with image separations $`<`$5″ associated with single galaxy lenses. The JVAS and CLASS lens surveys (King et al. 1999; Myers 1996) found a lensing rate of approximately 1 in 500 for lens systems with image separations between 0$`\stackrel{}{.}`$3 and 6″. This suggests that in a search for lens systems associated with clusters of galaxies, a background source population of at least $`10^3`$ is necessary before useful cosmological constraints can be set.

Previous searches for strong lensing events with image separations greater than a few arcseconds have started with quasar surveys because they provide a sample of objects with high mean redshift, which increases the lensing probability, and relatively low surface density, which reduces the number of false positives from random associations. Maoz et al. (1997) used a subsample of 76 quasars from the HST snapshot survey (Maoz et al. 1992) to look for lensing events with angular separations of 7″ to 50″. The fact that no lens systems were found allowed some extreme cosmological scenarios to be ruled out, for example that there exists a large population of ‘failed’ or ‘dark’ clusters. However, the small size of the parent sample limited the ability to place useful constraints on more conventional cosmological models.

In designing surveys to look for cluster-mass lenses, a basic problem arises when trying to recognise strongly-lensed images. As the image separation increases, the time delay between the images increases, and at a separation of tens of arcsec the characteristic delay can be hundreds of years. Basic properties of quasars such as the continuum emission, the emission line spectrum and the radio source nuclear structure may change significantly over this time. Due to their intrinsic variability, quasars are therefore not optimally-suited for use as probes for lensing by clusters of galaxies. A search method using background sources which are expected, a priori, to show no significant changes in physical properties over timescales of $``$1000 years is required to overcome this problem.

This paper is concerned with the Arcminute Radio Cluster-lens Search (ARCS), a project aimed at looking for strong gravitational lensing by galaxy clusters in a systematic and uniform manner. We use the fact that extended emission ($`>`$few kpc) from an extragalactic radio source is unlikely to vary over the timescales of interest. In order to maximise our parent sample of background sources, whilst minimising the observational effort in following up potential lens candidates, we started from the existing FIRST 1.4 GHz survey (Becker, White & Helfand 1995) and the Automatic Plate Measuring (APM) optical catalogue (for example Irwin, Maddox & McMahon 1994). The FIRST survey provides both high sensitivity and sufficient resolution to allow us to reject most multiple sources as candidate lens systems directly before any follow-up observations have taken place. Furthermore, demanding stellar optical identifications of the radio sources from the APM catalogue both increases the mean redshift of the parent sample and provides a filter against selecting intrinsic double sources as lens candidates.

We describe our selection algorithm for identifying gravitational lens candidates from the FIRST and APM catalogues in section 2. We then go on to describe the subsequent follow-up observations in section 3, including an analysis of the remaining CG-lens candidates at each observational stage. A review of the project and future work is laid out in section 4.

## 2 Lens Candidate Selection

Our aim is to maximise the probability of detecting gravitational multiple images produced by CGs while minimising the observational effort required to confirm or refute them. As we pointed out above, reliable searches for lensing events on these scales require us to select background parent sources which should not vary over the relevant timescales which could extend to hundreds of years. We must also maximise the redshift of these background sources which increases the lensing probability. The FIRST survey combined with the APM catalogue derived from the Palomar Observatory Sky Survey (POSS) plates enables us to achieve these goals.

### 2.1 The FIRST Survey

The Faint Images of the Radio Sky at Twenty centimetres (FIRST) survey will eventually cover 10 000 square degrees of the north galactic cap. The survey area will coincide with the Sloan Digital Sky Survey (SDSS), upon completion of which $``$50 per cent of all FIRST sources will be optically identified to the $`m(v)24`$ limit of SDSS. For our purposes, the north Galactic cap portion of the catalogue dated 98feb04 was used. This contains 382 892 sources covering 4150 square degrees of the sky in the region RA$`7^\mathrm{h}20^\mathrm{m}`$ to RA$`17^\mathrm{h}20^\mathrm{m}`$, and $`\delta `$22°12′to $`\delta `$57°36′. The lower flux density limit is 1 mJy with a typical rms noise of 0.15 mJy and a beam-size of 5$`\stackrel{}{.}`$4; source positions are tied to the JVAS calibrator reference frame (Patnaik et al. 1992; Browne et al. 1998; Wilkinson et al. 1998) with a systematic uncertainty of $`<`$0$`\stackrel{}{.}`$03, while the overall systematic positional errors are $`<`$0$`\stackrel{}{.}`$1 throughout the survey area. Source parameters were derived by fitting elliptical Gaussian models. A more complete description of the FIRST catalogue can be found at URL http://sundog.stsci.edu/.

Approximately 2.5 per cent of ‘sources’ in the FIRST catalogue are flagged as possible side-lobes of strong sources. Only about 15 per cent of these sources are likely to be real, equivalent to less than 0.5 per cent of the total. We therefore removed all possible side-lobes from our version of the catalogue prior to our analysis.

### 2.2 The APM Catalogue

The APM catalogue is derived from the UK Schmidt and POSS plates. For our present purposes we used APM data obtained from the POSS for the region covered by FIRST. The APM catalogue contains data from blue O plates and red E plates to limiting magnitudes of 21.5 and 20 respectively, with a positional accuracy of $``$1″. Each detection is parameterised in terms of position, magnitude, and a classification into a stellar or a non-stellar object based on the light profiles of a detection; a ‘merged’ classification exists where an object falls between these two types. A more complete description of the APM sky catalogues can be found at URL http://www.ast.cam.ac.uk/`~`apmcat/.

### 2.3 Selection of Primary Sources

In order to look for possible cases of multiple imaging we must define a parent sample of “primary” sources around which we will search for “secondary” images.

With an average beam size of $``$5″, the FIRST survey can confidently identify resolved structure down to an angular scale of $``$2$`\stackrel{}{.}`$5 (FWHM). For redshifts $`z0.2`$, this corresponds to a linear scale of $`7`$ kpc (for $`H_0=65`$ km s-1 Mpc-1). A region of this size will not vary in spectral or luminosity characteristics over periods of hundreds of years. We therefore placed a minimum angular size limit of 2$`\stackrel{}{.}`$5 along the deconvolved major axis of our primary sources. We chose the maximum angular size of our primary sources to be 15″ to avoid confusion with potentially-lensed secondary images with similar angular separations.

With the resolution of the FIRST survey many faint radio sources will be resolved into multiple components with angular separations less than 15″; these components are listed separately in the FIRST catalogue. We therefore constructed two classifications of primary sources, those with a single component and those with two components (see Figure 1). One could include primary sources with three or more components, but for simplicity we ignored such sources (inclusion of these sources would only increase our sample size by $``$10 per cent). In the case of a single component our maximum primary source size of 15″refers to the major axis given in the FIRST catalogue, in the case of two components it refers to the angular separation of the components. We defined the central position of a two-component primary source to lie on the line connecting the individual component positions, weighted by the integrated component flux densities (see for example Macklin 1981). The “fractional centroid distance” being the fractional distance to the centroid from the brighter component is defined as $`\beta =F_B/(F_A+F_B)`$ where $`F_A`$, $`F_B`$ are the integrated flux densities of the brighter and dimmer components respectively. We can then write the angular distance from the brighter component to the weighted central position as $`\theta \times \beta `$, where $`\theta `$ is the angular separation of the two components. Two extreme cases of composite primary sources are illustrated in Figure 2.

In order to select “clean” primary fields and to avoid double counting, we placed a window with a radius of 15″centred on the central position (defined above) of each primary source. We only included a source in our primary sample if no additional components were within the annulus. This will help prevent more complex (i.e. $`3+`$ component) sources from infiltrating our sample. We finally imposed the condition that the integrated flux density must be $``$35 mJy.

To increase the probability of lensing, it is desirable to select background sources with as high a mean redshift as possible. To this end we selected only those primary radio sources with stellar optical counterparts in the APM catalogue on the O and/or E plate. By choosing stellar objects with associated radio emission we expect to select a high proportion of quasars, and hence to maximise the mean redshift (see section 2.4.3). Radio sources with non-stellar i.e. galaxy identifications will tend to have redshifts $`0.2`$ and therefore have a greatly reduced probability of being lensed. Optical identifications with a limiting E magnitude brighter than 20 also make it possible to characterise the redshift distribution of the primary sample with optical spectroscopy using moderate-sized telescopes. In addition, an optical identification on a primary source will provide a filter against selecting classical double radio sources as lens candidates (see section 2.5).

A significant number of faint red galaxies are identified as non-stellar on the E plates, but have a stellar classification on the O plates due to the detection of only the central region of the galaxy. Figure 3 illustrates the three distinct populations present in a primary source sample when we allowed there to be a non-stellar identification on one of the POSS plates. To simplify the optical identification process for ARCS, an optical object with different classifications on the two plates is not acceptable as a quasar identification. In fact, as Figure 3 shows, those sources with stellar classifications on one plate and stellar or no classification on the other plate are predominantly blue, as would be expected for quasars (see for example the catalogue derived from four SDSS filters described by Newberg et al. 1999).

To take into account the positional accuracies of the APM catalogue, a 3″ difference between the APM and radio positions was tolerated (see section 2.4.2). Radio positions for the single component sources were taken directly from the FIRST catalogue, whereas the composite source positions were as defined above. After applying all these selection criteria, a final parent sample of 1023 primary sources remained (see Table 1).

### 2.4 Statistics of the Primary Source Sample

#### 2.4.1 Random Radio Associations Posing as Two-component Primaries

By selecting two-component primary sources from a catalogue with the depth of FIRST we will inevitably introduce a contamination to the primary sample in the form of pairs of sources which have no physical association. We expect this number to be small, but we have made a rough estimate of the number of chance associations of unrelated radio sources posing as two-component sources in the FIRST catalogue. Our selection routine found 5264 two-component primary sources from the radio selection criteria alone i.e. prior to the application of the optical criteria. From the known surface density of FIRST sources the number of random associations was estimated to be $``$100. This is a contamination of only $``$2 per cent amongst the two-component sources before our optical criteria have been applied (i.e. $`<`$1 per cent of all primaries – see Table 1). Within our final primary sample we expect this percentage to have been reduced significantly because randomly associated pairs of sources will be less likely to have stellar counterparts at their centroid position in comparison to a genuine physical pair. Our conclusion is that the level of contamination to the primary sample from the selection of false two-component primary sources will be negligible.

#### 2.4.2 Random Radio – Optical Associations

Since the error in the radio astrometry is negligible relative to the optical astrometric accuracy, and we assume that the optical positional error is purely gaussian, then the separation between the optical position and the true (i.e. radio) position, $`r`$, is described by a Rayleigh distribution, viz.

$$P(r)dr=\frac{r}{\sigma ^2}\text{exp}\left[\frac{r^2}{2\sigma ^2}\right]dr$$

where $`\sigma `$ is both the mean and the standard deviation of the distribution. However, there is a small but significant contamination within our primary sample of radio sources with chance associations with, for example, galactic stellar sources. This modifies the number of optical ‘identifications’ by

$$N_{ran}(r)dr=\mathrm{\Sigma }_{}N_{rad}2\pi rdr$$

where $`\mathrm{\Sigma }_{}`$ is the stellar surface density arcsec-2 in the APM catalogue from both POSS plates, $`N_{rad}`$ is the number of primary sources before optical identification, and $`r`$ is the distance from the radio source. Given the fact that $`\mathrm{\Sigma }_{}5\times 10^4`$ arcsec-2 within the FIRST region, the total number of random associations amongst the 1023 primary sources is estimated to be $``$180 to a radius of 3″. Figure 4 shows a histogram of the angular separation of optical and radio positions in our primary source sample. The distribution is well described by a Rayleigh-plus-random distribution, increasing our confidence that we are selecting extragalactic radio sources with quasar identifications. The best fit value of $`\sigma =`$0$`\stackrel{}{.}`$86 agrees well with the quoted error in the APM positions. However, to select as many primary sources as possible we retained our 3″ tolerance for FIRST–APM positional discrepancies. In an analysis of the lensing frequency within the sample, we must modify the primary sample size to $``$850 sources to take account of the expected number of random associations.

#### 2.4.3 Mean Redshift

We searched the NASA/IPAC Extragalactic Database (NED) for additional information on our primary sources. A total of 133 sources were found to have a measured redshift at the FIRST position. We calculated values of $`\overline{z}1.09`$, $`\overline{R}17.5`$ and $`\overline{B}18.1`$ for this sub-sample. The total primary sample is on average fainter with $`\overline{R}18.6`$ and $`\overline{B}19.5`$. Whilst the sub-sample obtained from NED will be unquantifiably contaminated by a multitude of selection effects, it is reasonable to assume that the total primary sample will have a higher value for $`\overline{z}`$. In fact, an initial analysis of spectroscopic data from a representative but limited sub-sample of the primary sample suggests that $`\overline{z}1.4`$ (Phillips et al., in preparation).

### 2.5 Initial Lens Candidate Selection

Having selected the primary sources the resolution and sensitivity of the FIRST survey allowed us to search for potential multiple imaging events directly. We looked in the catalogue near each primary source for any secondary sources within an angular separation range of 15″ and 60″ (Figure 5); any such secondary sources are putative lensed images of the primary sources. We did not impose any size selection criteria to the secondaries because a secondary lensed image, with a lower magnification than the primary, will not necessarily be resolved by the FIRST beam. We also omitted optical selection criteria; since the secondary will often be weaker, its optical emission could place it below the limit of the POSS/APM catalogue; in addition a stellar object may have varied in magnitude. However, we did impose a maximum flux ratio of 5:1 between the primary and the secondary. This places a lower limit on the integrated flux density of a secondary source in the weakest systems to 7 mJy. If a secondary source is indeed a lensed image of the primary, the conservation of surface brightness means that the secondary will not be resolved out i.e. it will be detectable in FIRST.

In section 2.3 we mentioned that a large proportion of extragalactic radio sources will be resolved into multiple components by the FIRST survey. Many such sources consist of two ‘lobes’ of emission, the so-called “classical doubles”. It is possible that in selecting primary sources with nearby secondaries, we are selecting a high proportion of classical double radio sources as lens candidates. However, the optical identification criterion applied during the primary source selection process provides a strong filter against these types of radio source posing as lens candidates. This is because the optical counterpart to such an object would be expected to lie between the lobes rather than on one or both of them as our selection criterion demands.

These selection criteria produced only 91 initial lens candidates i.e. over 90 per cent of the primary sample could be immediately rejected using the FIRST and POSS/APM information.

### 2.6 Further Rejection of Lens Candidates

Up until this point the lens candidate selection process was conducted automatically. We then went on to examine the FIRST and APM data individually by eye.

The FIRST radio images of the 91 initial lens candidates along with APM information for each field were combined using the AIPS (Astronomical Image Processing Software) package. Examples of typical fields are shown in Figure 6 and the top of Figures 7, 8 and 9. It is immediately clear that we can reject some more candidates by eye. The arguments for rejecting a system as a lens candidate can be broken down into two categories: (1) those where the primary and secondary sources are clearly components of a single larger source; and (2) those where lensing arguments can be invoked.
(1) Single source arguments:
The single source rejection arguments rely on the nature of known multiple-component radio sources. The morphologies of such sources can include specific features which would not be expected if the components were lensed images of each other. Examples are: extended emission bridging the region between the primary and secondary source, and ‘outer’ edges of a primary-secondary pair being brighter than the ‘inner’ edges. We also reject lens candidates where we find an optical counterpart between the primary-secondary pair coincident with an additional unresolved radio component; this is the position where we would expect an AGN to be present. Examples of sources rejected as single sources on the basis of the FIRST maps are 0713+369 and 1358+290 shown in Figure 6.
(2) Lens arguments:
We use the basic fact that the surface brightness of the background source is conserved in each lensed image. However, we can only apply this argument to the extended radio structure. The image flux ratio in the radio and optical regimes can be different due to the the variability arguments laid out in section 2.5.

From the original 91 sources with secondary components, 53 could be rejected following this individual examination leaving only 38 lens candidates which required further observations i.e. less than 4 per cent of the primary sample.

## 3 Radio Follow-up of Candidates

In order to test the lens hypothesis for the remaining candidates we made higher resolution observations to study candidate source morphology and surface brightness in more detail.

### 3.1 Stage 1 – VLA Observations

The sample of 38 lens candidates was observed with the VLA in A configuration at 6 cm and 20 cm in two 12-hour periods on 1 and 7 April 1998. The observing bandwidth was 50 MHz in each of two IF channels for both wavelengths, with observing periods of approximately 8 minutes per source per wavelength. The rms sensitivity was $`0.05`$ mJy at both wavelengths. The minimum beam sizes obtained were 0$`\stackrel{}{.}`$4 at 6 cm and 1$`\stackrel{}{.}`$4 at 20 cm when sources were observed near the zenith. The data were calibrated in a standard fashion using the AIPS package. Flux calibration was based on two 2-minute observations of 3C286 assuming a flux density of 7.360 Jy at 6 cm, and 14.554 Jy at 20 cm based on the expression derived by Baars et al. (1977). Polarisation calibration was based on three 1.5-minute observations of 1613+342, which was assumed to be unpolarised to correct for antenna effects, while 3C286 was used for phase offset corrections between the two measured orthogonal polarisation systems. Each source was CLEANed at both 6 cm and 20 cm for Stokes parameter I using the DIFMAP software package (Shepherd 1997) and where appropriate self-calibration was applied. Polarisation and spectral index data reduction was performed with the standard AIPS routines. In these cases a uniform convolving beam size of 1$`\stackrel{}{.}`$5 was used at both 6 cm and 20 cm for ease of image comparison. A Gaussian taper was applied in the uv plane to the 6 cm data to 30 per cent at $`200`$ k$`\lambda `$ in order to restrict the uv range to that of the 20 cm data.

We used similar arguments to those described in section 2.6 to classify each system. The 38 lens candidates fall into three distinct groups: single ‘classical double’ radio sources ($``$68.5 per cent); unrelated small angular separation pairs ($``$26.5 per cent); and 2 remaining lens candidates ($``$5 per cent). Figures 7,8 and 9 show examples from each of these groups.
Examples of single sources (Figure 7):
In both examples, 1248+393 and 1519+345, it can be seen from the A configuration 20 cm data that the primary and secondary components are outer edge-brightened. The 6 cm data show that the spectra of the outer edges are flattening, consistent with these regions being hot-spots embedded within diffuse lobes.
Examples of two unrelated double sources (Figure 8):
The 6 cm A configuration map of example 1000+417 shows a compact double primary (northerly) source with a central optical counterpart. The secondary source consists of a central unresolved flat-spectrum component and extended components to the north-east and south-west, with extended emission bridging the entire source. We can rule this example out as a lens candidate using surface brightness arguments.

In example 1034+379 the 20 cm A configuration data shows extended emission bridging two components in the primary (northerly) source, while the 6 cm data reveals an unresolved flat-spectrum component at the position of the optical counterpart. The secondary source at 20 cm contains two edge-brightened components with a non-stellar optical counterpart between, in a position that a host galaxy would be expected. This example can be ruled out as a lens candidate based on the dissimilar morphologies of the primary and secondary sources, and the orientations of each source relative to one another.
Remaining lens candidates (Figure 9):
The A configuration map of 1053+457 at 20 cm shows two compact double component sources in a plausible lens configuration. The 6 cm data show that the primary (western) source in fact consists of extended emission bridging two unresolved components, whereas the secondary (eastern) source has a single unresolved component with an extended component to the south-east. Despite the dissimilarities between sources, we cannot unequivocally rule out this system as a lens candidate since unresolved radio components can vary significantly over the time-scale of a time delay. The 1519+387 primary (western) source is a composite consisting of an extended western component and an unresolved eastern component in the A configuration 20 cm data. The secondary source appears largely unresolved in the 6 cm data with a possible detection of extended emission towards the unresolved component of the primary source. While these components could be interpreted as a single source with a central unresolved core at P1 (see Figure 9), high resolution MERLIN data are needed to confirm the extension in the secondary source. At this stage it remains as a lens candidate.

### 3.2 Stage 2 – MERLIN Observations

High resolution data on the two remaining lens candidates 1053+457 and 1519+387 were taken with MERLIN at 6 cm in two 16-hour long-track observations on 12 and 21 December 1998. The data were in 16$`\times `$1 MHz bands with an integration time of 8 seconds, giving an rms sensitivity of $`0.05`$ mJy; the minimum beam size was $``$0$`\stackrel{}{.}`$04. Flux calibration was based on 25 minute observations of 3C286 assuming a flux density of 7.09 Jy. Initial editing of the data was carried out with the standard MERLIN programs with calibration and mapping again performed in AIPS and DIFMAP.

The unresolved components in the VLA 6 cm data of the primary source of 1053+457 are both clearly resolved with MERLIN (Figure 10) and appear to be edge-brightened. Additionally a central unresolved component is revealed between the two extended components and is close to an optical object (see Figure 9). The unresolved component in the secondary source also appears to be slightly resolved while the extended component to the south east is almost resolved out. On the lens hypothesis we would expect the two resolved components in the primary source also to be detected in the secondary source. This is not the case. We detect only one resolved component in the secondary, which is also of lower surface brightness than either of the resolved components in the primary source. We can therefore reject 1053+457 as a lens candidate. The secondary source in 1519+387 is clearly resolved with MERLIN while the primary component remains unresolved. We can categorically reject this source as a lens candidate.

## 4 Discussion and Summary

We have completed what is currently the largest unbiased survey for multiple image gravitational lens systems with image separations $`>`$15″. Our primary sample of 1023 sources was selected from the FIRST and APM catalogues, and our search for lens candidates within this sample was achieved by looking for additional nearby images within the FIRST catalogue. As a result our follow-up required only 38 candidates be observed with the VLA, and just two of these necessitated higher resolution observations with MERLIN. No cases of gravitational multiple imaging with image separations in the range 15″to 60″were found. A synopsis of the numbers of sources at each stage of the selection process is given in Table 1.

Wambsganss et al (1995) made predictions of “wide-separation” lensing rates using $`N`$-body simulations with SCDM parameters (see also Wambsganss, Cen & Ostriker, 1998). However, the resolution and scale of the simulations were not adequate to provide a reliable estimate of the number of multiple imaging events that we would expect to see in an unbiased search. In particular, cluster mass density profiles must be simulated more accurately because the shape of the profile has a significant influence on a cluster’s ability to produce multiple images. In order to draw any significant cosmological conclusions from our result more realistic predictions of the number of wide-separation multiple imaging events in different cosmological scenarios need to be made. The results of such simulations are also necessary if we are to determine the sample size for any future expansion of the ARCS project.

To examine the practicability of increasing the size of the ARCS sample we have calculated the maximum size of a sample which can be selected from currently available FIRST data and the pending SDSS data. As described in section 2.1 FIRST has been designed to coincide with the SDSS survey region. It is estimated that $``$50 per cent of all sources in FIRST will be identifiable on SDSS to $`m(v)`$24 down to the 1 mJy limit of FIRST. Using our radio-only primary source sample of 12 932 (Table 1) this could result in a primary sample of $`>`$6500 if SDSS is used instead of the APM catalogue at the optical identification stage. Taking a more conservative approach by specifying an optical magnitude limit of $`m(v)`$23 we would still expect to increase our primary sample to $`>`$3000. It should be noted that these figures are firm lower limits on what could be achieved because our radio-only primary sample has a lower 1.4 GHz flux limit of 35 mJy (c.f. the 1 mJy limit of FIRST). The latest version of the FIRST catalogue dated 99jul21 contains 477 247 reliable sources over 5450 square degrees of the northern sky (c.f. 382 892 sources in the 98feb04 version we used). If we also lower the integrated flux density limit of a primary source to 5 mJy, thus giving a minimum secondary flux density at the detection limit of the FIRST survey i.e. 1 mJy, a total radio-only primary sample of 62 891 sources results. Applying the above optical arguments with a limiting magnitude of $`m(v)`$23, this should produce a primary sample of $`>`$15 000. It would appear that upon completion of SDSS there will be no problem constructing a sample of background sources at least an order of magnitude larger than that presented in this paper.

It is practicable to search such a primary source sample for instances of wide-separation lensing events. With the experience gained in this pilot study, we believe that the efficiency of our search technique could be improved by up to a factor $``$2. This is largely based upon our increased confidence in the FIRST survey images for use in rejecting false positives before any further observations of lens candidates are made at higher resolution. Selecting a primary sample of, say, 10 000 sources, a two-fold improvement in the initial candidate selection would result in only $``$200 candidates requiring VLA A configuration follow-up with $``$10 requiring further observations with MERLIN. This is equivalent to $``$3 days of VLA time and $``$3 days of MERLIN time.

In summary, we have established that it is possible to make a reliable unbiased search for lens systems with image separations in the range 15″to 60″. Furthermore it is possible to use existing radio and optical data to significantly reduce the required observing time for such a search. What we currently lack are reliable predictions of the lensing frequency for different cosmologies to compare with our null result.

## Acknowledgments

We wish to thank the VLA and MERLIN staff for making these observations possible. The VLA is the Very Large Array and is operated by the National Radio Astronomy Observatory which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. MERLIN is the Multi-Element Radio Linked Interferometer Network and is a national facility operated by the University of Manchester on behalf of PPARC. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research was supported in part by the European Commission, TMR Programme, Research Network Contract ERBFMRXCT96-0034 ‘CERES’. PMP would also like to thank PPARC for the support of a studentship award.