format: text

text:
# Melting and structure of the vortex solid in strongly anisotropic layered superconductors with random columnar pins

## I Introduction

The mixed phase of type-II superconductors with random pinning is generally regarded to be an archetypal test system for the study of the effects of quenched disorder on the structure and melting of solids. In this phase, magnetic flux penetrates the sample as quantized vortex lines which form a special physical system known as “vortex matter”. The fascinating equilibrium and dynamical properties of vortex matter in the mixed phase of high-temperature superconductors (HTSCs) have prompted considerable experimental and theoretical attention for more than a decade (see Ref. blat,  for an early review). The mixed phase of HTSCs is strongly affected by the presence of pinning centers, either intrinsic or artificially generated. Understanding the effects of pinning in these systems is basic for practical applications because pinning strongly influences the value of the critical current dg .

The vortex lines in a pure type-II superconductor form a triangular lattice (the Abrikosov lattice) at low temperatures. Because of enhanced thermal fluctuations in highly anisotropic, layered HTSCs, the Abrikosov lattice in very pure samples, in a magnetic field perpendicular to the layers (the field will be assumed to be in this direction throughout our discussion), undergoes a first-order melting transitionblat  into a resistive vortex liquid (VL) as the temperature $`T`$ increases. Random pinning destroys the long-range translational order of the Abrikosov lattice and leads to the occurrence of a variety of glassy phases at low $`T`$. It is now generally accepted that in systems with random point pinning, a topologically ordered low-temperature phase with quasi-long-range translational order (denoted as the “Bragg glass” (BrG) phase) occurs at low fields if the pinning disorder is sufficiently weak. This has been established theoretically nat ; gia ; ns  as well as experimentally (see e.g. Ref. kle, ). The possibility of an amorphous vortex glass (VG) phase, with nonlinear voltage-current characteristics and vanishing resistance in the zero-current limit, in systems with strong pinning (or at high magnetic fields where the effects of pinning disorder are enhanced) was suggested by Fisherfis  (see also Ref. ffh, ). However, in spite of extensive investigations, the existence of a true VG phase (i.e., an amorphous glassy phase thermodynamically distinct from the high-temperature VL) in systems with uncorrelated point pinning remains very controversial: different calculations bok ; ves  lead to different conclusions and the experimental situationpet ; str  is similarly contradictory. A variety of “glassy” behavior has been reported in different experiments (see Ref. ban1,  and references therein) near the first-order melting transition of the BrG phase of both conventional superconductors and HTSCs with point pinning. It has been suggestedban1 ; men  that these observations may be understood if it is assumed that the melting of the BrG phase occurs in two steps: the BrG first transforms into a “multidomain” glassy phase that melts into the usual VL at a slightly higher $`T`$.

Columnar pinning defects can be produced for example as damage tracks arising from heavy-ion bombardment. The technological importance of these extended defects oriented parallel to the direction of the external magnetic field has been long recognizedciv ; bud : they are highly effective in increasing the critical current by localizing vortex lines along their length. Heavy-ion irradiation produces a random array of parallel columnar defects each of which can trap one or more vortices at low temperatures. The effects of such random arrays of extended defects, oriented perpendicular to the superconducting layers, on the properties of the mixed phase of HTSCs have been extensively studied experimentallykha ; ban ; meng . The same question has also been examined theoreticallynel ; rad ; lar ; vin  and numericallywen ; sen ; tya ; non . When the columnar pinning is strong, and their concentration exceeds that of the vortex lines, a so-called “strong” Bose Glass (BoG) phase with nearly all the vortex lines localized at the pinning centers occursnel ; wen  at low $`T`$. This phase is strongly disordered with very short-range translational and bond-orientational correlations. The behavior in this regime is fairly well-understood in terms of a mappingnel  of the thermodynamics of the system to the quantum mechanical properties of a two-dimensional system of interacting bosons in a random potential.

Much less is known about the behavior in the dilute pin limit where the concentration of pins is much smaller than that of vortex lines. In such systems, one expectsrad  a “weak” BoG phase at low temperatures, with a small fraction of vortex lines localized strongly at the pinning centers and the remaining ones localized relatively weakly in the interstitial region between pinning centers. As the temperature is increased, this phase shouldrad ; vin  melt into an interstitial liquid (IL). In the IL phase, some of the vortices remain trapped at the pinning centers, while the other, interstitial ones, form a liquid. The pinned vortices are expectedrad ; vin  to delocalize, thereby forming the usual VL, at a crossover occuring at a higher temperature. The melting transition of the “weak” BoG phase into the IL is predictedrad ; lar  to be first-order for small pin concentrations, whereas the “strong” BoG to VL transition occurring for large pin concentrations is knownnel ; lar  to be a continuous one.

A first-order melting transition of the “weak” BoG phase has been observedmeng  in recent experiments on samples of $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_{8+\mathrm{x}}`$ (BSCCO) with a small concentration of columnar pins. The BoG phase in these systems is foundban ; meng  to have a polycrystalline structure with ordered vortex crystallites of different orientations embedded in the interstitial region between vortices pinned at the columnar defects. If the pin concentration is sufficiently small, the melting of this BoG phase upon increasing the temperature occurs into an IL phase with a fraction of vortex lines remaining pinned at the defects. These pinned vortices delocalize at a higher temperatureban . The temperature at which the first-order BoG to IL transition occurs iskha ; ban  very close to the melting temperature of the same system without the columnar pins, if the pin concentration is small. As the pin concentration is increased, the melting temperature increases and the transition eventually becomes continuous. The difference between the melting temperature and the temperature at which the pinned vortices in the IL phase delocalizes decreases and eventually goes to zero as the pin concentration increases. Another interesting feature, found experimentally for both pointsoi  and columnarban  pinning, is that the melting of the low $`T`$ “solid” phase is “inhomogeneous” in that it occurs locally over a range of temperatures: the local transition temperature, which can be measured from the local magnetization, is different in different regions of the sample. This inhomogeneity of the local melting temperature is believedsoi  to be closely related to the local arrangement of the specific pinning centers (the so-called “pinning landscape”) in each sample studied.

The polycrystalline nature of the BoG phase for small pin concentrations has also been observed in simulationssen . Recent numerical worknon  indicates that a BrG phase, with topological order, may also exist in such systems provided that the pin concentration is sufficiently small. In these simulations, the BrG phase is found to melt into the IL phase via a first-order transition as the temperature is increased. As the pin concentration is increased beyond a critical value, the BrG phase disappearsnon  and a low-temperature BoG phase with a continuous BoG to IL transition upon increasing the temperature is foundtya ; non .

In this paper, we approach these problems through a different numerical method. We consider a layered, strongly anisotropic superconductor (such as BSCCO) with a dilute random array of columnar defects oriented perpendicular to the layers, and with a magnetic field applied parallel to the columnar pins. We describe the equilibrium properties of this system in terms of vortex density variables, using a free energy functional of the Ramakrishnan-Yussouffry  (RY) form. We numericall minimize a spatially discretized version of this free energy functional and obtain the vortex density configuration at each local minimum. Analysis of these density configurations, both in terms of a variety of correlation functions and by direct visualization of the arrangement of the vortices (vortex positions are identified with those of local peaks of the vortex density), allows us to identify the nature of the phases corresponding to different local minima of the free energy. Comparison of the values of the free energy at the minima as the temperature varies yields the transition temperatures. A similar procedure has been successfully employed for the case of a regular arraydv1 ; dv2  of columnar pinning centers.

We find in our study that that a Bragg glass phase exists at low temperatures for samples with a small concentration of columnar pins. We will present evidence for this phase from our analysis of the vortex density and associated correlation functions. Upon warming, this phase melts into what we show to be an insterstitial liquid (IL) phase, but the melting is showndv3  to occur in two steps: the BrG and IL phases are separated, over a narrow temperature range, by an intermediate phase which exhibits a multi-domain polycrystalline structure. We show that this intermediate phase should be identified as a Bose glass (BoG). The temperature of the upper (BoG to IL) transition is approximately independent of the columnar pin concentration $`c`$, at the low $`c`$ values studied. Both the BrG to BoG and the BoG to IL transitions are found to be first-order. We also find that a local melting temperature can be suitably defined using a criterion based on the degree of localization of the vortices, and that its behavior (in particular, its spatial variation) is quite consistent with what is seen in experimentsban ; soi . The value of the local melting temperature is strongly correlated with the presence of topological defects (dislocations) in the vortex solid which, in turn, is correlated with the local arrangement of the pinning centers. We also show that the transition to the IL phase corresponds to a percolation of regions containing liquid-like (delocalized) vortices across the sample.

After this Introduction, we discuss our definitions, model and numerical methods in the next Section. Then in Section III we present our results. We discuss first the free energy minima and their study through the use of tools such as correlation functions, local peak-density plots, and Voronoi construction for the lattice formed by local density peaks. These tools allow the identification of the phases at each free energy minimum, as we shall show. Then we derive the phase diagram and show that indeed the BrG and IL phases are separated by a thin sliver of BoG. Finally, our results for the nature of the local melting are displayed and discussed. A brief conclusions Section recapitulates the main points of the paper and discusses them in the context of existing results.

## II Methods

The general procedures that we use are quite similar to those employed in previous work on a regular array of columnar pinsdv1 ; dv2 . We will therefore give here only a brief summary, emphasizing the details that are different in the random case considered here.

The system we study is a layered superconductor in the extreme anisotropic limit, that is, with vanishing Josephson coupling between layers, which are then coupled via the electromagnetic interaction only. This limit is appropriate for the Bi- and Tl- based HTSC compounds in a large region of the magnetic field ($`H`$)–$`T`$ plane. In this work, we will use material parameter values appropriate to BSCCO. With these assumptions, one can write the energy of a system of “pancake” vortices residing in the layers as a sum of anisotropic two-body interactions. For straight columnar pins normal to the layers, the pinning potential is the same in every layer. It is then possible to write the free energy as a functional of the time averaged areal vortex density $`\rho (𝐫)`$:

$$F[\rho ]F_0=F_{RY}[\rho ]+F_p[\rho ].$$
(1)
where $`F_0`$ is the free energy corresponding to a uniform vortex liquid of density $`\rho _0=B/\mathrm{\Phi }_0`$ ($`B`$ is the magnetic induction and $`\mathrm{\Phi }_0`$ the superconducting flux quantum). The first term in the right-hand side of Eq. (1) is the free energy functional in the absence of pinning. As explained above we use for this free energy the RYry  functional, which is knowndv1 ; dv2 ; seng ; men1  to give a quantitatively accurate description of the vortex-lattice melting transition in the absence of pinning. It is of the form:

$`\beta F_{RY}[\rho ]`$ $`=`$ $`{\displaystyle d^2r\left[\rho (𝐫)\{\mathrm{ln}(\rho (𝐫))\mathrm{ln}(\rho _0)\}\delta \rho (𝐫)\right]}`$ (2)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^2rd^2r^{}\stackrel{~}{C}(|𝐫𝐫^{}|)\delta \rho (𝐫)\delta \rho (𝐫^{})}.`$
Here $`\beta `$ is the inverse temperature, $`\delta \rho (𝐫)\rho (𝐫)\rho _0`$, and $`\stackrel{~}{C}(r)`$ is the usual direct pair correlation functionhan  which may be written as a sum over layers: $`\stackrel{~}{C}(r)_nC(n,r)`$, with $`C(n,r)`$ (where $`n`$ is the layer separation and $`r`$ the in-layer distance) being the corresponding direct pair correlation function of a layered liquid of pancake vortices. The direct correlation function, which is needed as input in our free energy, can be accurately calculated in a number of ways. We will use here the results of the hypernetted chain calculation of Ref. men1, . In general, the results in the limit considered depend on the values of the in-plane London penetration length $`\lambda (T)`$, the interlayer spacing $`d`$ and a dimensionless coupling parameter $`\mathrm{\Gamma }`$ given by:

$$\mathrm{\Gamma }=\beta d\mathrm{\Phi }_0^2/8\pi ^2\lambda ^2(T).$$
(3)
For BSCCO we will take $`d=15\AA `$, and we will assume a standard two-fluid temperature dependence for $`\lambda (T)`$ with $`\lambda (0)=1500\AA `$ and $`T_c=85`$K (at zero field).

The second term in the right-hand side of Eq. (1) represents the pinning and is of the form:

$$F_p[\rho ]=d^2rV_p(𝐫)\delta \rho (𝐫).$$
(4)
where $`V_p`$ is the pinning potential which can be written as

$$V_p(𝐫)=\underset{j}{}V_0(|𝐫𝐑_j|),$$
(5)
with the sum extending over the planar positions of the random pinning centers. We take the potential $`V_0`$ corresponding to a single pinning center to be of the usual truncated parabolic form:daf

$$\beta V_0(r)=\alpha \mathrm{\Gamma }[1(r/r_0)^2]\mathrm{\Theta }(r_0r)$$
(6)
where $`r_0`$ is the range. The basic length in the problem, which we will use as our unit of length unless otherwise indicated, is $`a_0`$, defined by the relation $`\pi a_0^2\rho _0=1`$. We choose $`r_0`$ to be $`r_0=0.1a_0`$ and take the dimensionless constant $`\alpha =0.05`$, for which value, as previously shown,dv2  each pinning center traps approximately one vortex in BSCCO, in the temperature range of interest. This range is determined by the following considerations. We will keep the field fixed at $`B=0.2`$T, and vary the temperature. The melting temperature of the unpinned lattice is then $`T_m^018.4`$Kdv2 , and we, therefore, consider the neighborhood of this temperature. Then one has, to a very good approximation $`\mathrm{\Gamma }=2650/T`$ (with the temperature in Kelvin), while $`\lambda \lambda (0)`$.

To carry out our numerical calculations, we discretize the density variables on a computational triangular grid of lattice spacing $`h`$, containing $`N^2`$ sites. We take the spacing $`h`$ to be $`h=a/16`$ in our calculations, where $`a=1.998a_0`$ is the equilibrium lattice constant of the vortex lattice in the temperature range considered at the indicated field. This value, as pointed out in Ref. dv2, , is slightly higher than that of a triangular vortex lattice of density $`\rho _0`$, which is $`(2\pi /\sqrt{3})^{1/2}a_0`$. We define at each site $`j`$ on this computational lattice a variable $`\rho _j`$, with $`\rho _j\rho (𝐫_j)v`$, where $`v`$ is the area of each computational unit cell. Results reported here are for $`N=1024`$, which, for the chosen value of $`h`$, corresponds to including $`N_v=`$ 4096 vortices in the calculation. Preliminary and confirmatory results at $`N=512`$ (1024 vortices) were also obtained. A number of pinning sites ($`N_p`$) are put, at random, on some of the computational lattice sites. The results presented here, which correspond to the dilute limit, will be primarily for a pin concentration $`c`$ of $`c=1/64`$, that is, a number of pins $`N_p=64`$ at $`N=1024`$ (or $`N_p=16`$ at $`N=512`$). Some preliminary results for $`c=1/32`$ ($`N_p=128`$ at $`N=1024`$, or $`N_p=32`$ at $`N=512`$) will also be discussed briefly. With random pins, the results depend on the particular random pin configuration and averaging over different such configurations is required. The dependence is however not strong: as we shall see, averaging over five to ten configurations is enough to make statistical errors sufficiently small for our purposes.

To perform our studies, we numerically minimize the discretized free energy with respect to the $`N^2`$ discretized density variables $`\{\rho _i\}`$. To do so, the interaction term in the right side of Eq. (2) must be repeatedly evaluated, and since this term is of a convolution form, this is most efficiently done in momentum space, through the use of efficient Fast Fourier Transform (FFT) routines. This avoids having to evaluate this term as a double sum which would be computationally much more cumbersome than performing the direct and inverse FFT’s. In performing the minimization, one must keep in mind that the variables $`\{\rho _i\}`$ must be nonnegative. This precludes the use of many efficient minimization algorithms. We use a procedurenum  that ensures that this constraint is satisfied. Numerical minimization is performed starting with an appropriate initial condition for the density variables. As explained in Sec. III below, different local minima may be found, at the same field and temperature, depending on the initial conditions. These minima correspond to different phases. The minimization procedure yields not only the value of the free energy at each minimum but also the detailed vortex density configuration at the minimum found, i.e. the values of the set $`\{\rho _i\}`$ at the minimum. It is then straightforward to analyze the actual density configuration in several ways. It is also possible to evaluate any desired density correlations. The nature of the phase corresponding to each local minimum of the free energy can be inferred from such analysis. The (mean-field) phase diagram is then obtained from a comparison of the free energies of the different minima as the temperature $`T`$ is varied. The results of carrying out this program are discussed below.

## III Results

### III.1 Free energy minima

As explained above, minima of the free energy are found by starting the minimization process with appropriate initial conditions. We use three different kinds of initial conditions. The first kind is a uniform density ($`\rho _i=\rho _0v`$ for all $`i`$), corresponding to a completely disordered liquid state. This is typically used to obtain, as we shall see, liquid-like states at temperatures near the equilibrium melting temperature of the pure vortex lattice, which isdv1 ; dv2  $`T_m^018.4`$K for the value of $`B`$ considered here. The second kind corresponds to a crystalline initial state, with values of $`\{\rho _i\}`$ as obtained in Ref. dv2,  for the pure vortex lattice. In our system with a spatially varying random pinning potential, the pinning energy of such a crystalline density configuration depends on the choice of the computational lattice sites at which the periodic local peaks of the density are located. We take, among all possible choices equivalent by symmetry operations in the computational lattice, the one for which the pinning energy is lowest for the specific pinning configuration under consideration. These initial conditions are used to obtain ordered states at lower temperatures. Finally, states originally obtained by either of these two procedures can be slowly warmed up, or cooled down: in this case the initial conditions are simply the final state obtained at the previous temperature. As one warms up or cools down a state, its nature in general changes: liquid configurations may become unstable upon cooling and ordered configurations upon warming.

Whichever initial conditions one uses, a local minimum of the free energy is eventually found. At a given temperature, in general several local minima with different density configurations, characterized by the values of the set of $`\{\rho _i\}`$ variables, are obtained. The minimum with the lowest free energy represents, at a given $`T`$, the true equilibrium state. It is obviously necessary to develop systematic procedures to identify the structure of each of these diverse minima. Since we have access to the full set of density variables at the computational lattice sites, we have at our disposal a variety of methods to achieve that goal. The first is to calculate the density correlations, e.g. the structure factor $`S(𝐤)`$:

$$S(𝐤)=|\rho (𝐤)|^2/N_v$$
(7)
where $`\rho (𝐤)`$ is the discrete Fourier transform of the set $`\{\rho _i\}`$. Equivalently, one can consider the Fourier transform of $`S(𝐤)`$, which is the two-point spatial correlation function $`g(𝐫)`$ of the time-averaged local density. We will present here results for $`S(𝐤)`$, considered as a function of the vector $`𝐤`$, and for the angularly averaged spatial density correlation function, $`g(r)`$.

It is also very useful to consider the spatial structure formed by the vortices at low temperatures. For this purpose we have extracted from $`\{\rho _i\}`$ the local peak densities. We say that the density locally peaks at site $`i`$ if the value of $`\rho _j`$ at $`j=i`$ is higher than that at any other computational lattice site within a distance $`a/2`$ from $`i`$, where, as defined above, $`a`$ is the equilibrium spacing of the unpinned vortex lattice. As expected, we find that at low-temperature solid-like minima where the vortices are strongly localized, the number of local density peaks matches the number of vortices $`N_v`$. The positions of these peaks determine what we will call the “vortex lattice”.

Much useful information about the spatial structure of this vortex lattice may be obtained by plotting directly the local density values at the vortex lattice points. An excellent, and complementary, alternative to help elucidate the degree of vortex lattice order is to carry out a Voronoi construction: we recall that in the Voronoi construction one determines cells around each lattice point by a Wigner-Seitz procedure. In a perfect crystal, all the resulting Wigner-Seitz cells are identical, while in the general case the cells can have different sizes and shapes. The number of sides of the Wigner-Seitz cell surrounding a lattice point is identified with the number of “nearest neighbors” of that lattice point, and the difference between this number and the average (six in our case) marks the position of defects (disclinations). We will make use of such plots below.

Since we will be concerned not only about translational, but also with orientational order, we will also examine the bond-orientational correlation function $`g_6(r)`$ in the vortex lattice defined as:

$$g_6(r)=\psi (𝐫)\psi (0)$$
(8a)
where the angular brackets denote overall average over the vortex lattice and the field $`\psi (𝐫)`$ is given by:
$$\psi (𝐫)=\frac{1}{n_n}\underset{j=1}{\overset{n_n}{}}\mathrm{exp}[6i\theta _j(𝐫)]$$
(8b)
with $`\theta _j(𝐫)`$ being the angle that the bond connecting a vortex at $`𝐫`$ to its $`j`$-th neighbor makes with a fixed axis, and $`n_n`$ is the number of neighbors of the vortex at $`𝐫`$.

It is also possible, and as we shall see, useful, to define what we will call a “translational correlation function” in the vortex lattice in a way that is quite analogous to the definition of $`g_6(r)`$. This function, which we denote as $`g_G(r)`$, is defined by an equation identical to the right-hand side of Eq. (8a), but with the field $`\psi `$ being defined as

$$\psi _G(𝐫)=\mathrm{exp}(i𝐆𝐫)$$
(9)
where $`𝐆`$ is a reciprocal lattice vector of the triangular vortex lattice in the absence of pinning. We will consider here only the case where $`𝐆`$ is a shortest nonzero reciprocal lattice vector and average over the results obtained for the six equivalent $`𝐆`$’s.

Our identification of the different phases is based on analysis of all this information. We show some of the results in the next few Figures. These all correspond to samples with 64 pins and 4096 vortices. First, in Fig. 1, we show the structure factor, as defined in Eq. (7). The three sets of results shown there are for local free energy minima at the same temperature, $`T=18.4`$K, obtained with different initial conditions of the three kinds described above, and the same pin configuration. Other pin configurations yield very similar results (a different example is shown in Fig. 1c of Ref. dv3, ). The green symbols correspond to the free energy minimum obtained by starting with uniform initial conditions. Clearly, the structure factor is completely featureless and liquid-like in this case – the absolute value of $`S`$ never exceeds five. The red circles are for a local minimum obtained with initial conditions corresponding to the best crystalline state, as explained above, for this pin configuration. The structure factor now exhibits typical ordered behavior, highlighted by six sharp Bragg-like peaks which are emphasized by the added vertical lines in the Figure. Finally, the blue triangles correspond to a minimum obtained by first “quenching” with uniform initial conditions to a relatively low temperature (16.8K) where the liquid-like state is found not to be stable, and then slowly warming, at 0.2K temperature intervals, back up to 18.4K. As one can see, the $`S(𝐤)`$ for this minimum has an intermediate structure, with several relatively well-defined peaks, more than six in number, whose heights are considerably lower than those of the peaks found for the ordered minimum.

In Fig. 2, results for the angularly averaged density correlation function $`g(r)`$ are displayed. These results are for minima obtained at the same temperature, and using the same initial condition procedures, as those for the data shown in Fig. 1. However, the results shown here are averages over five different pin configurations. The nearly flat (green) curve in the top panel corresponds to minima obtained from quenches with uniform initial conditions, the (red) curve with the well-defined peaks shown in the bottom panel corresponds to minima obtained with crystalline initial conditions, and the (blue) curve with the intermediate peaks shown in the top panel is for minima obtained by quenching with uniform initial conditions to a low temperature and subsequent slow warming (very similar results can alternatively be obtained by slowly cooling a high-temperature liquid-like minimum to a temperature where it is unstable). One can see then that the results for $`g(r)`$ are fully consistent with those found in Fig. 1: the minimum obtained from uniform initial conditions is fully disordered, while that obtained from crystalline initial conditions exhibits a large degree of crystalline order. For the the third kind of minima, we find some degree of intermediate range order.

We next examine the structure of the minima in terms of the “vortex lattice” found by the procedure explained above. Sample results are shown in Fig. 3, all for the same pin configuration, at the temperature $`T=18.2`$K. Results for three minima, obtained from the same procedures and initial conditions as the three minima in the previous two figures, are shown. The positions of the local density peaks and the values of the density at these peaks are displayed through a symbol and color coding scheme described in the caption of the Figure. Blank regions denote areas where no local peak was found using the algorithm described above: such areas are mostly found only in the most disordered case (uniform initial conditions). We see that the results are quite consistent with those obtained from the previous correlation functions: quenching with the appropriate crystalline initial conditions leads to a fairly well-ordered (but not perfectly ordered) lattice, while from uniform initial conditions, one obtains a rather disordered, liquid-like structure where the density is nearly uniform and close to $`\rho _0`$, except near the pinning centers at each of which one vortex is trapped. An intermediate result, on the whole more solid than liquid-like, but definitely disordered, is obtained through the slow warming or cooling scheme.

To better gauge the degree of disorder present in each case, we construct, in Fig. 4, the corresponding Voronoi plots. To display the behavior in a more obvious way, the temperature in these plots is lower by $`0.2`$K than that of the plots in Fig. 3. As mentioned above, sites with number of nearest neighbors ($`n_n`$) different from six represent the locations of disclinations. In the solid-like minima (top two panels of Fig. 4), nearly all the sites with $`n_n6`$ have $`n_n=5`$ or $`n_n=7`$ (the disclinations have unit “charge”). Also, six- and seven-coordinates sites, indicated respectively by (red) triangles and (green) solid circles in the plots, always occur in nearest-neighbor pairs. Such pairs correspond to dislocations. In the first panel, which displays the results for the most ordered state obtained from crystalline initial configurations, we see that the number of dislocations, while nonzero, is quite small and that they form fairly isolated small clusters, each of which has zero net Burgers vector. These defect clusters occur near the pinning centers and they represent the local disruption of crystalline order due to the pinning. In the second panel, which shows the results for the minimum with intermediate order, the dislocations are organized to form well-defined grain boundaries that separate crystallites with different orientations. As a result of this polycrystalline structure, both translational and bond-orientational correlations are short-range. Performing the same construction for other pin configurations (results for a different pin configuration are shown in Fig.1b of Ref. dv3, ), one finds that the crystallite arrangement depends on the pin configuration, as one would expect. The grain boundaries tend to lie away from pinning centers, which makes sense physically. Finally, in the last panel, which shows the results for the minimum obtained from uniform initial conditions, the high degree of disorder is evident. The Voronoi construction shown in this panel is not very meaningful because our method of obtaining the vortex positions is not reliable for liquid-like minima. The identification of the position of a vortex with that of a local peak of the density field is justified only when the peak is sharp. This is the case in the solid-like minima, but not so in the liquid-like minima with low peak densities. We have shown here the results for the liquid-like minimum only for the purpose of illustrating the difference between the structure of this minimum and those of the more ordered ones. For this minimum, we find a very large number of defects, and it is difficult to determine whether disclinations of opposite “charge” always pair up to form dislocations. Inspection of the defect distribution suggests the presence of free disclinations, but not conclusively. The total number of local peaks of the density field in liquid-like minima turns out to be substantially smaller than the expected number of vortices. This problem is not present at the solid-like minima where nearly all the vortices are strongly localized.

Next, we show in Fig. 5 examples of the bond-orientational and “translational” correlation functions, $`g_6(r)`$ and $`g_G(r)`$ as defined in Eqs. (8) and (9) respectively. All results shown in this figure are averages over five pin configurations at temperature $`T=17.6`$K. At this temperature the liquid is unstable, and in any case all liquid correlation functions are featureless, so that case is not shown. The (purple) triangles correspond to ordered minima obtained from quenching with crystalline initial conditions to $`T=16.8`$K and subsequent slow warming, and represent the bond-orientational correlation, $`g_6(r)`$. This function appears to saturate at a fairly large value as $`r`$ increases, indicating the presence of long-range bond-orientational order. The (red) open circles display $`g_6(r)`$ for the polycrystalline state with intermediate order, obtained through the procedures described above. This function decays exponentially with distance, as shown by the exponential fit (dashed line), over a distance of few tens in units of $`a_0`$. Finally, the (black) open squares show the results for the “translational” correlation function $`g_G(r)`$ for the ordered minima. The dotted line is an exponential fit to this translational correlation function for $`r/a_0<10`$, showing that the decay at longer distances is clearly slower than exponential (possibly a power-law).

¿From the overall examination of data such as those shown in these figures, we can reach the following conclusions: quenching with uniform initial conditions to temperatures above or slightly below $`T_m^0`$ leads to very disordered states, with very short range correlations and no structure, except for vortices at and near the pinning centers. These minima become unstable upon slow cooling: one then obtains minima of the third kind, discussed below. They clearly must be identified with the liquid phase, specifically the interstitial liquid (IL) phase discussed in the Introduction. The second kind of minima are obtained by quenching with crystalline initial conditions to any temperature below, or slightly above, $`T_m^0`$. These minima can then be cooled, without losing their character, but they melt into the IL upon sufficient warming. The states thus obtained are obviously nearly crystalline, with the reservation that the order is not truly long range, but has a slow decay. This is most evident in the results for $`g_G(r)`$ seen in Fig. 5. The bond-orientational order is nearly perfect. Defects are limited to isolated clusters. This state must therefore be identified as a Bragg glass (BrG). Finally, we have the states obtained as described above, by slow cooling of the IL state, which can also be obtained by quenching with uniform initial conditions to a temperature below the IL state stability limit. Either procedure yields, for any pin configuration, states that at a given $`T`$ differ only slightly in free energy or density configuration. The states do differ somewhat more, however, for different pin configurations. The Voronoi construction conclusively shows (see Fig. 4) that these states are polycrystalline. The grain boundaries are formed by dislocation chains. For any pin configuration, these grain boundaries survive thermal cycling across $`T_m^0`$, but if sufficiently warmed up, these minima melt, the melting beginning locally at the grain boundaries (see Sec III.2 below). The bond-orientational correlation function for these minima has an intermediate range. This, and also the observation that these states are not unique (as mentioned above, slightly different states are reached depending on the cooling or quenching protocol), clearly indicate a glassy state. We therefore identify it with the Bose glass (BoG) which has been experimentallyban ; meng  and numericallysen  shown to be polycrystalline.

### III.2 Properties of the condensed phases

In this section, we describe in detail how the properties of the low-temperature “solid-like” condensed phases represented by the BrG and BoG types of minima defined in the preceding section vary as the temperature is changed. The high-temperature liquid phase represented by the IL-type minima is not very interesting from the point of view of its temperature dependence: in the temperature range we have considered, this phase exhibits a liquid-like (nearly uniform) density distribution except in the immediate vicinity of the columnar pins, at each of which a vortex is trapped. These trapped vortices would eventually get delocalized at higher temperatures. This happens beyond the upper limit of the temperature range (16.8K $`T`$ 18.8K) we have considered here.

To study how the properties of a local minimum of the free energy change as the temperature is varied, we have “followed” minima as the temperature changes: for example, starting with a minimum obtained by quenching to the lowest temperature of 16.8K, we “follow” that minimum to higher temperatures by increasing the temperature in small steps (usually taken to be 0.2K) and finding a new minimum at the higher temperature by running the minimization routine with the configuration at the minimum obtained at the previous temperature as the initial state. This procedure leads to a new minimum of the same type as long as the minimum remains locally stable – if the temperature is increased to values substantially higher than $`T_m^0`$, where the BrG and BoG minima become unstable, this procedure leads then to the IL minima. We have also carried out “cooling” runs where the temperature was decreased in small steps, starting from a minimum obtained at a relatively high temperature. When the IL minima become unstable, BoG states are obtained. As long as the minimum under consideration did not become unstable in the range of temperatures considered in the runs, we did not find any substantial difference between the results obtained in the “heating” and ”cooling” runs

To characterize the density distribution at a minimum, we looked at “global” quantities such as the structure factor $`S(𝐤)`$, an example of which is shown above in Fig. 1. The $`T`$-dependent information contained in $`S(𝐤)`$ is, however, difficult to display in plots. We therefore introduce a closely related “order parameter” $`m`$ defined as

$$m\sqrt{(S_{max}/N_v)}$$
(10)
where $`S_{max}`$ is the largest value of the structure factor $`S(𝐤)`$ averaged over the six $`𝐤`$-vectors related by lattice symmetry to each $`𝐤`$. By definition, the order parameter $`m`$ is equal to unity in a state with perfect crystalline order (triangular lattice with $`\delta `$-function peaks). Since no long-range crystalline order is expected or found in either one of the BrG and BoG phases, the value of $`m`$ should go to zero at all temperatures in the thermodynamic limit. However, the infinite size limit is reached very slowly in the glassy phases, and for the finite-size systems considered here, this quantity provides a convenient measure of the degree of local order present at the minima.

We show in Fig. 6 the results for the order parameter $`m`$, averaged over five different pin configurations, for samples with 4096 vortices and 64 pins. The results for the BrG, BoG and IL minima are shown by (red) plus signs, (green) crosses and (blue) asterisks, respectively. Only two data points are shown for the IL minimum because it becomes unstable at lower temperatures. For all three kinds of minima, the value of $`m`$ decreases with increasing $`T`$, as expected. The rate of change is largest for the BrG minima, and smallest for the IL minima. However, the differences among the three kinds of minima in the degree of local order remain quite clear at all the temperatures considered. We have also obtained BrG-type minima for several ($`10`$) samples with 1024 vortices and 16 pins. The values of $`m`$ obtained from an average over these smaller samples are found to be only about 6% larger than those obtained for the 4096-vortex samples at the corresponding temperature. This result implies that the translational correlation function at the BrG minima falls off very slowly with distance. This is consistent with our interpretation of these minima as representing the BrG phase.

We also examined the temperature dependence of the detailed density distribution at the minima (obtained from “peak density” plots similar to those shown in Fig. 3) and their defect structure obtained from Voronoi plots such as those in Fig. 4. Examination of the peak-density plots reveals more information about how the density distribution at the minima changes with temperature. In Fig. 7, we have shown peak-density plots at two temperatures (17.4K in the top panel and 17.8K in the bottom panel) for the BoG minimum of the same sample for which a similar plot at 18.2K is shown in the middle panel of Fig. 3. The symbols and color scheme used in these plots are the same as those in Fig. 3. It is clear from these plots that the typical values of the local peak densities decrease (the vortices become less localized) as the temperature increases. The largest values of the local peak density occur near the pinning centers and the smallest values appear near the grain boundaries that separate different crystalline regions. This is physically reasonable – the disruption of local crystalline order near the grain boundaries should make the vortices in such regions more delocalized. At the relatively low temperature of 17.4K, there are very few local density peaks where the peak density is lower than $`3.5\rho _0`$ (such peaks are indicated by (red) circles in the plot). As the temperature is increased to 17.8K, the number of such peaks increases, and this trend continues as $`T`$ is increased further, as can be seen in the middle panel of Fig. 3 where the results for $`T=18.2K`$ are shown. It is also clear from these plots that the spatial regions where the low peak densities occur remain roughly unchanged as the temperature increases – as noted above, these regions are strongly correlated with the locations of the grain boundaries. These observations indicate that the vortices in the interior of the crystalline grains remain in a “solid” state as the temperature is increased to a value close to $`T_m^0`$, while those in the neighborhood of the grain boundaries begin to “melt” (get delocalized) at a lower temperature. This “inhomogeneity” of the melting process will be discussed in more detain in section III.4.

We also studied the temperature dependence of the number of topological defects present at the minima. This number increases slowly with increasing $`T`$, and then exhibits a sudden jump as liquid-like regions (characterized by low values of the local peak densities) appear near the grain boundaries at temperatures close to $`T_m^0`$. As noted in Sec. III.1 above, our method of obtaining the defect structure using the Voronoi construction becomes somewhat less reliable when liquid-like regions begin to appear. Due to this difficulty in obtaining reliable results for the defect structure and statistics at relatively high temperatures, we have not carried out a quantitative analysis of the dependence of these quantities on the temperature.

The gradual decrease in the values of the local peak densities with increasing temperature is also found in the BrG minima, but on the average the vortices remain more strongly localized at the BrG minima than at the corresponding BoG ones. The strong correlation between low values of the local peak density and the location of topological defects is found at the BrG minima too. The main difference between BrG and BoG minima is that the defects are more randomly distributed at the BrG minima: they do not line up along grain boundaries as they do at the BoG minima. For this reason, the regions of low peak density appear at fairly random locations at the BrG minima. The total number of topological defects at the BrG minima remains smaller than that at corresponding BoG minima at all temperatures. As the temperature approaches $`T_m^0`$, liquid-like regions appear in parts of the sample where the defect density is large. This is illustrated in the plots in the top panels of Figs. 3 and 4.

Thus, our investigation of the temperature dependence of the structure of the BrG and BoG minima confirms that these minima retain their topological structure as $`T`$ is increased toward $`T_m^0`$ and even beyond (provided they remain stable). The degree of localization of the vortices decreases with increasing $`T`$ and liquid-like regions with nearly delocalized vortices begin to appear at temperatures close to $`T_m^0`$. However, even at such temperatures, the BrG and BoG type minima are clearly distinguishable from one another and from the IL-type minima, as illustrated in Figs 3, 4 and 6.

### III.3 Two-step melting transition

As discussed in the preceding subsections, we find three kinds of coexisting, locally stable minima of the free energy at temperatures close to $`T_m^0`$. Two kinds (BrG and BoG) of local minima are found at temperatures substantially lower than $`T_m^0`$, and only the IL minimum is stable if $`T`$ is much higher than $`T_m^0`$. In our mean-field description, the thermodynamically stable phase at a particular temperature corresponds to the minimum with the lowest free energy. Therefore, crossings of the free energies of different kinds of minima correspond to first-order phase transitions in our description. The IL minima are expected to have the lowest free energy at high temperatures and the solid-like BrG or BoG minima should represent the globally stable phase at low temperatures. To determine how this “freezing” transition takes place, it is necessary to examine the dependence of the free energies of these different kinds of minima on the temperature $`T`$.

In Fig. 8, we have shown the results for the free energies of different minima of a sample with $`N_v=4096`$ and $`N_p=64`$. This sample is the same as the one for which results are shown in Figs. 1, 3, 4, and 7. A similar plot, showing the same general behavior for a different sample, may be found in Ref. dv3, . The BrG minimum clearly has the lowest free energy at low temperatures (while we have shown data for $`T17.8K`$, we have checked that this remains true at lower temperatures), while the IL minimum is the one with the lowest free energy at high temperatures. The interesting feature that emerges from data as that shown in this Figure is that the BoG minimum has the lowest free energy in an intermediate temperature range of small width (between 18.0K and 18.3K for this sample), indicating that the melting of the low-temperature solid phase to the high-temperature IL phase occurs in two distinct steps: the BrG phase that is the thermodynamically stable one at low temperatures undergoes a first-order transition into a BoG phase as the temperature is increased, and then this BoG phase melts into the IL phase via a second first-order transition at a slightly higher temperature. This two-step melting behavior is one of our main results.

The same qualitative behavior, indicating the occurrence of two separate first-order transitions, is found in all the 4096-vortex samples we have studied at pin concentration $`c=1/64`$ ($`N_p=64`$). The average value of the temperature interval in which the BoG phase has the lowest free energy is about 0.42K. This width exhibits fairly large sample-to-sample variations, ranging between about 0.1K in one sample to a maximum of 1.2K. As the system is cooled from the high-temperature liquid phase, it undergoes a first-order transition into a polycrystalline BoG phase. Such a transition has been observed meng  in experiments on BSCCO samples with a small concentration of columnar pins. The value of the upper (BoG to IL) transition lies between 18.2K and 18.3K in all the 4096-pin samples we have studied. These values are quite close to the first-order melting temperature $`T_m^0`$ of the same system in the absence of pins dv2 . Our results, thus, are consistent with the experimental observation kha ; ban  of a weak dependence of the freezing temperature of the vortex liquid on the pin concentration $`c`$ for small values of $`c`$. In addition, our work predicts a second first-order transition to a more ordered BrG phase at a slightly lower temperature.

For smaller samples at the same $`c=1/64`$ concentration ($`N_v=1024`$ with $`N_p=16`$) we did not always find the BoG minima: the minimization procedure that led to the BoG minima in the samples with $`N_v=4096`$ often converged to minima of the BrG type in the smaller samples. This is because the typical size of the crystalline grains at the BoG minima at this pin concentration is of the order of the sample size for 1024-vortex samples. This can be seen by comparing the middle panel of Fig. 9, where we have shown the Voronoi plot for the BoG minimum obtained for a $`N_v=1024`$ sample, with the Voronoi plot in the middle panel of Fig. 4, which is for the same $`T`$ and $`c`$, but a larger sample with 4096 vortices. One can plainly see that the domain size is comparable to the system size of the smaller samples. However, for the 1024-vortex samples where we found BoG-type minima, the behavior of the free energies was found to be very similar to that shown in Fig. 8. An example of such behavior is shown in the bottom panel of Fig, 9. The two-step melting transition found in the larger systems is found here also, indicating that this is the generic behavior. We have also carried out preliminary studies of the phase diagram for a slightly higher value of the pin concentration, $`c=1/32`$. We find the same qualitative behavior, namely the occurrence of two-step melting, with a larger difference between the two transition temperatures. We therefore conclude that the two-step melting transition we have found is a characteristic feature of our model of layered superconductors with a small concentration of columnar pins.

Since a first-order transition is not characterized by power-law behaviors of thermodynamic quantities with universal exponents, the usual finite-size scaling analysis of numerical data for continuous transitions does not apply to our work. However, we have examined the dependence of the magnitudes of various discontinuous changes at the transitions on sample size. From the results for the free energy as a function of the temperature, we have calculated the entropy jump per vortex at the first-order transitions. The entropy change at the BrG-BoG (lower temperature) transition is found to be $`0.1k_B`$ per vortex ($`k_B`$ is the Boltzmann constant), and the corresponding change at the BoG-IL (higher temperature) transition has a slightly higher value, $`0.15k_B`$ per vortex. The sum of these two entropy jumps is slightly smaller than the entropy change ($`0.29k_B`$ per vortex) at the single first-order melting transition found dv2  in the same system in the absence of pinning. All these results are physically reasonable. The difference between the degrees of order in the BrG and BoG minima is smaller than that between the BoG and IL minima, suggesting that the entropy jump at the lower transition should be smaller than that at the upper one. Also, the low-temperature BoG phase of the system in the presence of pinning is less ordered than the crystalline phase of the pure system, and the IL phase is slightly more ordered than the vortex liquid in the system without pinning (this is due to the local order of the vortices near the pinning centers which persists in the IL phase). Therefore, the net change of entropy in going from the BoG phase to the IL phase should be smaller than the entropy jump at the melting transition in the pure system. A comparison of our results for two sample sizes ($`N_v=4096`$ and 1024) does not show any appreciable dependence of the entropy changes on the sample size. We, therefore, conclude that there is no indication of any significant change in our phase diagram as the sample size is increased.

Since our mean-field treatment ignores the effects of fluctuations, it is important to address the question of whether one or both of the transitions found would become continuous if fluctuations were included. However, there are very few examples of such fluctuation-driven continuous transitions in three dimensions. In our calculations, the effects of the electromagnetic interaction among vortices on different layers are included. So, we expect our results regarding the nature of the transition to remain valid if fluctuations were included. This conclusion is supported by the result that our density functional theory provides a quantitatively a correct account dv2 ; seng ; men1  of the first-order melting of the vortex lattice in the pure system. Also, the transition in the presence of a small concentration of columnar defects is known to be first-order both experimentally kha ; ban ; meng  and from simulations non .

The appearance of a sliver of the BoG phase in our phase diagram may be qualitatively understood as arising from a competition between the elastic and pinning parts of the free energy. The elastic (free) energy plays a dominant role in the BrG minima: the effect of the randomly located pinning centers is accommodated in this structure by small displacements of the vortices from their ideal lattice positions towards the nearest pinning center. Where the occupation of a pinning center by a vortex would require a large displacement of the vortex, the pinning center is not fully occupied. For this reason, the number of pinning centers occupied by vortices at a BrG minimum is always slightly lower than the total number of pinning centers. At the corresponding BoG minima, on the other hand, the pinning centers are always occupied (except in very rare cases where two pinning centers are located very close to each other). Since the BoG minima are obtained from liquid-like initial conditions, the starting point has fully occupied pinning centers, and vortices localized in small crystalline patches in their immediate vicinity (see, for example, the bottom panel of Fig. 3 which shows such crystalline patches surrounding the pinning centers). The crystalline orientation of vortices around a vortex-binding pin site depends on the local arrangement of the pins.daf  This orientation is, in general, different in different regions of the sample. These crystalline patches grow as the temperature is reduced, and eventually meet one another at grain boundaries to form a polycrystalline BoG minimum. It is clear that the pinning centers are better accommodated in this structure than in the corresponding BrG structure. On the other hand, the creation of grain boundaries costs elastic energy. Thus, the BoG minima are expected to have lower pinning energy but higher elastic energy than the BrG ones for the same pin configuration. The elastic energy dominates over the pinning energy at low temperatures in our low pin-concentration samples. This is why the BrG phase is globally stable at low $`T`$. Both these components of the free energy decrease in magnitude as the temperature is increased. The softening of the lattice near the melting transition causes the elastic energy to decrease faster than the pinning energy. This makes the pinning component of the free energy more important than the elastic part near $`T_m^0`$, thereby making the total free energy of the BoG minima (which, as discussed above, have lower pinning energy) lower than that of the BoG minima. The BoG to IL transition at a slightly higher temperature is driven by the usual entropic mechanism.

We can substantiate this qualitative explanation by evaluating the pinning component of the free energy, as defined in Eq. (4), for the BrG and BoG minima. The results of our calculation, averaged over five samples (each with 4096 vortices and 64 pins), are shown in Fig. 10. The pinning energy is negative for both kinds of minima, as expected. The absolute value of the pinning energy of the BoG minima shows a smooth (nearly linear) decrease as the temperature is increased. This is because both the value of the parameter $`\mathrm{\Gamma }`$ (see Eq. (6)) that determines the depth of the pining potential, and the height of the local density peak at a pinning center decrease with increasing $`T`$. The plot for the BrG minima shows a similar but much slower behavior, still monotonic or nearly so: it is not clear whether the very shallow minimum near $`T`$ = 17.6K is significant. In some of the samples, the total number of occupied pining centers increases by a small amount near this temperature, thereby reducing the value of the total pinning energy. This may reflect a better accommodation of the pinning centers by the lattice, which softens as the temperature is increased. This plot clearly demonstrates that the pinning centers are better accommodated at the BoG minima. The difference in the pinning energies of the two kinds of minima decreases with increasing $`T`$, but the difference in the elastic component of the free energy decreases faster (this is clear from our results for the total free energy) due to a softening of the elastic constants with increasing temperature. The overall effect is a crossing of the free energies of the two kinds of minima near the melting transition.

Quantitatively, however, the values of the two parts of the free energy and their dependence on $`T`$ are determined by the material parameters of the superconductor and the properties of the pinning centers. Since several experimental studieskha ; ban ; meng  of the effects of irradiation-induced columnar pinning on the mixed state of BSCCO exist in the literature, we have used parameter values appropriate for this system. For other layered superconductors, the BoG phase might occur over a wider or narrower (even vanishing) range.

### III.4 Inhomogeneous melting of the vortex solid

As mentioned in section III.2, the density distribution in both BrG and BoG minima is very inhomogeneous: at temperatures close to $`T_m^0`$ there are liquid-like regions, characterized by low values of the local peak densities in parts of the sample. In this subsection, we show that this inhomogeneity leads to a spatial variation of a “local melting temperature”, defined below. As mentioned in the Introduction, such spatial variations of a local melting temperature have been deduced from measurements of the local magnetization in BSCCO samples with random pointsoi  and columnar pinningban .

Since we are dealing with time averaged densities, it is clear that the value of $`\rho _i`$ at a local peak of the density field provides a measure of the degree of localization of the vortex whose average position corresponds to the location of the density peak. A high (low) value of the local peak density implies strong (weak) localization. Smaller values of the local peak density imply mobile, liquid-like behavior. The value of the local density is, of course, equal to $`\rho _0`$ everywhere in the liquid state in the absence of pinning. The presence of pinning centers causes the local density to vary in space, but in the liquid state this variation does not lead to local peaks higher than about $`3\rho _0`$ (excluding the vortices localized at the pinning centers). A similar result was also obtained in previous studies dv2  of vortices in the presence of pinning, where it was found that values of the peak density lower than about $`3\rho _0`$ correspond to the liquid state. We, therefore, take the value $`3\rho _0`$ of the local peak density as separating solid- and liquid-like behaviors.

Using this criterion, we can determine whether a small region of the sample at a given minimum is in a locally “solid” or “liquid” state. To do this, we define a quantity $`\rho _{av}^p`$ as the average of the local peak densities in a small region containing $`100`$ vortices (see Fig. 11). The very high local density peaks representing vortices trapped at pinning centers are not included in this average. Values of $`\rho _{av}^p`$ substantially larger than $`3\rho _0`$ indicate solid-like behavior in the region under consideration, while substantially lower values of $`\rho _{av}^p`$ suggest a locally melted region. To determine how the local melting temperature varies from one part of the sample to another, we have studied the temperature dependence of $`\rho _{av}^p`$ for different regions of the sample. Typical results are shown in Fig. 11 for the pin configuration for which we have earlier shown detailed results in several Figures. In the top panel, where the positions of the local density peaks in the BoG minimum at $`T=18.0`$K and the locations of the pinning centers are shown, we have indicated three regions for which the average local peak density $`\rho _{av}^p`$ was calculated. The local density peaks in these regions are indicated by (red) plus signs (region A), (blue) crosses (region B) and (green) circles (region C), while the black dots represent the positions of the local density peaks in the other regions. The bottom panel of Fig. 11 shows the dependence of the $`\rho _{av}^p`$’s calculated for these three regions on the temperature $`T`$. At each value of $`T`$, the minimum with the lowest free energy at that temperature was used in the calculation of $`\rho _{av}^p`$.

It is clear from the plot and from similar ones for other pin configurations (one of which was shown in Ref. dv3, ) that the “local melting temperature”, operationally defined as the temperature at which $`\rho _{av}^p`$ drops from values well above $`3\rho _0`$ to values clearly below, varies from region to region. The range over which the local melting temperature varies is comparable to that found in the experiments ban ; soi . The values of the local melting temperature are strongly correlated with the local pin structure and the resulting defect structure of the vortex solid. Region B of the sample does not have any pinning center and is located across a grain boundary of the BoG minimum. As noted in section III.2, regions near grain boundaries appear to melt at temperatures that are slightly lower than the global transition temperature determined by the crossing of free energies. This is why the local melting temperature for region B is the lowest. This temperature corresponds to the transition from the BrG to the BoG state. In region A the vortices form a nearly perfect crystalline arrangement with no topological defects. The local melting temperature measured in this region is, as a consequence, close to that of the pure system. Region C contains a cluster of pinning centers, and the relatively large value of $`\rho _{av}^p`$ for this region at temperatures higher than the BoG to IL transition temperature reflects the local solid-like structure of vortices situated near pinning centers. Thus, the spatial variation of the local transition temperature is closely correlated with the “pinning landscape” of the sample.

We have found a second way of correlating the melting transition with the local density structure of the free energy minima. This second way is based on a previous studydv2  of the melting of the vortex lattice in the presence of periodic pinning, where it was found that the melting transition corresponds to the onset of percolation of liquid-like regions defined using the peak-density criterion mentioned above. We focus for this purpose on the upper transition between the BoG and IL phases. We classify vortices represented by local peaks of the density as solid-like or liquid-like, depending on whether the density at the local peak is higher or lower than $`3\rho _0`$. We then check whether the regions containing liquid-like vortices percolate across the sample. We find, as in the earlier study, that the melting transition coincides with the occurrence of percolation of the liquid-like regions. Typical results are shown in Fig. 12 where the top panel shown the locations of solid- and liquid-like local density peaks at a BoG minimum obtained at 18.2K (this minimum has the lowest free energy at this temperature), and the bottom panel shows a similar plot for the IL minimum of the same sample at 18.4K (the IL minimum is the one with the lowest free energy at this temperature, so that the melting transition in this sample occurs between 18.2K and 18.4K). The liquid-like regions do not percolate at 18.2K, whereas they do at the higher temperature of 18.4K. Thus, the melting transition in this sample also corresponds to the occurrence of percolation for the liquid-like regions. Similar behavior was found in the other samples we have studied. This observation provides a convenient way of approximately locating the transition point. This may be useful in other situations where the method of locating the transition temperature from a crossing of free energies may not be easily implementable (e.g. at higher pin concentrations where the melting transition is expected to become continuous).

## IV Conclusions

We have presented here the results of a detailed investigation of the structural and thermodynamic properties of a system of vortices in a highly anisotropic layered superconductor with a small concentration of randomly placed columnar pinning centers. Both the external magnetic field and the columnar pins are assumed to be perpendicular to the superconducting layers. Our method, based on numerical minimization of the appropriate free-energy functional, allows us to obtain very detailed information about the density distribution in the different free energy minima (which allows us to reliably identify the phases corresponding to these minima), and to map out the phase diagram of the system.

There are several salient results of our study. The first one is the occurrence of a topologically ordered BrG phase at low temperatures. While the occurrence of a low-temperature BrG phase in superconductors with a low concentration of random point pinning centers is well-established now, relatively little is known about the existence of such a phase in superconductors with random columnar pins. Our results are consistent with those of a recent numerical study non  of a similar system. It is clear from our work that the BrG minima represent a phase distinct from the polycrystalline BoG phase also found in our study. We can not conclusively rule out the possibility that free dislocations would appear at the nearly crystalline minima at length scales much longer than those considered in our numerical study. If this should happen, then the “hexatic glass” phase suggested in some earlier theoretical studies chud  would become a possible candidate for describing the BrG minima found in our study. This phase, however, would be distinct from the BoG.

Our second important result is the occurrence of a two-step melting transition: we find that the low-temperature BrG phase transforms into a polycrystalline BoG phase as the temperature is increased, and this BoG phase then melts into the high-temperature IL at a slightly higher temperature. This conclusion about the occurrence of two distinct transitions would remain valid even if the true nature of the low-temperature phase turns out to be slightly different from a Bragg glass: our work shows that the BrG and BoG minima are quite distinct from each other. The possibility of occurrence of a two-step melting transition of the vortex lattice in systems with random point pinning has been suggested earlier men . Our work provides support to this suggestion.

To our knowledge, the second (lower $`T`$) transition between the BoG and BrG phases has not been observed in experiments on layered superconductors with a small concentration of random columnar pins. This may be due to strong metastability: the BoG minimum into which the IL is expected to freeze as the temperature is decreased remains locally stable at temperatures lower than that at which its free energy crosses that of the BrG minimum, suggesting that it would be difficult to see in experiments the transition to the globally stable BrG phase. The situation here may be similar to that found in a recent experimental studyshobo  of a low-$`T_c`$ superconductor with weak point pinning which is expected to exhibit a BrG phase at low $`T`$. It is, however, found in the experiment that the vortex solid obtained by cooling the sample in the presence of the external magnetic field has a polycrystalline structure, indicating that the metastability of this disordered state prevents the system from reaching the more ordered (BrG) equilibrium state at low temperatures. In our numerical work, the BrG minima were obtained by performing the free-energy minimization from an initial configuration with crystalline order. It is not clear how a similar procedure can be adopted in experimental studies.

A low-temperature BrG phase has been observed in a recent simulationnon  of a model of HTSC with a small concentration of columnar pinning centers. However, this simulation finds a single first-order transition between the BrG and IL phases at low pin concentrations: the small intermediate region of BoG phase found in our study is not observed. This is probably due to the smallness of the system sizes ($``$ 100 vortex lines) used in the simulation. As discussed in Sec. III.3 above, the polycrystalline BoG minima are found only if the sample size is larger than the typical size of the crystalline domains: only the BrG minimum is found at low temperatures if this condition is not satisfied. Since the crystalline domains become large at small pin concentrations, it is quite likely that a simulation with small system sizes and low pin concentration would not see the BoG phase. Another possibility is that a narrow “two-phase” region found in Ref. non,  near the BrG melting transition actually corresponds to the intermediate BoG phase found in our study. A third possibility is that the intermediate BoG phase, which arises due to a competition between the elastic and pinning components of the free energy (see Sec. III.3 where it is shown that the pinning energy favors the BoG phase), does not appear in the simulation of Ref. non,  because the pins are assumed to be ”weak” in that work.

Our study also illustrates the spatial inhomogeneity of the melting process in the presence of disorder: we have found that a “local” melting temperature defined using a criterion based on the degree of localization of the vortices shows considerable spatial variation. This spatial variation is strongly correlated with the ”pinning landscape” associated with the random spatial location of the pinning centers. Our results about the spatial inhomogeneity of the local melting temperature are consistent with those of recent experiments ban ; soi  on superconductors with pinning disorder.

These results establish the usefulness of our numerical method in dealing with the problem of vortex matter in the presence of random pinning. It would be interesting to examine the extent to which our results depend on the values of the parameters appearing in the free energy functional. It would also be clearly useful to carry out similar studies of other related problems, such as the complete phase diagram of systems with random columnar pins in the $`Tc`$ plane, and the phase behavior of systems with random point pinning. Some of these investigations are currently in progress.

###### Acknowledgements.

We thank G. I. Menon for several useful discussions.


text:
# On the use of X-rays to determine dynamical properties of Elliptical galaxies

## 1 Introduction

Since its discovery in the late 1970’s (Fabbiano 1989) the presence of an extended and X-ray emitting interstellar medium (ISM) in early type galaxies raised the hope of having a new powerful and reliable tool to determine the dynamical properties of these galaxies, in particular their mass profile and total mass out to large radii ($`10`$ optical effective radii $`R_\mathrm{e}`$). In the simplest and commonly adopted approach, these determinations assume a spherical mass distribution. Then, the most general expression for the momentum equation of the gas is:

$$\frac{k}{\mu m_p}\frac{1}{\rho }\frac{d\rho T}{dr}+𝒱=g,$$
(1)
where $`T`$ and $`\rho `$ are the radial profiles of the ISM temperature and density, $`m_p`$ is the proton mass, $`\mu `$ is the mean molecular weight and $`g=GM(r)/r^2`$ is the total gravitational field. $`𝒱`$ includes all the non hydrostatic terms and, for simplicity, we do not consider other pressure terms such as the magnetic one (see however Sect. 4). Thus, in hydrostatic equilibrium, $`T(r)`$ and $`\rho (r)`$ – that in principle can be inferred from the X-ray surface brightness profile through deprojection, using also spectral information (e.g., Forman et al. 1985, Trinchieri et al. 1986) – are linked to the total mass profile by the well known formula (Fabricant, Lecar & Gorenstein 1980)

$$M(r)=r\frac{kT}{G\mu m_p}\left(\frac{d\mathrm{log}\rho }{d\mathrm{log}r}+\frac{d\mathrm{log}T}{d\mathrm{log}r}\right).$$
(2)
Forman et al. (1985) made the first total mass determinations using this approach, for a sample of early type galaxies with Einstein data extending out to radii of $`50`$ kpc. They obtained large mass values implying that dark matter should be present in a large quantity ($`M/L_B25100M_{}/L_{B,}`$). However, uncertainties in the $`\rho `$ and $`T`$ profiles were so large that mass values could be even as low as those derived from optical measurements (Trinchieri, Fabbiano & Canizares 1986; Fabbiano 1989), that typically give $`M/L_B12M_{}/L_{B,}`$ within $`R_\mathrm{e}`$ (Saglia, Bertin & Stiavelli 1992).

Mass profiles can in fact be derived also from the observed stellar velocity dispersion profiles $`\sigma _p(R)`$; these extend out to $`12R_\mathrm{e}`$ in few best cases (e.g., Kronawitter et al. 2000). In this context, a common approach is to use the Jeans equations; for a spherical model described by a phase-space distribution function $`f=f(E,L^2)`$, the relevant equation is:

$$\frac{1}{\rho _{}}\frac{d\rho _{}\sigma _r^2}{dr}+\frac{2\beta \sigma _r^2}{r}=g$$
(3)
(e.g., Binney & Tremaine 1987). Here $`\rho _{}`$ is the stellar density profile and $`\beta (r)=1\sigma _t^2(r)/2\sigma _r^2(r)`$, where $`\sigma _r`$ is the radial stellar velocity dispersion and $`\sigma _t^2=\sigma _\phi ^2+\sigma _\vartheta ^2`$ is the total tangential stellar velocity dispersion. $`\beta (r)`$ describes the radial trend of anisotropy in the stellar velocity dispersion tensor, which is radially anisotropic where $`0<\beta <1`$, tangentially anisotropic where $`\beta <0`$ and isotropic where $`\beta =0`$. For assumed visible and dark density profiles and a specified radial trend of $`\beta (r)`$, eq. (3) is solved to derive $`\sigma _r^2(r)`$; then by projection one obtains $`\sigma _p(R)`$ to be compared with the observed one. As well known, in absence of a physically motivated choice for $`\beta (r)`$, this method suffers from a “dark matter–$`\beta `$ degeneracy”: a sufficient amount of orbital anisotropy can mimic the presence of a massive dark matter halo (e.g., Tonry 1983). Observationally, $`\beta (r)`$ can be derived from the deviation of the line profiles from a Gaussian (e.g., Gerhard et al. 1998), but this process requires high quality data and care in analysis. In fact, the use of the hot ISM to derive dynamical quantities was very welcomed also because it is exempt from this degeneracy problem111Note that in other approaches (e.g., Bertin et al. 1994), the anisotropy profile is not a function to be specified by independent arguments, but derives from the physical picture on which the model is based..

In hydrostatic equilibrium, the two conditions (1) (with $`𝒱=0`$) and (3), holding for the gas and the stars lying in the same total gravitational field $`g`$, allow us to determine $`M(r)`$ and $`\beta (r)`$ uniquely, when coupled with the knowledge of the four observables $`T(r)`$, $`\rho (r)`$, $`\rho _{}(r)`$ and $`\sigma _p(R)`$. For example, by eliminating $`g`$ from (1) and (3) one obtains a relation linking gas and stellar properties only (Fabricant et al. 1989):

$$\frac{k}{\mu m_p\rho }\frac{d\rho T}{dr}=\frac{1}{\rho _{}}\frac{d\rho _{}\sigma _r^2}{dr}+\frac{2\beta \sigma _r^2}{r}.$$
(4)
Then, one can use the equation above to derive $`\beta (r)`$ using information coming from X-ray observations, as attempted recently by Mathews & Brighenti (2003a; hereafter MB03).

The lack of reliable ISM density and temperature profiles has been a severe limitation to the use of X-rays as a probe of the galactic potential and of stellar dynamical quantities. This was first overcome partially thanks to $`ROSAT`$ observations. Currently, quite accurate $`T`$ and $`\rho `$ profiles are being derived from $`Chandra`$ and $`XMMNewton`$ observations (e.g., Sun et al. 2003). In addition, in a few cases $`Chandra`$ has already shown signs of departure from the hydrostatic equilibrium in the hot gas coronae (e.g., Jones et al. 2002). Therefore a timely question is how much reliable the X-ray estimate of dynamical quantities can be.

In this paper we address the problem of the accuracy of the X-ray based total mass $`M(r)`$ and $`\beta (r)`$ determination when the hydrostatic equilibrium hypothesis is assumed but it does not hold. Our analysis relies upon the use of the relevant hydrodynamical equations, while we do not discuss technical problems related with the extraction of physical quantities from the observations, such as the well known ill-posed problem of Abel inversion for observed data (e.g., Tikhonov & Goncharsky 1987).

## 2 Estimates of mass and orbital anisotropy when $`𝒱0`$

In the most general case ($`𝒱0`$) the physically correct analog of eq. (4) is

$$\frac{k}{\mu m_p\rho }\frac{d\rho T}{dr}+𝒱=\frac{1}{\rho _{}}\frac{d\rho _{}\sigma _r^2}{dr}+\frac{2\beta \sigma _r^2}{r}.$$
(5)
Instead, by erroneously assuming hydrostatic equilibrium, one uses an equation involving $`\beta ^{est}`$:

$$\frac{k}{\mu m_p\rho }\frac{d\rho T}{dr}=\frac{1}{\rho _{}}\frac{d\rho _{}\sigma _r^2}{dr}+\frac{2\beta ^{est}\sigma _r^2}{r}.$$
(6)
Therefore the difference between the estimated and the true anisotropy profile is given by

$$\beta ^{est}\beta =\frac{r}{2\sigma _r^2}𝒱.$$
(7)
Following the same reasoning, the equation that should be used to derive the total mass profile $`M(r)`$ is eq. (1), while using it with $`𝒱=0`$, one derives an estimated mass $`M^{est}`$ given by:

$$M^{est}=M+\frac{r^2}{G}𝒱.$$
(8)
Note that, although $`M^{est}`$ and $`\beta ^{est}`$ derive from two independent procedures, the deviation in the mass estimate is linked to that in the $`\beta `$ estimate by a remarkably simple relation:

$$\frac{M^{est}(r)}{M(r)}1=\frac{2r\sigma _r^2(r)[\beta ^{est}(r)\beta (r)]}{GM(r)}.$$
(9)
Thus, for example, $`M^{est}<M`$ wherever $`\beta ^{est}>\beta `$.

In general, the modulating factor $`2r\sigma _r^2/GM`$ in eq. (9) depends on the exact slope of the stellar density profile, on the presence of dark matter and on $`\beta (r)`$ \[via $`\sigma _r(r)`$\]. The radial behavior of this factor for a few representative cases is given in Fig. 1. Within $`0.2R_\mathrm{e}`$ the factor is very sensitive to the shape of the density profile, as shown by the cases of the de Vaucouleurs (1948), Jaffe (1983) and Hernquist (1990) isotropic models (respectively the solid, long-dashed and dotted lines). In addition, Fig. 1 also shows that this dependence is not trivially linked to the radial behavior of the $`\rho _{}`$ profile, since $`\rho _{}(r)r^{3/4},r^2`$ and $`r^1`$ respectively for the de Vaucouleurs (Ciotti 1991), Jaffe and Hernquist models. The inner shape of the factor can be derived, however, from the asymptotic analysis discussed by Bertin et al. (2002).

The short-dashed line in Fig. 1 refers to an isotropic stellar de Vaucouleurs model coupled to a dark matter halo that is described by an Hernquist profile and is 5 times as massive as the stellar component; the mass enclosed within $`R_\mathrm{e}`$ is the same for the two components (which implies a core radius for the Hernquist model of $`2.47R_\mathrm{e}`$). Compared to the corresponding pure stellar model (solid line) there is a small increase in the modulating factor.

When $`\beta >0`$ the factor is larger than in the isotropic case; for example, when $`\beta `$ is radially constant222Any consideration on the physical consistency of this assumption is left apart here; see, e.g., Camm 1952., as adopted in some modelings (see, e.g., Sect. 3), and the density distribution is $`\rho _{}r^\gamma `$, the factor is $`1/(\gamma 10.5\beta )`$. In Fig. 1 the dot-dashed line refers to a de Vaucouleurs profile with $`\beta (r)`$ radially varying according to the Osipkov-Merritt parameterization (Osipkov 1979, Merritt 1985). In the central region, the effect of such a $`\beta (r)=r^2/(r^2+r_a^2)`$ (where $`r_a`$ is the anisotropy radius) is a small increase of the modulating factor with respect to the isotropic case (solid line). In the external regions an idea of the full range of the factor values derives by comparing the cases of isotropic and anisotropic orbital distributions for a de Vaucouleurs profile (solid and dot-dashed lines respectively). In fact, for this profile, the choice of $`r_aR_\mathrm{e}`$ adopted for Fig. 1 gives the maximum allowed level of anisotropy, based on stability arguments (Ciotti & Lanzoni 1997; see also Sect. 3).

In summary, outside the effective radius the modulating factor is everywhere $`<0.7`$; for the “maximally” anisotropic case, it is everywhere $`<1`$. Within $`R_\mathrm{e}`$ it is again everywhere $`\begin{array}{c}<\hfill \\ \hfill \end{array}1`$, except for the Hernquist model. Therefore for most models “errors” on $`M`$ (i.e., $`M^{est}/M1`$) are $`\begin{array}{c}<\hfill \\ \hfill \end{array}`$ than the $`\beta `$ deviations (i.e., $`\beta ^{est}\beta `$).

In the next sections we estimate the difference $`\beta ^{est}\beta `$ in eq. (7), first qualitatively (Sect. 2.1) and then quantitatively by using numerical models of hot gas flows for X-ray bright galaxies (Sect. 2.2).

### 2.1 Qualitative analysis

We assume that the hot ISM filling the galactic potential well has velocity, density and temperature $`u`$, $`\rho `$ and $`T`$. Usually, in the case of galactic flows, a mass source rate $`\alpha \rho _{}(r)`$ is present in the $`𝒱`$ term, due to the stellar mass losses (e.g., Ciotti et al. 1991). For an inviscid fluid, the $`𝒱`$ term is then written as

$$𝒱=\frac{u}{t}+\frac{1}{2}\frac{du^2}{dr}+\frac{\alpha \rho _{}u}{\rho }.$$
(10)
We consider first eq. (7) in the simplest case, namely the stationary case with $`\alpha =0`$; then

$$\beta ^{est}\beta =\frac{r}{4\sigma _r^2}\frac{du^2}{dr}.$$
(11)
This shows that when the velocity decreases with increasing radius in magnitude, i.e., when $`du^2/dr<0`$, then $`\beta ^{est}(r)>\beta (r)`$. This is the case of the central regions in cooling flows. If, instead, the velocity increases in modulus with increasing radius, then $`\beta (r)`$ is underestimated. This could be the case of galactic winds. Note that fairly large differences $`\beta ^{est}\beta `$ are expected when $`u^2`$ varies by $`\sigma _r^2`$ over a characteristic length of the order of $`r`$.

When $`\alpha >0`$, these results are reinforced in two cases: when the gas is inflowing ($`u<0`$) with a velocity decreasing with radius ($`du^2/dr<0`$), then even more radially anisotropic orbits are deduced; and when the gas is outflowing ($`u>0`$) with a velocity increasing with radius ($`du^2/dr>0`$), then even more tangentially anisotropic orbits are deduced. In the other cases the $`\alpha \rho _{}/\rho `$ term may have a competing effect on the $`du/dr`$ term that depends on the magnitude of the source function $`\alpha `$.

All the results above can be converted straightforwardly into results concerning mass estimates by using eq. (9).

The X-ray brightest galaxies, for which better quality data are available and dynamical studies can be attempted, contain large quantities of hot gas and are supposed to host inflows for which $`du^2/dr<0`$. In their inner regions (well within $`R_\mathrm{e}`$) the gas can reach high velocities, up to the local sound speed, while $`u0`$ outside $`R_\mathrm{e}`$ (e.g., Ciotti et al. 1991). Therefore, it is expected that more radially anisotropic orbits than true are deduced from X-ray observations of the central regions; $`\beta ^{est}`$ derived outside $`R_\mathrm{e}`$ should be reliable. More quantitative estimates follow below.

### 2.2 Quantitative results

Here we study the effect of the full $`𝒱`$ term given in eq. (10). For this investigation, two gas dynamical models are used, resulting from the numerical simulations of time evolving galactic gas flows by Pellegrini & Ciotti (1998) and representative of hot gas rich galaxies. In principle, one should use a numerically computed right hand side of eq. (10). However, in order to reduce the numerical noise due to spatial and temporal derivatives over a finite grid, the $`𝒱`$ term is obtained from eq. (1), so that:

$$\beta ^{est}\beta =\frac{r}{2\sigma _r^2}\left(g+\frac{k}{\mu m_p\rho }\frac{d\rho T}{dr}\right).$$
(12)
Similarly:

$$\frac{M^{est}}{M}=\frac{k}{g\mu m_p\rho }\frac{d\rho T}{dr}.$$
(13)
In fact, for the galaxy models we know $`g`$ and $`\beta `$. As representative cases we consider 1) a bright inflowing galaxy ($`L_B=7\times 10^{10}L_{}`$, $`R_\mathrm{e}=9.0`$ kpc, central $`\sigma =287`$ km s-1 and a dark halo nine times as massive as the stellar component); 2) a less massive galaxy in which the inflow is confined within 0.1 $`R_\mathrm{e}`$ ($`L_B=2.2\times 10^{10}L_{}`$, $`R_\mathrm{e}=4.4`$ kpc, central $`\sigma =210`$ km s-1 and a dark halo twice as massive as the stellar component). The underlying galaxy model is isotropic ($`\beta =0`$).

Fig. 2 shows in the upper panel the radial trend of the gas velocity at the end of the simulations (corresponding to an elapsed time of 15 Gyrs). The panel below shows whether any term in the right hand side of eq. (10) gives a dominant contribution to $`𝒱`$. Since the gas flows have reached a quasi steady state after 15 Gyrs, the $`u/t`$ term is negligible, and only the radial trend of the advective/source term ratio is plotted. For the largely inflowing galaxy the advective term is dominating in the central region and becomes comparable to the source term at $`R_\mathrm{e}`$, while for the partial wind model the two terms keep comparable. For both models the two terms produce a deviation in the same sense, since they are both negative in the inflowing region and positive in the outflowing region.

In the lowest two panels, Fig. 2 shows what $`\beta ^{est}`$ and $`M^{est}`$ would be derived from X-ray observations. Overall, the trends of $`\beta ^{est}`$ agree with the qualitative analysis of Sect. 2.1 and, as expected, $`\beta ^{est}`$ is $`>0`$ for the central gas inflowing region. For the galaxy hosting a global inflow, decidedly more radially anisotropic orbits than true would be deduced by neglecting the gas velocity; correspondingly, $`M^{est}<M`$ would be derived for the central galactic region \[eq. (9) and Fig. 2\]. Outside the inflow region, $`\beta ^{est}`$ is small and $`M^{est}M`$, consistently with the predictions for low gas velocities.

## 3 The case of NGC4472

The X-ray bright Virgo galaxy NGC4472 has been chosen since the earliest times for a detailed investigation of its dynamical properties using X-ray data (Thomas 1986; Trinchieri et al. 1986; Loewenstein 1992; Bertin, Pignatelli & Saglia 1993; Bertin & Toniazzo 1995; Irwin & Sarazin 1996). The study based on the most recent and extensive data set is that of MB03 (see also Mathews & Brighenti 2003b). MB03 employed eq. (4) to relate $`\beta (r)`$ and $`T(r)`$. They used $`\rho (r)`$ derived from Abel-inverted $`Chandra`$, $`ROSAT`$ and $`Einstein`$ surface brightness profiles, $`T(r)`$ derived from deprojection of $`ROSAT`$ and $`Chandra`$ data, the observed $`\sigma _p(R)`$ from Bender et al. (1994) and Fried & Illingworth (1994) and $`\rho _{}(r)`$ resulting from fitting the observed brightness profile with a de Vaucouleurs law (but in the central regions the density profile was assumed to steepen as $`r^{0.9}`$).

Adopting a radially constant $`\beta ^{est}`$ for the modeling, the value $`\beta ^{est}=0.71\pm 0.15`$ was derived, while the optically determined $`\beta (r)`$ within $`R_\mathrm{e}`$ is everywhere $`\begin{array}{c}<\hfill \\ \hfill \end{array}0.34`$ (Kronawitter et al. 2000). That is, $`\beta ^{est}`$ is about twice as large as the maximum $`\beta `$ found from optical data. Conversely, assuming the optically determined $`\beta (r)`$ of Kronawitter et al. (2000), $`T(r)`$ higher than observed is obtained. Commenting this result, MB03 state that “This strong preference for radial stellar orbits must be understood in terms of the formation history of massive elliptical galaxies. Conversely, if the smaller, optically determined $`\beta `$ is indeed correct, we are led to the important conclusion that the temperature profile $`T(r)`$ of the hot ISM in NGC 4472 must differ from that indicated by X-ray observations, or that the hot gas is not in hydrostatic equilibrium.

Here we argue that, if problems with the several steps involved in the modeling of the data can be excluded, the last possibility of lack of hydrostatic equilibrium is strongly favoured. In fact a simple stability argument rules out as unphysical such a high X-ray–derived $`\beta ^{est}`$ value. In case of a constant anisotropy parameter, the radial orbit instability indicator (Polyachenko & Shukhman 1981, Fridman & Polyachenko 1984) turns out to be

$$\frac{2K_r}{K_\vartheta +K_\phi }=\frac{1}{1\beta ^{est}}3.45,$$
(14)
where $`K_r`$ and $`K_\vartheta +K_\phi `$ are respectively the total radial and tangential kinetic energies. A galaxy is unstable due to radial orbit instability when this indicator is $`\begin{array}{c}>\hfill \\ \hfill \end{array}1.7\pm 0.25`$; therefore the value calculated in eq. (13) is well outside the stability regime. Note that even when inserting the maximum value of $`\beta (r)=0.34`$ derived by Kronawitter et al. (2000) in eq. (13), the indicator is $`1.5`$, safely in the stability regime. The criterion above has been proved to be obeyed essentially by all galaxy models to which it has been applied, even including massive dark halos (e.g., Stiavelli & Sparke 1991; Nipoti, Londrillo & Ciotti 2002; Bertin & Trenti 2003; see also Saha 1991).

Additional considerations argue in favour of the hypothesis of lack of hydrostatic equilibrium to explain the inconsistencies between X-ray and optically-derived orbital anisotropy. In fact, there is also a deviation between the X-ray and the optically–derived mass for NGC4472, and the two discrepancies (concerning $`\beta `$ and $`M`$) are related as predicted by eq. (9), as shown in the following. Having discarded the X-ray–derived $`\beta ^{est}`$ as unphysical, we assume that the optical estimate of $`\beta `$ is the correct one. Then, eq. (9) predicts the X-ray–derived $`M^{est}`$ to underestimate the true mass $`M`$. More specifically, the expected value of $`M^{est}/M`$ can be derived from the results for the de Vaucouleurs models in Fig. 1. For a range of values of the modulating factor within $`R_\mathrm{e}`$ of $`0.51`$ and for $`\beta ^{est}\beta =0.4`$, then $`M^{est}/M=0.80.6`$ is given by eq. (9). What $`M^{est}`$ derives from the X-ray data? The best fit total mass profile $`M^{est}(r)`$ determined from eq. (2), with the same $`T(r)`$ and $`\rho (r)`$ used by MB03, is a one-component de Vaucouleurs mass profile out to $`r=R_\mathrm{e}`$, with a total constant mass-to-light ratio of $`M^{est}/L_B=7`$ (Mathews & Brighenti 2003b). $`M/L_B`$ estimates were derived also from optical studies based on high quality photometrical and kinematical data extending out to $`0.5R_\mathrm{e}`$. The $`M/L_B`$ determined from one-component constant mass-to-light ratio axisymmetric dynamical models is 9.2, for the central galactic region, when rescaled to the distance adopted by MB03 (van der Marel 1991). Another analysis used spherical self-consistent two-component models (Saglia et al. 1992). The models selected at a 3 sigma confidence level have $`M/L_B=10.114.6`$ within $`R_\mathrm{e}`$, when rescaled to the MB03 distance333These selected models also have a velocity dispersion profile consistent with the data for the globular cluster system out to $`2.5R_\mathrm{e}`$.. These optically determined $`M`$ values, if taken as true values, imply that $`M^{est}`$ underestimates the true mass $`M`$. More precisely, $`M^{est}/M=0.76`$ and $`0.690.48`$ when adopting respectively the van der Marel’s and the Saglia et al.’s results, in remarkable agreement with the predictions above ($`M^{est}/M=0.80.6`$). Therefore, the hypotheses of i) correct optical estimates of $`\beta `$ and $`M`$ and ii) departure from hydrostatic equilibrium of the hot ISM, give a consistent explanation for the discrepancies of the X-ray–derived $`\beta ^{est}`$ and $`M^{est}`$ values. In particular, these hypotheses account for the fact that $`\beta ^{est}`$ is an overestimate and $`M^{est}`$ is an underestimate and also for the size of the discrepancies with the true $`\beta `$ and $`M`$.

## 4 Discussion and Conclusions

We have investigated the effect that deviations from hydrostatic equilibrium have on the X-ray estimates of the total mass $`M`$ and the orbital anisotropy parameter $`\beta `$ via the “standard” eqs. (2) and (4). We have found that:

$``$ in the most general case, the deviations of $`\beta ^{est}`$ and $`M^{est}`$ from the true values are linked by a remarkably simple relation. In particular, the mass is underestimated if $`\beta `$ is overestimated. Also, for the most common density and orbital anisotropy distributions, with or without a massive dark halo, the $`\beta ^{est}`$ deviations are larger than the “errors” on $`M`$ at all radii $`\begin{array}{c}>\hfill \\ \hfill \end{array}0.2R_\mathrm{e}`$.

$``$ in the stationary case without mass/momentum sources/sinks, more radially anisotropic orbits than true are deduced for the central regions of gas inflowing galaxies. Instead, $`\beta (r)`$ is underestimated for outflowing regions where the velocity increases with radius, as in galactic winds. When mass sources are present, the results above are reinforced in the inflow and in the wind cases. In the other cases the mass source term may have a competing effect.

$``$ the full time dependent problem is investigated with the aid of numerical simulations reproducing gas rich galaxies. The results are consistent with the simple qualitative predictions. $`\beta ^{est}`$ decidedly overestimates the true value (and $`M^{est}`$ underestimates it) if derived from the central region where gas is inflowing with a large velocity, while outside it $`\beta ^{est}\beta `$ and $`M^{est}M`$.

$``$ $`\beta ^{est}`$ and $`M^{est}`$ derived for the most thoroughly investigated bright elliptical galaxy (NGC4472), for which the best X-ray data are available, are considered as a representative application of our results. We show that the high $`\beta ^{est}`$ recently derived is unphysical. We also note a discrepancy in the X-ray and optically derived masses within $`R_\mathrm{e}`$. Assuming as true $`\beta `$ and $`M`$ the optically derived values, the $`\beta `$-overestimate implies (when applying our results for the non-hydrostatic equilibrium) an $`M`$-underestimate in remarkable agreement with $`M^{est}M`$. Of course, one cannot exclude that a different modeling of the X-ray and optical data would be able to produce a $`\beta (r)`$ profile consistent with the observed one, in the hydrostatic equilibrium hypothesis. However, this would reveal a degeneracy in the procedure and undermine the possibility to uniquely derive dynamical information from the X-ray data.

$``$ note that eqs. (7–9) describe the deviations of $`\beta ^{est}`$ and $`M^{est}`$ from the true values even when sources of pressure additional to the thermal one are effective (such as magnetic and radiation pressure). In fact, if these are accounted for explicitely in eqs.(2) and (4), then eqs. (7–9) remain unaltered. If these are not accounted for, again eq. (9) remains unaltered and eqs.(7–8) hold provided that the additional pressure terms are considered as included in the $`𝒱`$ term. The deviations of $`\beta ^{est}`$ and $`M^{est}`$ produced by these kinds of pressure may be confined to the central galactic regions, where magnetic fields and AGN-related phenomena more likely are concentrated.

Additional explanations can be figured out for the deviations of X-ray–derived $`\beta ^{est}`$ and $`M^{est}`$ from the true values. These include deviations from spherical symmetry, projection/deprojection effects not adequately taken into account, presence of a multiphase ISM at each radius, amplification of uncertainties in the process of extracting $`T(r)`$, $`\rho (r)`$, $`\rho _{}(r)`$ and $`\sigma _r(r)`$ from the data. However, it is yet to be investigated the possibility for these explanations to account at the same time for deviations in $`\beta `$ and $`M`$, and by precisely the observed amount, as demonstrated here for the hypothesis of lack of hydrostatic equilibrium in the case of NGC4472.

Our analysis casts doubts on the possibility to use X-ray information and the hydrostatic equilibrium assumption as robust tools to determine accurately the dynamical properties of X-ray emitting galaxies, at least within $`R_\mathrm{e}`$. Another warning to take with care X-ray mass estimates had already been made when it appeared that mass models constrained by accurate and extended stellar dynamical data, and including different amounts of dark matter, all lead to similar temperature and X-ray brightness profiles (Bertin et al. 1993), within the uncertainties in the X-ray measurements. In this sense, the situation might be improving with the new $`Chandra`$ and $`XMM`$ data.

In conclusion, instead of obtaining robust dynamical estimates from hot gas-based studies, it should perhaps be expected that stellar dynamical studies (possibly coupled with gravitational lensing studies) help understand the dynamical status of the hot gas.

## Acknowledgments

We would like to thank G. Bertin for useful discussions and the referee, M. Loewenstein, for comments that improved the paper.