text:
# Doppler Tomograms from Hydrodynamical Models of Dwarf Nova Disks

## 1 Introduction

With their typical dimensions of less than $`10^4`$ arcseconds, Dwarf Nova (DN) disks are far too small to be resolved directly. However, an indirect observational insight into their structure became possible already in mid-eighties, when a powerful technique of Doppler tomography was introduced (for a recent review see Marsh 2001). Nowadays tomographic observations are a standard, but the interpretation of Doppler tomograms in terms of theoretical models of DN disks is often problematic. This is particularly well visible in the case of spiral waves, which have been suggested as one of the agents responsible for the angular momentum transfer through the disk (for a recent review see Boffin 2001). The excellent discussion of the subject can be found in a recent paper by Smak (2001), and there is no need to repeat it here.

The main problem associated with the interpretation of Doppler tomograms is related to the nature of the data derived from the theory. Model calculations yield spiral patterns in the distribution of disk surface density, while patterns in observational tomograms are related to the distribution of the emissivity in specific spectral lines. While comparing these two sets of data one usually makes an implicit assumption that emissivity is proportional to surface density, which certainly is an oversimplification.

A more sophisticated approach was presented by Steeghs and Stehle (1999), who based their Doppler tomograms on emission line profiles calculated from 2D disk models. For the origin of the lines they adopted a purely thermal model with local Planckian source function. Such model was criticized by Smak (2001), who pointed out that it requires factor of 100 overabundance of helium in order to reproduce the observed intensities, while the relative brightness of the features it produces in the tomograms is incompatible with observations.

Smak himself proposed to base the theoretical tomograms on the distribution of velocity divergence. He argued that at a given radial distance from the disk center the compression regions with $`\stackrel{}{v}<0`$ would be distinguished by a higher-than-average disk thickness (because their density and temperature would be higher). As a result, the surface layer of the disk would be better exposed to the irradiating flux from the white dwarf and the boundary layer, and a local enhancement in line emission would be observed. Based on the three-body model of gas flow in a close binary he identified regions of maximum compression in the orbital plane and showed that they well reproduced shape, location and relative intensities of the arch-like structures observed in Doppler tomograms of DN disks.

Motivated by his paper, we obtained two- and three-dimensional hydrodynamical models of DN disks and calculated their tomograms. The details of the modelling procedure are given in Section 2. The results are presented in Section 3 and discussed in Section 4.

## 2 Numerical methods, input physics and initial conditions

All models presented here were obtained with the help of the ZEUS-3D code (Clarke & Norman 1994, Clarke 1996). The original code was modified to include conservative angular momentum transport (Kley 1998). Spherical coordinates $`(r,\theta ,\varphi )`$ centered on the primary and corotating with the system were used. The grid was extending from $`0`$ to $`2\pi `$ in $`\varphi `$, from $`0`$ to $`0.2\pi `$ in $`\theta `$, and from $`r_{\mathrm{in}}=0.1a`$ to $`r_{\mathrm{out}}=0.5a`$ in $`r`$, with $`a`$ standing for the orbital separation. Grid spacing was uniform in $`(\theta ,\varphi )`$ and logarithmic in $`r`$, resulting in zones of identical shape. After a few experiments we decided to limit the resolution to $`100\times 20\times 64`$ zones in $`r,\theta ,\varphi `$ directions, respectively (test runs, with resolution increased by a factor of 2 in $`r,\theta ,\varphi `$ consecutively, did not introduce any significant changes into the models).

A standard periodic boundary condition was imposed at $`\varphi =2\pi `$, and symmetry with respect to the orbital plane was assumed, implying a reflecting boundary condition at $`\theta =0`$. A free outflow was allowed for at $`r_{\mathrm{in}}`$ and $`r_{\mathrm{out}}`$. In the four innermost radial zones the radial velocity was reduced by 5% at every time-step, preventing the reflection of waves from the inner boundary of the grid. To check whether this damping procedure did not influence the structure of the disk or the shape of the spiral pattern, we calculated model B2 with $`r_{\mathrm{in}}`$ moved to $`0.045a`$ (see Table 1). It was found that shifting the inner grid boundary toward the white dwarf did not introduce any significant changes in the model.

The simulations did not include explicit viscosity (the von Neumann & Richtmyer and scalar linear artificial viscosities originally implemented in ZEUS were only used with coefficients $`C_1=0.5`$ and $`C_2=2.0`$, where $`C_1`$ is responsible for the magnitude of the artificial viscous pressure, and $`C_2`$ is a shock-spreading parameter, see the definitions in Stone and Norman (1992)). The energy equation was not solved; a polytropic equation of state ($`p=\kappa \rho ^\gamma `$) with $`\gamma =5/3`$ was employed instead. In a system of units in which gravitational constant, orbital separation, and primary’s mass are all equal to 1, the value of $`\kappa `$ was set to 6500. To mimic U Gem - like systems, all models had the same mass ratio, $`q=0.5`$. The stream flowing from the secondary through the $`\mathrm{L}_1`$ point was not included.

Every simulation consisted of three phases (relaxation, switch-on, and proper). The mass ratio was set to 0 in the relaxation phase and to 0.5 in the proper phase, while in the switch-on phase it was linearly increasing in time. At the beginning of each simulation ($`t=0`$) the grid was initialized with an exponential density distribution

$$\rho (r,\theta )=\mathrm{max}(\rho _0\mathrm{e}^{\alpha r^2\mathrm{sin}^2\theta },\rho _{\mathrm{min}}),$$
(1)
where

$$\alpha =\frac{GM_1}{2r^3}\frac{1}{c_{\mathrm{s},0}^2},$$
(2)
and

$$c_{\mathrm{s},0}^2=\frac{p}{\rho }|_{\theta =0}=\gamma \kappa \rho _0^{\gamma 1}.$$
(3)
In our system of units the value of the midplane density, $`\rho _0`$, was set to $`10^8`$, corresponding to $`210^8\mathrm{g}\mathrm{cm}^3`$ in U Gem (with U Gem parameters taken from Groot 2001). The azimuthal velocity of the disk was given a purely Keplerian pattern, and the remaining two velocity components were set to $`0`$. Since the content of the grid was not in hydrostatic equilibrium, we allowed it to relax for $`3.9`$ orbital periods ($`P_{\mathrm{orb}}`$). Throughout the relaxation phase the model was strictly axisymmetric, so that it was possible to speed the computations up by reducing the number of angular grid points to 2.

Because of extremely steep vertical gradients of density at the surface of the disk it was necessary to introduce a density limit. Every time step the grid was scanned for cells with $`\rho <\rho _{\mathrm{min}}`$, and whenever such a cell was found $`\rho `$ was reset to $`\rho _{\mathrm{min}}`$. We wanted $`\rho _{\mathrm{min}}`$ to be small enough to minimize side-effects caused by the newly-added matter falling onto the disk, and, simultaneously, large enough to avoid excessive computational slowdowns due to formation of strong shocks in the rarefied medium above the disk. After a few experiments $`\rho _{\mathrm{min}}`$ was set to $`10^{13}`$ for all models. At the end of the relaxation phase the midplane density of the disk increased, reaching up to $`410^8`$ (corresponding to $`810^8\mathrm{g}\mathrm{cm}^3`$ in U Gem).

At the beginning of the switch-on phase the relaxed model was mapped onto the standard grid, and the secondary’s gravity was “switched on”. The final value of $`q`$ was achieved at $`t=5.8P_{\mathrm{orb}}`$. At the end of the switch on phase the total mass contained in the grid, and scaled to U Gem, $`M_{\mathrm{disk}}`$, was equal $`10^{24}\mathrm{g}`$). The proper phase with $`q=0.5`$ lasted for another $`3.9P_{\mathrm{orb}}`$ so that the simulation was ended at $`t9.7P_{\mathrm{orb}}`$. By that time shape and location of the disk edge stabilized, and a stationary spiral pattern developed in the disk.

## 3 Results

The tidal forces affect the relaxed disk in two ways. First, some material is stripped from its outer edge and driven out of the grid through the outer grid boundary. Second, angular momentum is removed from the remaining material and transferred into the orbital momentum of the binary, causing the disk to shrink. In the three-body approximation, the radius of the disk cannot be larger than the radius of the largest non-intersecting orbit, $`r_{\mathrm{nis}}^{\mathrm{max}}`$. U Gem - like systems with $`q=0.5`$ have $`r_{\mathrm{nis}}^{\mathrm{max}}0.3`$, and in fact at the end of the simulation the disk barely extends beyond $`r0.3`$. The rest of the gas originally located at $`0.3<r<0.5`$ now resides in a ring-like density enhancement between $`r0.15`$ and $`r0.3`$.

The ring is markedly elliptical, but, when averaged over the azimuthal angle, it shows a well-defined density maximum at $`r0.19`$, i.e. slightly beyond the circularization radius ($`r_{\mathrm{circ}}=0.16`$ for $`q=0.5`$). The maximum density is factor of $`2.5`$ higher than the midplane density of the inner disk ($`r<0.15`$). The ratio $`h/r`$, where $`h`$ is the half-thickness of the disk, varies from $`0.1`$ in the inner disk to $`0.2`$ in the ring. The overall structure of the final model is reminiscent of the one expected at the early phase of the outburst, when the outer radius of the disk just begins to increase. However, because of the polytropic equation of state we employ, the interior of the disk is unrealistically hot (for a disk composed of pure hydrogen $`T`$ grows from $`10^4\mathrm{K}`$ at the surface of the disk to $`6.710^5\mathrm{K}`$ at the density maximum). We find that both 2-D and 3-D calculations produce disks of nearly the same shape and extent (Fig. 1). Below we shall argue, however, that fully three-dimensional models are needed to properly interpret the Doppler tomograms.

Both 2-D and 3-D models develop spiral shocks shown in Figs. 1 and 2. At a first glance, location and inclination of the shocks do not significantly depend on the number of dimensions of the model. Significant differences become visible when compression regions ($`\stackrel{}{v}<0`$) are compared: while model A has a clear two-armed pattern, three arms are present in models B1 and B2. We speculate that the third arm may originate due to tidal forcing of the disk matter in the direction perpendicular to the orbital plane (an effect entirely absent in 2D). Thus, to the long dispute on whether spiral shocks can exist in three dimensions, or rather their existence is limited to the two-dimensional world, we add a vote in favor of the first possibility. In 3D the shocks are definitely there, but their pattern is different than in 2D. Obviously, the validity of this conclusion is limited to hot polytropic disks with $`0.1h/r0.2`$.

Following the approach of Smak (2001), we obtained Doppler tomograms of the compressional power due to tidal forces, $`p\mathrm{d}V`$.The distribution of brightness on the $`(v_x,v_y)`$ plane was calculated from the integral

$$L_{p\mathrm{d}V}(v_x,v_y)=_V\mathrm{min}(p\stackrel{}{u},0)(u_x,v_x,\delta _{v_x})(u_y,v_y,\delta _{v_y})dV,$$
(4)
where ($`u_x,u_y`$) are velocity components of the volume element $`\mathrm{d}V`$, and the so-called boxcar function $``$ is defined as

$$(u,v,\delta )=\{\begin{array}{cc}1\hfill & |uv|<\frac{\delta }{2}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$
(5)
The resolution parameter $`\delta `$, related to the resolution of images ($`400\times 400`$ pixels), was set to $`6v_{\mathrm{orb}}/400`$.

In the three-body approximation of Smak (2001), the inner disk is essentially Keplerian, and its Doppler tomogram cannot show any nonaxisymmetric structures in the high-velocity range. However, in a more realistic polytropic disk the spiral shocks excited at its outer edge propagate deeply into the high-velocity regions. As a result, our tomograms (Fig. 3) show extended spirals instead of two crescent-shaped maxima reported by Smak (2001). Such spirals are not observed in real disks (Groot 2001). To account for the limited resolution of the observational data we blurred the tomograms to such a degree that the width of the brightest segments of the spiral became comparable to the width of the observed features ($`400÷500\mathrm{km}\mathrm{s}^1`$; see Groot 2001). However, the spiral pattern extending up to velocities of $`1000\mathrm{km}\mathrm{s}^1`$ was still clearly visible. Moreover, both location and relative intensity of the brightest segments of the spiral did not agree with those observed in U Gem.

Should this disagreement be regarded as an argument against the presence of spiral waves in CV disks? Certainly not. As we already indicated in the Introduction, the observed tomograms refer to the distribution of the emissivity in specific spectral lines rather than to the distribution of physical parameters directly obtainable from hydrodynamical simulations. The presently available hydrocodes are not sophisticated enough to predict the detailed spectrum of the disk, and the correspondence between those two sets of data is by no means clear. However, the models can yield data much more closely related to line emissivity than simple physical parameters or their combinations.

To obtain such data, we assume the line emission to originate mainly due to the irradiation of the surface layer of the disk by the central white dwarf and/or boundary layer (cf. Robinson et al. 1993, Smak 1991). Since our models are not detailed enough to resolve the surface layer, we assume that the lower boundary of the layer coincides with a constant density surface $`S_\mathrm{l}`$ at which $`\rho \rho _\mathrm{l}=10^8`$. Typical distance between $`S_\mathrm{l}`$ and the midplane of the disk, $`h_\mathrm{l}`$ was such, that $`h_\mathrm{l}/r=0.08`$ and $`h_\mathrm{l}/r=0.2`$ in the inner disk and in the outer ring, respectively. Further, we assume that the line flux from each element of the layer is proportional to the mass contained within that element, $`\rho \mathrm{d}V`$, multiplied by the irradiating flux. For simplicity, we also assume that all irradiating photons are emitted from a point source located at the centre of the white dwarf, so that the irradiating flux is given by $`/r^2`$, where $`r`$ is the radial coordinate of the volume element $`\mathrm{d}V`$, and $`=_{\mathrm{wd}}+_{\mathrm{bl}}`$ is the combined luminosity of the white dwarf and the boundary layer. The distribution of brightness on the $`(v_x,v_y)`$ plane is now given by the integral

$$L_{\mathrm{irr}}(v_x,v_y)=_V\rho \frac{}{r^2}𝒮(r,\theta ,\varphi )(u_x,v_x,\delta _{v_x})(u_y,v_y,\delta _{v_y})dV.$$
(6)
The function $`𝒮(r,\theta ,\varphi )`$ describes the shadow cast by $`S_\mathrm{l}`$, and it is given by

$$𝒮(r,\theta ,\varphi )=1\left(_0^r(\rho (r^{},\theta ,\varphi )\rho _\mathrm{l})𝑑r^{}\right),$$
(7)
where

$$(x)=\{\begin{array}{cc}1\hfill & x>0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$
(8)
is the Heaviside step function. The third velocity component, $`v_z`$, was neglected in (6) because nearly everywhere in the disk its value was smaller than $`5\%`$ of the local azimuthal velocity, $`v_\varphi `$. Formally, the integral in Eq. 6 subtends the whole space above $`S_\mathrm{l}`$. Practically, due to steeply falling density, only volume elements closest to $`S_\mathrm{l}`$ contribute to it significantly. The boundary layer above $`S_\mathrm{l}`$ has a mass of $`0.15M_{\mathrm{disk}}`$, but only about 25% of its volume is directly illuminated.

Location and relative intensity of the brightest areas on tomograms resulting from the irradiation approach (Fig. 4) agree rather well with those observed by Groot (2001) at an advanced outburst phase of U Gem (his Fig. 2, Episode 2). The major discrepancy is the bright ring visible in our tomograms at $`\sqrt{v_x^2+v_y^2}1000\mathrm{km}\mathrm{s}^1`$. In both models the ring is too bright compared to the observational data, but relative to the maximum intensity obtained in the model it is weaker in B2 where the disk extends down to $`r=0.045`$. We conclude that in B1 the ring is enhanced by a spurious contribution from the inner edge of the disk at $`r=0.1`$. The ring in B2 would be still weaker if absorption of the irradiating flux by the gas in the surface layer was taken into account in equation (6). Unfortunately, the present models are too crude for such an operation to be reliable. It is clear, however, that the effect of absorption should be particularly strong for $`r0.1`$, where $`h/r`$ is nearly constant (see Fig. 5), and the irradiating quanta propagate nearly parallel to $`S_\mathrm{l}`$. Further reduction in ring intensity could probably be achieved if the inner boundary of the grid was moved even closer to the white dwarf, as the shadow cast by $`S_\mathrm{l}`$ at $`r<0.045`$ might partly screen the region at $`r0.1`$ where the ring originates. On the other hand, in some phases of the activity cycle the intensity of the brightest areas in our irradiation tomograms is underestimated relative to the ring. This is because substantial $`p\mathrm{d}V`$ work is done by tidal forces directly on the gas in the surface layer of the outer disk, where the low-velocity emission originates.

In fact, the heating rate per unit mass, $`p\mathrm{d}V/\rho `$, reaches a clear maximum just below the surface of the outer disk (see Fig. 5). For the case of U Gem the height of this maximum is $`10^{11}\mathrm{erg}\mathrm{g}^1\mathrm{s}^1`$. With a typical outburst accretion rate of $`\dot{M}=310^{18}\mathrm{g}\mathrm{s}^1`$, and a white dwarf radius $`R_{\mathrm{wd}}=410^8\mathrm{cm}`$ we get $`L_{\mathrm{bl}}=\frac{1}{2}GM_1\dot{M}/R_{\mathrm{wd}}610^{35}\mathrm{erg}\mathrm{s}^1`$. The illuminated mass in the region of maximum $`p\mathrm{d}V/\rho `$ between $`r0.15a`$ and $`r0.25a`$ approaches $`0.025M_{\mathrm{disk}}`$, and it is distributed within a solid angle of $`0.3\pi `$. Assuming that the whole incident flux is absorbed there we obtain a radiative energy input of $`210^{12}\mathrm{erg}\mathrm{g}^1\mathrm{s}^1`$, and we see that during the outburst the contribution to the line flux from $`p\mathrm{d}V`$ heating is rather small. However, $`p\mathrm{d}V`$ dominates just before the outburst, when accretion rate is factor of $`100`$ lower.

The subsurface maximum of $`p\mathrm{d}V/\rho `$ in Fig. 5 originates mainly due to tidal forcing in the direction perpendicular to the orbit. It also contains a contribution from the ambient gas falling onto the disk; however the “rainfall” heating is much less efficient than the tidal one. We checked this by re-running simulation B1 with $`\rho _{\mathrm{min}}`$ reduced to $`310^{14}`$: at $`t8.0P_{\mathrm{orb}}`$ virtually no changes were seen in the distribution of $`p\mathrm{d}V/\rho `$.

## 4 Discussion

As discussed in Section 3, we find that spiral waves efficiently propagate from excitation regions at the outer edge of the disk toward the white dwarf, reaching to at least $`0.05a`$. This conclusion concerns both two-dimensional and three-dimensional models; it is however limited to hot polytropic disks presented in this paper. The main spiral features seen in the density distribution (two-dimensional case) and in the surface density distribution (three-dimensional case) are hard to distinguish. On the other hand, clear differences are visible in the distributions of the tidal heating rate, $`p\mathrm{d}V`$: in 3D the two main spiral arms are less tightly wound than in 2D, and a weaker third arm is excited. We suggest that the third arm may originate from tidal forcing in the direction perpendicular to the orbital plane. The effects of this forcing seem to be responsible for the origin of the clear maximum of tidal heating rate per unit mass, $`p\mathrm{d}V/\rho `$, which is located away from the midplane in subsurface layers of the outer disk.

Doppler tomograms of tidal heating rate derived from both 2D and 3D models correlate rather poorly with observed tomograms of U Gem (Groot 2001). A better agreement (but still not entirely satisfactory) is obtained for tomograms of the irradiation flux from the white dwarf through the surface layer of the disk. The brightest areas of such tomograms coincide with arches observed in U Gem at an advanced stage of the outburst. The irradiation tomograms can be derived from 3D models only, which indicates that fully three-dimensional modelling is needed for a reliable interpretation of the observational data on DN disks.

According to our results the arches originate in the outer part of the disk, fairly high above the midplane ($`h>0.1r`$). For this to happen, the outer disk would have to be substantially bulged. The bulging phenomenon may be explained within the following (not entirely new) scenario: prior to the outburst the gas transferred from the secondary mainly collects in a ring at the circularization radius, and only partly accretes through the disk onto the white dwarf. The ring expands as the gas flows in, but it remains cool until heating from tidal forcing in the orbital plane grows so strong that it cannot be balanced by radiative cooling. The ring begins to expand even more rapidly, and within it the gas located away from the midplane begins to receive additional internal energy from tidal forcing in the direction perpendicular to the orbital plane. Eventually, strong spiral shocks develop, and a dynamical instability of the kind described by Różyczka & Spruit (1993) sets in.

Obviously, such scenario cannot be the whole story, as it is not linked to the thermal instability believed to be at least partly responsible for eruptive phenomena in DN and other classes of cataclysmic variables. Nevertheless, it seems to indicate a promising direction of further research.

Acknowledgments. We are greatly indebted to Professor J. Smak for enlightening discussions and critical reading of the manuscript. This work was supported by the Committee for Scientific Research through the grant 2P03D 01419.

References

Boffin, H. M. J. 2001, “Astrotomography. Indirect Imaging Methods in Astronomy”, Eds. H. M. J. Boffin, D. Steeghs, J. Cuypers, Lecture Notes in Physics 573, p.69
Clarke, D. A. 1996, Astrophys. J., 457, 291
Clarke, D. A. and Norman, M. L., 1994, LCA preprint 7 (Urbana-Champaign: Univ. Ill.)
Groot, P. J. 2001, Astrophys. J., 551, L89
Kley, W. 1998, Astron. Astrophys., 338, L37
Marsh, T. R. 2001, “Astrotomography. Indirect Imaging Methods in Astronomy”, Eds. H. M. J. Boffin, D. Steeghs, J. Cuypers, Lecture Notes in Physics 573, p.1
Robinson, E. L., Marsh, T. R. and Smak, J. 1993, “Accretion Disks in Compact Stellar Systems”, Ed. J. C. Wheeler (Singapore: World Scientific), p. 75
Różyczka, M., and Spruit, H. 1993, Astrophys. J., 417, 677
Smak, J. 1991, “Structure and Emission Properties of Accretion Disks”, Eds. C. Bertout, S. Collin-Souffrin, J. P. Lasota, J. Tran Thanh Van (Paris: Editions Frontiéres), p. 247
Smak, J., 2001 Acta Astron., 51, 295
Steeghs, D. and Stehle, R. 1999, MNRAS, 307, 99


format: text

text:
# Isotropic microwave composites with negative refraction

## Abstract

The properties of artificial isotropic microwave composites which would possess simultaneously negative permittivity and permeability are studied theoretically. Four kinds of composites are considered. Two of them concern media with split-ring resonator particles with different particle arrangements. The other two are realized with Omega particles. An analytical antenna model of the electromagnetic behavior of a split-ring resonator (SRR) is suggested and verified by numerical simulations. Next, material parameters of composite media made with SRR or Omega particles are calculated. It is shown that the Omega-composites are more prospective for obtaining negative real values of permittivity and permeability than the split-ring resonators.

## 1 Introduction

V. G. Veselago studied in sixties  the electromagnetic properties of materials whose permittivity and permeability have negative real parts simultaneously. He showed that these media exhibit unusual effects like anomalous negative refraction of electromagnetic waves and existence of backward waves (when the group velocity and the Poynting vector are directed opposite with respect to the phase velocity). This concept implies many other interesting features of wave propagation that were reviewed in .

These Veselago media do not exist in nature. The interest in creating such media artificially, for optical applications, has been indicated by J. Pendry in . This kind of medium can be obtained with arrays of conducting inclusions of a special shape. It is obvious that the response of these materials will be resonant. Therefore, negative refraction is possible only within a more or less narrow frequency band (or bands) . This ”negative” band can be a part of a resonant band of inclusions that are used in the composite. This is quite clear analyzing the Lorentz model of the resonant dispersion.

At this stage of studies, since it is easier to use resonant inclusions with millimeter dimensions, it seems to be reasonable to focus on applications for the microwave region of the spectrum (optic composites will need micrometer particles) and to study their response experimentally. How to prepare these composites in a best way? Which shape of inclusions is more suitable? These are the general questions addressed by this article.

Works that have theoretically predicted negative permittivity for composite media have been reviewed in . However, publications containing the same results for the permeability are not so common. In  and  it was indicated that it was a very difficult problem to obtain a magnetic resonant response of inclusions (in order to get a negative real part of $`\mu _{eff}`$). Probably, it was the reason for which these authors have chosen a rather complicated particle named as split-ring resonator (SRR), see Fig. 1. They maintained in these works, that SRR has only a resonant magnetic polarizability and no resonant electric one.

In  it was suggested that to obtain negative permittivity simultaneously with negative permeability. SRR inclusions must be combined with a lattice of straight wires. This complex system is a very cumbersome technical solution. The only special case where wave propagates normally to the wires (which corresponds to ) is suitable to observe the negative refraction in such a structure. It is clear that the structures from  cannot be considered as an isotropic backward-wave composite.

In this paper we show that the SRR particles alone allow to obtain simultaneously negative permeability and permittivity. Indeed, these inclusions also possess a resonant electric polarizability that has not been discovered in  and . We also point out that the analytical model of SRR suggested in  does not give a satisfactory explanation of the frequency behavior of effective permeability. Moreover, it is possible to realize an isotropic composite with negative refraction choosing properly the geometry, arrangement and concentration of SRR particles in a host dielectric matrix.

However, we have found that SRR is not the best choice for obtaining a backward-wave medium at microwaves. Actually, the resonant electric response of SRR is weaker than the magnetic one. We show that it is better to use an inclusion known as Omega particle (or $`\mathrm{\Omega }`$-shaped conducting particle) suggested in 1992 by Saadoum and Engheta  since its magnetic and electric responses of are of the same order.

Numerical and analytical modeling of Omega particles have been recently made in - and the properties of uniaxial and planar Omega composite media have been derived. It appears that these media are prospective for creating low-reflecting shields and antenna radomes in microwave applications. With this kind of objectives, in 1995-1996, we studied theoretically uniaxial Omega composites and noticed that sometimes the real parts of both permittivity and permeability became negative, especially in the cases with high densities of particles. Thus, as it was not suitable for our goals, the experimental samples of Omega composites which were prepared under guidance of S.A. Tretyakov  were dilute mixtures in order to avoid the apparition of negative properties. Negative permittivity and permeability of Omega composites has been predicted in the paper . However, the model of the material parameters presented in this work and based on the transmission-line model of Omega particle was contradictory (particle size was taken of the order $`\lambda `$ but the constitutive relations were introduced as if the composites were continuous media) and violated some physical foundations (Kramers-Kronig relations, etc).

In the following we consider the effective properties of isotropic composites of both Omega and SRR particles embedded in a dielectric host medium. We study two different kinds of particles arrangement. The first one is a mixture with a random orientations of particles. The second one is obtained by positioning the particles at the sides of a cubic unit cell, like it is shown in Figs 2 and 3. Different inclusions and structures are finally compared.

## 2 Theory

Both SRR particles and Omega particles are bianisotropic particles. In other words, an external electric field produces a magnetic dipole in the particle and, at the same time, an external magnetic field produces an electric dipole. For SRR, this magnetoelectric effect is reduced by choosing the opposite positioning of the splits of both loops. However, it cannot cancel out exactly since the loops are of different sizes.

Let a SRR particle be located in an infinite and homogeneous dielectric medium with permittivity $`ϵ`$. Consider the electromagnetic response of an individual SRR particle to local magnetic and electric fields. The polarizabilities of a SRR particle are dyadic coefficients relating the local electric and magnetic fields $`𝐄^{\mathrm{loc}},𝐇^{\mathrm{loc}}`$ with the induced electric and magnetic dipole moments:

$`𝐩`$ $`=`$ $`\overline{\overline{a}}_{ee}𝐄^{\mathrm{loc}}+\overline{\overline{a}}_{em}𝐇^{\mathrm{loc}}`$ (1)
$`𝐦`$ $`=`$ $`\overline{\overline{a}}_{me}𝐄^{\mathrm{loc}}+\overline{\overline{a}}_{mm}𝐇^{\mathrm{loc}}`$ (2)
The electric and magnetic moments of the SRR as a whole are vector sum of the moments induced in the outer and inner loops (indices 1 and 2, respectively) :

$`𝐩`$ $`=`$ $`𝐩_1+𝐩_2`$ (3)
$`𝐦`$ $`=`$ $`𝐦_1+𝐦_2`$ (4)
In this section we evaluate analytically the polarizabilities of a SRR particle, assuming that the inclusion is fabricated from a round wire of radius $`r_0`$. The section of the wire is small enough to satisfy the inequality $`r_0a_{1,2}`$, where $`a_{1,2}`$ are the averaged radiuses of rings 1 and 2 (see Fig. 1a). Distance $`\delta `$ between the broken ends of each ring is assumed to be small compared to $`a_{1,2}`$. The analytical model of the Omega particle will not be reproduced here. Analytical antenna theory and its verification can be found in .

Finally, using the calculated polarizabilities we derive for each kind of particles the material parameters through the well-known Maxwell Garnett model.

### 2.1  SRR particle eigenfrequencies

We consider the SRR particle as a pair of two mutually coupled resonant scatterers. The complex amplitudes of induced currents at the splits of both rings can be written as follows:

$`I_1(Z_{C1}+Z_{L1})+I_2Z_{12}`$ $`=`$ $`_{H1}+_{E1},`$ (5)
$`I_2(Z_{C2}+Z_{L2})+I_1Z_{12}`$ $`=`$ $`_{H2}+_{E2}`$ (6)
Indices 1 and 2 refer to the outer and inner ring respectively. $`_H`$ and $`_E`$ are the electromotive forces induced by the external fields. $`I_{1,2}`$ are the currents at the splits of the rings. $`Z_{L1,L2}`$ represent the inductive loop impedances ($`Z_{L1,L2}=j\omega L_{1,2}`$), and $`Z_{C1,C2}`$ are the capacitive impedance ($`Z_{C1,C2}=1/(j\omega C_{1,2})`$) of the slits. For a magnetic field excitation we can neglect the individual loop capacitance which has only a small influence. For the study of electric polarization, influence of the proper capacitance of one loop must be taken into account. For that we use the known theory of a split loop antenna  and  which takes into account the capacitive properties of an individual loop scatterer and its radiation losses. $`Z_{12}`$ in (5),(6) is the mutual impedance of broken loops. We have chosen to include directly the capacitive mutual coupling in the term $`C_{1,2}`$. Otherwise the analytical model of mutual coupling would become too involved. With these considerations the mutual impedance $`Z_{12}`$ in (5),(6) is expressed only as a mutual inductance $`M`$ ($`Z_{12}=j\omega M`$).

The capacitances $`C_{1,2}`$ are composed of two contributions. The first part $`C_S`$ originates from the splits of each loop. The second one is due to capacitive coupling of loops and can be seen as a consequence of the following effect. Consider first as a reference the conductor of one loop (the outer loop of a SRR particle, for example). Next, place another conductor parallel to the first one (in our case, the inner loop, with the split at the opposite side). Then, assume that one half of the first loop (e.g. the upper part) contains positive charge $`+Q`$, and the second half of the first loop (e.g. the lower part) contains negative charge $`Q`$. These charges are distributed somehow along the loop perimeter. As the two loops are positioned closely, in the upper half of the second loop, a negative charge $`q`$ is induced. Hence, the distribution of the global charge $`q`$ around the second loop repeats the distribution of $`+Q`$ in the first loop. Note that we do not have necessarily the correspondence $`Q=q`$. Of course, the same thing happens in the lower part of the structure, where the induced charge $`+q`$ appears. So two equivalent capacitances $`C_{mut}`$ arise from this situation as shown in Fig. 4. Hence, we obtain an effective additional capacitance which is equal to $`C_{mut}/2`$, connected in parallel to the proper capacitance of the split. One can note that this equivalent scheme can be improved taking into account that there is an additional capacitance of the second slit. However, we are trying to derive a model of SRR as simple as possible and the complexity, increased by this additional capacitance, is not justified.

Values of inductances $`L_{1,2}`$, mutual inductance $`M`$ and total capacitances $`C_{1,2}`$ are analytically calculated in Appendix.

To obtain the eigenfrequencies of the SRR particle, we let the right-hand side of (5),(6) vanish and equate the determinant to zero:

$$(Z_{C1}+Z_{L1})(Z_{C2}+Z_{L21})Z_{12}^2=0$$
(7)
Equation (7) can be written as

$$\left(\frac{1}{j\omega C_1}+j\omega L_1+r_1\right)\left(\frac{1}{j\omega C_2}+j\omega L_2+r_2\right)+\omega ^2M^2=0$$
(8)
where $`r_{1,2}`$ are small resistances that represent the radiation and the Joule losses of both loops. They lead to a small imaginary correction term denoted as $`\kappa (\omega )`$ in equation (9):

$$(\omega ^2\omega _{}^2)(\omega ^2\omega _+^2)+j\kappa (\omega )=0$$
(9)
If we assume that the SRR particle is perfectly conducting, the term $`\kappa `$ describes only the radiation losses of a magnetic dipole and it can be derived from the energy balance condition, as it was done in  for arbitrary dipole scatterers. Comparing (8) with (9) we find the eigenfrequencies:

$$2\left(\frac{\omega _\pm }{\omega _m}\right)^2=\frac{\omega _m^2}{\omega _1^2}+\frac{\omega _m^2}{\omega _2^2}\pm \sqrt{\left(\frac{\omega _m^2}{\omega _1^2}+\frac{\omega _m^2}{\omega _2^2}\right)^24}$$
(10)
where we denote

$$\omega _m^2=\frac{1}{\sqrt{C_1C_2(L_1L_2M^2)}},\omega _{1,2}^2=\frac{1}{C_{1,2}L_{1,2}}$$

Frequencies (10) correspond to two eigenmodes. If the frequencies $`\omega _{}`$ and $`\omega _+`$ differ dramatically the two eigenmodes of the induced current are excited at different frequencies. This leads to two separate resonances (that magnetic and that electric ones) of the particle. However, we will see below that in practical cases $`\omega _+\omega _{}`$. The two resonance bands overlap, and we obtain a result close to that of the Lorentz model with only one resonance, as it goes for simple scatterers. As a consequence, it matches qualitatively the results of .

### 2.2  Magnetic polarizability of SRR particle

Next, we are going to derive the polarizabilities of SRR particle using (5),(6) and definitions (1) and (2). First, we study the magnetic polarizability. For a magnetic excitation ($`𝐄^{\mathrm{loc}}=0,𝐇^{\mathrm{loc}}=𝐳_0H^{\mathrm{loc}}`$, where the $`z`$-axis is perpendicular to the loop plane as it is shown in Fig. 2) we obtain:

$$_{H1,2}=j\omega \mu _0S_{1,2}H^{\mathrm{loc}},m_{1,2}=I_{1,2}\mu _0S_{1,2}$$
(11)
Here $`S_{1,2}=\pi a_{1,2}^2`$ are the effective loop areas.

The only component of the magnetic moment is vertical, i.e.

$$\overline{\overline{a}}_{mm}=𝐳_0𝐳_0a_{mm}^{zz}$$

Then, relations (11) together with (1)-(6) and (7)-(9) lead to the following expression:

$$a_{mm}^{zz}=D\omega ^2\frac{B_1(\omega ^2\omega _1^2)+B_2(\omega ^2\omega _2^2)}{\mathrm{\Delta }+j\kappa }$$
(12)
where we have denoted

$$\mathrm{\Delta }=(\omega ^2\omega _{}^2)(\omega ^2\omega _+^2)$$
(13)
$$B_1=S_1(S_1L_2MS_2),B_2=S_2(S_2L_1MS_1)$$

$$D=\frac{\mu _0^2}{L_1L_2M^2}$$

and

$$\omega _{1,2}^2=\frac{\omega _{1,2}^2}{\left(1\frac{MS_{2,1}}{L_{2,1}S_{1,2}}\right)}$$
(14)
Now let us find the term $`\kappa `$. The theory from  based on the balance of received and re-radiated power for any particle with magnetic dipole polarizability $`a_{mm}^{zz}`$ gives the equation

$$\mathrm{Im}\frac{1}{a_{mm}^{zz}}=\frac{k^3}{6\pi \mu _0}$$
(15)
Comparing this formula with (12) we come to the following expression for the term $`\kappa `$:

$$\kappa =\frac{\omega ^2Dk^3\left[B_1(\omega ^2\omega _1^2)+B_2(\omega ^2\omega _2^2)\right]}{6\pi \mu _0}$$
(16)
With an appropriate calculation of inductances and capacitances (see Appendix), it is now possible to completely evaluate the magnetic polarizability $`a_{mm}^{zz}`$.

### 2.3  Electric polarizability of SRR particle

Now we are interested in the electric polarizability. Consider the action of local electric fields on a SRR. Two possibilities exist to induce polarization on the SRR particle.

First, we study the case where the electric field is along the horizontal direction ($`x`$direction for an individual particle). The only electro-dipole mode of the current is excited in both loops in this case. The induced charges are concentrated at the splits and at the opposite points and the polarizability component $`a_{ee}^{xx}`$ is quasi-static. It does not depend on the presence or absence of the split and equations (5),(6) are useless in this case111 Because $`I_{1,2}=0`$ and the electromotive force $`_E`$ induced at the split is also equal to zero: charges are the same at the edges of each split..

In  we derived formula (28) for $`xx`$component of the polarizability of a circular loop with area $`S=\pi a^2`$. Denoting this polarizability as $`a_{QS}`$ we obtain from :

$$a_{QS}=\frac{4SJ_1^{}(ka)}{\omega \eta A_1(ka)}$$

In this equation, $`\eta =\sqrt{\mu _0/ϵϵ_0}`$ is the wave impedance of the host medium, $`k=\omega \sqrt{\mu _0ϵϵ_0}`$ is the wave number in the medium, $`J_1^{}(ka)`$ is the derivative of the Bessel function of the first order, $`A_1(ka)`$ is one of so-called King’s coefficient known from the theory of loop antennas ,. King’s coefficients with a very high accuracy are approximated in  and take the following form:

$$A_1(ka)=\left(\frac{ka\frac{1}{ka}}{\pi }\right)\left(\mathrm{log}\frac{8a}{r_0}2\right)\frac{0.667(ka)^30.207(ka)^5}{\pi }$$

$$j0.333(ka)^2j(0.133ka)^4+j0.026(ka)^6$$

$$A_0(ka)=\frac{ka}{\pi }\left(\mathrm{log}\frac{8a}{r_0}2\right)+\frac{0.667(ka)^30.267(ka)^5}{\pi }j0.167(ka)^4j0.033(ka)^6$$

$$A_2(ka)=\left(\frac{ka\frac{4}{ka}}{\pi }\right)\left(\mathrm{log}\frac{8a}{r_0}0.667\right)+\frac{0.4ka+0.21(ka)^30.086(ka)^5}{\pi }$$

$$j0.05(ka)^4j0.012(ka)^6$$

Since we have two loops, we must sum up their $`a_{QS}`$ to obtain the result:

$$a_{ee}^{xx}=\frac{4S_1J_1^{}(ka_1)}{\omega \eta A_1(ka_1)}\frac{4S_2J_1^{}(ka_2)}{\omega \eta A_1(ka_2)}$$
(17)
Now we consider the case where the electric field is directed along $`y`$axis. Here the resonant electric polarization appears in both loops of the SRR particle. To obtain the resulting polarizability, we combine the theory of a single loop with resonant electric polarization (see  and ) with equations (5) and (6).

In  the following formula has been derived for the $`yy`$component of the electric polarizability of a single loaded loop of radius $`a`$:

$$a_{ee}^{yy}=\frac{4\pi a^2J_1^{}(ka)}{\omega \eta A_1(ka)}\left(1+\frac{2j}{\pi \eta A_1(ka)}\frac{1}{Y_L+Y_{\mathrm{split}}}\right)$$
(18)
Here $`Y_{\mathrm{split}}`$ is the admittance that is connected to the loop split. For both loops we have $`Y_{\mathrm{split}}1/Z_{\mathrm{split}}=j\omega C_{1,2}`$ which include the mutual capacitance of the loops. $`Y_L=1/Z_L`$ is the proper loop admittance. In our case of coupled loops, we should replace $`Y_L`$ by the value $`Y_{1,2}`$ which takes into account their mutual inductance. For loop $`1`$, we can say that it is composed of the proper impedance of the loop $`Z_{L1}`$ and of the additional impedance $`Z_{12}I_2/I_1`$, which is the contribution of loop $`2`$ into loop $`1`$. Then, for loop $`1`$, we have to substitute into (18):

$$\frac{1}{Y_1}=Z_1=Z_{L1}+Z_{12}\frac{I_2}{I_1}$$

Similarly, for the second loop we have:

$$\frac{1}{Y_2}=Z_2=Z_{L2}+Z_{12}\frac{I_1}{I_2}$$

Denoting the ratio $`I_1/I_2`$ as $`\xi `$ and taking into account that $`Z_{12}=j\omega M`$, we come to the following formulas for the $`yy`$polarizabilities of both loops:

$$a_{e1}^{yy}=\frac{4S_1J_1^{}(ka_1)}{\omega \eta A_1(ka_1)}\left(1+\frac{2j}{\pi \eta A_1(ka_1)}\frac{1}{Y_1+j\omega C_1}\right)$$
(19)
$$a_{e2}^{yy}=\frac{4S_2J_1^{}(ka_2)}{\omega \eta A_1(ka_2)}\left(1+\frac{2j}{\pi \eta A_1(ka_2)}\frac{1}{Y_2+j\omega C_2}\right)$$
(20)
where

$$Y_1=\frac{Y_{L1}}{1+\frac{j\omega MY_{L1}}{\xi }}$$

$$Y_2=\frac{Y_{L2}}{1+j\omega MY_{L2}\xi }$$

The proper loop admittance is given by :

$$Y_L=\frac{1}{j\pi \eta }\left(\frac{1}{A_0(ka)}+\frac{2}{A_1(ka)}+\frac{2}{A_2(ka)}\right)$$
(21)
Expressions of $`Y_{L1,L2}`$ follow from (21) with substitutions $`a=a_{1,2}`$.

Equations (5),(6) allow to find $`\xi `$:

$$\xi =\frac{_{E1}(Z_{L2}+Z_{C2})_{E2}Z_{12}}{_{E2}(Z_{L1}+Z_{C1})_{E1}Z_{12}}$$
(22)
If a loop with proper impedance $`Z_L`$ loaded with the impedance of the split $`Z_{\mathrm{split}}`$ is submitted to an external electric field, an electromotive force is induced. Its value is given by formula (26) from :

$$_E=\frac{4aJ_1^{}(ka)(\frac{1}{Y_L}+\frac{1}{Y_{\mathrm{split}}})}{j\eta A_1(ka)}\left(1+\frac{j}{\pi \eta }\frac{1}{Y_L+Y_{\mathrm{split}}}\right)E_y^{\mathrm{loc}}$$
(23)
Substituting the values $`a_{1,2}`$ and $`r_0`$ into (21) and (23) allows to calculate the electromotive forces induced by an external electric field in each loop of the SRR particle. Then, substituting expressions of $`_{E1,E2}`$ into (22) with the help of relation (21) we can find $`\xi `$. Consequently, we derive an expression of $`Y_{1,2}`$ which is then substituted into (19) and (20). Finally, we evaluate $`a_{ee}^{yy}`$ as the sum

$$a_{ee}^{yy}=a_{e1}^{yy}+a_{e2}^{yy}$$
(24)
At this point, we have at our disposal a complete analytical model of the electric and magnetic polarizabilities of an individual SRR particle (other components of $`\overline{\overline{a}}_{ee}`$ and $`\overline{\overline{a}}_{mm}`$ are zeros).

### 2.4 Effective polarizabilities of a planar bianisotropic particle in isotropic mixtures

An isotropic mixture with a random orientation of particles can be considered as a mixture realised with isotropic particles which have averaged scalar polarizabilities :

$$a_{ee}^{av}=\underset{\alpha =1}{\overset{3}{}}\frac{a_{ee}^{\alpha \alpha }}{3},a_{mm}^{av}=\underset{\alpha =1}{\overset{3}{}}\frac{a_{mm}^{\alpha \alpha }}{3},a_{em}^{av}=a_{me}^{av}=0$$
(25)
Here we have designated Cartesian components with index $`\alpha =x,y,z`$. The particles under consideration do not have cross components: $`a_{ee}^{\alpha \beta }=0`$ for $`\alpha \beta `$. Moreover, we have seen with our analytical model, that the only polarizabilities that are non-zero are: $`a_{ee}^{xx},a_{ee}^{yy}`$ and $`a_{mm}^{zz}`$.

If SRR or Omega particles are arranged on the sides of cubic unit cells (see Figs. 2 and 3), the magnetoelectric effect also disappears and dyadics $`\overline{\overline{a}}_{em},\overline{\overline{a}}_{me}`$ of the individual particles do not contribute into material parameters. Furthermore, mutual coupling of six particles belonging to the same cell is weak enough if the cell size $`d_{cell}`$ is a bit more than the particle size $`d_p`$ (the particles do not touch one another). Now, let us evaluate the equivalent effective polarizability of each kind of unit cell ignoring the mutual coupling between the particles of a cubic cell.

#### 2.4.1 Effective polarizabilities of a cubic cell of SRR particles

The cubic cell under consideration is presented Fig. 2. An electric field directed along an arbitrary edge of a cubic cell excites the electric dipoles in four SRR particles. For two particles $`𝐄^{\mathrm{loc}}`$ is parallel to the line connecting the rings splits and corresponds to the polarizability $`a_{ee}^{xx}`$. For the other two particles $`𝐄^{\mathrm{loc}}`$ is perpendicular to this line and corresponds to $`a_{ee}^{yy}`$.

Consider now the magnetic field effect on the cell. A magnetic field directed along an arbitrary edge of a cell induces two equivalent magnetic moments in two opposite particles. So, the magnetic polarizability of a unit cell turns out to be the double of $`a_{mm}^{zz}`$.

To summarize, the cubic unit cell is equivalent to an isotropic particle with electric and magnetic polarizabilities given by

$$a_{ee}^{}=2a_{ee}^{xx}+2a_{ee}^{yy},a_{mm}^{}=2a_{mm}^{zz}$$
(26)
#### 2.4.2 Effective polarizabilities of a cubic cell of Omega particles

For a unit cubic cell of Omega particles (see Fig. 3) the situation is quite different. Here the resonant excitation by the electric field directed along the arms of $`\mathrm{\Omega }`$ (along $`y`$-axis in the model of an individual particle) is mainly due to the presence of the arms . Therefore the total electric polarization of two opposite omega particles keeps resonant behavior. Within the resonant band the quasi-static component of Omega particle polarizability $`a_{ee}^{xx\mathrm{\Omega }}`$ describing the response of a particle to the electric field applied normally to the arms (see Fig. 1) is small compared to $`a_{ee}^{yy\mathrm{\Omega }}`$. Then one obtains for a unit cell from Omega particles approximate relations

$$a_{ee}^{}2a_{ee}^{yy\mathrm{\Omega }},a_{mm}^{}=2a_{mm}^{zz\mathrm{\Omega }}$$
(27)
where $`a_{ee}^{yy\mathrm{\Omega }},a_{mm}^{zz\mathrm{\Omega }}`$ can be taken from .

#### 2.4.3 Material parameters of composites

The material parameters of the isotropic mixture are given by the standard Maxwell Garnett formulas :

$$ϵ_{eff}=ϵ+F\left(\frac{Na_{ee}^{}}{ϵ_0}\frac{N^2a_{ee}^{}a_{mm}^{}}{3ϵ_0\mu _0}\right)$$
(28)
$$\mu _{eff}=1+F\left(\frac{Na_{mm}^{}}{\mu _0}\frac{N^2a_{ee}^{}a_{mm}^{}}{3ϵ_0\mu _0ϵ}\right)$$
(29)
In these equations

$$F^1=1\frac{Na_{ee}^{}}{3ϵ_0ϵ}\frac{Na_{mm}^{}}{3\mu _0}+\frac{N^2a_{ee}^{}a_{mm}^{}}{9ϵ_0\mu _0ϵ}$$

These equations are general for the different kinds of composites we are studying. Depending on the composition of mixture, only expressions for polarizabilities change.

* for random media, which are realized by mixing particles in an isotropic and random way in a host medium, $`N`$ denotes the concentration of particles, and $`a_{ee}^{}=a_{ee}^{av}`$, $`a_{mm}^{}=a_{mm}^{av}`$ are given by (25)
* for cubic-cell media, which are obtained with previously mentioned cubic unit cells, $`N`$ denotes the concentration of cells, and $`a_{ee}^{}`$ and $`a_{mm}^{}`$ are given by (26) or (27).

## 3 Model validation for SRR particle

Our analytical model is below compared with numerical simulations. Numerical model which allows to simulate the averaged polarizabilities of an arbitrary bianisotropic particle with thin wire of round cross section, that allows to calculate also material parameters of a random composite has been derived from , . Since its publication this method has been widely used to simulate bianisotropic composites.

We compared the analytical and numerical models of the SRR prepared with wires of radius $`r_0=0.2`$ mm. Numerical example presented here corresponds to the following parameters. The averaged radiuses of the inner and outer loops are equal to $`a_2=2.1`$ and $`a_1=2.7`$ mm. The length of the broken part is taken $`\delta =0.1`$ mm or $`\delta =0.4`$ mm ($`\delta `$ is the same in both loops). The host medium permittivity $`ϵ=1.5`$.

In Fig. 5 the real part of $`a_{mm}^{av}`$ is shown calculated in both analytical and numerical ways for the case $`\delta =0.1`$ mm. The imaginary parts of $`a_{mm}^{av}`$ for this case are shown in Fig. 6. The real and imaginary parts of $`a_{ee}^{av}`$ calculated in both ways are shown in Figs. 7 and 8, respectively.

Comparing the resonant frequencies of $`a_{ee}^{av}`$ and $`a_{mm}^{av}`$, we can note that the difference is quite negligible. This means that the resonant bands overlap ($`\omega _+\omega _{}`$) as we mentioned previously.

We can observe also, that the analytical model shows a weak additional resonance for the electric polarization around 5.2 GHz. This is a defect of our analytical theory which is related to the adopted approximations. We think that the second resonance follows from our model of loops mutual coupling which is not accurate enough. Adding a very small imaginary value to the right-hand side of formula (31) we completely suppress this second resonance without changing the values of $`a_{mm}`$ and $`a_{ee}`$ within the first resonant band. Of course, this correction has nothing to do with the theory but it confirms that the model of $`Z_{1,2}`$ is not accurate enough and can be improved.

The electric polarizability calculated in both analytical and numerical ways for the case $`\delta =0.4`$ mm (the other parameters are the same as above) is presented in Fig. 9.

The agreement between two models is quite the same for $`\delta =0.1`$ and $`\delta =0.4`$ mm. The same remarks could be made for $`a_{mm}`$.

Polarizabilities of SRR particle are very sensitive to deviations of the particle parameters, and we can consider that the agreement between the two models is satisfactory. It makes our analytical theory suitable to calculate the effective composite parameters and, hence, allows to solve explicitly the inverse problem for such composites.

## 4 SRR and Omega isotropic composites

In this section we analyze the material parameters of effective media obtained with SRR particles or Omega particles. The mixture is fabricated by

* distributing the particles with random orientations in a host medium
* positioning the particles at the sides of a cubic unit cell (filled by the same medium as the matrix) as shown in Fig. 2 and 3.

### 4.1 SRR composites

In Figs. 10 and 11 both $`ϵ_{eff}`$ and $`\mu _{eff}`$ of a random SRR composite are shown (calculated in both analytical and numerical ways). Particles geometry and background permittivity has been taken the same as above (the split width is $`\delta =0.4`$ mm). The concentration of particles has been taken $`N=1`$ particle/cm3. For the particle maximal size $`5.8`$ mm this concentration allows to use the Maxwell Garnett averaging procedure.

We can observe on Fig. 10 that the permittivity is not becoming negative within the resonant band (whereas both models of the permeability allow the negative values for $`\mathrm{Re}(\mu _{eff})`$). Here, the permittivity of SRR composite remains positive because of the restriction of dilute mixtures (imposed by the validity of the Maxwell Garnett model).

In the case of SRR cubic media, since we apply the Maxwell Garnett model to the unit cells considering them as particles, we can take the concentration of SRR particles greater than in the case of a random mixture, and that, without violating the restriction of dilute mixtures.

As a consequence, in Fig. 12 we can see that $`\mathrm{Re}(\mu _{eff})`$ and $`\mathrm{Re}(ϵ_{eff})`$ are negative simultaneously. The concentration of cells has been calculated for the case of unit cells arrangement as presented in Fig. 3. The size of the cell is equal to $`6\times 6\times 6`$ mm3, and the averaged distance between the centers of two adjacent cells is close to $`11`$ mm. This allows to consider the composite with such cells as a dilute mixture.

### 4.2 Omega composites

In Fig. 13 and 14 we present the results calculated within the frame of the Maxwell Garnett model for an Omega composite prepared from cubic units as shown in Fig. 3. In that case, the maximal size of the Omega particle is equal to $`6.4`$ mm. The cell size is equal to $`6.5\times 6.5\times 6.5`$ mm3. For an illustration, we have chosen the following dimensions. The radius of the loop is $`3`$ mm and the radius of the wire is $`r_0=0.1`$ mm. The length of an arm is $`L=3`$ mm. The host medium permittivity is $`ϵ=2`$. Note that for cases when $`\delta a`$, parameter $`\delta `$ does not influence the polarizabilities at all (see ).

In Fig. 13, the real parts of the permittivity and permeability of the Omega composite versus the frequency are shown. In Fig. 14, the imaginary parts of permittivity and permeability are presented. We can see that $`\mathrm{Re}(ϵ_{eff})`$ and $`\mathrm{Re}(\mu _{eff})`$ reach $`1`$ practically at the same frequency. This is a very exciting result if we refer to paper  where the perfect lens has been predicted for the case $`Re(ϵ_{eff})=Re(\mu _{eff})=1`$.

The case of random Omega composite (Fig. 15) with density $`N=1`$ particle/cm3 only differs by the heights of the resonant peaks of the permittivity.

## 5 Conclusion

Based on the analysis of electromagnetic phenomena that appear in a split ring resonator, we have developed an analytical model of SRR particles. We have shown that the SRR particles can exhibit, on their own, simultaneously magnetic and electric polarizability. The proposed model has been successfully compared with a numerical approach. So, these results point out that previous works on SRR do not give a satisfactory explanation of the behavior of this particle. Next, we have calculated with the Maxwell Garnett approach the effective parameters of composites realized with SRR particles. We have shown that it is possible to obtain an isotropic composite with negative refraction choosing properly the geometry, arrangement and concentration of SRR particles in a host dielectric matrix. The way of making these composites has been presented. With this technique, we have obtained an effective material with simultaneously negative real part of permittivity and permeability.

However, we have shown that SRR is not the only and not the best choice for obtaining a backward-wave medium at microwaves. Actually, we can obtain the same kind of medium, using Omega particle instead of SRR particles. With these inclusions we have derived a material which gives $`\mathrm{Re}(ϵ_{eff})=\mathrm{Re}(\mu _{eff})=1`$ in the resonant band. This is the wished result exposed by the inventor of the SRR particle. Finally, we can conclude that the Omega particle is a more appropriate choice in order to create microwave composites with negative refraction.

We have to point out yet, that our model exhibits a high level of radiation losses within the resonant band. As a consequence, if material samples are so lossy, it would be problematic to measure the negative refraction on it. These high radiation losses are due to the Maxwell Garnett model which implies a random organization of particles and is consistent with the Kramers Kronig relations. However, the proposed cubic composite can be classified in a middle position between random media and lattices. The composite is constructed as a lattice, stacking up cubes in a tree dimensional chessboard assembly, but it is not a perfect lattice. There is indeed, a certain random organization in the positions of the particles on the sides of the cubic cells (as it can be seen on figures 2 and 3). So, we can reasonably expect that losses given by the Maxwell Garnett model will be overestimated.

## 6 Appendix

Inductance of a round loop of radius $`R`$ made with a wire of cross section radius $`r_0`$ can be found e.g. in :

$$L=\mu _0R\left(\mathrm{log}\left(\frac{8R}{r_0}\right)2\right)$$
(30)
Inductances of both rings of SRR can be obtained with (30) with the substitution $`R=a_{1,2}`$.

Formula (5.23) from  for two concentric loops whose radii $`a_1,a_2`$ are close to one another, gives the following approximation for the mutual inductance:

$$M=\mu _0a_1\left[(1\chi )\mathrm{log}\left(\frac{4}{\xi }\right)2+\xi \right]$$
(31)
where $`\chi =(a_1a_2)/2a_1`$ is a small parameter ($`0<\chi <0.25`$). In (31) we have neglected the terms of order $`\xi ^2`$.

Now, consider the capacitances included in the loops splits. We can write in accordance with Fig. 4

$$C_1=C_{s1}+\frac{C_{mut}}{2},C_2=C_{s2}+\frac{C_{mut}}{2}$$
(32)
For $`C_s`$ we have the simple equation of a conventional (parallel-plate) capacitor

$$C_s=\frac{\pi r_0^2ϵϵ_0}{\delta }$$
(33)
if the split width $`\delta `$ is much smaller than the wire cross section diameter. Otherwise, if $`\delta r_0`$ one has $`C_{s1,2}C_{mut}`$ and the capacitance of the split can be ignored. We do not study the cases when $`r_0\delta `$.

Next, we calculate $`C_{mut}`$. The capacitance per unit length of two parallel wires of diameter $`2r_0`$ with the distance between the wire axes $`d=a_2a_1`$ denoted below as $`C_0`$ is known :

$$C_0=\frac{\pi ϵϵ_0}{\mathrm{arccosh}\left(\frac{d}{2r_0}\right)}$$
(34)
Roughly (but successfully) we can consider $`C_{mut}`$ as the capacitance between two equivalent parallel straight wires of length $`\pi a\pi (a_1+a_2)/2`$ with distance $`d`$ between their axes. However, we would obtain a bad result if we simply multiply $`C_0`$ by $`\pi a`$ in order to find $`C_{mut}`$. We should take into account the non-uniformity of the charge distribution along these equivalent wires and then write $`C_{mut}=C_0P`$ (where the effective length of two parallel charged wires $`P`$ is smaller than $`\pi a`$ due to this non-uniformity).

We use the King theory of loop antennas  in order to find effective length $`P`$. For it we write the Fourier expansion for the current induced in the reference loop by a magnetic field (with respect to the azimuth angle $`\varphi `$):

$$I(\varphi )=I^{(0)}+I^{(1)}\mathrm{cos}\varphi +I^{(2)}\mathrm{sin}\varphi +\mathrm{}$$
(35)
It corresponds to the following expansion for the induced charges per unit angle:

$$\tau (\varphi )=\tau _1\mathrm{sin}\varphi +\tau _2\mathrm{cos}\varphi +\mathrm{}$$
(36)
In the King theory it was proven that (for the case of small loop antenna ($`ka_{1,2}1`$) with a finite impedance of the split) the omitted terms in expansions (35) and (36) are of the next order of smallness compared to the remained terms . The uniform term $`I^{(0)}`$ represents the magneto-dipole mode of the induced current and has nothing to do with induced charges. Two other terms in (35) are responsible for the electric polarization of the loop and are related with the charge distribution. We can rewrite (36) as:

$$\tau (\varphi )=\tau _{max}\mathrm{sin}(\varphi +\varphi _0)$$
(37)
where $`\tau _{max}\sqrt{\tau _1^2+\tau _2^2}`$. So, at it was already mentioned, half of the loop ($`\varphi _0<\varphi <\pi +\varphi _0`$) contains positive charge $`Q`$, and the other half ($`\pi +\varphi _0<\varphi <\varphi _0`$) contains negative charge $`Q`$. The distribution of these charges is sinusoidal. The induced charges in the other loop have the opposite sign but the same sinusoidal distribution (see also Fig. 4).

We obtain that the distribution of the charge along two effective equivalent parallel wires of length $`\pi a`$ is sinusoidal with zeros at the ends of these wires. It leads to the result $`P=2a`$ which gives together with (34) the approximate relation for the capacitance $`C_{mut}`$:

$$C_{mut}=\frac{2\pi ϵϵ_0a}{\mathrm{arccosh}\left(\frac{d}{2r_0}\right)}$$
(38)
If the loop is excited by $`y`$polarized electric field, the formula of the current distribution (35) keeps valid . Coefficients $`I^{(0,1,2)}`$ change in this case, but the charges $`+Q`$ and $`Q`$ remain sinusoidally distributed on the two halves of each the loop. Therefore the formula (38) keeps valid, too.


text:
# Constraints on the cosmic neutrino background

## 1 Introduction

The quest for cosmological parameters is the major goal of cosmology, and we are experiencing an exciting era in which newer and better data allow for the precise determination of an increased number of parameters. In this paper we focus on one specific parameter, the amount of relativistic energy density at recombination/equivalence, and we discuss the implications of the current CMB and LSS observations. The relativistic energy density is often parametrized by the equivalent number of standard model neutrino species $`N_{eff}`$. The standard model predicts three neutrino species, plus a correction imposed by neutrinos not being completely decoupled during electron-positron annihilation ($`\mathrm{\Delta }N_{eff}=0.03`$) (Dolgov et al. 1997, Gnedin & Gnedin 1998) and by the finite temperature QED correction to the electromagnetic plasma ($`\mathrm{\Delta }N_{eff}=0.01`$) (Lopez et al. 1999, Mangano et al. 2002).

Departures from the standard model which imply a variation in $`N_{eff}`$ may be due to decaying dark matter (Bonometto & Pierpaoli 1998, Hannestad 1998, Lopez et al 1998, Kaplinghat & Turner 2001 or quintessence (Bean et al 2001). In these scenarios $`N_{eff}`$ can be either smaller or greater than 3.04, and the constraints on $`N_{eff}`$ implied by CMB and LSS need not to agree with the ones implied by Big Bang Nucleosyntesis (Kaplinghat & Turner 2001). We will analyze and discuss the CMB+LSS results on $`N_{eff}`$. Unlike previous work on this subject, we allow for a cosmology with general curvature. We will use the CMB data from the WMAP satellite (Bennet et al. 2003), and complement them with the CBI ones (Pearson et al. 2002). As for the LSS, we use the 2dF data (Percival et al 2002). We use a general cosmology in our analysis because at the present time we do not have an independent probe of the flatness of the Universe other than the CMB. It is therefore fictitious and misleading to restrict ourselves to flat cosmologies in the so called “precision cosmology” era. We will show that allowing for non–zero curvature changes quite significantly the results on $`N_{eff}`$, despite the fact that the inferred constraints on curvature are still compatible with a flat Universe.

The paper is organized as follows: in section 2 we discuss the effects of $`N_{eff}`$ on the matter and radiation power spectra, in section 3 we discuss the CMB+LSS constraints, while section (4) is dedicated to conclusions.

## 2  Effects of a relativistic background on the matter and radiation power spectra

An excess of energy in relativistic particles with respect to the one predicted by the standard model has an impact on both the radiation and the matter power spectrum, for different reasons. The increased relativistic energy delays equivalence, and scale entering the horizon at equivalence is bigger. Since the amplitude of dark matter perturbations on smaller scales is frozen during the radiation dominated era, the matter power spectrum turnover is shifted toward lower $`k`$’s. The effect can be quantified via the shape parameter $`\mathrm{\Gamma }\mathrm{\Omega }_mh(g_{}/3.36)^{1/2}`$ (White, Gelmini & Silk 1995), where $`g_{}`$ represents the relativistic degrees of freedom ($`g_{}=2+0.454N_{eff}`$). For a given power at large scales, the small scale power s therefore reduced, and a lower value of $`\sigma _8`$ is implied.

An increased radiation density implies a smaller conformal time and therefore a smaller sound horizon. As a consequence, the peaks of the CMB power spectrum, which reflects the peaks/valleys of the oscillations at recombination, are shifted toward smaller scales (bigger $`l`$’s) (Pierpaoli & Bnometto 1999). In addition, since recombination occurs when the Universe is not completely matter dominated, the early Integrated Sachs-Wolfe (ISW) effect causes an enhanced first peak.

A description of parameter degeneracies in a flat Universe is given by Bowen et al. (2002), where they find that $`N_{eff}`$ is mainly degenerate with $`\mathrm{\Omega }_m`$ and the spectral index $`n`$. They show that almost complete degeneracy can be obtained by keeping $`z_{eq}`$, $`\mathrm{\Omega }_bh^2`$ and $`R`$ (the position of the first acoustic peak with respect to a reference model) fixed while varying $`N_{eff}`$ produces almost degenerate power spectra. Here we extend the analysis to the real data and consider general curvature.

Allowing for possible curvature affects the radiation power spectrum. In particular, the position of the peaks is strongly dependent on curvature. A given comoving scale at the last scattering surface subtends an angle on the sky that is curvature–dependent. If the Universe is open ($`\mathrm{\Omega }_k>0`$) \[ or closed ($`\mathrm{\Omega }_k<0`$)\] then the angle associated with a particular scale at last scattering is smaller \[larger\] than the corresponding flat Universe one. Consequently, the peaks/valleys of the CMB power spectrum will be located at larger \[smaller\] $`l`$’s. A slightly close Universe can compensate the shift of the peaks toward higher $`l`$’s implied by an increased radiation background. We therefore expect some degeneracy between curvature and $`N_{eff}`$, and if we aim to determine both from the CMB measurements, we ought to treat both of them as free parameters in the data analysis.

## 3  CMB and LSS constraints

We used a modified version of the COSMOMC (Lewis & Bridle 2002) package to compute the likelihoods also including the WMAP new results (Verde et al. 2003). We considered only CDM adiabatic perturbations, and the following set of parameters, with usual definitions: $`\mathrm{\Omega }_{dm}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_k`$, $`H_0`$, $`N_\nu `$, $`n_s`$, $`A_s`$, $`z_{re}`$ (reionization redshift). In our analysis we use the WMAP and CBI data together with the 2dF power spectrum, and we discuss the effect of adding other priors.

First we compare the results for a flat model ($`\mathrm{\Omega }_k=0`$) with those for a general Universe when $`H_0`$ to vary only in the range $`64H_080`$ corresponding to the 1$`\sigma `$ allowed range of the HST result (Freedman et al. 2001). In fig.1 we present the likelihoods for each parameter, after marginalization over all the others, for flat and generally curved cases. For the number of effective neutrinos we find a best fit $`N_{eff}=4.08`$ ($`N_{eff}=3.70`$) for the general curvature (flat) case, with $`2.23N_{eff}6.13`$ ($`1.82N_{eff}5.74`$) at the 95% C.L. The assumption of a flat Universe biases the result toward low $`N_{eff}`$ and shrinks the inferred error bars. Our results are in agreement with what found by Crotty et al. (2003) with the assumption of flat geometry. Notice that the general curvature analysis is important at the present time because of the great improvement that the WMAP results have implied on the CMB power spectrum. Constraints from CMB and LSS on flat Universes previous to the WMAP were much weaker: $`N_{eff}15`$ at 95% C.L. (Hannestad 2001, see also Hansen et al., 2002).

In fig.1 the flat and curve case present the similar upper limit for $`N_{eff}`$ mainly because we imposed a tight top-hat prior on $`H_0`$ ($`<80`$). The likelihood indeed seems to prefer high values for $`H_0`$, which may allow for extra radiative component.

For the general curvature case, we now proceed to discuss the effect of various priors on the determination of $`N_{eff}`$.

In fig. 2 we present the results obtained allowing for a wider $`H_0`$ (and $`z_{re}`$) range. The best fit for $`N_{eff}`$ in this case is $`N_{eff}=4.31`$, with a 95 per cent C.L. range : $`1.6N_{eff}7.1`$ from CMB and 2dF only, and $`N_{eff}=4.08`$ with $`1.9N_{eff}6.62`$ when the $`H_0`$ prior from the HST project is included.

As for the other parameters, we note that the inclusion of extra relativistic energy causes a higher $`n_s`$ values than in the standard case (Spergel et al 2003), a higher $`z_{re}`$ and a slightly closed Universe ($`\mathrm{\Omega }_k=0.013\pm 0.015`$). We argue that, since a high $`N_{eff}`$ boosts the first peak through the early ISW effect, a higher $`n_s`$ combined with an early reionization of the Universe can still ensure a good fit to the CMB data. Matter power spectrum data on smaller scales than the ones probed by 2dF (e.g. Ly-$`\alpha `$ forest data) may be used in breaking this degeneracy, because the matter power spectrum is only sensitive to $`n_s`$ and not to $`z_{re}`$. Data of the matter power spectrum at very small scales would be incredibly sensitive to small increases in $`n_s`$. However, there is still much debate on small scale data, since their interpretation may be complicated by the non–linear growth of the fluctuations. We adopted here a conservative approach and chose not to include them in the analysis. In fig. 2 we note that the derived value of $`\sigma _8`$ is quite high ($`0.97\pm 0.1`$) yet because of the favoured high $`n_s`$ values. Actual quotations of the $`\sigma _8`$ value range between 0.75 and 1 ( Pierpaoli et al. 2003, and references therein), so that a broad prior on $`\sigma _8`$ would already decrease the allowed $`n_s`$ and improve the errors on $`N_{eff}`$. Consistent 5 per cent accuracy measurements of $`\sigma _8`$ from both weak lensing and clusters would greatly improve the constraint on $`N_{eff}`$. As an example, in fig. 2 we plot the curves obtained with an hypothetical prior on $`\sigma _8`$ with the typical scaling derived from cluster and weak lensing: $`\sigma _8(\mathrm{\Omega }_m/0.3)^{0.6}=0.85\pm 0.05`$. Such prior would imply a lower $`\mathrm{\Omega }_m`$ and $`N_{eff}=3.8\pm 1.6`$.

As for the other parameters, we note that $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are well constrained in the range: $`\mathrm{\Omega }_m=0.29\pm 0.06`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.72\pm 0.05`$ (1 $`\sigma `$ error).

Hannestad (2003) while analyzing flat Universes uses a tight prior on $`\mathrm{\Omega }_m`$ derived from the Supernova Cosmology Project (Perlmutter e al. 1999). As a consequence he derives a small $`N_{eff}`$ value with small error bars. We run a separate chain treating CMB+LSS+SN data in a general cosmology with COSMOMC, and found that the inclusion of the Supernovae in the analysis doesn’t particularly improve the results on $`N_{eff}`$ (see fig. 2). Including the SN in the analysis slightly improves the $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ determination, but the errors on $`N_{eff}`$ are practically unchanged. We conclude that the strong constraints obtained by Hannestad (2003) are mainly due to the restriction to a flat Cosmology and to the artificially small prior assumed on $`\mathrm{\Omega }_m`$ as a way of introducing the SN constraint.

In table 1 we plot the correlation matrix for the parameters under consideration. Note that $`N_\nu `$ is mostly correlated with $`\mathrm{\Omega }_{dm}h^2`$ and $`\mathrm{\Omega }_k`$. A better independent estimate of the redshift of equivalence (typically probing $`\mathrm{\Omega }_mh`$), combined with a better independent measurement of $`H_0`$, may help in breaking the degeneracy.

## 4 Conclusions

The new CMB data (WMAP and CBI) together with the matter power spectrum derived by the 2dF galaxy survey can constrain the effective number of neutrino species much more precisely than previous experiments. Previous estimates of $`N_{eff}`$ were derived under the assumption of nul curvature. Since we don’t have any independent confirmation of the flatness of the Universe other than the CMB itself, we argue that the curvature should be kept as a free parameter in the estimation of $`N_{eff}`$. We compare the results derived from the two different hypothesis. Applying a top hat prior for $`H_0`$ ($`64H_080`$), we find for a flat Universe $`1.82N_{eff}5.74`$ at 95 per cent C.L., with a best fit of $`N_{eff}=3.70`$, while with general curvature $`2.23N_{eff}6.13`$ with a best fit $`N_{eff}=4.08`$ . Allowing for general curvature shifts the acceptance range of $`N_{eff}`$ toward higher values mainly because the curvature tends to compensate the effect of $`N_{eff}`$ on the peak locations in the CMB power spectrum.

In the case of general curvature, we explored a wider range in $`H_0`$ and $`z_{re}`$ and applied different priors. We find $`1.6N_{eff}7.1`$ (best fit $`N_{eff}=4.31`$) at 95% C.L. from CMB and 2dF only, and $`N_{eff}=4.08`$ with $`1.9N_{eff}6.62`$ when the $`H_0`$ prior from the HST project is included as a proper Gaussian prior. No significant modifications derive from the inclusion of the SN constraint.

By looking at the likelihood distribution for each parameter after marginalization over all the others we conclude that the inclusion of an extra relativistic component would suggest a higher expansion rate $`H_0`$, a higher spectral index $`n_s`$ and $`\mathrm{\Omega }_k`$ slightly negative, if compared to the standard analysis with three neutrinos (Spergel et al. 2003).

We analyze the correlations between the various parameters and conclude that $`N_{eff}`$ is most degenerate with $`\mathrm{\Omega }_{dm}h^2`$ and $`\mathrm{\Omega }_k`$. We argue that other independent measurements of the matter power spectrum, like precise determinations of $`\sigma _8`$ from clusters and lensing or probes at smaller scales from the Ly–$`\alpha `$ forest, would help in constraining the epoch of equivalence and therefore would improve the results on $`N_{eff}`$. Moreover, it would greatly improve the constraints on the large $`n_s`$ now allowed.

## 5 Acknowledgments

EP is supported by NASA grant NAG5-11489. The author thanks Antony Lewis for assistance with the use of COSMOMC; Rachel Bean and Olivier Doré for useful discussions, and the anonymous referee for useful comments.


text:
# Weak Universality of Spin Glasses in Three Dimensions

## 1 Introduction

A spin glass (SG) phenomenon has been attracting great interests for thirty years both theoretically and experimentally sgreview . It is also one of the toughest subject in a field of the computational physics. Though many theoretical investigations have been done, it is still an unsettled issue whether or not a simple random bond model explains the SG phase transition. Spins of many SG materials are well-approximated by the Heisenberg spins. However, numerical studies have suggested that there is no finite-temperature SG transition in the Heisenberg model mcmillan ; olive . This is quite different from the Ising SG case: many numerical investigations maricampbell  consistently explain corresponding experimental results isingsgexp . Kawamura chirality  proposed the chirality mechanism in order to settle the discrepancy in the Heisenberg spin glass. The chiral-glass transition occurs without the spin-glass order. A small but finite random anisotropy in the real materials mixes the chirality and the spin. This effect induces the observed SG transition. The scenario sounds nice. However, Matsubara matsubara1 ; matsubara2 ; matsubara3  recently recalculated the domain-wall excess energy and the spin-glass susceptibility and suggested that the finite SG transition occurs. The methods are quite similar to the previous ones mcmillan ; olive . The final conclusion drastically changes by a subtle difference in the analyses of the obtained data.

Computer simulations on SG models suffer from serious slow dynamics. The system sizes treated up to now are accordingly limited to very small ones, e.g., mostly a linear size $`L=20`$ or less in three dimensions. The contradiction of the numerical results mentioned above possibly originates in the smallness of the system. This is a limit of the conventional equilibrium simulational study. One realizes the equilibrium state in finite size systems. Then, the thermodynamic properties are obtained by the finite-size scaling analysis.

We take an opposite approach to the thermodynamic limit and reexamine the SG transition. First, we take the infinite size limit by dealing with a very large system within a finite time range before the finite-size effect appears. Then, the finite-time scaling analysis is performed to obtain the thermodynamic properties. This method is known as the nonequilibrium relaxation (NER) method ner ; huse ; blundell . In our previous paper totasg1  it is clarified that the SG transition occurs at the same finite temperature as the chiral-glass transition occurs. Of course, the chirality freezes if the spin freezes. The estimated critical exponent $`\gamma `$ is consistent with the corresponding experimental result heisenexp . The chirality mechanism is not necessary. In this article we expand the same analysis to the Ising SG model and the XY SG model in three dimensions. It is found that the Ising, the XY, and the Heisenberg models undergo finite-temperature SG transitions with the weak universality.

## 2 Method

The model we work with is the nearest-neighbor $`\pm J`$ random bond model. For each random bond configuration we prepare eight real replicas with different paramagnetic spin configurations and with different random number sequences. Spins are updated by the single-spin-flip algorithm. We calculate at each time step $`t`$ the spin-glass susceptibility $`\chi _{\mathrm{sg}}`$ through overlaps $`q_{\mu ,\nu }^{\alpha \beta }`$ between two replicas ($`\alpha `$ and $`\beta `$) out of eight:

$$\chi _{\mathrm{sg}}=\frac{1}{N}\underset{i,j}{}\left[𝑺_i𝑺_j_{\mathrm{therm}}^2\right]_{\mathrm{conf}}=N\left[\frac{2}{m(m1)}\underset{\alpha >\beta }{\overset{m}{}}\underset{\mu ,\nu }{\overset{x,y,z}{}}(q_{\mu ,\nu }^{\alpha \beta })^2\right]_{\mathrm{conf}}$$

Indices $`\mu `$ and $`\nu `$ denote three components of spin degrees of freedom. The thermal average is replaced by the summation over the real replicas, whose number is denoted by $`m`$. It determines the statistics of the thermal average. Thus, a relaxation function $`\chi _{\mathrm{sg}}(t)`$ for one bond configuration is obtained. A configurational average is taken over $`\chi _{\mathrm{sg}}(t)`$ with different random bond configurations, different initial paramagnetic states, and different random number sequences.

The spin-glass susceptibility is expected to diverge at the transition temperature ($`T_{\mathrm{sg}}`$) as $`\chi _{\mathrm{sg}}(t)t^{\gamma /z\nu }`$, where $`z`$ is the dynamic exponent. We obtain $`T_{\mathrm{sg}}`$, $`\gamma `$, and $`z\nu `$ by the finite-time scaling analysis. The NER of the Binder parameter $`g(t)`$ is also calculated. Since the quantity is related to the fourth-order cumulant, quite a lot of samples are necessary to obtain meaningful data. Numbers of the bond configurations to obtain the results in this article is summarized in Table 1. The Binder parameter is expected to behave at $`T_{\mathrm{sg}}`$ as $`g(t)t^{d/z}`$, by which $`z`$ is independently obtained. Then, $`\nu `$ is estimated from a value of $`z\nu `$ obtained by the scaling. The NER of the spin-glass susceptibility was originally applied to the Ising SG model by Huse huse  and that of the Binder parameter was applied by Blundell et al. blundell .

## 3 Results

Numerical results suggesting the weak universality is shown in Fig. 1. The NER data of the spin-glass susceptibility exhibit an algebraic divergence with the exponent $`\gamma /z\nu =0.38`$. The NER function of the Ising model and that of the Heisenberg model are not distinguishable. That of the XY model differs by a factor three. An amplitude of the $`\chi _{\mathrm{sg}}`$ is dependent on the model. If we take into account a correction-to-scaling term with its exponent $`w/z=0.5`$ maricampbell , the relaxation functions fit to the expression $`At^{\gamma /z\nu }[1Bt^{w/z}]`$ from the first step. A deviation of the XY data from the diverging line after $`t=10^5`$ is possibly due to the finite-size effect. The NER of the Binder parameter shows an algebraic divergence with $`z=6.2`$ for both models. Therefore, we obtain $`\gamma /\nu =2.4`$ which is common to the Ising and the Heisenberg model. This suggests the weak universality.

Results on the Heisenberg model are briefly reported in totasg1 . We have obtained $`T_{\mathrm{sg}}=0.20(2)`$, $`\gamma =1.9(5)`$, and $`z\nu =4.7(12)`$. Since $`z=6.2`$, $`\nu `$ is estimated as $`\nu =0.8(2)`$.

An example of a finite-time scaling analysis of the Ising model is shown in Fig. 2. This is an only preliminary result with a small system size ($`L=19`$) and a short time range. The error bars are consequently large. Nevertheless, the obtained results are consistent with previous numerical investigations maricampbell  and the corresponding experiment isingsgexp  as summarized in Table 2. The exponent estimated by the scaling ($`\gamma /z\nu =0.39`$) is consistent with that of the raw data at $`T_{\mathrm{sg}}`$ as shown in Fig. 1 ($`\gamma /z\nu =0.38`$).

Time evolution of the distribution function of the overlap $`P(q,t)`$ is investigated at $`T=T_{\mathrm{sg}}=1.17`$ (figures not shown). It exhibits a single Gaussian shape with a peak at $`q=0`$ before the size effect appears. As the time increases, the width of the Gaussian distribution grows in accordance with the divergence of the spin-glass susceptibility. It is possible to scale $`P(q,t)/t^{\gamma /2z\nu }`$ plotted versus $`qt^{\gamma /2z\nu }`$ for various time steps from $`t=10`$ to $`t=10^4`$. That is, the NER function knows the critical phenomenon from the 10th step. The distribution changes its shape to having double peaks at $`\pm q_{\mathrm{eq}}`$ after the finite-size effect appears.

It is also possible to investigate the low-temperature phase. The data suggest that the $`\chi _{\mathrm{sg}}`$ diverges algebraically with an exponent linearly dependent on the temperature. Details will be reported elsewheretotaisingsg .

## 4 Concluding Remarks

By applying the NER method it is clarified that the $`\pm J`$ spin glass models in three dimensions undergo finite-temperature SG transitions. There is possibly the weak universality with a ratio of the critical exponents $`\gamma /\nu =2.4`$. We compare our results with other numerical results and the experimental results in Table 2. They agree well within the numerical errors. Critical exponents $`\gamma `$ and $`\nu `$ in the Heisenberg model are almost one half of those in the Ising model. Reexamination of the experiment on the critical exponents is strongly encouraged.

The NER method is particularly effective in slow-dynamic systems. Applications to various complex systems should be fruitfultotaner . It potentially bring a breakthrough to a system stagnant due to the difficulty in computer simulations.

###### Acknowledgements.

The authors thank F. Matsubara for fruitful discussions, and N. Ito and Y. Kanada for providing them with a random number generator. Computations were partly done at the Supercomputer Center, ISSP, The University of Tokyo. The author TN thanks a financial support from The Sumitomo Foundation.