text:
# Dielectric photonic crystal as medium with negative electric permittivity and magnetic permeability

## Abstract

We show that a two-dimensional photonic crystal (PC) made from a non-magnetic dielectric is a left-handed material in the sense defined by Veselago. Namely, it has negative values of both the electric permittivity $`ϵ`$ and the magnetic permeability $`\mu `$ in some frequency range. This follows from a recently proven general theorem. The negative values of $`ϵ`$ and $`\mu `$ are found by a numerical simulation. Using these values we demonstrate the Veselago lens, a unique optical device predicted by Veselago. An approximate analytical theory is proposed to calculate the values of $`ϵ`$ and $`\mu `$ from the PC band structure. It gives the results that are close to those obtained by the numerical simulation. The theory explains how a non-zero magnetization arises in a non-magnetic PC.

About 35 years ago Victor VeselagoVeselago (1967) considered theoretically propagation of the electromagnetic waves (EMW’s) in a hypothetical left-handed medium (LHM), where both the electric permittivity $`ϵ`$ and the magnetic permeability $`\mu `$ are negative in some frequency range. Since the square of the speed of light $`c^2=c_0^2/ϵ\mu `$ is positive, the EMW’s propagate but they have unusual properties, like anomalous Doppler and Cherenkov effects, negative light pressure, and negative refraction (NR) at the interface with a regular medium. The latter property enables three-dimensional (3D) imaging. The origin of all these properties is the fact that the energy flux, S and the wave vector, k are antiparallel in the LHM. An important assumption of Veselago’s theory is the absence of spatial dispersion, which means that both $`ϵ`$ and $`\mu `$ are $`𝐤`$-independent. However both of them are functions of frequency. A generalization of Veselago’s theory for the case of spatial dispersion would not be simple because in this case the expression for the Poynting vector includes an extra term Landau and Lifshitz (1960).

To create the LHM Veselago proposed to combine two independent different subsystems Veselago (1967), one with a negative $`ϵ`$ and another with a negative $`\mu `$. As an example, he considered a ferromagnetic metal at frequency close to the ferromagnetic resonance and below the plasma frequency of the metal. Following this line, the San Diego groupShelby et al. (2001) has observed a NR in an artificial composite system and attributed this effect to a two-dimensional (2D) LHM. Their system consists of two elements: 2D metallic photonic crystal (PC) creating negative $`ϵ`$ below a cutoff frequency and the so-called split-ring resonators creating negative $`\mu `$ due to anomalous dispersion near the resonance. The experiments of the San Diego group are very important because they have ignited the whole field, started by Veselago about 35 years ago but completely forgotten since then. However they caused new hot debates concerning different aspects of their experiments and interpretationsPokrovsky and Efros (2002a); Valanju et al. (2002); Smith et al. (2002a, b); Garcia and Nieto-Vesperinas (2002).

Many groupsKosaka et al. (1998); Notomi (2000); Luo et al. (2002a, b); Cubukcu et al. (2003) predicted and observed NR in PC’s using for the prediction the dispersion relation of EMW’s and methods of geometrical optics. However, the “effective refractive index” they found was, in general, a function of $`𝐤`$. We think that it originates from the spatial dispersion of the dielectric constant for EMW with $`𝐤`$ deep in the Brillouin zone. Several authorsKosaka et al. (1998); Notomi (2000); Luo et al. (2002b) found that near the $`\mathrm{\Gamma }`$ point this index becomes $`k`$-independent. As we show below, together with the negative group velocity this indeed is a manifestation of the LHM.

The main idea of this work is based upon the theorem proved in Ref. Pokrovsky and Efros (2002b). It reads that the group velocity $`𝐯_g=\omega /𝐤`$ in an isotropic medium is positive in the regular media (RM)($`ϵ>0`$, $`\mu >0`$) and negative in the LHM. This means that vectors $`\omega /𝐤`$ and $`𝐤`$ are parallel in the RM and antiparallel in the LHM. This theorem follows from the thermodynamics requirement that the energy of electromagnetic fields is positive.

Now we generalize this theorem to a 2D PC illustrated in the left bottom inset of Fig. 1. The PC represents a square lattice of period $`d`$, comprising cylindrical holes in a dielectric medium with the axes of cylinders along $`z`$-direction. We consider propagation of EMW’s with the electric field $`𝐄`$ in this direction (s-wave). The calculated weu  low frequency spectrum of EMW’s in this PC is shown in Fig. 1 for the case when the dielectric material of the crystal is characterized by $`ϵ_m=12`$, $`\mu _m=1`$. The fields in the PC are represented by the Bloch functions with a 2D wave vector $`𝐤`$=($`k_x`$, $`k_y`$). One can see from Fig. 1 that in the frequency range a little below the top of the band $`1`$ at the $`\mathrm{\Gamma }`$ point, the group velocity $`\omega /𝐤`$ is negative and there are no other bands in this frequency range. It is also important that band $`1`$ is not degenerate at $`\mathrm{\Gamma }`$-point. Near the $`\mathrm{\Gamma }`$-point the wavelength is much larger than the lattice period of the PC and one can introduce $`𝐤`$-independent $`ϵ_{ik}`$ and $`\mu _{ik}`$. Since $`z`$-axis is the axis of $`C_4`$ symmetry, we have a uniaxial crystal and both tensors in the principal axes have only two components $`ϵ_{zz}`$, $`ϵ_{}=ϵ_{xx}=ϵ_{yy}`$ and $`\mu _{zz}`$, $`\mu _{}=\mu _{xx}=\mu _{yy}`$ respectively. The electromagnetic energy density of the s-wave can be represented in the formLandau and Lifshitz (1960)

$$\overline{U}=\frac{1}{16\pi \mu _{}\omega }\frac{d[\omega ^2n^2]}{d\omega }|E|^2,$$
(1)
where $`n^2=ϵ_{zz}\mu _{}`$. In the rest of the paper we use notation $`ϵ_{zz}=ϵ`$, $`\mu _{}=\mu `$. Since $`\overline{U}>0`$ we get $`\mu (𝐯_g𝐤)>0`$, where 2D group velocity $`𝐯_g=2c_0^2𝐤/[d(\omega ^2n^2)/d\omega ]`$. Thus, negative $`𝐯_g𝐤`$ means negative $`\mu `$ and vise versa. For the propagating waves it also means $`ϵ<0`$. It follows that this PC is a LHM in the frequency range under study.

Note, that the p-wave (magnetic field in $`z`$ direction) is characterized by $`ϵ_{}`$ and $`\mu _{zz}`$. These components might also be negative in the vicinity of the $`\mathrm{\Gamma }`$ point but within a different frequency range.

We find $`ϵ`$ and $`\mu `$ by microscopic calculations and show that they are indeed both negative in the frequency range under study. For this purpose we consider a PC slab infinite in $`z`$-direction and embedded in a homogeneous medium with $`ϵ^{}`$, $`\mu ^{}`$. A 2D point source (an infinite dipole along $`z`$-axis) is located in the homogeneous medium to the left of the slab. If we now choose $`ϵ^{}`$ and $`\mu ^{}`$ so that the cylindrical wave in the system is the same as it would be in a completely homogeneous medium with $`ϵ^{}`$ and $`\mu ^{}`$ (though inside the slab there is a microscopic field), the values of $`ϵ^{}`$ and $`\mu ^{}`$ obtained in this way are the macroscopic $`ϵ`$ and $`\mu `$ of the slab at a given frequency. The best choice for $`ϵ^{}`$ and $`\mu ^{}`$ at the working frequency $`\omega =0.3(2\pi c_0/d)`$ is obtained by a fitting procedure and the result is shown in Fig. 2. It yields $`ϵ=1.2`$, $`\mu =0.096`$. Thus, at this frequency the PC slab is a LHM and this result is in accordance with the above theorem. Note, that in the case of the LHM the vectors S and k are in the opposite directions (see animation http://www.physics.utah.edu/ efros ).

One can calculate independently the refractive index of the slab $`n^2=ϵ\mu `$ from the dispersion relation of the band $`1`$. At small $`k`$ it can be written in a form $`\omega ^2=\omega _1^2\alpha c_0^2k^2`$. Our calculation of spectrum gives $`\alpha =0.94`$, $`\omega _1=0.314(2\pi c_0/d)`$. Using the alternative definition of $`n^2`$ namely, $`\omega ^2n^2=c_0^2k^2`$ one gets

$$n^2=\frac{1}{\alpha }\left(1+\frac{\omega _1^2}{\omega ^2}\right).$$
(2)
Note, that $`n0`$ as $`\omega \omega _1`$. At the working frequency $`n^2=(0.32)^2`$, which agrees with $`ϵ`$ and $`\mu `$ found above.

If we choose $`ϵ^{}=|ϵ|`$ and $`\mu ^{}=|\mu |`$ as the parameters of the homogeneous medium surrounding the PC, we obtain a unique device proposed by VeselagoVeselago (1967) that we call the Veselago lens. It is shown schematically in Fig. 3a. The computer simulation of the Veselago lens is shown in Fig. 4 (also see animation http://www.physics.utah.edu/ efros ). In the above simulations the PC is cut in the direction $`[10]`$. We have checked that the results do not change if the PC is cut in the direction $`[11]`$. This confirms the absence of spatial dispersion at the working frequency. The size of the focus is of the order of the wavelength as predicted by the diffraction theory Pokrovsky and Efros (2002c).

Now we explain how the starting theorem that follows from Eq.(1) and is formulated in terms of macroscopic electrodynamics can be interpreted microscopically. The negative group velocity shown at Fig. 1 results from the repulsion of bands $`1`$ and $`2`$, $`3`$. Therefore a microscopic explanation of negative $`\mu `$ should be found in terms of the interaction between these bands. We propose a quantitative explanation by writing Maxwell’s equation for the microscopic field $`e_z(𝐫)`$ in a matrix formJoannopoulos et al. (1995); Sakoda (2001) in the basis of Bloch functions at $`k=0`$. Taking into account that near the $`\mathrm{\Gamma }`$-point the bands $`1`$, $`2`$, $`3`$, $`4`$ are far from the other bands, we truncate the matrix keeping only these four bands. This approximation is analogous to the Kane model in semiconductor physicsKane (1957). By diagonalizing the truncated matrix we find the spectra of all four bands and the correct combinations of basis functions for each band at finite $`k`$. We have checked that the dispersion law for all bands $`1`$-$`4`$ as obtained from the truncated model does not differ by more than $`4\%`$ from the result presented in Fig. 1 at $`kd/2\pi <0.1`$. At small $`k`$ the microscopic electric field for the band $`1`$ has the form

$$e_z(𝐫)=e^{i𝐤𝐫}\left(u_1\frac{ik_0k_xc_0^2}{\omega _2^2\omega _1^2}u_2\frac{ik_0k_yc_0^2}{\omega _2^2\omega _1^2}u_3\frac{ik_0k_1k^2c_0^4}{(\omega _2^2\omega _1^2)(\omega _4^2\omega _1^2)}u_4\right),$$
(3)
where $`𝐫`$ is the 2D-vector, $`u_1`$, $`u_2`$, $`u_3`$ and $`u_4`$ are the Bloch functions at $`k=0`$, the positive parameters $`k_0`$, $`k_1`$ have the dimensionality of a wave vector and they are integrals from these functions. The symmetry of the function $`u_1`$ is $`X^2Y^2`$. The function $`u_4`$ has $`X^2+Y^2`$ symmetry. The functions $`u_2`$ and $`u_3`$ have the symmetries $`X`$ and $`Y`$ respectively. The parity of $`X`$ and $`Y`$ functions within the elementary cell is shown in the central inset of Fig. 1.

If the material of the PC were metallic, the electric field in the bands $`2`$ and $`3`$ would create two currents in opposite directions. These currents would induce a magnetic moment in $`x`$ and $`y`$ directions. In case of the dielectric PC these currents are imaginary displacement currents and they do not create any magnetic moment in the modes $`2`$ and $`3`$. However, as follows from Eq.(3) the admixture of bands $`2`$ and $`3`$ to band $`1`$ is also imaginary. Therefore the band $`1`$ has a real magnetic moment.

The magnetization can be calculated using the expression

$$𝐌=\frac{1}{2c_0d^2}𝐫\times \frac{𝐏}{t}𝑑s,$$
(4)
where $`𝐏=𝐳e_z(𝐫)[ϵ(𝐫)1]/4\pi `$, $`𝐳`$ is a unit vector in $`z`$ direction, $`ϵ(𝐫)=1`$ inside the holes and $`ϵ_m`$ otherwise, the integration is performed over the unit cell, and $`e_z(𝐫)`$ is given by Eq.(3). Only the second and the third terms from Eq.(3) contribute to the magnetization. The result can be expressed through the integrals from $`u_2`$ and $`u_3`$. These functions have been taken from the computation of the spectrum shown at Fig. 1. The macroscopic vectors $`𝐁`$ and $`𝐄`$ are connected by the relation $`𝐁=c_0𝐤\times 𝐄/\omega `$, where $`𝐄=𝐳<e_z(𝐫)>`$ and $`<\mathrm{}>`$ means averaging over the unit cell. Only the function $`u_4`$ gives a non-zero contribution after the averaging. It gives the main contribution to $`𝐄`$. Other basis functions give non-zero contribution only when multiplied by $`𝐫`$-dependent terms in the expansion of the $`\mathrm{exp}(i𝐤𝐫)`$ within the unit cell. Their contributions to $`𝐄`$ are of the same order in $`k`$ but they are much smaller than the contribution from $`u_4`$ since $`(\omega _4\omega _1)/\omega _11`$, where $`\omega _1`$ and $`\omega _4`$ are the frequencies of the bands $`1`$ and $`4`$ at $`k=0`$. In fact, the latter strong inequality permits truncation of the matrix in our analytical model.

Now we can find $`\mu `$ from the definition $`\mu ^1=14\pi |𝐌|/|𝐁|`$. The result is

$$\mu =0.89(1\omega _1^2/\omega ^2).$$
(5)
Comparing Eq.(2) and Eq.(5) one can see that the root of $`\mu `$ at $`\omega =\omega _1`$ is responsible for the root of $`n^2`$. This means that $`ϵ`$ is non-zero and negative at $`\omega =\omega _1`$. On the other hand, $`\mu `$ changes its sign at $`\omega =\omega _1`$. Thus, the gap at $`\omega >\omega _1`$ is due to $`\mu ϵ<0`$. At $`\omega =0.3(2\pi c_0/d)`$ one gets $`\mu =0.084`$ and that is very close to the value obtained above by the fitting procedure. Thus, we have found analytically the value of $`\mu `$.

The 3D LHM can be obtained in the same way, for example, using the PC with the cubic symmetry proposed by Luo et al.Luo et al. (2002b). This PC also has a negative group velocity near the $`\mathrm{\Gamma }`$ point with no other bands at this frequency.

Thus, we have shown that the LHM can be readily obtained using only one dielectric PC that creates both negative $`ϵ`$ and $`\mu `$. We have presented the physical explanation of this phenomenon in a 2D case. The working frequency can be easily changed in a wide range by changing the lattice constant and $`ϵ_m`$ of the PC. It is important for applications that the losses in a dielectric PC are much smaller than in a metallic one. We have found not only $`n`$ but also $`ϵ`$ and $`\mu `$ that are necessary to describe the dynamical properties of the interfaces. This is especially important for the creation of the Veselago lens, which requires no reflection at interfaces, thus imposing two matching conditions on $`ϵ`$ and $`\mu `$. By matching $`n`$ only one gets a modification of the Veselago lens with multiple foci (Fig. 3b)Pokrovsky and Efros (2002d). We realize, that we have not presented a method of matching the LHM with the surrounding medium. This matching is very easy in computational approach but it must be handled creatively on the way toward real devices.

###### Acknowledgements.

We are grateful to Serge Luryi and Valy Vardeny for helpful discussions. The work has been funded by the NSF grant DMR-0102964.


format: text

text:
# Theory of fusion hindrance and synthesis of the superheavy elements

## 1 INTRODUCTION

Since the discovery of the periodicity in chemical elements by Mendelejeff, our knowledge has been expanded. The heaviest element in nature is known to be Uranium with the atomic number Z=92. Heavier elements than that have been synthesized artificially. Thus, it is an intriguing question how many elements can exist or what is the heaviest element we can synthesize. Larger $`Z`$ values result in the instability due to the larger Coulomb repulsion which would dominate over nuclear attraction. But it is also well-known that closed shells of nucleonic structure give rise to an extra-binding which contributes to the stabilization of atomic nucleus. Many attempts have been made to predict double closed shell nuclei heavier than 208Pb. The magic number for proton next to 82 is predicted to be 114, 120 or 126, depending on the nuclear models employed, while that for neutron next to 126 is commonly predicted to be 184. These results suggest that there is a stable region in the nuclear chart far away from that of the known isotopes, which is sometimes called an island of the superheavy elements (SHE). Naturally, enormous experimental efforts have been devoted to answer to the question last few decades. Until now, the elements with up to Z=112 have been synthesized by heavy ion fusion experiments with 208Pb target, and about even heavier elements with Z=114 and 116 are reported indicative experimental observations with 48Ca beams on the actinide targets. But reaction mechanisms of nuclear fusion of massive systems are not well known, and thus the experiments have been performed following the systematics of the experimental data available so far. Therefore, it is a very important and urgent subject to develop a theory which permits us to predict which combination of projectile and target is favorable for synthesis of an element under consideration, and at what incident energy residue cross sections are optimized. Of course, in order to predict residue cross sections quantitatively, we have to treat not only fusion processes, but also cooling processes precisely, because compound nuclei formed by fusion reactions are excited. More explicitly, the residue cross section is given by the following formula, presuming the statistical theory of the compound nucleus,

$$\sigma _{\text{res}}(E_{\text{c.m.}})=\pi \lambda _{}^{}{}_{}{}^{2}\mathrm{\Sigma }_J(2J+1)P_{\text{fusion}}^J(E_{\text{c.m.}})P_{\text{surv}}^J(E^{}),$$
(1)
where $`E^{}`$ is the excitation energy of the compound nucleus and is equal to $`E_{\text{c.m.}}+Q`$ with $`Q`$ being the $`Q`$-value of the reaction, and $`\lambda ^{}`$ is the wave length divided by $`2\pi `$. The symbol $`J`$ is the quantum number of total angular momentum of the system, as usual. $`P_{\text{fusion}}^J`$ and $`P_{\text{surv}}^J`$ denote the fusion and the survival probabilities, respectively. The latter is the probability for surviving against fission decay and charged particle emissions, which are well described as statistical decays. A new computer program KEWPIE (Kyoto Evaporation Width calculation Program with tIme Evolution) is constructed so as to minimize the ambiguities, which is briefly described later. The remaining unknown is the nuclear masses of the superheavy isotopes, or in other words, their shell correction energies, to which $`xn`$ residue cross sections are extremely sensitive.

On the other hand, as stated above, there is no commonly accepted theory for the fusion probability, though it is well-known experimentally that there is a strong hindrance in fusion reactions of massive systems. That is, it does not reach 1/2 at the incident energy equal to an expected barrier height, and in order for that, a large amount of energy is required additionally, so-called extra-push energy. Two possible explanations have been proposed. One is due to the loss of the incident energy by friction which is supposed to act between the ions of the entrance channel. This is inferred from the studies on so-called Deep-Inelastic Collisions (DIC). Gross and Kalinowski proposed surface friction model (SFM) which successfully explained the characteristic aspects of DIC. The other one is due to dissipation of energy of the collective motion which leads the system to the spherical compound nucleus, starting from the pear-shaped configuration formed by the sticking of the ions of the entrance channel. This is conceivable, because there is the strong dissipation in the fissioning motion, i.e., in an “inverse” process. This is confirmed from the studies of anomalous multiplicities of neutrons, gamma rays, etc. emitted prior to fission of excited nuclei, with the strong friction given by the one-body wall-and window formula (OBM). Thus, it is orthodox and comprehensive to consider that both mechanisms are working in the fusion reactions of massive systems. As schematically shown in Fig. 1, two processes undergo in time-sequence, i.e., collision processes are up to the contact point of the projectile and the target to stick each other, and subsequently from the sticked configuration, shape evolution of the united system starts toward the spherical shape over the conditional saddle point. Therefore, the fusion process is described by the two steps, and then the fusion probability $`P_{\text{fusion}}^J`$ is given by a product of the sticking probability $`P_{\text{stick}}^J`$ and the formation probability $`P_{\text{form}}^J`$;

$$P_{\text{fusion}}^J(E_{\text{c.m.}})=P_{\text{stick}}^J(E_{\text{c.m.}})P_{\text{form}}^J(E_{\text{c.m.}})$$
(2)
This is a new two step model. The method of connection is not like that from a diabatic to an adiabatic treatments, but should be called “statistical”, because the first step generally results in statistical distributions of physical quantities of the united system and they should be used as initial conditions for the second step, as will be detailed below.

In order to realize the treatment of the first step, SFM is employed, but in an extended form so as to include the fluctuation forces associated with the friction forces in accord with the dissipation-fluctuation theorem. As stated above, the model is used up to the contact point, not for the whole processes like in DIC. Then, the results of the calculations provide the sticking probability $`P_{\text{stick}}^J`$ for the projectile and the target to stick each other, as well as information on kinetic energy of the radial motion, i.e., information on the distribution of the radial momentum at the contact point. Since they were already reported in detail elsewhere, we do not repeat it here, but summarize the results.

1. An energy dependence of the sticking probability shows that a barrier height is effectively shifted to an energy about 10 MeV higher than the original one, which would already provide an explanation of the extra-push energy.
2. A distribution of the radial momentum at the contact point has almost exactly a Gaussian form whose center, i.e., the mean value is equal to zero. Its variance is consistent with the temperature calculated from the energy conservation in average.
3. The orbital angular momentum is found to reach the dissipation limit.

These results indicate the fact that the projectile and the target stick each other at the contact point to form a united system. Therefore, with this distribution as the initial values, we solve shape evolution of the united mono-nuclear system to obtain the formation probability, i.e., the probability that the system overcomes the saddle point to reach the spherical shape, by the use of the same type of Langevin equation which was used for description of fission of excited nuclei.

Recently it has been shown analytically that the dissipation manifests its effects in an extremely enhanced way in over-saddle-point problem with a schematic one-dimensional model. This is a modeling of the fusion reactions of massive systems where the sticked configurations are located outside of the conditional saddle point and thus the system has to overcome it in order to reach the spherical shape. Furthermore, we can obtain a simple expression for the extra-push energy. Since they have been reported elsewhere, in the present paper we concentrate ourselves on realistic calculations of the formation and then the fusion probabilities by the use of the liquid drop model (LDM). In the next section, we describe dynamical evolution of nuclear shapes over the ridge line with Langevin equations for the collective variables. In section 3, we discuss about 48Ca- induced fusion reactions on the actinide targets, on which experimental fusion excitations are available for comparisons. By the use of KEWPIE, we calculate the survival probability, which is combined with the fusion probability to give $`xn`$ residue cross sections for SHE with Z=114, 116 and 118. The results are compared with the experimental data available for the first two elements.

## 2 FORMATION PROBABILITY OF THE SPHERICAL COMPOUND NUCLEUS

Shapes of the united nucleus are described by the two-center parameterization (TCP) with the distance between two mass-centers, the mass-asymmetry, and the neck parameter, while the deformations of fragments are neglected. In view of the extremely strong friction for the neck degree of freedom on the basis of OBM, the neck motion is expected to be very slow, which permits us to freeze it in a good approximation. Thus, the parameter $`ϵ`$ is taken to be 0.8. It should be worth noticing that we are here interested only in fusion, but if we are interested in properties of fission-like fragments, we have to describe dynamics of the neck degree of freedom as well as those of deformations of the fragments along decaying paths. The shell correction energy is also neglected in formation process, since the united system is already well excited in the preceding processes of overcoming the Coulomb barrier. Therefore, we employ LDM for the potential as well as for the inertia-mass tensor for the collective motion. As for friction, we employ OBM for the calculation of the friction tensor. Now, we need only initial values to start calculations of trajectories with the Langevin equation. As stated in the introduction, our calculations of shape evolution start at the contact point. That is, the initial value of the distance is equal to the distance between the mass centers of the configuration, and that of the mass-asymmetry is equal to the mass-asymmetry of the entrance channel. As for the conjugate momenta, that for the mass-asymmetry is taken to be zero, while that for the radial motion can have various values, due to the Gaussian distribution given by the collision processes over the Coulomb barrier. In order to obtain the formation probability, we have to calculate probabilities $`F^J(p_0,T)`$ for the system to overcome the ridge line among numerous trajectories, starting with various initial momenta $`p_0`$, and then, to make an average of $`F^J(p_0,T)`$ with the weight of the Gaussian distribution of the initial momenta $`g^J(p_0,T_c)`$ obtained by SFM. Therefore, the formation probability is given as follows,

$$P_{\text{form}}^J(E_{\text{c.m.}})=𝑑p_0F^J(p_0,T)g^J(p_0,T_c),$$
(3)
where $`T_c`$ denotes the temperature of the system at the contact point. Actually, values of $`T_c`$ and $`T`$ are very close each other. In Fig. 2, the formation probability calculated for 48Ca+244Pu system is shown as a function of excitation energy $`E^{}`$. It is remarkable that it does not increase quickly as transmission coefficient usually used for fusion probability, but remains to be very small over the wide range of energy considered. This is in accord with the typical feature of the fusion hindrance observed.

## 3 FUSION CROSS SECTIONS AND RESIDUE CROSS SECTIONS OF SHE

Fusion probabilities are obtained with Eq. (2). Then, fusion excitation functions are calculated as usual,

$$\sigma _f(E_{\text{c.m.}})=\pi \lambda _{}^{}{}_{}{}^{2}\mathrm{\Sigma }_J(2J+1)P_{\text{fusion}}^J(E_{\text{c.m.}})$$
(4)
The results are shown in Fig. 3 for 48Ca+ 238U, +244Pu, +248Cm and 252Cf systems, together with the available data. Apparently, it is seen that the calculations reproduce the experimental feature of the energy-dependence, and furthermore, reproduce the absolute values of the first three systems very well systematically. It is worth emphasizing here that there is no adjustable parameter employed all through the calculations, except the probable choice of $`ϵ`$. Experiments on the last system are strongly desired for confirmation of validity of the present model.

From the comparisons made above, it is reasonable to consider that the fusion probabilities calculated by the present model are accurate, and therefore, they can be used for residue calculations with Eq. (1) in a high reliability. The remaining is to calculate accurately the other factor, i.e., the survival probability $`P_{\text{surv}}^J`$, which is made by KEWPIE. There, it is aimed to reduce ambiguities often introduced arbitrarily in analyses.

1. The level density parameter a is calculated with Töke and Swiatecki’s formula for the spherical and the saddle-point shapes. For particle evaporations, it is taken to accommodate the shell correction energy with Ignatyuk’s prescription with the standard damping energy of 18.5 MeV.
2. As for the saddle point and fission barrier height, we refer to Cohen-Plasil-Swiatecki’s LDM. But we also adopt the barrier height systematics by Dahlinger et al. for SHE, which is supposed to be more realistic.
3. Kramers factor for dynamical fission width is included with the reduced friction being $`510^{20}`$/sec which is consistent with OBM, together with Strutinski’s correction factor from inclusion of the collective degree of freedom around the spherical shape.

The frequencies of the potential at the spherical shape and the saddle point are temporally taken to be the same, the corresponding energy quantum taken equal to 1 MeV. Of course, one could use more realistic values, but they may give rise to a factor 2 or less. Thus, there is no free parameter again in the statistical calculations, except unknown masses or shell correction energies of isotopes of SHE. As stated in the introduction, there are many structure calculations, most of which predict shell correction energies of a few to several MeV. Møller et al.’s are larger ones among them, so we use their values with reduction factor 1/2 in keeping their tendency in mass-dependence. Calculated excitation functions of $`xn`$ reactions are shown in Fig. 4 for Z=114, 116, and 118, together with available experimental points. For Z=114 and 116, the calculations appear not to be inconsistent with the data. But the tendency of the predicted masses is not satisfactory. Nevertheless, it seems that in 48Ca+248Cm system a lower energy is more favourable. For the last system, the model provides predictions, on which experiments are being waited for.

In brief, the present two-step model has turned out to reproduce the experimental fusion excitation functions systematically for 48Ca induced reactions, as well as the measured $`xn`$ residue cross sections for SHE with the shell correction energies of Møller et al. with the reduction factor, the absolute values of which appear to be close to some of the mean-field calculations. Furthermore, we have made predictions on the fusion excitation function and the $`xn`$ cross sections for Z=118, which would stimulate experiments. A systematic study is now being made with typical theoretical predictions of nuclear masses and/or of shell correction energies for SHE. The present model will be applied to other elements and/or other entrance channels soon. Then, we will be able to answer the question raised in the introduction.

One of the authors (B.B) appreciates a post-doctoral fellowship provided by JSPS which enables the present collaboration. The authors acknowledge communications and encouragements by the experimental groups at Dubna, GSI, GANIL, RIKEN and JAERI. The work is partially supported by a Grant-in-Aids by JSPS (no. 13640278).


text:
# PSR J0609+2130: A disrupted binary pulsar?

## 1 Introduction

Radio pulsars are categorized as either “normal” or “recycled” objects. Normal pulsars are observationally more numerous and are predominantly young isolated objects with spin periods in the range 30 ms $`<P<8`$ s, strong inferred surface dipole magnetic field strengths ($`B10^{12}`$ G) and characteristic ages in the range $`10^3<\tau <10^7`$ yr. Recycled pulsars have shorter spin periods (1.5 ms $`<P<60`$ ms), weaker magnetic fields ($`10^8<B<10^9`$ G), larger ages ($`10^8<\tau <10^{10}`$ yr) and are usually members of binary systems with either white dwarf or neutron star companions.

The normal pulsars are thought to be born in supernovae and spin down due to magnetic dipole braking (see e.g. Gunn & Ostriker 1970). The favoured theory for the origin of recycled pulsars involves an old neutron star that is spun up (recycled) through the accretion of matter from a binary companion which overflows its Roche lobe (Bisnovatyi-Kogan & Komberg 1974). During the spin-up phase, the system is expected to be visible as an X-ray binary. The duration of the X-ray phase and amount of spin up depends upon the mass of the companion star. For those companions massive enough to explode as a supernova, the X-ray lifetime is relatively short ($`10^{67}`$ yr) and the most likely outcome is the disruption of the binary. Those systems fortunate enough to survive the explosion are the double neutron star (hereafter DNS) binaries (e.g., the original binary pulsar B1913+16). For less massive companions, where the period of spin up is longer ($`10^8`$ yr), the subsequent collapse leaves a white dwarf star in orbit around a rapidly spinning millisecond pulsar. For a detailed review of these evolutionary scenarios, see Bhattacharya & van den Heuvel (1991) and references therein.

In order to understand the population of recycled pulsars in detail, a statistically significant sample is required. Because of their short spin periods and, often, binary orbits, detecting these elusive objects is a computationally intensive process requiring state-of-the-art data acquisition systems, significant data storage and computational power for post processing. In 1994 the Penn State Pulsar Machine (PSPM), a new analogue filterbank spectrometer (Cadwell 1997), was installed at the 305-m Arecibo telescope and has been in regular use as a pulsar search and timing instrument ever since. In March 2002, we began using part of a high-speed computer cluster to process various sets of PSPM data taken during the latter stages of the Arecibo upgrade (1996–8). To date, we have discovered a dozen new pulsars and re-detected over 40 previously known ones. Preliminary details of the search were given by McLaughlin et al. (2003).

Here we report on an isolated 55.7-ms pulsar J0609+2130, one of the first new discoveries from our searches. While the spin parameters of the new pulsar are similar to other recycled pulsars, the lack of a binary companion poses interesting questions as to its origin. The layout of the rest of this paper is as follows: in §2 we briefly describe the survey observations, data acquisition system and analysis pipelines. The results of the analysis and subsequent timing observations of J0609+2130 are given in §3. Finally, in §4, we discuss the new pulsar’s properties in the context of its likely origin and evolution.

## 2 Survey observations and analysis

The data presented here were taken with the Arecibo radio telescope during a 9.5-hr maintenance period in October 1998 when upgrade operations required that the telescope was parked at $`35.7^{}`$ azimuth. The 430-MHz line-feed system, traditionally used for pulsar observations, was not available for use on this occasion. We therefore took advantage of the newly commissioned 430-MHz receiver in the Gregorian dome which was in a fixed position at $`3.9^{}`$ from the zenith. Although the 430-MHz Gregorian receiver111see http://www.naic.edu/~astro/RXstatus+ for full details of currently available Arecibo receiving systems illuminates less of the dish than the line-feed so that the forward gain is lower (11 K/Jy cf. 18 K/Jy for the line feed), the lower system temperature (45 K cf. 120 K) and reduced levels of interference make it a very attractive receiver for pulsar searching. In this setup, sources at declination $`\delta `$ drift through the 11-arcmin primary beam in $`44/\mathrm{cos}\delta `$ s. The total area covered was $`22`$ deg2 along an 11-arcmin slice through the Galactic anti-centre in the right ascension range $`1.5^\mathrm{h}<\alpha <10^\mathrm{h}`$ centred at $`\delta 21.5^{}`$.

The incoming signals from the telescope were passed to the PSPM which summed the two independent polarizations before 4-bit sampling the band into $`128\times 60`$-kHz channels every 80 $`\mu `$s. For this choice of sampling interval, a $`2^{19}`$-pt Fourier transform is well matched to the 44-s transit time of sources through the 430-MHz beam. The resulting data were then written directly to magnetic tape for offline analysis.

Searches for periodic and transient radio signals in the data were carried out using COBRA, a 180-processor Beowulf cluster at Jodrell Bank Observatory. Deferring a more complete discussion for a future paper, we now briefly outline the main procedures. To correct for the dispersive effects of the interstellar medium, the raw PSPM data were dedispersed at 392 different trial dispersion measures (DMs) in the range $`0`$ DM $`491.2`$ cm-3 pc using freely available analysis tools (Lorimer 2001). The resulting dedispersed time series were then passed to two different analysis pipelines. In the first, the data were Fourier transformed and the resulting amplitude spectra were searched for significant features indicating the presence of a periodic signal. To increase the sensitivity to narrow duty cycle pulses, spectra summing 2, 4, 8 and 16 harmonics were also searched. Candidates with signal-to-noise (S/N) ratios above 8 were then analysed in the time and radio frequency domain to produce diagnostic plots which were saved for later visual inspection. Further details can be found in Lorimer et al. (2000).

In the second analysis pipeline, the dedispersed data were searched in the time domain for pulses with S/N $`>5`$ using software developed by Cordes & McLaughlin (2003). This search was aimed at finding individual and giant pulses from some pulsars where the periodic signal may be below the threshold of the Fourier analysis (see e.g. Nice 1999). To maximize sensitivity over a variety of pulse widths, each time series was smoothed by up to 64 adjacent samples.

Given the measured system equivalent flux density of the 430-MHz Gregorian receiver (3.3 Jy), we estimate the sensitivity of the Fourier analysis to be $`0.5`$ mJy to long-period pulsars observed away from the Galactic plane. For millisecond pulsars, with typical periods of 5 ms and pulse duty cycles of order 60%, we estimate the sensitivity to be closer to 3 mJy. The detection threshold of the single-pulse search is about 0.5 Jy to the narrowest pulses. Full details of the search sensitivity will be presented elsewhere.

## 3 Results and timing observations

Two pulsar-like signals were detected in the periodicity search. One of these was the 33-ms Crab pulsar, B0531+21 (also seen in the single-pulse analysis) detected in several adjacent beams with S/N $`18`$. The remaining candidate, a 55.7-ms signal with DM $`40`$ and S/N $`9`$, was re-observed and confirmed as the new pulsar J0609+2130 in April 2002. Only one other previously known pulsar, B0525+21, lies within the region covered by this search. This was not detected in either the periodicity or the single-pulse search. A closer scrutiny of the raw search data revealed no detectable signal. Since the data were not obviously affected by radio-frequency interference, two possible effects are responsible for the non-detection: either (a) since this pulsar is known to null for approximately 25% of the time (Biggs 1990) it is quite likely that the survey observations caught this pulsar in a null state; or (b) the pulsar’s flux density was significantly reduced by diffractive interstellar scintillation. Given the normally high mean flux density of this pulsar (57 mJy at 408 MHz; Lorimer et al. 1995), nulling seems the most likely explanation for the non-detection.

To determine the detailed spin and astrometric parameters of J0609+2130, we have been carrying out regular 430-MHz timing observations since June 2002. These observations use the PSPM in timing mode, where the incoming signals are folded modulo the period predicted from an ephemeris initially derived from the confirmation observation. In timing mode, folded 1024-bin pulse profiles for each of the 128 frequency channels of the PSPM are written to disk every 180 s. We initially collected at least three of these 180-s observations per epoch in order to enable high-precision estimates of the pulse period. The timing analysis proceeded by forming the dedispersed profile across the band of the PSPM for each 180-s observation and estimating the pulse phase by cross correlating each profile with a high S/N template using procedures identical to those described in detail by Lorimer, Camilo & Xilouris (2002).

The topocentric pulse arrival times were fitted to a simple timing model involving spin and astrometric parameters using the TEMPO222http://pulsar.princeton.edu/tempo+ software package. Once a preliminary phase-connected timing solution was obtained, a better template was formed by phase-aligning the pulse profiles according to the best-fit model. The final template profile is shown in Fig. 1. As part of this refinement procedure, individual 180-s observations on a given day were phase aligned and averaged together to produce one pulse arrival time per observing session. The estimate of DM from the search analysis was improved by analysing PSPM observations with the 327-MHz and 610-MHz Gregorian receivers on February 17, 2003, and using TEMPO to fit for the dispersion delay between the frequency bands. The best-fit timing model for 30 arrival times spanning the MJD range 52427–52796 results in featureless timing model residuals with a post-fit RMS of 44 $`\mu `$s. The resulting parameters are given in Table 1.

## 4 Discussion

### 4.1 Basic properties of J0609+2130

PSR J0609+2130 is a very weak radio source. Using the calibration procedures described by Lorimer, Camilo & Xilouris (2002), we find the mean 430-MHz flux density of the pulsar $`S_{430}`$ to be only $`0.8\pm 0.1`$ mJy at 430 MHz. Given the distance estimate $`d1.2`$ kpc from the new electron density model of Cordes & Lazio (2003), the 430-MHz radio luminosity $`L=Sd^21.15`$ mJy kpc2. Only 10 out of the 612 pulsars currently in the online catalogue333http://www.atnf.csiro.au/research/pulsar/psrcat+ with measured luminosities are fainter sources than PSR J0609+2130.

The pulsar’s flux density is clearly affected by interstellar scintillation. We observe variations between 0.1 and 5 mJy as shown in Fig. 1. Indeed, for about 75% of the observations we have taken, the flux density is below our nominal survey threshold of about 1 mJy for this pulsar. Although our observations at other frequencies are not yet numerous enough to characterize the radio spectrum of J0609+2130 with any certainty, our preliminary flux estimate of 2 mJy at 327-MHz implies that it is a steep-spectrum source that is most readily detectable in low-frequency surveys.

Our timing measurements of PSR J0609+2130 reveal a period derivative $`\dot{P}=3\times 10^{19}`$. PSR J0609+2130 therefore appears to be a an old neutron star with characteristic age $`\tau =P/(2\dot{P})=2.8`$ Gyr and a relatively weak surface dipolar magnetic field $`B=3.2\times 10^{19}\sqrt{P\dot{P}}=4.2\times 10^9`$ G. Of all the sources in the current pulsar catalogue, the one which has properties most similar to J0609+2130 is PSR J2235+1506, a solitary 57.9-ms pulsar with $`\dot{P}=1.7\times 10^{19}`$ (Camilo, Nice & Taylor 1993, 1996). The positions of both these pulsars on the $`BP`$ diagram are shown in Fig. 2. As can be seen, J0609+2130 and J2235+1506 are the only two solitary pulsars in this part of the diagram. We shall return to the positions in the $`BP`$ plane later in the discussion.

### 4.2 The Galactic population and evolution of weakly magnetized solitary pulsars

Based on the close proximity and low radio luminosities of J0609+2130 and J2235+1506 we expect pulsars like them to be quite common in the Galaxy. Indeed, Kalogera & Lorimer (2000) estimated the Galactic population of pulsars like J2235+1506 to be $`5000/f`$, where $`f`$ is the fraction of $`4\pi `$ sr covered by the radio beam. For a beaming fraction of 30% (i.e. $`f=0.3`$) the Galactic population of these objects is likely to be of order 15,000.

What is the origin of weakly magnetized isolated pulsars like J0609+2130 and J2235+1506? One possibility is that they are simply part of the same population as the normal pulsars. To investigate this idea quantitatively, we can take the best-fitting initial magnetic field distribution in a recent population study of isolated pulsars by Arzoumanian, Chernoff & Cordes (2002). For their distribution, only about one in the expected Galactic population of $`10^9`$ neutron stars is born with $`B<10^{10}`$ G, clearly inconsistent with the population estimates made above. Unless there is a weak-field component to the neutron star magnetic field distribution that has been overlooked by population studies so far, a more efficient means of forming these weakly-magnetized neutron stars is required.

A pulsar born with a stronger magnetic field could, in principle, evolve to resemble PSR J0609+2130 through long-term decay of the magnetic field. While the existence of field decay for the pulsar population at large remains controversial (see e.g. Bhattacharya et al. 1992), available evidence suggests that, for objects comparable in age and magnetic-field strength to J0609+2130, any decay must act on extremely long timescales. Where such pulsars are found in binary systems, studies of their white-dwarf companions provide an independent age estimate through cooling models, affirming the long-lived nature of the pulsars and constraining the existence of field decay (Kulkarni 1986). Of course, J0609+2130 is an isolated pulsar, but it seems unlikely that non-binarity would somehow encourage rapid field decay.

Camilo, Nice & Taylor (1993) suggested that J2235+1506 might be the remains of a high-mass binary system that disrupted at the time of the second supernova explosion. In this scenario, PSRs J0609+2130 and J2235+1506 are the first-born neutron stars in their respective binaries and initially spin down as regular neutron stars. During Roche-lobe overflow from the subsequently evolved secondary star, mass accretion reduces the magnetic field (see e.g. Romani 1990), and also spins up the neutron star in the process. If the secondary is sufficiently massive to undergo a supernova explosion, the most likely outcome is disruption of the binary, releasing a newly-born neutron star and a recycled pulsar. Those binary systems more likely to survive this explosion as DNS binaries are the tightly-bound systems with short orbital periods. The wider binaries, which are much less likely to survive the explosion, release recycled pulsars with spin properties akin to J0609+2130 and J2235+1506 into the Galactic field. We hereafter refer to these as “disrupted recycled pulsars” (DRPs).

Using this evolutionary scenario for DRPs as a working hypothesis, we now return to the $`BP`$ diagram in Fig. 2 to investigate the post-accretion spin periods of DRPs after the accretion phase. In their analysis of DNS binaries, Arzoumanian, Cordes & Wasserman (1999; hereafter ACW) discuss two equilibrium spin period relations for the spin-up of a neutron star in a binary system. The first of these is the familiar limit for spherical accretion of matter in which the limiting spin period is set by the corresponding Keplerian orbital period at the so-called Alfvén radius (see e.g. Ghosh & Lamb 1992). On the $`BP`$ plane, this translates to a relationship of the form $`BP^{5/3}`$. The constant of proportionality in this expression depends on a number of factors such as the mass accretion rate, the radius and moment of inertia of the neutron star and the opacity of the accreting material. The uncertainties in these parameters result in a family of $`BP^{5/3}`$ curves. Two such curves, assuming accretion at the Eddington rate and at one-fifth of this value, are plotted on Fig. 2 using the relationship quoted by ACW. For any pulsar, in the absence of magnetic field decay, we can estimate the spin period at the end of the accretion phase by tracing a horizontal line from the current position on Fig. 2 to the spin-up line. Under these assumptions, we infer that the intial spin period of J0609+2130 and J2235+1506 would have been in the range 10–20 ms.

ACW pointed out that, for high-mass binary systems, the above picture is oversimplified since the lifetime of the donor star will be much less than the accretion time scale at which the Alfvén limit is reached. They show that this condition produces a second family of limiting period curves which, on the $`BP`$ plane, take the form $`BP^{7/2}`$. As noted by ACW, these curves are much more sensitive to uncertainties in the accretion parameters. In Fig. 2 we plot two such spin-up lines, again assuming accretion at the Eddington rate and one-fifth of this limit. It is clear that, for sub-Eddington accretion, the initial spin period of these DRPs would be much longer, perhaps very similar to their currently observed spin periods.

As noted by ACW, since pulsar characteristic ages $`\tau =P/(2\dot{P})`$ assume a negligible post-accretion spin period, one consequence of these accretion scenarios is that DNS binaries are younger objects than their characteristic ages suggest. If DRPs are formed by a similar process, the same applies to them. For example, a more likely “post-spinup” age for J0609+2130 based on a limiting spin period of 40 ms would be about half its characteristic age (i.e. about 1.5 Gyr).

### 4.3 Testing the DRP hypothesis

Is the DRP hypothesis consistent with the observations? One parameter that can be confronted by both theory and observation is the “survival probability” of the binary system, hereafter $`\eta `$, defined as the fraction of binary systems that remain bound after the second supernova explosion. Numerous authors have followed the orbital evolution of a wide variety of binary systems containing neutron stars using detailed Monte Carlo simulations (see e.g. Lipunov & Prokhorov 1984; Dewey & Cordes 1987; Fryer & Kalogera 1997). We use the relatively up-to-date simulations of Portegies Zwart & Yungelson (1998) who include a Maxwellian kick velocity distribution with a mean of 415 km s-1 in their simulations. From the results presented in Table 1 of their paper, for this model, we infer $`\eta 4\%`$. Based upon this low survival probability, we therefore expect DRPs to be much more common in the Galaxy than the DNS binaries.

To estimate $`\eta `$ from an observational perspective, we note that the expected spin parameters of the recycled pulsar we observe in a DNS binary are the same as those of the putative DRPs. We can therefore reasonably expect the pulsars in DNS binaries to have the same radio lifetimes as the DRPs. Hence, if $`N_{\mathrm{DNS}}`$ and $`N_{\mathrm{DRP}}`$ are the total numbers of DNS binaries and DRPs in the Galaxy, it follows that

$$\eta =\frac{N_{\mathrm{DNS}}}{N_{\mathrm{DNS}}+N_{\mathrm{DRP}}}.$$

Using the results of a recent population study of DNS binaries by Kim, Kalogera & Lorimer (2003), we take $`N_{\mathrm{DNS}}`$ population to be $`1500`$ (a factor of 2 higher than quoted by Kim, Kalogera & Lorimer 2003, to crudely account for the population of DNS binaries that do not coalesce in a Hubble time). Using the above estimate of $`N_{\mathrm{DRP}}=15,000`$, we find $`\eta 1500/(15000+1500)=9\%`$.

One outstanding issue is the observed numbers of DRPs. Given the above estimate for $`\eta `$, we would expect to see roughly 10 DRPs for each DNS binary. So far, however, we have identified only two DRPs: J0609+2130 and J2235+1506, compared to four DNS binaries in Fig. 2. Three main possibilities which might explain this apparently significant deficit of DRPs are: (a) DRPs are present in the observed sample, but are hard to identify since they have no binary companion. Although there are perhaps a further four DRP-like pulsars known (see Table 1 of Kalogera & Lorimer 2000), a more thorough analysis of the sample would be worthwhile. (b) The neutron star kick velocity distribution could be lower than assumed in the population synthesis. This would lead to an increase in $`\eta `$ and reduce the discrepancy. A population synthesis dedicated to this question would certainly be useful. (c) The above DRP population estimates are somewhat crude and are based currently on J2235+1506 which is a nearby, low-luminosity object. A revision of these estimates to better account for small-number statistics is strongly recommended.

In summary, while it seems currently plausible that J0609+2130 and J2235+1506 are examples of disrupted X-ray binary systems, there are a number of open issues linking the two populations. Careful population studies to properly identify and estimate the DRP population and binary population syntheses to predict the survival rate are both strongly recommended. Further observational input should come both in the substantial numbers of pulsars now being found in the Parkes Multibeam survey (Manchester et al. 2001), and also proper motion measurements. For J2235+1506, the implied transverse speed is $`V_t=100\pm 40`$ km s-1 (Camilo, Nice & Taylor 1996). If J0609+2130 has a similar proper motion, it should be detectable in our on-going Arecibo timing observations by early 2005.

## Acknowledgments

The Arecibo observatory, a facility of the National Astronomy and Ionosphere Center, is operated by Cornell University in a co-operative agreement with the National Science Foundation (NSF). We thank Alex Wolszczan for making the PSPM freely available for use at Arecibo. Without this superb instrument, the results presented here would not have been possible. We also thank Paulo Freire for observing assistance and the referee, Simon Johnston, for many useful comments. DRL is a University Research Fellow funded by the Royal Society. MAM is an NSF MPS-DRF Fellow.