text:
# Compressibility Consideration in the Boundary of a Strongly Collapsing Bubble

## I Introduction

The problem of non-linear radial oscillations of a gas bubble in a liquid, when it experiences a high amplitude spherical sound field, is an old challenging problem. Several complications are present in this problem arising from the effects of heat conduction, mass diffusion, compressibility, chemical reactions, surface tension, and viscosity. Many authors have studied different aspects of this matter, but a rather complete description has not been presented yet.

The radial dynamics of a bubble in an incompressible liquid is described by the well-known incompressible Rayleigh-Plesset equation Rayleigh:1917 ; Noltingk:1951 . The extension of this equation to the bubble motion in a compressible liquid has been studied by many authors; for example Herring Herring:1941 , Trilling Trilling:1952 , Gilmore Gilmore:1952 , Keller and Kolodner Keller:1956 , Flynn Flynn:1975 , Lastman and Wentzell Lastman:1979 , Löfstedt et al. fstedt:1993 , and Nigmatulin et al. Nigmatulin:2000 . In these works, several different forms of the compressible bubble dynamics equations have been presented, using different approximation for consideration of the compressibility effects. On the other hand, the influences of heat conduction and mass diffusion in the bubble motion have been reported by Hickling Hickling:1963 , Fujikawa and Akumatsu Fujikawa:1980 , and Yasui Yasui:1997 .

In a generalized approach, Prosperetti and Lezzi Prosperetti:1987  used a singular-perturbation method of the bubble-wall Mach number and derived a one-parameter family of equations describing the bubble motion in the first order approximation of compressibility. This family of equations are written as:

$`(1(\eta +1){\displaystyle \frac{\dot{R}}{C}})R\ddot{R}+{\displaystyle \frac{3}{2}}(1{\displaystyle \frac{1}{3}}(3\eta +1){\displaystyle \frac{\dot{R}}{C}})\dot{R}^2=`$
$`{\displaystyle \frac{R}{\rho C}}{\displaystyle \frac{d}{dt}}\left(P_lP_a\right)+(1+(1\eta ){\displaystyle \frac{\dot{R}}{C}})\left({\displaystyle \frac{P_lP_aP_0}{\rho }}\right),`$
where, $`R`$, $`C`$, $`P_0`$, $`P_a`$, and $`\rho `$ are bubble radius, liquid sound speed, ambient pressure, driving pressure, and density of the liquid, respectively. Also, $`\eta `$ is an arbitrary parameter. Equation (I) must be supplemented by a boundary condition equation at the bubble interface to relate the liquid pressure, $`P_l`$, to the gas pressure inside the bubble. Like all previous authors, Prosperetti and Lezzi Prosperetti:1987  used the following incompressible equation for this purpose:

$$P_l=P_g4\mu \frac{\dot{R}}{R}\frac{2\sigma }{R},$$
(2)
where, $`P_g`$, $`\mu `$, and $`\sigma `$ are gas pressure at the bubble interface, liquid first viscosity coefficient, and surface tension, respectively. As Prosperetti and Lezzi showed Prosperetti:1987 , different existing forms of the bubble dynamics equations belong to this single parameter family of equations, corresponding to different values of $`\eta `$. Especially, $`\eta =1`$ and $`\eta =0`$ correspond to Keller and Miksis equation Keller:1980  and Herring-Trilling equation Herring:1941 ; Trilling:1952 , respectively.

Two specific approximations have been considered in the derivation of Eq’ns. (I) and (2). First, Eq’n. (I) has been derived from the Euler equation, in which the liquid viscosity effects have been neglected. However, these effects are not important in the usual applications of the bubble dynamics equations. It can be shown Moshaii:2002  that these effects are significant when $`\mu \rho CR_0`$, where $`R_0`$ is the bubble ambient radius. Therefore, for micron size and larger bubbles the elimination of these effects is thoroughly justified.

The second approximation, which is more important, is the incompressibility assumption of the liquid and the gas at the bubble interface, which has been used in the derivation of Eq’n. (2). Note that, all of the effects of the liquid compressibility in the work of Prosperetti and Lezzi, as well as in all previous works, have been resulted from the liquid motion around the bubble, but not from the bubble boundary condition equation. In fact, all previous authors, on one hand take into account the compressibility of the liquid motion around the bubble, but on the other hand neglect its consideration at the bubble interface. Although Eq’n. (2) has widely been used in all old and recent works Rayleigh:1917 ; Noltingk:1951 ; Herring:1941 ; Trilling:1952 ; Gilmore:1952 ; Keller:1956 ; Flynn:1975 ; Lastman:1979 ; fstedt:1993 ; Nigmatulin:2000 ; Hickling:1963 ; Fujikawa:1980 ; Yasui:1997 ; Prosperetti:1987 ; Keller:1980 , but its applicability for the bubble motion needs to be clarified (especially at the end of the collapse, where the bubble motion is significantly compressible).

In this paper, we have modified Eq’n. (2) considering the effects of compressibility at the bubble interface for both the liquid and the gas. The modified equation has new terms resulted from the effects of two coefficients of viscosity of the liquid and the gas. This work shows that the liquid compressibility effects in the bubble dynamics originate from two sources. First, they originate from the liquid motion around the bubble, which has already been presented in the generalized Eq’n. (I). Second, they are initiated through the liquid compressibility consideration at the bubble interface as presented in this work. We have performed numerical analysis to investigate the influence of the modification of Eq’n. (2).

## II Compressible Bubble Boundary Equation

To derive the compressible bubble boundary equation, we assume that, the motions of the bubble interface and the surrounding liquid are always spherically symmetric. The continuity equation and the radial component of the stress tensor can be written as:

$$\frac{1}{\rho }\left[\frac{\rho }{t}+u\frac{\rho }{r}\right]=\frac{u}{r}\frac{2u}{r}=\mathrm{\Delta },$$
(3)
$$T_{rr}=p+\lambda \mathrm{\Delta }+2\mu \left(\frac{u}{r}\right).$$
(4)
where, $`\rho `$, $`u`$, $`p`$, and $`\mathrm{\Delta }`$ are density, velocity, pressure, and divergence of the velocity, respectively. Also, $`\lambda `$ is the second coefficient of viscosity. Inserting $`u/r`$ from Eq’n. (3), into Eq’n. (4) yields:

$`T_{rr}`$ $`=`$ $`p+(\lambda +2\mu )\left({\displaystyle \frac{u}{r}}+{\displaystyle \frac{2u}{r}}\right)4{\displaystyle \frac{\mu u}{r}}`$ (5)
$`=`$ $`p+(\lambda +2\mu )\mathrm{}4{\displaystyle \frac{\mu u}{r}}.`$
The velocity divergence, $`\mathrm{}`$, can be written as:

$$\mathrm{}=\frac{1}{\rho }\frac{d\rho }{dt}=\frac{1}{\rho c^2}\frac{dp}{dt},$$
(6)
where, the sound speed, $`c`$, is defined as $`c^2=dp/d\rho `$. The boundary continuity requirement at the bubble interface is:

$`T_{rr}(liquid)_R=T_{rr}\left(gas\right)_R+2{\displaystyle \frac{\sigma }{R}}.`$ (7)
Applying Eq’n. (5) for the gas and the liquid parts of Eq’n. (7) leads to:

$`P_l+4{\displaystyle \frac{\mu \dot{R}}{R}}`$ $``$ $`\left(\lambda +2\mu \right)\mathrm{}_l=P_g+4{\displaystyle \frac{\mu _g\dot{R}}{R}}`$ (8)
$``$ $`\left(\lambda _g+2\mu _g\right)\mathrm{}_g2{\displaystyle \frac{\sigma }{R}},`$
where, $`\mu _g`$ and $`\lambda _g`$ are the first and the second coefficients of viscosity of the gas at the bubble interface, respectively. Also, $`\mathrm{\Delta }_l`$ and $`\mathrm{\Delta }_g`$ are the divergence of velocity of the liquid and the gas, respectively. Substituting the divergence of velocity for the liquid and the gas from Eq’n (6) into Eq’n. (8) yields:

$`P_l+4{\displaystyle \frac{\mu \dot{R}}{R}}`$ $`+`$ $`\left({\displaystyle \frac{\lambda +2\mu }{\rho C^2}}{\displaystyle \frac{dP_l}{dt}}\right)=P_g+4{\displaystyle \frac{\mu _g\dot{R}}{R}}`$ (9)
$`+`$ $`\left({\displaystyle \frac{\lambda _g+2\mu _g}{\rho _g}}{\displaystyle \frac{d\rho _g}{dt}}\right)2{\displaystyle \frac{\sigma }{R}},`$
where $`\rho _g`$ is the gas density at the bubble interface. Equation (9) represents the bubble boundary equation containing all effects of the compressibility and the viscosity of both the liquid and the gas. Comparison of Eq’ns. (2) and (9) indicates the existence of three new terms in Eq’n. (9) due to the liquid and the gas compressibility and viscosity effects. Here, we concentrate on the effects of the new term arising from the liquid compressibility. Therefore, the gas viscosity effects are neglected due to their smallness relative to the liquid viscosity as in previous works Rayleigh:1917 ; Noltingk:1951 ; Herring:1941 ; Trilling:1952 ; Gilmore:1952 ; Keller:1956 ; Flynn:1975 ; Lastman:1979 ; fstedt:1993 ; Nigmatulin:2000 ; Hickling:1963 ; Fujikawa:1980 ; Yasui:1997 ; Prosperetti:1987 ; Keller:1980 . Under this circumstance, Eq’n. (9) becomes:

$$P_l+\left(\frac{\lambda +2\mu }{\rho C^2}\frac{dP_l}{dt}\right)=P_g4\frac{\mu \dot{R}}{R}2\frac{\sigma }{R}.$$
(10)
To generalize the argument, we express Eq’ns. (I) and (10) in dimensionless forms. The dimensionless variables are:

$`R^{}={\displaystyle \frac{R}{R_0}},`$ $`\dot{R^{}}={\displaystyle \frac{\dot{R}}{C}},t^{}={\displaystyle \frac{tC}{R_0}},P_l^{}={\displaystyle \frac{P_l}{\rho C^2}},`$ (11)
$`P_g^{}={\displaystyle \frac{P_g}{\rho C^2}},P_a^{}={\displaystyle \frac{P_a}{\rho C^2}},`$
where, $`R_0`$ is the ambient radius of the bubble. Substituting the dimensionless variables into Eq’ns. (I) and (10) results the dimensionless equations:

$`\left(1(\eta +1)\dot{R}^{}\right)`$ $`R^{}\ddot{R^{}}+{\displaystyle \frac{3}{2}}\left(1{\displaystyle \frac{1}{3}}(3\eta +1)\dot{R}^{}\right)\dot{R}^2`$ (12)
$`=`$ $`(1+(1\eta )\dot{R}^{})\left(P_l^{}P_a^{}P_0^{}\right)`$
$`+`$ $`R^{}{\displaystyle \frac{d}{dt^{}}}\left(P_l^{}P_a^{}\right),`$
$$P_l^{}=P_g^{}4\frac{\mu ^{}\dot{R}^{}}{R^{}}2\frac{\sigma ^{}}{R^{}}\left(\lambda ^{}+2\mu ^{}\right)\frac{dP_l^{}}{dt^{}}.$$
(13)
The quantities $`\sigma ^{}`$, $`\lambda ^{}`$, and $`\mu ^{}`$ are the dimensionless form of surface tension, first, and second viscosity coefficients of the liquid and are defined as: $`\sigma ^{}=\sigma /\rho R_0C^2`$, $`\lambda ^{}=\lambda /\rho R_0C`$, and $`\mu ^{}=\mu /\rho R_0C`$, respectively. These dimensionless quantities are basically the inverse of Weber Number and Reynolds Number.

It should be mentioned that, although the effects of compressibility in Eq’n. (12) are in the first order approximation, but these effects have been introduced completely in Eq’n. (13). To quantify the effects of our modification, we have performed our calculations within the framework of Eq’n. (12), using $`\eta =1`$. Of course, similar results are obtained using the other values of $`\eta `$ in Eq’n. (12), because of the similarity of the order of approximation in all equations corresponding to different values of $`\eta `$.

## III Numerical Analysis

To quantify the effects of the new term in Eq’n. (13), numerical analysis was carried out for the conditions of Single Bubble Sonoluminescence (SBSL) Brenner:2002 ; Barber:1997 . The driving pressure in its dimensionless form was $`P_a^{}(t)=P_a^{}\mathrm{sin}\left(\omega ^{}t^{}\right)`$, where $`\omega ^{}`$ is dimensionless angular frequency and is defined as $`\omega R_0/C`$. To have a well-posed problem, the value of the gas pressure at the bubble interface, $`P_g^{}`$, must be specified. It can be determined, in the most complete approach, from the numerical solution of the conservation equations for the bubble interior along with Eq’ns. (12) and (13) simultaneously. Also, the bubble content changes during the bubble evolution because of the chemical reactions at the end of the collapse. In addition, the mass exchange and the heat transfer between the bubble and the surrounding liquid affect the bubble interior.

During the past ten years, several different approaches have been presented to describe the real state of the gas and its evolution considering the above-mentioned complexities. These approaches, which at first assumed inviscid and without consideration of chemical reactions, heat transfer, and mass exchange C.C.WU ; Moss:1994 , became gradually more complex by including dissipating effects of radiative transfer Kondict:1995 , heat transfer, and viscous gas dynamics Voung:1996 ; Yuan:1998 . Recent gas dynamic model of Storey and Szeri Storey:2000  accounts for all effects of chemical reactions and water vapor evaporation and condensation.

During the bubble motion, strong spatial inhomogeneities inside the bubble are not remarkably revealed unless at the end of the collapse. Therefore, the uniformity assumption for the bubble interior seems to be useful and provides many features of the bubble motion Barber:1997 ; Brenner:2002 . By this assumption, the gas pressure inside the bubble is obtained assuming polytropic evolution for the gas bubble as:

$$P_g^{}=\left(P_0^{}+2\sigma ^{}\right)\left(\frac{1a^3}{R^3a^3}\right)^\gamma ,$$
(14)
where, Van der Waals equation of state has been used. Here, $`a^{}`$ is the dimensionless Van der Waals hard core radius; $`a^{}=a/R_0`$. The bubble evolution is assumed to be isothermal ($`\gamma =1.0`$) for the radii larger than $`R_0`$, when the bubble moves relatively slow. While, adiabatic assumption ($`\gamma =\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$ is the ratio of the specific heats of the bubble interior) is applied for the smaller radii, when the bubble experiences rapid changes Barber:1997 .

Hilgenfeldt et al. Hilgenfeldt:1999  extended this model considering $`\gamma `$ to be a function of the Peclet number of the bubble. Of course, the Peclet number changes with the variations of the bubble radius and the bubble velocity. They showed that, many features of SBSL can be explained by their model. However, in practice there is very small difference (less than 5%) between the Hilgenfeldt et al.’s model and the isothermal-adiabatic model Putterman:2001 . Although, these two models can not illustrate production of discontinuities and shock waves at the end of the collapse, but they are really useful for our purpose to investigate the importance of the new term during the collapse. It should be mentioned that, these models have been used in several recent works in sonoluminescence studies ISO-ADI-model .

We have used isothermal-adiabatic model in this paper. Also, we have argued about the extension of the results for more complicated models. Under these circumstances, the time variations of the bubble properties under the framework of Eq’n (12) (for $`\eta =1`$) have been numerically calculated for two cases: (a) compressible boundary condition (Eq’n. 13) and (b) incompressible boundary condition (Eq’n. 2). We used Runge-Kutta method for our numerical analysis. The constants and the parameters used in the calculations were set for an air bubble in water at room temperature, $`T_0=293.0K`$, and atmosphere pressure, $`P_0=1.0atm`$, CRC:1991 ; $`\rho =998.0kg/m^3`$, $`C=1483.0m/s`$, $`\mu =1.01\times 10^3kg/ms`$, $`a^{}=1.0/8.745`$, $`\sigma =0.0728kgs^2`$, and $`\mathrm{\Gamma }=1.4`$. The second coefficient of viscosity of water was set to be $`\lambda =2.23\mu `$. It was derived from the value of the bulk viscosity of water at room temperature that, $`\mu _b=\lambda +(2\mu /3)=2.9\mu `$ Pierce:1991 . Also, the angular frequency of the deriving pressure was $`\omega =2\pi \times 26.5kHz`$.

## IV Results

The results of our calculations have been illustrated in Figs. (1-5) for $`P_a=1.6atm`$ and $`R_0=2.0\mu m`$. Similar values for these parameters were reported in the results of recent experimental works of Ketterling and Apfel Ketterling:2000 , Simon et al. Simon:2001 , and Vazquez et al. Vazquez:2002 . Figure (1) shows the variations of the bubble radius and the gas temperature in one period of the applied pressure field using the compressible boundary condition. The bubble motion has incompressible characteristics during a period, except for an infinitesimal time interval at the end of the collapse. Basically, the new term in Eq’n. (13) originates from the compressibility of the liquid motion. Therefore, it is expected that, its effects not to be revealed, until the bubble motion becomes significantly compressible. The results of our calculations clearly confirm this point. Except for the end of the collapse, the differences between the bubble characteristics resulting from the compressible and incompressible boundary conditions are less than $`1.0\%`$. This result thoroughly justifies the elimination of the new viscous term in Eq’n. (13) for all times, but not for the end of the collapse. Indeed, in Fig. (1), the results for the bubble properties of the two boundary conditions completely coincide. However, the maximum bubble temperature in Fig. 1(b) has a considerable increase relative to that of the incompressible case. This discrepancy instigated us for further concentration of the bubble properties around the minimum bubble radius.

In Figs. (2-4), we have presented the evolution of the bubble characteristics using compressible boundary condition (a), and incompressible boundary condition (b), around the minimum radius time. The time interval of these figures is 40 ps. Since, the bubble experiences its maximum compression at the end of the collapse, the effects of the new term in Eq’n. (13) are more significant during this time interval.

Figure (2) shows the bubble radius evolution for the two boundary conditions. This figure illustrates a considerable difference between the compressible and incompressible cases. After the bubble reaches its minimum radius, a number of small bouncing oscillations appear in the graph of the compressible case, which do not occur in that of the incompressible case. The period of these oscillations is nearly 3 ps. The details of our calculations show that, the times of minimum radius for the two cases are the same. While, the minimum radius for the new Eq’n. (13) is 1% less than that of the old Eq’n. (2).

Figure (3) represents the time variations of the gas temperature near the minimum radius time. It evidently represents that, introducing the new viscous term in Eq’n. (13) strongly affects the gas temperature evolution at the end of the collapse. After the minimum radius time, the incompressible case illustrates smooth behavior for the bubble temperature. While, remarkable sharp peaks appear in the temperature evolution using the compressible boundary condition. Also, as the bubble goes away from the minimum radius, the peaks become weaker. Note that, the value of the maximum temperature increases considerably (nearly 2.4 times) for the new boundary equation.

Comparison of the gas pressure evolution for the two cases is shown in Fig. (4). Similar behaviors as in Fig. (3) are observed in this figure for the difference between the new and the old boundary equations. Note that, the gas pressure is much more sensitive to the presence of the new viscous term compared to the gas temperature. In fact, the presence of the new term increases the maximum pressure more than one order of magnitude.

In Fig. (5), the time variations of the three different pressure terms in Eq’n. (13) have been illustrated near the minimum radius. These terms are due to the effects of surface tension and viscosity, i.e. $`P_\sigma =2\sigma ^{}/R^{}`$, $`P_\mu =4\mu ^{}\dot{R}^{}/R^{}`$, and $`P_{\lambda \mu }=(\lambda ^{}+2\mu ^{})(dP_l^{}/dt^{})`$. This figure clearly shows that, at the end of collapse, the collective effects of the viscous terms are by far greater than the surface tension term. Moreover, the new viscous term, $`P_{\lambda \mu }`$, is the dominant term in this time interval. Note that, the order of the maximum values of the three pressure terms are completely different. $`P_{\lambda \mu }`$ is the highest (up to $`10^3`$) and $`P_\sigma `$ is the lowest (less than $`10^2`$). These results emphasize that, the elimination of $`P_{\lambda \mu }`$ is not reasonable, when the bubble evolves near the minimum radius.

Figures (6) and (7) illustrate the dependence of the peak pressure and the peak temperature of the bubble to the driving pressure and to the ambient radius, for the two boundary equations. Figure (6) is for the case in which, the ambient radius is constant ($`R_0=1.5\mu m`$), while the amplitude of the driving pressure is increased from $`P_a=1.4atm`$ to $`1.8atm`$. On the other hand, Fig. (7) shows the results for the constant amplitude of the driving pressure, ($`P_a=1.7atm`$), with varying $`R_0`$ from $`1.0\mu m`$ to $`5.0\mu m`$. Different values of $`R_0`$ can be experimentally adapted for a specific value of $`P_a`$, adjusting the concentration of the dissolved gas in the liquid Barber:1997 .

The dependence of the maximum pressure and the maximum temperature of the bubble to $`P_a`$ is represented in Fig. (6). There are significant differences between these bubble properties resulted from the new boundary condition with those of the old one. The differences are more considerable for the higher driving pressures. Moreover, the differences are more remarkable for the maximum pressure than the maximum temperature. Note that, for $`P_a=1.8atm`$, the maximum temperature using the new equations is about 8 times greater than that of the old ones. While, the increase of the maximum pressure in this case is about 1400 times. These results indicate that, the compressibility effects should be much higher than what is considered in the incompressible case.

The effects of the variation of $`R_0`$ on the mentioned bubble properties are illustrated in Fig. (7). Similar considerable differences between the new and the old cases as in Fig. (6) are also present in this figure. The differences are more significant for the smaller ambient radii.

Time grid resolution study of the problem shows that, the values of the bubble properties at the end of collapse are sensitive to the step size. But, this sensitivity disappears as the resolution is reduced sufficiently. Figure (8), presents the step size dependence of the value of maximum gas pressure for the compressible case, where $`P_a=1.8atm`$ and $`R_0=1.5\mu m`$. It shows a diminishing step size dependence bellow $`10^{17}s`$.

It should be mentioned that, although the results of the last figures were acquired by assumption that $`\mathrm{\Gamma }=1.4`$ (the bubble content was assumed to be a diatomic gas), however the appearance of the new bubble behavior is independent of the selected value of $`\mathrm{\Gamma }`$ (and gas content). Figure (9) shows the evolution of the bubble temperature for two different states; $`\mathrm{\Gamma }=1.67`$ (monoatomic gas) and $`\mathrm{\Gamma }=1.33`$ (polyatomic gas). The comparison of Figs. (3) and (9) shows that, the new bubble behavior will be established even when the gas content changes. Of course, as Figs. (3) and (9) indicates, the configuration and the number of the new peaks depend on the value of $`\mathrm{\Gamma }`$. The number of peaks is more for the smaller $`\mathrm{\Gamma }`$. While, the values of the peak temperatures are greater for higher $`\mathrm{\Gamma }`$, due to the decrease of the number of the degrees of freedom.

## V Discussion

Although, the isothermal-adiabatic model used in this study does not account for the effects of gas dynamics, chemical reactions and water vapor exchange, but the consideration of these effects can not cover the importance of the new term at the end of the collapse. In this section, the influences of these corrections on the new bubble behavior are argued:

(i) Gas dynamics: Considerations of the gas dynamics effects inside the bubble have been presented by several different approaches C.C.WU ; Moss:1994 ; Kondict:1995 ; Voung:1996 ; Yuan:1998 . These approaches show that, the overall bubble temperature and pressure at the collapse in gas dynamics models are at least no less than those of the isothermal-adiabatic model. In fact, in the gas dynamics models the bubble experiences more compression. Therefore, the pressure upon the liquid will be greater in the gas dynamic models. This means that, the effects of the new viscous term should even be more remarkable in the gas dynamics models.

(ii) Chemical reactions: The temperature of the bubble at the end of the collapse is high enough (higher than 10000 K) to destroy the chemical bonds of N<sub>2</sub> and O<sub>2</sub> molecules of the air bubble. The chemical reactions between the water vapor and the existing oxygen and nitrogen atoms mainly produce very soluble substances in water (HN, NH3, and HNO3), which are completely absorbed. Therefore, a sonoluminescence air bubble mainly contains inert gases. This idea, which was initially presented by Lohse et al. Lohse:1997  is known as dissociation hypothesis (DH) and has been confirmed very well by experimental reports exp . The simulation of chemical reactions at the collapse in a gas dynamics model has been recently presented by Storey and Szeri Storey:2000 . Their work shows that, a considerable decrease appears in the bubble peak temperature due to the consideration of chemical reactions. Since most of the reactions at the collapse are endothermic, their influences can be considered as addition of some extra degrees of freedom Brenner:2002 . This means that, the bubble actually evolves near the minimum radius by an effective exponent, $`\mathrm{\Gamma }_{eff}`$, which is less than the monoatomic exponent ($`\mathrm{\Gamma }=5/3`$). However, according to our results, the appearance of the new behaviors is independent of the values of $`\mathrm{\Gamma }`$. Therefore, the new bubble behaviors should be established even when chemical reactions are introduced in the model.

(ii) Water vapor: Evaporation and condensation of the water vapor between the bubble and the liquid occur during the expansion and compression of the bubble. Recent simulations Yasui:1997 ; Storey:2000 ; Togel:2003  show that, a large amount of water vapor evaporates into the bubble during the expansion. Indeed, at maximum radius, about 90% of the bubble content is water vapor. During the collapse, the water vapor molecules condense to the liquid so that, near the minimum radius only a small fraction of a sonoluminescence bubble is water vapor and the remaining is mainly inert gases Yasui:1997 ; Storey:2000 ; Togel:2003 . Presence of these water vapor molecules inside the bubble decreases the maximum temperature due to the increase of the bubble’s total degrees of freedom. This effect causes more decrease of the effective exponent Brenner:2002 . Therefore, similar to our argument in the last paragraph, the effects of the water vapor can not also cover the strong effects of the compressibility consideration of this work.

## VI Conclusions

The modification of the bubble dynamics equation to account for the viscosity of a compressible liquid was performed by deriving a new equation for the bubble boundary. This equation includes a new term, which has been resulted from simultaneous effects of viscosity and compressibility of the liquid. The new term is the prominent term at the end of the collapse, where the bubble is highly compressed. It exhibits its role by intensifying the strength of the collapse. The more intense the collapse is, the more significant the role of the new term is. Moreover, the new effects can not completely be covered by the dissipating effects of the water vapor and the chemical reactions.

The results of this work strongly indicate that, the neglect of the new term at the end of the collapse in the previously derived equations is not reasonable. Especially, it is more remarkable for high amplitudes single bubble sonoluminescecne. It is expected that, these new theoretical results can be confirmed in the experiments, if resolution of the bubble motion measurement at the end of an intense collapse is less than 0.1 ns. Of course, a stable high amplitude sonoluminescence bubble can be produced, if the concentration of the dissolved gas in the liquid is sufficiently small Simon:2001 .

## ACKNOWLEDGEMENTS

This work was supported by Sharif University of Technology and Bonab Research Center. Partial support of this work by Institute for Studies in Theoretical Physics and Mathematics is appreciated. The authors would like to thank Prof. Andrea Prosperetti for his helpful comments.


format: text

text:
# 1 Introduction

## 1 Introduction

The interest of semi-inclusive hadron production in Deep Inelastic Scattering (DIS) on nuclear targets is at least twofold. First, and on rather general grounds, this reaction should help us to understand how a given QCD medium affects the non-perturbative hadronization mechanism. In particular, the use of nuclei — with a size $`R`$ of a few fm — provides a direct information on the soft time scales $`\mathrm{\Lambda }1/R`$ involved in the dynamics of hadronization. In the second place, such a process also allows one to gain a quantitative insight into the scattering properties of cold nuclear matter. These may serve as a baseline to which properties of hot QCD media should be contrasted. This is clearly an important issue in high energy heavy ion physics where the quark gluon plasma is expected to be formed.

Measurements of semi-inclusive hadron electroproduction on nuclei were performed a decade ago by the European Muon Collaboration (EMC) experiment who reported a significant attenuation of charged hadron spectra in Copper as compared to Deuterium targets . This was a strong indication that quark fragmentation functions get modified in a nuclear environment111The attenuation of hadron electroproduction on a Carbon target was actually first observed at SLAC . However, the SLAC measurements were not normalized to the number of inclusive events $`N_A^e`$ (see Eq. (1)) but the atomic mass number, hence rather reflected the nuclear modifications of parton densities (see discussion in Section 2). For this reason, these data are not considered in the present study.. More recently, the HERMES collaboration at DESY reported on high statistics measurements of charged hadron ($`h^\pm `$), pion ($`\pi ^\pm `$), kaon ($`K^\pm `$), proton and antiproton ($`p`$, $`\overline{p}`$) electroproduction on several nuclei (Deuterium, Nitrogen, Neon, Krypton) and on a wide kinematic acceptance . These two data sets make quantitative analyses possible.

Several suggestions have been advanced so far to account for the experimental results. Models based on the nuclear absorption of the produced hadron , gluon bremsstrahlung , or partial deconfinement in nuclei , were shown to describe the trend of the data fairly well. Following Wang and Wang , we shall assume here that the observed hadron quenching actually comes from the energy loss incurred by hard quarks traveling through the nuclear medium. Although the effects of quark energy loss in DIS on nuclear targets have already been investigated in Refs. , a more complete analysis of EMC and HERMES data is still lacking. This is the aim of the present study.

The outline of the paper is as follows. Hadron electroproduction in nuclei is computed to leading order in the coupling constant (Section 2) using nuclear fragmentation functions modeled in Section 3. Predictions on hadron attenuation are then systematically compared to EMC and HERMES preliminary data as a function of $`\nu `$ and $`z`$ (Section 4). Calculations in Xenon targets in the HERMES kinematic acceptance will also be presented. Finally, we shall discuss and summarize our results in the last Section.

## 2 Hadron production in DIS on nuclei

Hadron production in semi-inclusive DIS can be computed in perturbation theory within the QCD improved parton model. To leading order (LO) in the strong coupling constant $`\alpha _s`$, the multiplicity $`N_A^h`$ of a specific hadron species $`h`$ in lepton-nucleus reactions is expressed in terms of nuclear parton densities $`f_q^A`$ ($`f_{\overline{q}}^A`$) and the $`qh`$ ($`\overline{q}h`$) fragmentation functions, $`D_q^h`$ ($`D_{\overline{q}}^h`$). Normalizing it to the number of DIS inclusive events $`N_A^e`$ , it is given by

$`{\displaystyle \frac{1}{N_A^e}}{\displaystyle \frac{dN_A^h(\nu ,z)}{d\nu dz}}`$ $`=`$ $`{\displaystyle 𝑑x\underset{q,\overline{q}}{}e_q^2\left[Zf_q^{p/A}(x,Q^2)+(AZ)f_q^{n/A}(x,Q^2)\right]\sigma ^{\gamma ^{}q}(x,\nu )D_q^h(z,Q^2,A)}`$ (1)
$`/`$ $`{\displaystyle 𝑑x\underset{q,\overline{q}}{}e_q^2\left[Zf_q^{p/A}(x,Q^2)+(AZ)f_q^{n/A}(x,Q^2)\right]\sigma ^{\gamma ^{}q}(x,\nu )}`$
where $`\nu `$ and $`Q`$ are respectively the energy and the virtuality of the hard photon, and $`z=E_h/\nu `$ is the energy fraction of the virtual photon carried away by the leading hadron. The nucleus is characterized by its atomic mass number $`A`$ and its number of protons $`Z`$. The LO $`\gamma ^{}q`$ cross section appearing in (1) is given by

$$\sigma ^{\gamma ^{}q}(x,\nu )=\frac{4\pi \alpha _{em}Mx}{Q^4}\times \left[1+\left(1\frac{Q^2}{xs}\right)^2\right],$$
(2)
where $`\alpha _{em}`$ is the fine structure constant, $`M`$ the nucleon mass, and $`s`$ the center-of-mass energy squared of the lepton-nucleon reaction. The integral over Bjorken $`x=Q^2/2M\nu `$ in Eq. (1) is given by the $`Q^2`$ acceptance of the experiment.

At small $`x0.05`$ and low $`Q^2`$, the parton distribution functions (PDF) in nuclei $`f_q^A`$ differ significantly to those measured in a free proton, $`f_q^p`$ . These so-called shadowing corrections may therefore affect the hadron production cross section (1). In the following, however, we shall not consider such corrections and assume $`f_q^{p/A}=f_q^p`$ for two reasons. First, the data we would like to discuss cover a Bjorken $`x`$ range, $`x0.1`$, where the effects of nuclear shadowing prove rather small. Furthermore, the normalized hadron multiplicity Eq. (1) is rather insensitive to the absolute magnitude of the parton densities. As noted in , the shadowing corrections mostly cancel out in the ratio (1) and finally affect hadron production by a percent at most.

The parton densities in a proton $`f_q^p`$ we use in the present analysis are given by the MRST 2001 LO parameterization . It has been checked however that similar results were obtained using either the GRV LO  and CTEQ5L  parton distribution functions. The neutron PDF $`f_q^n`$ are directly deduced from $`f_q^p`$ by isospin conjugation, i.e. $`u^p=d^n`$, $`d^p=u^n`$, $`\overline{u}^p=\overline{d}^n`$, $`\overline{d}^p=\overline{u}^n`$, and $`\overline{s}^p=\overline{s}^n`$.

In this paper, we would like to discuss the nuclear dependence of semi-inclusive hadron production. In particular, the theoretical calculations will be compared to the hadron production ratio in a heavy nucleus $`A`$ with respect to a (light) Deuterium target,

$$R_A^h(z,\nu )=\frac{1}{N_A^e}\frac{dN_A^h(z,\nu )}{d\nu dz}/\frac{1}{N_D^e}\frac{dN_D^h(z,\nu )}{d\nu dz},$$
(3)
which has been reported experimentally. When $`x`$ is not too small, hadron production (1) is dominated by the fragmentation of valence quarks. Assuming for simplicity that only the up quark channel contributes to the process, the double production ratio Eq. (3) approximately reduces to

$$R_A^h(z,\nu )D_u^h(z,Q^2,A)/D_u^h(z,Q^2,D).$$
(4)
Therefore, hadron production in DIS on nuclei is mostly sensitive to the nuclear dependence of the fragmentation functions. In particular, the reported depletion of hadron multiplicity on heavy nuclei $`R_A^h(z,\nu )<1`$ on a wide kinematic range is a clear hint that fragmentation functions are somewhat modified in the medium.

The starting point to determine the hadron attenuation $`R_A^h(z,\nu )`$ is therefore the computation of nuclear fragmentation functions $`D(z,Q^2,A)`$. This task is carried out in the following section.

## 3 Nuclear fragmentation functions

### 3.1 Model for medium-modified fragmentation functions

To leading order in the coupling and at leading twist, the virtual photon picks up a quark222In the following discussion, we shall ignore the fragmentation of antiquarks as $`x`$ is not too small. The theoretical computations Eq. (1) consider of course this channel as well. in the nucleus which subsequently hadronizes into the observed leading hadron. In presence of a QCD medium, however, the hard quark suffers multiple scattering and radiates soft gluons all along its path. Due to this medium-induced gluon radiation, the quark energy is reduced from $`E=\nu `$ to $`E=\nu ϵ`$ at the time of hadronization. This picture is schematically depicted in Figure 1.

How to relate these final state interactions to nuclear fragmentation functions remains however unclear. Such an attempt has been performed recently in a series of papers  in which a higher-twist perturbative framework has been applied successfully to describe hadron production in DIS and heavy ion reactions.

Here, we shall rather adopt the effective model suggested in Ref. . Within this approach, the quark energy shift leads to a rescaling of the momentum fraction in presence of a QCD medium,

$$z=\frac{E_h}{\nu }z^{}=\frac{E_h}{\nu ϵ}=\frac{z}{1ϵ/\nu },$$
(5)
where $`E_h`$ is the measured hadron energy and $`ϵ`$ the energy lost by the hard quark while going through the nucleus. Consequently, the fragmentation functions in nuclei $`D_q^h(z,Q^2,A)`$ may simply be expressed as a function of the standard (vacuum) fragmentation functions $`D_q^h(z,Q^2)`$ through

$$zD_q^h(z,Q^2,A)=_0^{\nu E_h}𝑑ϵD(ϵ,\nu )z^{}D_q^h(z^{},Q^2).$$
(6)
Here, $`D(ϵ,\nu )`$ denotes the probability for a quark with energy $`E=\nu `$ to lose an energy $`ϵ`$. This probability distribution, or quenching weight, will be made explicit hereafter. Since fragmentation functions fall steeply at large $`z`$, even the small shift $`\mathrm{\Delta }z=z^{}zzϵ/\nu `$ in (5) may substantially suppress the hadron yields in nuclear targets.

Note that in the Bjorken limit $`\nu ϵ`$, hence $`Qϵ`$ at finite $`x`$, the effects of the final state interactions become negligible, $`zz^{}`$. No more nuclear dependence of the fragmentation functions is observed $`D_q^h(z,Q^2,A)=D_q^h(z,Q^2)`$, in agreement with QCD factorization theorems .

The vacuum fragmentation functions $`D_q^h(z,Q^2)`$ into various hadron species (charged hadrons, pions, kaons) which enter Eq. (6) are here given by the leading order analysis of $`e^+e^{}`$ data performed by Kretzer .

Let us now turn to the computation of the quenching weights for hard quarks going through nuclear matter.

### 3.2 Quenching weight

As apparent in Eq. (6), the knowledge of the probability distribution $`D(ϵ,\nu )`$ is essential to determine the nuclear fragmentation functions. This point has been raised recently by Baier, Dokshitzer, Mueller, and Schiff (BDMS) in Ref.  in which they express $`D(ϵ,\nu )`$ as a Poisson series

$$D(ϵ,\nu )=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n!}\left[\underset{i=1}{\overset{n}{}}𝑑\omega _i\frac{dI(\omega _i,\nu )}{d\omega }\right]\delta \left(ϵ\underset{i=1}{\overset{n}{}}\omega _i\right)\mathrm{exp}\left[_0^+\mathrm{}𝑑\omega \frac{dI(\omega ,\nu )}{d\omega }\right],$$
(7)
assuming gluon emissions to be independent.

The expression Eq. (7) explicitly relates the probability distribution in the energy loss to the gluon spectrum $`dI/d\omega `$ radiated by hard quarks produced in QCD media, which has been computed perturbatively by Baier, Dokshitzer, Mueller, Peigné, and Schiff (BDMPS). To first order in the quark energy $`𝒪\left(1/\nu \right)`$, it reads

$$\frac{dI(\omega ,\nu )}{d\omega }=\frac{\alpha _sC_F}{\pi \omega }\left(1\frac{\omega }{\nu }\right)\mathrm{ln}\left[\mathrm{cosh}^2\sqrt{\frac{\omega _c}{2\omega }}\mathrm{sin}^2\sqrt{\frac{\omega _c}{2\omega }}\right]\mathrm{\Theta }(\nu \omega ),$$
(8)
where $`C_F=4/3`$ is the Casimir operator in the fundamental representation and $`\alpha _s=g^2/4\pi 1/2`$ the strong coupling constant. Note that $`\alpha _s`$ in (8) is evaluated at a soft scale given by the virtuality of the radiated gluons, while the renormalization scale appearing in the parton densities and fragmentations functions in (1) is given by the photon virtuality, $`Q`$. The appearance of these two scales follows our implicit assumption on the factorization between the production process and the final state interactions. The typical gluon energy in the spectrum is set by the relevant scale

$$\omega _c=\frac{1}{2}\widehat{q}L^2$$
(9)
which characterizes the QCD medium. The gluon transport coefficient $`\widehat{q}`$ reflects the medium gluon density (say, its scattering property) while $`L`$ is the length of matter covered by the hard quark. Note that $`\widehat{q}`$ can be very large in practice in a hot pion gas or quark gluon plasma . Energy loss of partons may thus be used to probe such dense media expected to be formed in high energy heavy ion reactions . The absolute values for $`\omega _c`$ in nuclei will be discussed below.

Using Mellin transform techniques, BDMS resumed the Poisson series to eventually give the quenching weight Eq. (7) a simple integral representation. This allows for the numerical computation of $`D(ϵ,\nu )`$ from the BDMPS gluon spectrum (8) performed in . The distribution has been given a simple analytic parameterization which we shall use in the present analysis.

*Transport coefficient for nuclear matter $`\widehat{q}`$*

The quenching of hadron spectra depends crucially on the transport coefficient for nuclear matter. The higher the $`\widehat{q}`$, the larger the typical quark energy loss $`ϵ\widehat{q}`$, hence the stronger the depletion in large nuclei. Based on the (not too) small $`x`$ gluon distribution in a nucleus, BDMPS provide a perturbative estimate for the transport coefficient in nuclear matter, $`\widehat{q}0.25`$ GeV/fm2 .

A slightly larger — although consistent — value has been obtained recently from our LO analysis of Drell-Yan production data  measured in $`\pi ^{}`$-A collisions by the NA3 collaboration,

$$\widehat{q}=0.72\mathrm{GeV}/\mathrm{fm}^2.$$
(10)
This value will be taken for the numerical applications to come. It should already be kept in mind, however, that this value could not be determined with a high accuracy ($`1\sigma =0.54`$ GeV/fm2) because of the large statistical errors in the Drell-Yan measurements.

As we shall discuss later (Section 5), hadron depletion in DIS on light nuclei (say, Nitrogen and Neon targets) may provide an independent estimate for the transport coefficient $`\widehat{q}`$. In particular, the high statistics HERMES data should prove of considerable help to constrain further the estimate (10) from Drell-Yan production.

*Length $`L`$*

Once $`\widehat{q}`$ is fixed, only the mean length $`L`$ covered by the hard quark in the nuclear medium remains to be determined. Assuming that hadronization occurs well after the hard quark has escaped from the nucleus, in other words when the hadron formation time $`t_f`$ is large, $`L`$ may simply be estimated for any given density profile $`\rho (𝐫)`$. Taking for simplicity a sharp sphere density $`\rho (𝐫)=\rho \times \mathrm{\Theta }(R|𝐫|)`$, it is given by

$$L_{\mathrm{max}}=\frac{3}{4}R,$$
(11)
where the nuclear radius $`R=(3/(4\pi \rho ))^{1/3}A^{1/3}1.17A^{1/3}`$ fm with $`\rho =0.15`$ fm-3.

When $`t_f`$ gets smaller than twice the nuclear radius, however, quarks produced on the “back” side of the nucleus hadronize while still inside the medium. The length (11) given above is therefore no longer correct. Integrating over the nuclear volume, we now obtain

$$L(t_f)=t_f\times \left[1\frac{3}{8}\frac{t_f}{R}+\frac{1}{64}\left(\frac{t_f}{R}\right)^3\right]\mathrm{if}t_f2R.$$
(12)
Of course, $`L(t_f)=L_{\mathrm{max}}`$ at $`t_f=2R`$. For small formation times (large nuclei), $`t_f/R1`$, energy loss effects are strongly reduced as the quark has propagated over shorter distances than (11), $`L(t_fR)t_fL_{\mathrm{max}}`$. It is implicitly assumed here that a color singlet hadron ($`C_A=0`$ in (8)) does not lose any energy by gluon radiation. The produced hadron in the medium may nevertheless interact and be absorbed in the nucleus. In this approach — that we want to keep as simple as possible —, this mechanism is not taken into account. However, we shall come back to this point in the discussion section.

The soft scales involved in the hadronization process prevent any perturbative calculation of the formation time $`t_f`$. Therefore, the absolute value for $`t_f`$ as well as its kinematic dependence is not clearly known. We shall take the Bialas-Gyulassy estimate  obtained in the framework of the Lund model333It has been stressed in  that Eq. (13) slightly overestimates the Lund model predictions with standard parameters at small $`z`$. We checked however that this does not strongly affect our results.

$$t_f=\left(\frac{\mathrm{ln}(1/z^2)1+z^2}{1z^2}\right)\times \frac{z\nu }{\sigma },$$
(13)
where $`\sigma `$ is a non-perturbative scale in the order of the string tension, although slightly reduced in a nuclear environment . We take $`\sigma =0.75`$ GeV/fm throughout this study. The formation time (13) exhibits a maximum around $`z0.3`$ and decreases like $`t_f(1z)`$ at large $`z`$, observed as well within gluon bremsstrahlung models .

Plugging the Lund formation time (13) into Eq. (12), the mean (normalized) length $`L/R`$ is plotted as a function of the quark energy $`\nu `$ and the hadron momentum fraction $`z`$ in Figure 2. Also shown as a dashed line is the asymptotic length $`L_{\mathrm{max}}/R`$ when hadronization occurs outside of the nucleus. We observe that $`L`$ gets much smaller than $`L_{\mathrm{max}}`$ when $`z`$ approaches 0 or 1, as $`t_f`$ vanishes. We shall thus expect in these specific kinematic limits a dramatic discrepancy with the prediction $`L=L_{\mathrm{max}}`$ assuming hadronization outside of the nucleus. Because of the Lorentz boost $`L\nu `$ in (13), highly energetic quarks $`\nu \sigma R`$ will fragment into hadrons after escaping the nuclear medium. As discussed in the next section and already apparent in Figure 2, finite formation time effects become irrelevant to describe high energy data.

To give the reader a feeling for such formation time effects, hadron attenuation will be computed in the following account assuming both the asymptotic length $`L_{\mathrm{max}}`$ (11) (dashed lines) as well as the more realistic estimate Eq. (12) (solid lines).

## 4 Results

Hadron production in semi-inclusive DIS on nuclear targets may now be determined from the perturbative expression (1) using the nuclear fragmentation functions derived in the latter section. Calculations are performed using the acceptance of the experiments. We have used the mean kinematic variables $`\nu `$, $`z`$, and $`Q^2`$ measured experimentally whenever they were available, as it was shown in  that these were slightly underestimated (by less than 10%) by the LO prediction (1).

The quenching of hadron multiplicity $`R_A^h(\nu )`$ and $`R_A^h(z)`$ is compared to the EMC  and HERMES  data on Cu as well as on N and Kr nuclei, respectively. We shall also compare in this section our results to the HERMES preliminary data on Neon . Finally, predictions on a Xenon target in the HERMES kinematic window are presented.

### 4.1 $`\nu `$ dependence

We start the section by first considering the quark energy dependence $`\nu `$ of the hadron quenching. As previously emphasized, the energy loss effects should vanish as the energy $`\nu `$ gets large as compared to the typical scale $`\omega _c`$ of the medium. Naturally, the quenching factor $`R_A^h(\nu )`$ should therefore increase with $`\nu `$ and eventually approaches unity at asymptotic energies.

This is what is observed experimentally. In Figure 3 is plotted the charged hadron attenuation $`R^{h\pm }(\nu )`$ in a Cu nucleus on a wide range, $`15\nu 200`$ GeV. The calculations are found to be in good agreement with the measurements performed by the EMC collaboration. Notice that these high energy data do not exhibit any formation time effects in such a small nucleus, $`\nu /\sigma R6`$, as shown from the tiny difference between the dashed and solid lines in Figure 3.

Let us turn to lower energy data. The HERMES collaboration measured semi-inclusive hadron production on Nitrogen and Krypton targets. What is more, the experimental apparatus allows for a particle identification of various hadron species, such as pions, kaons, protons and anti-protons. Since no fragmentation functions are provided in the (anti-)proton channel, only the quenching of charged hadron, pion, and kaon spectra will be presented.

Figure 4 (upper panel) indicates that the quenching of charged hadrons on Nitrogen is extremely well reproduced by our predictions, while the calculations slightly overestimate the attenuation in the pion (isospin averaged) sector. It should be noted that the stronger pion quenching $`R^\pi (\nu )<R^h(\nu )`$ reflects only the different acceptance444More precisely, it comes from the different cuts in the $`z`$ acceptance. Pions were measured with $`z>0.5`$ while a cut $`z>0.2`$ has been applied to the charged hadron channel. in the HERMES experiment. In the Kr target, the suppression of charged hadrons is again well described by the calculations. The trend predicted in the pion channel in Kr is close to what is observed experimentally, although the absolute magnitude of the pion quenching slightly underestimates the data ($`58\%`$).

Perhaps more interesting are the HERMES preliminary measurements on kaon quenching in krypton, $`R_{\mathrm{Kr}}^{K\pm }(\nu )`$ (figure 4, lower panel). A significantly larger depletion is observed in $`K^{}`$ production on the whole energy range, $`R^K^{}(\nu )<R^{K^+}(\nu )`$. Furthermore, this noticeable isospin dependence in the kaon sector is apparent in the theoretical predictions as well. It is indeed clear from Eqs. (4) to (6) that the quenching is directly related to the slope of the vacuum fragmentation functions555We omit, for clarity, the integral $`𝑑ϵD(ϵ)`$ in Eq. (14).

$`R_A^h(\nu )`$ $``$ $`1+{\displaystyle \frac{D_u^h}{z}}{\displaystyle \frac{1}{D_u^h(z)}}{\displaystyle \frac{zϵ}{\nu }}`$ (14)
$``$ $`1\eta _u^h\times {\displaystyle \frac{zϵ}{\nu (1z)}}`$
where $`\eta _u^h`$ is the slope of the fragmentation functions at large $`z`$, $`D_u^h(z,Q^2)(1z)^{\eta _u^h}`$. Schematically, the up quark will preferably fragment into a $`K^+`$ by picking up easily a surrounding $`\overline{s}`$ quark in the sea (valence-type fragmentation channel). Inversely, the difficulty to produce a negative kaon from a hard up quark (sea-type) will translate into a steeper slope, $`\eta _u^K^{}=\eta _u^{K^+}+1`$ and thus in a stronger quenching Eq. (14). At very small $`x1`$, semi-inclusive production originates from the sea quark fragmentation, $`u=\overline{u}`$ 666Furthermore, the sea is almost flavor symmetric $`\overline{u}\overline{d}`$ at small $`x`$. However, down quarks will be disfavored due to their electric charge.. We therefore predict that no more isospin effect should be observed in this region, and thus $`R^K^{}(\nu )R^{K^+}(\nu )`$

Formation time effects start to become visible in Krypton targets when $`\nu `$ is not too large, $`\nu 12`$ GeV (see Figure 4). Nevertheless, the rather small differences between the dashed and solid curves indicate that the krypton nucleus (A=84) is not too large to probe the hadronization in the medium. More pronounced effects are expected in heavier nuclei.

Finally, hadron attenuation has been computed on a Neon target (Figure 5, upper panel) measured very recently by the HERMES collaboration . While the observed pion quenching is well reproduced by the calculations, the predicted attenuation noticeably overestimates the preliminary data in the kaon sector, and particularly in the $`K^{}`$ channel. This feature was already present, although less pronounced, in the comparison with Krypton measurements. The underestimate of the slopes of the flavor-tagged $`uK^+(K^{})`$ fragmentation functions, poorly constrained by $`e^+e^{}`$ data, may be at the origin of this discrepancy . In particular, the assumption made in Ref.  that the sea type channels (e.g. $`D_u^K^{}`$) get suppressed by one extra power in $`(1z)`$ (see discussion above) seems to be disfavored by experimental data on kaon production .

Figure 5 (lower panel) also displays our predictions on Xenon nuclei (A=131), in which formation time effects get larger. As later stressed in the discussion section, measurements in such a heavy nucleus will be of crucial importance to disentangle the possible various mechanisms responsible for the observed hadron attenuation.

### 4.2 $`z`$ dependence

Another way to probe the dynamics of hadronization in a nuclear medium is to look at the $`z`$ dependence of hadron quenching. In our calculations, we expect a stronger hadron quenching as $`z`$ increases for two reasons. First, as can been seen from the rough estimate Eq. (14), the slope of the fragmentation functions increases with $`z`$, and thus so does the attenuation $`R_A^h(z)`$.

The second reason that leads to significant hadron depletion at large $`z`$ actually comes from the restricted phase space. Measuring hadrons with large momentum fraction $`z`$ selects on those (rare) fragmenting quarks that have lost a tiny amount of energy. Assuming for simplicity that the fragmentation functions are flat, $`z^{}D(z^{},Q^2)zD(z,Q^2)`$, the quenching factor may be approximated by

$$R(\nu ,z)_0^{(1z)\nu }𝑑ϵD(ϵ,\nu ),$$
(15)
which vanishes as $`z`$ approaches 1. Let us note however that the ratio (15) will tend to a constant $`p_0>0`$ in the parton energy loss frameworks developed by Gyulassy, Lévai, Vitev  and Wiedemann  in which the quenching weight $`D(ϵ,\nu )`$ contains a discrete part, $`p_0\delta (ϵ)`$ . Similar behavior is also observed within the BDMPS formalism when ultra soft gluons are removed from the perturbative spectrum .

Looking at Figure 6, our predictions on the $`z`$ dependence of the charged hadron quenching $`R^{h^\pm }(z)`$ in Copper targets describe well the EMC measurements, with the exception of the highest $`z`$ bin. Again, high energy quarks hadronize well after escaping the nucleus in this whole $`z`$ range.

On the contrary, we do observe huge formation time effects in the HERMES data which are taken at much lower energy $`\nu 10`$$`20`$ GeV (Figure 7). The predictions fairly reproduce the HERMES measurements for all hadron species when $`z`$ is not too high. At large $`z0.7`$, nevertheless, the discrepancy between the predictions assuming $`L=L_{\mathrm{max}}`$ (dashed lines) and the data is striking. As discussed above, the hadron attenuation $`R(z)`$ falls too steeply as $`z`$ goes to one, unlike the experimental results. This does not come as a surprise as we have already anticipated in the latter section that this assumption is no longer justified since $`Lt_f(1z)\nu /\sigma L_{\mathrm{max}}`$.

When the realistic length Eq. (12) is used, the typical scale $`\omega _cL^2(1z)^2`$ is getting smaller and weakens the quark energy loss — and subsequently the hadron attenuation — as $`z`$ increases. This effect somehow compensates the quenching due to the slope of the fragmentation functions (14) and to phase space restriction (15) at moderate $`z`$, and eventually takes over at $`z0.9`$ (see Figure 7). This non-monotonic behavior of the hadron depletion is actually a manifestation of the quadratic behavior of the parton energy loss in QCD media. We may indeed convince ourselves of this by replacing the energy loss $`ϵ\omega _c(1z)^2`$ in the rough estimate (14). It is however doubtful whether this could seen in the data as the final state interaction of the produced hadron may somehow modify this picture. As shown in Figure 7, a nice agreement is reached between these calculations (solid lines) and the HERMES data for all $`z`$.

The predicted hadron attenuation in Ne and Xe targets is also plotted in Figure 8 as a function of $`z`$. The preliminary HERMES measurements in Neon are well described by these predictions. Once more, let us emphasize that measurements in Xenon targets will prove essential to further constrain the hadron formation time estimates at large $`z`$.

Let us finally stress that theoretical uncertainties might get large in the high $`z1`$ region. In particular, important higher twist corrections may come into play in this restricted region of phase space , given moreover the fact that the typical values $`Q^224`$ GeV2 measured in these reactions are pretty small. Moreover, the resummation of soft and collinear gluons leads to large contributions $`\alpha _s\mathrm{log}(1z)/(1z)`$ in the fragmentation functions  which are not taken into account in the leading order analysis  we use here.

## 5 Discussion and summary

The nuclear dependence of semi-inclusive hadron production in DIS comes in our calculations from the sole effect of the quark energy loss in nuclei. Of course, many other mechanisms, which have not been considered here, could play a role in the hadron quenching observed experimentally. The absorption of the (pre)hadrons produced in the medium may for instance substantially affect its production . This will be particularly true at large $`z`$ and/or small $`\nu `$ when the formation time is small. It was also stressed that the partial deconfinement which could occur in large nuclei rescales the fragmentation functions and eventually suppress the hadron yields .

However, although such effects have not been taken into account here, a very nice description of the EMC and more recent HERMES preliminary data for several hadron species and various nuclei is still achieved. The already good agreement assuming only energy loss effects in our approach is therefore a hint that the above mentioned competing processes play a rather modest role. Let us moderate this statement. The strength of each of these mechanisms is governed by one parameter which is in general poorly constrained. Take the quark energy loss for instance. We have used here the transport coefficient for nuclear matter $`\widehat{q}=0.72`$ GeV/fm2 extracted from Drell-Yan data, whose absolute magnitude is not known with a great accuracy. Similarly, the nuclear absorption strongly depends on the pre-hadron cross section $`\sigma ^{}`$ about which little is known.

It is therefore a delicate matter to disentangle these various processes. We suggest a couple of observables which may help us in this prospect. First, we have seen that the mean length $`Lt_f`$, hence the scale $`\omega _c`$, becomes independent of $`R`$ when the nuclear radius is large, $`R\nu /\sigma `$. The effects of quark energy loss on hadron quenching thus somehow saturates when the nuclei are sufficiently heavy. This can already be seen from the comparison of the rather similar hadron quenching in Kr and Xe targets, $`R_{\mathrm{Kr}}R_{\mathrm{Xe}}`$ (Figures 7 and 8, solid lines). It is clearly at variance with the absorption models which predict of course a much stronger attenuation as the nucleus gets larger. We therefore predict that the ratio of hadron depletion in two large nuclei (say, Xe/Kr) will provide a good measurement of hadron absorption, only. Inversely, hadron quenching in small nuclei could help to constrain the transport coefficient in nuclear matter, without being spoiled by hadronic absorption777Let us however mention that rescaling models predict a non-negligible suppression in such small nuclei .. Fitting $`\widehat{q}`$ to the charged hadron and pion quenching in Nitrogen and Neon data, we found $`\widehat{q}=0.59\pm 0.18`$ GeV/fm2 and $`\widehat{q}=0.70\pm 0.06`$ GeV/fm2, respectively, in excellent agreement with the Drell-Yan estimate888This consistency between Drell-Yan and DIS data was discussed in Ref. , Eq. (10), which suffered from a larger statistical uncertainty.

We also think that the isospin dependence of the kaon quenching could be important. As already mentioned, a similar attenuation ratio $`R^K^{}(\nu ,z)R^{K^+}(\nu ,z)`$ is expected at small $`x`$ if quark energy proves the most relevant mechanism at work. On the contrary, a stronger $`K^{}`$ quenching should still be observed in the framework of absorption models from the larger interaction of negative kaons. The $`Q^2`$ dependence (at a given $`\nu `$) of the ratio $`R^{K^+}/R^K^{}`$ may thus clarify the relevant underlying processes.

To conclude, a leading order calculation of semi-inclusive hadron production in DIS on nuclei has been performed. The fragmentation functions were modified to take into account the energy loss incurred by hard quarks while they propagate through the nuclear medium. The quenching of hadron spectra was then compared to the EMC and HERMES preliminary data in several nuclear targets, as a function $`\nu `$ and $`z`$. Using the transport coefficient determined from Drell-Yan production, a nice agreement between our predictions and the experimental measurements is observed for the various hadron species. In particular, we emphasized that hadron production in small nuclei may help us to determine and constrain what the transport coefficient (hence the quark energy loss) of nuclear matter is. Measurements of semi-inclusive hadron production in heavy nuclei are highly desirable since, in addition to be a sensitive probe of formation time effects, these prove crucial to disentangle quark energy loss on the one hand from hadron absorption on the other hand.

Finally, it should be interesting to constrain, from the HERMES and EMC data, the rapidity dependence of the quenching of charged hadron and pion hadroproduction inclusive spectra in deuteron-gold collisions at RHIC energy .

### Acknowledgments

I would like to thank Néstor Armesto, Valeria Muccifora and Jianwei Qiu for discussions and many useful comments about this study, and Garegin Elbakyan for providing me with the HERMES preliminary Neon data.


text:
# Thermal fluctuations in ultrasmall intrinsic Josephson junctions

## I Introduction

In the strongly anisotropic cuprate superconductors, such as Bi2Sr2CaCu2O8+δ (BSCCO) or Tl2Ba2CaCu2O8+δ (TBCCO), the CuO2 layers together with the intermediate material form a stack of Josephson junctions. In the presence of a bias current perpendicular to the layers each junction of the stack is either in the resistive or in the superconducting state leading to the well-known multibranch structure of the IV-curves, see Ref. Kleiner et al., 1992; Kleiner and Müller, 1994; Müller, 1994; Yurgens, 2000 and references therein.

The effect of thermal noise on the properties of Josephson junctions made from low-temperature superconductors, and in particular on the critical current and voltage-current curves, has been investigated experimentally and theoretically in many publications, see Ref. Fulton and Dunkleberger, 1974; Iansiti et al., 1989; Kautz and Martinis, 1990; Martin et al., 1993. For high-temperature superconductors considerable interest has been attracted to the problems of thermal fluctuations in the vortex state, see Ref. Glazman and Koshelev, 1991; Blatter et al., 1994; Bulaevskii et al., 1996a; Koshelev, 1996; Koshelev and Bulaevskii, 1999, but, to our knowledge, no detailed study of the influence of thermal fluctuations on the critical current and current-voltage curves of mesoscopic intrinsic Josephson junctions has been carried out. Some fluctuation effects are reported recently in Ref. Warburton et al., 2003.

In our paper we present the results of experimental and theoretical investigations of small area intrinsic Josephson junctions, when the Josephson energy $`E_J(T)=J_c(T)/2e`$ is of the order of $`kT`$ ($`k`$ is the Boltzmann constant). For temperatures $`TT^{}`$ thermal fluctuations become important. This cross-over temperature can be estimated with help of Fig. 1, where the Ambegaokar-Baratoff temperature dependence of the critical current density $`j_c(T)`$ is plotted together with the linear functions $`2ekT/\mathrm{}A`$ (marked as $`kT=E_J`$ in the figure) for mesas with different area $`A`$. $`T^{}`$ can be read-off from the crossing points of the curves.

We discuss in the following junctions with cross-over temperatures in the range of 10-60 K ($`A0.5\times 0.5\mu `$m). Note that for all our samples the charging energy $`E_c=e^2/2C`$ is still small compared to $`E_J`$, such that quantum effects (Coulomb blockade and quantum tunneling) are not essential. From the behaviour of the current-voltage characteristics of the junctions three temperature regimes can be distinguished. At low temperatures, $`T<T^{}`$, straight lines with hysteretic jumps to other resistive branches are observed (see Fig. 3 and Fig. 8). At temperature $`TT^{}`$ the current voltage curves are bended near the critical current (like in Fig. 4 and Fig. 11), but still some hysteresis is found. At high temperatures, $`T>T^{}`$, no more jumps in the current voltage curve occur and no real critical current exists, but still a maximum in the curvature can be distinguished (see Fig. 5), which can be regarded as maximum of the supercurrent. In order to describe all temperature regimes together we introduce the concept of a ”switching” current, which will be defined below. It replaces the critical current in the non-hysteretic high-temperature regime. A plot of the switching current as function of temperature shows a rather non-monotonic behaviour, which is very different from the AB-curve found for the critical current of large junctions. We will show below that the temperature dependence of the switching current, in particular the sharp drop, can be explained by phase diffusion processes. This will be done with help of an extended RSJ model containing thermal current noise and the coupling to an external RC circuit (Fig. 2). The latter is necessary to obtain the observed hysteretic behaviour for the phase diffusion branch at intermediate temperatures $`TT^{}`$. In our case the external RC circuit is formed most probably by a large capacitance between the ground superconducting layer and the upper electrode together with the resistance of the contact between mesa and upper electrode.

The paper is organized as follows. In section II we discuss typical voltage-current curves at different temperatures and give a brief summary of the theory used in calculations. In section III sample preparation and experimental set-up are described. The main section IV is devoted to experimental results and their discussion. Here also a discussion of magnetic field effects is included.

## II Outline of the theory

### II.1 RSJ + Johnson noise

The simplest model to describe thermal phase fluctuations in Josephson systems is the usual RSJ model with Johnson noise

$$\frac{\mathrm{}C}{2e}\frac{d^2\phi }{dt^2}+\frac{\mathrm{}}{2eR}\frac{d\phi }{dt}+J_c\mathrm{sin}\phi =J+J_T,$$
(1)
$$J_T(t)J_T(t^{})=\frac{2kT}{R}\delta (tt^{}).$$
(2)
Here $`\phi `$ is the Josephson phase difference, $`C`$, $`R(T)`$, $`J_c(T)`$ are junction capacitance, resistance, and critical current, $`J`$ is the bias current, $`J_T(t)`$ is a random current noise. This equation can be written in dimensionless form as

$$\beta \frac{d^2\phi }{d\tau ^2}+\frac{d\phi }{d\tau }+\mathrm{sin}\phi =j+j_T,$$
(3)
$$j_T(\tau )j_T(\tau ^{})=2\gamma \delta (\tau \tau ^{}),$$
(4)
where $`\tau =\omega _ct`$, $`j=J/J_c`$, and $`\omega _c=2eRJ_c/\mathrm{}`$ is the characteristic frequency.

Here we introduce two main parameters:
(i) the McCumber parameter

$$\beta =\frac{\omega _c^2}{\omega _p^2}=\frac{2eR^2CJ_c}{\mathrm{}},$$
(5)
where $`\omega _p^2=2eJ_c/\mathrm{}C`$ is the Josephson plasma frequency;
(ii) the fluctuation parameter (normalized temperature)

$$\gamma =\frac{kT}{E_J}=\frac{2ekT}{\mathrm{}J_c}.$$
(6)
### II.2 Extended RSJ model

As mentioned above, for the description of our experiments we need a frequency dependent damping which is achieved by the coupling to an external circuit as shown in Fig. 2. This model was already considered in detail by Kautz and MartinisKautz and Martinis (1990) and is described by the following set of equations

$$\beta \frac{d^2\phi }{d\tau ^2}+\frac{d\phi }{d\tau }+\alpha \left(\frac{d\phi }{d\tau }v\right)+\mathrm{sin}\phi =j+j_T+j_{Ts},$$
(7)
$$\beta \frac{dv}{d\tau }=\rho \left(\frac{d\phi }{d\tau }v\frac{j_{Ts}}{\alpha }\right),$$
(8)
where

$$j_{Ts}(\tau )j_{Ts}(\tau ^{})=2\alpha \gamma \delta (\tau \tau ^{}).$$
(9)
Here the parameters $`\alpha =R/R_s`$ and $`\rho =\alpha C/C_s`$ are introduced, $`v(t)`$ is the voltage across the external capacitance.

From the geometry of our samples we assume that $`C_s`$ is the capacitance between the top and bottom contacts of the mesa and consequently $`C_sC`$. The external resistance $`R_s`$ is determined by a contact resistance and the resistance of the leads. For small mesas the single junction resistance $`R`$ is very large, up to 10-100 k$`\mathrm{\Omega }`$ and $`R_sR`$ can be assumed. Therefore for small junctions the coupling to the external circuit should be included in the model calculations. In the classical fluctuation regime discussed in our paper it leads to the existence of a phase-diffusion branch instead of a zero-voltage superconducting branch.

Using this model as basis for numerical simulations we discuss in the following the typical behaviour of current-voltage curves in the different temperature regimes. In Section IV the model will be refined, in order to describe the experimental results in detail.

### II.3 Premature switching at low temperatures ($`T<T^{}`$).

At low temperatures a typical $`I`$-$`V`$ curve as shown in Fig. 3 contains a superconducting branch (S-state) with zero voltage and a hysteretic jump to the resistive state (R-state). Due to fluctuations the critical (or, better to say switching) current $`J_s`$ is random: the first strong enough thermal kick leads to a premature transition to the R-state, the probability of a back-transition to the S-state is very small. The average switching current $`J_s(T)`$ depends on the sweeping velocity $`\dot{I}`$ of the current. It can be calculated numerically and analytically. We found that the results of our numerical simulations are in good agreement with the well known analytical result $`J_c(T)`$ at low temperaturesTinkham (1996)

$$J_s(T)=J_c(T)\left[1\left(\frac{\gamma }{2}\mathrm{ln}\frac{\omega _p\dot{I}}{2\pi I}\right)^{2/3}\right],$$
(10)
Typically the logarithm is of the order of 10.

### II.4 Phase diffusion branch at intermediate temperatures ($`TT^{}`$)

At intermediate temperatures the S-branch shows a bending with a finite voltage close to the switching current (see Fig. 4). As mentioned above this behaviour can be explained by phase diffusion within a model with frequency-dependent damping. Let us discuss briefly the properties of this model.

In the limit of zero frequency dissipation is determined by the junction resistance $`R`$ alone, but at frequencies of the order of the plasma frequency $`\omega _p`$ due to the small impedance $`1/\omega _pC_s`$ of the external capacitance the dissipation is determined by the external resistance $`R_sR`$ in parallel with $`R`$. Therefore for small junctions with $`R_sR`$ dissipation at the plasma frequency is significantly larger then at zero frequency. For large-$`\beta `$ hysteretic junctions it leads to the coexistence of a finite-voltage phase-diffusion S-branch and resistive branch. Thermal current-kicks stimulate jumps of the Josephson phase difference from one potential minimum to the next and therefore diffusion. On the other hand, the high dissipation at the plasma frequency – the characteristic frequency of these jumps, prevents from a fast transition to the resistive state with a running phase.

In Fig. 4 the phase diffusion branch in the model with frequency-dependent damping is shown (circles) together with the usual RSJ curve at similar parameters (triangles).

The phase diffusion branch can be found numerically, and has been observed in our experiments, only at intermediate temperatures. At low temperatures the voltage in the S-state is exponentially small. At high temperatures the usual phase diffusion takes place with a single-valued voltage-current curve.

### II.5 Switching current at high temperatures ($`T>T^{}`$).

At high temperatures the hysteresis disappears and only one single-valued $`I`$-$`V`$ curve with a finite voltage at all currents exists (Fig. 5), but the average supercurrent, which is large at low voltages leads to a pronounced nonlinearity of $`I`$-$`V`$ curve. We can extract approximately the contribution of the supercurrent to the $`I`$-$`V`$ curve by subtracting the quasi-particle current, which in this model is a linear function of voltage. The current $`J`$ at which the difference $`JV/R`$ reaches its maximum (see Fig. 5) defines the switching current $`J_s`$. Note that this procedure to determine a switching current can also be applied at low temperatures and will be generally used later when experimental and theoretical results are presented and compared.

## III Sample preparation and experimental set-up

To measure the properties of ultrasmall intrinsic Josephson junctions we used a mesa geometry as sketched in Fig. 6.

Bi2Sr2CaCu2O8+δ (BSCCO) single crystals and Tl2Ba2CaCu2O8+δ (TBCCO) thin films were used to prepare the samplesKoval et al. (2000). Fig. 7 shows an SEM picture of a TBCCO sample. With our technique we were able to structure samples with lateral dimensions as small as $`0.2\times 0.2\mu m^2`$. The height of the mesas was between 60 Å and 1500 Å, which is equivalent to 4 - 100 intrinsic Josephson junctions.

Several small mesas with lateral dimensions of a few $`\mu m`$ or lower were arranged on one chip. On the same chip also some larger mesas with lateral dimension of several $`\mu m`$ were prepared. They serve as electrical contacts to the bulk superconductor. All mesas were contacted with Au leads. Current flows through one large mesa into the bulk superconductor and is extracted through the small mesa or vice versa. For voltage measurements one of the other large mesas is used. Due to the dimensions of the small mesa it was not possible to attach two isolated Au leads. Thus transport measurements were performed in a three-terminal configuration.

The Au electrode connected to the mesa leads to a suppression of the critical current and critical temperature of the uppermost intrinsic Josephson junctionKim et al. (1999); Rother et al. (2003). This degeneration might be explained by the assumption that the first superconducting layer is in close proximity with the normal Au electrode.

One sample was prepared in a step like geometry as proposed by Wang et al.Wang et al. (2001). Here we are able to record the $`I`$-$`V`$-characteristics in a true four-point geometry (sample s7n4 jj1 Table 1).

The bias current was provided by a battery powered current source. The $`I`$-$`V`$-characteristics of the samples were recorded by applying dc currents and recording the voltages across mesas by digital voltmeters.

Measurements were done either in a standard helium dewar or in a magnet cryostat equipped with a 5 T Helmholtz solenoid. Temperatures could be varied between 4.2 K and room temperature. The temperature was measured with a platinum resistor for higher temperatures and a cernox resistor for low temperatures.

To reduce external noise low pass filters were used. In the helium dewar both cold filters and room temperature filters were used.

## IV Experimental results and discussion

### IV.1 Experimental voltage-current curves and $`J_s(T)`$

We systematically measured voltage-current curves of mesas with different zero-temperature critical currents $`J_c(0)`$ at different temperatures (only a part of the sample parameters is presented in Table. I). Transition from high-hysteretic curves at low temperatures to phase-diffusion branches at crossover temperatures, and to a single resistive curve at high temperatures was clearly observed. Examples of experimental current-voltage curves at low (squares) and intermediate (filled circles) temperatures are shown in Fig. 8 (here the contact resistance is not subtracted). The well-known multibranch structure is observed.

The measurements of critical current were done mostly for the first (superconducting) branch, such that all junctions are in the S-state at low current and only one junction with the lowest critical current switches to the R-state. For this branch the average switching current $`J_s`$ was determined. Its temperature dependence is shown in Fig. 9 for five samples with different values of $`J_c(0)`$ (open symbols) together with the results of the theoretical calculation (filled symbols, discussed in the next section). This picture is the main result of our work.

The main feature of all curves is a kink (or even minimum) in $`J_s(T)`$ at some temperature $`T_{exp}^{}`$. This kink is absent in the $`J_s(T)`$ curves of mesas with large values of $`J_c(0)`$ and is shifted to lower temperatures with decreasing $`J_c(0)`$. $`T_{exp}^{}`$ corresponds approximately to the condition $`E_J(T^{})=kT^{}`$, which strongly suggests thermal fluctuations as an origin of the observed features.

### IV.2 Theoretical voltage-current curves and $`J_s(T)`$

To describe the experimental result for the switching current $`J_s(T)`$ in detail it is necessary to take into account the temperature dependence of the parameters $`\beta `$ and $`\gamma `$. These parameters are known functions of:

(i) the critical current $`J_c(T)`$ without fluctuations, which is determined (qualitatively also in HTSC) by the Ambegaokar-Baratoff relation.

(ii) the quasiparticle resistance $`R(T)`$, which we determined experimentally from current-voltage curves at small voltages (crosses in Fig. 10). We also calculated $`R(T)`$ at $`V0`$ from microscopic theory for different models of interlayer tunneling (coherent and incoherent) and intra-layer scattering (Born and unitary). The best fit to the experimental results was obtained for incoherent tunneling in the unitary limit of intra-layer scattering with an elastic scattering frequency $`\nu =1`$ meV. For the numerical simulations the function $`g(T)=R(0)/R(T)`$ is approximated by the fit (circles in Fig. 10).

The results of the theoretical calculations are presented in Fig. 9 for three different zero-temperature critical currents $`J_c(0)`$, the zero-temperature McCumber parameter $`\beta (0)`$ is in the range of $`10^3`$ to $`10^4`$. For each temperature we calculated $`J_c(T)`$, $`\beta (T)`$, $`\gamma (T)`$, and then produced current-voltage curves by numerical simulations starting from zero external current $`J`$. From these results the current $`J`$ was determined, at which the average supercurrent $`J_{sup}=JV/R`$ reaches its maximum. This current was defined previously as switching current $`J_s`$. The experimental switching currents $`J_s`$ shown in Fig. 9 are defined in the same way. We find good agreement between experimental and theoretical curves. Also the current-voltage curves calculated at different temperatures show the same behaviour as the experimental curves.

From these results we conclude that the non-monotonic temperature dependence of $`J_s(T)`$ observed in experiment can well be explained by the model of large-$`\beta `$ Josephson junction with thermal noise and coupling to an external impedance, if the proper temperature dependence of junction resistance is taken into account. In particular the strong suppression of $`J_s`$ with temperature at low temperatures, $`T<T^{}`$, is due to the premature switching effect in highly-hysteretic junctions. At high temperatures, $`T>T^{}`$, the McCumber parameter $`\beta `$ is smaller, the hysteresis disappears and the contribution of the average supercurrent to the current at finite voltages becomes larger, which leads to the plateau in the observed $`J_s(T)`$ dependence. Further discussions of the origins of non-monotonic temperature dependence of $`J_s`$ can be found in the papers of Iansiti et al.Iansiti et al. (1989) and Kautz and MartinisKautz and Martinis (1990).

### IV.3 Phase-diffusion at low currents

For the junctions with the smallest critical current $`J_c(0)1`$ $`\mu `$A and highest quasiparticle resistance $`RR_s`$ we investigated the low-current part of the current-voltage curves more precisely and found clear evidence of the phase-diffusion behaviour induced by the external circuit.

In the case $`R_sR`$ and for $`E_J<kT`$ the phase-diffusion branch can be described by the analytical expressionIngold et al. (1994); Tinkham (1996)

$$J(V)=\frac{4eTJ_s}{\mathrm{}}\frac{Z_1V}{V^2+(2eZ_1T/\mathrm{})^2},$$
(11)
and the small finite resistance in the S-state at very low currents and at any temperature is given byIngold et al. (1994)

$$R(J0)=\frac{Z_1}{I_0^2(1/\gamma )1},$$
(12)
where $`Z_1R_s`$ is the impedance at low frequencies (of the order of the plasma frequency $`\omega _p`$) and $`I_0(x)`$ is a Bessel function.

The experimental points presented in Fig. 11 can be fitted by Eq. (11). The best fitting result was achieved for a temperature dependent impedance $`Z_1`$ as shown in the insert. The origin of this temperature dependence is still unclear. It may be a contribution of the temperature dependent quasiparticle conductivity of other junctions in the mesa: when the temperature increases this conductivity also increases and the impedance decreases.

### IV.4 $`J_s(T)`$ in external magnetic field

Finally, we present here some results of measurements in an external magnetic field. Magnetic fields were applied by a 5 T superconducting split coil magnet. To measure $`J_s(B)`$ for parallel fields, the CuO2-layers must be aligned with high precision parallel to the magnetic field. If the magnetic field has a component perpendicular to the layers pancake vortices enter the Josephson junctions and suppress the critical current. The orientation of samples in a magnetic field follows the procedure described in Ref. Irie et al., 2000. The junction is biased at a fixed current and a fixed magnetic field is applied. Then the voltage change is monitored while the angle $`\phi `$ between the external field and the CuO2-layers was varied with a precision rotation stage. When measured at high enough temperatures the field-induced changes of the $`I`$-$`V`$-characteristic were reversible and could be used for alignment. Thus, the field orientation relative to the CuO2-layers could be adjusted to an accuracy of $`0.01^{}`$.

The main and most general experimental result is presented in Fig. 12. In both perpendicular and parallel magnetic fields $`J_s(T)`$ is not changed for $`T>T^{}`$, while for $`T<T^{}`$ one finds a field dependent suppression of $`J_s(T)`$.

To understand this behaviour qualitatively, we explore the well-known sin-Gordon equation, which describes long Josephson junctions (see Ref. Barone and Paterno, 1982; Likharev, 1996)

$$\beta \frac{^2\phi }{d\tau ^2}\frac{^2\phi }{dx^2}+\frac{\phi }{d\tau }+\mathrm{sin}(\phi +B_{}x+\delta \phi (x,\tau ))=j+j_T(x,\tau ).$$
(13)
Here the coordinate $`x`$ is in units of the Josephson length $`\lambda _J`$, $`B_{}`$ is the component of the external magnetic field parallel to the layers, and $`\delta \phi (x,\tau )`$ is a random phase distribution produced by pancake-vortices in perpendicular magnetic fieldBulaevskii et al. (1996b); Bulaevskii et al. (1997); Koshelev et al. (2000). We include the phase variation due to the external magnetic field in the $`sin`$-function instead of using the more familiar procedure, where the magnetic field is taken into account by boundary conditions (both methods are equivalent). The random functions $`j_T(x,\tau )`$ and $`\delta \phi (x,\tau )`$ describing noise are assumed to be Gaussian and $`\delta `$-correlated in space

$$j_T(x,\tau )j_T(x,\tau ^{})=2\gamma \delta (\tau \tau ^{})\delta (xx^{}),$$
(14)
$$\delta \phi (x,\tau )\delta \phi (x,\tau ^{})B_z\delta (xx^{}).$$
(15)
In the last expression we accept a linear magnetic field dependence of the correlation function and neglect the time-dependence of phase shifts assuming that the pancakes are frozen or move slowly enough. In a more realistic approach the space and time dependence of correlation functions of the pancake disorder depend on a state of the vortex matter, is it liquid, glass-like, or solid.

Results of our numerical simulations with magnetic field are presented in Fig. 13 for several values of perpendicular and parallel fields. The main result of our calculations (as well as similar experimental results) is the sensitivity to a magnetic field at low temperatures and the absence of considerable influence at high temperatures.

Qualitatively this result can be understood in the following way. The main effect of the magnetic field is to produce phase shifts, regular for the parallel field and random for the perpendicular field. At low temperatures, when $`J_s`$ is determined by the premature switching mechanism, the reduction of $`J_c`$ in the magnetic field is important because the energy barrier becomes smaller. But in the high-temperature regime where phase diffusion is strong, additional phase shifts $`\delta \phi (x,\tau )`$ produced by the magnetic field cannot compete with the existing large thermal phase fluctuations.

## V Summary

In this paper we discuss new experiments on the influence of thermal noise on current-voltage curves and the switching current of small area intrinsic Josephson junctions. Two main effects of thermal fluctuations: premature switching and phase-diffusion, are observed. The non-monotonic temperature dependence of the switching current is explained theoretically with a help of an extended RSJ model (including thermal noise and a frequency dependent external impedance) taking into account the proper temperature dependence of parameters. The different behaviour of the switching current at low and high temperatures in the presence of magnetic fields can be explained in a similar way.

## Acknowledgements

We thank G. Blatter and L.N. Bulaevskii for valuable discussions.

This work was supported by the German Science Foundation (D.R.), and by the Swiss National Center of Competence in Research ”Materials with Novel Electronic Properties-MaNEP” (C.H.).