text:
# Halo and Relic Sources in Clusters of Galaxies

## 1 Introduction

Diffuse radio sources in clusters of galaxies are the probe of the existence of non-thermal emission from the cluster intergalactic medium (IGM). These sources have no optical identification, steep radio spectrum and size of $``$ 1 Mpc. Diffuse sources with regular shape, permeating the cluster center are called radio halos, the prototype being the radio halo Coma C in the Coma cluster (Giovannini et al. 1993). Other similar sources, with irregular morphologies, located in cluster peripheral regions are named relics. In a few clusters, both a central halo and a peripheral relic are present.

The difficulty in explaining these diffuse sources arises from the combination of their large size and the short lifetimes of radiating electrons. Several theoretical models have been proposed: (i) in-situ acceleration by turbulent gas motion or shocks produced in the intergalactic medium during the merger of a subgroup into the main cluster, (ii) diffusion of relativistic electrons out of the present head-tail radio galaxies, (iii) secondary particle production by hadronic interaction of relativistic protons with the background gas of the IGM.

The knowledge of the physical properties of these sources, and of their origin and evolution is limited by the low number of halos and relics well studied up to now. To search for new sources of this class and to test the possible correlation between the presence of tailed radio galaxies and diffuse sources suggested by Giovannini et al. (1993), we have obtained observations with the Very Large Array (VLA 111The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under contract with the National Science Foundation.) of 21 clusters with the following selection criteria: a) 10 clusters where the possible presence of a diffuse source was suggested in the literature, but the angular resolution and sensitivity of the available images was not good enough to confirm its existence and derive the source parameters; b) 11 clusters showing at least a tailed radio galaxy within 300 kpc (0.1 Abell radius) from the cluster center in the sample of O’Dea and Owen (1985).

To properly map these extended sources one needs large sensitivity to the extended features, but also with the high resolution necessary to distinguish a real halo from the blend of unrelated sources. Observations were carried out using different VLA configurations at 20 cm for more distant clusters whereas the nearest clusters were observed at 90 cm. We present here sensitive images of the detected diffuse sources. For the computation of intrinsic parameters, a Hubble constant H0=50 km s-1 Mpc-1 and a deceleration parameter q0=0.5 are assumed.

## 2 Observations and data reduction

The list of observed cluster of galaxies is given in Table 1 where the observing frequency and the VLA configurations are also given.

The bandwidths in the 1.4 GHz observations are of 25 MHz for the B array, and of 50 MHz for the C and D arrays. At 0.3 GHz, the bandwith is of 3.125 MHz.

The observing frequency and array were chosen in order to have a good sampling of short spacings, ensuring that a halo-type source similar to Coma C could be easily detected and imaged. We note that Coma C is completely resolved out from the VLA at 20 cm also in its more compact array (D), because of missing short spacings. For the most distant clusters in the present sample, the VLA in D configuration at 20 cm was found to be suitable, whereas for the nearest sources observations at 90 cm were necessary. The signal to noise ratio achieved by our observations should be enough to map halo sources with a surface brightness 3-5 times lower than Coma C. Moreover, each cluster was observed with higher angular resolution to allow the separation of the discrete sources from the diffuse cluster emission.

The data were calibrated and reduced with the Astronomical Image Processing System (AIPS), following the standard procedure: Fourier-Transform, Clean and Restore. Self-calibration was applied to minimise the effects of amplitude and phase fluctuations. At each frequency, the data of different arrays were combined to obtain the final maps. Images were produced at different angular resolutions to study both the extended structures and the discrete sources. In the caption of each figure, we give the Half Power Beam Width (HPBW) and the noise level of the maps. A comparison with available maps from the NRAO VLA Sky Survey (NVSS, Condon et al. 1998) and the Faint Images of the Radio Sky at Twenty-cm (FIRST, Becker et al. 1995) Survey was also made.

## 3 Results

Among the 10 clusters where the presence of a diffuse sources was suggested by previous authors, we have detected diffuse sources in 7 clusters. The parameters of these sources are given in Table 2. We have not detected halos or relics in the sample of clusters selected because of the presence of tailed radio galaxies. In Table 3 we give the surface brightness limit for the undetected extended sources and the largest size detectable in our observations. In fact, we cannot exclude the presence of a source with size larger than the limit, because of the missing short spacings. Notes on the clusters hosting diffuse sources, and on some other interesting clusters are given in the following. We present images of all the diffuse sources, as well as a few maps of extended and/or peculiar radio galaxies detected in this study.

### 3.1 Notes on individual clusters of galaxies with diffuse sources

A 85. This cluster was studied in detail by Bagchi et al. (1998). In Fig. 1 we present the image of the central cluster region obtained at 90 cm, with angular resolution of 25′′. A small size radio emission is detected from the central cD galaxy (A), and two head-tail radio galaxies are visible (B,C).

Moreover, two extended regions of radio emission with no obvious optical identification are present: they are labelled ’R’ and ’D’ in Fig. 1. Source ’R’ is the well known relic source of steep spectrum studied by Joshi et al. (1986) and Bagchi et al. (1998). In this image, the relic has an irregular structure. An image with slightly better angular resolution (see Fig. 2) shows that the relic has a complex structure with evidence of internal substructures and filaments. This is confirmed by the high resolution image presented by Andernach et al. (2000). The total flux density (2.74 Jy) is lower than the total flux density measured with the Ooty Synthesis Radio Telescope at the same frequency (3.15 Jy, Joshi et al. 1986). We refer to Bagchi et al. (1998) for the total radio spectrum of this source. We obtained a spectral index map between 90 and 20 cm, combining our image with the 20 cm image retrieved from the NVSS, presented in Giovannini et al. 1999.) Because of the low resolution of the NVSS image (45′′), the substructures are smoothed in the maps used for the spectral comparison. The spectrum in very steep ($`\alpha `$2.5 - 3.0) with no evidence of substructure.

The sources ‘B’ and ’D’ are blended together in the 90 cm image by Bagchi et al. (1998), where the extended feature formed by these two sources is suggested to be another possible diffuse source. From the present high resolution image, we distinguish the tailed radio galaxy ’B’ and remark the existence of the extended feature labelled ’D’, which has a steep spectrum ($`\alpha _{0.3}^{1.4}22.5`$, obtained by comparison with the NVSS) and a total extent of about 2, corresponding to $``$175 kpc. Its shape and the presence of a nearby galaxy with a possible faint radio emission at about the 3$`\sigma `$ level suggests that this structure could be a dying tailed radio galaxy, where the nuclear emission has almost completely ceased.

A 610. We detect a faint diffuse emission (Fig. 3), coincident with the halo source detected by Valentijn et al. (1979). It is elongated in shape and displaced from the cluster center, thus we classify it as a relic source. The high resolution image given in Fig. 4 confirms that the source is diffuse. Using the fux at 0.6 GHz from Valentijn et al. (1979), S$`_{0.6GHz}`$ = 59 mJy, a spectral index $`\alpha _{0.6}^{1.4}`$ = 1.4 is derived.

A 665. Moffet & Birkinshaw (1989) detected a diffuse emission coincident with the cluster center. Our image (Fig. 5) shows that this halo source is very extended and elongated in the SE-NW direction. The highest resolution map (not shown here) indicates that 2 faint discrete sources are present within the diffuse radio emission (crosses in Fig. 5) and that the radio emission is asymmetric with respect to the cluster center, being brighter and more extended toward NW.

A 2142. The presence of a halo in this cluster was suggested by Harris et al. (1977). We confirm the presence of diffuse emission, located around the brightest cluster galaxy (Fig. 6), as also detected in the NVSS (Giovannini et al. 1999). The radio emission could be more extended in S - SW direction up to 27 10 in declination (2 sigma level), but more sensitive observations are necessary to confirm the reality of this feature. In the FIRST survey, the diffuse emission is completely resolved out, and no emission is detected at a flux level of 0.4 mJy/beam (3$`\sigma `$) from any of the the bright galaxies present in the region of the extended radio emission. The radio source could therefore be classified as a cluster halo. We note, however, that this source is much smaller than radio halos commonly found in clusters. We cannot rule out the possibility that it is the remnant of a single galaxy which was active in the past.

A 2218. We confirm the existence of the diffuse source found by Moffet & Birkinshaw (1989). It is slightly displaced from the cluster center, toward East (Fig. 7). In a deep image obtained by Zwaan et al. (2000) the diffuse source is more extended, with radio emission of very low surface brightness also to the West of the cluster center. This halo source is smaller than typical cluster-wide radio halos. It is similar, both in structure and in size, to the halo in A 1300 (Reid et al. 1999). As in the case of A 1300, the diffuse radio source is located in the same direction as the extension of the cluster X-ray emission detected by the ROSAT HRI (Markevitch 1997).

At 5 GHz, the flux density of the halo is reported to be 0.6 mJy (Partridge et al. 1987). Comparison of this flux with the present measurement at 1.4 GHz provides a spectral index $`\alpha _5^{1.4}`$ $``$ 1.6, which should however be taken cautiously since the two maps at the two frequencies do not have matched Fourier coverage.

CL 0016+16. We confirm the presence of a diffuse emission near the cluster center found by Moffet & Birkinshaw (1989). Unfortunaly, our high resolution observations (VLA configurations B and C) are strongly affected by interferences therefore we cannot produce a higher resolution map. The diffuse radio emission (Fig. 8) is located at the cluster center, and is extended more than 1 Mpc, therefore it is a typical radio halo. We note the presence of a fainter radio emission west of the halo, at around RA2000 = 00h 18m 25s, DEC2000 = 16 26 48′′, coincident with a galaxy group visible in the DPSS (see Fig. 8).

0917+75. This source was studied in detail by Dewdney et al. (1991) and Harris et al. (1993). Our image (Fig. 9) shows that the diffuse source is elongated in the E-W direction, with flux density and morphology in very good agreement with their data.

As discussed in the above mentioned papers, this source cannot be obviously associated with a well defined cluster of galaxies. In the region, there are the 3 clusters A 786, A 787 and A 762, at redshift larger than 0.1, belonging to the Rood Group of clusters of galaxies #27, but the source is located at very large distance from their centers (respectively 30, 45, and 44, corresponding to $``$ 5 Mpc, $``$ 8 Mpc and $``$ 8 Mpc). Thus, although the radio emission is reminiscent of the giant radio relic in A 3667 (Röttgering et al. 1997), its very large distance from the cluster center would favour the hypothesis that this diffuse source is related to a possible group near the source center (see the optical DPSS image in Fig. 9) or to the supercluster. This makes this source quite unique.

### 3.2 Notes on individual clusters without diffuse sources

A 84. An extended head-tail radio galaxy is present in the cluster center region, as reported by Rudnick & Owen (1976) and O’Dea & Owen (1985). In our map, the tailed radio galaxies is detected to a much larger extent than previously believed, and for the first time a prominent bend of the tail is visible (Fig. 10).

Moreover, a faint extended radio emission is detected at about 2 SE of the tailed radio galaxy. This feature has no obvious optical counterpart, and is completely resolved out in a higher resolution map. Its structure and the presence of a galaxy at the boundary of the radio emission might support the hypothesis that it is the remnant of a tailed radio galaxy, where the nuclear emission has ceased, but it could also be a distant unrelated source with no optical counterpart.

A 119. This cluster was studied in detail by Feretti et al. (1999). The present data at 90 cm confirm the absence of diffuse emission.

A401 \- Our image at 90 cm shows the same features detected by Harris et al. (1980) and Roland et al. (1981), in particular the two components A and B (see Fig. 11). On the contrary, the faint extended emission detected near the central cD galaxy in the NVSS (see Giovannini et al. 1999) at 20 cm is not seen here. With the present data, the existence of diffuse emission in this cluster remains quite dubious, as well as the classification of source B as a cluster halo.

A 629. No diffuse emission is detected in this cluster, which is characterized by the presence of an extended wide-angle-tail radio galaxy and of a head-tail radio galaxy at the center (see Owen & Ledlow 1997).

A 754. The extended emission at the cluster center, reported by Harris et al. (1980) and also visible in the NVSS image (Giovannini et al. 1999) is at the limit of significance in our map at 90 cm. This implies that this feature cannot have a steep spectral index between 90 and 20 cm ($`\alpha _{90}^{20}`$ $`<`$ 0.7) and confirms the suggestion of Giovannini et al. (1999) that the radio emission might be the blend of many bright cluster galaxies in this region, rather than a diffuse halo. The diffuse emission detected in the NVSS in the eastern peripheral region (Giovannini et al. 1999) cannot be seen in the 90 cm image because of the too low angular resolution.

A 1609. The radio halo suggested by Valentijn (1979) has not been detected in our image, which shows several radio sources. The extended radio galaxy with a tailed structure (Fig. 12) is actually the blend of two tailed radio galaxies (Rudnick, private communication).

A 1775. The head-tail radio galaxy in this cluster (B1339+266B) shows a very narrow tail (Fig. 13), extended $``$ 8 (850 kpc), i.e. twice more than in the image of Owen & Ledlow (1997). A similar head-tail radio galaxy, but with a shorter extent ($``$550 kpc) is found in the cluster A 2256.

A 2250. The head-tail radio galaxy near the cluster center (B1709+397) is much more extended than in the previous map by O’Dea and Owen (1985) ($``$ 6, corresponding to $``$ 600 kpc), but no diffuse halo is visible (Fig. 14).

## 4 Discussion and conclusions

The present paper increases the number of halo and relic sources with detailed radio information. The new data on the two clusters A 665 and CL 0016+16 are of great importance in the study of cluster-wide radio halos, since they represent two more examples of giant sources of this class. Both the halos in A 665 and in CL0016+16 are more powerful than the Coma halo.

The cluster A 665 has a bolometric X-ray luminosity of 4.2 1045 erg s-1 (Wu et al. 1999), a temperature of 8.3 keV (Wu et al. 1999) and no cooling flow (White et al. 1997). From X-ray data, it is suggested to be in a postmerger state, with an asymmetric X-ray brightness distribution, extended in the NW direction (Markevitch 1996, Gomez et al. 2000). This is also the orientation of the radio structure (Fig. 15 and Sect. 3), indicating a spatial correlation between the radio and X-ray structure, as found in other clusters (see e.g. Feretti 2000).

The cluster CL 0016+16 has a bolometric X-ray luminosity of 2.81 1045 erg s-1 (Wu et al. 1999). Its X-ray emission has been studied by Neumann & Böhringer (1997) who obtain a temperature of 8.2 keV, conclude that the cluster has not yet formed a steady cooling flow, and report strong evidence of substructure which is possibly due to merging of a galaxy group with the main cluster. The radio-X-ray overlay, presented in Fig. 16, shows that the radio halo is located at the cluster center, at the position of the X-ray peak.

In the clusters A 2142 and A 2218, diffuse halos of about 500 kpc in size, or less, are detected. Both these clusters are characterized by signatures of recent merger events (Markevitch 1997, Markevitch et al. 1998). No cooling flow is detected in A 2218, whereas the cluster A 2142 has a massive cooling flow (Peres et al. 1998) centered on the diffuse radio source. The presence of a radio halo in a cooling flow cluster is quite uncommon (Feretti 2000). The fact that the halo in A 2142 has a small size could be related with the presence of the cooling flow. Future deep observations of more clusters of galaxies are necessary to understand if these sources have a different nature with respect to the giant cluster-wide halos, or if they are the low power-small size tail of the same class of radio sources.

We have detected no radio halos in the 11 clusters showing tailed radio galaxies at the cluster center. In the Coma cluster, Giovannini et al. (1993) suggested that the tailed radio galaxy NGC4869, which is orbiting around the cluster center, could be responsible for the supply of relativistic electron radiating within the radio halo Coma C. Also in the clusters A 2255 (Feretti et al. 1997a) and A 2319 (Feretti et al. 1997b), tailed radio galaxies are found within the diffuse radio emission, or close to it, supporting the hypothesis that the tailed radio galaxies could provide relativistic electrons to the radio halo. In the sample of tailed radio galaxies presented by O’Dea and Owen (1985), there are 15 clusters with at least a tailed radio galaxy within 300 kpc from the cluster center (Giovannini et al. 1993). These are the 11 clusters studied here, plus A 85, A 1656 (Coma), A 2255 and A 2256. The last 3 clusters show a halo source, while A 85 hosts a relic radio emission (Sect. 3). For the 11 clusters observed here, we only give a limit in size and brightness for the presence of a possible extended source. The lack of detection of bright halos in these clusters would imply that the connection between halos and tailed radio galaxies seems not relevant. On the other hand, we cannot exclude that these tailed radio galaxies are not really at the cluster center but simply projected onto it, or that they are crossing the cluster center and not orbiting at the cluster center, as needed to fuel a diffuse source. From rotation measure studies of the sources in A 119 (Feretti et al. 1999), where no halo is detected, it was derived that the two powerful tailed radio galaxies are genuinely located at the cluster center. We note, in addition, that no tailed radio galaxy is found to be present in A 665, hosting a cluster-wide radio halo. We conclude that the tailed radio galaxies could be a source of relativistic electrons, but they are not a necessary ingredient for the halo formation. These results are in line with the model developed by Brunetti et al. (2000a), who suggest that relativistic particles radiating within the halos were injected in the cluster volume by past major merger events, starbust and AGN activity. Following a recent suggestion by Brunetti et al. (2000b), the relativistic particles deposited by the tailed radio galaxies could account for the emission of extreme ultraviolet radiation.

According to Bliton et al. (1998), the clusters containing narrow-angle tailed radio galaxies are in general dynamically complex systems udergoing merger events. Actually, structure and other merger signatures have been detected in all the clusters of Table 3 (e.g. A 754, Henry & Briel 1995; A 119, Markevitch et al. 1998) which have been adequately observed in optical or X-ray, except for A 401, which is a rather unusual cD cluster with neither a merger activity nor a cooling flow. Therefore, we confirm that the presence of a recent merger is not sufficient for the formation of a radio halo.

Recent studies have shown that radio halos are found in merging clusters which are also characterized by high X-ray luminosity, i.e. high temperature and large mass (Giovannini et al. 1999, Feretti 2000). The halos in A 665 and CL 0016+16 are fully confirming the above result. The connection between halos and mergers, and the presence of halos in hot, massive, X-ray bright clusters can be interpreted in the framework of the two-phase model suggested by Brunetti et al. (2000a) and successfully applied to Coma C. According to this model, relativistic particles produced in the early cluster history, are reaccelerated by the energy supplied from shocks and turbulence in a recent merger event. Alternative models and the theoretical implications of non-thermal emission from the intracluster medium are reviewed by Enßlin (2000).

We thank the referee, Dr. T. Enßlin, for helpful comments. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with the National Aeronautics and Space Administration. We acknowledge partial financial support from the Italian Space Agency (ASI).
References

Abell G.O., Corwin H.G. Jr., Olowin R.P., 1989, ApJS 70, 1

Andernach H., et al., 2000, The Universe at low Radio Frequencies, Symposium IAU 199, Pune (India), in press

Bagchi J., Pislar V., Lima Neto G.B, 1998, MNRAS 296, L23

Becker, R.H., White, R.L., Helfand, D.J., 1995, ApJ 450, 559

Bliton M., Rizza E., Burns J.O., Owen F.N., Ledlow M.J., 1998, MNRAS 301, 609

Brunetti G., Setti G., Feretti L., Giovannini G., 2000a, MNRAS Submitted

Brunetti G., Setti G., Feretti L., Giovannini G., 2000b, in preparation

Condon J.J., Cotton W.D., Greisen E.W. Yin Q.F., Perley R.A., Taylor G.B., Broderick J.J., 1998, AJ 115, 1693

Dewdney P.E., Costain C.H., McHardy I., Willis A.G., Harris D.E., Stern C.P., 1991, ApJS 76, 1055

Enßlin T.A. 2000, In: The Universe at low Radio Frequencies, Symposium IAU 199, Pune (India), in press, astro-ph/0001433

Feretti L., Böhringer H., Giovannini G., Neumann D., 1997a, A&A 317, 432

Feretti L., Giovannini G. Böhringer H., 1997b, New Astronomy 2, 501

Feretti L., Dallacasa D., Govoni F. Giovannini G., Taylor G.B., Klein U., 1999, A&A 344, 472

Feretti L., 2000, In: The Universe at low Radio Frequencies, Symposium IAU 199, Pune (India), in press, astro-ph/0006379

Giovannini G., Feretti L., Venturi T., Kim K.-T., Kronberg P.P., 1993, ApJ, 406, 399

Giovannini G., Tordi M., Feretti L., 1999, New Astronomy 4, 141

Gomez P.L., Hughes J.P., Birkinshaw M, 2000, ApJ, in press, astro-ph/0004263

Harris D.E., Bahcall N.A., Strom R.G., 1977, A&A 60, 27

Harris D.E., Kapahi V.K., Ekers R.D., 1980, A&AS 39, 215

Harris D.E., Stern C.P., Willis A.G., Dewdney P.E., 1993, AJ 105, 769

Henry J.P., Briel U.G., 1995, ApJ 443, L9

Joshi M.N., Bagchi J., Kapahi V.K., 1986, In: Radio continuum processes in clusters of galaxies, NRAO, p. 73

Markevitch M., 1996, ApJ 465, L1

Markevitch M., 1997, ApJ 483, L17

Markevitch M., Forman W.R., Sarazin C.L., Vikhlinin A., 1998, ApJ 503, 77

Moffet A.T., Birkinshaw M., 1989, AJ 98, 1148

Neumann, D.M., Böhringer H., 1997, MNRAS 289, 123

O’Dea C.P., Owen F.N., 1985, AJ 90, 927

Owen F.N., Ledlow M.J., 1997, ApJS 108, 41

Partridge R. B., Perley R. A., Mandolesi N., Delpino F., 1987, ApJ 317, 112

Peres C.B., Fabian A.C., Edge A.C., Allen S.W., Johnstone R.M., White D.A., 1998, MNRAS 298, 416

Reid A.D., Hunstead R.W., Lémonon L., Pierre M.M., 1999, MNRAS 302, 571

Roland, J. Sol, H. Pauliny-Toth, I. Witzel, A. 1981, A&A 100, 7

Röttgering H.J.A., Wieringa M.H., Hunstead R.W., Ekers R.D., 1997, MNRAS 290, 577

Rudnick L., Owen F.N., 1976, ApJ 203, L107

Valentijn E.A., 1979, A&AS 38, 319

White D.A., Jones C., Forman W., 1997, MNRAS 292, 419

Wu Xiang-Ping, Xue Yan-Jie, Fang Li-Zhi, 1999, ApJ 524, 22

Zwaan M.A., van Dokkum P.G., Verheijen M.A.W., Briggs F., 2000, poster presented to the conference Gas and Galaxy evolution, Socorro, NRAO, in press


format: text

text:
# Attractive Forces Between Global Monopoles and Domain Walls

## I Introduction

In a previous paper in collaboration with Andrew Sornborger, we investigated the interaction of global topologically stable monopoles with embedded domain walls in a $`SO(3)\times Z_2`$ linear sigma model . We were motivated by Dvali et al.’s proposal  that monopoles could be swept up by domain walls which, if true, provides a potential solution for the monopole problem. In our previous analysis the $`Z_2`$ symmetry factor corresponding to the domain wall was a “semi-discrete” symmetry, i.e. part of the connected symmetry group. Hence, the walls corresponding to this $`Z_2`$ were not stable; they were embedded walls . However, we confirmed that both monopoles and anti-monopoles are attracted to the domain wall and that after colliding with the wall they unwind before the wall itself decays. Furthermore, we confirmed Dvali et al.’s conjecture that the monopole charge dissipates on the surface of the wall by studying the evolution of both the winding number $`Q_w`$, and the winding number density $`\rho _w`$.

The question immediately prompted by our previous results is whether the same situation will be found when global monopoles interact with topologically stable domain walls. Recent work  provides numerical evidence that the attraction-unwinding mechanism extends to models with both topologically stable local monopoles and domain walls. In this paper, we ask whether or not the attraction-unwinding mechanism carries over to models with stable global monopoles and $`Z_2`$ walls.

In our previous work , we argued that the attractive force between the monopole and the wall would also be found in systems with general global and gauge monopoles and domain walls because the underlying physics is generic, and rests primarily on the minimization of the gradient energies of the field configurations. In this paper we present evidence supporting this conjecture by studying a $`G=O(3)\times Z_2`$ linear sigma model spontaneously broken to $`H=O(2)`$. The vacuum manifold is $`=𝒢/=\mathrm{S}^2\times \mathrm{Z}_2`$ and thus admits both topologically stable domain walls and monopoles. Our numerical simulations confirm the existence of the attractive force between monopole and domain wall. However, in contrast to the interaction between global monopoles and embedded global walls, we find that the monopole fails to unwind on the surface of the domain wall. Finally, we provide an analytic explanation of the results seen in the numerical simulations.

## II The Model

Topologically stable monopoles are ensured  when the second homotopy group of the vacuum manifold is nontrivial, i.e. $`\mathrm{\Pi }_2()0`$. Similarly, stable domain walls require the existence of a nontrivial $`\mathrm{\Pi }_0()`$. Therefore, to get both stable domain walls and monopoles one needs $`Z_2`$ and $`S^2`$ factors to be present in the vacuum manifold. In this investigation we consider a field theory with a novel potential which admits just such a vacuum manifold. The symmetry breaking pattern is

$$G=O(3)\times Z_2H=O(2).$$
(1)
Therefore, the vacuum manifold is

$$=\frac{G}{H}=\frac{O(3)\times Z_2}{SO(2)}S^2\times Z_2.$$
(2)
To realize a symmetry breaking pattern which yields both stable monopoles and domain walls, we use  the following Lagrangian:

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}_\mu \mathrm{\Phi }^a^\mu \mathrm{\Phi }^a{\displaystyle \frac{\lambda _1}{4}}\varphi _4^2(\varphi _1^2+\varphi _2^2+\varphi _3^2\eta _1^2)^2`$ (4)
$`{\displaystyle \frac{\lambda _2}{4}}(\varphi _4^2\eta _2^2)^2`$
where $`\mathrm{\Phi }^a=(\varphi _1,\varphi _2,\varphi _3,\varphi _4)`$ is a four component scalar field.

This model has an exact $`Z_2`$ symmetry corresponding to the discrete transformation $`\varphi _4\varphi _4`$. Furthermore, the Lagrangian is invariant under an $`O(3)`$ rotation of the fields $`\varphi _1,\varphi _2,\varphi _3`$. Hence the Lagrangian has an $`O(3)\times Z_2`$ symmetry. The set of vacuum states is

$`\varphi _1^2+\varphi _2^2+\varphi _3^2`$ $`=`$ $`\eta _1^2`$ (5)
$`\varphi _4^2`$ $`=`$ $`\eta _2^2`$ (6)
and the vacuum manifold hence is $`S^2\times Z_2`$.

Since

$$\mathrm{\Pi }_2()=\mathrm{\Pi }_2(G/H)=Z\times Z_21,$$
(7)
the existence of stable monopoles is ensured. A monopole configuration is

$`\varphi `$ $`=`$ $`(\varphi _1,\varphi _2,\varphi _3)=f(r,t)({\displaystyle \frac{x}{r}},{\displaystyle \frac{y}{r}},{\displaystyle \frac{z}{r}})\eta _1`$ (8)
$`\varphi _4`$ $`=`$ $`\pm \eta _2`$ (9)
where $`r^2=x^2+y^2+z^2`$ and where the profile function $`f(r,t)`$ obeys $`f(r,t)=f(r)0`$ as $`r0`$ and $`f(r,t)=f(r)1`$ as $`r\mathrm{}`$. A reasonable ansatz for the profile function at the initial time $`t_0`$ is

$$f(r,t_o)=(1e^{\frac{r}{r_c}}),$$
(10)
where $`r_c`$ is the core radius of the monopole which is determined by balancing gradient and potential energies, which yields $`r_c\lambda ^{1/2}\eta _1^1`$.

Likewise, since

$$\mathrm{\Pi }_0()0,$$
(11)
our Lagrangian also supports the existence of stable domain walls. A domain wall is constructed by picking one spatial direction (e.g. the $`z`$ direction), and setting up a kink configuration in the $`\varphi _4`$ field:

$$\varphi _4(z)=\eta _2\mathrm{tanh}\left(\frac{z}{w}\right),$$
(12)
and letting the other three fields be at a fixed point in the $`S^2`$ subspace of $``$. Here, $`w`$ is the thickness of the wall, $`\sqrt{2/\lambda _2}/\eta _2`$.

This model has several advantages compared to the model with three scalar fields studied in . In particular, it admits stable domain walls rather than unstable embedded walls. Moreover, it is much easier to set up a configuration corresponding to a monopole displaced a certain distance $`z_0`$ from a domain wall as one does not have to rely on a product ansatz to simulate the monopole / domain wall configuration. The $`\varphi _4`$ field independently generates a pure domain wall while the other three fields will support the monopole solution without affecting its winding number. Specifically, we can start with the $`\varphi _4`$ field in the configuration (12) and set up the other three fields in the configuration (9) centered not at $`(x,y,z)=(0,0,0)`$, but at $`(x,y,z)=(0,0,z_0)`$. The mixing of the winding number due to the product ansatz was a concern in our previous work , and now we can address it explicitly.

We also need to study the monopole winding number, $`Q_w`$, in order to determine whether the monopole unwinds. We can compute both the winding number and the winding number density. The topological charge associated with $`\mathrm{\Pi }_2()`$ (which in our case is $`𝒵`$) is the winding number. This number quantifies how many times the field configuration $`\varphi (x)`$ wraps around the vacuum manifold as $`x`$ ranges over a sphere $`S^2`$ in coordinate space. Hence the winding is defined as the homotopy class of the map $`\widehat{\varphi }=\frac{\varphi ^a}{|\varphi |}`$ from coordinate space to the vacuum manifold, known as isospace: $`\widehat{\varphi }:S_{space}^2S_{iso}^2`$.

For our O(3) theory, the winding number in a volume enclosed by the surface $`S_k`$ is:

$$N=\frac{1}{8\pi }𝑑S_kϵ^{ijk}ϵ_{abc}\widehat{\varphi ^a}_i\widehat{\varphi ^b}_j\widehat{\varphi ^c},$$
(13)
where the indices $`a,b,c`$ run from 1 to 3.

Since we consider the time evolution of the winding over the entire coordinate space, useful information is obtained from the winding number density. Using Stokes’ theorem in Eq. 13 by performing a total derivative on the surface flux in the integrand, the winding number becomes:

$$N=\frac{1}{8\pi }d^3xϵ_{abc}ϵ^{ijk}_i\widehat{\varphi ^a}_j\widehat{\varphi ^b}_k\widehat{\varphi ^c}.$$
(14)
By visualizing the evolution of the integrand in Eq. 14 we obtain information about the topological charge density over the whole space. In principle, we could track the charge by studying the surface flux alone, but to do this one has to know where the winding is in order to choose a correct surface of integration, or compute the surface integral over many different trial boundaries. Consequently, we have found the winding number density to be a more useful variable.

The one drawback of using the winding density to track the evolution of the monopole winding is that the density is initially a delta function, since $`ϵ_{abc}ϵ^{ijk}_i\widehat{\varphi ^a}_j\widehat{\varphi ^b}_k\widehat{\varphi ^c}\delta (r)`$ and the delta function is ill defined after we discretize the equations of motion. Consequently, while the observed peak in the winding density will give an accurate qualitative picture of the monopole’s location, the volume integral Eq. 14 is not conserved during the numerical evolution, even when we are looking at an isolated monopole.

## III Results and Discussion

Partly as a test of the code, we first investigated the monopole and domain wall in isolation. In both cases we found that the solitons are stable. Note that if the initial thickness of the solitons is not chosen ideally, there will be excess energy which can and will radiate away. Due to the reflective boundary conditions implemented, this energy reflects back towards the core. Nonetheless, in our studies of monopole - wall interactions this was not an obstruction, since the interaction occurred before the boundary effects became significant.

Our primary goal was to investigate the attractive force between the monopole and wall. We did this by choosing initial conditions corresponding to a monopole that was initially at rest adjacent to a domain wall, and then numerically evolving the three dimensional equations of motion. The evolution of the energy and winding densities for a typical set of parameters is displayed in Figs 1 and 2. The parameter choices for the plotted solution are $`\lambda _1=\eta _1=\eta _2=1`$ and $`\lambda _2=2`$, and the equations of motion were solved on a $`192^3`$ point grid.

The attractive force between the monopole and wall is obvious from the simulation, and does not depend on any special parameter choices. In Fig. 1, we see that the monopole starts moving towards the domain wall, but does not pass through it. At later times (not shown in the plots) the distance between the monopole center and the wall does not decrease uniformly, but oscillates. This is apparently due to the natural oscillations of the monopole that are observed in isolated monopoles if the initial profile is not a time independent solution, and may possibly be accentuated by the edge effects due to the finite size of the simulation volume. However, there is no ambiguity about the existence of the attractive force that leads to the initial displacement of the monopole.

We see from Fig. 2 that the monopole winding is initially localized at the core of the monopole. The charge remains localized on the wall at the point of impact for late times and does not dissipate (unwind) on the wall. This differs from the evolution seen when a monopole interacts with an embedded wall , in which case the monopole can unwind on the wall. Consequently, we will now examine the differences between the two models so as to understand the lack of unwinding with a stable domain wall, but in order to do this we must understand the mechanism which attracts the domain wall to the monopole.

The mechanism underlying the attraction between general monopoles and domain walls is based on the assumption that the fields evolve to minimize the sum of gradient and potential energies. As shown in , the gradient energy decreases as the monopole approaches the wall. Another way to understand the origin of the attraction is to focus on the change in potential energy and to associate it with a force between the monopole and the wall. If the force is attractive, it behaves like a tension which induces an acceleration along the line that connects the solitons.

From the form of the potential in (4) it is obvious that the potential energy associated with the monopole part of the configuration ($`\varphi _1\varphi _3`$) decreases with the distance between the monopole and the domain wall, as the monopole’s contribution to the potential energy being well localized. The potential energy is proportional to $`\varphi _4^2`$ and thus decreases as the distance between the monopole and the wall decreases. However, this is a short range force since the monopole profile function $`f(r)`$ approaches 1 at large $`r`$. The decrease in the gradient energy as the distance between monopole and domain wall shrinks is a longer range effect (for global symmetries - for local symmetries the gradients are screened by the gauge fields and hence also the force stemming from minimizing the gradient energy is short range). This effect comes from the fact that the volume of the region over which large gradients are localized shrinks as the defects approach. Note that this argument works for general classes of defects, and adds support to the first crucial assumption in the DLV  mechanism for solving the monopole problem.

We can understand the lack of unwinding in the present model by studying the energetics associated with a monopole unwinding on the domain wall. Let us assume that the monopole winding can spread out to a distance $`l`$ in the plane of the domain wall. When the monopole hits the wall its potential energy vanishes. What is left over is the gradient energy of both the monopole configuration (since the domain wall is topological, there is no reason to expect the $`\varphi _4`$ configuration to change). As the monopole winding tries to spread out, the $`\varphi _1\varphi _3`$ hedgehog configuration will be distorted so that for most values of $`x^2+y^2<l^2`$ the $`\varphi _3`$ component dominates. Thus, in order to estimate the energetics associated with a speading monopole charge, let us foucs on the contribution of the $`\varphi _3`$-component of the hedgehog field to the gradient energy. Making use of the symmetry in the plane of the wall, the energy $`E`$ can thus be approximated by integrating the gradient energy density of the $`\varphi _3`$-component of the monopole integrated over all $`z`$ out to a distance $`l`$ in the plane of the wall.

$$E>\pi l^2_w^{\mathrm{}}(\frac{\varphi }{z})^2𝑑z\pi l^2_w^{\mathrm{}}\frac{1}{z^2\pi ^2}𝑑z$$
(15)
where $`w`$ is the thickness of the wall along the $`z`$ axis. Thus,

$$E>\frac{1}{\pi }l^2\eta ^2w^1.$$
(16)
Hence we see that the energy grows as $`l^2`$. We can obtain an effective force by taking a derivative with respect to $`l`$. What we see is a restoring force $`F`$ which acts to keep the monopole winding charge confined to a localized region of the domain wall.

$$F=kl\mathrm{with}\mathrm{k}=\frac{1}{\pi }\eta ^2\mathrm{w}^1.$$
(17)
We conclude that this “restoring force” is the basis of the non-unwinding of the monopole on the domain wall.

## IV Conclusions

We have studied the interaction between monopoles and domain walls in a scalar field theory with symmetry group $`G=O(3)\times Z_2`$ symmetry group spontaneously broken to $`H=O(2)`$, a theory which admits both stable monopoles and domain walls. We demonstrated the existence of one of the two mechanisms which was crucial for the proposed solution of the monopole problem of , namely the attractive form between monopoles and domain walls. However, in contrast to the case of global monopoles interacting with embedded walls  we found that the charge of the global monopole will not unwind on the surface of the domain wall. The gradient energy associated with delocalizing the monopole topological charge yields a tension which keeps the winding confined to a local region on the wall. Hence, we conclude that global monopoles can not unwind unless there are gauge fields to cancel the tension associated with the gradients of the scalar fields. Therefore, the lack of unwinding of global defects on domain walls is a universal behavior.

## Acknowledgements

We thank Bill Unruh for a very useful discussion. Computational work in support of this research was performed at the Theoretical Physics Computing Facility at Brown University. This work has been supported in part by the US Department of Energy under Contract DE-FG02-91ER40688, TASK A, and one of us (SA) acknowledges support from the US Department of Education under the GAANN program.


text:
# Magnetic Fields and the Triaxiality of Molecular Cloud Cores

## 1. Introduction

Observations establish that magnetic fields play an important role in shaping the structure and dynamics of molecular clouds. Polarization maps reveal that there is a magnetic field threading molecular clouds that is ordered over large scales (e.g., Vrba, Strom, & Strom 1976; Goodman et al. 1990; Schleuning 1998) and Zeeman measurements establish that measured field strengths are invariably strong enough to yield a mass-to-flux ratio very close to the critical value (Crutcher 1999; hereafter C99), i.e., the magnetic energy approximately equals the gravitational binding energy. A problem with absorption maps of the polarization of background starlight is that the densest regions (the cloud cores) are not sampled (Goodman et al. 1995). However, new maps of polarized submillimeter emission (Matthews & Wilson 2000; Ward-Thompson et al. 2000) do sample these regions and provide a crucial test of theoretical scenarios involving magnetic fields; in particular, the relative alignment of the projected magnetic field with the projected minor axis of the core density distribution is an important diagnostic.

When making the above comparison, two important factors must be kept in mind: 1) the measured relative alignment is only between two projected quantities in the plane of the sky, and care must be taken to draw inferences about the relative orientation in the actual three dimensional object, and 2) theoretical models of magnetically dominated clouds (e.g., Mouschovias 1976; Lizano & Shu 1989; Fiedler & Mouschovias 1993), while predicting the important effect of flattening along the mean magnetic field direction, are usually restricted (for numerical convenience) to the physically unnecessary assumption of axisymmetry. In fact, it has been known for a long time (Mestel 1965; Lin, Mestel, & Shu 1965) that gravity amplifies anisotropic structure, so that an initially unstable fragment will first collapse along one dimension, forming a sheet, which will subsequently break into elongated filamentary structures. The presence of a large scale magnetic field can enhance this effect, as the first stage of collapse is preferentially along the mean field direction. The subsequent contraction lateral to the field will occur in qualitatively the same manner as described by Mestel (1965), particularly if the cloud is magnetically supercritical, but even if it is subcritical, since the evolution is still driven by gravity (Mouschovias 1978; Langer 1978), but on the ambipolar diffusion (rather than free-fall) time scale. Thus, the most general formation mechanism for a dense core in the magnetic field scenario yields an object flattened along the mean field direction and also having an elongated structure in the lateral direction; such an object is triaxial.

The properties of the projected shapes of triaxial bodies have been investigated in the context of galaxy studies. Stark (1977) has proven the important result that the isocontours of a triaxial body seen in projection are elliptical. Binney (1985) has extended this analysis, and his results can be applied to find the distribution of angular offsets $`\mathrm{\Psi }`$ between the projected field direction and the projected minor axis of a triaxial body, given that the body is primarily flattened along the direction of the mean magnetic field.

In § 2, we review observational evidence to establish that molecular cloud fragments are ubiquitously flattened along the mean magnetic field direction, and in § 3 calculate the distribution of $`\mathrm{\Psi }`$ for various triaxial bodies. In § 4 we discuss how the results relate to current observations and to other theoretical scenarios presented in the literature.

## 2. Cloud Flattening along Magnetic Field Lines

The magnetic field strength data of C99 also reveals that the correlation of the line-of-sight field strength $`B_{\mathrm{los}}`$ with the density $`\rho `$ is in apparent agreement with models of preferential flattening along the magnetic field ($`B\rho ^{1/2}`$; Mouschovias 1976) rather than isotropic contraction ($`B\rho ^{2/3}`$). Figure 1a (equivalent to Fig. 1 of C99) plots (in logarithmic form) the confirmed magnetic field strength detections of C99 against the measured number densities $`n`$. The solid line is the least squares best fit, with slope $`0.47\pm 0.08`$, and chi-squared error statistic $`\chi ^2=2.36`$. The correlation coefficient is 0.84.

We seek an even better correlation for the C99 data as follows. Consider a cloud that is flattened along the mean magnetic field direction. In this direction, force balance between an isotropic pressure and self-gravity yields

$$\rho _0\sigma _v^2=\frac{\pi }{2}G\mathrm{\Sigma }^2,$$
(1)
where $`\rho _0`$ is the density at the midplane, $`\sigma _v`$ is the total (thermal and non-thermal) one-dimensional velocity dispersion, $`\mathrm{\Sigma }`$ is the column density of the cloud, and $`G`$ is the gravitational constant. Here we assume that the effect of an external confining pressure is negligible. Furthermore, the relation

$$\frac{\mathrm{\Sigma }}{B}\mu (2\pi G^{1/2})^1$$
(2)
defines the dimensionless ratio $`\mu `$ of the mass-to-flux ratio to the critical value $`(2\pi G^{1/2})^1`$ for a disk (Nakano & Nakamura 1978). Combining equations (1) and (2) yields

$$B=(8\pi c_1)^{1/2}\frac{\sigma _v\rho ^{1/2}}{\mu }$$
(3)
(see Mouschovias 1991 for a similar relation), where $`c_1(1)`$ is an undetermined proportionality constant between $`\rho _0`$ and the mean density $`\rho `$. Clearly, the correlation in Figure 1a is imperfect because variations in $`\sigma _v`$ and $`\mu `$ are unaccounted for, and the observed quantity $`B_{\mathrm{los}}=B\mathrm{cos}\theta `$ (where $`\theta `$ is the angle between the field direction and the line of sight) does not correspond exactly to $`B`$. We improve the correlation by accounting for variations in the observable quantity $`\sigma _v`$. Figure 1b plots $`B_{\mathrm{los}}`$ versus $`\sigma _vn^{1/2}`$ in logarithmic form; equation (3) predicts a slope of 1 to the extent that $`\mu `$ is constant from cloud to cloud. We use values of $`\mathrm{\Delta }v=\sigma _v(8\mathrm{ln}2)^{1/2}`$ listed in Table 1 of C99. The correlation coefficient in this plot is 0.95. The solid line is the best least squares fit and has slope $`1.00\pm 0.09`$, and $`\chi ^2=0.79`$. The dashed line is the theoretical relation (3) for $`\sqrt{c_1}/\mu =1`$. The two lines coincide if $`\sqrt{c_1}/\mu =0.8`$. This exceptional correlation leads us to the conclusion that clouds (and cloud fragments) are flattened (at least moderately) along the magnetic field direction, even in the presence of turbulent motions, due to the dynamically important roles of self-gravity and the mean magnetic field. We note that the correlation found in Figure 1b is equivalent to that presented by Myers & Goodman (1988) in their Figure 1 (see also Mouschovias & Psaltis 1995) based on earlier magnetic field data. We also note that our fit to the magnetic field data is quite general; it does not require axisymmetry, or that the Alfvén speed $`v_A=B/(4\pi \rho )^{1/2}`$ be constant from one cloud to another. However, equation (3) implies that $`\sigma _v/v_A=(2c_1)^{1/2}\mu `$, so that the Alfvén mach number may be approximately constant, if non-thermal motions dominate thermal motions and $`\mu `$ does not vary much from one cloud to another.

Fig. 1  a) Plot of $`\mathrm{log}B_{\mathrm{los}}`$ versus $`\mathrm{log}n`$ from the C99 data set. The solid line is the best fit (see text). The correlation coefficient is 0.84. b) Plot of $`\mathrm{log}B_{\mathrm{los}}`$ versus $`\mathrm{log}\sigma _vn^{1/2}`$ for the same data. The correlation coefficient is 0.95. The solid line is the best fit, and the dashed line is a theoretical relation (eq.  for $`\sqrt{c_1}/\mu =1`$).

## 3. Projection of Triaxial Bodies

Let $`z`$ be the direction of the mean magnetic field. Since the mean field contributes cloud support against gravity only in the lateral ($`x`$ and $`y`$) directions, but the cloud need not be axisymmetric in the $`xy`$ plane, the most general cloud configuration is a triaxial body with the shortest dimension along the $`z`$axis. The simplest mathematical description of such a body is

$$x^2+\frac{y^2}{\zeta ^2}+\frac{z^2}{\xi ^2}=a^2,$$
(4)
where $`a`$ is a constant, and $`1\zeta \xi `$. Here we follow the geometric analysis of Binney (1985) to calculate projected quantities when the object is viewed from an observing angle $`(\theta ,\varphi )`$ (where the angles are defined on an imaginary viewing sphere and have their usual meaning in a spherical coordinate system). Defining the quantities

$$A\frac{\mathrm{cos}^2\theta }{\xi ^2}\left(\mathrm{sin}^2\varphi +\frac{\mathrm{cos}^2\varphi }{\zeta ^2}\right)+\frac{\mathrm{sin}^2\theta }{\zeta ^2},$$
(5)
$$B\mathrm{cos}\theta \mathrm{sin}2\varphi \left(1\frac{1}{\zeta ^2}\right)\frac{1}{\xi ^2},$$
(6)
and

$$C\left(\frac{\mathrm{sin}^2\varphi }{\zeta ^2}+\mathrm{cos}^2\varphi \right)\frac{1}{\xi ^2},$$
(7)
it can be shown that the projection in the sky of the triaxial body is an elliptical body of apparent axial ratio

$$q=\left(\frac{A+CD}{A+C+D}\right)^{1/2},$$
(8)
where $`D\sqrt{(AC)^2+B^2}`$. Furthermore, the quantity

$$\psi =\frac{1}{2}\mathrm{arctan}\left(\frac{B}{AC}\right)$$
(9)
can be used to find the absolute value of the angular offset between the projected $`z`$ axis and the apparent minor axis of the projected ellipse. Letting $`p=(AC)\mathrm{cos}2\psi +B\mathrm{sin}2\psi `$, the required positive definite angle is

$$\mathrm{\Psi }=\{\begin{array}{cc}|\psi |\hfill & \text{for }p0\hfill \\ \pi /2|\psi |\hfill & \text{for }p>0.\hfill \end{array}$$
(10)
Since $`\mathrm{\Psi }`$ corresponds to the offset between the projected magnetic field $`B_{}`$ and the projected minor axis of a core, we calculate the distribution of observed offsets $`\mathrm{\Psi }`$ using a Monte Carlo simulation. A large number $`N=10^6`$ of randomly oriented viewing angles $`(\theta ,\varphi )`$ are used to generate a distribution of $`\mathrm{\Psi }`$ for a given set of intrinsic axial ratios $`\zeta `$ and $`\xi `$ of a triaxial body. Figure 2 shows, for each of three cases $`\xi =0.1,0.3,`$ and 0.5, the probability distribution function $`f`$ versus $`\mathrm{\Psi }`$ for three different values of $`\zeta `$, as well as the dependence of the mean value $`\mathrm{\Psi }`$ of the distribution and the mean axial ratio $`q`$ on $`\zeta `$. The lower limit of $`\zeta =\xi `$ corresponds to a prolate object and the upper limit $`\zeta =0.9`$ is very nearly an oblate object.

Our calculated distributions show the robust result that there is always a peak in $`f`$ at $`\mathrm{\Psi }=0`$, so that it is the most probable value, but that there is a long tail towards all nonzero angles, including $`\mathrm{\Psi }=90^{}`$. Figures 2b, 2d, and 2f show how $`\mathrm{\Psi }`$ varies for clouds of different triaxial ratios; it becomes smaller as a cloud goes toward the oblate limit. For a perfectly oblate cloud $`(\zeta =1)`$, there would be no distribution, as all viewing angles yield $`\mathrm{\Psi }=0`$. The distributions for prolate clouds ($`\zeta =\xi `$)111Such clouds are distinct from existing models of prolate clouds (Tomisaka 1991; Fiege & Pudritz 2000) in which the mean magnetic field lies along the major axis; here the mean field lies along the minor axis. The cloud can be prolate ($`\zeta =\xi `$) by chance, but there is no intrinsic symmetry that favors such a configuration, and the pressure exerted by the mean field makes it likely that $`\zeta >\xi `$ for most clouds. have the largest possible values for $`\mathrm{\Psi }`$. The values of $`q`$ illustrate which triaxial models may be most likely. The crosses highlight the values of $`\zeta `$ for which $`q`$ is in the range $`0.50.65`$ that is consistent with observations of dense cores (e.g., Myers et al. 1991; Jijina, Myers, & Adams 1999).

While data sets that measure $`q`$ and $`\mathrm{\Psi }`$ for a large number of cores will ultimately settle the values (and distributions) of $`\zeta `$ and $`\xi `$, we do not currently favor highly flattened cases $`\xi 0.1`$ due to 1) the observational result that highly elongated cores are rarely observed, and 2) the theoretical result that flattening along the magnetic field yields cores of axial ratio $`0.250.33`$ (see Ciolek & Basu 2000 and references therein). Hence, the $`\xi =0.3`$ curves represent likely scenarios, although the $`\xi =0.5`$ curves cannot be ruled out observationally. The observational constraint on $`q`$ also allows the widest plausible range of $`\zeta `$ values when $`\xi =0.3`$, rather than imply near-oblate (for $`\xi =0.1`$) or near-prolate (for $`\xi =0.5`$) configurations. For $`\xi `$ in the range $`0.30.5`$, clouds with $`q`$ in the range $`0.50.65`$ have $`\mathrm{\Psi }`$ in the range $`10^{}30^{}`$.

To illustrate the point that no single observation can establish the relative orientation of the field and the minor axis of the intrinsic three-dimensional object, we present in Figure 3 a view of a triaxial body with $`(\xi ,\zeta )=(0.3,0.6)`$ from three positions on the viewing sphere. This illustrates that all possible orientations of $`B_{}`$ with the projected minor axis are possible, although the first panel shows the single most likely orientation and the middle panel is closest to the average orientation $`\mathrm{\Psi }=20^{}`$ for this case.

## 4. Discussion

Our predicted probability distributions, with their correlation toward $`\mathrm{\Psi }=0`$, are in good agreement with early submillimeter polarimetry of dense regions. Matthews & Wilson (2000) measure three distinct regions in OMC-3, of which two imply near perfect alignment of $`B_{}`$ with the projected minor axis, and one has $`\mathrm{\Psi }30^{}`$. Ward-Thompson et al. (2000) infer $`B_{}`$ in three dense cores, L1544, L183, and L43, and also find it to be correlated with the projected minor axes. They measure $`\mathrm{\Psi }=29^{},34^{}`$, and $`44^{}`$ for the three cores, respectively, although the last measurement is complicated by a nearby outflow and a weak polarization signal toward the core center. For L1544, the nonzero $`\mathrm{\Psi }`$ implies that the axisymmetric assumption in the magnetic model for that core developed by Ciolek & Basu (2000) needs to be relaxed, at least moderately. However, since that model does already account for flattening along the field lines, the magnetic field strength predictions would likely not change substantially.

Future large data sets of the offset angle $`\mathrm{\Psi }`$, as submillimeter polarization measurements become common, will properly test the distributions of $`\mathrm{\Psi }`$ presented in this paper. Such data will provide a clear observational test between the scenario presented here (which is already supported by the Zeeman data), in which gravity and the mean magnetic field play an important role in shaping core structure (even in the presence of turbulence) and those scenarios in which turbulent compression is dominant, leading to the claim (Ballesteros-Paredes, Vazquez-Semadeni, & Scalo 1999) that there is no correlation between the magnetic field direction and the cloud elongation.

I thank the anonymous referee for comments which helped improve the paper. This research was supported by a grant from the Natural Sciences and Engineering Research Council of Canada.

Fig. 2  Probability distribution function $`f`$ versus angle $`\mathrm{\Psi }`$ for various values of $`(\xi ,\zeta )`$, as well as the dependence of mean values $`\mathrm{\Psi }`$ and $`q`$ on $`\zeta `$ for individual values of $`\xi `$. a) $`f`$ versus $`\mathrm{\Psi }`$ for $`\xi =0.1`$ and $`\zeta =0.1,0.5,0.9`$. b) $`\mathrm{\Psi }`$ (solid line) and $`q`$ (dashed line) versus $`\zeta `$ for $`\xi =0.1`$. c) Same as a) but $`\xi =0.3`$ and $`\zeta =0.3,0.6,0.9`$. d) Same as b) but $`\xi =0.3`$. e) Same as a) but $`\xi =0.5`$ and $`\zeta =0.5,0.7,0.9`$. f) Same as b) but $`\xi =0.5`$. In b), d) and f), the crosses highlight the region in which $`q=0.50.65`$, consistent with observations.

Fig. 3  Simulated contours (solid lines) and mean polarization direction (dashed lines) for a triaxial body with axial ratios $`(\xi ,\zeta )=(0.3,0.6)`$ seen from three sets of viewing angles $`(\theta ,\varphi )`$. The The $`x^{}`$ and $`y^{}`$ axes are in the plane of the sky and are chosen to lie along the projected major and minor axes, and $`a^{}`$ is the apparent semi-major axis. The dashed lines lie along the projected direction of the intrinsic $`z`$ axis. a) $`(\theta ,\varphi )=(\pi /2,0)`$, yielding an apparent axis ratio $`q=0.3`$ and offset angle $`\mathrm{\Psi }=0^{}`$. b) $`(\theta ,\varphi )=(\pi /4,\pi /4)`$, yielding $`q=0.57`$ and $`\mathrm{\Psi }=28^{}`$. c) $`(\theta ,\varphi )=(\pi /6,0)`$, yielding $`q=0.68`$ and $`\mathrm{\Psi }=90^{}`$.


text:
# Instabilities in Josephson Ladders with Current Induced Magnetic Fields

## Abstract

We report on a theoretical analysis, consisting of both numerical and analytic work, of the stability of synchronization of a ladder array of Josephson junctions under the influence of current induced magnetic fields. Surprisingly, we find that as the ratio of the mutual to self inductance of the cells of the array is increased a region of unstable behavior occurs followed by reentrant stable synchronization. Analytic work tells us that in order to understand fully the cause of the observed instabilities the behavior of the vertical junctions, sometimes ignored in analytic analyses of ladder arrays, must be taken into account.

Ladder arrays of Josephson junctions are intriguing systems for a wealth of reasons: the possibility of phase locking, or synchronizing, a maximal subset of junctions suggests their use as microwave sources; they offer rich dynamical behavior, accessible to both theorists and experimentalists, in the field of coupled nonlinear oscillators (with recent interest in the prediction and observation of discrete rotobreathers); their complexity is between that of better understood one dimensional serial and parallel arrays and full 2D arrays (e. g. square arrays) and thus they offer a nice link between the two geometries; and ladder arrays can, under circumstances that are partially understood, be modeled by the discrete sine-Gordon (DSG) equation, which is itself a source of research interest among many. With a desire to understand better the conditions under which stable synchronization can occur, we consider ladder arrays biased with uniform dc bias currents $`I_B`$ greater than the critical currents of the junctions and include the effects of current induced magnetic fields (CIMF’s) via self and mutual inductances of the cells of the array. (see Fig. 1) Since the array is current-biased above the critical current there will be a nonzero voltage across some subset of junctions in the array. These “active” junctions are synchronized, or phase locked, if their voltages, after some initial transients, are identical functions of time. Furthermore, previous workers have established that the effects of CIMF’s may be important in determining the static and dynamics properties of arrays, and so it is natural to consider the effects of CIMF’s on synchronization as well.

In this Letter we present numerical and analytical evidence that mutual inductance between cells of a ladder array can lead to rich dynamical behavior, including destabilization of synchronization and reentrant synchronization as the relative size of the mutual to self inductance of the cells is increased. The instability results for a finite range of values of the mutual inductance and occurs when, for the resistively and capacitively shunted junction (RCSJ) model, the numerical solutions to the model equations diverge with time. In addition to our numerical solutions we also investigate the dynamics analytically. Namely, the two coupled RCSJ equations for the horizontal and vertical junctions can be approximately decoupled. The simplified equations allow us to calculate analytically the Floquet exponents, which measure the degree of stability of the synchronization. Then, by direct comparison of the analytic results for the Floquet exponents, based on the simplified equations, and the numerical results for the Floquet exponents, based on the full RCSJ equations, we learn valuable information about the roles of the horizontal and vertical junctions in the array. This technique appears to be a powerful way to analyze the relative importance of subsets of junctions that experience different local conditions.

Josephson junctions consist of superconducting islands separated by a thin layer of nonsuperconducting material. In the superconductors, the coherent motion of the paired electrons, or Cooper pairs, leads to a wavefunction of the form $`\mathrm{\Psi }=|\mathrm{\Psi }|e^{i\theta }`$, where $`\theta `$ is the macroscopic phase of the superconductor. The equations describing the dynamics of a single junction depend on the gauge-invariant phase difference, or Josephson phase, across the junction, $`\varphi =\theta _1\theta _2(2\pi /\mathrm{\Phi }_0)_1^2𝐀𝑑𝐥`$, where $`𝐀`$, the vector potential due to an external magnetic field, is integrated along a path from one side of the junction to the other. $`\mathrm{\Phi }_0=\mathrm{}/2e`$ is the magnetic flux quantum, where $`\mathrm{}`$ is Planck’s constant divided by $`2\pi `$ and $`e`$ is the electronic charge. We assume in this work that $`𝐀=0`$.

Consider a ladder array of underdamped junctions with $`N`$ cells and periodic boundary conditions, as shown in Fig. 1. Each junction has a McCumber parameter $`\beta _c2\pi CI_{cx}R^2/\mathrm{\Phi }_0`$, where $`C`$($`R`$) is the junction’s capacitance (resistance), and $`I_{cx}`$ is the critical current of a “horizontal” junction. We assume each cell of four junctions is described by a self-inductance $`L`$ and also has a mutual inductance $`M`$ ($`M>0`$) with each of its two nearest neighbors, where $`M<L`$. As an adjustable parameter in the theory, $`M/L`$ is allowed to range between zero and one. Basic physics arguments for a simple ladder with nearest neighbor mutual inductance precludes $`M/L`$ from exceeding $`0.5`$. Nevertheless, it is informative and worthwhile to study the behavior of the model equations over the range $`0<M/L<1`$. Furthermore, such theoretical work could pique the interest of others to try, using modern fabrication techniques, to enhance the mutual over the self inductance of the cells and thereby probe the broader range of $`M/L`$ values studied here. Alternatively, one could point to these results for $`M/L>0.5`$ as an interesting portent to the behavior that might be observed in arrays with overlapping cells or in three dimensional arrays.

An analysis based on the RCSJ equations and including the effect of induced magnetic flux leads to a pair of coupled equations for the Josephson phases of the horizontal and vertical junctions,

$$\beta _c\varphi ^{\prime \prime }+\varphi ^{}+\mathrm{sin}\varphi +\frac{1}{\beta _L}\underset{¯}{Z}^{\text{Tr}}\underset{¯}{X}^1\left(\underset{¯}{Z}\varphi +2\psi \right)=0,$$
(1)
$$\beta _c\psi ^{\prime \prime }+\psi ^{}+\alpha \mathrm{sin}\psi +\frac{1}{\beta _L}\underset{¯}{X}^1\left(\underset{¯}{Z}\varphi +2\psi \right)=0,$$
(2)
where $`\beta _L2\pi LI_{cx}/\mathrm{\Phi }_0`$ is the dimensionless self-inductance of a cell, and $`\alpha I_{cy}/I_{cx}`$ is the critical current anisotropy. The prime symbols denote differentiation with respect to dimensionless time, $`\tau t/t_c`$ where $`t_chbar/2eI_{cx}R`$. These equations are compactly represented in matrix notation, where $`\varphi `$ and $`\psi `$ are $`N`$-dimensional vectors representing the the Josephson phases of the horizontal and vertical junctions, respectively. $`\underset{¯}{Z}`$ is an $`N\times N`$ matrix that depends on geometry, and $`\underset{¯}{X}`$ is the dimensionless inductance matrix, also $`N\times N`$ in size. The diagonal terms of $`\underset{¯}{X}`$ represent the self-inductance of a given cell, i.e. $`X_{jj}=+1`$, while the mutual inductance of the nearest neighbor cells is represented by the terms $`X_{j,j\pm 1}=\mu _L`$, where $`\mu _LM/L`$ is the dimensionless mutual inductance. All other elements of $`\underset{¯}{X}`$ are zero. Eqs. 1 and 2 were solved numerically for $`\varphi `$, $`\psi `$, $`\varphi ^{}`$, and $`\psi ^{}`$ via a fourth-order Runge-Kutta algorithm as a function of the parameters $`N`$, $`\beta _c`$, $`\beta _L`$, $`\mu _L`$, and $`i_BI_B/I_{cx}`$ (the dimensionless bias current). The starting configuration consisted of random voltages and zero Josephson phases.

As described elsewhere for a similar system a stability analysis of the solutions to Eqs. 1 and 2 follows by letting $`\varphi =\varphi _0+\eta `$ and $`\psi =\psi _0+\delta `$, where $`\varphi _0`$ and $`\psi _0`$ are solutions to these equations. Eqs. 1 and 2 are linearized with respect to $`\eta `$ and $`\delta `$. The time dependence of the perturbations has the form $`\eta e^{\lambda t_c\tau }`$ and $`\delta e^{\mathrm{\Lambda }t_c\tau }`$, where $`\lambda `$ ($`\mathrm{\Lambda }`$) represent the Floquet exponents for the horizontal (vertical) junctions. If $`\text{Re}(\lambda t_c)>0`$ or $`\text{Re}(\mathrm{\Lambda }t_c)>0`$ we expect unstable behavior of the array for the given set of circuit parameters. In fact, we are interested in that exponent whose magnitude is closest to zero, as that describes the stability of the longest-lived mode of the array.

Figure 2(a) shows the minimum Floquet exponent as a function of $`\mu _L`$ for a 5-cell ladder. The three different sets of symbols represent three different sets of initial values for the voltages, $`\varphi ^{}`$ and $`\psi ^{}`$. The data show that the exponents do depend, to some degree, on the values of the starting voltages, although all the runs follow the same trend. As $`\mu _L`$ is increased from zero towards $`0.5`$, the stability of phase-locking increases, as shown by the negative exponent of growing magnitude, while the degree of stability decreases for increasing $`\mu _L`$ greater than approximately $`0.6`$. Even more interesting is the behavior of the ladder in the range $`0.5\mu _L0.6`$. For these values of the mutual inductance the ladder is actually unstable! This is evidenced by very rapidly growing phases and voltages with time as Eqs. 1 and 2 are numerically integrated. For $`N=5`$, the lower limit of this instability region is $`\mu _L^{(1)}=0.5`$ independent of other circuit parameters such as $`\beta _c`$ and $`\beta _L`$. The upper limit of this region, which we denote by $`\mu _L^{(2)}`$, depends on such quantities as the value of the starting voltages as well as on the value of $`\beta _L`$. For example, for a fixed set of starting voltages, we find that $`\mu _L^{(2)}`$ is a decreasing function of $`\beta _L`$. Also, it is interesting to note that this instability region does not appear at all if both the phases and the voltages are initialized to zero! (See discussion below for the reason for this behavior.)

Physically, this instability is due to a competition between the self-inductance of a given loop (say loop $`j`$), which wishes to have a current with a given sense of circulation, and the mutual inductance of the two neighboring loops ($`j\pm 1`$), which wish to have the current in loop $`j`$ flow in the opposite sense. We have also looked at ladders with $`N=6`$, $`7`$, $`8`$, and $`9`$. All show this instability in the vicinity of $`\mu _L=0.5`$. Indeed, we would expect this competition-induced instability to be independent of ladder size for the case of nearest-neighbor mutual inductance in that the onset of the instability should always occur at $`\mu _L=0.5`$.

Interestingly, ladders with $`N=7`$, $`8`$, and $`9`$ also exhibit a second instability region that the 5-cell ladder does not exhibit. Fig. 2(b) shows the Floquet exponents for $`N=7`$. This second instability region has an onset at a value of $`\mu _L^{(3)}>\mu _L^{(2)}`$ that is dependent on ladder size. (See cross-hatched region in Fig. 2(b).) We now turn to an analytic calculation of the Floquet exponents, which helps us understand the source of these instabilities.

A reasonable starting approximation is to ignore the coupling between the horizontal and vertical junctions but otherwise not to ignore the effects of the vertical junctions. That is, we let $`\psi =0`$ in eq. 1 and $`\varphi =0`$ in eq. 2. An analysis like that described in Ref.  applied to eq. 1 leads to a Mathieu equation describing the time dependence of the perturbations to the horizontal Josephson phases. In such a case, the corresponding Floquet exponents for the horizontal junctions can be calculated analytically. The result is

$$\text{Re}\left(\lambda _mt_c\right)=\frac{1}{2\beta _c}\pm \frac{1}{2\beta _c}\text{Re}\sqrt{14\omega _m^{(N)}\left(\frac{\beta _c}{\beta _L}\right)},$$
(3)
where we can think of $`\omega _m^{(N)}`$ as the effective normal-mode frequencies of the ladder. For example, $`\omega _m^{(5)}=(4\mathrm{sin}^2(m\pi /N)+2\mu _L\{\mathrm{cos}(2\pi m/N)\mathrm{cos}(4\pi m/N)\})/(\mu _L^2\mu _L1)`$, where $`0mN1`$. Note that $`\omega _m^{(N)}`$ is a function of $`\mu _L`$.

Eq. 3 was used to produce the solid curves in Figs. 2(a) and (b). Note that if, for particular values of $`\omega _m^{(N)}`$, $`\beta _c`$ and $`\beta _L`$, the argument of the square root in Eq. 3 is positive and larger than one, then at least one of the Floquet exponents will be positive, signaling unstable phase-locked solutions. In fact, since $`\beta _c>0`$ and $`\beta _L>0`$ such an instability will occur if $`\omega _m^{(N)}<0`$! Plots of $`\omega _1^{(N)}`$ versus $`\mu _L`$ for $`N=5`$ and $`7`$ (not shown here) demonstrate that $`\omega _1^{(5)}>0`$ for $`0\mu _L1`$, but $`\omega _1^{(7)}`$ is negative for $`\mu _L>0.8`$. We have checked that $`\omega _m^{(N)}>0`$ for $`m1`$ and for $`N=5`$ and $`7`$. Thus the cause of this instability in the 7-cell ladder for $`\mu _L>\mu _L^{(3)}`$ is the $`m=1`$ normal mode. That is, this instability is a geometrical effect, in that it does not occur for $`N=5`$ for example, and it is triggered by an effective normal mode frequency of the horizontal junctions becoming negative.

This analytic work, however, does not point to the horizontal junctions as the cause of the instability near $`\mu _L=0.5`$. To appreciate this behavior it is crucial to look to the vertical junctions. A procedure similar to that which led to eq. 3 leads to a set of effective Floquet exponents for the vertical junctions

$$\text{Re}\left(\mathrm{\Lambda }_mt_c\right)=\frac{1}{2\beta _c}\pm \frac{1}{2\beta _c}\sqrt{14\beta _c\left[\alpha \mathrm{cos}\psi _0+\frac{2\gamma _m^{(N)}}{\beta _L}\right]},$$
(4)
where the geometrical factor $`\gamma _m^{(N)}`$ is also a function of $`\mu _L`$ and is similar but not identical to $`\omega _m^{(N)}`$. In this case, the vertical junctions will exhibit an exponentially growing Josephson phase if $`\gamma _m^{(N)}<(\alpha \beta _L\mathrm{cos}\psi _0)/2.`$ Now a plot of $`\gamma _0^{(5)}`$ versus $`\mu _L`$ (see inset of Fig 2(a)) shows that the function abruptly becomes negative at $`\mu _L=0.5`$ and asymptotically approaches zero from the negative side as $`\mu _L`$ is increased further. (We have checked that $`\gamma _m^{(5)}>0`$ for $`m0`$. Also, we see similar behavior for the 7-cell ladders.) If we assume that $`\mathrm{cos}\psi _0>0`$, then the vertical junctions will be unstable for $`\gamma _m^{(N)}<0`$. Based on the behavior of $`\gamma _0^{(5)}`$ an instability region will exist for a range of $`\mu _L`$ values, $`\mu _L^{(1)}\mu _L\mu _L^{(2)}`$ where $`\mu _L^{(1)}=0.5`$ and $`\mu _L^{(2)}`$ will depend on $`\alpha `$, $`\beta _L`$, and $`\mathrm{cos}\psi _0`$. For example, as $`\beta _L`$ increases we expect that $`\mu _L^{(2)}`$ will decrease, i.e. approach a value of $`0.5`$. This agrees with the behavior of the numerical results for the Floquet exponents. Also, the inequality $`\gamma _m^{(N)}<(\alpha \beta _L\mathrm{cos}\psi _0)/2`$ suggests that the value of $`\mu _L^{(2)}`$ should depend on the value of $`\mathrm{cos}\psi _0`$. This is relevant to the numerical results in Fig. 2, where we see that the value of $`\mu _L^{(2)}`$ does indeed on the choice of the starting configuration of phases and voltages. In general, then, it is clear that the instability near $`\mu _L=0.5`$ originates with the vertical junctions and would thus be missed by an analysis that was based solely on the horizontal junctions. It is also clear why this instability does not appear numerically when both the Josephson phases and the voltages across the junctions are initialized to zero. In such a scenario, although the horizontal junctions may be active, the only possible solution for the vertical junctions is to keep zero voltages and Josephson phases for all times. Since we know this instability region is triggered by the vertical junctions, the vertical junctions have no chance to “go unstable” and thus the instability never appears.

We conclude that mutual inductance between cells of an underdamped ladder array has the effect of destabilizing synchronization for ranges of values of $`\mu _L`$, the ratio of the mutual to self inductance. These specific ranges of $`\mu _L`$ that lead to unstable behavior are geometry dependent. An analytic calculation of the Floquet exponents based on the horizontal junctions agrees with the numerical exponents, based on the full RCSJ equations, for those values of $`\mu _L`$ for which stable phase locking occur. To understand the cause of all the observed instabilities, however, it is crucial in the analytic work to consider the behavior of the vertical junctions.

Although some values of the mutual inductance used in these simulations can not be obtained in simple ladder arrays, this work suggests that experimentalists may wish to attempt fabrication of arrays that enhances the mutual over the self-inductance, perhaps making it possible to look for the rich dynamical behavior predicted here. Certainly researchers working on the problem of coherent emission from Josephson junction arrays should be aware this potential for unstable behavior exists.

The authors wish to thank Barbara Andereck, Tom Dillman, Steve Herbert, Mark Jarrell, and David Stroud for useful conversations. This research was funded by the Howard Hughes Medical Institute Undergraduate Biological Sciences Education Program grant #71196-529503 to Ohio Wesleyan University.


text:
# Precise calculation of the Casimir force between gold surfaces

## Abstract

We analyse the result of precise measurement of the Casimir force between bodies covered with gold. The values of the parameters used to extrapolate the gold dielectric function to low frequencies are very important and discussed in detail. The finite temperature effect is shown to exceed considerably the experimental errors. The upper limit on the force is found which is smaller than the measured force. Many experimental and theoretical uncertainties were excluded with gold covering and we conclude that, possibly, a new force has been detected at small separations between bodies.

The Casimir force  attracts considerable attention nowadays. At early stage it has been measured between two conducting plates  with large errors and much more accurately for plate and sphere covered with chromium . A new wave of interest raised when the force between metallized bodies has been measured with high precision . Theoretical prediction for the force is not so straightforward as expected. In the original Casimir work  the force $`F^{pl}(a)=\frac{\pi ^2}{240}\frac{\mathrm{}c}{a^4}`$ has been found as change in the zero point energy between two parallel perfectly conducting plates separated by the distance $`a`$. For sphere and plate the expression above has to be modified with the proximity force theorem (PFT)  to the following

$$F_c^0(a)=\frac{\pi ^3\mathrm{}c}{360}\frac{R}{a^3},$$
(1)
where $`R`$ is the sphere radius. Real metals are not ideal conductors and there is an important correction to (1) taking into account their real properties. In the Lifshitz theory , which generalizes the Casimir approach, the force is represented in terms of the metal dielectric function $`\epsilon \left(i\omega \right)`$ at imaginary frequencies. This theory was applied for description of the experiments in . Different but equivalent method has been developed in , where detailed calculations were made for $`Al`$, $`Au`$, and $`Cu`$ metallization of the bodies using the handbook optical data for these metals. The result for $`Au`$ is in good agreement with the torsion pendulum experiment .

For the atomic force microscope (AFM) experiment  $`Al`$ metallization has been used. This metal is oxidized easily and charges can be trapped in the oxide. To prevent influence of these charges on the force, $`Al`$ was covered with a thin layer ($`20`$ or $`8nm`$ thick) of $`Au/Pd`$. The conclusion has been made  that the result of the AFM experiment was in agreement with the theory. However, in all calculations  the role of $`Au/Pd`$ layer was ignored. Importance of this layer was stressed in  (see also ), where the upper limit on the force has been found using parameters of single-crystalline materials. This limit was smaller than the measured force in both experiments  and the discrepancy far exceeded the experimental errors. To exclude uncertainties connected with aluminum oxidation, $`Au/Pd`$ layer and clarify the origin of the discrepancy, it was proposed  to use $`Au`$ metallization instead of $`Al`$. Very recently the AFM experiment with $`Au`$ metallization has been carried out . In this experiment the residual potential was reduced to a negligible level so as the surface roughness. The force was measured for smaller separations up to $`62nm`$ with uncertainty $`3.5pN`$. In the previous experiments the photodiods response on the cantilever deflection was calibrated assuming that the plate-sphere system on the contact behaves as a rigid body. As the new experiment showed this natural assumption was not quite correct and different calibration method was proposed . To all appearance the excessive force revealed in the experimental data  originates mainly from the calibration effect.

The new experiment gives us a solid ground for verification of the theory. Though good agreement with the theory is reported , this conclusion cannot be considered as final because of the following reasons. Deviations between theory and experiment were significant at small separations but the authors averaged them over the whole range of separations. An error was admitted in extrapolation of the optical data for $`Au`$ to low frequencies. At last, the finite temperature correction was neglected but it is important at small separations . These are the main points we will discuss in this paper.

We start from the general expression for the Casimir force given by Lifshitz  modified for the case of sphere-disk geometry:

$$F(a)=\frac{kTR}{c^2}\underset{n=0}{\overset{\mathrm{}}{}}^{}\zeta _n^2\underset{1}{\overset{\mathrm{}}{}}𝑑pp\mathrm{ln}\left[(1G_1)(1G_2)\right],$$
(2)
where

$$G_1=\left(\frac{ps}{p+s}e^{p\zeta _na/c}\right)^2,G_2=\left(\frac{\epsilon ps}{\epsilon p+s}e^{p\zeta _na/c}\right)^2,$$

$$s=\sqrt{\epsilon \left(i\zeta _n\right)1+p^2},\zeta _n=2\pi nkT/\mathrm{}.$$
(3)
$`\epsilon \left(i\zeta _n\right)`$ is the dielectric function of $`Au`$ at imaginary frequencies. The prime over the sum sign indicates that the first $`n=0`$ term has to be taken with the factor $`1/2`$. Special care needs to treat this term. The formal reason is that $`\zeta _n^2`$ becomes zero but the integral over $`p`$ diverges. The physical reason is that this term corresponds to the static limit when for metallic bodies we should take $`\epsilon \mathrm{}`$. The resulting contribution in the force corresponds to the classical limit $`F_{cl}\left(a\right)`$ for metals

$$F_{cl}\left(a\right)=\frac{kTR}{4a^2}\zeta \left(3\right),$$
(4)
where $`\zeta \left(N\right)`$ is the zeta-function. If the metal is ideal, this term is canceled by the rest ones in the sum leaving only negligible temperature correction. The cancelation is incomplete for real metals  and the correction becomes important. To evaluate the Casimir force, the integral instead of the sum in (2) was used in . The difference between the sum and integral is just the temperature correction $`\mathrm{\Delta }_TF`$ which is nearly linear in $`T`$. The relative value of $`\mathrm{\Delta }_TF`$ becomes smaller for smaller separations but its absolute value increases. To have an estimate for the correction, we used the Drude approximation for $`\epsilon \left(i\zeta \right)`$ with the parameters of perfect single-crystalline gold  $`\omega _p=1.3710^{16}s^1`$, $`\omega _\tau =3.710^{13}s^1`$ and found for $`\mathrm{\Delta }_TF`$ the values $`4`$ and $`14pN`$ at $`a=100`$ and $`60nm`$, respectively. Obviously this correction cannot be ignored in precise calculations.

The function $`\epsilon \left(i\zeta _n\right)`$ cannot be measured directly for a given material but can be expressed via imaginary part of the dielectric function on the real axis $`\epsilon ^{\prime \prime }\left(\omega \right)`$ with the dispersion relation

$$\epsilon \left(i\zeta \right)1=\frac{2}{\pi }\underset{0}{\overset{\mathrm{}}{}}𝑑\omega \frac{\omega \epsilon ^{\prime \prime }\left(\omega \right)}{\omega ^2+\zeta ^2}.$$
(5)
Actual uncertainty in the force calculations arise from uncertainties in the function $`\epsilon ^{\prime \prime }\left(\omega \right)`$. Part of the available data on $`\epsilon ^{\prime \prime }\left(\omega \right)`$ is shown in Fig. 1. The solid line represents the data taken from  which were extrapolated to low frequencies as discussed in  (see below). The points marked with the open squares and triangles are taken from . The circles represent the data from . The points marked with the crosses  are the only data in the figure for the evaporated gold film while all the rest are for the bulk material. The data can be separated conventionally in two frequency regions with the boundary at the curve minimum $`\omega _1=3.210^{15}s^1`$. In the region $`\omega <\omega _1`$ the electron scattering on defects dominates in the material absorptivity. That is why in this range one can find different values for $`\epsilon ^{\prime \prime }\left(\omega \right)`$ corresponding to the samples prepared in different conditions. At $`\omega \omega _1`$ the dielectric function is well described by the Drude model

$$\epsilon \left(\omega \right)=\epsilon ^{}+i\epsilon ^{\prime \prime }=1\frac{\omega _p^2}{\omega \left(\omega +i\omega _\tau \right)}.$$
(6)
The parameters one can find by fitting the data for the real and imaginary parts of $`\epsilon \left(\omega \right)`$ with (6). It has been done  for the points marked by squares and triangles (rows 2 and 3 in Table I). Different procedure has been used for the solid line . The data in  are available for frequencies $`\mathrm{}\omega >0.1eV`$. They have been extrapolated to lower frequencies according to (6). The plasma frequency was estimated using its relation with the concentration $`n`$ of free electrons and supposing that every atom produce one electron: $`\omega _p=\sqrt{e^2n/m_e^{}\epsilon _0}`$. Here $`e`$ is the electron charge, $`\epsilon _0`$ is the free space permittivity, and $`m_e^{}`$ is the effective electron mass which is close but larger than $`m_e`$. Then at fixed $`\omega _p`$ the data  have been fitted with (6) to find $`\omega _\tau `$ (row 1 in Table I). The $`\omega _p`$ found in this way has actually the largest of possible values and, of course, the real material has smaller $`\omega _p`$ due to presence of the defects.

In contrast to low frequencies all the points at $`\omega >\omega _1`$ describe practically the same curve. That is because at high frequencies the interband absorption dominates and the defects become unimportant. Therefore, one can consider this part of the curve as the universal one (sample independent). We estimated uncertainty in $`\epsilon \left(i\zeta \right)`$ due to a reasonable variation of this part of the curve. It is less than 1% for important frequencies $`\zeta c/2a`$. The resulting uncertainty in the force will be much less than the experimental errors. One can separate contributions in $`\epsilon \left(i\zeta \right)`$ from different frequency regions and write $`\epsilon \left(i\zeta \right)=1+\epsilon _1+\epsilon _2+\epsilon _3`$. The contribution of low frequencies $`0<\omega <\omega _0`$, where the optical data are unavailable, can be expressed analytically using (6) in the dispersion relation:

$$\epsilon _1\left(i\zeta \right)=\frac{2}{\pi }\frac{\omega _p^2}{\zeta ^2\omega _\tau ^2}\left[\mathrm{tan}^1\left(\frac{\omega _0}{\omega _\tau }\right)\frac{\omega _\tau }{\zeta }\mathrm{tan}^1\left(\frac{\omega _0}{\zeta }\right)\right].$$

$`\epsilon _2`$ and $`\epsilon _3`$ can be found numerically using the data in the ranges $`\omega _0<\omega <\omega _1`$ and $`\omega >\omega _1,`$ respectively. To get an idea of the contributions of different intervals, we have found for the solid line in Fig. 1 $`\epsilon _1=26.6`$, $`\epsilon _2=8.3`$, and $`\epsilon _3=5.3`$ at $`\zeta =c/2a`$ with the smallest separation $`a=63nm`$ (the first point presented in ). Let us stress once more that $`\epsilon _3`$ does not depend on the defect concentration in the material while $`\epsilon _1`$ is quite sensitive to this parameter via actual values of $`\omega _p`$ and $`\omega _\tau `$ and its contribution dominates in $`\epsilon \left(i\zeta \right)`$.

Proceeding exactly as in  we successfully reproduced $`\epsilon \left(i\zeta \right)`$ given there in Fig. 2. However, we failed to reproduce the same function in  though the authors claimed to use the same data and the same values for $`\omega _p`$ and $`\omega _\tau `$. The difference between the curves presented in  and  is 46% at $`\zeta =10^{12}s^1`$ and decreases with frequency increase. To understand the origin of this difference, consider the limit $`\zeta 0`$ when the dielectric function behaves as $`\epsilon \left(i\zeta \right)`$ $`\left(\epsilon _0\rho \zeta \right)^1`$. Here $`\rho `$ is the static resistivity defined via the Drude parameters as $`\rho =\omega _\tau /\epsilon _0\omega _p^2`$. The parameters $`\omega _p`$ and $`\omega _\tau `$ used in  correspond to $`\rho =3.2\mu \mathrm{\Omega }cm`$. For unknown reason the function $`\epsilon \left(i\zeta \right)`$ presented in  at low frequencies agrees much better with $`\rho =2.3\mu \mathrm{\Omega }cm`$ that has no relation with the data in .

The Casimir force was evaluated in  and  without the temperature correction using the integral instead of the sum in (2). To check the procedure, we repeated the calculations there at $`a=100nm`$ and found for the reduction factor $`\eta =F/F_c^0`$ the values $`0.547`$ and $`0.559`$ to be compared with $`0.55`$  and $`0.56`$ , respectively. The small difference in $`\eta `$ results from different resistivities used to extrapolate the optical data. At $`a=63nm`$ this difference is not so harmless: we found for the force $`F=467pN`$ and $`477pN`$, respectively. The latter value is in agreement with the theoretical curve in  since the calculations have been done as in . The force has to be compared with the experimental value of $`491pN`$. The difference is 4 times larger than the experimental errors although the rms deviation of the experiment from theory for the whole range of separations is only $`3.8pN`$ .

Now let us take into consideration the temperature correction which is positive and will improve agreement with the experiment. We make calculations using Eq.(2) with $`\epsilon ^{\prime \prime }\left(\omega \right)`$ given by the solid line, squares, and triangles in Fig. 1. For the squares there is a gap in the interval $`6.310^{14}s^1<\omega <\omega _1`$ which is filled by the linear interpolation on a log-log scale. Different ways to fill the gap lead to uncertainty in the force $`\pm 1pN`$. The results of calculations are presented in Table I. Note that the force given in  accidentally coincides with the right one (row 1) because the error in $`\rho `$ was compensated by the temperature correction. The parameters in row 2 are quite close to those of perfect single-crystal. Therefore, the force has to be maximal in this case  due to maximal reflectivity. However, we see that it is a bit larger for row 1. The reason is that $`\omega _p`$ in  was not extracted from the optical data but set to be maximal by hand. Real material will have smaller $`\omega _p`$ that will lead to smaller force. The residual force defined as $`\mathrm{\Delta }F\left(a_i\right)=F_{exp}\left(a_i\right)F_{theor}(a_i)`$, where $`a_i`$ are the experimental separations, is shown in Fig. 2 for $`a_i<103nm`$. It has been calculated for the squares in Fig. 1 and can be considered as the minimal residual force. As we can see $`\mathrm{\Delta }F`$ exceeds 4 standard deviations in a few points of the closest approach. Parameters of evaporated gold films will deviate significantly from those of the perfect single-crystal because the films contain much more defects than the bulk material . Therefore, the residual force for real body covering will be definitely larger and will be significant at larger separations. One can easily improve the upper limit on the predicted Casimir force by measuring the resistivity of the gold films used in the experiment.

The experiment  has been substantially improved in comparison with  and many uncertainties have been excluded. If there is no an additional systematic error at small separation, the deviation between theory and experiment has to be considered seriously. In this connection it is very important to reproduce the result in independent measurements. Good potential has the micromechanical resonator experiment  which already has the necessary sensitivity but an additional work has to be done to control the separations $`100nm`$.

One can speculate that the residual force is explained by a new Yukawa force mediated by a light scalar boson. In this case interaction of two atoms is described by the potential

$$V_Y\left(r\right)=\alpha N_1N_2\frac{\mathrm{}c}{r}\mathrm{exp}(r/\lambda ),$$
(7)
where $`\alpha `$ is the dimensionless interaction constant, $`\lambda `$ is the Compton wavelength of a particle responsible for the interaction, and $`N_{1,2}`$ is the number of nucleons in atoms of the interacting bodies. With 95% cl the residual force is larger than $`1723.5=10pN`$. Integrating (7) for the bodies in the experimental configuration one finds the lower limit on the interaction constant

$$\alpha >\frac{6.2710^{25}e^{a/\lambda }}{11.74e^{h/\lambda }+0.75e^{2h/\lambda }}\left(\frac{100nm}{\lambda }\right)^3,$$
(8)
where $`a=63nm`$ is the minimal separation and $`h=96nm`$ is the gold film thickness. The interaction constant is restricted from the He burning stars  as $`\alpha <1.510^{22}`$. Then the allowed region for the wavelength is $`\lambda >33nm`$ or for the scalar boson mass $`m<38eV`$.

In  it is stated that the sphere-plate geometry cannot be used for verification of the Casimir effect since a rigorous theoretical calculation has never been carried out for this configuration. In this respect we would like to note that such a logic does not allow to verify the effect for finite-size plates, either, because the plates always can be considered as parts of spherical surfaces. From physical point of view it is obvious that the force for infinite plates will be a good approximation if $`La`$, where $`L`$ is the plate size. The same is true for a sphere and plate if $`Ra`$. One can discuss the applicability of PFT but there is no any reason to worry (see discussion in  and ).

When this work was finished we become aware about the paper , where the authors claimed that the Schwinger prescription  for the $`n=0`$ term in (2) cannot be applied for nonideal metals. As the result they found the force in the classical limit which is 2 times smaller than (4). We think that the statement needs further investigation (see, for example, ), but in any case our conclusion will not change since the Casimir force will be only smaller for the smaller $`n=0`$ term.

In conclusion, we found the upper limit on the Casimir force between gold covered plate and sphere using the handbook data for $`\epsilon \left(\omega \right)`$. The effect of nonzero temperature is shown to be considerably larger of the experimental errors. The excessive force found at small separations exceeds 4 standard deviations. The lower limit on the interaction constant is found assuming that the residual force is the manifestation of a new Yukawa interaction.

We thank R. Onofrio for attracting our attention to the papers ,  and for stimulating discussion. We thank also N.V. Mikheev for discussion and comments.