text:
# Phase diagram of a hard-sphere system in a quenched random potential: a numerical study

## I Introduction

The equilibrium phase diagram of a classical system of interacting particles in a quenched, random, pinning potential is an important subject on which much effort is currently being spent . There are several experimentally studied systems, such as vortices in the mixed phase of high-Tc superconductors , fluids confined in porous media , magnetic bubble arrays , and Wigner crystals , which provide physical realizations of a collection of interacting classical objects in the presence of an external, time-independent, random potential. In the absence of such a potential, systems of this kind are expected to crystallize at sufficiently low temperatures. Several years ago, Larkin  showed that the presence of arbitrarily small amount of random pinning disorder destroys long-range translational order in all physical dimensions $`d<4`$. However, recent theoretical studies  suggest that weak disorder distorts the crystalline state only slightly, leading to a phase with perfect topological order and logarithmic fluctuations of the relevant displacement field. This phase, with quasi-long-range translational order and power-law Bragg peaks in the structure factor, is called a “Bragg glass” . The transition point between a Bragg glass and the high-temperature liquid phase is likely to be shifted with increasing disorder, but the transition is believed to remain first order as long as the disorder is weak. A question of obvious interest is how this transition temperature and the nature of the transition depend on the strength of the random potential.

As the relative strength of the disorder is increased, so that the week-disorder situation described above no longer applies, the Bragg glass phase is expected to undergo a transition to a topologically disordered amorphous phase with only short-range translational correlations. It is not yet clear whether this phase is thermodynamically distinct from the high-temperature liquid. An interesting possibility is that it is analogous to the glassy phase obtained by supercooling a liquid below the structural glass transition temperature in the absence of external quenched disorder . If this is so, then the phase diagram of such systems would contain three phases: a Bragg glass phase obtained at low temperature and weak disorder, an amorphous (without quasi-long-range translational order) glassy phase at low temperatures and strong disorder, and a weakly inhomogeneous (because of the presence of the external random potential) liquid phase at high temperatures. The glassy phase would be a thermodynamically stable one in these systems. This is different from the situation in the absence of external disorder where the crystalline state is the true equilibrium state near the structural glass transition and both the supercooled liquid and the glass are metastable. In other words, the presence of external disorder may lead to the possibility of occurrence of a true, thermodynamically stable, glassy phase.

The phase diagram  in the temperature ($`T`$) – magnetic field ($`H`$) plane of layered, highly anisotropic, type-II superconductors such as $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_8`$ in a magnetic field perpendicular to the layers is a credible candidate to exhibit these three phases. For a wide range of values of $`H`$, the flux lines in these materials may be regarded as columns of interacting “pancake” vortices  residing on the layers, and the properties of the mixed phase may be described in terms of the classical statistical mechanics of these point-like objects. In these compounds, at low enough fields, a flux-lattice melting transition separates a nearly crystalline state of the flux lines from a disordered “vortex liquid” state. The first-order character of this transition has been carefully documented . When the magnetic field is increased, the transition becomes continuous , and the nearly crystalline state appears to be replaced by an amorphous state called “vortex glass”  which is endowed with glassy properties such as non-ohmic current-voltage characteristics . It is generally assumed  that the vortex glass phase owes its existence to the presence of point-like pinning disorder. Observation of Bragg peaks in neutron scattering experiments  confirms that the phase at low $`H`$ and $`T`$ is a Bragg glass. As the effective strength of the disorder is increased, either indirectly by increasing $`H`$ (which is believed to increase  the effective strength of the disorder), or directly by increasing the amount of defects in the sample , the Bragg glass phase changes over to the vortex glass. The latter is separated from the liquid by a continuous transition . This phase diagram, thus, suggests that the first-order liquid-to-crystal transition in a three-dimensional (3d) system of point-like objects may be driven by the pinning disorder into a continuous liquid-to-glass transition.

The formation of a glassy phase at strong disorder was investigated recently  in an analytic study of the phase diagram of a system of hard spheres in a random pinning potential of arbitrary strength. This work used a combination of two “mean-field”- type approaches based on the “replicated liquid formalism” : the replica method  was used for treating the effects of quenched disorder, and the hypernetted chain approximation  to calculate the equilibrium correlation functions in the liquid in the presence of the pinning potential. These correlation functions were then the input in a replicated density functional  of the Ramakrishnan-Yussouff (RY) form  from which the location of the freezing transition of the liquid into a nearly crystalline (Bragg glass) phase was obtained. The possibility of a liquid-to-glass transition was investigated using the phenomenological approach of Mézard and Parisi . The resulting phase diagram in the density – disorder plane (the density, rather than the temperature, is the appropriate control parameter for a hard-sphere system) shows three thermodynamic phases: a nearly crystalline Bragg glass, an amorphous glassy phase, and a low-density liquid. It is consistent with the expectation (from earlier work  and the Lindemann criterion ) that the density at which the Bragg glass to liquid transition occurs should move to higher values as the strength of the disorder is increased. The first-order crystallization transition is replaced by a continuous glass transition as the disorder strength is increased above a threshold value. This phase diagram is thus qualitatively similar to that proposed for some layered type-II superconductors if, as noted above, the density is replaced by the temperature $`T`$ and the disorder strength by the magnetic field $`H`$.

In the present paper, we report the results of a numerical investigation of the phase diagram of the same system, i.e. a hard-sphere fluid in the presence of a random pinning potential with short-range spatial correlations. Our work involves the use of direct numerical minimization to study the effects of the presence of an external random potential on the minima of a discretized version of the RY free energy functional for the hard-sphere system. It is known  that in the absence of external disorder, this model free energy functional exhibits, at sufficienly high densities, a large number of “glassy” local minima characterized by inhomogeneous but aperiodic density distributions. In addition, a global minimum corresponding to the crystalline solid is also found at high densities if the sample size and the discretization scale are commensurate with the crystal structure. If they are incommensurate, only the glassy minima are present. We have carried out extensive numerical investigations of the resulting free-energy landscape  in the absence of disorder. In the present study, and with the physical situations described above in mind, we develop similar numerical methods to find the location and structure of the local minima of the same model free energy with the addition of the presence of a time-independent, random, one-body potential.

Using these numerical methods we investigate how the uniform liquid, crystalline solid and glassy minima of the free energy in the absence of the random potential evolve as the strength of the random potential is increased from zero. We also examine the dependence of the free energy of these minima, and of the density structure (as given by the two-point density correlation function) of the system at these minima, on the strength of the disorder. In this picture, a transition from one phase to another is signalled by the crossing of the free energies of the corresponding minima of the free energy. By monitoring where these crossings occur as the density and the strength of the disorder are varied, we are able to map out the phase diagram in the density – disorder plane. This phase diagram is qualitatively very similar to the one obtained in the analytic work. For weak disorder we find, in the commensurate case as described above where a crystalline minimum exists, a first-order liquid-to-crystal transition that moves to higher density as the strength of the disorder is increased. In the metastable “supercompressed” regime (i.e. at a density higher than the value at which equilibrium crystallization takes place for the commensurate case), we find in all cases a liquid-to-glass transition. The density at which this transition occurs decreases (very slowly for the largest systems studied, which are incommensurate, and more rapidly for the smaller, commensurate systems) as the strength of the disorder is increased. The nature of this glass transition depends on the strength of the disorder: it is first order when the disorder is weak, but it changes to second order beyond a certain critical value of the disorder strength. For the commensurate case, the crystallization line crosses the glass transition line at or very near the same critical value of the disorder strength, so that the system at stronger disorder then undergoes a liquid-to-glass transition (instead of the liquid-to-crystal transition found for weak disorder) as the density is increased from a low value. The continuous nature of the glass transition in the large disorder regime is in contrast with the first order transition from the liquid to a crystalline or glassy state (depending on the commensurability) at small values of the disorder strength. Thus, this work provides support to the prediction that the first-order liquid-to-crystal (Bragg glass) transition should change over to a continuous liquid-to-glass transition as the strength of the pinning disorder is increased above a critical value.

The rest of the paper is organized as follows. In section II, we define the model studied here and outline the numerical procedure used in this study. The numerical results obtained for the different transition lines in the density–disorder plane are described in detail in section III. Section IV contains a summary of our main results and a few concluding remarks.

## II Methods

### A The Free Energy functional

As discussed in the Introduction, our starting point is the free energy as a functional of the time-averaged local density $`\rho (𝐫)`$ at each point $`𝐫`$. We write this free energy functional in the form:

$$F[\rho ]=F_{RY}[\rho ]+F_s[\rho ]$$
(1)
where the first term in the right-hand side is the RY free enrgy functional  for hard spheres in the absence of disorder, and the second is the contribution arising from the presence of a quenched random potential. Thus we have:

$`\beta F_{RY}[\rho ]`$ $`=`$ $`{\displaystyle 𝑑𝐫\{\rho (𝐫)\mathrm{ln}(\rho (𝐫)/\rho _0)\delta \rho (𝐫)\}}`$ (2)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}C(|𝐫𝐫^{}|)\delta \rho (𝐫)\delta \rho (𝐫^{})}.`$ (3)
Here, we have defined $`\delta \rho (𝐫)\rho (𝐫)\rho _0`$ as the deviation of $`\rho (𝐫)`$ from $`\rho _0`$, the density of the uniform liquid, and taken the zero of the free energy at its uniform liquid value. In Eq.(3), $`\beta =1/(k_BT)`$, $`T`$ is the temperature and the function $`C(r)`$ is the direct pair correlation function  of the uniform liquid at density $`\rho _0`$, which can be analytically expressed in terms of the usual dimensionless density for hard spheres of diameter $`\sigma `$, $`n^{}\rho _0\sigma ^3`$, by making use of the Percus-Yevick approximation . This approximation is sufficiently accurate in the density ranges ($`n^{}1.0`$) considered in this paper. We write also:

$$\beta F_s[\rho ]=𝑑𝐫\delta \rho (𝐫)V_s(𝐫)$$
(4)
where $`V_s(𝐫)`$ is an external potential (in dimensionless form) representing the random, quenched disorder. We will assume that $`V_s`$ has zero mean and short-range Gaussian correlations as detailed below.

In order to carry out numerical work, we discretize our system. We introduce for this purpose a simple cubic computational mesh of size $`L^3`$ with periodic boundary conditions. On the sites of this mesh, we define density variables $`\rho _i\rho (𝐫_i)h^3`$, where $`\rho (𝐫_i)`$ is the density at site $`i`$ and $`h`$ the spacing of the computational mesh. It is known from previous work  that in the absence of any random potential, this discretized system crystallizes at sufficiently high densities if the quantities $`h`$ and $`L`$ are commensurate with a fcc structure with appropriate lattice spacing, whereas no crystalline state exists when the computational mesh is incommensurate with a fcc structure. Both commensurate and incommensurate systems exhibit  many glassy (inhomogeneous but aperiodic) minima of the free energy at densities higher than the value at which crystallization occurs in commensurate samples.

To model the random potential $`V_s(𝐫)`$, we introduce random variables $`\{V_i\}`$ defined at the sites of the computational mesh. These variables are uncorrelated with one another, and distributed according to a Gaussian probability distribution with zero mean and variance $`s`$. Thus, $`s`$ represents the dimensionless strength of the disorder. In terms of these quantities, the dimensionless free energy of our discretized system has the form

$`\beta F`$ $`=`$ $`{\displaystyle \underset{i}{}}\{\rho _i\mathrm{ln}(\rho _i/\rho _{\mathrm{}})(\rho _i\rho _{\mathrm{}})\}`$ (5)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}C_{ij}(\rho _i\rho _{\mathrm{}})(\rho _j\rho _{\mathrm{}})+{\displaystyle \underset{i}{}}V_i(\rho _i\rho _{\mathrm{}}),`$ (6)
where the sums are over all the sites of the computational mesh, $`\rho _{\mathrm{}}\rho _0h^3`$, and $`C_{ij}`$ is the discretized form of the direct pair correlation function $`C(r)`$ of the uniform liquid.

Our objective is to study the phase diagram of this system in the $`(n^{},s)`$ plane, in which in principle crystalline, liquid, and glassy phases may be found. The thermodynamics of hard spheres in the clean limit is determined by the dimensionless density $`n^{}`$ only. Our rescaling of the potential $`V_s`$ by $`\beta `$ (see Eq.(4)) ensures that $`s`$ is now the only additional relevant variable. We wish to locate various transition lines in this $`(n^{},s)`$ plane, that is, we wish to determine which phase (crystal, glass or liquid) is the thermodynamically stable one at different points in this plane. We also wish to know when and how a certain phase becomes metastable or unstable as we move around in the $`(n^{},s)`$ plane. In our mean-field description, different phases are represented by different minima of the free energy. If several local minima of the free energy are simultaneously present, then the minimum with the lowest free energy represents the thermodynamically stable phase and the other local minima correspond to metastable phases. A crossing of the free energies of two different minima signals a first-order phase transition. The point where a minimum becomes locally unstable (i.e. changes from a true minimum to a saddle point or disappears altogether) corresponds to a mean-field spinodal point representing the limit of metastability of the corresponding phase. A merging of the transition point with the spinodal points of the two phases signals a continuous phase transition in this description. Thus, a study of how the minima of the free energy of Eq.(6) evolve as $`n^{}`$ and $`s`$ are changed is sufficient for mapping out the mean-field phase diagram in the $`(n^{},s)`$ plane.

We locate the minima of the free energy by using a numerical procedure generalized from that originally developed for the clean case . This procedure works by changing the local density variables $`\{\rho _i\}`$ in a way that ensures that these changes always decrease the free energy. Given an initial configuration of the variables $`\{\rho _i\}`$, this procedure finds, by constantly moving downhill on the free-energy surface in the multidimensional configuration space spanned by the $`L^3`$ variables $`\{\rho _i\}`$, the local minimum whose basin of attraction contains the intial state. Thus, different local minima of the free energy can be located by using this minimization procedure for different, appropriately chosen, initial configurations.

As noted earlier, there are in our system three different kinds of free-energy minima: liquid, crystalline and glassy. In the clean limit ($`s=0`$), it is easy to distinguish among them: the liquid minimum has uniform density ($`\rho _i=\rho _{\mathrm{}}`$ for all $`i`$), the crystalline minimum has a periodic distribution of the density variables ($`\rho _i`$ is close to unity at mesh points corresponding to the sites of a fcc lattice, and close to zero at all other mesh points), and a glassy minimum exhibits a strongly inhomogeneous nonperiodic density distribution (some of the $`\rho _i`$’s are close to unity and the others are close to zero). This symmetry-based distinction among minima of different kind becomes less clear when the external random potential is turned on: for $`s0`$, the density distribution in the liquid phase is not completely homogeneous, and the crystalline state is not strictly periodic.

We use here, therefore, a procedure of “adiabatic continuation” to distinguish among the liquid, crystalline and glassy minima in the presence of the disorder. This procedure works as follows: We start with a minimum of a particular kind obtained at $`s=0`$ for a given value of $`n^{}`$. There is no difficulty in generating the liquid (and if appropriate the crystalline) configuration for the pure system. Glassy states at $`s=0`$ are easily obtainable also, in the right density ranges, by the procedures described in Ref. . Indeed, we have used in many cases the same density configurations obtained there, which were available as computer files. After thus choosing the initial state, we generate a set of uncorrelated random numbers $`r_i,i=1,\mathrm{},L^3`$, distributed according to a Gaussian with unit variance. A “realization” of the random potential $`\{V_i\}`$ is obtained by multiplying these random numbers by the strength parameter $`s`$. The initial $`s=0`$ minimum is then “followed” to finite $`s`$ by increasing $`s`$ in small steps $`\delta s`$ . After each step increase, the minimization routine is run, to find the nearest local minimum. Thus, in the first step of this process, the initial configuration is that of the minimum obtained at $`s=0`$ and the values of $`V_i`$ are set at $`(\delta s)r_i`$. The resulting density configuration obtained from the minimization routine is then used as the initial configuration for the next step, with the values of $`V_i`$ incremented to $`2(\delta s)r_i`$. During this process, the random variables $`\{r_i\}`$ are held fixed – only the strength parameter $`s`$ in increased in steps of $`\delta s`$. By iterating this procedure, minima of different kinds obtained at $`s=0`$ for a certain $`n^{}`$ are “followed” at constant density to the desired value of $`s`$. We use the terms “liquid”, “crystalline” and “glassy” to denote the continued $`s0`$ minima obtained from a $`s=0`$ minimum of the corresponding kind by using this continuation procedure without crossing transition lines. We will see below that even at large $`s`$, there are distinguishable differences in the structure of the different kinds of minima.

Once a minimum of the desired kind is obtained at a particular point in the $`(n^{},s)`$ plane, the free energy at the minimum reached, as well as the entire density configuration of the system at the minimum are obtained and can be analyzed. The translational correlations can be quantified by the two-point correlation function $`g(r)`$ of the density variables $`\{\rho _i\}`$. This function is defined as

$$g(r)=\underset{i>j}{}\rho _i\rho _jf_{ij}(r)/[\overline{\rho }^2\underset{i>j}{}f_{ij}(r)],$$
(7)
where the distance $`r`$ is measured in units of $`\sigma `$, $`\overline{\rho }_i\rho _i/L^3`$ is the average value of the $`\rho _i`$ variables at the minimum under consideration, and $`f_{ij}(r)=1`$ if the separation between mesh points $`i`$ and $`j`$ lies between $`r`$ and $`r+\mathrm{\Delta }r`$ ($`\mathrm{\Delta }r`$ is a suitably chosen bin size), and $`f_{ij}(r)=0`$ otherwise. This function represents the spatial correlation of the time-averaged local density, and is distinct from the equal-time, two-point density correlation function which is often called $`g(r)`$ in the literature. We also calculate $`\rho _{max}`$, the maximum value of the $`\rho _i`$ variables at the minimum, which gives additional information about the inhomogeneity when contrasted with $`\overline{\rho }`$ or its rescaled equivalent $`\rho _{av}\overline{\rho }(\sigma /h)^3`$ at the minimum.

In addition to examining the transitions by looking at discontinuities in $`F`$, $`g(r)`$ and the density configurations, we also directly check on the stability of the corresponding minima. The stability of a local minimum requires that all the eigenvalues of the Hessian matrix $`𝐌`$ whose elements are given by

$$M_{ij}\frac{^2(\beta F)}{\rho _i\rho _j}=\frac{1}{\rho _i}\delta _{ij}C_{ij}$$
(8)
evaluated at the minimum must be positive. This matrix is difficult to handle numerically if the minimum under consideration is strongly inhomogeneous, with some of the $`\rho _i`$’s very close to zero. In such cases, the $`1/\rho _i`$ in the first term on the right-hand side of Eq.(8) causes numerical difficulties. To avoid this problem, we consider instead the closely related matrix $`𝐌^{}`$ whose elements are given by

$$M_{ij}^{}\sqrt{\rho _i}M_{ij}\sqrt{\rho _j}=\delta _{ij}C_{ij}\sqrt{\rho _i\rho _j},$$
(9)
evaluated at the minimum under consideration. It is easy to show that an instability of the minimum corresponds to the vanishing of the smallest eigenvalue $`\lambda `$ of this matrix. In our numerical work, we calculate the value of $`\lambda `$ in order to check whether the minimum under study becomes unstable as $`n^{}`$ or $`s`$ is varied.

In our computations we have included the density range from $`n^{}=0.65`$ to $`n^{}=0.95`$, and values of $`s`$ from zero to about two. These are sufficient to encompass the phenomena that we wish to study. We have used three lattice sizes, $`L=12,15`$, and 25. For the last two we have used an incommensurate ratio $`h/\sigma =1/4.6`$, whereas for the smallest lattice we have taken the commensurate value $`h/\sigma =0.25`$.

## III Results

### A General considerations: Phase diagram

At a general point in the $`(n^{},s)`$ plane, the system may have a number of local minima, one of which is the absolute minimum while the others are metastable. Possible minima are that corresponding to the liquid (the one with uniform density at $`s=0`$ and its continuation to finite $`s`$), the crystalline minimum, by which we similarly mean the one with a periodic structure at $`s=0`$, and its continuation as described in the preceding section, and a large number of glassy minima. As the values of $`n^{}`$ or $`s`$ change, free energy minima may in general appear and disappear, and the free energy values of those that remain change. There will therefore be a number of instabilities and transitions, which are the main subject of our study.

Consider first the previously studied  case of the disorder-free system ($`s=0`$ line). There, only the uniform liquid minimum is present at low densities. As $`n^{}`$ increases, a crystalline minimum appears if the computational mesh is commensurate. When $`n^{}`$ is further increased, a density is reached at which the crystal becomes thermodynamically stable, that is, its free energy becomes lower than that of the liquid state. We will denote this density as $`n_D^{}`$. Regardless of commensurability, many glassy minima appear as the density is further increased. We denote by $`n_C^{}`$ the density at which the first glassy minimum makes its appearance. Alternatively, one may consider the evolution of the glassy minima as $`n^{}`$ is decreased from a large initial value, and define $`n_C^{}`$ as the density at which the last remaining glassy minimum becomes locally unstable and disappears: the free energy of this last remaining glassy minimum crosses that of the liquid at a density $`n_B^{}`$ which is somewhat higher than $`n_C^{}`$. This density corresponds to a liquid to glass transition. In the commensurate case, the density $`n_C^{}`$ is above the crystallization density $`n_D^{}`$, and the free energy of the crystalline minimum is lower than that of the glassy minima. Thus, the glass transition in the pure system occurs in a “supercompressed” regime where the crystalline state is the thermodynamically stable one.

When we include the effects of the disorder ($`s>0`$), we find yet another density, $`n_A^{}`$, at which the liquid minimum becomes locally unstable (i.e. ceases to exist as a local minimum of the free energy). For weak disorder, the value of $`n_A^{}`$ is large (substantially higher than $`n_B^{}`$) so that the four densities $`n_A^{}`$, $`n_B^{}`$, $`n_C^{}`$, and $`n_D^{}`$ are in decreasing order. Thus, we have four (three in the incommensurate case where the crystalline state is absent) functions $`n_X^{}(s)`$ with $`X=A,B,C,D`$, representing precisely the four transitions or instabilities defined above. We denote the corresponding lines in the $`(n^{},s)`$ plane as the $`A,B,C,D`$ lines. The determination of the location of these lines is one the main results of our work. These results will be discussed below, but to fix ideas and to make this discussion easier to follow, we show in Fig. 1 these four lines for the $`L=12`$ commensurate case. There,     FIG. 1.: The overall phase diagram of the hard sphere system in the density ($`n^{}`$) – disorder ($`s`$) plane, obtained for the $`L=12`$ commensurate sample. The meaning of the line labels is explained in the text. The results shown are averages over 5 realizations of the disorder.  the general structure of the phase diagram, including the general shape of the four lines $`n_X^{}(s)`$ can be seen. Similarly, we show in Fig. 2 the three lines $`n_X^{}(s),X=A,B,C`$ (from top to bottom) found in the incommesurate, $`L=25`$ system. The similarities and differences between the comensurate and incommensurate cases are discussed below. The lines in the phase diagram for the incommensurate $`L=15`$ case are within error bars the same as those shown in Fig. 2, so that the differences between Figs. 1 and 2 must be attributed to different commensurability rather than to different sample size.

There are certain trends that can be easily discerned   FIG. 2.: The overall phase diagram for the incommensurate case at size $`L=25`$ as explained in the text. The diamonds represent the $`A`$ line, the crosses the $`B`$ line and the dashed line is the $`C`$ line. Sample error bars have been indicated. They reflect sample-to-sample variations for six to twelve (the number increases with $`s`$) realizations of the disorder.  when one follows a free energy minimum as $`s`$ is increased at constant $`n^{}`$. If one starts from the uniform liquid minimum at $`s=0`$ and a relatively small value of $`n^{}`$, the free energy value at the minimum (initially zero according to our convention) decreases steadily with increasing $`s`$. For the case in Fig. 2 at $`n^{}=0.66`$, for example, $`\beta F`$ is close to $`180`$ at $`s=1.8`$. The density distribution becomes progressively less uniform, with $`\rho _{max}`$, which at $`s=0`$ equals the average value $`\overline{\rho }=\rho _{\mathrm{}}`$, rising by more than one order of magnitude as $`s`$ increases from zero to one. For a deep glassy state at a relatively large value of $`n^{}`$, the free energy is strongly negative even at $`s=0`$, and its value decreases further as $`s`$ increases. For example, at $`n^{}=0.78`$ the glassy minimum for $`L=25`$ with $`\beta F=63`$ at $`s=0`$ can be continued to a minimum with $`\beta F`$ equal to $`231`$ at $`s=1.8`$. The density distribution at a glassy minimum is considerably more inhomogeneous than that of the liquid minimum continued to the same value of $`s`$ and it is less sensitive to the value of $`s`$: the quenched disorder has less effect on a state that is inhomogeneous and disordered to begin with.

These trends in the behavior of liquid and glassy minima as $`s`$ is increased from zero are clearly illustrated by examining the pair correlation function $`g(r)`$, defined in   FIG. 3.: Liquid phase correlations. The pair correlation function $`g(r)`$ as defined in Eq. (7) plotted as a function of $`r`$, the distance in units of $`\sigma `$, for the liquid-like minimum at density $`n^{}=0.66`$. The curves shown, in order of increasing peak height at $`r=1`$, correspond to $`s=0.2,0.6,1.0,1.4,1.8`$. The system size is $`L=25`$.    FIG. 4.: The pair correlation function $`g(r)`$ for a glassy minimum. The curves shown correspond to the same values of $`L`$ and $`s`$ as in Fig.3, but for a glassy minimum at $`n^{}=0.78`$ as discussed in the text.  Eq.(7), at each minimum. In Fig. 3, we show $`g(r)`$ computed for the liquid minimum at size $`L=25`$ and density $`n^{}=0.66`$. The curves shown, in order of increasing value of the peak near $`r=1`$, correspond to increasing values of $`s=0.2,0.6,1.0,1.4,1.8`$. There is a clearly visible increase in structure, which becomes more evident as the value of $`s`$ increases beyond unity. However, this level of structure is still quantitatively different from that found for glassy minima at relatively high densities. This can be seen by comparing Fig. 3 with Fig. 4 where we plot $`g(r)`$ for a $`L=25`$ glassy minimum continued from $`s=0`$ to $`s=1.8`$ at $`n^{}=0.78`$, as mentioned in the preceding paragraph. We see that the $`s`$-dependence of the structure is now much weaker, and the heights of the peaks at $`r=0`$ and near $`r=1`$ are much larger than those in Fig. 3. These results can be compared to those found in the replica calculation . To make contact with those results, our $`g(r)`$ for the liquid-like minimum should be compared with the function $`g_0(r)`$ of the replica symmetric solution, and our $`g(r)`$ for a glassy minimum with the function $`g_1(r)`$ of the replica-symmetry-broken solution. Although, due to differences in the modeling of the random potential and effects of discretization in the present study (some of these effects are discussed in section IV), a detailed, quantitative comparison of our results with those of Ref. is not possible, it is clear that the main features we have discussed are qualitatively similar.

The crystalline minimum obtained for $`s=0`$ in commensurate systems at sufficiently high densities shows very little change in structure as it is followed to non-zero values of $`s`$. Any effects of weak pinning disorder on the crystalline order may be too subtle  to show up at the system sizes and discretization scales used here in the commensurate case.

### B Instability of the liquid minimum

We consider first the $`A`$ line, that is, the density at which the liquid minimum becomes locally unstable as $`n^{}`$ is increased from a low initial value, keeping $`s`$ fixed. This transition is detected at any desired value of $`s`$ in the following way. At a density previously determined to be well below the value of $`n_A^{}(s)`$ (this determination is easily performed by trial and error), one “follows” the $`s=0`$ liquid minimum, as previously explained, to the value of the disorder strength being studied. The density configuration at this minimum is the initial condition. Then, one proceeds to increase $`n^{}`$ by small intervals, thus moving up along a vertical line in the phase diagram. At every value of $`n^{}`$ that is reached, we run our minimization routine (using the configuration at the minimum obtained at the previous step as the starting point) to locate the nearest minimum. The density configuration at the minimum is analyzed and then used as the initial condition to study the next higher density.

In the initial stages of this process, the system remains in the liquid-like minimum, with little change in its properties. However, as $`n^{}`$ reaches the value $`n_A^{}(s)`$, discontinuities are found. These are more prominent for the larger system sizes (Fig. 2) and particularly dramatic for values of $`s`$ not too large. The free energy drops abruptly, as the liquid minimum becomes unstable and the system has to find some other nearby minimum (our numerical minimization procedure is designed to converge only to stable local minima of the free energy). Computationally, this is heralded by a very sharp and obvious increase in the number of iterations required by our numerical procedure to find the free energy minimum nearest to the starting configuration. This new minimum is invariably glassy, as one might expect, since a considerable number of glassy minima are close in configuration space to the liquid-like minimum . The value of the free energy at the minimum that the system has reached drops sharply as the $`n_A^{}(s)`$ value is crossed, because the free energies of glassy minima are considerably lower in the region of the $`(n^{},s)`$ plane being considered. Also, every measure of structure in the system inceases abruptly, since, as discussed above in connection with Figs. 3 and 4, glassy states are much more inhomogeneous than the liquid-like ones in this region of the $`(n^{},s)`$ plane.

An example of the behavior found is displayed in Figs. 5 and 6. In the main part of Fig. 5, we show the evolution of the free energy as $`n^{}`$ is increased     FIG. 5.: Discontinuities at the $`A`$ line. In the main plot, the free energy in dimensionless form for a $`L=12`$ sample ($`s=0.6`$) is plotted as a function of density. A sharp drop in the free energy is seen as the liquid minimum becomes unstable and the system switches to a glassy minimum. As shown in the inset, this switch is also reflected in the discontinuity in $`\lambda `$, the smallest eigenvalue of the matrix $`𝐌^{}`$ defined in Eq.(9).  in steps of 0.001, keeping $`s`$ fixed at 0.6 for a $`L=12`$ sample. One can clearly see that $`\beta F`$ varies little while the system remains in the liquid minimum and jumps abruptly as this minimum becomes unstable near $`n_A^{}0.78`$. The behavior for the larger incommensurate samples is quite similar, the main difference being that the drop in $`\beta F`$ is much larger, and that the transition occurs, at this value of $`s`$, at $`n_A^{}0.84`$ for both $`L=15`$ and $`L=25`$. The value of $`n_A^{}`$ can readily be found to very high precision and it varies little as one averages over different realizations of the quenched disorder, for the same $`s`$. The error bars shown in Fig. 2 correspond to an average over six to twelve realizations (the larger number at larger $`s`$.) The results in Fig. 1 are averages over five realizations. In the inset, we show that the smallest eigenvalue, $`\lambda `$, of the matrix $`𝐌^{}`$ defined in Eq.(9) approaches zero as $`n^{}`$ approaches $`n_A^{}`$ from below. This is as would be expected – as noted in section II, the instability of a local minimum is signalled by the vanishing of $`\lambda `$.     FIG. 6.: Example of how the system becomes more structured as the $`A`$ line is crossed for an $`L=12`$ sample at $`s=0.6`$. The height $`g_{max}`$ of the first finite-$`r`$ peak in $`g(r)`$ increases discontinuously, the density nonuniformity represented by $`\rho _{max}`$ exhibits a large increase, and the average density $`\rho _{av}`$ shows a small discontinuous increase. In order to be able to use a single vertical scale, we have displayed $`\rho _{av}`$ rather than $`\overline{\rho }`$ (see text).  In Fig. 6, we show three quantities which characterize the nature of the density distribution at a minimum. These are: $`g_{max}`$, the value of the pair correlation function $`g(r)`$ at its first finite $`r`$ maximum (near $`r=1`$); $`\rho _{max}`$, the maximum value of the $`\rho _i`$; and the dimensionless average density $`\rho _{av}`$ defined in section II. All these quantities exhibit discontinuous changes as the system switches minima at $`n^{}=n_A^{}0.78`$, (or at $`n_A^{}0.84`$ in the incommensurate case). $`g_{max}`$ remains close to unity as long as the system stays in the liquid state, and then jumps to a substantially larger value consistent with the increased short-range order present in a glassy minimum. This can also be seen from Figs. 3 and 4. The value of $`\rho _{max}`$ also increases by a considerable amount, indicating the increased inhomogeneity of a glassy minimum relative to the liquid-like one. The small increase in the value of $`\rho _{av}`$ reflects that the average density at a glassy minimum is slightly higher than that at the liquid-like minimum.

The behavior discussed above changes as $`s`$ is increased. The change occurs near $`s=1`$ for $`L=12`$, and at somewhat larger $`s`$ for the other system sizes, as the $`A,B,C`$ lines become very close. The results obtained in the larger-$`s`$ regime are described in the next subsection.

### C Instability of glassy minima and the liquid to glass transition

To find the $`B`$ and $`C`$ lines, we must start with a carefully chosen glassy configuration at a relatively high $`n^{}`$ and fixed $`s`$, and then follow this configuration to lower densities by decreasing $`n^{}`$ in small steps ($`\delta n^{}0.001`$), keeping the value of $`s`$ unchanged. This is continued until the minimum becomes unstable and the minimization routine converges to a new minimum which, if the starting minimum is chosen as described below, turns out to be the liquid-like one. The density at which this occurs defines the value of $`n_C^{}`$. Then, comparing the free energy of the glassy minimum with that of the liquid minimum obtained for the same realization of the disorder, it is easy to determine the value of $`n_B^{}`$ – this is the value of $`n^{}`$ at which the two free energies are equal.

The determination of the appropriate starting glassy minimum is nontrivial. Glassy minima for $`s0`$ are obtained by continuation from those of the pure system ($`s=0`$). One may think that the best choice would be to take the glassy minimum with the lowest free energy at the starting ($`n^{},s`$) point. In practice, this is difficult to implement because an exhaustive enumeration of all the glassy minima is computationally very hard. The glassy minimum with the lowest free energy at a particular point in the $`(n^{},s)`$ plane does not in general continue to have the lowest free energy as the values of $`n^{}`$ and $`s`$ are changed. Also, in the pure system, all the configurations obtained by applying one of the symmetry operations of the computational mesh to the density configuration at a particular glassy minimum also correspond to local minima with exactly the same free energy. For $`s0`$, all these symmetry-related minima have to be considered separately because the presence of the random potential destroys the symmetries present in the pure limit.

We have not found a rigorous solution to this problem. Instead, we first carried out an exploratory study of how the locations of the $`B`$ and $`C`$ lines in the phase diagram depend on the choice of the initial glassy minimum. The following choices, among others, were considered in our initial exploration. (a) One of the low-lying $`s=0`$ glassy minima, continued to finite $`s`$. (b) Beginning with the same starting configuration as in (a) and a specific realization of the random variables $`\{V_i\}`$, minimize the random potential energy (the last term in Eq.(6)) with respect to all symmetry operations of the computational mesh. This attempts to find the configuration that minimizes, among all the symmetry related ones, the contribution of the random potential to the free energy but it is not quite rigorous because the minimization is performed using the values of $`\{\rho _i\}`$ at the $`s=0`$ minimum. (c) The glassy minima to which the system moves when the density is increased above the $`A`$ line, as discussed in the preceding subsection.

The outcome of this study is that the locations of the $`B`$ and $`C`$ lines in the $`(n^{},s)`$ plane are not sensitive to the choice of the glassy minimum as long as it is one of the low-lying minima. (Even when we have deliberately or accidentally chosen a “wrong”, non-low-lying, minimum, we have found that the system often spontaneously makes a glass-to-glass transition  to a low-lying minimum as one decreases $`n^{}`$ above the $`B`$ line.) The variation of the values of $`n_B^{}`$ and $`n_C^{}`$ for different choices of the glassy minimum is comparable to the uncertainty of these values arising from sample-to-sample variations caused by differences in the realization of the disorder. The results described below were obtained (unless otherwise indicated) from runs in which a low-lying glassy minimum obtained from continuation of one at $`s=0`$ was taken to be the initial state for the density-lowering run. As explained at the beginning of this subsection, the initial configuration is followed to lower densities at fixed $`s`$ and the $`n_C^{}(s)`$ and $`n_B^{}(s)`$ points are found for that value of $`s`$. For relatively small values of $`s`$, the signatures of the $`C`$ instability are very easy to detect: they are similar to the discontinuities shown in Figs. 5 and 6. At larger values of $`s`$ more care is required.

For small values of $`s`$, the $`A`$, $`B`$ and $`C`$ lines are well separated from one another. However, as the value of $`s`$ is increased, these three lines begin to approach each other. As shown in Figs. 1 and 2, the separation between lines $`A`$ and $`B`$ decreases rather rapidly with increasing $`s`$, while the separation between lines $`B`$ and $`C`$ decreases more slowly. Finally, near $`s=1`$, these three lines appear to merge with one another for the $`L=12`$ system. For $`L=25`$ (and also for $`L=15`$), the separation between them does not exceed the combined error bars, but separate $`B`$ and $`C`$ transitions can be detected in most (not all) “runs” (i.e. realizations of the disorder) as explained in detail below. At larger $`s`$, it becomes increasingly difficult to resolve these three lines as they come close to each other. Since lines $`A`$ and $`C`$ represent, respectively, the limits of stability of the liquid and glassy minima and line $`B`$ corresponds to the first-order liquid-glass transition, a merging of these three lines suggests that this transition becomes continuous as $`s`$ is increased beyond a “tricritical” value which would be close to unity for the $`L=12`$ commensurate sample and somewhat larger for the incommensurate samples. Another possibility is that the first-order liquid-glass transition disappears beyond a critical point near $`s=1`$.

To examine the behavior in this region more closely, we have carried out several numerical experiments in which the value of $`n^{}`$ is “cycled” through the liquid-glass transition, keeping $`s`$ fixed at values close to unity. In this way, the three lines are detected in the same “run”. These numerical experiments are similar to simulations of hysteresis in magnetic phase transitions. We start with the liquid minimum at a low value of $`n^{}`$ (below line $`C`$), and increase $`n^{}`$ in small steps, keeping $`s`$ fixed. The liquid minimum is thus followed to higher densities until it undergoes a rapid change signalling a possible instability. The process of increasing $`n^{}`$ in small increments is continued for a few more steps, and then the local minimum so obtained is followed to lower densities by decreasing $`n^{}`$ is small steps. This is continued until the starting value of $`n^{}`$ is reached. If the liquid-glass transition at the chosen value of $`s`$ is first-order with the three densities $`n_A^{}`$, $`n_B^{}`$, and $`n_C^{}`$ separated from one another, then the cycling experiment described above should exhibit clear evidence of hysteresis. This is indeed what we find, for     FIG. 7.: Hysteresis and discontinuities across the liquid-glass transition at small values of $`s`$. In the main plot, the dimensionless free energy of the stable minimum is plotted vs. density as one cycles across the $`A`$, $`B`$ and $`C`$ lines as explained in the text. Hysteresis is clearly observed. In the inset, the quantity $`g_{max}`$ is shown. The results shown are at $`s=0.8`$ for a commensurate $`L=12`$ system, but the same behavior is found in this range of $`s`$ for incommensurate systems.  all system sizes and at every run, if the value of $`s`$ is lower than a certain critical value. A typical example is shown in Fig. 7 which shows the results for a $`L=12`$ sample at $`s=0.8`$. The hysteresis in the free energy and $`g_{max}`$ (shown in the inset) is evident: the liquid minimum becomes unstable at $`n_A^{}0.735`$ as $`n^{}`$ is increased from a low initial value, while the glassy minimum found for $`n^{}>n_A^{}`$ can be continued all the way down to $`n_C^{}0.720`$ before it becomes unstable. The liquid-glass transition occurs at $`n_B^{}0.725`$ where the two branches of the free energy cross. The same situation occurs for the incommensurate $`L=25`$ system except that the values of the transition points are $`n_A^{}0.79`$, $`n_B^{}0.73`$ and $`n_C^{}0.71`$ for $`s=0.8`$. The results at $`L=15`$ are, within error bars, the same as those for $`L=25`$ at this value of $`s`$.

The behavior in Fig. 7 is to be contrasted with that shown in Fig. 8 which displays the results of the cycling experiment on a $`L=12`$ sample at $`s=1`$. The distribution of the random variables $`\{r_i\}`$ in this sample is the same as that of Fig 7 – only the strength of the disorder is changed. In this figure, there is no evidence of     FIG. 8.: Cycling across the liquid-glass transitions for $`s`$ = 1 in a $`L=12`$ commensuarte system. The same quantities are plotted as in Fig. 7, and now no hysteresis is seen. Incommensurate systems exhibit the same behavior at somewhat larger values of $`s`$, but not in all runs.  hysteresis in the free energy. The plot of $`g_{max}`$ shown in the inset exhibits a sharp change near $`n^{}=0.706`$ for both inreasing-$`n^{}`$ and decreasing-$`n^{}`$ runs, and the results for the two runs are nearly identical. Given the rounding off errors associated with the numerical procedures we use, the small differences between the increasing-$`n^{}`$ and decreasing-$`n^{}`$ values of $`g_{max}`$ are likely to be insignificant. We, therefore, conclude that at least within the resolution of our numerical procedures, there is no hysteresis at $`s=1.0`$ for this $`L=12`$ sample. This implies that the first order transition found in this sample for $`s=0.8`$ either becomes a continuous one or disappears as the value of $`s`$ is increased to 1.0. The sharp change in the value of $`g_{max}`$ near $`n^{}=0.706`$ suggests that the transition persists as a continuous one. To investigate this further, we have calculated the derivatives of $`g_{max}`$, $`\rho _{max}`$, and $`\rho _{tot}_i\rho _i`$ with respect to $`n^{}`$ in the region where these quantities change rapidly. We have also examined the behavior of $`\lambda `$ as a function of $`n^{}`$ in this region.     FIG. 9.: Derivatives with respect to $`n^{}`$ of the quantities $`\rho _{tot}`$, $`\rho _{max}`$ and $`g_{max}`$, as defined in the text, plotted as functions of $`n^{}`$ across a putatively continuous liquid-glass transition in a $`L=12`$ sample with $`s=1.0`$. The three quantities have sharp peaks at $`n^{}=0.706`$. The eigenvalue $`\lambda `$, also defined in the text, shows a pronounced dip at the same point.  Results for these quantities are shown in Fig. 9 for the same sample as that of Fig. 8. All the derivatives exhibit sharp peaks at $`n^{}=0.706`$, and the value of $`\lambda `$ goes through a minimum that is very close to zero at the same point. These results strongly suggest the occurrence of a continuous phase transition at $`n^{}=0.706`$. However, due to the limited resolution of our numerical calculations and the smallness of sample size, we can not rule out the possibility that the observed behavior reflects a sharp crossover rather than a true phase transition. Similar results are found for larger values of $`s`$. The continuation of the “transition line” beyond the point where the lines $`A`$, $`B`$ and $`C`$ come together is determined by locating the value of $`n^{}`$ at which the eigenvalue $`\lambda `$ reaches a minimum. The value of $`s`$ at which the $`A`$, $`B`$ and $`C`$ lines merge and the hysteresis in the cycling experiment disappears is found to be weakly dependent on the realization of the disorder – it varies between 1.0 and 1.2 for the five different $`L=12`$ samples studied.

For the incommensurate samples, the situation is somewhat more ambiguous. For $`L=25`$, the same cycling procedure shows that the transition is clearly hysteretic for all runs with $`s1.1`$. For larger values of $`s`$, an increasingly larger percentage of the runs is non-hysteretic (i.e. the results for $`\beta F`$ look like those in Fig. 8), while the other runs display a behavior similar to that in Fig. 7 but with much smaller discontinuities. As $`s`$ is increased beyond $`s=1.8`$, it becomes, in most of the “runs”, impossible to distinguish the discontinuities, if any, from computer noise. Thus, it is possible in this case to plot separate $`A`$, $`B`$, and $`C`$ lines all the way up to $`s=1.8`$. This accounts for the obvious difference in this respect between Figs. 1 and 2. The results for $`L=15`$ are quite consistent with those for $`L=25`$, but the smaller system size makes all interpretations more difficult. Thus, it is more difficult to identify the precise position of any well-defined tricritial point (or a critical point) from the results for the incommensurate samples. One might alternatively say that these incommensurate results are indicative of a crossover. It is not possible to completely rule out that the behavior is different for the commensurate and incommensurate samples, or that the poorer resolution of the smaller samples masks discontinuous behavior in some of the larger $`s`$ runs.

### D Crystallization

To study how the crystallization density $`n_D^{}`$ changes as $`s`$ is increased from zero (i.e. the location of the line $`D`$), we start with the crystalline minimum obtained for a commensurate sample at a large value     FIG. 10.: Free energy crossing at the crystallization transition. The solid and dotted lines represent, respectively, the free energies of the crystalline and liquid-like minima of a $`L=12`$ sample with $`s=0.6`$. Their crossing point is the density $`n_D^{}(s)`$.  of $`n^{}`$. We then find the symmetry related configuration that minimizes the random potential energy for a particular realization of the disorder and continue this configuration to the desired value of $`s`$. This configuration is then continued to smaller values of $`n^{}`$ by decreasing $`n^{}`$ in small steps. The crystalline minimum turns out to be quite robust under changes of the density and the strength of the disorder – the minimization routine converges rather quickly to the new minimum as the value of $`n^{}`$ or $`s`$ is changed by a small amount. While decreasing the value of $`n^{}`$, we keep track of the free energy of the crystalline minimum and find the value of $`n^{}`$ at which this free energy crosses that of the liquid minimum for the same realization of the disorder. For relatively small values of $`s`$, the crossing point determines the value of $`n_D^{}`$ for the chosen value of $`s`$. Typical results for the crossing of these two free energies are shown in Fig. 10. Our results for line $`D`$, averaged over five realizations of the disorder, are shown in Fig. 1. The crystallization transition is strongly first order for all values of $`s`$. In the small $`s`$ regime, the crystalline minimum has the lowest free energy for all densities above line $`D`$. Therefore, the lines $`A`$, $`B`$ and $`C`$ do not have any equilibrium thermodynamic significance in this regime for a commensurate system: the liquid-glass transition at line $`B`$ can be observed only if the crystallization transition at line $`D`$ is avoided, e.g. by rapid compression.

As shown in Fig. 1, the crystallization line crosses the liquid-glass transition line at a point which is very close to that where the lines $`A`$, $`B`$ and $`C`$ seem to come together. Beyond this point, line $`D`$ is determined by the crossing of the free energies of glassy and crystalline minima. The procedure is quite analogous to that shown in Fig. 10. This line, therefore, represents a first order transition between crystalline and glassy states in this regime. The phase diagram of Fig. 1 implies that the system undergoes a first order liquid-to-crystal transition for small values of $`s`$ as the density is increased from a low initial value. However, as the value of $`s`$ is increased above a critical value (which is close to unity for the $`L=12`$ system), the transition as $`n^{}`$ is increased becomes a continuous liquid-to-glass transition (or perhaps a sharp crossover). The glassy state then undergoes a first order transition to the crystalline state as the density is increased further. The observed curvature of line $`D`$ for large $`s`$ also implies that the system would undergo a first order crystal-to-glass transition as the strength of the disorder is increased at constant density.

## IV Summary and Discussion

We have mapped out the mean-field phase diagram of a hard sphere system in the presence of a quenched random potential by numerically studying the evolution of the minima of a model free energy as a function of the density $`n^{}`$ and the strength $`s`$ of the disorder. The phase diagram in the $`(n^{},s)`$ plane exhibits liquid, glassy and crystalline (for commensurate samples) phases. We find that the standard first order crystallization transition which occurs at $`s=0`$ upon increasing $`n^{}`$ retains its character at small $`s`$ as a first order transition from a weakly inhomogeneous liquid phase to a nearly crystalline state. The density at which this transition occurs increases somewhat with the strength of the disorder. We also find for all samples a liquid to glass transition in the metastable, “supercompressed” regime. This transition is first-order for small $`s`$, but within the resolution of our results it appears to become continuous as $`s`$ is increased above a critical value. This critical value is larger for incommensurate samples. The crystallization line in the $`(n^{},s)`$ plane crosses the glass transition line near the point where the glass transition becomes continuous. Thus, the first order crystallization transition found for small $`s`$ as the density is increased from a small initial value is replaced, for sufficiently large $`s`$, by a continuous liquid to glass transition. The phase diagram also shows, at larger $`s`$, a first order crystal to glass transition as the strength of the disorder is increased at constant density.

All the qualitative features of our phase diagram (i.e. its topology, the shapes of the transition and instability lines, and the nature of the transitions) are identical to those found in the replica-based analytic study  of the same syatem . Two of our most important results, namely the change in the nature of the liquid to glass transition as the strength of the disorder is increased above a critical value, and the crossing of the crystallization and glass transition lines above this critical value of the disorder strength, were also found in the analytic calculation. This similarity between the results of two studies using extremely different methodologies strongly suggests that the qualitative features of our phase diagram are correct, at least at the mean-field level. The considerable quantitative differences that exist between the numerical and replica results, i.e. that all the transition and instability lines in our phase diagrams lie at substantially lower densities than those obtained in the analytic study, have the same origin as the discrepancy between our $`s=0`$ results and those of molecular dynamics simulations  of the pure hard sphere system. As noted in our earlier work , these differences result from the discretization of the free energy functional. The use of a simple cubic mesh of spacing $`h0.2\sigma `$ in the discretization procedure increases the relative stability of inhomogeneous local minima of the free energy and thus leads to substantially lower values for the densities at which crystallization and the glass transition occur. The quantitative differences between our results in Fig. 1 and those in Fig. 2, on the other hand, appear to arise chiefly from the incommensurability of the latter sample, rather than from the slight difference in the values of $`h`$, or even from that in the values of $`L`$: we have found negligible sample-size effects in comparing the $`L=15`$ results to those at $`L=25`$ at the same value of $`h`$. The effects of discretization would presumably disappear for $`h`$ much smaller than the width of the approximately Gaussian density distributions near the points where the particles are localized at an inhomogeneous minimum of the continuum free energy functional. Unfortunately, a numerical calculation with such small values ($`0.01\sigma `$) of $`h`$ would require dealing with a very large number (of the order of $`10^6`$) of variables $`\{\rho _i\}`$. This appears to be computationally intractable.

Our phase diagram is a mean-field one – possible effects of fluctuations are not included in our calculation. The first-order crystallization transition is not expected to be strongly affected by fluctuations. The situation is more complex for the glass transition because there are a large number of glassy local minima. When the effects of fluctuations are included, the system might visit a large number of different glassy minima during its evolution over a long period of time, and thus behave like a liquid in the sense that the particles would no longer be localized in space and the time-averaged local density would be only weakly inhomogeneous. A true thermodynamic glass transition would occur only if the characteristic time scale for transitions between different glassy minima diverges in the thermodynamic limit. Whether this happens in the pure system is still a highly controversial issue. Further investigations of this question for systems of particles in the presence of quenched disorder would be very worthwhile. Also, the presence of multiple low-lying glassy minima of the free energy is expected to lead to slow relaxation even if no thermodynamic glass transition is present. Therefore, signatures of the mean-field glass transition found in our study should show up in the dynamics of the system even if no thermodynamic glass transition occurs when fluctuations are taken into consideration.

Our density–disorder phase diagram exhibits qualitative similarities with the field–temperature phase diagram of some high-$`T_c`$ superconductors in the presence of random point pinning. For a system of vortices in the mixed phase of type-II superconductors, the temperature $`T`$ plays the role of the density $`n^{}`$ of the hard sphere system – increasing $`T`$ is analogous to decreasing $`n^{}`$. As pointed out in the Introduction, increasing the magnetic field $`H`$ is believed  to increase the effective strength of the pinning disorder. Using these analogies, one can translate, in a very crude and qualitative sense, our phase diagram in the $`(n^{},s)`$ plane to a phase diagram for superconductors in the $`(T,H)`$ plane. Then, our result that the crystallization transition at weak disorder is replaced by a continuous glass transition as the strength of the disorder is increased translates into the statement that for superconductors, the first-order liquid to Bragg glass transition at low fields should change over to a continuous glass transition as the field is increased. As noted in the Introduction, this is precisely the behavior found in experiments on a family of high-$`T_c`$ superconductors. The experimentally obtained phase diagram of these superconductors also exhibits a Bragg glass to amorphous solid transition as the field $`H`$ is increased at low temperatures. This is analogous to the crystal to glass transition found in our phase diagram as the strength of the disorder is increased at constant density. Further evidence in support of this analogy is provided by a recent numerical study  that suggests that the high-field, low-temperature phase of high-$`T_c`$ superconductors (the so-called vortex glass phase) is very similar to a structural glass. In view of these similarities, an extension of our calculation to a system of pancake vortices in layered superconductors with random point pinning, using the appropriate form of the free energy, would be of obvious interest.

We are not aware of any experimentally studied system that provides a direct and precise physical realization of the model studied here. Colloidal suspensions in the presence of a time-independent, spatially random external potential (produced, for example, by suitably configured laser fields ) would probably provide a close approximation to our model. Since simple liquids with short range pair potentials which are strongly repulsive at short distances behave in many ways like a hard sphere liquid, our calculation is expected to apply, at least qualitatively, to such systems also.

###### Acknowledgements.

We are grateful to F. Thalmann, G. I. Menon, A. K. Sood, S. Bhattacharya, G.F. Mazenko and T. Witten for helpful discussions or comments.


format: text

text:
# Compact Stellar Systems in the Fornax Cluster: Super-massive Star Clusters or Extremely Compact Dwarf Galaxies?

## 1 Introduction

In cold dark matter (CDM) galaxy formation, small dense halos of dark matter collapse at high redshift and eventually merge to form the large virialised galaxy clusters seen today. The CDM model is very good at reproducing large-scale structure, but only very recently have the best numerical simulations (Moore et al. 1998) had the resolution to trace the formation of small halos within galaxy clusters, with masses $`10^9`$M. We do not yet know what the lower mass limit is for the formation of halos in the cluster environment: determining the lower limit of galaxy mass in clusters will provide an important constraint on these models. Most of the smallest cluster galaxies are low surface brightness dwarfs for which mass estimates are very difficult, though comparison with field low surface brightness dwarfs would suggest that they may be dark matter dominated (e.g. Carignan & Freeman 1988).

In this paper we describe a population of small objects we have found in the Fornax Cluster (see also Minniti et al. 1998 and Hilker et al. 1999) which have high surface brightness. The origin and nature of these objects is not yet clear, but if they are a product of the galaxy formation process in clusters, their high surface brightness will make it possible to probe the low-mass limit discussed above. They may represent extreme examples of compact low luminosity (“M32-type”) dwarf ellipticals. Alternatively, these objects may be super-massive star clusters—there is a very large population of globular clusters associated with the central galaxy of the Fornax Cluster, NGC 1399 (Grillmair et al. 1994). These objects are generally similar to Galactic globular clusters with similar colours and luminosities (Forbes et al. 1998). There is evidence that they are not all bound to the NGC 1399 system. Kissler-Patig et al. (1999) show that the kinematics of 74 of the globular clusters indicate that they are associated with the cluster gravitational potential rather than that of NGC 1399. They infer that the most likely origin of these globular clusters is that they have been tidally stripped from neighbouring galaxies. This has also been suggested by West et al. (1995), although the effect would be diluted by the large number of halo stars that would presumably be stripped at the same time.

Bassino, Muzzio & Rabolli (1994) suggest that the NGC 1399 globular clusters are remnants of the nuclei of dwarf nucleated galaxies that have survived the disrupture of being captured by the central cluster galaxy. A related suggestion is a second model proposed by West et al. (1995) that intra-cluster globular clusters could have formed in situ in the cluster environment. Bassino et al. (1994) conclude their discussion by noting that remnant nuclei an order of magnitude larger (and more luminous) than standard globular clusters would also be formed in significant numbers, but that existing globular cluster searches would not have included them. In Section 2 of this paper we describe how the observations of our Fornax Spectroscopic Survey have sampled this part of the cluster population by measuring optical spectra of all objects brighter than $`B_J=19.7`$ in the centre of the Fornax Cluster. In Section 3 we describe the properties of a new population of compact objects found in the cluster that appear to be intermediate in size between globular clusters and the smallest compact dwarf galaxies. We discuss the nature of these objects in Section 4 and show that higher resolution observations will enable us to determine if they are more like globular clusters or dwarf galaxies.

## 2 Discovery Observations: The Fornax Spectroscopic Survey

Our Fornax Spectroscopic Survey, carried out with the 2dF multi-object spectrograph on the Anglo-Australian Telescope (see Drinkwater et al. 2000), is now 87% complete in its first field to a limit of $`B_J=19.7`$. The 2dF field is a circle of diameter 2 degrees (i.e. $`\pi `$ square degrees of sky). We have measured optical spectra of some 4000 objects (some going fainter than our nominal limit) in a 2dF field centred on the central galaxy of the Fornax Cluster (NGC 1399). This survey is unique in that the targets (selected from digitised UK Schmidt Telescope photographic sky survey plates) include all objects, both unresolved (“stars”) and resolved (“galaxies”) in this large area of sky. The resolved objects measured are mostly background galaxies as expected with a minor contribution from Fornax Cluster members. The unresolved objects are mostly Galactic stars and distant AGN, also as expected, but some are compact starburst (and post-starburst) galaxies beyond the Fornax Cluster (Drinkwater et al 1999a).

Finally, in addition to the dwarf galaxies already listed in the Fornax Cluster Catalog (FCC: Ferguson 1989) which we have confirmed as cluster members, we have found a sample of five very compact objects at the cluster redshift which are unresolved on photographic sky survey plates and not included in the FCC. These new members of the Fornax Cluster are listed in Table 1 along with their photometry measured from the UKST plates.

Our 2dF measurements of unresolved objects are 80% complete in the magnitude range of these objects ($`17.5<b_j<20.0`$). There is therefore about one more similar compact object still to be found in in our central 2dF field. The number density of these objects is therefore $`6\pm 3`$ per 2dF field ($`\pi `$ square degrees). Two of the objects (the two brightest) were also identified as cluster members by Hilker et al. (1999). Hilker et al. measured spectra of about 50 galaxies brighter than $`V=20`$ in a square region of width 0.25 degrees at the centre of the Fornax Cluster. In the “galaxies” they include objects which were very compact, but still resolved. By contrast our own survey covers a much larger area and also includes all unresolved objects.

## 3 Properties of the compact objects

### 3.1 Sizes

These object images are unresolved and classified “stellar” in our UKST plate data, although imaging with the CTIO Curtis Schmidt shows that the brightest two objects have marginal signs of extended structure. In Fig. 1 we present R-band (Tech Pan emulsion + OG 590 filter) photographic images of these compact objects from the UKST. These were taken in seeing of about 2.2 arcseconds FWHM and the third object (Thales111Thales of Miletus was the first known Greek philosopher and scientist and possibly the earliest astronomer. 3) is resolved with a 3.2 arcsecond FWHM. Applying a very simple deconvolution of the seeing this corresponds to a scale size (HWHM) of about 80 pc. This is much larger than any known globular cluster, so this object, at least, is not a globular cluster. The other objects are all unresolved, so must have scale sizes smaller than this.

### 3.2 Luminosity and Colours

These new objects have absolute magnitudes $`13<M_B<11`$, based on a distance modulus of 30.9 mag to the Fornax Cluster (Bureau et al. 1996). These values are at the lower limit of dwarf galaxy luminosities (Mateo 1998), but are much more luminous than any known Galactic globular clusters (Harris, 1996) and the most luminous of the NGC1399 globulars (Forbes et al. 1998) which have $`M_B10`$. The luminosities of the compact objects are compared to several other populations of dwarf galaxy and star cluster in Fig. 2. We note that the magnitude limit of the 2dF data corresponds to an absolute magnitude of $`M_B11`$ here. In order of decreasing luminosity the first comparison is with the dwarf ellipticals listed in the FCC as members of the Fornax Cluster. The possible M32-type galaxies in the FCC are not included as none of them have yet been shown to be cluster members (Drinkwater, Gregg & Holman 1997). The Figure shows that the Fornax dEs have considerable overlap in luminosity with the compact objects, but morphologically they are very different, being fully resolved low surface brightness galaxies. Recently, several new compact dwarf galaxies have been discovered in the Fornax Cluster (Drinkwater & Gregg 1998) but these are all brighter than $`M_B=14`$ and do not match any of the objects we discuss here. Binggeli & Cameron (1991) measured the luminosity function of the nuclei of nucleated dwarf elliptical galaxies in the Virgo Cluster. The Figure shows that this also overlaps the distribution of the new compact objects. In this case the morphology is the same, so the compact objects could originate from the dwarf nuclei. The Figure also shows the luminosity functions of both the NGC 1399 globular clusters (Bridges, Hanes & Harris, 1991) and Galactic globular clusters (Harris 1996). These are quite similar and have no overlap with the compact objects.

For completeness we note that the luminosities of the compact objects have considerable overlap with the luminosities of dwarf galaxies in the Local Group (Mateo 1998), but even the most compact of the Local Group dwarfs, Leo I ($`M_B=11.1`$) would be resolved ($`r_e3^{\prime \prime }`$) in our images at the distance of Fornax. The only population they match in both luminosity and morphology is the bright end of the nuclei of nucleated dwarf ellipticals.

### 3.3 Spectral Properties

The 2dF discovery spectra of these compact objects are shown in Figure 3. They have spectra similar to those of early type dwarf galaxies in the sample (two are shown for comparison in the Figure) with no detectable emission lines. As part of the spectral identification process in the Fornax Spectroscopic Survey, we cross-correlate all spectra with a sample of stellar templates from the Jacoby et al. (1984) library. The spectra of the new compact objects were best fit by K-type stellar templates, consistent with an old (metal-rich) stellar population. The dE galaxies observed with the same system by contrast are best fit by younger F and early G-type templates. This gives some indication in favour of the compact objects being related to globular clusters, although we note that two of them were analysed by Hilker et al. (1999) in more detail without any conclusive results. We do not have the spectrum of a dE nucleus available for direct comparison, but since our 2dF spectra are taken through a 2 arcsec diameter fibre aperture, the spectrum of FCC 211, a nucleated dE, is dominated by the nucleus. This spectrum was fitted by a younger F-type stellar template, again suggestive of a younger population than the compact objects. We cannot draw any strong conclusions from these low-resolution, low signal-to-noise spectra.

### 3.4 Radial Distribution

The main advantage of our survey over the previous studies of the NGC 1399 globular cluster system (e.g. Grillmair et al. 1994, Hilker et al. 1999) is that we have complete spectroscopic data over a much larger field, extending to a radius of 1 degree (projected distance of 270 kpc) from the cluster centre. This means that we can determine the spatial distribution of the new compact objects. In Figure 4 we plot the normalised, cumulative radial distribution of the new compact objects compared to that of foreground stars and cluster galaxies. This plot is the one used to calculate Kolmogorov-Smirnov statistics and allows us to compare the distributions of objects independent of their mean surface densities. It is clear from the Figure that the new compact objects are very concentrated towards the centre of the cluster, at radii between 5 and 30 arcminutes (20–130 kpc). Their distribution is more centrally concentrated than the King profile fitted to cluster members by Ferguson (1989) with a core radius of 0.7 degrees (190 kpc). The Kolmogorov-Smirnov (KS) test gives a probability of 0.01 that the compact objects have the same distributions as the FCC galaxies: they are clearly not formed (or acreted) the same way as average cluster galaxies. To test the hypothesis that the compact objects are formed from nucleated dwarfs, we also plot the distribution of all the FCC nucleated dwarfs, as these are more clustered than other dwarfs (Ferguson & Binggeli 1994). However in the central region of interest here the nucleated dwarf profile lies very close to the King profile of all the FCC galaxies, so this does not provide any evidence for a direct link with the new compact objects.

West et al. (1995) suggest that a smaller core radius should be used for intra-cluster globular clusters (GCs). This profile, also shown in the Figure, is more consistent with the distribution of the new objects: the KS probability of the compact objects being drawn from this distribution is 0.39.

We also note that the radial distribution of the compact objects is much more extended than the NGC 1399 globular cluster system as discussed by Grillmair et al. and extends to three times the projected radius of that sample. It unlikely that all the compact objects are associated with NGC 1399. This is emphasised by a finding chart for the central 55 arcminutes of the cluster in Fig. 5 which indicates the location of the compact objects. They are widely distributed over this field and Thales 3 in particular is much closer to NGC 1404, although we note that its velocity is not close to that of NGC 1404 (see below).

### 3.5 Velocity Distribution

We have some limited information from the radial velocities of the compact objects. The mean velocity of all 5 ($`1530\pm 110\mathrm{km}\mathrm{s}^1`$) is consistent with that of the whole cluster ($`1540\pm 50\mathrm{km}\mathrm{s}^1`$) (Jones & Jones (1980). However, given the small sample, it is also consistent with the velocity of NGC 1399 ($`1425\pm 4\mathrm{km}\mathrm{s}^1`$) as might be expected for a system of globular clusters. Interestingly, the analysis of the dynamics of 74 globular clusters associated with NGC 1399 by Kissler-Patig et al. (1999) notes that their radial velocity distribution has two peaks, at about 1300 and 1800$`\mathrm{km}\mathrm{s}^1`$. Our sample is far too small to make any conclusions about the dynamics of these objects at present.

## 4 Discussion

We cannot say much more about the nature of these objects on the basis of our existing data. In ground-based imaging, they are intermediate between large GCs and small compact dwarf galaxies, so it becomes almost a matter of semantics to describe them as one or the other. The most promising way to distinguish between these possibilities is to measure their mass-to-light (M/L) ratios. If they are large, but otherwise normal, GCs, they will be composed entirely of stars giving very low M/L. If they are the stripped nuclei of dwarf galaxies we might expect them to be associated with some kind of dark halo, but we would not detect the dark halos at the small radii of these nuclei, so we would also measure small M/L values. Alternatively, these objects may represent a new, extreme class of compact dwarf elliptical (“M32-type”) galaxy. These would presumably have formed by gravitational collapse within dark-matter halos, so would have high mass-to-light-ratios, like dwarf galaxies in the Local Group (Mateo 1998). One argument against this interpretation is the apparent lack of M32-like galaxies at brighter luminosities (Drinkwater & Gregg 1998). If the compact objects are dwarf galaxies, they will represent the faintest M32-like galaxies ever found. They may also fill in the gap between globular clusters and the fainter compact galaxies in the surface brightness vs. magnitude distribution given by Ferguson & Binggeli (1994).

A further possibility is that these are small scale length ($`100`$ pc) dwarf spheroidal galaxies of only moderately low surface brightness. While Local Group dSphs of equivalent luminosities generally have substantially larger scale sizes (and consequently lower surface brightnesses) (Mateo, 1998), Leo I for example has $`M_B=11.0`$, and a scale length of only 110 pc (Caldwell et al., 1992), but as we discuss above this would be resolved in our existing imaging.

Our existing data will only allow us to estimate a conservative upper limit to the mass of these objects. If we say that the core radii of the objects are less than 75 pc and the velocity dispersions are less than 400$`\mathrm{km}\mathrm{s}^1`$ (the resolution of our 2dF spectra) we find that the virial mass must be less than $`10^{10}\text{M}\text{}`$. For a typical luminosity of $`M_B=12`$ this implies that $`M/L<2\times 10^3`$. This is not a very interesting limit, so we plan to reobserve these objects at higher spectral resolution from the ground and higher spatial resolution with the Hubble Space Telescope (HST) in order to be sensitive to $`M/L100`$. This will allow us to distinguish globular clusters from dwarf galaxies.

In order to demonstrate what we could measure with high-resolution images, we present two extreme possibilities in Figure 6: a very compact Galactic globular cluster and a dwarf galaxy with an $`r^{1/4}`$ profile ($`r_e=0.2`$ arcsec), both normalised to magnitudes of $`B=19`$ ($`V=18.4`$) and the Fornax cluster distance. We also plot the PSF of the Space Telescope Imaging Spectrograph (STIS) in the Figure for reference. The globular cluster profile is that of NGC 2808 (Illingworth & Illingworth 1976) with the radius scaled to the distance of the Fornax Cluster and the surface brightness then scaled to give the desired apparent magnitude. The globular cluster profile is very compact and will only just be resolved with HST, but it will clearly be differentiated from the dwarf galaxy profile.

In addition to measuring the size of these objects for the mass measurement, the radial surface brightness profiles may also give direct evidence for their origin and relationship to other kinds of stellar systems. For example, if they are the stripped nuclei of galaxies, the remnants of the outer envelope might show up in the HST images as an inflection in the surface brightness profile at large radius.

## 5 Summary

We have reviewed the observed properties of these new compact objects discovered in the Fornax Cluster. Their luminosities are intermediate between those of known globular clusters and compact dwarf galaxies, but they are consistent with the bright end of the luminosity function of the the nuclei of nucleated dwarf ellipticals. The 2dF spectra are suggestive of old (metal-rich) stellar populations, more like globular clusters than dwarf galaxies. Finally the radial distribution of the compact objects is more centrally concentrated than cluster galaxies in general, but extends further than the known globular cluster system of NGC 1399.

These objects are most likely either massive star clusters (extreme globular clusters or tidally-stripped dwarf galaxy nuclei) or very compact, low-luminosity dwarf galaxies. In the latter case these new compact objects would be very low-luminosity counterparts to the peculiar compact galaxy M32. This would be particularly interesting given the lack of M32-like galaxies at brighter luminosities (Drinkwater & Gregg 1998). With higher resolution images and spectra we will be able to measure the mass-to-light ratios of these objects and determine which of these alternatives is correct.

## Acknowledgements

We thank the referee for helpful suggestions which have improved the presentation of this work. We wish to thank Dr. Harry Ferguson for helpful discussions and for providing the STIS profile. We also thank Dr. Trevor Hales for assistance in the naming of the objects. MJD acknowledges support from an Australian Research Council Large Grant.

## References

Caldwell, N., Armandroff, T.E., Seitzer, P., Da Costa, G.S., 1992, AJ, 103, 840

Bassino, L.P., Muzzio, J.C., Rabolli, M. 1994, ApJ, 431, 634

Binggeli, B., Cameron, L.M., 1991, A&A, 252, 27

Bridges, T.J., Hanes, D.A., Harris, W.E., 1991, AJ, 101, 469

Bureau, M., Mould, J.R., Staveley-Smith, L., 1996, ApJ, 463, 60

Carignan, C., Freeman, K.C. 1988, ApJ, 332, L33

Drinkwater, M.J., Gregg, M.D., 1998, MNRAS, 296, L15

Drinkwater, M.J., Gregg, M.D., Holman, B.A., 1997 in Arnaboldi M., Da Costa G.S., Saha P., eds, ASP Conf. Ser. Vol. 116, The Second Stromlo Symposium: The Nature of Elliptical Galaxies. Astron. Soc. Pac., San Francisco, p. 287

Drinkwater, M.J., Phillipps, S., Gregg, M.D., Parker, Q.A., Smith, R.M., Davies, J.I., Jones, J.B., Sadler, E.M., 1999a, ApJ, 511, L97

Drinkwater, M.J., Phillipps, S., Jones, J.B., Gregg, M.D., Deady, J.H., Davies, J.I., Parker, Q.A., Sadler, E.M., Smith, R.M. 2000, A&A, submitted

Ferguson H.C., 1989, AJ, 98, 367

Ferguson H.C., Binggeli, B., 1994, A&ARv, 6, 67

Forbes, D.A., Grillmair, C.J., Williger, G.M., Elson, R.A.W., Brodie, J.P. 1998, MNRAS, 293, 325

Grillmair, C.J., Freeman, K.C., Bicknell, G.V., Carter, D., Couch, W.J., Sommer-Larsen, J., Taylor, K. 1994, ApJ, 422, L9

Harris, W.E., 1996, AJ, 112, 1487

Hilker, M., Infante, L., Vieira, G., Kissler-Patig, M., Richtler, T., 1999, A&AS, 134, 75

Illingworth, G., Illingworth, W. 1976, ApJSup, 30, 227

Jacoby, G.H., Hunter, D.A., Christian, C.A., 1984, ApJSup, 56, 257

Jones, J.E., Jones, B.J.T. 1980, MNRAS, 191, 685

Kissler-Patig, M., Grillmair, C.J., Meylan, G., Brodie, J.P., Minniti, D., Goudfrooij, P., 1999, AJ, 117, 1206

Mateo, M., 1998, Ann. Rev. Astron. Astrophys., 36, 435

Minniti, D., Kissler-Patig, M., Goudfrooij, P., Meylan, G., 1998, AJ, 115, 121

Miller, L. A., Cormack, W., Paterson, M., Beard, S., Lawrence, L., 1992, in ‘Digitised Optical Sky Surveys’, eds. H.T. MacGillivray, E.B Thomson, Kluwer Academic Publishers, p. 133

Moore, B., Governato, F., Quinn, T., Stadel, J., Lake, G. 1998, ApJ, 499, L5

West, M.J., Cote, P., Jones, C., Forman, W., Marzke, R.O. 1995 ApJ 453 L77