text:
# Quantum effects in small-capacitance single Josephson junctions

## I Introduction

Small-capacitance superconducting tunnel junctions provide an ideal system for studying the interplay between quantum mechanically conjugate variables: the Josephson-phase difference across the junction and the charge on the junction electrode. Their behavior is influenced by dissipation, such as discrete tunneling of quasiparticles, coupling to an electromagnetic environment, etc. The simplest and the most fundamental example is the single junction. Current-voltage ($`I`$-$`V`$) characteristics of single junctions have been the subject of extensive theoretical investigations.Averin and Likharev (1991); Schön and Zaikin (1990) A small capacitance $`C`$ gives a large charging energy $`e^2/2C`$, and thus the theory predicts that at low enough temperatures, a superconducting small-capacitance tunnel junction behaves like an insulator so long as the junction is isolated from the electromagnetic environment. Experimentally, however, the observation of this charging effect (Coulomb blockade) in a single junction is not straightforward, because the measurement leads attached to the single junction should have a high enough impedance to isolate the junction.Averin and Likharev (1991) For this reason, thin-film resistors were employed for the leads, and Coulomb blockade was successfully observed in single Josephson junctions.Haviland et al. (1991) Tunnel-junction arrays were also tried for the leads, and an increase of differential resistance around $`V=0`$ was reported.Geerligs (1990); Shimazu et al. (1997)

We use one-dimensional (1D) arrays of dc superconducting quantum interference devices (SQUIDs) for the leads.Watanabe and Haviland (2001) The advantage of this SQUID configuration is that the effective impedance of an array can be varied in situ by applying an external magnetic field perpendicular to the SQUID loop. Thus, the zero-bias resistance of the SQUID arrays at low temperatures can be controlled over several orders of magnitude. This phenomenon has been extensively studied in terms of the superconductor-insulator (SI) transition.Chow et al. (1998); Haviland et al. (2000, 2001); Watanabe and Haviland (2002) The single junction in our samples, on the other hand, does not have a SQUID configuration, and therefore its parameters are practically independent of the external magnetic field. This enables us to study the same single junction in different environments. We show that the $`I`$-$`V`$ curve of the single junction is indeed sensitive to the state of the environment. Furthermore, we can induce a transition to a Coulomb blockade of the single junction by increasing the zero-bias resistance $`R_0^{}`$ of the SQUID arrays. When $`R_0^{}`$ is much higher than the quantum resistances: $`R_Kh/e^226`$ k$`\mathrm{\Omega }`$ for quasiparticles and $`R_Qh/(2e)^26.5`$ k$`\mathrm{\Omega }`$ for Cooper pairs, single-junction $`I`$-$`V`$ curve has a region of negative differential resistance. The negative differential resistance is a result of coherent single-Cooper-pair tunneling from the viewpoint of the theory for current-biased single Josephson junctions whose electromagnetic environment has a sufficiently high impedance.Averin and Likharev (1991); Schön and Zaikin (1990)

In the context of the theory, we calculate the $`I`$-$`V`$ curves numerically by following Ref. Geigenmüller and Schön, 1988, and compare with the experimental ones at $`R_0^{}R_K`$. The numerical results are in qualitative agreement with the experiments,Watanabe et al. (2001a) however, some quantitative discrepancies suggest that even at $`R_0^{}R_K`$, the single-junction $`I`$-$`V`$ curve would be influenced by the environment. According to the perturbation theoryFalci et al. (1991); Ingold and Nazarov (1992); Ingold and Grabert (1999) that deals with the effect of a linear-impedance environment $`Z(\omega )`$ on single Josephson junctions, $`\text{Re}[Z(\omega )]`$, or the dissipation, plays a key role. An insulating behavior characterized by the Coulomb blockade of Cooper-pair tunneling is expected for the single junction when $`\text{Re}[Z(\omega )]R_Q`$, i.e., the dissipation is sufficiently small, because the large $`\text{Re}[Z(\omega )]`$ suppresses the charge fluctuation. The dissipative dynamics of single Josephson junctions has also been discussed by studying resistively shunted junctions both theoreticallySchön and Zaikin (1990); Kato and Imada (2000); Chung (2002) and experimentally.Yagi et al. (1997); Penttilä et al. (2001) As the dissipation is increased (the shunt resistance is decreased), the tunneling of the Josephson phase is suppressed, and as a result the Josephson junction exhibits a transition from an insulator to a superconductor. The theoretical phase diagramSchön and Zaikin (1990) for this SI transition was supported by the experimentsYagi et al. (1997); Penttilä et al. (2001) qualitatively. The superconducting behavior has been observed experimentally in single Josephson junctions not only by shunting with a resistor but also by biasing with a low-impedance circuit,Steinbach et al. (2001) i.e., by achieving a low-impedance electromagnetic environment. Thus, we expect that there should be a SI transition driven by the electromagnetic environment. We discuss the environment-driven SI transition in our system, where the environment is tunable with the SQUID arrays.

## II Theory

The Hamiltonian of a single Josephson junction in an environment with sufficiently high impedance is written as

$$H=\frac{Q^2}{2C}E_J\mathrm{cos}\varphi ,$$
(1)
where $`Q`$ is the charge on the junction electrode, $`C`$ is the capacitance of the junction, $`E_J`$ is the Josephson energy, and $`\varphi `$ is the Josephson-phase difference across the junction. The charge $`Q`$ and $`\mathrm{}\varphi /2e`$ are quantum-mechanically conjugate valuables, and a set of the eigenfunctions are Bloch waves of the form

$$\psi (\varphi )=u(\varphi )\mathrm{exp}(i\varphi q/2e),$$
(2)
where $`q`$ is called quasicharge and $`u(\varphi )`$ is a periodic function,

$$u(\varphi +2\pi )=u(\varphi ).$$
(3)
The energy eigenvalue $`E`$ plotted as a function of $`q`$ has a band structure, and in all the allowed bands, it is $`2e`$ periodic. An example of the energy diagram for $`E_J/E_C=0.3`$, where $`E_Ce^2/2C`$ is the charging energy, is shown in Fig. 1a.

Under constant current bias $`I_\mathrm{x}`$, in the absence of quasiparticle or Cooper-pair tunneling, $`q`$ increases uniformly in time according to

$$\frac{dq}{dt}=I_\mathrm{x},$$
(4)
so that the state of the system advances toward higher $`q`$ within a given band as time goes on. The average voltage is given by

$$V=\underset{i_b,q}{}P(i_b,q)\frac{dE(i_b,q)}{dq},$$
(5)
where $`i_b`$ is the band index and $`P(i_b,q)`$ is the probability that the system is in the state $`(i_b,q)`$. The probability $`P(i_b,q)`$ can be calculated by solving a set of coupled differential equations of the form

$$\frac{dP(i_b,q)}{dt}=\underset{i_b^{},q^{}}{}A(i_b,q,i_b^{},q^{})P(i_b^{},q^{})=0,$$
(6)
where the matrix element $`A(i_b,q,i_b^{},q^{})`$ describes the rate of transition between the states $`(i_b,q)`$ and $`(i_b^{},q^{})`$. The dominant process for $`A(i_b,q,i_b^{},q^{})`$ depends on the magnitude of $`I_\mathrm{x}`$.

An example of the theoretical $`I`$-$`V`$ curve is shown in Fig. 1b. For sufficiently small $`I_\mathrm{x}`$ (region CB), the dominant process is stochastic quasiparticle tunneling, where $`q`$ changes by $`e`$. This tunneling always occurs along the energy parabola $`q^2/2C`$, such that $`i_b`$ changes by 0 or $`+1`$ if the initial state is in the lowest band ($`i_b=1`$) and by $`\pm 1`$ for all the other initial states.Schön and Zaikin (1990) The rate for the quasiparticle tunneling is given by

$$\mathrm{\Gamma }(\mathrm{\Delta }E)=\frac{\mathrm{\Delta }E/e^2R_{\mathrm{qp}}}{\mathrm{exp}(\mathrm{\Delta }E/k_BT)1},$$
(7)
where $`R_{\mathrm{qp}}`$ is the quasiparticle resistance and $`\mathrm{\Delta }E`$ is the difference in energy between the initial $`(i_b,q)`$ and final $`(i_b^{},q^{})`$ states,

$$\mathrm{\Delta }EE(i_b^{},q^{})E(i_b,q).$$
(8)
At sufficiently low temperatures, the tunneling with $`\mathrm{\Delta }E>0`$ is extremely unfavorable, and the $`I`$-$`V`$ curve is highly resistive (Coulomb blockade, region CB in Fig. 1b). For larger $`I_\mathrm{x}`$ (region BO in Fig. 1b), the quasicharge is frequently driven to the boundary of the Brillouin zone, $`q=e`$, then taken to $`e`$ as a Cooper pair tunnels (Bloch oscillation). This process decreases $`V`$, and as a result the $`I`$-$`V`$ curve has a region of negative differential resistance, or “back bending,” in the low-current part. For still larger $`I_\mathrm{x}`$ (region ZT in Fig. 1b), Zener tunneling becomes important, and $`V`$ increases again. In Zener tunneling, no quasicharge is transferred but the state of the system jumps from one band to another as it passes by the gap between the bands. The probability of Zener tunneling from band $`i_b`$ to $`i_b+1`$ or vice versa is given by

$$P_Z=\mathrm{exp}\left[\frac{\pi }{8}\frac{(\mathrm{\Delta }E)^2}{i_bE_C}\frac{e}{\mathrm{}I_\mathrm{x}}\right].$$
(9)
Following Ref. Geigenmüller and Schön, 1988 which takes into account the above tunneling processes (quasiparticle, Cooper pair, and Zener), we have calculated the $`I`$-$`V`$ curve numerically.111The numerical code was written by R. L. Kautz. The parameters for the calculation are $`E_J/E_C`$, $`k_BT/E_C`$, and

$$\alpha \frac{h}{\pi ^2e^2R_{\mathrm{qp}}}.$$
(10)
The current and the voltage are in units of $`e/R_{\mathrm{qp}}C`$ and $`e/C`$, respectively.

As we have seen in Fig. 1b, a typical $`I`$-$`V`$ curve consists of three regions, so that it is characterized by the local voltage maximum, or blockade voltage $`V_b`$, and the local current minimum, or crossover current $`I_{\mathrm{cr}}`$. (Here, we have to mention that the back-bending feature is blurred out if $`k_BT/E_C`$ or $`\alpha `$ is increased considerably.) An analytical expression of $`V_b`$ and $`I_{\mathrm{cr}}`$ has been obtained theoretically for limiting cases.Schön and Zaikin (1990) The value of $`V_b`$ is a function of $`E_J/E_C`$, and given by

$$V_b\{\begin{array}{cc}0.25e/C\hfill & \text{for }E_J/E_C1\hfill \\ \delta _0/e\hfill & \text{for }E_J/E_C1\text{,}\hfill \end{array}$$
(11)
as $`T0`$, where

$$\delta _0=\frac{e^2}{C}\mathrm{\hspace{0.17em}8}\left(\frac{1}{2\pi ^2}\right)^{1/4}\left(\frac{E_J}{E_C}\right)^{3/4}\mathrm{exp}\left[\left(8\frac{E_J}{E_C}\right)^{1/2}\right]$$
(12)
is the half-width of the lowest energy-band. As for $`I_{\mathrm{cr}}`$,

$$I_{\mathrm{cr}}\left(I_Z\frac{e}{R_{\mathrm{qp}}C}\right)^{1/2}$$
(13)
is expected for $`\alpha (E_J/E_C)^21`$ and $`T0`$, where

$$I_Z\frac{\pi }{8}\frac{eE_J^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}}}{\mathrm{}E_C}$$
(14)
is the Zener breakdown current. Note that $`I_{\mathrm{cr}}`$ is much smaller than $`I_Z`$. When we compare our experimental results with the theory, we need theoretical prediction for finite $`k_BT/E_C`$, and arbitrary $`E_J/E_C`$ and $`\alpha `$. For this reason we have done the numerical calculation. The measured $`V_b`$ and $`I_{\mathrm{cr}}`$ are compared with the calculation in Sec. IV.3.

## III Experiment

### III.1 Sample preparation and characterization

A schematic diagram of the Al/Al2O3/Al samples is shown in Fig. 2a.

A single Josephson junction is at the center of Fig. 2a. The junction area is 0.1$`\times `$0.1 $`\mu `$m2. For a scanning electron micrograph of the single junction, see Fig. 1 of Ref. Watanabe and Haviland, 2001. On each side of the single junction, there are two leads enabling four-point measurements of the single junction. A part of each lead close to the single junction consists of a 1D array of dc SQUIDs. We show in Fig. 2b a scanning electron micrograph of a part of the SQUID array. Each of the electrodes in the array is connected to its neighbors by two junctions in parallel, thus forming a dc SQUID between nearest neighbors. The area of each junction in the SQUID array is 0.3$`\times `$0.1 $`\mu `$m2 and the effective area of the SQUID loop is 0.7$`\times `$0.2 $`\mu `$m2. All of the samples have the same configuration, except for the number $`N`$ of junction pairs in each SQUID array. The samples are characterized by the normal-state resistance $`R_n`$ of the single junction, the normal-state resistance $`r_n^{}`$ per junction pair of the SQUID arrays, and $`N`$. The parameters of the samples are listed in Table 1.

It is difficult to determine the uniformity of these 1D SQUID arrays when only the series resistance of all junctions can be measured. Previous workHaviland et al. (2000) showed that the variations of the normal-state resistance per junction pair were less than 6% for six arrays having nominally identical junctions made simultaneously on the same chip, with the number of junction pairs ranging from 7 to 255. Furthermore, critical-current switching at zero magnetic field for these arrays indicated a high degree of uniformity (Fig. 4 of Ref. Haviland et al., 2000), and the magnetic-field dependence of the arrays was smooth. For the samples of this work, the uniformity of the SQUID arrays could be estimated by comparing the two arrays on the right- and left- hand sides of the single junction. On an average, the normal-state resistance of the arrays agreed to within 5%. However, if the tunnel-barrier thickness was identical for all junctions, one would expect $`R_n/r_n^{}6`$ from the junction area for our samples. In Table 1, we see a large variation in $`R_n/r_n^{}(=117)`$, which would be mainly due to the variation of $`R_n`$ for the following reasons: In order to obtain $`r_n^{}`$, we do not examine one junction pair in the arrays, but measure two arrays in series and divide by $`2N`$. This kind of “averaging” does not occur when we determine $`R_n`$. Moreover, the single junction is much smaller than the junctions in the SQUID arrays, and small variations in the lithography and evaporation angle cause larger changes in the area of the smaller junction, thus leading to larger variations in the junction area and thereby the resistance and capacitance.

The tunnel junctions were fabricated on a SiO2/Si substrate with electron-beam lithography and a double-angle shadow evaporation technique.Haviland et al. (1996) We employed a double-layer resist (ZEP520/PMGI) method in order to achieve a large enough undercut profile for the desired shift, which in our case was 0.2 $`\mu `$m. The electron-gun evaporation of Al was done at a rate of $`0.10.4`$ nm/s in a turbo-molecular-pumped vacuum system with the base pressure of $`<2\times 10^6`$ Pa. The vacuum system is equipped with an angle stage, and with the stage tilted to an angle $`+\theta 13^{}`$, Al was deposited to a thickness of 20 nm. After deposition of the base electrode, O2 was introduced into the chamber to a pressure of $`1.34.0`$ Pa for a period of $`210`$ min. The oxidation conditions, or the thickness of the Al2O3 layer, determines $`R_n`$ and $`r_n^{}`$. After oxidation, the top layer of Al was deposited to a thickness of 30 nm with the stage tilted to an angle $`\theta `$.

In the above process, bonding pads ($`0.3\times 0.35`$ mm2) were also fabricated simultaneously, i.e., the pads are 20+30=50 nm-thick Al. Onto the pads, 25 $`\mu `$m-diameter Al wires were wedge bonded.

### III.2 Measurements

The samples were measured in a 3He-4He dilution refrigerator (Oxford Instruments, Kelvinox AST-Minisorb) at $`T=0.020.6`$ K. The temperature was determined by measuring the resistance of a ruthenium-oxide thermometer fixed at the mixing chamber. To measure the resistance of the thermometer, we used an ac resistance bridge (RV-Elekroniikka, AVS-47). Magnetic fields on the orders of $`110`$ mT were applied by means of a superconducting solenoid. In this magnetic-field range, the temperature error of ruthenium-oxide thermometers due to the magnetoresistance is negligibly small (e.g., a typical value of the error at $`T=0.05`$ K is less than 0.1%).Watanabe et al. (2001b) The samples were placed inside a copper rf-tight box that was thermally connected to the mixing chamber. All the leads entering the rf-tight box were low-pass filtered by 1 m of Thermocoax cable.Zorin (1995)

We measured the $`I`$-$`V`$ curve of the single junction in a four-point configuration (see Fig. 2a). The bias was applied through one-pair of SQUID-array leads, and the potential difference was measured through the other pair of SQUID-array leads with a high-input-impedance instrumentation amplifier based on two operational amplifiers (Burr-Brown OPA111BM, $`10^{13}10^{14}`$ $`\mathrm{\Omega }`$ input impedance according to the data sheet) and a preamplifier \[Stanford Research Systems (SRS) SR560\]. The current was measured with a current preamplifier (SRS SR570). When the voltage drop at the SQUID arrays was much larger than that at the single junction, the single junction was practically current biased.

The SQUID arrays could be measured in a two-point configuration (same current and voltage leads) on the same side of the single junction. Note that the two arrays are connected in series and that current does not flow through the single junction. In order to obtain $`R_n`$ and $`r_n^{}`$, we measured at $`T=24`$ K (above the superconducting transition temperature of Al).

## IV Results and discussion

### IV.1 SQUID arrays: tunable electromagnetic environment for the single junction

The effective Josephson energy $`E_J^{}`$ between adjacent islands in the SQUID arrays is modulated periodically by applying an external magnetic field $`B`$ perpendicular to the SQUID loop,

$$E_J^{}\left|\mathrm{cos}\left(\pi \frac{BA}{\mathrm{\Phi }_0}\right)\right|,$$
(15)
so long as $`B`$ is sufficiently smaller than the critical field, where $`A`$ is the effective area of the SQUID loop and $`\mathrm{\Phi }_0=h/2e=2\times 10^{15}`$ Wb is the superconducting flux quantum. Figure 3 shows how the $`I`$-$`V`$ curve of two SQUID arrays on the same side of the single junction changes depending on $`B`$.

The $`I`$-$`V`$ curve is Josephson-like at $`B=0`$, while a Coulomb blockade is developed at $`B=7`$ mT, where the normalized flux

$$f\frac{BA}{\mathrm{\Phi }_0}$$
(16)
is close to 1/2 for our samples with $`A=0.14`$ $`\mu `$m2. Note that the current scale for $`B=0`$ is $`10^4`$ times larger than that for $`B=7`$ mT. In this manner, the $`I`$-$`V`$ curve is varied periodically until $`B`$ becomes comparable to the critical field. At $`B=50`$ mT, the superconductivity is suppressed considerably, and the $`I`$-$`V`$ curve is almost linear. The magnetic-field dependence of the SQUID arrays is summarized in Fig. 4, where the zero-bias resistance $`R_0^{}`$ of the same SQUID arrays is plotted vs. $`B`$.

As expected from Eq. (15), $`R_0^{}`$ oscillates in the low-field regime. In the shaded regions, a Coulomb blockade is clearly observed and the maximum value of $`R_0^{}`$ can be as high as $`10^910^{10}`$ $`\mathrm{\Omega }`$ (see Fig. 2 of Ref. Watanabe and Haviland, 2001 for $`R_0^{}`$ of Sample B at $`0B9`$ mT). At $`B>40`$ mT, the oscillation is suppressed and $`R_0^{}`$ approaches the value in the normal state, $`2Nr_n^{}`$.

In summary of this subsection, the SQUID arrays, which act as an electromagnetic environment for the single junction in this work, can be tuned over a wide range in situ with $`f`$. Henceforth, we use $`f`$ defined by Eq. (16) rather than $`B`$ for the purpose of indicating the external field.

### IV.2 Coulomb blockade induced in the single Josephson junction by tuning the environment

In this subsection, we take Sample G for example. We look at the $`I`$-$`V`$ curve and the influence of the electromagnetic environment on it. (The $`I`$-$`V`$ curves of Sample I are shown in Figs. 3–5 of Ref. Watanabe and Haviland, 2001.) Figure 5 shows the $`I`$-$`V`$ curves for the single junction at $`f=0`$. The vertical broken lines represent the superconducting-gap voltage $`\pm 2\mathrm{\Delta }_0/e`$, which is $`\pm 0.34`$ mV for Al according to Ref. Kittel, 1996. In a conventional Josephson junction ($`E_JE_C`$), a finite current flows at $`V=0`$ as expected within the RCSJ (resistively and capacitively shunted junction) model.Tinkham (1996) This supercurrent has also been observed in small-capacitance Josephson junctions ($`E_JE_C`$), and when a small-capacitance junction was biased with a carefully filtered low-impedance circuit,Steinbach et al. (2001) the magnitude of the supercurrent reached the value predicted by the theory. Our single Josephson junction is biased with SQUID arrays, whose $`R_0^{}`$ is much higher than the impedance ($`24`$ $`\mathrm{\Omega }`$) of the biasing circuit of Ref. Steinbach et al., 2001 even at $`f=0`$. For the case of Fig. 5 (Sample G at $`f=0`$), $`R_0^{}=0.21`$ k$`\mathrm{\Omega }`$. Ingold and GrabertIngold and Grabert (1999) considered how phase fluctuations suppress the supercurrent in a small-capacitance Josephson junction. The small capacitance of a junction is a source of phase fluctuations because the charge on the capacitance is the conjugate variable to the Josephson-phase difference. In addition, there is a strong relationship between the fluctuations and the electromagnetic environment. They assumed a purely Ohmic environment $`R`$, and showed that the supercurrent is rounded as $`R`$ is increased. In such a context, the shape of the $`I`$-$`V`$ curve in Fig. 5 may be explained for the region in the vicinity of $`V=0`$, although our SQUID arrays are not represented by a purely Ohmic impedance. At higher voltages, finite currents still flow in Fig. 5 even within the superconducting gap. Theoretically, at a finite voltage $`|V|<2\mathrm{\Delta }_0/e`$, current can flow when the tunneling Cooper pairs can release their excess energy $`2eV`$, e.g., to the degrees of freedom in the electromagnetic environment. Quasiparticle excitations would also be responsible for the current especially at $`|V|2\mathrm{\Delta }_0/e`$. Below we will concentrate on the immediate vicinity of $`V=0`$ (typically $`\pm 30`$ $`\mu `$V), where the $`I`$-$`V`$ curve is most sensitive to the state of the electromagnetic environment.

Figure 6 shows how the single-junction $`I`$-$`V`$ curve depends on the environment.

As $`f`$ is varied, the curve develops a Coulomb blockade. As we have mentioned in Sec. I, the Josephson energy of the single junction is independent of $`f`$ because it does not have a SQUID configuration and the field $`f\mathrm{\Phi }_0/A<7`$ mT applied here is much smaller than the critical field. The electromagnetic environment for the single junction (the SQUID arrays), however, is strongly varied with $`f`$. The behavior of the single junction demonstrated in Fig. 6 does not result from the magnetic-field influence on the single-junction $`I`$-$`V`$ curve, but rather from an environmental effect on the single junction. This experiment demonstrates in a direct way that the single-junction $`I`$-$`V`$ curve is indeed sensitive to the electromagnetic environment.

Figure 7 shows the $`I`$-$`V`$ curves of the two SQUID-array leads connected in series at the same $`f`$ as in Fig. 6.

The $`I`$-$`V`$ curves of the leads are nonlinear, and in general the SQUID array is not described by a liner impedance model.Haviland et al. (2000) However, the environment may be characterized by $`R_0^{}`$. From (a) to (d) in Fig. 7, $`R_0^{}=0.21`$ k$`\mathrm{\Omega }`$, 11 k$`\mathrm{\Omega }`$, 0.33 M$`\mathrm{\Omega }`$, and 6.5 M$`\mathrm{\Omega }`$, respectively. Coulomb blockade is distinct only when $`R_0^{}R_K`$, which is consistent with the theoretical conditions for the clear observation of Coulomb blockade in single junctions.Ingold and Nazarov (1992) For an arbitrary linear environment characterized by $`Z(\omega )`$, $`\text{Re}[Z(\omega )]R_K`$ is required for the Coulomb blockade of single-electron tunneling and $`\text{Re}[Z(\omega )]R_Q`$ for that of Cooper-pair tunneling.Ingold and Nazarov (1992) It is interesting to note that while the $`I`$-$`V`$ curve of the leads is “Josephson-like” (differential resistance is lower around $`V=0`$) at all $`f`$ shown in Fig. 7, that of the single junction is “Coulomb-blockade-like” at $`f=0.40`$ and exhibits a clear Coulomb blockade at $`f0.45`$.

### IV.3 Comparison with the numerical calculation

The region of negative differential resistance seen at the bottom curve of Fig. 6 is related to coherent tunneling of single Cooper pairs according to the theoryAverin and Likharev (1991); Schön and Zaikin (1990) of a current-biased single Josephson junction in an environment with sufficiently high impedance. Following Ref. Geigenmüller and Schön, 1988, we have calculated the blockade voltage $`V_b`$ numerically as a function of $`E_J/E_C`$.Watanabe et al. (2001a) The numerical results are compared with the measured $`V_b`$ in Fig. 8.

The data from Ref. Haviland et al., 1991 for the samples with the same nominal junction area ($`0.1\times 0.1`$ $`\mu `$m2) as in this work are also plotted. For Sample G we used the data at $`f=0.46`$ (the bottom curve in Fig. 6) in order to obtain $`V_b`$. At $`f=0.46`$ the voltage drop at the SQUID arrays is $`10^2`$ times larger than that at the single junction, and the single junction is therefore considered to be current biased. Compare the voltage scale of Figs. 6 and 7d. We calculated $`E_J`$ from $`R_n`$,

$$E_J=\frac{h\mathrm{\Delta }_0}{\mathrm{\hspace{0.17em}8}e^2R_n}.$$
(17)
For $`E_C`$ we employed $`c_s=130`$ fF/$`\mu `$m2, and with this value the experimental data, especially those for this work (the circles in Fig. 8), agree with the numerical calculation. Actually, a smaller value, $`c_s=45\pm 5`$ fF/$`\mu `$m2 (Ref. Lichtenberger et al., 1989), which was obtained for the junctions with $`3\times 28`$ $`\mu `$m2 and $`7\times 54`$ $`\mu `$m2, has been frequently employed.Chow et al. (1998); Haviland et al. (2000, 1991); Shimazu et al. (1997); Haviland et al. (2001) Our apparently large $`c_s`$ may be partly explained by distributed capacitance of the SQUID arrays or by the non-infinite impedance of the environment (SQUID arrays), which would have the effect of shunting the single junction with an additional capacitance (see Fig. 11). We also note that the uncertainty in $`c_s`$ seems to be large when the junction area is on the order of 0.01 $`\mu `$m2 or smaller, as we have discussed in Ref. Watanabe and Haviland, 2001.

Fixing $`C=1.3`$ fF, we tried to reproduce the measured $`I`$-$`V`$ curve by optimizing the other parameters of the numerical calculation, $`\alpha `$ and $`k_BT/E_C`$. Examples are shown in Fig. 9, where we adjusted $`\alpha `$ roughly by considering the current scale.

For both Samples G and I, the slope of the experimental data at $`|V|<V_b`$ and $`I0`$ can be fit by adjusting $`\alpha `$, so that it is consistent with the theoretical curve for $`k_BT/E_C=0.05`$ ($`T=0.036`$ K for $`C=1.3`$ fF), which is close to the temperature of the mixing chamber, $`0.020.04`$ K. At higher bias currents, however, the theory predicts for $`k_BT/E_C=0.05`$ a “back bending” much bigger than that in the experimental data. In terms of the magnitude of “back bending,” the curves for higher temperatures, $`k_BT/E_C=0.30.6`$ ($`T=0.210.43`$ K), are more similar to the data. “Back bending” in the theoretical curve disappears when the temperature is raised to $`k_BT/E_C=0.8`$ ($`T=0.57`$ K). Although the theory explains the general shape of the measured $`I`$-$`V`$ curve, it does not yield a complete agreement. It would be reasonable to assume that the effective temperature of the single junction depends on the bias current, and increases with increasing current due to heating. This assumption reduces the discrepancy between the experimental data and the theory. Another possible origin of the discrepancy would be the imperfect isolation of the single junction. While the theory assumes that the single junction is completely isolated from the environment, we have biased the single junction with the SQUID arrays having non-infinite $`R_0^{}`$. Moreover, the SQUID arrays are a nonlinear environment.

Finally in this subsection, we estimate $`R_{\mathrm{qp}}`$ from $`I_{\mathrm{cr}}`$ based on the theory. We plot the measured $`I_{\mathrm{cr}}`$ as a function of $`E_J/E_C`$ in Fig. 10 together with some theoretical curves based on our numerical calculation.

For the estimate of $`R_{\mathrm{qp}}`$ we used the numerical results for $`k_BT/E_C=0.3`$ and 0.5 ($`T=0.21`$ and 0.36 K), because with these values the numerical calculation gave an $`I`$-$`V`$ curve similar to the experimental one around $`I_{\mathrm{cr}}`$ as we have seen in Fig. 9. The estimated $`R_{\mathrm{qp}}`$ is listed in Table 1. The order of $`R_{\mathrm{qp}}`$ ($`R_{\mathrm{qp}}=10^010^1`$ M$`\mathrm{\Omega }`$) and the ratio to $`R_n`$ ($`R_{\mathrm{qp}}/R_n=10^210^3`$) are within a reasonable range.

### IV.4 Superconductor-insulator transition in single Josephson junctions driven by the environment

The numerical calculation in the preceding subsection agrees with the experiment qualitatively, however, the apparently large $`c_s`$ and the discrepancies in Fig. 9 suggest that the single junction would still be dependent on the environment even at $`R_0^{}R_K`$. The influence of a linear electromagnetic environment $`Z(\omega )`$ on small-capacitance Josephson junction has been studied in the framework of the perturbation theoryFalci et al. (1991); Ingold and Nazarov (1992); Ingold and Grabert (1999) by modeling the whole circuit as in Fig. 11a.

By assuming an Ohmic environment $`Z(\omega )=R`$, the Cooper-pair current-voltage characteristics have been calculated for $`E_J/E_C1`$. The theory predicts Coulomb blockade of Cooper-pair tunneling for $`R_Q/R1`$, and supercurrent at or in the vicinity $`V=0`$ for $`R_Q/R1`$. By comparing Fig. 2a with Fig. 11a, we realize that the total impedance $`Z_t`$ seen by the junction in Fig. 2a corresponds to $`Z`$ in Fig. 11a. If we assume in Fig. 2a that the SQUID arrays alone determine the effective environment for the single Josephson junction at relevant frequencies, where all the four leads are shunted at the far ends of the SQUID arrays, and that the SQUID array is described by a linear impedance, we obtain $`Z_t=\left(Z_1^1+Z_2^1\right)^1+\left(Z_3^1+Z_4^1\right)^1Z_1`$ (all the SQUID arrays are nominally the same). Note that the dc $`I`$-$`V`$ curve is the result of high-frequency quasicharge dynamics of the single junction, and thus one must consider what the junction sees at high frequencies. The large capacitance introduced by the cryostat leads will effectively shunt the sample with a low impedance at the relevant frequencies. This is why we fabricated the SQUID arrays in the immediate vicinity of the single junction in order to avoid the lead capacitance shunting the single junction, and this is why we assume that the SQUID arrays alone determine the environment for the single junction. We have been characterizing the environment by the zero-bias resistance of two arrays measured in series, $`R_0^{}`$. Since $`Z_tZ_1`$, it would be more appropriate to use $`R_0^{}/2`$. The factor of 2, however, is not very important in this work, and thus we keep using $`R_0^{}`$. We should note here that $`ZR_0^{}`$ is a bold approximation. The $`I`$-$`V`$ curves of the SQUID arrays are nonlinear as we have seen in Figs. 3 and 7, and in general the SQUID array is not described by a liner impedance. The frequency dependenceHaviland et al. (2000) is also neglected in the approximation.

In order to examine the $`R_0^{}`$ dependence further, we plot in Fig. 12 the zero-bias resistance $`R_0`$ of the single junction vs. $`R_0^{}`$ at three different temperatures for Sample G.

At all the temperatures $`R_0`$ is a smooth function of $`R_0^{}`$ and increases with increasing $`R_0^{}`$, which supports the notion that the approximation $`ZR_0^{}`$ is adequate. Even at the largest $`R_0^{}`$, $`R_0`$ is still increasing with no sign of saturation. The top left point in Fig. 12 for $`T=0.04`$ K and $`R_0^{}=6.5`$ M$`\mathrm{\Omega }`$ ($`f=0.46`$) corresponds to the bottom $`I`$-$`V`$ curve in Fig. 6. Thus, $`R_0`$ vs. $`R_0^{}`$ also suggests that the single junction would still be influenced by the environment at $`R_0^{}R_K`$. The change of $`R_0`$ is weak at $`R_0^{}<10^310^4`$ $`\mathrm{\Omega }`$. This fact would be understood in the following context. In the limit of $`R_0^{}0`$, the impedance of the electromagnetic environment for the Josephson junction at relevant frequencies is determined by the microwave impedance of the leads, which is usually of the order of the free-space impedance $`Z_0377`$ $`\mathrm{\Omega }`$. Thus, when $`R_0^{}`$ is of the order of $`Z_0`$ or smaller, our assumption that the SQUID arrays determine the environment is not valid, and $`ZZ_0`$ irrespective of the value of $`R_0^{}`$. In Ref. Steinbach et al., 2001, special care was taken in order to achieve $`Z/Z_0<0.1`$ at all relevant frequencies.

The interaction with a dissipative environment, or the dissipative dynamics of single Josephson junctions has also been studied theoreticallySchön and Zaikin (1990); Kato and Imada (2000); Chung (2002) in the circuit shown in Fig. 11b with $`Z(\omega )=R`$. Note that the two circuits in Fig. 11 are equivalent when $`I=V/Z(\omega )`$. As the dissipation is increased, i.e., $`R`$ is decreased, the Josephson junction undergoes a transition from an insulator to a superconductor. The phase diagram for this SI transition has been derivedSchön and Zaikin (1990) and qualitatively supported by the experiments.Yagi et al. (1997); Penttilä et al. (2001) The phase (“superconductor” or “insulator”) is usually determined by the $`I`$-$`V`$ curve (or $`dV/dI`$ vs. $`I`$ curvePenttilä et al. (2001)) at the lowest temperatures and/or the temperature dependence of the zero-bias resistance.Yagi et al. (1997) We also find evidence of the SI transition in the nonlinearity of the $`I`$-$`V`$ curve. For Sample G, shown in Fig. 6, we see that the $`I`$-$`V`$ curve is “Josephson-like” when the leads have $`R_0^{}=0.21`$ k$`\mathrm{\Omega }`$ ($`f=0`$), whereas at $`R_0^{}0.33`$ M$`\mathrm{\Omega }`$ ($`f0.45`$), there is a distinct Coulomb blockade. The temperature dependence of $`R_0`$ also supports the existence of the SI transition. We see in Fig. 12 that at $`R_0^{}<10^310^4`$ $`\mathrm{\Omega }`$, $`R_0`$ has little temperature dependence, which is similar to the behavior of a 1D systemChow et al. (1998); Haviland et al. (2000, 2001) in the superconducting phase. At $`R_0^{}>10^410^5`$ $`\mathrm{\Omega }`$, on the other hand, we see insulating behavior, where $`R_0`$ increases rapidly as the temperature is lowered.

We show in Fig. 13 the phase diagram for the SI transition determined by the $`I`$-$`V`$ curves at the lowest temperatures of all the samples studied in this work.

Note that $`R_0^{}`$ was varied over as large as nine orders of magnitude in this work. In Samples F–H, we observed the SI transition. The squares and the solid lines to the left of the squares indicate that insulating behavior (Coulomb blockade or “Coulomb-blockade-like” $`I`$-$`V`$ curve) was observed. The circles and the solid lines to the right of the circles represent the superconducting phase (“Josephson-like” $`I`$-$`V`$ curve). In the range denoted by broken lines, the $`I`$-$`V`$ curve was almost linear, and it was hard to determine the phase. For Samples I–L (Sample E), the insulating (superconducting) phase alone was identified. Samples A–D have a large $`E_J/E_C`$ ($`5`$), and we expect from Fig. 8 that their $`V_b`$ is too small ($`1`$ $`\mu `$V) to be detected in our dc measurements, i.e., the insulating phase is never identified. Indeed, the $`I`$-$`V`$ curve of Samples A–D is almost independent of $`R_0^{}`$. As an example, we show in Fig. 14 the $`I`$-$`V`$ curve of Sample D at $`R_0^{}=0.14`$ k$`\mathrm{\Omega }`$ ($`f=0`$) and $`R_0^{}=2`$ M$`\mathrm{\Omega }`$ ($`f=0.5`$).

It is Josephson-like and in a contrast with the curve in Fig. 5. We do not discuss the SI transition for Samples A–D, and just indicate the $`R_0^{}`$ range covered in the measurements by dot-dashed lines in Fig. 13. According to the theorySchön and Zaikin (1990) of single Josephson junctions with an Ohmic shunt $`R`$, the phase boundary is located at $`R_Q/R=1`$ irrespective of the ratio $`E_J/E_C`$ when the quasiparticle tunneling is not important ($`\alpha 1`$) and the transition is driven solely by $`R`$. In Samples F–L, $`\alpha =10^410^21`$ from Table 1. It is interesting to see that Fig. 13, the phase diagram for our system with non-Ohmic environment, is consistent with the theoretical phase diagram where an Ohmic shunt is assumed.

## V Conclusion

We have studied the current-voltage ($`I`$-$`V`$) characteristics of single Josephson junctions biased with a tunable electromagnetic environment composed of SQUID arrays. We have demonstrated that the single junction is indeed sensitive to the state of the environment, and a Coulomb blockade is induced in the single junction by increasing the zero-bias resistance $`R_0^{}`$ of the SQUID arrays. When $`R_0^{}`$ is much higher than the quantum resistance $`R_Kh/e^226`$ k$`\mathrm{\Omega }`$, a region of negative differential resistance has also been observed in the $`I`$-$`V`$ curve of the single junction. The negative differential resistance is evidence of coherent single-Cooper-pair tunneling according to the theory of current-biased single Josephson junctions. Based on the theory we have calculated the $`I`$-$`V`$ curves numerically. The calculation has reproduced the measured $`I`$-$`V`$ curves at $`R_0^{}R_K`$ qualitatively. We have also discussed the superconductor-insulator (SI) transition in single Josephson junction driven by the environment. We have characterized the environment by $`R_0^{}`$ and determined the phase diagram for the transition. This phase diagram is consistent with the theoretical one for the SI transition in single Josephson junctions driven by an Ohmic shunt alone.

###### Acknowledgements.

We are grateful to R. L. Kautz for great help in the numerical calculation, and to T. Kato and F. W. J. Hekking for fruitful discussions. This work was supported by Swedish NFR, and Special Postdoctoral Researchers Program and President’s Special Research Grant of RIKEN. M. W. would like to thank the Japan Society for the Promotion of Science (JSPS) and the Swedish Institute (SI) for financial support.


format: text

text:
# Melting of ‘porous’ vortex matter

## Abstract

Bitter decoration and magneto-optical studies reveal that in heavy-ion irradiated superconductors, a ‘porous’ vortex matter is formed when vortices outnumber columnar defects (CDs). In this state ordered vortex crystallites are embedded in the ‘pores’ of a rigid matrix of vortices pinned on CDs. The crystallites melt through a first-order transition while the matrix remains solid. The melting temperature increases with density of CDs and eventually turns into a continuous transition. At high temperatures a sharp kink in the melting line is found, signaling an abrupt change from crystallite melting to melting of the rigid matrix.

Melting of heterogeneous systems, and in particular of nanocrystals embedded in porous rigid matrices, is a complex process with many uncontrolled parameters. Metal and semiconductor nanocrystals with free surfaces, for example, usually show a decrease in their melting temperature with decreasing size , whereas nanocrystals encapsulated in a porous matrix often display an increase in melting temperature . Although the contribution of the different factors is still a matter of debate, the melting process is known to depend on the size, dimensionality, material properties of the nanocrystals and the matrix, as well as the interface energies between the materials . In this work we investigate an analogous, but a more controllable composite system, which is a ‘porous’ vortex matter consisting of vortex nanocrystals encapsulated in a matrix of strongly pinned vortices. As shown below, this system is present in the commonly heavy-ion irradiated superconductors when the vortices outnumber the columnar defects (CDs). The rigid matrix is created by vortices localized on the network of random CDs, while the softer nanocrystals are formed within the ‘pores’ of this matrix by the interstitial vortices. The size of the nanocrystals can be readily varied from several hundred down to a few vortices by changing the applied field or the density of CDs. We find that this composite vortex matter reveals a number of intriguing mechanisms: Similarly to the metallic nanocrystals in a matrix, we observe for the first time a pronounced upward shift in the vortex melting temperature $`T_m`$, while *preserving* the first-order nature of the transition (FOT). With increasing density of CDs, the size of the pores decreases, resulting in a larger shift in $`T_m`$. We also find a critical point at which the FOT changes into a continuous melting. Moreover, the crystallites can melt while the matrix remains rigid. As a result, at high temperatures we find an abrupt breakdown in the upward shift of $`T_m`$ and a sharp kink in the FOT line, which apparently result from the collapse of the matrix due to vortex depinning from the CDs.

The reported findings were obtained using Bitter decoration and differential magneto-optical (MO)  techniques. High quality Bi2Sr2CaCu2O8 (BSCCO) crystals ($`T_c89`$ K) were covered by various patterned masks and irradiated at GANIL by 1 GeV Pb ions with doses corresponding to matching fields of $`B_\varphi =`$ 5, 10, 20, and 50 G. Figure 1a shows schematically one of these masks which results in the formation of CDs in BSCCO crystals only within the circular apertures of about 90 $`\mu `$m diameter. This patterning allows very sensitive simultaneous comparison of the vortex structure and the local melting processes in adjacent irradiated and pristine regions, which is not possible by other methods.

In the absence of CDs the vortices form the Bragg glass phase  which has quasi-long-range order with no topological defects, as seen in the pristine part of the magnetic decoration image in Figs. 1b and 1c. What happens to this phase when sparse CDs are added? The irradiated region in Fig. 1b shows that it is no longer Bragg glass since it has significant amount of topological defects (solid circles in Fig. 1c) and no orientational long-range order . It is also not an amorphous or glass phase in the usual sense, nor is it a simple polycrystal as discussed below. In the presence of CDs the Bose glass (BG) theory  is usually applied, which describes the vortex matter in terms of anisotropic *homogeneously* pinned medium. Such a description is adequate for the common situation of high irradiation doses ($`B_\varphi >B`$), in which all the vortices reside on CDs and the vortex pinning energies are comparable. In contrast, we investigate here mainly the opposite extreme of low doses in which the vortices greatly outnumber the CDs, i.e., $`BB_\varphi `$. In this case the system is inherently *heterogeneous*, consisting of two vortex populations with *well separated characteristic energies*: The vortices residing on CDs are strongly pinned and form a rigid network or matrix, whereas the interstitial vortices are localized by significantly weaker elastic interactions and form relatively soft crystallites within the pores of the matrix. In the following we refer to this state (upper parts of Figs. 1b and 1c) as porous vortex matter in order to emphasize the important consequences of the heterogeneous structure.

Figure 2 shows several frames from a ‘movie’  of the melting process as a function of temperature $`T`$ at two fields, 30 and 60 G. Each frame is obtained by taking the difference between the MO images at $`T+0.15`$ K and $`T0.15`$ K and averaging a large number of such differential images, as described previously . The bright features show the regions in the sample that undergo a FOT within the temperature interval of 0.3 K at the indicated $`T`$, and the intensity of this bright paramagnetic signal is the equilibrium magnetization step $`\mathrm{\Delta }B`$ at the transition . Figure 2a shows the nucleation of the liquid phase in the form of bright inclined strips in the central pristine region which arise from intrinsic sample disorder (see Ref. ). With increasing $`T`$ (Fig. 2b) the liquid expands, remarkably *avoiding* the irradiated apertures. In Fig. 2c the entire central pristine part of the sample is liquid, while the apertures with $`B_\varphi =20`$ G are still solid. In Fig. 2d the central apertures melt at 82.45 K, which is about 1 K above $`T_m`$ of the adjacent surrounding pristine regions in Fig. 2b. The apertures closer to the sample edges begin to melt in Fig. 2e. The FOT is equally strong in the irradiated and pristine regions: The $`\mathrm{\Delta }B`$ step derived from the paramagnetic melting signal  is the same in Figs. 2d and 2b. Also the width of the local melting transition is the same within our resolution, i.e., each point in the sample melts within 0.3 K or less. This is the first direct observation of an upward shift of the FOT by correlated disorder. Note that for small $`B_\varphi `$ the melting in the irradiated apertures occurs while the melting in remote pristine regions is still in progress, and therefore the quantitative measurement of the shift in $`T_m`$ has to be determined by comparing neighboring irradiated and pristine regions and cannot be readily detected by global techniques.

The melting process at 60 G (Fig. 2, second row) reveals two important differences. First, the shift of the melting temperature, $`\mathrm{\Delta }T_m`$, is about 4 K (difference between Figs. 2i and 2f) which is much larger as compared to 1 K at 30 G. In Fig. 2h, for example, the entire pristine sample has melted while the irradiated apertures are still solid. Second, the brightness of the paramagnetic melting signal $`\mathrm{\Delta }B`$ in the apertures in Fig. 2i is much lower than in the pristine sample, and moreover, the melting in each aperture is broadened over several frames, as exemplified in Fig. 2j. At still higher fields, above 100 G, no paramagnetic FOT signal is detected in the irradiated apertures with $`B_\varphi =20`$ G.

In order to investigate the phase transition over a wider range of fields we have used differential MO imaging in which the field is modulated by 1 G (Fig. 3) instead of the $`T`$ modulation. In addition to detecting the FOT , this method provides a very sensitive measurement of the irreversibility line at a very low effective frequency of about 0.1 Hz. Figure 3 demonstrates the determination of the irreversibility line in the region where no FOT is present. In Fig. 3a the entire pristine sample is in the liquid state, while the apertures are still solid. The apertures appear black, showing that the external field modulation is shielded due to the enhanced pinning. Upon increasing the field the black apertures disappear (Figs. 3b and 3c) revealing the value of the local irreversibility field.

Figure 4 shows the location of the onset of the FOT for $`B_\varphi =5`$, 10, 20, and 50 G obtained by $`T`$ modulation (solid symbols) and of the irreversibility line obtained by field modulation (open symbols). The solid lines are guides to the eye for the FOT lines which terminate at the corresponding critical points. The irreversibility data coincide with the FOT line below the critical point and smoothly extrapolate the location of the transition line to higher fields. The first interesting observation here is that although the structures of the porous vortex matter and of the Bragg glass are very different (Fig. 1), their phase diagrams for $`B_\varphi =`$ 5 G in Fig. 4 are almost identical. The melting remains FOT in most of the field range and the melting line is shifted only slightly. This brings us to an important conclusion that the quasi-long-range order that characterizes the Bragg glass is not an essential requirement for the existence of a FOT , and the presence of the short-range order within the crystallites is apparently sufficient.

We can understand the upward shift in $`B_m(T)`$ in Fig. 4 by generalizing the concept of the cage model. In a pure system each vortex is confined in a potential cage arising from the elastic interactions with its neighbors. Melting occurs when the transverse thermal fluctuations of vortices $`u_T`$ reach a certain fraction $`c_L`$ of the lattice spacing $`a_0`$. The pores of the matrix are vertical cylinders which provide an additional confining cage potential to the interstitial vortices. This enhanced rigidity reduces $`u_T`$, and hence stabilizes the solid crystallites within the pores, thereby shifting $`B_m(T)`$ upwards. The shift in $`B_m(T)`$ grows with $`B_\varphi `$ due to the decrease in the size of the pores. Since the pinning energy of the CDs is usually significantly larger than the elastic energy, the melting of the crystallites may occur without the destruction of the matrix, which remains solid up to a higher temperature as described below. The rigid matrix may thus coexist with an interstitial liquid as shown theoretically  and in numerical simulations .

Upon increasing the density of CDs the FOT is weakened and it eventually transforms into a continuous transition in Fig. 4. This transformation is seen more clearly by varying the temperature along a given $`B_m(T)`$ line. The inset to Fig. 4, left axis, shows the height of the equilibrium magnetization step $`\mathrm{\Delta }B`$ vs. $`T`$ for $`B_\varphi =20`$ G. As the temperature is decreased $`\mathrm{\Delta }B`$ drops and eventually vanishes at a critical point (CP) below which no discontinuity is found and the melting becomes continuous. In addition, the width of the FOT, $`\delta T_m`$, shows a significant broadening on approaching the CP (Fig. 4 inset, right axis). This indicates that it is a true CP at which the FOT transforms into a continuous transition as expected theoretically . This CP is unlike the apparent point-disorder-induced CP in BSCCO at lower $`T`$ where no appreciable broadening was observed , and where the FOT was found to be obscured by pinning-induced irreversibility in the liquid phase . Here, in contrast, the liquid is fully reversible and no such experimental limitations exist. Figure 4 shows that the CP is shifted to lower fields with $`B_\varphi `$ . It is interesting to note that in YBCO the CP was observed to shift upward with CDs . There, however, the CP is induced by the point disorder in the pristine crystals and the CDs just modify slightly its position by reducing the point-disorder-induced vortex meandering. In our case, in contrast, the CDs apparently create a new critical point that does not exist in pristine crystals.

We now discuss another key finding which is a kink in the $`B_m(T)`$ lines marked by the arrows in Fig. 4. This kink becomes very prominent when the temperature shift $`\mathrm{\Delta }T_m`$ between the irradiated and pristine melting lines is plotted as shown in Fig. 5. The inset displays the corresponding upward shift in field $`\mathrm{\Delta }B_m`$. In BG theory a kink in the transition line is expected to occur at $`B_kB_\varphi `$ . We argue that the kink in Fig. 5 is a new feature of completely different nature, which reflects the collapse of the rigid matrix. The BG kink is experimentally found to occur at $`B_k/B_\varphi `$ of 1/6 to 1 for large $`B_\varphi `$, and $`B_k`$ usually scales linearly with $`B_\varphi `$ . The kink in Fig. 5 occurs in the opposite regime of $`B_k/B_\varphi 2`$ to 8, and it does not scale with $`B_\varphi `$. At high $`B_\varphi `$ the BG kink is a broad feature that was found to occur only along a glass transition line . In contrast, here the kink is extremely sharp, and it is the first observation of a kink that occurs along a FOT line. Finally, the BG kink is a transition from a strong influence of the CDs below the kink to a much weaker effect of CDs above the kink . Here, in contrast, the situation is just the opposite, as seen in the inset to Fig. 5. The effect of CDs is large at fields above the kink, where the shift $`\mathrm{\Delta }B_m`$ is almost constant and approximately equal to $`B_\varphi `$, and it collapses rapidly at $`T>T_k`$. This collapse occurs at $`T_k75`$ K, with little dependence on $`B_\varphi `$. It is also interesting to note that for $`B_\varphi 20`$ G the kink occurs along the FOT line, while for $`B_\varphi =50`$ G it falls in the region of continuous transition.

In Fig. 4 three regions are enumerated with respect to $`B_\varphi =50`$ G curve. In region 1 the crystallites in the pores are stabilized by the rigid matrix and melt at $`B_m^{por}`$, well above the pristine $`B_m(T)`$. In region 2 the pores are liquid while the matrix remains intact, and hence $`B_m^{por}`$ reflects a melting transition of the softer of the two substances in a heterogeneous medium. This unconventional melting cannot be described by a single set of parameters since it depends on the properties of both the pores and the matrix, each described by a different set of parameters with different field and temperature dependencies. As a result, $`B_m^{por}(T)`$ does not extrapolate to $`T_c`$ as commonly expected, but rather well above it. To the best of our knowledge this is the first observation of a melting or irreversibility line that has such an unusual property at high temperatures. This heterogeneous process collapses sharply at $`T_k`$, above which the matrix rapidly delocalizes, resulting in a homogeneous liquid in region 3. The $`B_m^{mtx}`$ line thus describes the melting process of the matrix due to depinning from the CDs, which leads to immediate melting of the superheated crystallites. Thus in contrast to the common homogeneous BG description, the $`B_m^{por}(T)`$ and $`B_m^{mtx}(T)`$ lines originate from *two independent heterogeneous processes*, melting of crystallites within a rigid matrix and the collapse of the matrix itself, resulting in a sharp kink at the intersection point. We expect the $`B_m^{mtx}(T)`$ line to extend also into the liquid phase , like the dashed line in Fig. 4, separating the interstitial liquid within a solid matrix in region 2 from the homogeneous liquid 3. Our MO measurements cannot detect this line since both these regions are fully reversible. In region 2, however, we expect a higher c-axis correlation than in region 3. Recent Josephson plasma resonance studies of BSCCO with low $`B_\varphi `$ indeed find a recouping transition within the liquid phase at which an enhancement in the interlayer phase coherence is observed . This extension of the $`B_m^{mtx}`$ line should also be detectable by transport measurements since the flux-flow resistance of the interstitial liquid should be lower than of the homogeneous liquid . These studies will be the subject of future work.

In summary, when vortices outnumber CDs, heterogeneity rather than the average properties of the lattice has to be taken into account for a proper description of the structure and the thermodynamic behavior of the vortex matter. We have found evidence for two mechanisms: melting of superheated crystallites within the pores of a solid matrix and the destruction of the rigid matrix. The intersection point of these two independent processes results in a sharp kink in the observed melting line. The heterogeneous melting can be either first-order or continuous, depending on temperature and the density of CDs. The porous vortex matter may thus provide a tunable model system for general comprehension of melting of nanocrystals in porous materials.

This work was supported by the Israel Science Foundation and Center of Excellence Program, by the German-Israeli Foundation G.I.F., by the Minerva Foundation, Germany, and by the Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture, Japan. EZ and FC acknowledge the support by the Fundacion Antorchas - WIS collaboration program.

FIGURE CAPTIONS

Fig. 1. (a) Schematic of BSCCO crystal irradiated through a mask with an array of circular apertures. (b) Bitter decoration image ($`B=40`$ G, $`T=4.2`$ K) showing pristine region at the bottom (Bragg glass) and a section of the irradiated aperture ($`B_\varphi =10`$ G) on top where porous vortex matter consisting of ordered crystallites embedded in a rigid matrix is formed. (c) Corresponding locations of six-fold (open circles) and of five- and seven-fold coordinated vortices (solid circles) obtained by Delaunay triangulation.

Fig. 2. Melting process as a function of $`T`$ at fields of 30 and 60 G in BSCCO crystal $`0.65\times 0.45\times 0.01`$ mm3 using differential MO imaging with $`T`$ modulation of 0.3 K. The bright regions are the areas that undergo melting within the 0.3 K interval. (a) $``$ (c) and (f) $``$ (h): melting of the pristine regions while the irradiated apertures are still solid. (d), (e), (i), and (j): melting of the irradiated apertures ($`B_\varphi =20`$ G). The area outside the crystal is blackened. Color movies of the melting process are available at http://www.weizmann.ac.il/home/fnsup/.

Fig. 3. Differential MO images using field modulation of 1 G in BSCCO crystal of Fig. 2 at three fields at $`T=68.4`$ K. (a) All the pristine regions are liquid while the irradiated apertures ($`B_\varphi =20`$ G) are still solid and irreversible, and hence appear black due to shielding. (b) Partial reversibility of the central apertures. (c) Central apertures are liquid while the apertures closer to the edges begin to melt.

Fig. 4. The melting lines $`B_m(T)`$ of the pristine and irradiated regions with indicated $`B_\varphi `$. Solid (open) symbols are temperature (field) modulation data showing the location of the FOT (irreversibility line). Solid (dotted) lines are guides to the eye of the first-order (continuous) transitions that terminate at the critical points $``$. Inset: The height of the FOT equilibrium magnetization step $`\mathrm{\Delta }B`$ which vanishes at the CP, and the local width $`\delta T_m`$ of the FOT vs. $`T`$.

Fig. 5. The shift $`\mathrm{\Delta }T_m`$ in the melting temperature for different $`B_\varphi `$ with respect to the pristine $`T_m`$ at various fields $`B`$. Inset: The upward shift $`\mathrm{\Delta }B_m`$ in the melting field vs. $`T_m`$.


text:
# Ram pressure stripping in Planetary Nebulae

## 1 INTRODUCTION

The process or processes that may cause the departure of spherical symmetry in PNs have been the subject of many studies e.g. the hydrodynamic collimation of the fast wind by an equatorial density enhancement (Frank, & Mellema, 1994), the rotation of a single star with a magnetized wind (Różyczka & Franco, 1996; García-Segura, 1997; García-Segura et al., 1999) and the evolution of a PN in a close binary system (Soker, 1996). These mechanisms are able to account for the formation of elliptical, point-symmetric and bipolar PNs. However, several PNs have been observed in which mainly the external shell shows departure from sphericity (e.g. MA 3, PM 1-295, M 2-40). In these PNs, the presence of asymmetries cannot be explained by the physical processes related to the central star (CS), such as rotation, magnetic fields or binarity. Instead, the interaction with the ISM was proposed to be the cause of the observed asymmetries.

The interaction of a PN with the ISM was first suggested by Gurzadyan (1969) and the first theoretical study was made by Smith (1976); he assumed a thin shell approximation and the ‘snow-plough’ momentum-conserving model of Oort (1951). Later on, Isaacmann (1979) used the same approximation, but considering higher velocities and ISM densities. Both theoretical works arrived at a similar conclusion, namely that the nebula fades away before any disruption of the nebular shell becomes noticeable. Borkowski, Sarazin, & Soker (1990) analytically, and Soker, Borkowski, & Sarazin (1991) numerically confirmed the thin-shell approximation originally developed by Smith (1976), concluding that the interaction becomes dominant when the density of the nebular shell drops below a critical limit, due to the nebular expansion. As a consequence, the interaction should be easily observable during the late stages of the nebular evolution. Moreover, high ISM densities, large relative velocities or both were needed to explain PN asymmetries due to interaction with the ISM.

This is the scenario commonly accepted nowadays and as a result, the search for interaction of PNs with the ISM has been mostly restricted to old PNs with large angular extent (Tweedy & Kwitter, 1996; Xilouris et al., 1996). Evidence for the interaction has been also reported in the study of some individual objects (e.g. Borkowski, Tsvetanov, & Harrington 1993; Tweedy, Martos, & Noriega-Crespo 1995; Soker & Zucker 1997) and the observations have been interpreted in the light of the available models. As a consequence high spatial velocities were invoked to explain the morphologies of PN at large distances from the Galactic plane.

The problem is that in the available models, the PN-ISM interaction was studied by considering the relative movement when nebular shell was already formed. In this letter, we demonstrate that when the PN formation is followed by considering the evolution of the stellar wind from the AGB, high velocities, high ISM densities or the presence of a magnetic field are not necessary to explain the observed asymmetries in the external shells of PNs. Moreover, in contrast to what is commonly believed, the interaction is observable even during the very early stages of the PN phase.

## 2 THE NUMERICAL METHOD

We have carried out the numerical simulations with the fluid solver ZEUS-3D (Stone & Norman 1992a, b; Stone, Mihalas, & Norman 1992), developed by M. L. Norman and the Laboratory for Computational Astrophysics. The computations have been performed using a 2D spherical polar grid with the angular coordinate ranging from 0 to 180 and a physical radial extension of 2 pc which gives a total grid size of 4$`\times `$2 pc. We have adopted a resolution of 400$`\times `$360 zones in the radial and angular coordinates of the grid respectively. We consider the photoionization of the gas during the PN phase by using a simple approximation to derive the location of the ionization front for arbitrary density distributions (see Bodenheimer, Tenorio-Tagle, & Yorke 1979; Franco, Tenorio-Tagle, & Bodenheimer 1990). This is done by assuming that ionization equilibrium holds at all times and that the gas is fully ionized inside the H II region. The models include the Raymond & Smith (1977) cooling curve above 104 K. For temperatures below 104 K, the unperturbed gas is treated adiabatically but the shocked gas region is allowed to cool down with the radiative cooling curves given by Dalgarno & McCray (1972) and MacDonald & Bailey (1981). The photoionized gas is always kept at 104 K, so no cooling curve is applied to the H II region unless there is a shock.

### 2.1 The Boundary Conditions

We assume that the ISM moves relative to the star perpendicular to the line of sight by fixing the position of the star at the center of the grid and allowing the ISM to flow into the grid at the outer boundary from 0 to 90. On the other hand, from 90 to 180 we set an outflow boundary condition. The temporal evolution of the stellar wind for a 1 M star during the AGB phase (taken from Vassiliadis & Wood 1993) and the PN stage (Vassiliadis & Wood, 1994) has been set within a small (five radial zones) spherical region centered on the symmetry axis, where reflecting boundary conditions are used. The gas evolution in a static configuration as well as the stellar wind inputs for the models are described in detail in Villaver, García-Segura, & Manchado (2002a) and Villaver, Manchado & García-Segura (2002b).

In order to study the interaction process we have adopted a purely hydrodynamical scheme, in which we assume that the ISM pressure is simply the standard gas thermal pressure, $`P=nKT`$ (where n is the number of particles per unit volume, K is the Boltzman’s constant and T is the gas temperature). Since the strength of the interaction depends linearly on the ISM density and quadratically on the relative velocity of the star, we have adopted two ISM densities and a very conservative value for the star’s velocity, 20 km s-1. The ISM is assumed to be homogeneous with the characteristics of the warm neutral medium (WNM) which is the main observed constituent of the ISM, with temperatures estimated to be in the range 5000-8000 K (Kulkarni & Heiles, 1988) and neutral atomic hydrogen densities between 0.1 and 1 $`\mathrm{cm}^3`$ (Burton, 1988). We use an ISM temperature of 6000 K and densities of 0.1 and 0.01 $`\mathrm{cm}^3`$, the latter to study the interaction with a very low density ISM.

## 3 RESULTS

### 3.1 The Evolution during the AGB through a Low Density Medium

The logarithm of the gas density during the AGB phase is shown in Figure 1 for a star moving with 20 km s-1 through a low density ISM (no = 0.01 $`\mathrm{cm}^3`$). Figure 1 shows only half of the $`r\theta `$ plane, the star is fixed in the grid and the ISM flows in from the top to the bottom.

During the very early stages of the evolution of an AGB star a free-streaming steady stellar wind with a velocity of 2 km s-1 and mass-loss rate of Ṁ = 10-8 M yr-1 reaches ram pressure equilibrium with the ISM at a stand-off distance from the star of 0.13 pc (first panel on the left of Fig. 1). After this, the mass-loss rate and wind velocity are changing continuously in the inner boundary giving rise to the formation of shocks that develop inside the bow-shock cavity formed by the early stationary wind. From then onwards the stand-off distance cannot be computed from a simple ram pressure balance argument between a free streaming stellar wind and the ISM. The ram pressure of the ISM is balanced by the ram pressure of the stellar wind inside the bow-shock cavity, making it a time dependent problem.

As a consequence of the interaction process a highly asymmetric shell is formed that grows in size mainly in the downstream direction. The outer bow shock is radiative for the stellar velocity considered, so the density of the swept-up shell is proportional to the square of the upstream Mach number and proportional to the ISM density. As the mass-loss rate continues (highly modulated by the thermal pulses in the star, see Fig. 1 of Villaver et al. 2002a) new shells are subsequently formed inside the structure generated by the previous stellar wind. During most of the AGB evolution (355 000 yr of the $``$ 495 000 yr that the AGB phase lasts) the leading part of the gas is not able to expand far away from the CS. Only the high mass-loss rate periods associated with the last two thermal pulses during the AGB have enough momentum to compete efficiently with the previously developed density structure. It is at this moment that the leading parts of the outer shell growin size. The last two panels show the density at the end of the second thermal pulse and in the middle of the third one respectively. Because of the relative movement, mass is constantly flowing away from the head of the bow shock forming cometary structures behind the star.

### 3.2 The Evolution during the AGB through a Higher Density Medium

In Figure 2 we show the gas evolution for an ISM density of no = 0.1 cm-3, lower than the value in the Galactic plane ($`1`$ $`\mathrm{cm}^3`$). The snapshots of the evolution have been selected at the same times as those shown in Figure 1. As expected, the effects of the stellar movement are more noticeable, since the ram pressure exerted by the ISM depends linearly on the ISM density. The compression and deformation of the upstream side of this shell is higher than in the case shown in Figure 1. As for the low density case, the stagnation radius can be computed from a ram pressure balance argument only when the wind is stationary in the early evolution (see the first panel of Fig. 2 where the stagnation radius is a factor 10-1/2 smaller than in the first panel of Fig. 1). In the present case, the upstream side of the shell does not increase in size until very late in the evolution and the size is always at least 1.5 times smaller than the one formed when the star evolves in a lower density medium.

### 3.3 The PN stage

During the PN phase the stellar wind increases its velocity, generating an adiabatic shock that forms the brightest shell of the nebula. In principle, this shell should not be affected by the interaction with the ISM, unless other mechanisms such as ambipolar diffusion under the presence of an external magnetic field are invoked to propagate the asymmetry to the innermost regions of the nebula.

In the right panel of Figure 3 we show the gas density when the PN is 2 500 yr old111That is 2 500 yr after the photoionization has started at the end of the AGB phase. and when the ISM has a density of 0.1 cm-3. The simulated region has been reflected around the symmetry axis in order to more clearly show the structure. The left panel shows the PN at the same age, but when the star is evolving at rest with respect to the ISM (see Villaver et al. 2002b). Two bright shells are clearly distinguished in both panels of Figure 3. These are a spherical innermost shell shaped by the fast wind and an external shell which shows a characteristic bow-shock configuration in the right panel Figure 3.

The first effect of the interaction, apart from that on the morphology, is that the total size of the outer PN shell (halo) is reduced considerably. The radius of the halo in the left panel of Figure 3 is 1.8 pc, whereas the outer shell of the interacting nebula has a radius of 0.7 pc. The displacement of the CS with respect to the geometrical center of the nebula is 0.1 pc, which corresponds to 20 ″ at a distance of 1 kpc. A less obvious effect is that the formation of the outer shell of the PN takes place in a lower density environment, since most of the mass ejected during the AGB phase is stripped by the ram pressure of the ISM and left in the downstream direction of the stellar movement. The model at rest has 0.7 M of gas above the ISM value in the circumstellar envelope (0.43 M that were lost by the star and 0.27 M of ISM matter swept up by the wind), only 0.24 M of gas are recovered when the star is moving.

## 4 DISCUSSION

We find that the movement of the CS with respect to its surrounding medium considerably alters the circumstellar gas structure formed during the AGB phase and hence the PN formation. Up to now, it has been generally accepted that the deceleration of the nebular shell due to its interaction with the ISM can be more likely observed during the late stages of the nebular evolution. We find that when the evolution of the stellar wind is properly considered in the simulations, the interaction with the ISM cannot be studied using simple ram pressure balance arguments and plays a major role even during the early AGB evolution.

From an observational point of view, the presence of asymmetries in the haloes of PN that could be related to the interaction with the ISM is a common feature222Multiple shells appear in 24% of the spherical and elliptical PN sample of the northern hemisphere (Manchado et al., 1996); 40% of those show asymmetries in the haloes that could be related to the interaction with the ISM.. The high rate of observed asymmetries cannot be explained if all of these PNs are moving at high velocities through high density media. Instead, it can be understood when the evolution of the stellar wind is taken into account, producing a much higher efficiency in the interaction. Moreover, in agreement with our results, Guerrero, Villaver, & Manchado (1998) show kinematical evidences that the interaction of the external shells of PN with the ISM are found at all evolutionary stages.

According to our simulations it may not be necessary to invoke the presence of a magnetic field in the ISM to explain the asymmetries found in the slow moving (10-20 km s-1) PNs Sh 2-216 (Tweedy, Martos, & Noriega-Crespo, 1995) and NGC 6894 (Zucker & Soker, 1997). For the model parameters assumed in this paper, no instabilities appear in the nebula, which is in agreement with the analytical results of Dgani & Soker (1994, 1998).

Another important aspect is that the interaction with the ISM considerably reduces the mass of the circumstellar envelope during the AGB and PN phases due to ram pressure stripping. Stellar evolution calculations predict that stars with initial masses in the range of $``$1–5 M will end as PN nuclei with masses around 0.6 M. Most of the mass is lost on the AGB phase and should be easily observable as ionized mass during the PN stage. However, observations of Galactic PNs reveal on average only 0.2 M of ionized gas. The mass stripped away by the ISM when the star is moving could solve the problem of the missing mass in PNs.

## 5 SUMMARY

We have studied the effects the relative movement of the CS with respect to the ISM has on the evolution of the stellar ejecta during the AGB phase, and on shaping the PN. Under the very conservative conditions assumed for the ISM density and the CS’s velocity (20 km s-1) we find that the PN-ISM interaction provides an adequate mechanism to explain the high rate of observed asymmetries in the external shells of PNs. According to our simulations, spherical haloes in PNs are expected only if the star is at rest with respect to the ISM or if the star is moving at low angles with respect to the line of sight, otherwise, bow-shock like and cometary structures should be present.

We are grateful to the referee for very constructive suggestions. We thank M. L. Norman and the Laboratory for Computational Astrophysics for the use of ZEUS-3D. EV would like to thank Tariq Shahbaz and Mónica Sánchez Cuberes for fruitful discussions. GGS is partially supported by DGAPA-UNAM grant IN117799 and CONACyT.


text:
# Distant globular clusters with anomalously small masses

(received July 29, 2002 accepted October 13, 2002)

## Abstract

We found that 10 metal-poor globular clusters are greately distinguished for anomalously small masses on the ”destruction rate–mass” plain. As it turned out, these poor clusters, situated 15 kpc farther from the Galactic centre, are somewhat younger than the bulk of the metal-poor globulars, and have anomalously red horizontal branches. All these clusters are supposed to belong to the ”young halo” subsystem, i.e. they are supposed to have been captured by the Galaxy at different stages of their evolution. We discovered a significant correlation between the ages found from isochrones and masses of globulars which lie at galactocentric distances greater than the radius of the solar orbit. At the same time, deficiency of distant massive clusters is noticeable with increasing distance from the galactic centre. So we see unticorrelation between the galactocentric distance and masses of distant clusters. Both relations are negligible for the inner clusters of the Galaxy. We assume that favourable conditions for violent dissipation with considerable loss of mass are realised inside the protoglobulars which formed far from the Galactic centre.

Van den Berg (1998, 2000) payed attention to the fact that all known to him ”young” globular clusters with ages by $`3`$ billion years younger than the oldest clusters of the same metallicity lie at galactocentric distances $`R>15`$ kpc and have masses smaller than the average. That is why he assumed that the initial mass spectra of globular clusters formed in the outer halo of the Galaxy displace with time towards smaller masses. Independently, we (Borkova, Marsakov, 2000) have found that the average mass of clusters located within the radius of the solar circle somewhat increases, while in the outer halo it decreases noticeably with growing galactocentric distance. The clusters of the outer halo have markedly smaller masses than those of the thick disk, even among the youngest objects (see Fig. 12 and Table 2 in the papers by Borkova and Marsakov, 2000). In the present paper this phenomenon is discussed in more details and an attempt to explain it is made.

As initial data we have used the compiled catalogue of the measured characteristics of the globular clusters of the Galaxy (Harris, 1996) and the summary catalogue of basic parameters for 145 clusters compiled on its basis by Borkova and Marsakov (2000) with involvement of other sources. For the present paper we have recalculated the parameters on the basis of a new version of the catalogue by Harris (1996) compiled for 1999 June 22, where the distances and luminosities of clusters of the inner region of the Galaxy were mostly revised. We found the masses from integral stellar magnitudes of clusters using a constant mass-to-light ratio $`M/L_V=3`$, where the mass $`M`$ and luminosity $`L_V`$ are expressed in solar units. It should be borne in mind, however, that the constancy of the ratio $`M/L_V`$ is not satisfied in all the cases. In particular, the drop in luminosity of clusters with a given mass will be observed in the case of poor population of their red giant branch. Similarly, the mass of a cluster distant from the galactic centre may be underestimated due to a considerable fraction of small-mass stars preserved from tidal destruction. Fig. 1a shows a relationship between the galactocentric distance and the mass of clusters. The dashed line restricts the region of slow evolution of globular clusters the existence of which is theoretically grounded in the paper by Surdin and Arkhipova (1998). The upper line corresponds to the critical value of mass resistent to the effects of dynamical friction that leads to deceleration of a massive stellar cluster moving through the field stars and to its destruction at the centre of the Galaxy under the action of tidal forces, while the lower one corresponds to the effects of dissipation from tidal ”shocks” when the cluster flies through the galactic plane. It can be stated with high probability that the clusters lying on the diagram outside of this region are at the end of their lifetime. The solid lines in Fig. 1a represent straight regressions for the regions inside and outside the solar circle ($`R_{}=8.5`$ kpc). The corresponding correlation coefficients are shown in the figures. It can be seen that in the inner region the inspected data do not reveal variation of the average mass of clusters with increasing distance from the galactic centre. However, in the outer region of the Galaxy the observed anticorrelation is different from zero beyond the errors ($`r=0.3\pm 0.1`$). The application of the $`t`$ criterion has shown that correlation between the mass and the galactocentric distance of globular clusters in the outer halo takes place at a confidence level of about 95 % despite the small value of correlation coefficient. The slope of the regression does not decrease, and the correlation coefficient remains beyond the errors even if six the most remote points are excluded from the diagram. According to Surdin and Arkhipova (1998), such distant clusters did not undergo a considerable dissipation and dynamical friction. This is why the initial mass distribution in the outer halo remained almost unaffected. Thus, a question arises why the deficiency of massive globular clusters increases in the outer halo with increasing galactocentric distance (the lower right angle in the diagram is practically empty).

The present-day theory of dynamical evolution of globular clusters offers a possibility of estimating the mass loss rate of a large number of clusters. In particular, Gnedin and Ostriker (1997) calculated the destruction rate of 119 clusters of galaxies caused by the combined action of two–body relaxation, tidal destruction, and shocks with the disk and bulge. The authors disregarded the energy increase of the clusters caused by their interaction with giant molecular clouds, following the conclusions made by Chernoff et al. (1986). The tangential components of velocities of each cluster, which were lacking for the calculation of the galactic orbits were derived by Gnedin and Ostriker (1997) statistically from the kinematic model for the galactic globular cluster system. In the recently published paper by Dinesku et al. (1999) the destruction rate for 38 globular clusters with the measured tangential and radial components of space velocities were computed. On the whole, the authors note a good agreement with the results of Gnedin and Ostriker (1997). That is why, in order to use as many clusters as possible, we considered to make use of the calculation results of destruction rates derived by Gnedin and Ostriker (1997). Fig. 1b presents the relation between the mass of clusters and the total destruction rates $`\nu `$. The light symbols denote the clusters situated at a distance of $`R>15`$ kpc from the galactic centre. It can be seen that a sufficiently reliable decrease of the average mass with increasing destruction rate is observed. However, rather a narrow sequence consisting of 10 distant clusters parallel to the main one separates sharply. Their masses are approximately an order of magnitude smaller than the average masses of clusters at the same values of $`\nu `$. The 10 clusters are shown in Fig. 1 by open triangles. Only one cluster with an anomalously small mass (Pal 3) fell in the sample of Dinesku et al. (1999), and in the diagram destruction rate — mass, plotted from the data of this paper, it lies lower than the main relation by the same value.

It is interesting that three of the isolated clusters of anomalously small mass are usually believed to belong to the Sagittarius dwarf spheroidal galaxy. These are the clusters Arp 2, Ter 7 and Pal 12. They are closely spaced and located at a distance of about 18 kpc from the centre of the Galaxy. Apart from them, there are located another two clusters from the same galaxy, Ter 8 and NGC 6715. They did not fall into our group of clusters of anomalously small mass. Although the mass of Ter 8 is the same as that of the three isolated clusters, but its destruction rate has not been found. The mass of the cluster NGC 6715 is two orders of magnitude larger. A group of five clusters with $`R100`$ kpc and $`log(M/M_{})4.5`$ is distinguished in Fig. 1, but in fact, they are spaced by a few dozen kiloparsec. Such a distant, but massive cluster NGC 2419 has a very low central density which is more characteristic of small-mass clusters, ($`log(\rho )=1.5`$), with an average density logarithm for all clusters of the Galaxy approximately equal to three. For the clusters of anomalously small mass that we have selected this value is $`0.2`$.

Check the ages of the clusters included in this group and the existence of clusters with masses above the average among the ”young” clusters in the entire sample. It is of interest to see if there is any time variation of the average mass of the clusters inside and outside the solar circle separately. For this purpose, we make use of the ages from the paper by Borkova, Marsakov (2000) where, based on the ages of 47 sources taken from literature and 336 individual determinations reduced to the unified scale, the weighted mean estimates of ages for 63 globular clusters are computed. In averaging, a two-stage iteration procedure with assigning weight both to each of the used source of ages of the clusters and to each individual determination of individual age. The internal accuracy of the obtained estimates is equal to $`\pm 0.9`$ billion years. In the present paper, in addition to the original list, we have computed the ages of another 14 clusters. The diagrams $`tlog(M/M_{})`$ are displayed in Fig. 2. Mass variation as a function of age is not observed in the inner region of the Galaxy. The correlation coefficient within the errors is equal to zero, see the figures in Fig. 2a. In the outer region the correlation turns out to be rather high ($`r=0.6\pm 0.1`$), it remains significant ($`r=0.4\pm 0.1`$) even if four the youngest clusters are disregarded (i.e. the points far away from the centre of distribution). It should be noted that for the most frequently used method, when the ages of globular clusters are estimated from the luminosity of the stars of the turnoff points, the difference in ages between the clusters with different content of heavy elements depends strongly on the adopted relation between luminosity and metallicity of the horizontal branch stars. This is why, the younger age of comparatively metal-rich globular clusters is still in question. The method of estimation of ages from the effective temperatures of the turnoff point stars, which assumes coincidence of the positions of the red giant branches of clusters of the same metallicity, can be used to obtain only the difference in ages between the clusters being compared. In this case, the relative ages prove to be more reliable. The ages of the clusters of anomalously small mass distinguished in Fig. 1b have been obtained by both methods in about ten and a half papers. And both the methods suggest their age to be smaller as compared with the majority of metal-poor clusters.

The presence of an explicit correlation between mass and age in the outer halo suggests a time shift of the initial mass spectrum of globular clusters towards small masses. Indeed, massive enough clusters in Fig. 2b appear only at $`t12`$ billion years. However, attention should be drawn to the fact that in the range of ages (12–14 billion years) clusters are observed in a wide range of masses and all the clusters of this age with relatively great masses (NGC 362, 1851, 3201, 6981, IC 4499, and Rup 106) lie in the range of galactocentric distances ($`9÷18`$) kpc, whereas small-mass clusters (Pal 3, 4, 5, 13, 14, Arp 2, and Pyxis) are situated at far greater distances ($`18÷102`$) kpc. It means that at the same age the decisive factor that causes the small mass of clusters is their distance from the galactic centre.

So, the decreasing of the average mass of globular clusters with diminishing age is observed only in the outer regions of the Galaxy. The inner clusters do not reveal such an effect. Note that 10 distant clusters of anomalously small mass selected by us proved to be younger and with an average metallicity $`<[Fe/H]>=1.36\pm 0.14`$ possess anomalously red horizontal branches (except for Arp 2) non-typical for globular clusters genetically related with the Galaxy. All of them belong to the subsystem of the ”young halo”, i.e. they were assumingly captured by the Galaxy at different stages of their evolution (see Zinn (1993), Da Costa and Armandroff (1995)).

It looks as if only far from the centre of the Galaxy, conditions are created that lead to formation of globular clusters with anomalously small masses. The farther from the centre of the Galaxy the more frequently these conditions are realized: almost all the (five out of six) clusters with $`R>60`$ kpc have masses several times smaller than the average over the Galaxy. Probably, here the predicted by Agekian (1962) process of violent dissipation of globular clusters turns to be efficient, as a result of which they may lose the greater part of their masses during a short lifetime. It is caused by the fact that at the stage of condensation of a massive diffuse cloud of a protocluster, every star formed inside it tends to its centre due to the cloud attraction. When passing by the centre of inertia of the cluster, part of the stars have time to accelerate so much that they leave it. In this case the percentage of stars ”evaporated” from the cluster depends on the degree of original sphericity of the protocluster. Only the protoclusters, formed at a considerable distance from any massive objects, including those formed at a considerable distance from the Galaxy centre, can get a ”regular” shape. The existence in the outer halo of older clusters quite rich in stars can be explained in this case by the fact that these clusters (located, mainly, much closer to the galactic centre than the younger ones) formed under the condition of close interaction with other massive objects and at the stage of the protocloud could not acquire a regular spherical shape. (The known very distant massive cluster NGC 2419 is not only an older one, but its central density, as has already been noted, more corresponds to a density of low-mass clusters.) In contrast to them, younger clusters were captured by the Galaxy from a relatively poorly populated, more distant region of the intergalactic space, where they lost a considerable part of mass at the initial stage of evolution. Within such an explanation, the deficiency of massive clusters at large galactocentric distances becomes clear. This explanation is somewhat hindered by the large range of metallicities of our clusters of anomalously small mass ($`1.8<[Fe/H]<0.6`$). However, two clusters with $`[Fe/H]>1.0`$ belong to Saggitarius and one more, the youngest one having at the same time the smallest mass, Pal 1, can be referred with the same degree of probability both to globular clusters and open clusters (Van den Berg, 2000). The remaining six small-mass clusters have metallicities close to $`[Fe/H]1.6`$ with $`\sigma _{[Fe/H]}0.1`$.

Certainly, each of the correlations discussed here turns out to be of little evidence because of the small number of the observed clusters and large errors in determination of their parameters. Nevertheless, since the distances, masses, ages, radial velocities and morphological composition of the horizontal branches of clusters are determined independently, the assembly of all the results suggests that in the outer halo an isolated group of mass globular clusters really exists. To draw a reliable conclusion if the deficiency of small-mass young globular clusters at the Galaxy periphery is related to the initial mass function variation of these objects with the age, or to the existence of violent dissipation in their evolution, it is necessary to measure proper motions of distant clusters. The elements of galactic orbits computed on their basis make it possible to define the supposed sites of their formation and estimate reliably the destruction rate of every cluster.

Acknowledgements: The work was done under the support of the RFBR through grants 00–02–17689 and 02–02–06911.