text:
# High-field side of Superconductor-Insulator Transition

## I Introduction

The zero-temperature superconductor-insulator transition (SIT) is a phase transition of the second order and is driven by fluctuations of a quantum nature. The superconducting phase is considered to be a condensate of Cooper pairs with localized vortices, and the insulating phase is a condensate of vortices with localized Cooper pairs. The theoretical description based on this assumption was suggested in Ref. . According to the theory: (i) the film resistance $`R`$ near the field-induced SIT at low temperature $`T`$ in the vicinity of the critical field $`B_c`$ is a function of one scaling variable $`\delta =(BB_c)/T^{1/\nu z}`$, with $`\nu z`$ being a parameter; (ii) at the transition point the film resistance has the universal value $`h/4e^2`$ (the quantum resistance for Cooper pairs). Although much work has been done, and in many systems the scaling relations hold, the SIT in disordered films remains a controversial subject, especially concerning the insulating phase and the bosonic conduction on the high-field side of the SIT at $`B>B_c`$. In some works the behavior of the resistance in this region can be explained by weak-localization theory. In Ref.  it has been found that, on the high-field side of the transition, the magnetoresistance reaches a maximum and the phase can be insulating as well as metallic, with at high fields the zero-temperature asymptotic value of the resistance is approximately equal to $`R_c`$ (the value of the resistance at $`B=B_c`$). Here we present a careful examination of the presence of the insulating phase and its evolution on the high-field side of the SIT.

## II Experiments and discussion

A TiN film with a thickness of 5 nm was formed on 100 nm of SiO2 grown on top of $`<`$100$`>`$ Si substrate by atomic layer chemical vapour deposition at 400C . Analysis by atomic force spectroscopy shows that the film exhibits low surface roughness and consists of a dense packing of grains, with a rather narrow distribution of grain size and the average size is roughly $`30`$ nm.

The samples for the transport measurements were fabricated into standard Hall bridges using conventional UV lithography and subsequent plasma etching. The channel length and width were 250 and 50 $`\mu `$m, respectively. The magnetoresistance measurements were performed in a temperature-stabilized dilution refrigerator. The magnetic field was applied perpendicular to the film. Four terminal transport measurements were performed using standard low frequency techniques. The resistance data were taken at a measurement frequency of 10 Hz with an ac current $`1`$ nA.

A typical set of resistance per square vs magnetic field data is given in Fig. 1. The main features of these results are the presence of the point of intersection ($`B_c,R_c`$) and the negative magnetoresistance in high fields. Using the $`B_c`$, we plot the same data against the scaling variable $`|BB_c|/T^{1/\nu z}`$, and adjust the power $`\nu z`$ to obtain the best visual collapse of the data. Such behavior, previously regarded as the main evidence of the existence of SIT, is actually not incontestable proof of the presence of the insulating phase at $`B>B_c`$. Following the approach of the authors in , we analyze the temperature dependence of the film resistance on the high-field side of the SIT in terms of 3D “bad” metals in the vicinity of the metal-insulator transition :

$$G(T)=a+bT^{1/3}.$$
(1)
The sign of the parameter $`a`$ discriminates between a metal and an insulator at $`T0`$. If $`a>0`$, it yields the zero temperature conductance $`G(0)a`$, whereas negative $`a`$ points to activated conductance at lower temperatures. The temperature dependence of the high-field conductance at different magnetic fields is shown in Fig. 2. It is well described by Eq. (1). As $`G(0)`$ is negative in fields higher than the critical field we can conclude that this phase is insulating. With further increase of $`B`$ the resistance decreases. This behavior is in agreement with the concept of localized Cooper pairs in the insulating phase, and that of field-induced Cooper-pair breaking. Extrapolation to $`T=0`$ allows us not only to determine the field of the insulator-metal transition ($`B_{IM}`$), but also the magnetic field dependence of the zero temperature conductance. The result of this procedure is presented in the lower inset to Fig. 2. The zero temperature conductance at $`B>B_{IM}`$ is well described by an empirical expression $`G(B)=e^2/h(1\mathrm{exp}[(B_{IM}B)/B^{}])`$ shown by the solid line in the lower inset to Fig. 2. In contrast to the statement of the authors in , the resistance at $`T0`$ on the high-field side of the SIT approaches the value of $`h/e^2`$ rather than $`R_c`$. The exponential dependence of $`G(T=0,B)`$ may result from a broad dispersion of the binding energies of localized Cooper pairs.

This work was supported by Grants RFBR 00-02-17965, RFBR 02-02-16782 and by the programs “Low-dimensional and Mesoscopic Condensed Systems” of the Russian Ministry of Science, Industry, and Technology, and “Quantum Macrophysics” of the Russian Academy of Science.


format: text

text:
# Liquid-Gas Phase Transition and Instabilities in Asymmetric Two Component Systems

## Abstract

The liquid-gas phase transition and associated instability in two component systems are investigated using a mean field theory. The importance of the role of both the Coulomb force and symmetry energy terms are studied. The addition of the Coulomb terms bring asymmetry into a mean field and thus result in important differences with previous approaches which did not include such terms. The Coulomb effects modify the chemical instability and mechanical instability domains shifting many features away from proton fraction point $`y=1/2`$ to a value closer to the valley of $`\beta `$ stability. These features are discussed in detail.

Thermodynamic properties of two and, more generally, multi-component systems are of interest in many areas of physics. For example, binary systems of two components are encountered in nuclei where the two components are neutrons and protons. While light nuclei have $`N=Z`$, heavy nuclei have $`N>Z`$ and systems with large neutron excess such as neutron stars exist in astrophysics. In condensed matter physics two component systems appear in liquid 3He where the two components are spin-up and spin-down fluids. Such systems contain a more complex phase structure than one component systems. This letter is an extension of previous work concerning the nuclear case which was discussed from both the fragmentation point of view  and from a mean field picture . In the mean field picture we studied the phase surface associated with the liquid-gas phase transition incorporating the role of surface and Coulomb effects . Previous mean studies  of two component systems were investigations of the role of the symmetry energy on the phase structure of the liquid-gas phase transition. While the symmetry energy plays an important role in governing the proton/neutron ratios in the liquid and gas phases, favoring more symmetric systems in the liquid phase than the gas phase, the effects of both the Coulomb interaction and surface energry were also shown to be important. This was first discussed in Ref. and more recently in Ref.. Here we will focus on the stability, metastability and instability of asymmetric systems of protons and neutrons. Such systems are of interest for radioactive beams studies which can be used to explore nuclei at the limits of isospin asymmetry. In a recent article, V. Baran etal  have pointed out some new effects associated with two component systems where a new kind of liquid-gas phase transition associated with chemical instability may appear. This new type of transition may manifest itself as an isospin distillation or fractionization as reported in Ref. where the gas phase is enriched in neutrons relative to the liquid phase. Such observations suggest that neutron diffusion occurs at the same rate or is more rapid than fragment production. These observations are based solely on symmetry energy effects, which make the liquid phase more symmetric, with the role of the Coulomb interaction not included. The importance of Coulomb effects driven by $`Z^2`$ terms versus symmetry energy effects arising from $`(NZ)^2`$ terms can be seen by just looking at the most stable nucleus for a given $`A`$. Using a liquid drop model, the $`Z`$ with the smallest mass or greatest binding energy is given by $`Z=(A/2)/(1+B_cA^{2/3}/B_s)`$ where $`B_c`$ is the Coulomb coefficient $`B_cZ^2/A^{1/3}`$ and $`B_s`$ is the symmetry energy coefficient $`B_s(NZ)^2/A`$. When the Coulomb interaction is turned off $`Z`$ is simply $`A/2`$. Moreover, the Coulomb interaction was shown  to have a very dramatic effect on the specific heat associated with a first order liquid-gas phase transition, substantially reducing the peak previously reported in Ref. where surface energy effects associated with one component systems played the dominant role. Moreover, the Coulomb force introduced another pair of coexistence densities having larger proton concentration in gas phase than in liquid phase due to its asymmetry in the mean field .

In symmetric systems, mechanical instability from $`dP/d\rho <0`$ and chemical instability from $`d\mu /d\rho <0`$ at fixed $`T`$ are simply related since

$`Nd\mu =SdT+VdP`$ (1)
where $`\mu `$ is the chemical potential, $`S`$ the entropy, $`P`$ the pressure, $`N`$ the number of particles and $`T`$ the temperature. Thus $`\rho d\mu /d\rho =dP/d\rho `$ at fixed $`T`$ and the boundaries of both instabilities are common with $`d\mu /d\rho =dP/d\rho =0`$. The coexistence curve of the Maxwell pressure $`P`$ versus density $`\rho `$ or volume $`V`$ is determined by the equality of the pressure and chemical potentials between the liquid and gas phase at the same $`T`$. In the region between the spinoidal curve determined by the instability line and the coexistence line are metastable regions of supercooled liquid and superheated gas where nucleation processes can occur. These processes for one component systems were studied in Ref. , where nucleation processes gave a dynamical picture of the liquid-gas phase separation. Specifically, classical nucleation theory  looks at the change in the Gibbs free energy with and without a liquid drop made of $`A`$ particles. This change is written as $`\mathrm{\Delta }G=(\mu _l\mu _g)A+4\pi R^2\sigma `$ where $`\mu _l`$ is the liquid chemical potential of the drop, $`\mu _g`$ is the gas chemical potential from which the drop is formed, $`R`$ the radius of the drop and $`\sigma `$ is the surface tension. A superheated gas has $`\mu _g>\mu _l`$, while equilibrium states have $`\mu _g=\mu _l`$. This change in Gibbs free energy has a maxima at a critical size drop whose radius is given by $`R_c=2\sigma /(\mu _g\mu _l)\rho _l`$, where $`\rho _l`$ is the liquid drop density. The $`A=\rho _lV`$ with $`V`$ the volume of the drop. Drops with smaller $`R`$ evaporate particles thereby lower their free energy, while drops with bigger $`R`$ grow by accumulating particles and also lowering their free energy. Such processes in two component systems of neutrons and protons with additional Coulomb and symmetry energy effects have also been partly explored. Some preliminary studies of Coulomb effects in nucleation were done in Ref.. For a two component drop, the $`\mathrm{\Delta }G`$ is

$`\mathrm{\Delta }G=(\mu _{lp}\mu _{gp})Z+(\mu _{ln}\mu _{gn})N+4\pi R^2\sigma +\delta E_C+\delta E_{sym}`$ (2)
The $`\delta E_C`$ is the change in Coulomb energy with and without the drop and $`\delta E_{sym}`$ is the change in symmetry energy.

In two component systems the coexistence line of a one component system now becomes a surface in $`(P,y,T)`$, with $`y`$ the proton fraction. This surface has several new elements associated with it that are not present in one component systems. For example, the isolated critical point of a one component system now becomes a line of critical points with $`y=1/2`$ the critical point of a one component system. Associated with this surface is also a line of maximal asymmetry in the proton/neutron ratio and a line of equal concentration in liquid and gas phase . The intersection of this $`(P,y,T)`$ with a plane at fixed $`T`$ leads to loops of $`P`$ versus $`y`$ as shown in Fig. 4 of Ref.. The chemical potentials $`\mu _n`$ and $`\mu _p`$ versus $`y`$ at fixed $`T`$, $`P`$ also have interesting behavior as shown in Fig. 3 of Ref. . For two component system, Eq.(1) is replaced by

$`y(d\mu _p/d\rho )+(1y)(d\mu _n/d\rho )=(1/\rho )(dP/d\rho )`$ (3)
$`y(d\mu _p/dy)+(1y)(d\mu _n/dy)=(1/\rho )(dP/dy)`$ (4)
In the first equation, $`y`$ and $`T`$ are held fixed, and in the second equation $`\rho `$ and $`T`$ are held fixed. We will use these figures to study the spinoidal instability region and its associated surface in $`(P,y,T)`$.

Some properties of the instability region can best be analyzed using some thermodynamic identities. We first note that the pressure and the chemical potentials are usually written as functions of $`\rho `$, $`y`$ and $`T`$; $`P(\rho ,y,T)`$ and $`\mu _q(\rho ,y,T)`$. One important identity arises from the equation of state $`P=P(\rho ,y,T)`$ which gives

$`\left({\displaystyle \frac{dy}{d\rho }}\right)_{P,T}=\left({\displaystyle \frac{dP}{d\rho }}\right)_{y,T}/\left({\displaystyle \frac{dP}{dy}}\right)_{\rho ,T}`$ (5)
for the relation of $`y`$ and $`\rho `$ at fixed $`P`$ and $`T`$. Since $`P`$, $`T`$, $`y`$ and $`\rho `$ are connected by an equation of state $`P(\rho ,y,T)`$, we can write

$`\left({\displaystyle \frac{d\mu _q}{dy}}\right)_{P,T}=\left({\displaystyle \frac{d\mu _q}{dy}}\right)_{\rho ,T}+\left({\displaystyle \frac{d\mu _q}{d\rho }}\right)_{y,T}\left({\displaystyle \frac{d\rho }{dy}}\right)_{P,T}`$ (6)
Combining Eqs.(5) and (6), we can write $`(dP/d\rho )_{y,T}`$ as

$`\left({\displaystyle \frac{dP}{d\rho }}\right)_{y,T}={\displaystyle \frac{(dP/dy)_{\rho ,T}(d\mu _q/d\rho )_{y,T}}{(d\mu _q/dy)_{P,T}(d\mu _q/dy)_{\rho ,T}}}`$ (7)
Since $`(d\mu _q/dy)_{\rho ,T}`$ is finite, singularites in $`(d\mu _q/dy)_{P,T}`$ (chemical instability) give rise to zeroes in $`(dP/d\rho )_{y,T}`$. The zeroes of $`(dP/d\rho )_{y,T}`$ form the spinodal surface, boundary for mechanical instability region. In two component systems, mechanical and chemical instability regions are separated for asymmetric cases. Using Eq.(5), the chemical instability boundary ($`(d\mu _q/dy)_{P,T}=0`$), determined by setting Eq.(6) equal to 0, is obtained from

$`\left({\displaystyle \frac{dP}{d\rho }}\right)_{y,T}=\left({\displaystyle \frac{dP}{dy}}\right)_{\rho ,T}{\displaystyle \frac{(d\mu _q/d\rho )_{y,T}}{(d\mu _q/dy)_{\rho ,T}}}`$ (8)
The following features should be noted. First, mechanical equilibrium ($`(dP/d\rho )_{y,T}=0`$) is determined by the condition that the left hand side of this equation is 0. Thus the mechanical and chemical instability boundaries are the same when $`(dP/dy)_{\rho ,T}=0`$ since $`(d\mu _q/dy)_{\rho ,T}0`$ here. Due to Eq.(5), $`(d\rho /dy)_{P,T}=0`$ too when $`(dP/dy)_{\rho ,T}=0`$ except on the mechanical instability boundary. This condition is also related to the equal concentration $`y_E`$ of liquid and gas phase on the coexistance surface . For systems without Coulomb terms this occurs at $`y=1/2`$. Secondly, for chemical equilibrium the right hand side can be obtained from the equation of state $`P=P(\rho ,y,T)`$ and the behavior of the proton or neutron chemical potentials $`\mu _q=\mu _q(\rho ,y,T)`$. The $`q`$ will distinguish neutron results from proton results. We will use Eq.(8) to explore the boundary of chemical equilibrium.

We begin our study by first considering a simpler system where Coulomb and surface effects are turned off. Using Eq.(54) of Ref. with $`C_s=0`$ (surface term), $`C=0`$ (Coulomb term), $`x_3=1/2`$ and $`\alpha =1`$, then $`(dP/d\rho )_{y,T}`$ is quadratic in $`\rho `$ and is given by

$`\left({\displaystyle \frac{dP}{d\rho }}\right)_{y,T}`$ $`=`$ $`T+\left({\displaystyle \frac{3}{4}}\right)t_0\rho +\left({\displaystyle \frac{3}{8}}\right)t_3\rho ^2+T\left({\displaystyle \frac{1}{4\sqrt{2}}}\right){\displaystyle \frac{\lambda ^3}{\gamma }}\rho `$ (10)
$`\left[t_0(1/2+x_0)T\left({\displaystyle \frac{1}{2\sqrt{2}}}\right)\left({\displaystyle \frac{\lambda ^3}{\gamma }}\right)\right]\left({\displaystyle \frac{\rho }{2}}\right)(2y1)^2`$
The last term involving $`(2y1)^2`$ does not contribute for symmetric systems. The terms involving $`\lambda ^3`$ arise from degeneracy corrections. The mechanical instability points can easily be obtaind by solving the quadratic equation $`(dP/d\rho )_{y,T}=0`$ to obtain the two densities that define the instability region for each $`T`$ and $`y`$. These two roots become a single solution at a specific $`T_c(y)`$ at each $`y`$. This occurs when $`(d^2P/d\rho ^2)_{y,T}=0`$. In one component systems and in symmetric systems with $`y=1/2`$, this is the critical point $`T_c`$.

The calculations using Eqs.(53) and (54) of Ref.  with degeneracy correction are shown in Fig. 1 and Fig. 2. The mean field force parameters are the same as used in Ref. except for using a larger radius $`R=8`$ fm of a uniform sphere for a larger Coulomb force. These figures show that the mechanical instability boundary (thin solid line) is inside of the chemical instability region (thick solid line) and these coincide at the proton concentration of $`y_E=0.4256`$ at $`T=10`$ MeV (Fig.1 and upper middle box of Fig.2). The curves of $`(d\mu _q/d\rho )_{y,T}=0`$ also coincide with the instability boundary at $`y_E`$. The coexistance boundary coincides with the chemical instability boundary at the line of critical points (peaks in Fig.1 for a $`T=10`$ MeV plane cut) but not with the mechanical instability boundary . We can see from Fig.2 that the chemical and mechanical instability boundaries coincide at the proton concentration $`y_E(T)`$ which is on the line of equal concentration on the coexistance surface . While all the curves coincide at for $`y=1/2`$ (without Coulomb) as shown in the lower left box of Fig.2, the top middle box of Fig. 2 with Coulomb exhibit a separation between these curves due to the temperature dependence of $`y_E(T)`$. All curves in the upper middle box of Fig. 2 coincide at $`T=10`$ MeV since $`y_E=0.4256`$ at this temperature. From Fig.2 and the left column of Fig.1, we can see that all the four surfaces in $`(\rho ,y,T)`$ space coincide at the line of temperature dependent $`y_E(T)`$ and that the mechanical instability surface is inside of the chemical instability surface. Upper Fig. 2 and upper left box of Fig.1 also shows that the role of $`(d\mu _p/d\rho )_{y,T}=0`$ and $`(d\mu _n/d\rho )_{y,T}=0`$ are exchanged when $`y`$ crosses $`y_E`$.

To keep equations simple in the dicussions of various new features in two component system with Coulomb force, we will drop the degeneracy corrections. However, they are included in the actual calculations done and in the fiqures (Figs. 1 and 2) presented here. In this limit $`P=\rho Ta_0(y)\rho ^2+2a_3\rho ^3`$ where we have followed the notation of Ref.  with $`P`$ now generalized to the case of asymmetric systems. The generalization changes $`a_0`$ to $`a_0(y)=a_0[1(2/3)(1/2+x_0)(2y1)^2]`$ with $`a_0=(3/8)t_0`$ and $`a_3=t_3/16`$. Then $`dP/d\rho =0`$ at

$`\rho _\pm =\left(a_0(y)\pm \sqrt{a_0^2(y)6a_3T}\right)/(6a_3)`$ (11)
The two roots are equal at $`T_c(y)=a_0^2(y)/(6a_3)`$ and are $`\rho _c(y)=\rho _+=\rho _{}=a_0(y)/(6a_3)`$. At this point the pressure is $`P_c(y)=a_0^3(y)/(108a_3^2)`$. Under the same limits the chemical potentials are given by

$`\mu _q=T\mathrm{ln}[(\lambda ^3/\gamma )(\rho /2\pm (2y1)\rho /2)]2a_0\rho \pm (4/3)(1/2+x_0)a_0(2y1)\rho +3a_3\rho ^2`$ (12)
with the upper sign for protons and the lower sign for neutrons and $`\gamma =2`$ for the spin up and down degeneracy. This equation for $`\mu _q=\mu _q(\rho ,y,T)`$ when combined with the equation of state $`P(\rho ,y,T)`$ can be used to calculate the chemical potentials as a function $`\mu _q=\mu _q(P,y,T)`$. Boundaries of the chemical instability region can also be determined using the condition $`d\mu _q/dy=0`$ at constant $`T`$, $`P`$. Substituting the simplified equation of state and expression for the chemical potential into Eq.(8) leads to a cubic equation $`c_3\rho ^3+c_2\rho ^2+c_1\rho +c_0=0`$ with $`c_3=32(1/2+x_0)a_0a_3(1y)y`$, $`c_2=[6a_3T(32/3)(1/2+x_0)a_0^2(1y)y]`$, $`c_1=(4/3)(1x_0)a_0T`$, $`c_0=T^2`$. Chemical instability occurs when three real roots exist, two of which will be postive, one negative. One negative root is present because $`c_0/c_3`$ is positive. At $`y=1/2`$, the two positive solutions are the same as that for mechanical instability. When $`y`$ is no longer 1/2, the mechanical instability and chemical instability regions are no longer the same. This result is shown in the lower part of Fig.2 for various choices of $`y`$ with the degenercy correction of Ref. . It should also be noted that the coefficients of the cubic equation have a certain symmetry arising from the form in $`y`$ which is $`y(1y)`$, i.e., interchanging $`y`$ (protons) and $`1y`$ (neutrons) leaves the solution invariant (see curves for $`(d\mu _q/d\rho )_{y,T}=0`$ in Figs. 1 and 2). This result is no longer true when Coulomb interactions are included. We now turn our attention to the simple EOS and $`\mu _q`$, but now with Coulomb effects included in a simplified form which will illustrate its qualitative effects.

When Coulomb effects are included in the approximation of a uniformly charged sphere of fixed radius $`R`$, the pressure equation contains an additional term $`Cy^2\rho ^2`$ and the proton chemical potential contains an extra term $`2Cy\rho `$. The $`C`$, which depends on R, is taken to be a constant. For this choice, we can consider R as an effective range of the Coulomb force which is now independent of the density. This simplification reduces the complexity of many results. To make the Coulomb effects more visible in our figures (Figs. 1 and 2), we use a large value of $`C=231`$ MeV/fm3 which is appropriate for systems near U. These additions to $`P`$ and $`\mu _p`$ change the equation for chemical and mechanical instability. Specifically the coefficients of the cubic equation for chemical instability for both protons and neutrons are now $`c_3=[32(1/2+x_0)a_0+12C]a_3y(1y)`$, $`c_2=6a_3T[(32/3)(1/2+x_0)a_0^2+(8/3)(1x_0)a_0C]y(1y)`$, $`c_1=[(4/3)(1x_0)a_0+2Cy]T`$, $`c_0=T^2`$.

Some qualitative features of the effects of such a Coulomb behavior are the following (see Fig. 1 and Fig.2 also). 1) The mechanical and chemical instability curves are the same when $`(dP/dy)_{\rho ,T}=0`$ which occurs at $`y_E=0.5/(1+C/[(8/3)(1/2+x_0)a_0])`$ and thus no longer at $`y_E=0.5`$ when $`C0`$. This equation for $`y`$ has a similar structure to the equation that gives the $`Z`$ which has the maximum binding energy at fixed $`A`$ (valley of $`\beta `$ stability) obtained from the liquid drop mass formulae or from Eq.(55) of Ref.. With higher order $`T`$ dependence coming from degeneracy corrections, $`y_E`$ depends on $`T`$ and $`y_E=0.4256`$ for $`x_0=1/2`$ at $`T=10`$ MeV. Explicitly, with first order degeneracy corrections, $`y_E=0.5/(1+C/[(8/3)(1/2+x_0)a_0+T\lambda ^3/(2\sqrt{2}\gamma )])`$. 2) For both larger and smaller values of $`y`$ from $`y_E`$, the chemical instability curve is outside the mechanical instability curve as shown in Fig.1 and 2. The mechanical instability has its peak in $`T`$ versus $`\rho `$ at $`T_c(y)=[a_0(2/3)(1/2+x_0)a_0(2y1)^2Cy^2]^2/(6a_3)`$ which is maximum at $`y=y_E`$ and coincide with chemical instability at $`y_E`$. 3) The maximum critical $`T`$ obtained from the peak in the chemical equilibrium surface of $`T`$ versus $`\rho `$, occurs at $`y0.5`$, but is now shifted to $`y=y_E(T)`$. 4). At $`y=y_E(T)`$, the equilibrium curve or binodal curve of $`P`$ versus $`y`$ obtained by cutting the binodal surface at a given $`T`$ meets at a point, i.e., the liquid and gas phase have the same proton concentration (see Fig.4 of Ref.). For $`y>y_E`$ the gas phase is proton richer than the liquid phase when compared to a situation without Coulomb effects . When Coulomb effects are turned off the binodal section meets at $`y=1/2`$. At $`y=y_E1/2`$ with Coulomb effects included a two component system will behave as a one component system, keeping the proton/neutron ratio the same in both liquid and gas phases.

In summary, we have presented calculations based on a mean field theory with asymmetry which include the role of Coulomb force on the liquid-gas phase instability regions for both the mechanical (spinodal) surface and the chemical (diffusive) surface. The Coulomb force modifies these surfaces so that their peaks are no longer at proton fraction $`y=1/2`$ but shifted to a smaller value called $`y_E`$ which is close to the proton fraction in the valley of $`\beta `$-stability. The mechanical instability surface is always inside the chemical instability surface and the two surfaces are the same at the equal concentration $`y_E`$ (allowing for $`T`$ dependence in $`y_E`$). At $`y_E`$, the peak in the chemical instability surface and binodal equilibrium surface also coincide. For systems with $`yy_E`$, the peak in the $`T`$ versus $`\rho `$ curve for chemical instability and chemical equilibrium – binodal or coexistance curve – coincide, but the peak for mechanical or spinodal instability is below that of the other two. At $`y=y_E`$, a phase transition in a two component system behaves the same as a one-component system, keeping the proton fraction at $`y_E`$ in both phases during the entire transformation.

This work was supported in part by the DOE Grant No. DE-FG02-96ER-40987 and in part by Grant No. R05-2001-000-00097-0 from the Basic Research Program of the Korea Science and Engineering Foundation.


text:
# Dependence of Nuclear Level Density on Vibrational State Damping

## 1  Introduction

Level density $`\rho `$ is one of the main quantities to define characteristics of nuclear decay. The collective states can strongly effect on level density, specifically, at low excitation energies -. Calculation of the variation factor $`K`$ of level densities is the simplest method to estimate effect of the vibrational states on level densities. The factor $`K`$ is the ratio of level densities with and without allowing for vibrational states. The level densities can be calculated within framework of statistical approach by standard saddle-point method or modificated one. It was found that factor $`K`$ is practically independent of type of statistical approach. Therefore we will use below standard saddle-point method as a simplest way for calculation of the variation factor of level density.

In this contribution work the effect of the vibration state damping on nuclear level density is studied on the base of response function approach in wide range of excitation energies. The dependence of variation factor $`K`$ on damping of vibrational state is investigated.

## 2  Response function method for calculation of the level density variation

Variation factor of level density is defined as

$$K=\rho /\rho _0,$$
(1)
where $`\rho `$ and $`\rho _0`$ are level densities with and without allowing for vibrational states.

The level density has the following form within framework of statistical approach in standard saddle-point method

$$\rho (U,A)=\left(4\pi ^2D\right)^{1/2}\mathrm{exp}S(\alpha _0,\beta _0),$$
(2)
where $`S(\alpha _0,\beta _0)=\alpha _0A+\beta _0E+\mathrm{ln}Z(\alpha _0,\beta _0)`$ is entropy of nucleus with mass number $`A`$ at excitation energy $`U=EE_{g.s}`$ with $`E_{g.s}`$ for ground state energy; $`Z(\alpha ,\beta )=\mathrm{Tr}\left[\mathrm{exp}\left(\beta \widehat{H}+\alpha \widehat{A}\right)\right]`$ \- partition function; $`\widehat{H}`$ \- nuclear Hamiltonian; $`\widehat{A}`$ \- operator of particle number; $`D`$ \- determinant of second partial derivatives of logarithm of partition function with respect to parameters $`\alpha `$ and $`\beta `$. The saddle-points $`\alpha _0,\beta _0`$ define temperature ($`T=1/\beta _0`$) and chemical potential ($`\mu =\alpha _0/\beta _0`$). They are solutions of the system of thermodynamic state equations $`A=\mathrm{ln}Z/\alpha `$, $`E=\mathrm{ln}Z/\beta `$.

We use conception of random phase approximation and suppose that collective vibrational states of the multipolarity $`L`$ are formed by the separable residual interaction $`V_{res}^k=k\widehat{Q}_L\widehat{Q}_L^+`$, where $`\widehat{Q}_L`$ is one-body operator of multipole momentum with multipolarity $`L`$, $`\widehat{Q}_Lr^LY_{LM}`$, and $`k`$ is a coupling constant.

The partition function can be presented in the form $`Z=Z_0\mathrm{\Delta }Z,`$ where $`Z_0`$ is partition function of the independent particle model. The quantity $`\mathrm{\Delta }Z`$ can be expressed by variation of thermodynamic potential $`\mathrm{\Delta }\mathrm{\Omega }`$ in the form $`\mathrm{\Delta }Z=\mathrm{exp}(\beta \mathrm{\Delta }\mathrm{\Omega })`$ with the following equation for $`\mathrm{\Delta }\mathrm{\Omega }`$

$$\mathrm{\Delta }\mathrm{\Omega }=\underset{0}{\overset{k}{}}𝑑k^{}\widehat{Q}_L\widehat{Q}_L^+/2=\frac{2L+1}{2\pi }\underset{0}{\overset{k}{}}𝑑k^{}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{1}{e^{\beta \epsilon }1}Im\mathrm{X}_L^k^{}𝑑\epsilon .$$
(3)
Here, $`\mathrm{X}_L^k^{}`$ is nuclear response function on external field proportional to $`\widehat{Q}_L`$ with frequency $`\omega =\epsilon /\mathrm{}`$; symbol $`<\mathrm{}>`$ denotes averaging over total Hamiltonian with residual interaction $`V_{res}^k^{}`$. The final expression in right-hand side of the eq.(3) is obtained after transformation of relationship for $`\mathrm{\Delta }\mathrm{\Omega }`$ with the use of Green function method. Finally variation of the partition function is presented in the following form

$$\mathrm{\Delta }Z\left(\beta \right)\mathrm{\Delta }Z_{vib}=\underset{j}{}\left|\left[\frac{1\mathrm{exp}\left(\left(E_j^0+i\mathrm{\Gamma }_j^0\right)\beta \right)}{1\mathrm{exp}\left(\left(E_j+i\mathrm{\Gamma }_j\right)\beta \right)}\right]^{2L+1}\right|$$
(4)
with the use of transformation similar to that one from Ref. Here, $`E_j`$, $`\mathrm{\Gamma }_j`$ are real and imaginary parts of solutions of the equation

$$k1/\mathrm{X}_L^0(\epsilon )=0,$$
(5)
and $`E_j^0`$, $`\mathrm{\Gamma }_j^0`$ are real and imaginary parts of roots of the equation $`1/\mathrm{X}_L^0(\epsilon )=0`$ with $`\mathrm{X}_L^0`$ for response function of independent particle model ($`k=0`$).

## 3  Calculations and conclusion

Semiclassical Vlasov-Landau approach for response function was used in the calculations. The nuclear mean field and relaxation processes can be easy taken into account in this approach. Collision integral was considered in relaxation time approach. Three different parametrizations of relaxation time $`\tau _c`$ were adopted with taken into account retardation effects.

1) Expression within framework of transport theory

$$\mathrm{}/\tau _c(\epsilon ,T)=\left(\left(\epsilon /2\pi \right)^2+T^2\right)/4.07\mathrm{MeV}.$$
(6)
2)Relaxation time due to model of doorway state

$$\mathrm{}/\tau _c(\epsilon ,T)=2\pi \left(\frac{K_m}{A^3}\right)\frac{g^3}{2}\frac{\left(\epsilon +U\right)^2}{3+\overline{n}},$$
(7)
where $`\overline{n}=0.843aT`$; $`K_m`$ is defined from agreement between relaxation times given by eqs.(6) and (7) at $`T=0`$ and at GQR energy; $`g=g\left(\epsilon _F\right)=6a/\pi ^2`$ with level density parameter $`a=A/8`$MeV-1.

3) Model of doorway state

$$\mathrm{}/\tau _c(\epsilon ,T)=2\pi \left(\frac{K_B}{A^3}\right)\frac{g^3}{8}\left(\epsilon +U\right),$$
(8)
where $`K_B`$ is defined from agreement with relaxation times given by eqs.(6) and (8) at $`T=0`$ and at GQR energy.

The dependence of variation factor $`K`$ on excitation energy is shown on fig.1 for nucleus 56Fe. Low-lying quadrupole state with energy $`E_{2^+,exp}=0.847`$MeV was taken into account; histogram is microscopic calculation and experimental values. Solid line shows calculation with relaxation time given by eq.(6), dash - with relaxation time according to eq.(7), dotted - eq.(8) for $`\tau _c`$.

Parametrization of relaxation time given by eq. (8) leads to small value of variation factor in comparison with other parametrization of relaxation times. Level density rather strongly depends on a relaxation times. It can give an additional possibility for investigation of the relaxation time dependence on temperature and collective state energy.


text:
# Deeply Virtual Compton Scattering at Jefferson Lab, Results and Prospects

## 1 INTRODUCTION

The recently developed formalism of “Generalized Parton Distributions” (GPDs)  showed that information on quark-quark correlations, the transverse quark momentum distribution, and contributions of correlated quark-antiquark pairs (mesons) to the nucleon wave function can be obtained in hard exclusive leptoproduction experiments. GPDs provide a unifying framework for the interpretation of an entire set of fundamental quantities of hadronic structure, such as, the vector and axial vector nucleon form factors, the polarized and unpolarized parton distributions, and the spin components of the nucleon due to orbital excitations.

Deeply Virtual Compton Scattering (DVCS) is one of the key reactions to determine the GPDs experimentally, and it is the simplest process that can be described in terms of GPDs. The dominance of the handbag diagram and the behavior of the reduced forward cross section as $`1/Q^4`$ (scaling regime) is expected to be reached at lower $`Q^2`$ than in the case of deeply exclusive meson production. This is supported by measurements of the $`\gamma ^{}\gamma \pi ^0`$ form-factor in $`e^{}e^+`$ collisions .

One of the first experimental observation of DVCS was obtained from the recent analysis of CLAS data with a $`4.2`$ GeV polarized electron beam in a kinematical regime near $`Q^2`$ = 1.5 GeV2 and $`x_B=0.22`$ . New measurements at higher energies are currently being analyzed, and dedicated experiments are planned. The high luminosity available for these measurements will make it possible to determine details of the $`Q^2`$, $`x_B`$, and $`t`$ dependences of GPDs.

## 2 FIRST OBSERVATION OF EXCLUSIVE DVCS WITH THE CLAS DETECTOR

The DVCS/Bethe-Heitler (BH) interference has recently been measured using the CEBAF Large Acceptance Spectrometer in Hall B at Jefferson Lab . The data were collected as a by-product of the 1999 run with a 4.25 GeV polarized electron beam. At energies above 4 GeV, the CLAS acceptance covers a wide range of kinematics in the deep inelastic scattering domain ($`W`$ 2 GeV and $`Q^2`$ 1 GeV2). The open acceptance of CLAS and the use of a single electron trigger ensures event recording for all possible final states. This experiment measures DVCS via the interference with the Bethe-Heitler (figure 1.

At beam energies accessible at Jefferson Lab, the BH contribution in the cross section is predicted to be several times larger than the DVCS contribution in most regions of the phase space. The dominant BH process can be turned into an advantage by using a longitudinally polarized electron beam: one can measure the helicity-dependent interference term that is proportional to the imaginary part of the DVCS amplitude. In this case the pure real BH contribution is subtracted out in the cross section difference.

For the present DVCS analysis, electron and proton were both detected in the CLAS detector, the reaction $`\stackrel{}{e}pepX`$was studied and the number of single photon final states was extracted by fitting the missing mass ($`M_X^2`$) distributions. The beam spin asymmetry was then calculated as:

$`BSA={\displaystyle \frac{1}{P_e}}{\displaystyle \frac{\left(N_\gamma ^+N_\gamma ^{}\right)}{\left(N_\gamma ^++N_\gamma ^{}\right)}}`$ (1)
Here $`P_e`$ is the beam polarization and, $`N_\gamma ^{+()}`$ is the extracted number of $`\stackrel{}{e}pep\gamma `$events at positive (negative) beam helicity. The resulting $`\varphi `$-dependence is shown in figure 2. A fit to the function

$`F\left(\varphi \right)=A\mathrm{sin}\varphi +B\mathrm{sin}2\varphi `$ (2)
yields $`A=0.217\pm 0.031`$ and $`B=0.027\pm 0.022`$. If the handbag diagram dominates, as expected in the Bjorken regime, $`B`$ should vanish and only the contribution from transverse photons should remain, described by the parameter $`A`$. The GPD analysis including twist-3 contribution shows sensitivity of these data to $`\overline{q}Gq`$ correlations .

## 3 DEDICATED DVCS EXPERIMENTS AT JEFFERSON LAB

There are two dedicated DVCS experiments planned to run using the 6 GeV polarized electron beam. Both experiments plan to detect all three particles in the final state, the scattered electron, the recoil proton, and the photon.

The first experiment E00-110  is a Hall A experiment which is expected to run in 2003. The DVCS beam spin asymmetries and cross section differences will be measured at three $`Q^2`$ intervals, for a fixed interval of $`x_B`$. The experiment will provide a precise check of the $`Q^2`$ dependence of the $`epep\gamma `$ cross section differences (for different beam helicities).

The second experiment E01-113  is a dedicated CLAS DVCS experiment. The main goal of this experiment is to measure the $`t`$, $`\varphi `$, and $`x_B`$ dependence of the beam spin asymmetry for several fixed $`Q^2`$ bins. This quantity is sensitive to the model description of the GPDs. This will be the first time this dependence will be studied with high sensitivity using the DVCS process. A second goal will be to extract the helicity-dependent cross section difference, which directly determines the imaginary part of the DVCS amplitude.

Figure 3 shows examples of the expected CLAS results, the $`t`$ and $`\varphi `$ dependence of the cross section and of the beam spin asymmetry, for three bins in $`Q^2`$ and $`x_B`$. Expected data points are shown only for $`Q^2=2`$ GeV2 and $`x_B=0.35`$. Comparison with different models for the $`\xi `$-dependence of the GPDs is made, together with a first estimation of the twist-3 effects for this process.

The results of these two experiments, Hall A and CLAS, on the $`\stackrel{}{e}pep\gamma `$cross section will allow tests of the $`Q^2`$-dependence to check the scaling behavior. As CLAS covers a broad kinematic range, we will be able to test the $`Q^2`$ dependence for different $`x_B`$. This will verify that we are in a regime where a direct interpretation of the results in terms of GPDs is possible. Observation of significant scaling violations would provide important input for the analysis in terms of higher twist effects.

## 4 CEBAF 12 GeV UPGRADE AND DVCS PROGRAM

The study of GPDs via deeply exclusive reactions is one of the main research programs driving the CLAS upgrade to 12 GeV .

The DVCS cross section and spin asymmetries will be measured in a large range of kinematical bins, simultaneously. Figure 4 illustrates the power of such measurements, showing the expected statistical accuracy of the data, and an example of their binning. In the figure, simulated data on the beam spin asymmetry are shown in the 56 bins of $`x_B`$, $`Q^2`$, and $`t`$ (a total of 1064 points). The outer horizontal scale corresponds to $`x_B`$, divided into 4 bins, shown by solid lines. The outer vertical scale represents $`Q^2`$. The $`Q^2`$ bins shown as rows of asymmetry distributions. Larger bins at higher $`Q^2`$ are to compensate for the fast drop of the cross section. In each row, several beam spin asymmetry distributions as a function of $`\varphi _{\gamma \gamma }`$ (running from $`0^{}`$ to $`360^{}`$) correspond to different $`t`$ bins. The average value of $`t`$ is shown in the upper corner of each plot.

Figure 5 shows more clearly the accuracy of such a measurement and its sensitivity to the models in 2 selected bins, $`x_B=0.35`$ and $`t=0.3`$ (GeV/c)2, for $`Q^2=2.75`$ (GeV/c)2 in figure 5.a, and for $`Q^2=5.4`$ (GeV/c)2 in figure 5.b. As an input to the simulation cross sections calculated in Ref. based on Ref.  have been used. Different curves in the figure correspond to different model assumptions. One sees sufficient sensitivity for the separation of models at both $`Q^2`$.

## 5 SUMMARY

A first measurement of the beam spin asymmetry in the exclusive electroproduction of real photons in the deep inelastic regime was presented. We see a clear asymmetry, as expected from the interference of the DVCS and BH processes. It has been shown that our results can be accomodated within a GPD analysis . This supports the expectations that DVCS will allow access to GPDs at relatively low energies and momentum transfers. This opens up a new avenue for the study of nucleon structure which is inaccessible in inclusive scattering experiments. Dedicated DVCS experiments at 6 GeV electron beam energy are planned, which will allow significant expansion of the $`Q^2`$ and $`x_B`$ range covered in these studies. The high luminosity available for these measurements will make it possible to map out details of the $`Q^2`$, $`x_B`$, and $`t`$ dependences of GPDs.

The CEBAF 12 GeV upgrade will allow a breakthrough program, in the study of nucleon structure via deeply virtual exclusive reactions, to be carried out with unprecedented precision.


text:
# Implications of Variability Patterns observed in TeV Blazars on the Structure of the Inner Jet

## 1 INTRODUCTION

Blazars are active galactic nuclei exhibiting the most rapid and largest amplitude variability of all AGN. Historically, radio observations first revealed that the emission was luminous and rapidly variable. With that, assuming that the radio emission was due to synchrotron radiation, calculations of Compton up-scattering predicted much higher X–ray fluxes than the observed values unless the radio emission was relativistically beamed (Hoyle et al., 1966; Jones, O’Dell, & Stein, 1974). This led to our current model for blazars where the entire electromagnetic emission arises in a relativistic jet pointing close to the line of sight (Blandford & Rees, 1978). Subsequent radio observations using the Very Long Baseline Interferometry (VLBI) showed superluminal motion in many sources, which served as the direct evidence for the relativistic motion.

The broadband spectra of blazars consist of two peaks, one in the radio to optical–UV range (and in some cases, reaching to the X–ray band), and the other in the X–ray to $`\gamma `$–ray region. From the high polarization of the radio to optical emission, the lower energy peak is best interpreted as produced via the synchrotron process by relativistic electrons in the jet. The higher energy peak is believed to be due to Compton up-scattering by the same population of relativistic electrons. Several possibilities exist for the source of the seed photons; these can be the synchrotron photons internal to the jet (Jones, O’Dell, & Stein, 1974; Ghisellini & Maraschi, 1989), but also external, such as from the broad emission line clouds (Sikora, Begelman, & Rees, 1994) or from the accretion disk (Dermer, Schlickeiser, & Mastichiadis, 1992; Dermer & Schlickeiser, 1993).

Blazars are commonly detected as $`\gamma `$–ray sources. A number of them with peak synchrotron output in the X–ray range also have been detected in the TeV range with ground-based Cherenkov arrays. These are the so-called “TeV blazars.” In TeV blazars, X–rays provide the best means for studying variability properties: this is because X-ray flux is presumably produced by electrons that are accelerated to the highest energies (where the cooling time scales are most rapid), and thus the dilution by the non-varying components is the smallest.

Variability studies of blazars have entered a new stage after a number of continuous long-look X–ray observations conducted with the ASCA satellite. One such observation, of Mrk 421, conducted in 1998, showed for the first time that flaring is actually occurring on a daily basis, and that long-duration flares detected in previous observations were probably unresolved superpositions of multiple, more rapid flares (Takahashi et al., 2000). Such excellent data provided the new knowledge about the radiation processes, allowing an exploration of the actual dynamics of the particles that are accelerated in the jets.

In §2, we present the properties of X–ray variability observed from TeV blazars, and summarize the issues we need to explain. In §3, we consider the internal shock scenario, which involves shells propagating rapidly along the jet. Using this scenario, we simulate expected X-ray light curves, and study how various quantities that can be measured from observations depend on the input parameters to the model. In §4, we compare the results of simulations to the data to determine whether the observed variability features can be reproduced, and if so, what are the implications on model parameters. We briefly consider the possibility of the external shock scenario in §5, and give a summary in §6.

## 2 VARIABILITY PROPERTIES OBSERVED IN TEV BLAZARS

One of the most surprising results from the long look observations of TeV blazars was the repeated occurrence of flares with a time scale of $``$ one day. This was first observed in the 7–day observation of Mrk 421 in 1998 (Takahashi et al., 2000), and was confirmed via the 10–day observations of both Mrk 501 and PKS 2155–304 in 2000 (Tanihata et al., 2001). The latter two sources were in a relatively low flux state compared to previous observations, which implied that a high state is not a requirement for rapid variability. In other words, this feature indicates that the time scale of the rise and decay of the flares are similar to the time scale of the repeating of the flares.

Another observational fact is that the X–ray flares always appear to lie on top of an underlying offset-like component (see, e.g., light curves in Urry et al. 1997; Tanihata et al. 2001; Zhang et al. 2002). At this point, we cannot distinguish whether this is due to superpositions (or pile-ups) of flares, or whether there is a separate, non-varying offset component, but the observations indicate that there is always some component such that a flare does not start from zero, but from a level comparable to the flare amplitude.

An important advantage of long continuous observations is that the variability can be treated statistically. Recent structure function analysis have showed that there is clearly a characteristic time scale $`t_{\mathrm{chr}}`$ of an order of a day for each of three TeV blazars, Mrk 421, Mrk 501, and PKS 2155–304, and also that the variability power at shorter time scales than $`t_{\mathrm{chr}}`$ is strongly suppressed (Kataoka et al., 2001; Tanihata et al., 2001). Together with the near-symmetry of the flares and also by comparing the various possible time scales, it was suggested that this $`t_{\mathrm{chr}}`$ is determined by some dynamical time scale, instead of energy dependent cooling or acceleration times (Kataoka, 2000).

X–ray emission from TeV blazars shows spectral variability, where the common trend indicates a harder X–ray spectrum for higher intensity states. Furthermore, the improved spectral coverage allows a measurement of the exact location of the synchrotron peak. Perhaps the most striking case is Mrk 501, where the peak frequency shifted from below 1 keV to over 100 keV during the very high flux state (Pian et al., 1998; Tavecchio et al., 2001). Likewise, Fossati et al. (2000b) has shown using the BeppoSAX observations of Mrk 421 in 1997 and 1998 that the peak frequency shifted to the higher energy for higher intensity states; this is particularly apparent during the flare observed in 1998 April. This was also shown in the long look observation of Mrk 421 in 1998 using the combined ASCA and EUVE spectrum (Tanihata et al., 2002). Through a detailed analysis of the spectral evolution at the rise of the flare, we found several flares that start to appear from the higher energy, which strongly suggests that an appearance of a new harder (i.e., higher synchrotron peak frequency) emission component is associated with the generation of the flare.

Summarizing, the following features are the issues concerning the variability derived from observations, that we need to explain in considering any model: (1) Daily-flares (i.e. $`t_{\mathrm{cycle}}t_{\mathrm{chr}}`$), (2) the offset component, (3) the structure function, and (4) the energy dependence.

## 3 THE INTERNAL SHOCK SCENARIO

Among many studies addressing the mechanism of particle acceleration in jets, acceleration via shocks appears to be the most viable. Such shocks can efficiently accelerate particles to very high energies (e.g. Longair 1994; Bell 1978; Drury 1983; Blandford & Eichler 1987; Jones & Ellison 1991), and since it is highly probable that shocks form inside jets, we consider it in more detail here.

In order to form a shock, there must be a large velocity difference between the colliding parcels of matter. The key idea of the internal shock model is a shock-in-jet scenario, where the central engine injects energy into the jet in a discontinuous manner, producing individual shells having slightly different bulk Lorentz factors and energies. If this occurs, there will be collisions by a faster shell catching up to a slower one, forming a shock. This internal shock scenario is among the most promising models to explain the emission of gamma–ray bursts (e.g. Sari & Piran 1995), although it was originally suggested in reference to AGN more than 20 years ago by Rees (1978). It has been recently suggested that this model could be successfully applied to blazars, as it can explain some of their basic properties such as the low efficiency (Ghisellini, 2001; Spada et al., 2001; Sikora, Błażejowski, Begelman, & Moderski, 2001).

### 3.1 The Model

In order to investigate whether the observed variability properties can be reproduced in the internal shock scenario, we developed a simulation code. Recently, Spada et al. (2001) has given a detailed calculation covering the formation, propagation, and collision of such shells. Included are hydrodynamic calculations to determine the structure of the shock fluid, and the full radiation spectrum is derived by summing up all the locally produced spectra from the electrons accelerated by the shocks. The model and simulations presented here are a much simplified version, to be compared specifically against the actual measurements of TeV blazars.

Due to the difference of the initial velocity of the shells ejected from the base of the jet, collisions occur when a faster shell catches up to a slower one. This is where the shock is formed, electrons are accelerated, and lose energy through radiation. Since we consider only the colliding shells, in our simulations we generate pairs of shells with one having a bulk Lorentz factor (BLF) of $`\mathrm{\Gamma }_1`$ ejected from the base of the jet, and a following shell with a BLF of $`\mathrm{\Gamma }_2`$ ejected after an initial separation distance of $`D_0`$. If $`\mathrm{\Gamma }_2`$ is larger than $`\mathrm{\Gamma }_1`$, the latter shell will catch up to the former one, and the two will merge into a single shell producing a shock. Each collision then generates radiation (called hereafter a “shot”) for a duration of $`t_{\mathrm{shot}}`$. The time profile of each shot is assumed to have a symmetric linear rise and decay. The collisions are distributed randomly in time following a Poisson distribution, and superposition of these individual shots results in the output light curve.

For simplicity, the rest mass of all shells is assumed to be the same, and the shell thickness $`l`$ is assumed to be equal to the initial separation $`D_0`$. The average frequency of the collision is set to be consistent with the separation of the shells (i.e. $`F_{\mathrm{col}}=c/2(D_0+l)`$), and only the first collision is considered. The initial BLFs are assumed to be distributed around an average value $`\mathrm{\Gamma }_{\mathrm{avrg}}`$, following a Gaussian distribution with its width described by the sigma, $`\sigma _\mathrm{\Gamma }`$. In this case, there are only three input parameters: two of them describe the distribution of the initial BLFs of the ejected shells ($`\mathrm{\Gamma }_{\mathrm{avrg}}`$, $`\sigma _\mathrm{\Gamma }`$), and one describes the initial separation of the two colliding shells ($`D_0`$). We also define the fractional width $`\sigma _\mathrm{\Gamma }^{}`$ $`\sigma _\mathrm{\Gamma }/\mathrm{\Gamma }_{\mathrm{avrg}}`$.

With the assumptions described above, using the momentum and energy conservation laws in an inelastic collision, the BLF of the merged shell is

$$\mathrm{\Gamma }_\mathrm{m}=(\mathrm{\Gamma }_1\mathrm{\Gamma }_2)^{\frac{1}{2}}.$$
(1)
The newly generated total internal energy is

$$E_\mathrm{m}=Mc^2(\mathrm{\Gamma }_1\mathrm{\Gamma }_\mathrm{m})+Mc^2(\mathrm{\Gamma }_2\mathrm{\Gamma }_m),$$
(2)
and thus the dynamical efficiency is given by

$$\eta =\frac{Mc^2(\mathrm{\Gamma }_1\mathrm{\Gamma }_\mathrm{m})+Mc^2(\mathrm{\Gamma }_2\mathrm{\Gamma }_m)}{Mc^2\mathrm{\Gamma }_1+Mc^2\mathrm{\Gamma }_2}.$$
(3)
We assume that all of the newly generated internal energy is converted to the random energy of the electrons.

Each collision will take place at distance

$$D\frac{2\mathrm{\Gamma }_1^2\mathrm{\Gamma }_2^2}{\mathrm{\Gamma }_2^2\mathrm{\Gamma }_1^2}D_0$$
(4)
from the core, where in making the above approximation, we used $`\mathrm{\Gamma }_1^2`$, $`\mathrm{\Gamma }_2^21`$. We assume that the jet is collimated into a cone with an opening angle $`\theta 1/\mathrm{\Gamma }`$. The radius of the shell at the location of the collision can thus be written as

$$R=D\mathrm{tan}\theta \frac{D}{\mathrm{\Gamma }}.$$
(5)
The time duration of each shot $`t_{\mathrm{shot}}`$ results from several competing time scales: the acceleration and cooling time scales, and the dynamic time scale, which includes the hydrodynamic time and the angular spreading time. It has been shown for TeV blazars that the cooling or acceleration times are significantly shorter than the observed rise and decay time scales of the X–ray and gamma–ray flares (Kataoka, 2000; Tanihata et al., 2001), and thus the dynamical time scale most likely determines $`t_{\mathrm{shot}}`$, and also filters out any faster variability.

The hydrodynamic time scale is determined by the time that the shock takes to cross the shell. Using the conservation of mass, energy, and momentum at the shock, and the equality of pressure and velocity along the contact discontinuity, the Lorentz factors of the forward and reverse shocks $`\mathrm{\Gamma }_{\mathrm{fs}}`$ and $`\mathrm{\Gamma }_{\mathrm{rs}}`$ are given by (Kobayashi, Piran, & Sari, 1997)

$`\mathrm{\Gamma }_{\mathrm{fs}}`$ $``$ $`\mathrm{\Gamma }_\mathrm{m}\sqrt{(1+{\displaystyle \frac{2\mathrm{\Gamma }_\mathrm{m}}{\mathrm{\Gamma }_1}})/(2+{\displaystyle \frac{\mathrm{\Gamma }_\mathrm{m}}{\mathrm{\Gamma }_1}})}`$ (6)
$`\mathrm{\Gamma }_{\mathrm{rs}}`$ $``$ $`\mathrm{\Gamma }_\mathrm{m}\sqrt{(1+{\displaystyle \frac{2\mathrm{\Gamma }_\mathrm{m}}{\mathrm{\Gamma }_2}})/(2+{\displaystyle \frac{\mathrm{\Gamma }_\mathrm{m}}{\mathrm{\Gamma }_2}})}.`$ (7)
We note that Equations (6) and (7) are for the case of relativistic shocks, where the adiabatic index is 4/3. For non-relativistic shocks, the adiabatic index should be 5/3, but since this does not make a large difference, we use this formula for all shocks in this paper. Following Kobayashi, Piran, & Sari (1997), we estimate the shock crossing time $`t_{\mathrm{crs}}`$ by the longer time scale of the two shocks to cross the shell; this happens to be the time which the reverse shock takes to cross the faster shell. Because the emitting region is moving with a relativistic speed towards the observer, the observed duration of the flare will be shortened by a factor of $`(1\beta \mathrm{cos}\theta )1/\mathrm{\Gamma }_\mathrm{m}^2`$, and thus,

$$t_{\mathrm{crs}}=\frac{l}{c(\beta _2\beta _{\mathrm{rs}})\mathrm{\Gamma }_\mathrm{m}^2}\frac{2}{c\mathrm{\Gamma }_\mathrm{m}^2}\left(\frac{1}{\mathrm{\Gamma }_{\mathrm{rs}}^2}\frac{1}{\mathrm{\Gamma }_2^2}\right)^1D_0$$
(8)
where $`\beta _2c`$ and $`\beta _{\mathrm{rs}}c`$ are respectively the velocity of the catching up shell and the reverse shock. The angular spreading time is given by

$$t_{\mathrm{ang}}=\frac{R}{c\mathrm{\Gamma }_\mathrm{m}}.$$
(9)
The average time cycle of the collisions is determined from the frequency of the ejection of the shells, and thus

$$t_{\mathrm{col}}=\frac{4}{c}D_0.$$
(10)
For the $`t_{\mathrm{shot}}`$ in the simulations, we use the longer one of the $`t_{\mathrm{crs}}`$ and $`t_{\mathrm{ang}}`$ for each collision. As a result, since $`t_{\mathrm{crs}}`$ is always longer than $`t_{\mathrm{ang}}`$ for all cases, thus $`t_{\mathrm{shot}}`$=$`t_{\mathrm{crs}}`$.

### 3.2 Results and the Dependence on Each Parameter

We first consider the case of $`\mathrm{\Gamma }_{\mathrm{avrg}}`$=10, $`\sigma _\mathrm{\Gamma }^{}`$=0.05, and $`D_0`$=$`3\times 10^{13}`$ cm. The distribution of the collision distances in this case is shown in Figure 1(a). Our simulations show that the collisions take place at distances ranging from $`D10^{17}`$ cm up to $`D10^{20}`$ cm. Figure 1(b) shows the amount of the newly generated internal energy for each collision $`E_\mathrm{m}`$ plotted against the collision distance $`D`$, and Figure 1(c) is the time scale of each generated flare $`t_{\mathrm{shot}}`$ plotted against $`D`$. It is apparent that $`E_\mathrm{m}`$ is larger and $`t_{\mathrm{shot}}`$ is shorter for collisions which occur at smaller $`D`$.

A portion of the simulated light curve for this case is shown in Figure 1(d). It can be seen that the overall light curve is characterized by repeating flares having time scales of $`50100`$ ks, resembling the light curves observed from TeV blazars. This results from the fact that the shells which collide at smaller distances have larger $`E_\mathrm{m}`$ and shorter $`t_{\mathrm{shot}}`$; with this, the amplitude of the emission becomes much larger compared to the shots generated from collisions at larger distances. Accordingly, only the shots produced by collisions at the smallest distances will be apparent as flares in the observed light curve. The average frequency of the collision in this case is $`F_{\mathrm{col}}=0.24`$ mHz, corresponding to one collision per $``$4 ks on average, while the number of flares which are observed is only about 6% of the total number of shots. The distribution of the resolved flares are shown as the shadowed area in Figure 1(a).

The calculated structure function for the simulated light curve is shown in Figure 1(e). A clear break, indicating a characteristic time scale $`t_{\mathrm{chr}}`$ is seen, and the slope at the shorter time scales is $`\beta `$2. This indicates that the values of $`t_{\mathrm{shot}}`$ of observable flares are restricted to a rather narrow range, which determines $`t_{\mathrm{chr}}`$, and that there is very little variability power below $`t_{\mathrm{chr}}`$ – exactly what we have observed from TeV blazars.

Another remark regards the offset-like component. This is due to the other 94% of the collisions that generate the longer, smaller amplitude shots, which overlap each other, resulting in the observed offset. In the following, we show how each of the parameters affects on these results.

#### 3.2.1 Dependence on $`\mathrm{\Gamma }_{\mathrm{avrg}}`$

We first consider the effect of the change in the average BLF $`\mathrm{\Gamma }_{\mathrm{avrg}}`$, while the relative width $`\sigma _\mathrm{\Gamma }^{}`$, and $`D_0`$ are kept constant. When $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ increases by a factor of $`a`$, Equation 4 implies that the collisions will occur at a distance $`D`$ which is greater, by a factor of $`a^2`$. On the other hand, Equation 8 shows that the the duration of each shot will be conserved, and thus only the amplitude will be different in the resulting light curve.

#### 3.2.2 Dependence on $`D_0`$

We then consider the case where the initial separation of the colliding shells $`D_0`$ changes, while $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ and $`\sigma _\mathrm{\Gamma }`$ are kept constant. When $`D_0`$ is larger by a factor of $`a`$, Equations 8 and 10 imply that both $`t_{\mathrm{crs}}`$(=$`t_{\mathrm{shot}}`$) and $`t_{\mathrm{col}}`$ will be longer by a factor of $`a`$, and thus the simulated light curve will simply be a stretched-in-time version of the original one – the amount of relative offset components will be identical. The collision distance $`D`$ will become larger by a factor of $`a`$.

#### 3.2.3 Dependence on $`\sigma _\mathrm{\Gamma }^{}`$

This turns out to be the most important parameter. Here, we fix $`\mathrm{\Gamma }_{\mathrm{avrg}}`$=10 and $`D_0`$=$`3\times 10^{13}`$ cm, and simulate light curves with different $`\sigma _\mathrm{\Gamma }^{}`$. First, it can be shown from Equation 4 that as $`\sigma _\mathrm{\Gamma }^{}`$ becomes larger, the collisions start to take place at shorter distances from the core.

Since the characteristic time scale $`t_{\mathrm{chr}}`$ is always determined by the shots due to collisions which took place at the smallest distances, the time scales of the observable flares in the light curve become shorter when $`\sigma _\mathrm{\Gamma }^{}`$ is larger. This is shown in the simulated light curve for different $`\sigma _\mathrm{\Gamma }^{}`$ ($`\sigma _\mathrm{\Gamma }^{}`$=0.001, 0.005, and 0.05) in Figure 2. Here, the flares appear more spiky as $`\sigma _\mathrm{\Gamma }^{}`$ becomes larger. Note that since all 3 simulations assume the same $`D_0`$, the number of collisions per unit time is the same for all 3 light curves. On the other hand, the number of visible flares in the light curves is clearly different. In order to see the difference in $`t_{\mathrm{chr}}`$, we calculated the structure function for the simulated light curves. This is shown in Figure 3. The break, indicating $`t_{\mathrm{chr}}`$, is clearly seen to shift to the shorter time scales as $`\sigma _\mathrm{\Gamma }^{}`$ becomes larger.

The differences in the spikiness in the light curve can be regarded as the differences in the relative amplitude of the flare and offset components. Importantly, what we have shown is that this relative amplitude, $`r_{\mathrm{fo}}`$, changes with the value of $`\sigma _\mathrm{\Gamma }^{}`$. This indicates that, assuming that all of the offset component is generated by the internal shocks, we can compare the observed $`r_{\mathrm{fo}}`$ to the that derived through simulations. Indeed, the actual observations have clearly shown a larger offset component than, for instance is apparent in Figure 2(c), which suggests that $`\sigma _\mathrm{\Gamma }^{}`$ must be relatively small. We will quantify this in section 4.

#### 3.2.4 Effect of Fluctuation in $`l`$

We have so far fixed the shell thickness $`l`$ to be equal to the initial separation of the two shells $`D_0`$. Here at the end of this section, we consider the effect when $`l`$ also fluctuates. We assume that $`l`$ is distributed around an average $`D_0`$, following a Gaussian function with its width described by the relative sigma, $`\sigma _l^{}`$ ($`\sigma _l/D_0`$). We also make an assumption that the mass of the ejected shells scales with the shell thickness (i.e. density of shell is constant).

The effect of the fluctuation is demonstrated in the plotted correlation of the duration of the radiation ($`t_{\mathrm{shot}}`$) and the dissipated energy ($`E_\mathrm{m}`$) for each collision in Figure 4. For the case where $`l`$ is fixed (a), $`t_{\mathrm{shot}}`$ and $`E_\mathrm{m}`$ become simply respectively shorter and larger as $`\mathrm{\Delta }\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }_2\mathrm{\Gamma }_1`$) of the colliding shells increases, indicating the clear trend in the correlation. When $`l`$ fluctuates, the spread becomes wider, as shown in Figures 4(b)(c).

This indicates that there will be increasing power in the faster variability time scales when $`\sigma _l^{}`$ becomes larger. For a shell with smaller $`l`$, $`E_\mathrm{m}`$ will be reduced because of the smaller mass of the shell. However, since the shock crossing time $`t_{\mathrm{crs}}`$ (=$`t_{\mathrm{shot}}`$) also becomes shorter, the amplitude of the shot will be as large as that of a longer shot due to a larger $`l`$. This is demonstrated in the calculated structure function from the light curves simulated for different $`\sigma _l^{}`$, shown in Figure 5. While the break is seen to stay nearly at the same value, the slope at the faster time scale appears to flatten as $`\sigma _l^{}`$ increases.

### 3.3 Summary of the Simulations

The simulation results can be summarized as follows:

1. Collisions of two shells which had the largest relative velocity, and accordingly collided at the shortest distances ($`D`$), are the shots which appear as the strongest observable flares in the light curve.
2. These shots determine the characteristic time scale ($`t_{\mathrm{chr}}`$) of the variability. There is very little variability power on time scales shorter than this $`t_{\mathrm{chr}}`$.
3. An offset component will arise from emission due to overlapping shots produced by collisions at larger distances. The relative amplitude of flare to offset ($`r_{\mathrm{fo}}`$) is a function of the initial width of the bulk Lorentz factor ($`\sigma _\mathrm{\Gamma }^{}`$).
4. The dependences of each parameter are:

* $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ determines $`E_\mathrm{m}`$: higher $`E_\mathrm{m}`$ yields higher $`\mathrm{\Gamma }_{\mathrm{avrg}}`$.
    * $`D_0`$ determines the normalization of the time series: longer $`t_{\mathrm{chr}}`$ is a result of larger $`D_0`$.
    * $`\sigma _\mathrm{\Gamma }^{}`$ has an effect on $`t_{\mathrm{chr}}`$ and $`r_{\mathrm{fo}}`$: larger $`\sigma _\mathrm{\Gamma }^{}`$ results in smaller $`t_{\mathrm{chr}}`$ and larger $`r_{\mathrm{fo}}`$.
    * $`\sigma _l^{}`$ has effect on the slope $`\beta `$ of the structure function: $`\beta `$ flattens as $`\sigma _l^{}`$ increases.

## 4 APPLICATION OF THE INTERNAL SHOCK MODEL TO THE X-RAY DATA FOR TEV BLAZARS

### 4.1 Light Curves

In the previous section, we have shown that the internal shock scenario naturally predicts the main features of blazar light curves described in §2, under the condition of $`\sigma _\mathrm{\Gamma }^{}0.01`$. One prediction is that the typical time scale of the observed flares ($`t_{\mathrm{chr}}`$) always becomes similar to the time scale of the flare $`cycles`$. This is very much consistent with the actual observations, where the “day-scale” flares are observed “daily.”

The next prediction concerns the features observed in the structure function analysis. The observations show that all TeV blazars show a break in the structure function, with $`t_{\mathrm{chr}}1`$ day. The slope below this $`t_{\mathrm{chr}}`$ is steep ($`>`$1), suggesting that very little variability power is below this $`t_{\mathrm{chr}}`$. We have shown that the structure functions calculated from the simulated light curves show the same features. In fact, this provides an explanation for the non-existence of shorter time scale variability – there are no collisions until a certain distance $`D`$ from the central core.

In order to see the actual $`\sigma _\mathrm{\Gamma }^{}`$ dependence on $`t_{\mathrm{chr}}`$, we simulated light curves for a series of $`\sigma _\mathrm{\Gamma }^{}`$ for the case of $`\mathrm{\Gamma }_{\mathrm{avrg}}`$=10, and $`D_0`$=$`3\times 10^{13}`$ cm as in §3.2.3. We calculated the structure function for each simulated light curve, and assumed that it is well described by two power laws. We then fitted for the two slopes, and estimated the characteristic time scale as the point where the two slopes cross (as shown as the dotted lines in Figure 3). The derived values of $`t_{\mathrm{chr}}`$ plotted as a function of $`\sigma _\mathrm{\Gamma }^{}`$ are shown in Figure 6. It appears that $`t_{\mathrm{chr}}`$ becomes shorter with larger $`\sigma _\mathrm{\Gamma }^{}`$, and the best-fit power-law gives a slope of $`0.9`$. Note that we have shown in §3.2.1 that the time scales do not depend on $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ if $`\sigma _\mathrm{\Gamma }^{}`$ is constant.

Another point regards the presence of the offset component. As discussed in §2, the observed flares in TeV blazars always appear to lie on top of an underlying offset component. The data do not tell us directly whether this offset is due to flares overlapping each other, or due to some steady emission component. What we have shown here is that considering the internal shock scenario, there will always be many overlapping flares which will appear as the offset component.

The amount of this offset component should also be useful in modeling of an actual observation, as this suggests that $`\sigma _\mathrm{\Gamma }^{}`$ can be estimated from an observed light curve if it has a sufficient length to estimate the relative amplitude of the flare and offset component. In an attempt to quantify the offset component, here we use the parameter $`r_{\mathrm{fo}}`$, describing the relative amplitude of the flare and offset component as follows. We first generate a histogram of the count rates which form a peak. Since the lower and higher end of the peak represent the minimum and maximum count rates in the light curve, these two can be considered as an indicator of the amplitude of the offset component, and the offset-plus-flare component. As there are fluctuations, we define the offset amplitude as the point where 10% of total counts is reached, and the offset-plus-flare amplitude as the point where 90% of the total counts is included. The amplitude of the flare component is estimated by subtracting the offset amplitude from the offset-plus-flare amplitude.

We calculate this flare-to-offset ratio $`r_{\mathrm{fo}}`$ for the same set of simulated light curves used in Figure 6. The result is shown in Figure 7, which shows that $`r_{\mathrm{fo}}`$ increases with increasing $`\sigma _\mathrm{\Gamma }^{}`$. As we showed in §3, $`r_{\mathrm{fo}}`$ does not depend on $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ nor $`D_0`$, which means that $`\sigma _\mathrm{\Gamma }^{}`$ can be directly estimated from a $`r_{\mathrm{fo}}`$ given by an observed light curve.

Finally, we note that the $`t_{\mathrm{chr}}`$ plotted in the Figure 6 is for the case of $`D_0`$=$`3\times 10^{13}`$ cm. Given that the time series scales linearly with $`D_0`$, the vertical axis in Figure 6 can be regarded as $`t_{\mathrm{chr}}/D_{0,3\times 10^{13}}`$, where $`D_{0,3\times 10^{13}}`$ denotes $`D_0`$ in units of $`3\times 10^{13}`$ cm. Accordingly, this indicates that $`D_0`$ can also be estimated if $`r_{\mathrm{fo}}`$ and $`t_{\mathrm{chr}}`$ can be measured. This is interesting, given that the initial separation is a value determined from the frequency of the emitted shells. Thus, $`D_0`$ should reflect the frequency of the activity of the central engine, that generates and ejects individual shells.

### 4.2 Variation of the Synchrotron Peak Frequency

The important observational result from the spectral analysis was that the peak synchrotron frequency increases during both high intensity states, and also within the daily flares. In this section, we discuss whether this can be interpreted within the same model. We start with formulating the peak energy, by generally following the prescriptions of Inoue & Takahara (1996) and Kirk, Rieger & Mastichiadis (1998).

The peak frequency of the synchrotron spectrum reflects the maximum energy of the accelerated electrons, which is determined from where the cooling and the acceleration time scale becomes comparable. The cooling time of the electrons due to synchrotron and inverse Compton emission can be written as (Rybicki & Lightman, 1979),

$$\tau _{\mathrm{cool}}(\gamma )=\frac{3m_\mathrm{e}c}{4(U_\mathrm{B}+U_{\mathrm{soft}})\sigma _\mathrm{T}\gamma },$$
(11)
where $`U_\mathrm{B}`$ and $`U_{\mathrm{soft}}`$ are the energy densities of the magnetic field and the soft photons to be upscattered in the Inverse Compton process, $`m_\mathrm{e}`$ the rest mass and and $`\gamma `$ is the random Lorentz factor of the electron, and $`\sigma _\mathrm{T}`$ is the Thomson cross section.

The acceleration time is not as well understood as the cooling time. Perhaps the most promising theory suggests that it is determined from the time scale of the 1st order Fermi acceleration process operating in a shock. In this case, the acceleration time can be approximated by considering the mean free path $`\lambda (\gamma )`$ for the scattering of the electrons with the magnetic disturbances. Taking the mean free path to be proportional to the Larmor radius by introducing another parameter $`\xi `$ ($`\lambda \xi \frac{\gamma mc^2}{eB}`$), the acceleration time is given by,

$$\tau _{\mathrm{acc}}(\gamma )=\frac{20\lambda (\gamma )c}{3v_\mathrm{s}^2}=\frac{20m_\mathrm{e}c^3}{3eB}\frac{\xi }{v_\mathrm{s}^2}\gamma ,$$
(12)
where $`v_\mathrm{s}`$ is the shock velocity. By equating the radiative cooling and acceleration times in Equations 11 and 12, the maximum energy of the electrons is,

$$\gamma _{\mathrm{max}}=\frac{v_\mathrm{s}}{c}\left(\frac{9eB}{80(U_\mathrm{B}+U_{\mathrm{soft}})\sigma _\mathrm{T}\xi }\right)^{\frac{1}{2}}.$$
(13)
Furthermore, in the case of TeV blazars, the Compton cooling for the energy range of electron emitting X–rays is strongly suppressed by the Klein-Nishina cutoff, and thus $`U_{\mathrm{soft}}/U_\mathrm{B}1`$ (e.g. Kataoka 2000; Li & Kusunose 2000). Thus for this case, the observed maximum synchrotron frequency can be written as

$`\nu _{\mathrm{max}}`$ $`=`$ $`1.2\times 10^6B\delta \gamma _{\mathrm{max}}^2`$ (14)
$``$ $`2.49\times 10^{21}{\displaystyle \frac{\delta _{10}}{\xi }}\left({\displaystyle \frac{v_\mathrm{s}}{c}}\right)^2\mathrm{Hz},`$ (15)
where $`\delta _{10}`$ denotes the beaming factor in units of 10. Note that the above implies that when synchrotron losses dominate (as is likely to be the case for TeV blazars), the observed peak frequency is independent on $`B`$.

We will first consider the differences in the synchrotron spectrum for the flare components and the underlying offset component. Above, we have shown that in the framework of the internal shock scenario, the observed flares can be regarded as being due to collisions which had the largest $`\mathrm{\Delta }\mathrm{\Gamma }`$ within the distribution, while the offset component results from emission due to fainter, overlapping flares produced by collisions that had smaller $`\mathrm{\Delta }\mathrm{\Gamma }`$.

Since the collisions which generate the observed flares occur at shorter distances $`D`$, the magnetic field $`B`$ is expected to be higher. Accordingly, the radiative cooling times of the electrons are shorter there, preventing electrons from being accelerated to higher energies. On the other hand, Equation 12 suggests that the acceleration time will also be shorter, due to the larger $`B`$, and also due to larger shock velocity $`v_\mathrm{s}`$. Equation 15 shows that concerning the observed synchrotron peak frequency, the changes in $`B`$ cancel out and only values of $`v_\mathrm{s}`$, $`\delta `$, and $`\xi `$ affect $`\nu _{\mathrm{max}}`$.

In Figure 8(a), we show the calculated shock velocity $`v_\mathrm{s}`$ for each collision (in units of $`c`$; $`\beta _\mathrm{s}`$=$`v_\mathrm{s}/c`$) as a function of $`D`$, for the case of $`\mathrm{\Gamma }_{\mathrm{avrg}}`$=10, $`\sigma _\mathrm{\Gamma }^{}`$=0.05, and $`D_0`$=$`3\times 10^{13}`$ cm (similar to the ones used in Figure 1). It shows that $`v_\mathrm{s}`$ decreases as $`D`$ increases, following a relation of $`\beta _\mathrm{s}D^1`$. Accordingly, if $`\xi `$ were to be similar for all collisions, Equation 15 shows that $`\nu _{\mathrm{max}}`$ will be higher for collisions at smaller $`D`$. Since these are the shots that appear as flares, this is consistent with the actual observations, where the flare spectrum components show higher $`\nu _{\mathrm{max}}`$ than the offset component. The calculated $`\nu _{\mathrm{max}}`$, considering the Bohm limit (i.e. $`\xi `$=1) is shown in Figure 8(b). The dependence of $`\nu _{\mathrm{max}}`$ on $`D`$ is $`\nu _{\mathrm{max}}D^2`$.

Finally, in addition to the spectral variation during the day-scale flares, the same trend regarding the relationship between the flux and synchrotron peak frequency is also observed in comparing different observations. The same discussions as above will hold for the case when $`\sigma _\mathrm{\Gamma }^{}`$ becomes larger. The shock velocity will become higher for the collisions which generate the observed flares, and accordingly the average $`\nu _{\mathrm{max}}`$ will be higher. Another way to generate a change in the average $`\nu _{\mathrm{max}}`$ is the change in the $`\mathrm{\Gamma }_{\mathrm{avrg}}`$. If $`\mathrm{\Gamma }_{\mathrm{avrg}}`$ increases, it can be shown from equation 15 that $`\nu _{\mathrm{max}}`$ will also increase.

### 4.3 The Case of Mrk 421

Here we consider the observed light curve of Mrk 421 during the ASCA long look in 1998 April. In modeling the observed light curve, we wish to reproduce the observables determined from these data such as the characteristic time scale $`t_{\mathrm{chr}}`$ and the flare-to-offset ratio $`r_{\mathrm{fo}}`$ (as defined in § 4.1).

The first step is to estimate $`\sigma _\mathrm{\Gamma }^{}`$ from the observed $`r_{\mathrm{fo}}`$. The normalized light curve (in counts s<sup>-1</sup>) and the calculated histograms of the count rate are shown in Figure 9. The calculated $`r_{\mathrm{fo}}`$=0.7 is compared with the derived relation of $`\sigma _\mathrm{\Gamma }^{}`$ and $`r_{\mathrm{fo}}`$ from the simulations shown in Figure 7, which gives an estimate of $`\sigma _\mathrm{\Gamma }^{}`$0.001. Then, from the observed characteristic time scale $`t_{\mathrm{chr}}40`$ ks (Tanihata et al., 2001), $`D_0`$ can be estimated from Figure 6 as

$$D_0\frac{40}{120}\times 3\times 10^{13}1\times 10^{13}\mathrm{cm}.$$
(16)
From Equation 10, this would give $`t_{\mathrm{col}}1300`$ seconds, which indicates that the shells are ejected from the central engine on average every $``$10 minutes.

On the other hand, we also remark here that the $`r_{\mathrm{fo}}`$ derived in Figure 7 is for the variations of the dissipated energy. Concerning the $`r_{\mathrm{fo}}`$ for a particular energy, this will have an energy dependence, which requires a proper treatment of synchrotron formulae, and will also depend on the energy dependence of the efficiency of the energy that is to be converted to radiation. This might be further investigated by studying the difference of $`r_{\mathrm{fo}}`$ between different observations. In particular, a similar analysis as above, applied to the BeppoSAX observation of Mrk 421 in 2000 May, shows a larger value of $`r_{\mathrm{fo}}`$=1.7, whereas the value is somewhat similar to the ASCA observations for the BeppoSAX observations in 2000 April. This is shown in Figure 10 and 11. An examination of the light curve during 2000 May reveals that flares with shorter time scales are dominant in the light curve. These differences in the characteristic time scales can be seen clearly in the structure functions calculated for each observation, shown in Figure 12. The break in the structure function is at $``$10 ks for the 2000 May data, whereas it is at $``$30–40 ks for the other two. While this might be only a coincidence, the 3 observations show a trend that $`r_{\mathrm{fo}}`$ is larger when $`t_{\mathrm{chr}}`$ is shorter: this is actually the trend expected from our model – when $`\sigma _\mathrm{\Gamma }^{}`$ is larger, $`r_{\mathrm{fo}}`$ is larger and $`t_{\mathrm{chr}}`$ is shorter.

### 4.4 The Efficiency

One concern is that the dynamical efficiency $`\eta `$ of the collisions is rather small. For instance, for the parameters derived above for Mrk 421, $`\eta `$ given by Equation 3 is $`10^6`$ for the collisions that occur at the shortest $`D`$. For the more distant collisions, this becomes even smaller.

It has been claimed that the jet cannot radiate all of its energy in the sub-parsec region considered here, since a substantial power must be transported to the kpc-scale radio lobes. With this, the efficiency must be low. However, this seem too low – as this means that only 10<sup>-6</sup> of the total energy is radiated at the base of the jet, which indicates that the jet is 10<sup>6</sup> times energetic than is observed via the blazar phenomenon. On the other hand, since there is no measure of the energetics of the entire jet, such low radiative efficiency near the core cannot be simply ruled out. Nevertheless, as it does seem rather low, we discuss how this could be increased.

The very low efficiency mainly results from the fact that we intend to reproduce the offset component, which reduces the $`\sigma _\mathrm{\Gamma }^{}`$, and thus reduces the velocity difference of the two colliding shells. Thus, the easiest way to increase the efficiency would be to have larger $`\mathrm{\Gamma }_2\mathrm{\Gamma }_1`$, but this would make the flare time scales shorter than the flare cycles, which conflicts with the observed flares that occur repeatedly (daily-flares).

Part of the low efficiency also is due to the fact that we are assuming a Gaussian distribution for the initial BLF, which emphasizes the effect. The origin of the modulation of the Lorentz factor is not known; it can be due to any physical condition of the central engine that is not stable, such as instabilities in the innermost parts of an accretion disk or magnetic eruptions in the corona (e.g., Sikora & Madejski 2000). Thus one way to increase the efficiency would be to assume a broader distribution, such as a flat distribution within a certain range (e.g., Spada et al. 2001; Kobayashi, Piran, & Sari 1997), which would increase the number of efficient collisions. According to our simulations, this was shown to increase the efficiency by roughly one order of magnitude. We can consider even more extreme distributions where the values of the BLFs are concentrated at the low and high end of the distribution. At most, this could improve the efficiency by another order of magnitude, but not any more, since $`\mathrm{\Gamma }_2\mathrm{\Gamma }_1`$ is always small. If this is still too inefficient, this would indicate that there is a problem assuming random distribution of BLFs for the shells in the internal shock scenario. In this case, we would probably need to consider another possible origin of the offset component. For instance, if there is emission from significantly larger distances, such as by reconfinement shocks or by shocks due to collisions with inhomogeneities that have much larger radial extensions, this may produce a low-amplitude relatively steady component. For this case, the efficiency can be much larger. On the other hand, we emphasize here that the internal shock scenario alone reproduces successfully many of the observed features such as the daily-flares, the existence of the characteristic time scale $`t_{\mathrm{chr}}`$ and the non-existence of variability power below $`t_{\mathrm{chr}}`$.

## 5 THE EXTERNAL SHOCK SCENARIO

In the previous sections we have discussed how the internal shock model successfully reproduces the observed variability of the both the fluxes and spectra of blazars. While the this model works rather nicely naturally explaining the observed properties, we briefly consider the alternative model, the external shock scenario.

In the internal shock scenario, a shock is generated when a faster shell catches up to a slower shell. This is where electrons are accelerated to relativistic energies and radiate through synchrotron and inverse Compton emission. The alternative scenario is that these shells do not collide with each other but run into either an external material or external field, where the shock emerges (e.g. Dermer et al. 1999). This is somewhat similar to the acceleration mechanism which is usually considered in supernova remnants or the afterglows of gamma-ray bursts. In these two examples, the external material is provided by the interstellar medium, while in the case of blazars, the promising candidates for the external medium are broad line clouds and/or intercloud material. We will discuss below how this may work in blazars; below, we will call such external material “clumps.”

Here, we assume a shell ejected with a Lorentz factor $`\mathrm{\Gamma }`$, which runs into a broad line cloud “clumps” (or any other external material) at a distance $`D`$. In similarity to the case of the internal shock scenario, given that the radiation cooling time is much faster than the observed variability time scales, the observed variability time scales are most likely determined by the dynamical time scale. There are two time scales which may affect the dynamical time scale: the time for the shock to cross the region, and the angular spreading time. Here, in contrast to the internal shock scenario, the shock can be considered as relativistic.

For this case, the time for the shock to cross the shell is given by

$$t_{\mathrm{crs}}=\frac{l}{2c\mathrm{\Gamma }^2},$$
(17)
where $`l`$ is the shell thickness. The angular spreading time is similar to the case of the internal shock model, and thus

$$t_{\mathrm{ang}}=\frac{R}{c\mathrm{\Gamma }}=\frac{D}{c\mathrm{\Gamma }^2},$$
(18)
where $`R`$ is the radius of the shell at distance $`D`$.

If the shock crossing time $`t_{\mathrm{crs}}`$ were to determine the observed time scale of the flares, which is of the order of a day, Equation 17 would suggest that the thickness of the shell must be as large as $`l`$$`3\times 10^{17}`$ cm. This means that the central engine must be continuously ejecting material for 100 days to produce a single shell – which is in severe conflict with the observed rate of flares which occur daily. Thus, in contrast to the case of the internal shocks, the angular spreading time $`t_{\mathrm{ang}}`$ should be the dominating time scale in the case of external shocks. For this case, to obtain flares with time scales of $``$ 1 day, the location of the shock is calculated to be at $`D`$$`3\times 10^{17}`$$`\mathrm{\Gamma }_{10}^2`$ cm, where $`\mathrm{\Gamma }_{10}`$=$`\mathrm{\Gamma }/10`$.

Because the emission is beamed, the observed time scale of $`1`$ day reflects a $`10`$ day emission in the jet co-moving frame. The daily occurrence of flares indicates that there must be at least 10 shells radiating at the same time, which requires at least 10 separate “clumps.” The external fields must also be restricted in a rather limited range, as to keep the flare time scales similar. Figure 13 shows a simulated light curve for the case of $`\mathrm{\Gamma }`$=15, and $`D`$=$`5.4\times 10^{17}`$ cm (calculated for $`t_{\mathrm{shot}}=t_{\mathrm{ang}}=80`$ ks). Here, the average flare cycle is set to 5 ks, showing that it is possible for the offset component to be produced only if the flare cycle is high (i.e. substantial number of “clumps” of the external material available for collisions with the jet at the same time).

We also remark that for the case of the external shock model, the collisions are likely to occur farther away from the central source than for the internal shock model. A concern arises to whether the broad line clouds are still present at such distances. From observations of the time delay between the continuum and line variations in AGNs, the broad line regions are suggested to extend to distances around 10<sup>16-18</sup> cm, although there are no measurements for blazars. If this is of the same order for blazars, the distance calculated above (for $`\mathrm{\Gamma }`$=15, $`t_{\mathrm{shot}}`$=1 day) is actually at the right location, but may be close to the limit. Summarizing, we do not claim that the discussion above rules out the external shock scenario, but we remark that more fine-tuning of the parameters is necessary as compared with the internal shock model.

## 6 SUMMARY

The high-quality light curves and spectra obtained via the recent long look observations of TeV blazars provided the first opportunity to use the variability as a new tool to study the structure of jets in blazars in more detail than was previously possible. Our approach is unique as it investigates not only a single flare, but considers a series of flares resulting from well-sampled, long-duration time series. We summarized the observed variability properties and suggested a physical interpretation to explain these features based on the framework of the internal shock scenario.

We presented a simplified model applied specifically for TeV blazars, and investigated through simulations how each of the parameters would affect the observed light curve or spectrum. In particular, we showed that the internal shock scenario naturally accounts for the observed variability properties, but it requires a condition that the fluctuation of the initial bulk Lorentz factors are small, $`\sigma _\mathrm{\Gamma }^{}0.01`$. Remarkably, this explains both features observed in the flux but also in the spectral variability – the repeatedly occurring flares, the offset component, the structure function, and the shift of the synchrotron peak frequency. We also showed that several useful observational quantities can be used to probe the physical parameters of the inner jet: the characteristic time scale, the flare-to-offset ratio, and the slope of the structure function. We applied this model to the ASCA X–ray light curves of the TeV blazar Mrk 421, which allows us to determine several physical parameters of the jet such as the frequency of ejection of shells (on average one shell every 10 minutes), the average thickness of the shell ($`10^{13}`$ cm), and the location of their collisions (typically $`10^{17}`$ cm away from the central source). We also briefly commented on the external shock scenario, and claimed that this scenario is viable, but requires rather detailed fine-tuning to the parameters.

We thank Marek Sikora and Fumio Takahara for valuable discussions on this manuscript. We would also like to thank the anonymous referee for constructive comments to improve this paper. Support for this work was provided by the Fellowship of Japan Society for Promotion of Science for Young Scientists, and by NASA via Chandra grant no. GO0-1038A from SAO to Stanford University.