text:
# Dimensional crossover and hidden incommensurability in Josephson junction arrays of periodically repeated Sierpinski gaskets

## I Introduction

A wide variety of disordered materials, including superconductors, is known to exhibit geometrical inhomogeneities over a broad range of length scales. The properties of such systems can be conveniently described in terms of percolation , the simplest idea to understand randomness. Percolation can be regarded as a geometrical phase transition taking place at a ”critical concentration” $`p_c`$ separating a phase of finite clusters ($`p<p_c`$) from a phase where an infinite cluster is present ($`p>p_c`$). Like other critical phenomena, it is characterized by a correlation length $`\xi _p(p)`$ which diverges at the percolation threshold $`p_c`$. Right at $`p_c`$, a system with percolative disorder exhibits a natural self-similar structure at all length scales and can therefore be modeled by a family of scale-invariant fractal lattices, such as the Sierpinski gasket (SG), which has been suggested  to mimic the essential geometrical features of the percolating cluster’s backbone. In the critical region above and below $`p_c`$, where $`\xi _p`$ is finite, the nature of the geometry depends on the length scale l at which one is probing the system: if $`\text{l}<\xi _p`$, its structure is fractal, whereas it can be regarded as homogeneous with conventional euclidian features for $`\text{l}>\xi _p`$.

Allowing an accurate control of both the nature and the amount of disorder and exhibiting properties quite sensitive to dimensionality, Josephson junction arrays and wire networks prepared with modern micro- and nanofabrication techniques provide ideal model systems to investigate the dimensional crossover from the euclidian to the fractal regime. The first step in this direction was made by Gordon et al. , who investigated the superconducting-to-normal phase boundary $`T_c(B)`$ of aluminum wire networks formed by periodically repeated SGs and of analogous networks with percolative geometry exposed to a perpendicular magnetic field $`B`$. In those experiments the magnetic length, $`\text{l}(B)(\varphi _0/B)^{1/2}`$, which is a measure of the typical nearest-neighbor distance between the vortices present in the system, was shown to be the relevant length scale to explore the euclidian-fractal (EF) crossover ($`\varphi _0`$ is the superconducting flux quantum). While the scaling behavior of the phase boundary of the SG networks was found to exhibit a crossover from the euclidian to the fractal regime consistent with theoretical predictions  based on extensions  of the Ginzburg-Landau theory and allowing for comparison  with the anomalous diffusion exponent, no EF crossover was observed in percolative networks which, surprisingly, were found to behave like a homogeneous system at all length scales $`\text{l}(B)`$ probed in the experiment.

More recently, compelling evidence for the EF crossover in a percolative system emerged from ac conductance measurements  performed in zero magnetic field on (unfrustrated) site-diluted Josephson junction arrays with site occupation probabilities $`p`$ very close to $`p_c`$. In these experiments the crossover was controlled by the driving angular frequency $`\omega `$, which determines the ratio of the impedances associated with the two types of links forming the random network. According to Efros and Shklovskii , the increase of this ratio with decreasing $`\omega `$ also leads, if $`p`$ is sufficiently close to $`p_c`$, to a crossover from the fractal to the euclidian regime, as confirmed by the experiments of Ref. .

In this article we report the results of experimental and theoretical studies of a model superconductor sharing essential geometrical features with a percolative system. The samples we have investigated are arrays of proximity-effect coupled SNS junctions (where S stands for superconductor and N for normal metal) consisting of SGs connected to each other at the vertices in such a way as to form a regular triangular lattice. As can be seen from Fig. 1(a), where a part of an array of second-order gaskets is shown, in these systems the linear size $`L_N=2^Na`$ of an individual gasket of order $`N`$ can be regarded as playing the role of $`\xi _p`$ ($`a`$ is the length of an elementary link of the gasket). As shown in detail in this work, a remarkable feature of the SG arrays is that in these systems the EF crossover is clearly manifest, in contrast to truly percolative systems, where its signatures are elusive.

The quantity at the heart of the present study is the magnetoinductance $`L(B)`$ of the SG arrays, extracted from measurements of their ac impedance. Its interest resides in the observation that, being inversely proportional to the areal superfluid density, it provides a tool to appreciate how the degree of superconducting phase coherence in the system changes with $`B`$ or, equivalently, with the level of frustration imposed by $`B`$. Previous impedance measurements performed on weakly frustrated arrays similar to those studied in this work focused merely on the fractal regime and demonstrated, in particular, the unusual scaling properties  of the vortex energy  as well as the asymptotic ($`B0`$) scaling behavior of the field-induced correction to the array inductance resulting from the hierarchical structure of the gaskets .

Although some preliminary evidence for a dimensional crossover was already reported in Ref. , the phenomenon was not exhaustively investigated. In the present work we rely on high-resolution studies of the complex fine structure of $`L(B)`$ to explore in detail both the euclidian and the fractal regimes of the SG arrays. Reflecting flux quantization phenomena occurring in loops with a hierarchical distribution of sizes up to the gasket size $`L_N`$, the fine structure provides a unique tool to reveal how the geometrical properties of the system change as the magnetic length $`\text{l}(B)`$ is swept through $`L_N`$. We show that, in terms of the frustration parameter $`f`$ expressing the magnetic flux threading an elementary triangular cell of a gasket in units of $`\varphi _0`$, the EF crossover occurs at $`f_{cN}=1/(2\times 4^N)`$. In the fractal regime ($`f>f_{cN}`$) the most relevant contributions to the fine structure of $`L(f)`$ are shown to arise from a particular set of ground states defined by $`f=Mf_{cN}`$, where $`M`$ is an integer. Corresponding to vortex configurations where the vortex lattice is strongly pinned by the hierarchical potential landscape created by the gaskets, these states are particularly robust against thermal fluctuations and are therefore quite prominent in the fine structure of $`L(f)`$.

A very interesting aspect emerged from the study of the magnetoinductance in the fractal regime. The analysis of the data revealed anomalous features (specifically, the suppression of the periodicity corresponding to a shift of $`f`$ by $`1`$ and of the symmetry with respect to $`f=1/2`$) inconsistent with theoretical predictions based on the description of the system in terms of a uniformly frustrated $`XY`$ model. We suggest that, because of the asymmetric shape of the superconducting islands forming the junction pattern of a gasket \[see Fig. 1(a)\], the screening currents flowing in the islands create a distortion of the magnetic field distribution in the array such that the fluxes threading the various loops slightly deviate from being proportional to their areas. This introduces an effective incommensurability (which we call ”hidden” to distinguish it from the ”geometric” one studied earlier  in wire networks with incommensurate cells) and perturbs the self-similarity of the gaskets. As a result, the system is no longer uniformly frustrated. We demonstrate that the anomalous features mentioned above can be quantitatively accounted for by a simple model, in which the area of the different plaquettes of a gasket is changed according to an appropriate deformation scheme.

Compelling evidence for the existence of the euclidian regime is provided by the array’s magnetoinductance for $`f<f_{cN}`$. Besides the absence of the power-law scaling behavior characteristic of the fractal regime, $`L(f)`$ contains structures which reflect the presence of ground states corresponding to vortex configurations commensurate with the underlying triangular lattice formed by the largest triangular loops of the array (below $`f_{cN}`$, it is energetically unfavorable for vortices to penetrate loops of smaller size), thereby allowing an unambiguous identification of the euclidian regime.

In general, these vortex configurations have the same structure as in the $`XY`$ model on a honeycomb lattice, but correspond to values of $`f`$ reduced by a factor of $`1/f_{cN}=2\times 4^N`$. However, the case $`f=(1/2)f_{cN}`$ requires special attention. In a honeycomb lattice the ground state of the corresponding (fully frustrated) $`XY`$ model is characterized by an accidental degeneracy, which (to lowest order) survives even in the presence of thermal fluctuations . Owing to the more complex structure of the system, however, this peculiar degeneracy is removed in our SG arrays, allowing to identify the vortex configuration in the ordered phase at $`f=(1/2)f_{cN}`$.

The paper is organized as follows. Experimental details are given in Sec. II. Sec. III is devoted to the fractal regime ($`f>f_{cN}`$). Relying on the methods developed in Ref., in Sec. III A we present the calculation of the magnetoinductance of a frustrated SG array for the particular set of frustrations $`f=Mf_{cN}`$ corresponding to the sequence of the most stable states, which are characterized by a relatively compact structure (with a unit cell not larger than a single gasket). Experimental data for the magnetoinductance in the same regime are presented and discussed in Sec. III B, where we show, in particular, how the anomalous features of $`L(f)`$, revealing the presence of hidden incommensurability, can be accounted for by a simple model, in which the areal changes of different plaquettes (related to the redistribution of the magnetic field) are determined by only one adjustable parameter. In Sec. IV we provide experimental evidence for the existence of the EF crossover and in Sec.V the fine structure of $`L(f)`$ in the euclidian regime below $`f_{cN}`$ is shown to be consistent with the existence of vortex states commensurate with the periodic lattice formed by the largest triangular cells of the arrays. A few concluding remarks are given in Sec. VI.

## II Experimental aspects

The samples studied in this work consist of second- ($`N=2`$) and fourth-order ($`N=4`$) gaskets sitting on the sites of, respectively, a $`313\times 313`$ and a $`78\times 78`$ triangular lattice and connected to each other at the vertices \[see Fig. 1(a)\]. Each gasket contains, respectively, $`3^{2+1}=27`$ and $`3^{4+1}=243`$ SNS Josephson junctions consisting of superconducting Pb islands proximity-effect coupled to each other by an underlying normal Cu layer. The geometrical and physical parameters of the junctions are almost identical to those of the array studied in Ref. . Most of the data presented below have been obtained in experiments performed on the array of second-order gaskets.

The sheet magnetoinductance $`L(f)`$ was inferred from measurements of the array’s linear sheet impedance $`Z=R+i\omega L`$ performed with a very sensitive SQUID-operated two-coil mutual inductance technique  at driving frequencies typically in the range 0.1-1.0 kHz. With this method we were able to resolve inductance changes of the order of 10 pH in swept-frustration measurements, in which $`f`$ could be tuned with a precision better than $`10^3`$. The experimental data are presented and analyzed in terms of $`L^1(f)`$, the quantity measuring the degree of superconducting phase coherence in the samples. When needed, the resistive component $`R(f)`$, related to dissipative vortex motion, is also shown for completeness. Additional details concerning the samples and the measuring technique can be found in Ref. .

In the following, temperatures are expressed in terms of the reduced temperature relevant for the statistical mechanics of the system, $`\tau k_BT/J(T)`$, where $`J(T)`$ is the temperature-dependent Josephson coupling energy. At temperatures well below the zero-field critical temperature $`\tau _{cN}`$, $`J(T)`$ was deduced from measurements of the ”bare” sheet kinetic inductance $`L(T)=(\varphi _0/2\pi )^2(5/3)^N/\sqrt{3}J(T)`$ of the unfrustrated samples . Extrapolation to higher temperatures was then achieved by fitting the low-temperature data to theoretical expressions  for $`J(T)`$.

Because of their two-dimensional (2D) nature at length scales larger than $`L_N`$, both samples are expected to exhibit, at zero frustration ($`f=0`$) and in the limit $`\omega 0`$, a Berezinskii-Kosterlitz-Thouless (BKT) phase transition . A sharp depression of $`L^1(0)`$, which can be associated with the BKT transition, has been indeed observed at, respectively, $`\tau _{c2}0.57`$ and $`\tau _{c4}0.23`$, in good agreement with the theoretical prediction  $`\tau _{cN}=(3/5)^N\tau _{c0}`$, where $`\tau _{c0}1.5`$ is the reduced temperature of the BKT transition of a regular triangular Josephson junction array  with the same $`J(T)`$.

## III The fractal regime

### A The ground states of a regular array of Sierpinski gaskets and their sheet inductance

We start by recalling that, within the framework of an approximation ignoring thermal fluctuations, a Josephson junction array behaves, with respect to an external (dc) current source, like a network of inductors {$`L_{ij}`$}, whose inductances are given  by:

$$L_{ij}(\theta _{ij})=\frac{(\varphi _0/2\pi )^2}{J\mathrm{cos}\theta _{ij}},$$
(1)
where $`\theta _{ij}`$ is the gauge-invariant phase difference across the link $`ij`$. As required by fluxoid quantization, the sum of $`\theta _{ij}`$ around a lattice cell is equal to $`2\pi (fSm)`$, where $`f`$ is the frustration parameter expressing the magnetic flux (in units of $`\varphi _0`$) threading an elementary triangular cell of a gasket \[$`f=Ba^2\sqrt{3}/(4\varphi _0)`$\], $`S`$ the area of the cell (expressed in units of the area of an elementary triangular loop) and $`m`$ the number of flux quanta (vortices) penetrating the cell under consideration.

In writing Eq. (1) we have assumed that the proximity-effect coupled SNS junctions forming the arrays studied in this work have a sinusoidal current-phase relation at the temperatures of interest . It clearly follows from Eq. (1) that, even if all the junctions are identical, their effective inductances in a frustrated system may differ substantially from each other on account of the nonuniform distribution of $`\{\theta _{ij}\}`$. Since the array magnetoinductance $`L(f)`$ can be found by applying Kirchhoff’s laws to the inductor network {$`L_{ij}`$}, it is evident that, at any frustration $`f`$, $`L(f)`$ will be completely determined once the distribution of $`\{\theta _{ij}\}`$ is known.

The ground state of a uniformly frustrated array of periodically repeated gaskets of order $`N`$ \[see Fig. 1(a)\] can be constructed by a simple juxtaposition after finding the ground state of an isolated $`N`$th-order gasket only if the constraints of fluxoid quantization imposed on the largest triangular loops (located between the $`N`$th-order gaskets) are automatically fulfilled. Recalling the definition of $`f`$ given above, it can be shown  that this condition is satisfied only for a particular set of frustrations given by:

$$f=\frac{M}{2\times 4^N},$$
(2)
where $`M`$ is an integer corresponding to the total number of vortices in the rhombohedral unit cell of the SG array composed, as shown in Fig. 1(b), by a single $`N`$th-order gasket and the adjacent ”empty” triangular loop. Thus, in order to calculate the magnetoinductance of the system at the values of $`f`$ given by Eq. (2), all we need is to determine the ground-state distribution of $`\{\theta _{ij}\}`$ in one of its constituent gaskets.

Assuming that the ground state of a gasket is the one having the highest possible symmetry consistent with its reflection and third-order rotation symmetries  (for wire networks with the geometry of a third-order SG this conjecture has been confirmed by numerical calculations ), it can be shown  that the number of independent bond variables $`\{\theta _{ij}\}`$ is equal to $`(3^N+1)/2`$. They can be found from an equal number of constraints imposed by current conservation at the nodes and fluxoid quantization in the loops of the gasket. For example, in the second-order gasket of Fig. 2 the links sharing the same symbol have the same values of $`\theta _{ij}`$, thereby showing that there are only 5 independent bond variables.

Imposing current conservation at one of the 6 equivalent nodes (denoted by $`K`$ in Fig. 2) and fluxoid quantization in the 4 non-equivalent loops $`a`$, $`b`$, $`c`$, and $`d`$ (with, respectively, vortex occupation numbers $`m_a`$, $`m_b`$, $`m_c`$, and $`m_d`$) leads to a system of 5 equations (only one of which, describing current conservation, is nonlinear) which have to be solved numerically under the constraint that the distribution of the quantum numbers $`\{m_\alpha \}`$ ($`\alpha =a,b,c,d`$) is such that the gasket energy:

$$E=J\underset{<ij>}{}(1\mathrm{cos}\theta _{ij})$$
(3)
is minimized. Quite remarkably, for a given frustration satisfying Eq. (2), the distribution of $`\{m_\alpha \}`$ corresponding to the lowest energy turns out to be identical  to that emerging from a calculation based on junctions with a linear current-phase relation, for which a fully analytical treatment is possible . The result is illustrated in Fig. 3, where the ground-state vortex configurations for the rhombohedral unit cell of a regular array of second-order gaskets are shown for $`M`$ in the range \[$`0,40`$\].

Inspection of Fig. 3 reveals characteristic features, which are valid for arbitrary gasket order. One first observes that, with increasing frustration, vortex nucleation spreads from the largest to smallest loops , a property reflecting the hierarchical character of the energy needed to create a vortex excitation . Next, one notices that the vortices penetrate the gaskets only for $`M>1`$, thereby implying that a SG array will exhibit fractal behavior only for $`f>1/(2\times 4^N)`$. One further recognizes that, since for $`M=1`$ the rhombohedral unit cell contains just one single vortex sitting in the largest triangular loop, the ground state of the array at $`f=1/(2\times 4^N)`$ corresponds to a triangular lattice of vortices with a nearest-neighbor distance equal to the gasket size $`L_N`$. Recalling that $`L_N`$ plays the role of $`\xi _p`$, one expects that for $`f<1/(2\times 4^N)`$ the system will behave like a regular Josephson junction array with conventional euclidian geometry. Thus, we identify

$$f_{cN}=\frac{1}{2\times 4^N}$$
(4)
as the frustration at which the EF crossover occurs. The nature of the ground states of the SG array in the euclidian regime ($`0<f<f_{cN}`$) is discussed in Sec. V.

Having shown how the structure of the ground state can be determined for the particular set of frustrations (2), we can now proceed with the calculation of the sheet magnetoinductance $`L(f)`$ of the SG array, the quantity measured in our experiments. We first notice that, for this particular set of frustrations, $`L(f)`$ is proportional to the magnetoinductance of a constituent gasket. Therefore, if we calculate the inductance of a single gasket and normalize it to its value at $`f=0`$ to eliminate the trivial dependence on the gasket size, we obtain a result also expressing the normalized sheet magnetoinductance $`L(f)/L(0)`$ of the composite periodic system. In Ref.  we have pointed out that for a given distribution of $`\{\theta _{ij}\}`$ the calculation of the inductance of a single gasket can be performed by successive application of the triangle-star transformation well known in the theory of electric networks . The result of such a calculation for a regular array of second-order gaskets is shown in Fig. 4 for multiples of $`f=1/32`$ in the range $`[11,43]`$. In order to compare this calculation with the experimental data presented in the following subsection, we plot the normalized inverse magnetoinductance $`L^1(f)/L^1(0)`$. Notice that, as expected for a uniformly frustrated Josephson junction array, $`L^1(f)`$ is symmetric with respect to $`f=1/2`$ and periodic with period $`f=1`$. We also recall that, although hardly visible in the linear plot of Fig. 4, in the fractal regime the frustration-induced correction $`\mathrm{\Delta }L(f)=L(f)L(0)`$ to the array inductance is predicted to scale , in the limit of small frustrations and low temperatures, as $`f^\nu `$ with $`\nu =\mathrm{ln}(125/33)/\mathrm{ln}40.96`$. Obviously, the power-law behavior of $`\mathrm{\Delta }L(f)`$ should no longer persist in the euclidian regime below $`f_{cN}`$.

So far we have considered only the frustrations given by Eq. (2), for which the ground state of a regular array of gaskets can be regarded as a periodic replication of the ground state of a single gasket. In order to determine the array ground state at rational frustrations differing from those given by Eq. (2), one should consider a ”supercell” comprising more than one gasket. The analysis of the ground states based on such supercells rapidly becomes cumbersome and (for $`f>f_{cN}`$) is not pursued in this work. However, as supercells imply that superconducting phase coherence extends at larger length scales, one can expect the corresponding ground states to be more vulnerable to thermal fluctuations \[and, consequently, less prominent in the fine structure of $`L(f)`$\] than those at $`f=Mf_{cN}`$. At these particular values of frustration the vortex configurations (shown in Fig. 3) are strongly pinned by the hierarchical potential landscape provided by the gaskets , thereby making these ground states particularly robust against thermal fluctuations.

### B Comparison with experiment and effects of hidden incommensurability

Focusing on the fine structure of $`L(f)`$ we now compare the theoretical predictions of Sec. III A with high-resolution magnetoinductance measurements performed on the array of second-order gaskets. Fig. 5(a) shows the normalized inverse sheet magnetoinductance $`L^1(f)/L^1(0)`$ of the array of second-order gaskets (measured at 1 kHz) at three different (reduced) temperatures. We first observe that, although the overall shape of the magnetoinductance curves looks roughly similar, the fine structure becomes richer and much sharper with increasing temperature, thereby revealing very clearly almost all the ”superfluid” peaks corresponding to the states with a unit cell not larger than a single gasket and belonging to the sequence $`f=M/32`$ given by Eq. (2). Notice that, to make the identification of the structures easier, the frustration unit on the horizontal axis of Fig. 5(a) is chosen to be equal to $`f_{cN}=1/32`$. The striking evolution of the fine structure with temperature suggests that the motion of vortices due to thermal fluctuations plays a major role in the dynamic response of these arrays. We interpret the behavior shown in Fig. 5(a) as clear evidence that, at sufficiently high temperatures, superconducting phase coherence in the neighborhood of the ground states at $`f=M/32`$, for which the vortex lattice is pinned, is drastically disrupted by vortex-lattice defects, created by excess or missing vortices, moving almost freely on the ”frozen” vortex background. This process dramatically sharpens the fine structure, thereby enhancing the amplitude of the oscillations. Similar behavior was also observed in experiments performed on wire networks of interconnected gaskets  and in regular triangular Josephson junction arrays  as well as in numerical simulations .

In sharp contrast to the theoretical prediction for a uniformly frustrated SG array (see Fig. 4), the $`L^1(f)/L^1(0)`$ curves of Fig. 5(a), although still symmetric with respect to $`f=0`$, are no longer periodic with period $`1`$, a behavior leading unavoidably to the suppression of the symmetry with respect to $`f=1/2`$. We attribute these anomalous features to the inhomogeneous frustration resulting from the change in the effective areas of different plaquettes caused by the asymmetric (with respect to the direction of the links) diamagnetic response of the truncated-star-shaped superconducting (Pb) islands \[see Fig. 1(a)\]. Because of this particular geometrical form, the screening currents flowing in these grains create a distortion of the current patterns associated with the individual loops which leads to a redistribution of the magnetic field and perturbs the self-similarity of the gaskets. In the temperature range ($`5.5K<T<6.4K`$) of the data shown in Fig. 5(a), the magnetic penetration depth of the Pb islands \[as estimated from the zero-temperature bulk Pb value ($`\lambda (0)40nm`$) and the (proximity-effect reduced) transition temperature ($`T_c=6.9K`$) of the Pb grains\] is at least 25 times smaller than their smallest planar geometrical dimension, which corresponds to the width of the junctions ($`2\mu m`$). Thus, one expects the distortion of the current patterns and, consequently, the nonuniformity of the frustration to have a considerable effect on $`L^1(f)`$.

The origin of the phenomenon being intimately related to the geometry of the superconducting islands rather than to the physical properties of the junctions, it seems plausible to describe the effect of the inhomogeneous frustration by changing merely the effective area of the different plaquettes according to a prescribed rule. This is illustrated in Fig. 6(a) for a second-order gasket, whose distortion is modeled by shifting the vertices of the triangular loops toward the ”centers of mass” of the corresponding superconducting islands, as shown in Fig. 6(b). Since the distribution of the screening currents in the islands is shifted in the same direction, this deformation scheme appears to be a reasonable approach offering, above all, the advantage of a simple description in terms of a single paramater, the (small) areal change $`ϵ`$ (expressed in units of the area of an elementary triangular cell) defined in Fig. 6(a). Notice, incidentally, that this procedure does not alter the frustrations \[Eq. (2)\] for which the ground state of the array can be constructed by a simple juxtaposition of independent gaskets.

As the resulting deformation preserves the reflection and third-order rotation symmetries of the gasket \[see Fig. 6(a)\], the determination of the vortex and $`\{\theta _{ij}\}`$ configurations in the highly symmetric ground states at $`f=Mf_{cN}`$ can be carried out by following again the procedure described in Sec. III A. The only modification appears in the formulation of the fluxoid-quantization constraints for the deformed loops, whose areas (and, therefore, the associated magnetic fluxes) change. For example, for the second-order gasket of Fig. 6(a) the areas of the deformed plaquettes $`a`$, $`b`$, $`c`$, and $`d`$ are equal, in an approximation linear in $`ϵ`$, to $`1ϵ`$, $`1(3/2)ϵ`$, $`13ϵ`$, and $`4+3ϵ`$, respectively. An analogous problem, however with a quite different deformation scheme, was considered by Ceccatto et al.  for a system with a linear current-phase relation.

Ground-state vortex configurations for the rhombohedral unit cell of an array of second-order gaskets with a deformation parameter $`ϵ=7.4\%`$ (the reason for the choice of this value is explained below) are shown in Fig. 7 for $`M`$ in the range \[$`0,40`$\]. A comparison with Fig. 3 reveals that some of the ground states of the undistorted array have been replaced by new ones. More precisely, while for the regular SG lattice the symmetry with respect to $`f=1/2`$ implies that the ground state for $`(1f)`$ may be conceived as resulting from the superposition of the ground states for $`f=1`$ and $`f`$, for the distorted system this property no longer holds. The absence of periodicity (corresponding to the $`f1+f`$ symmetry) and of the related $`f1f`$ symmetry are clearly reflected in the normalized inverse magnetoinductance shown in Fig. 5(b). This $`L^1(f)/L^1(0)`$ curve was calculated for an array of second-order gaskets using the method outlined in Sec. III A and was fitted to the low-temperature ($`\tau =0.05`$) data of Fig. 5(a) using $`ϵ`$ as an adjustable parameter. The best fit was obtained for $`ϵ=7.4\%`$. In this connection, it should be noticed that, although weak, thermal fluctuations still affect the data at $`\tau =0.05`$ and tend to enhance the amplitude of the oscillations with respect to that predicted by our calculation, which neglects thermal fluctuations. Nevertheless, the agreement is quite remarkable, especially if one considers that it involves only one fitting parameter. Moreover, the vertex displacement \[denoted by the arrow in Fig. 6(b)\] defining $`ϵ`$ turns out to be $`70\%`$ of the distance to the ”center of mass” of the truncated-star-shaped superconducting islands, thereby demonstrating the basic validity of our interpretation.

Similar behavior was observed in the magnetoinductance of the array of fourth-order gaskets, whose fine structure was found to be richer than that of the array of second-order gaskets, as demonstrated by the incipient splitting of some of the structures at $`f=M/32`$. However, because of the strong overlap resulting from the nonvanishing width of the superfluid peaks, only a fraction of the states at multiples of $`f=1/512`$ could be resolved and unambiguously identified in a plot at (relatively) low frustration resolution like that of Fig. 5 (for high-resolution data at very small $`f`$, see Fig. 9 in Sec. V).

Before closing this section, we would like to notice that when $`ϵ`$ is a rational number, the periodicity of $`L^1(f)`$ is restored, however with a period larger than 1.

## IV The crossover between the two regimes

Below $`f_{cN}`$ \[see Eq. (4)\] the SG arrays are expected to be in the euclidian (or homogeneous) regime. In order to provide preliminary evidence for the existence of the EF crossover at $`f_{cN}`$, in Fig. 8 we compare, in a log-log plot, the quantity $`\mathrm{\Delta }L^1(f)/L^1(0)[L^1(0)L^1(f)]/L^1(0)`$, which measures the relative change in superfluid density caused by frustration, for the two SG arrays studied in this work. Both curves were taken at 160 Hz and at temperatures such that the structures corresponding to the ordered states are emphasized by thermal fluctuations. With decreasing frustration the data for the sample of fourth-order gaskets exhibit, down to $`f_{c4}=1/512`$, clear fractal features: specifically, four self-similar stages (the number of stages being consistent with the order of the gaskets) reflecting the dilational symmetry of the gaskets and an overall scaling with $`f`$ which, in spite of the indisputable evidence for fluctuation effects, follows the asymptotic prediction $`\mathrm{\Delta }L^1(f)f^\nu `$. For comparison, the result of a calculation for an undistorted infinite gasket based on the methods discussed in Sec. III A is also shown in Fig. 8. Below $`f_{c4}`$ the data tend to flatten out, thereby signaling a possible change of regime. However, considering the fact that this change sets in almost at the limit of our inductance resolution ($`\mathrm{\Delta }L/L1\%`$), we refrain from drawing a firm conclusion with regard to the existence of a dimensional crossover in the array of fourth-order gaskets.

On this subject, the data for the array of second-order gaskets convey a much stronger message. In Sec. III B the magnetoinductance of this sample was shown to obey the predictions for the fractal regime whose signatures, although less pronounced than in the $`N=4`$ case, can also be identified, above $`f_{c2}=1/32`$, in the log-log plot of Fig. 8 \[notice that the self-similar and scaling properties of $`\mathrm{\Delta }L^1`$ are expected to become clearly manifest only in the asymptotic limit ($`f0`$) of higher-order gaskets \]. Below $`f_{c2}`$, however, the behavior of $`\mathrm{\Delta }L^1(f)/L^1(0)`$ drastically changes. Besides the loss of self-similarity and scaling, the data reveal, by closer inspection, the presence of structures \[the ”dips” in $`\mathrm{\Delta }L^1(f)`$\] corresponding to new commensurate states, which can not be ascribed to the fractal regime.

To strengthen the evidence for a regime crossover, in Fig. 9 we present the results of sheet inductance measurements performed at high frustration resolution in the range $`|f|f_{c2}`$. Once again, to promote structures corresponding to ordered states we tuned the effect of thermal fluctuations by increasing the temperature and reducing the measuring frequency (10 Hz) in an appropriate way. Moreover, in Fig. 9 we also include the dissipative component $`R(f)`$, whose remarkably well resolved fine structure provides additional evidence for the EF crossover. Notice that maxima in the (normalized) superfluid density, $`L^1(f)/L^1(0)`$, correspond to minima in $`R(f)`$, as it should be. For the array of fourth-order gaskets both the superfluid and the dissipative components display clear fractal features with marked structures at multiples of $`f_{c4}=1/512`$ corresponding to the frustration unit on the $`f`$-axis. In sharp contrast with this behavior, the dynamic response of the array of second-order gaskets is, within experimental accuracy, symmetric with respect to $`f=(1/2)f_{c2}`$ and contains well resolved structures at $`|f|=(1/4)f_{c2}`$, $`|f|=(1/3)f_{c2}`$ and $`|f|=(1/2)f_{c2}`$. In particular, the presence of prominent structures at $`|f|=(1/3)f_{c2}`$ provides compelling evidence for the existence of a new regime, since no commensurate state corresponding to this frustration can ever emerge from the characteristic sequence $`f=Mf_{cN}`$ of the fractal regime for gaskets of any order. We interpret these features as an unambiguous signature of the euclidian regime, in which the fine structure of the array’s impedance reflects, as shown in the following section, the existence of vortex ground states with a unit cell larger than a second-order gasket.

## V The euclidian regime

In order to understand how our arrays of gaskets behave in the euclidian regime, we first recall that the energy of a vortex in a Sierpinski gasket decreases with increasing size of the cell in which the vortex core is localized . Accordingly, a single vortex in a lattice of periodically repeated SGs can be considered as interacting with an external potential whose minima coincide with the centers of the largest triangular cells located between the $`N`$th-order gaskets and, therefore, form a triangular lattice. For $`f<f_{cN}`$ the number of vortices in the system is smaller than the number of the largest triangular cells, so that it is energetically favorable for the vortices to occupy only these largest cells, the concentration of the occupied ones being equal to the ratio $`f/f_{cN}`$. Therefore, one can expect that, for $`f<f_{cN}`$, the behavior of an array of periodically repeated SGs resembles that of a uniformly frustrated $`XY`$-model on a honeycomb lattice (with a reduced value of frustration $`f_H=f/f_{cN}<1`$), whose ground states can also be thought of as formed by vortices occupying the sites of a triangular lattice with the same concentration $`f_H`$.

Another approach leading to the same conclusion relies on the iterative procedure, based on successive applications of the triangle-star transformation , developed in Ref. . At low frustrations, such that there are no vortices in the gaskets, this method allows to replace each gasket by an elementary ”star” consisting of three links with a modified interaction. This transforms a lattice of periodically repeated SGs into a honeycomb lattice (see Fig. 10). The difference with respect to the conventional $`XY`$-model on a honeycomb lattice is that in the iteration process the interaction becomes almost harmonic , the anharmonic corrections becoming smaller and smaller as the number of iterations, which is set by the gasket order $`N`$, increases. Moreover, for $`f0`$ the resulting effective interaction is no longer an even function of $`\theta _{ij}`$ .

After establishing that the ground states in the euclidian regime ($`f<f_{cN}`$) are formed by vortices occupying the sites of a triangular lattice with concentration $`f_H=f/f_{cN}<1`$, it is natural to expect that the states which are particularly robust against thermal fluctuations correspond to values of $`f_H`$ allowing the formation of an undistorted triangular vortex lattice (analogous to the Abrikosov lattice in bulk superconductors) commensurate with the underlying lattice provided by the array. The energy required to create a defect in these highly symmetric states is larger than for frustrations requiring the vortex lattice to be distorted. Therefore, these Abrikosov-like states will be less vulnerable to thermal fluctuations. For $`1/2<f_H<1`$ almost equally stable states can be constructed when the vacancies in the densely packed triangular vortex lattice corresponding to $`f_H=1`$ also form an undistorted triangular lattice. This leads to the symmetry $`f_H1f_H`$, although it is not rigorous.

It is readily seen that triangular vortex lattices commensurate with the underlying lattice can be constructed for $`f_H=1/(m^2+mn+n^2)`$, where $`m`$ and $`n`$ are integers ($`m1`$, $`0nm`$). In particular, $`m=n=1`$ gives $`f_H=1/3`$, $`m=2`$ and $`n=0`$ give $`f_H=1/4`$, $`m=2`$ and $`n=1`$ give $`f_H=1/7`$, etc. The states for $`f_H=1/3`$ and $`f_H=1/4`$ are shown, respectively, in Figs. 11(a) and 11(b). Within this family, the most ”dense” state, i.e. the state for $`f_H=1/3`$, can be expected to be, on account of its stronger vortex-vortex interaction, the most stable one, followed, in order of decreasing stability, by those for $`f_H=1/4`$ and $`f_H=1/7`$. This is precisely what we observe in the magnetoimpedance data of Fig. 9 for the array of second-order gaskets where, focusing on the interval $`|f_H|<1/2`$), we find that the most prominent structures, signaling a very stable commensurate state, appear at $`|f_H|=1/3`$. Weaker structures corresponding to the next-stable triangular vortex lattice are also well resolved at $`|f_H|=1/4`$, whereas structures at $`|f_H|=1/7`$ are barely visible and present only in the inverse magnetoinductance data. It is worth mentioning that the observation of the particularly prominent vortex state at $`|f_H|=1/3`$ is entirely consistent with theoretical predictions  as well as with results of Monte Carlo simulations  for the frustrated $`XY`$ model on the honeycomb lattice.

The magnetoimpedance data of Fig. 9 show weak structures also at $`|f_H|=1/2`$. The ground states of the $`XY`$-model on a honeycomb lattice at this particular frustration and the corresponding vortex configurations were studied in Ref.  and shown to possess a so-called accidental (i.e., not related to the symmetry) degeneracy, which can be discussed in terms of zero-energy domain walls. Fig. 11(c) shows an example of a ground state for $`f_H=1/2`$ characterized by a regular 1D superlattice structure. In this state all the variables $`\{\theta _{ij}\}`$ take the values $`\theta _{ij}=0,\pm \pi /4`$. Redistributing the same set of variables in a different way among the bonds of the lattice allows to transform the state of Fig. 11(c) into another one with the same energy, shown in Fig. 11(d), where a domain wall separates two ”striped” states of the type shown in Fig. 11(c). An infinite family of states with the same energy can be constructed by creating sequences of such zero-energy domain walls parallel to each other. Fig. 11(e) shows another example of a periodic ground state, which can be obtained by introducing the largest possible number of domain walls into the state shown in Fig. 11(c).

It is known  that in systems allowing the formation of zero-energy domain walls the stabilization of long-range order at nonzero temperatures is achieved if the accidental degeneracy is removed by thermal fluctuations. This mechanism (the so-called ”order by disorder” ) is relatively weak and, therefore, the ordering in frustrated $`XY`$-models with accidental degeneracy should be less stable (and destroyed at lower temperatures) than in other $`XY`$ models. This effect should be even more pronounced for the fully frustrated (i.e. for $`f_H=1/2`$) model on a honeycomb lattice, which, in contrast to other 2D $`XY`$-models with accidental degeneracy , does not display any difference in the spin-wave free energy, calculated in the harmonic approximation, of different periodic ground states . The Monte-Carlo simulations of Shih and Stroud  have indeed demonstrated that the phase transition of the frustrated $`XY`$-model on a honeycomb lattice at $`f_H=1/2`$ takes place at a much lower temperature than at $`f_H=1/3`$ or $`f_H=1/4`$, where the accidental degeneracy is absent.

Returning to the system of periodically repeated SGs at $`f=f_{cN}/2`$ ($`f_H=1/2`$), it is also possible to start the search for its ground state by finding the structure of the state which minimizes the energy (3) for the periodic vortex configuration shown in Fig. 11(c). Obviously, this state is periodic, and its unit cell comprises two gaskets. By redistributing the same set of $`\{\theta _{ij}\}`$ among the bonds of the lattice in a different way, one can construct the state with the same energy corresponding to the vortex configuration shown in Fig. 11(d). However, in contrast to the fully frustrated $`XY`$-model on a honeycomb lattice, the state obtained in this way will not be an extremum of the Hamiltonian and, therefore, a slight readjustment of $`\{\theta _{ij}\}`$ can further decrease its energy. This means that in a system of periodically repeated SGs the domain wall shown in Fig. 11(d) has a negative energy. Accordingly, the state with the lowest energy corresponds to the periodic vortex configuration shown in Fig. 11(e), which is characterized by the highest possible density of such domain walls. The unit cell of this state comprises four gaskets. In the fully frustrated $`XY`$-model on a honeycomb lattice a ground state with the same structure is selected if one takes into account the interaction, of arbitrary sign, with the second-nearest neighbors .

Thus, the phase transition taking place, with decreasing temperature, at $`f_H=1/2`$ should be related to the appearance of long-range order corresponding to the vortex configuration shown in Fig. 11(e). The selection of this state relies on a weak mechanism, which loses its efficiency with increasing $`N`$, since the effective interaction becomes almost harmonic under decimation. One can therefore expect the ordered phase at $`f_H=1/2`$ to be again rather vulnerable to thermal fluctuations. This explains why the structures at $`|f_H|=1/2`$ in Fig. 9 are much weaker  than those at $`|f_H|=1/3`$. Their strength is at most comparable to that of the structures corresponding to the triangular vortex lattice at $`f_H=1/4`$, whose weaker vortex-vortex interaction makes the ordering less robust than at $`f_H=1/3`$. Numerical simulations  show that in the conventional $`XY`$ model on a honeycomb lattice the ordered phases appear, with decreasing temperature, in the same order: first at $`f_H=1/3`$, then at $`f_H=1/4`$, and only further down at $`f_H=1/2`$.

## VI Conclusion

Our magnetoinductance measurements on Josephson junction arrays of periodically repeated Sierpinski gaskets have clearly demonstrated the existence of two regimes. In agreement with theoretical analysis, in one of them (the fractal regime) the peaks observed in the inverse sheet magnetoinductance, reflecting those states which are the most stable against thermal fluctuations, are equally spaced. Neighboring states in this sequence differ from each other by the penetration of an additional vortex into each unit cell of the array. In the other, euclidian, regime the sequence of the observed stable states corresponds to periodic lattices of vortices occupying the largest triangular cells of the array.

The agreement with theory is achieved not only for the positions of the different peaks, but also for their relative strengths. In the euclidian regime the structure of the ordered states is analogous to that on a honeycomb lattice with reduced frustration and the relative stability of the different states can be understood in terms of vortex lattice disordering. However, the amplitudes of the peaks observed in the fractal regime can be quantitatively explained only if the redistribution of the magnetic field in the array due to the asymmetric shape of the superconducting islands is taken into account. In this connection, we would like to mention the recent work by Park and Huse , who compared the energies of different states in a wire network with a kagomé lattice geometry at full frustration. These authors came to the conclusion that the effects related to the finite width of the wires can be compensated by bending the wires. This is equivalent to our conjecture that the influence of the asymmetry associated with the screening currents can be reduced to a redistribution of the magnetic field with respect to an ideal system.

The results of our magnitoinductance measurements show that this phenomenon is more pronounced at high frustration levels, an observation consistent with the analysis presented in Sec. III B. Indeed, in our model, the triangular cells exhibiting the largest relative areal changes are the smallest (elementary) ones, which are precisely those providing the dominant contribution to $`L^1(f)`$ at high values of $`f`$. Data taken at very small frustrations, like those shown in Fig. 9, are practically unaffected by the nonuniform frustration resulting from the asymmetry of the screening currents.

Early studies  of a large-order ($`N=10`$) single gasket of superconducting aluminum wires revealed that the period associated with adjacent minima of the superconducting-to-normal phase boundary $`T_c(f)`$ was larger than that extracted from adjacent maxima. We have found a similar perturbation of the periodic field dependence in our calculations of $`L^1(f)`$ based on the model for the magnetic field redistribution proposed in Sec. III B. The same mechanism may be responsible for the anomalous feature observed in Ref. . It should be noticed, however, that the experiments of Ref.  were performed at temperatures very close to $`T_c(0)`$, where the estimated magnetic penetration depth of the Al wires turns out to be comparable to their width. Accordingly, the effective gasket distortion resulting from screening effects is expected to be weaker than in our arrays. A quantitative verification is therefore needed before drawing a conclusion as for the ability of our model to explain the anomaly observed in Ref. .

The influence of the incommensurability of different lattice cells on the magnetoresistance of Josephson junction arrays was discussed by Kosterlitz and Granato  in relation to experiments performed on periodic arrays with a complex unit cell . However, quantitative agreement between experiment and theory in treating incommensurability effects in Josephson junction arrays has been demonstrated only by the present work . In our system the phenomenon of incommensurability is present in an hidden form and its manifestations appear as a consequence of the asymmetric shape of the superconducting islands.

###### Acknowledgements.

We are grateful to H. Beck for several interesting discussions. We also thank R. Théron for his assistance in preparing the figures. P.M. would like to thank Ø. Fischer and J.M. Triscone for the hospitality extended to him during an eight-month visit at the DPMC of the University of Geneva, where part of this article was written. This research was supported by the Swiss National Science Foundation, the Swiss Federal Office for Education and Science within the framework of the TMR-project ”Superconducting Nanocircuits” of the European Union, and the ”Vortex Scientific Program” of the European Science Foundation. S.E.K. additionally acknowledges the support of the Program ”Quantum Macrophysics” of the Russian Academy of Sciences and of the Program ”Scientific Schools of the Russian Federation” (grant No. 00-15-96747).


format: text

text:
# The Phase Diagram of Four Flavor SU(2) Lattice Gauge Theory at Nonzero Chemical Potential and Temperature

## I Introduction

Recently there has been a resurgence of interest in QCD at nonzero chemical potential for quark number. Arguments based on instantons  and phenomenological gap equations  support the old idea that diquark condensation and a color superconductivity phase transition  will occur at a critical chemical potential somewhat greater than that at which nuclear matter forms. Nuclear matter is expected to be the ground state of QCD above a chemical potential $`\mu _N`$, which is slightly less than one-third of the proton’s mass. These arguments neglect the forces of confinement and may only be reliable at asymptotically large chemical potentials. Since diquarks carry color in real QCD, a quantitative understanding of confinement and screening is needed to estimate the energies of the states which control the phases of the system. Unfortunately, brute force lattice simulation methods, which are so useful at $`T0`$ and $`\mu =0`$, do not yet provide a reliable simulation algorithm for these environments in which the fermion determinant becomes complex.

Given these problems, theorists have turned to simpler models which can be analyzed and simulated by current methods. One of the more interesting is the color $`SU(2)`$ version of QCD which addresses some of the issues of interest , , . In this model diquarks do not carry color, so their condensation does not break color symmetry dynamically. Instead, the diquark condensed phase resembles a superfluid rather than a superconductor. The critical chemical potential is also expected to be one-half the mass of the lightest meson, the pion, because quarks and anti-quarks reside in equivalent representations of the $`SU(2)`$ color group. Hence hadron flavor multiplets include both mesons and diquarks. Chiral Lagrangians can be used to study the diquark condensation transition in this model because the critical chemical potential vanishes in the chiral limit, and the model has a Goldstone realization of the spontaneously broken quark-number symmetry . The problem has also been studied within a Random Matrix Model at non-zero $`\mu `$ and $`T`$ . Lattice simulations of the model are also possible because the fermion determinant is real and non-negative for all chemical potentials. One hopes that these developments will uncover generic phenomena that will also apply to QCD at nonzero chemical potential. Furthermore, most of the theories that aim at a description of the phase diagram of 3-color QCD at nonzero baryon chemical potential and temperature can be implemented to describe the phase diagram of 2-color QCD at nonzero baryon chemical potential and temperature. Therefore, an important motivation to conduct Lattice simulations is to check the validity of these theories in the 2-color case. If a theory is not correct in the two-color case, it is very doubtful that the predictions of the corresponding 3-color theory can be trusted.

Preliminary lattice simulations of the $`SU(2)`$ model with four species of quarks, simulation data and an Effective Lagrangian analysis of aspects of the $`T`$-$`\mu `$ phase diagram were recently published . It is the purpose of this article to present the continuation of that work on larger lattices, closer to the theory’s continuum limit with a focus on the system’s phase diagram at nonzero $`\mu `$ and $`T`$. Very early work on this model at finite $`T`$ and $`\mu `$ was performed by . A simulation study of the spectroscopy of the light bosonic modes will be presented elsewhere .

In our exploratory study  of this model’s phase diagram, we found an unanticipated result. At fixed but small quark mass, which insures us that the pion has a nonzero mass and chiral symmetry is explicitly broken, we found that at high $`T`$ and $`\mu `$ there is a line of transitions separating a quark-gluon plasma phase from one with a diquark condensate. Along this line there is a tricritical point, labeled $`2`$ in Fig. 1, where the transition switches from being second order (and well described by mean field theory) at relatively low $`\mu `$ to a first order transition at an intermediate $`\mu `$ value. We will present evidence here for metastability along the line of high $`\mu `$ and $`T`$. We have argued elsewhere  that this simulation result, the existence of a tricritical point $`2`$, has a natural explanation in the context of chiral Lagrangians. Following the formalism of  we argued that trilinear couplings among the low lying boson fields of the Lagrangian become more significant as $`\mu `$ and $`T`$ increase and they can cause the transition to become first order at a $`\mu `$ value in the vicinity of the results found in the simulation. Using hindsight, this behavior should not have come as a surprise. $`\mu `$ plays the rôle of a second ‘temperature’ in this theory in that it is a parameter that controls diquark condensation, but does not explicitly break the (quark-number) symmetry. The existence of two competing ‘temperatures’ is a characteristic of systems which exhibit tricritical behavior. A tricritical point is also found in Chiral Perturbation Theory, as well as in a new Random Matrix Theory model .

The conjectured phase diagram of Fig. 1 also has a critical point labeled $`1`$ which connects the line of first order transitions to the dashed line of crossover phenomena extending to the $`\mu =0`$ axis where we expect a finite T chiral transition between conventional hadronic matter and the quark-gluon plasma. For sufficiently light quarks in the four flavor SU(2) model, the $`\mu =0`$ transition is known to be first order both from theoretical arguments  and simulations . For quark mass $`m=0.05`$ we find that the transition is smoothed out to a crossover region as $`\mu `$ is increased until the dynamics favors diquark condensation. Clearly, the results will depend on the precise value of the quark mass $`m`$. At $`m=0.05`$ we find that the transition at $`\mu =0.10`$ is as smooth as that at $`\mu =0.0`$, while the transition at $`\mu =0.20`$ is noticeably sharper. This result suggests that if $`m`$ were chosen somewhat smaller than $`m=0.05`$, the scenario with the critical point 1 and the first order line might emerge. However, this critical point 1 might be absent so the crossover line would intersect the curve separating the diquark condensation phase from the normal phase. We will be pursuing this suggestion in separate simulations. The simulations done here are not accurate enough to probe the details of the region of the phase diagram Fig. 1 where the two first order lines split apart. Simulations at lower quark masses on larger lattices might be needed.

In QCD at nonzero baryon chemical potential Fodor and Katz  have estimated the position of the critical point 1 using a modified, multivariable Glasgow algorithm . It will be interesting to find point 1 in the four flavor SU(2) model at smaller $`m`$ values using our conventional simulation methods and then checking whether the Fodor-Katz method can find it with competitive accuracy. If point 1 is absent, the Fodor-Katz method should be able to find the intersection of the crossover line with the curve delimiting the diquark phase.

We also present extensive and accurate simulation results on a $`16^4`$ lattice, essentially at vanishing $`T`$. The critical indices of the transition at nonzero $`\mu `$ are measured and favor mean-field scaling in agreement with Chiral Perturbation Theory, analyzed through one loop corrections, and simulations of 2-color QCD in the strong coupling limit . The critical indices $`\beta _{mag}`$ and $`\delta `$ are emphasized. Equation of State fits, suggested by Chiral Perturbation Theory, also support these conclusions. As a third method to determine the critical indices, we use the static scaling hypothesis. This analysis is also compatible with mean field, although this approach is unable to predict the index $`\delta `$ accurately because of its sensitivity to estimates of the critical chemical potential, $`\mu _c`$. Unfortunately, because simulations at small $`\lambda `$ are extremely costly in computing requirements, we have not been able to approach the critical point closely enough to make definitive measurements, insensitive to the detailed nature of the critical scaling. In a recent paper in which we study quenched versions of this and closely related theories, we have indicated how sensitive measurements of critical scaling are to the form of the scaling variable . This paper does indicate, however, that there are classes of scaling variables for which there exist extended regions outside of the true scaling window, where the order parameter scales with the correct critical exponents.

## II Simulation Method and The Algorithm

The lattice action of the staggered fermion version of this theory is:

$$S_f=\underset{sites}{}\left\{\overline{\chi }[D/(\mu )+m]\chi +\frac{1}{2}\lambda [\chi ^T\tau _2\chi +\overline{\chi }\tau _2\overline{\chi }^T]\right\}$$
(1)
where the chemical potential $`\mu `$ is introduced by multiplying links in the $`+t`$ direction by $`e^\mu `$ and those in the $`t`$ direction by $`e^\mu `$ . The diquark source term (Majorana mass term) is added to allow us to observe spontaneous breakdown of quark-number on a finite lattice. The parameter $`\lambda `$ and the usual mass term $`m`$ control the amount of explicit symmetry breaking in the lattice action. We will be particularly interested in small values of $`\lambda `$ and the extrapolation $`\lambda 0`$ for a given $`m`$ to produce an interesting, realistic physical situation. This is the system that has been studied analytically using effective Lagrangians at non-zero chemical potential $`\mu `$ , and non-zero temperature .

Integrating out the fermion fields in Eq.1 gives:

$$pfaffian\left[\begin{array}{cc}\lambda \tau _2& 𝒜\\ 𝒜^T& \lambda \tau _2\end{array}\right]=\sqrt{\mathrm{det}(𝒜^{}𝒜+\lambda ^2)}$$
(2)
where

$$𝒜D/(\mu )+m$$
(3)
Note that the pfaffian is strictly positive, so that we can use the hybrid molecular dynamics  method to simulate this theory using “noisy” fermions to take the square root, giving $`N_f=4`$.

For $`\lambda =0`$, $`m0`$, $`\mu 0`$ we expect no spontaneous symmetry breaking for small $`\mu `$. For $`\mu `$ large enough ($`\mu >m_\pi /2`$ according to most approaches including Chiral Perturbation Theory ) we expect spontaneous breakdown of quark number and one Goldstone boson – a scalar diquark. (The reader should consult  for a full discussion of the symmetries of the lattice action, remarks about spectroscopy, Goldstone as well as pseudo-Goldstone bosons, and for early simulations of the 8 flavor theory at $`\lambda =0`$.)

## III Simulation Results and Analysis

We now consider the simulation results for the $`N_f=4`$ theory on $`16^4`$, $`8^4`$ and $`12^3\times 6`$ lattices, measuring the chiral ($`\overline{\chi }\chi `$) and diquark condensates ($`\chi ^T\tau _2\chi `$), the fermion number density, the Wilson/Polyakov line, etc.

Let us first briefly mention our $`8^4`$ simulations. These were performed at $`\beta =1.5`$, a relatively strong coupling, and $`m=0.025`$. These extend our earlier work with the same $`\beta `$ but $`m=0.1`$. Using a smaller quark mass increases the validity of chiral perturbation analyses. In addition, we are performing simulations on a $`12^3\times 24`$ lattice with the same parameters to examine the spectrum of Goldstone and pseudo-Goldstone bosons of this theory. For this purpose it is necessary to have the pion mass significantly lighter than that of non-Goldstone, non-pseudo-Goldstone excitations. We have simulated at $`\lambda =0.0025`$, and $`0.005`$. Our measurements of the diquark condensate are given in figure 2. Each ‘data’ point represents a 2000 molecular-dynamics time-unit simulation. The fact that a simple linear extrapolation of our measurements to $`\lambda =0`$ is very small, not only at $`\mu =0`$, where we know that the condensate vanishes, but also at $`\mu =0.1`$ leads us to believe that the condensate vanishes in this limit, probably up to $`\mu 0.2`$. (Note that is is known that there is some small but finite value of $`\mu `$, below which this condensate vanishes .) For $`\mu `$ much greater than $`0.2`$ but less than the saturation value it is reasonably clear that the condensate does not vanish in the limit as $`\lambda 0`$ (on an infinite lattice). Hence there is a phase transition for some value of $`\mu `$ to a state where quark number is spontaneously broken by a diquark condensate. The transition occurs over a relatively small range of $`\mu `$, but does not appear to be first order. Because of the relatively small scaling window this implies, we will not be able to analyze critical scaling until we perform measurements at more $`\mu `$ values within this scaling window.

The chiral condensate remains approximately constant up to $`\mu =\mu _c`$ after which it falls rapidly approaching zero for large $`\mu `$. The fermion density rises slowly and approximately linearly from zero for $`\mu \mu _c`$. Above $`\mu 0.6`$, it increases more rapidly approaching its saturation value of $`2`$ at $`\mu 1.0`$.

We now turn to measurements on a $`16^4`$, ‘zero temperature’ lattice. We simulated the $`SU(2)`$ model at a relatively weak coupling $`\beta =1.85`$, within the theory’s apparent scaling window, in order to attempt to make contact with the theory’s continuum limit. The quark mass was $`m=0.05`$, as in , and a series of simulations were performed at $`\lambda =0.0025`$, $`0.005`$, and $`0.01`$ so that our results could be extrapolated to vanishing diquark source, $`\lambda =0`$. Simple linear and spline extrapolations were used and these procedures appeared to be sensible away from the transition. None of the conclusions to be drawn here will depend strongly on the limit. In fact, everything we say could be gathered from our data at a fixed (small) $`\lambda `$ value, $`0.005`$ or $`0.0025`$, say. However, only in the limit of vanishing $`\lambda `$ do we expect real diquark phase transitions (at least when the transitions are second order), so it is particularly interesting and relevant to investigate the $`\lambda 0`$ limit. The quantitative agreement between fits at fixed, small $`\lambda `$ values and those using extrapolated data enforce the prevalent idea that simple extrapolation procedures are reliable except in the immediate vicinity of the critical point where non-analyticities are inevitable.

Because we are forced by practical considerations to work at $`\lambda `$ values large enough that the condensate is relatively smooth over the transition region, we have to deduce critical scaling from the curvature of the condensate further from $`\mu _c`$ than we would like. Thus we run the risk that we are outside the scaling window and the apparent scaling we are seeing does not represent true critical scaling. However, we saw for quenched theories  at relatively weak couplings, that the form of the natural scaling variable is such that it has an inflection point when expressed in terms of the simplest scaling variable $`\mu \mu _c`$. This leads to an extended domain where the order parameter scales with the true critical exponents. It is reasonable to assume that dynamical theories behave similarly.

In Table I we present the raw data on the $`16^4`$ lattice at the three $`\lambda `$ values. Table II exhibits the results of extrapolating this data to $`\lambda =0`$. Typically, one thousand time units of the Hybrid Molecular Dynamics algorithm were accumulated to generate these data sets. The error bars in the data sets of Table I account for the correlations in the raw data sets. The low $`\lambda `$ runs were our most CPU intensive.

Since the Hybrid algorithm suffers from systematic errors proportional to the square of the discrete time step used in integrating its stochastic differential equations forward in Monte Carlo time , we were forced to a small discrete time step. The simulations reported here used time steps as small as $`dt=0.0025`$ to control these errors. Even then, because our statistical errors are so small, these systematic errors might not be negligible. We will later indicate cases where we used $`dt`$ values which were larger than optimal and have seen signs of such errors.

Now let us discuss some results. In Fig. 3 we show the diquark condensate linearly extrapolated to $`\lambda =0`$ plotted against the chemical potential $`\mu `$. The quark mass $`m`$ was fixed at $`0.05`$, the coupling was $`\beta =1.85`$ and the linear extrapolation used the raw data at $`\lambda =0.0025`$ and $`0.005`$.

We see evidence to a quark-number violating second order phase transition in this figure. The dashed line is a power law fit from $`\mu =0.30`$ to $`0.35`$ which predicts the critical chemical potential of $`\mu _c=0.2860(2)`$. The power law fit is good, its confidence level is $`48`$ percent, and its critical index is $`\beta _{mag}=0.54(3)`$ which is consistent with the mean field result $`\beta _{mag}=1/2`$, predicted by chiral perturbation theory ,  including one loop corrections. Note that the quark mass is fixed at $`m=0.05`$ throughout this simulation, so chiral symmetry is explicitly broken and the transition of Fig. 2 is due to quark number breaking alone.

In figure 4 we show the diquark data extrapolated to $`\lambda =0`$ quadratically using the raw data at $`\lambda =0.010`$, $`0.005`$ and $`0.0025`$. The fit to the data at $`\mu `$ ranging from $`0.30`$ through $`0.35`$ is also shown. The fit predicts a critical point $`\mu _c=0.2870(2)`$, has a confidence level of $`47`$ percent and has a critical exponent of $`\beta _{mag}=0.54(4)`$, again in agreement with mean field theory.

Further evidence for $`\beta _{mag}=1/2`$ can be obtained directly from the raw data without any extrapolations. For example, in figure 5 we show the $`\lambda =0.0025`$ data and a power law fit ranging form $`\mu =0.30`$ through $`0.38`$. We find $`\mu _c=0.2749(2)`$, $`\beta _{mag}=0.54(5)`$ with a confidence level of $`24`$ percent.

We also considered extrapolations of the form $`\chi ^T\tau _2\chi =\rho `$ where $`\rho `$ is the solution of an equation of the form $`\rho ^3=a\rho +b\lambda `$ with $`a`$ and $`b`$ determined as functions of $`\mu `$ by performing an independent fit at each $`\mu `$. This form is suggested my mean-field scaling forms, which would predict that $`a`$ should be proportional to $`\mu \mu _c`$ close to the transition. However, these fits were so poor that this approach was abandoned.

In addition, we considered the fermion number density which is predicted to rise linearly from $`\mu _c`$ and found good agreement with this expectation, as shown in figure 6.

The fit predicts the critical point $`\mu _c=0.2765(5)`$, with a critical index $`1.17(7)`$ and a confidence level of $`13.1`$ percent. Plots and fits of the fermion number density at the other $`\lambda `$ values were equally good. In fact the fermion number density shows relatively little $`\lambda `$ dependence. This is in agreement with Chiral Perturbation Theory .

Finally, we considered the $`\lambda =0.005`$ simulations and extended the measurements up to large $`\mu 1.2`$, shown in figure 7. We see the curious phenomenon observed in the past: as $`\mu `$ becomes large, the condensate falls to zero. The point is that once $`\mu `$ becomes of order unity, the discreteness of the lattice comes into play and suppresses the density of states which forces the condensate to fall toward zero. The quark-number density approaches 2 (per lattice site), the maximum allowed by Fermi statistics. This is a lattice artifact. Nonetheless, we can try a power law fit to the diquark condensate over the entire region of $`\mu `$ from $`0.30`$ to $`0.60`$ where the curve is increasing. The result is shown in the figure and the fit predicts $`\beta _{mag}=0.49(7)`$ and $`\mu _c=0.2651(2)`$. Since $`\chi ^2/dof5`$, the fit is poor, but the deviations from this fit are small enough that this fit should be considered a reasonable zeroth order approximation.

In summary, $`\beta _{mag}`$ measurements are relatively easy and self-consistent. For all the parameters and procedures we tried, we found good agreement with mean field theory. It was important, however, to be working at relatively weak coupling, $`\beta =1.85`$, in the theory’s scaling window where such an extended scaling window exists. Experience has shown that at stronger coupling, the true scaling is described by effective Lagrangians closer to the non-linear sigma model form, where the scaling window in $`\mu \mu _c`$ is small and the extended scaling regime exhibits scaling in $`\mu \mu _c`$ but with an apparent critical exponent $`\beta _m`$ which is half the true value .

A second critical index in which we are interested is $`\delta `$, which controls the response of the system to symmetry breaking fields at the critical point $`\mu _c`$,

$$\chi ^T\tau _2\chi =A\lambda ^{1/\delta }$$
(4)
in the limit of small $`\lambda `$. Mean field theory predicts $`\delta =3`$. This scaling law Eq. 4 is usually difficult to use in practice because first it requires accurate estimates of $`\mu _c`$, and second it generally suffers from a small scaling window in $`\lambda `$. Not having a precise estimate of $`\mu _c`$ is probably our biggest problem, since as we have shown in , fits to the extended scaling regime can underestimate $`\mu _c`$ by $`8\%`$. In order to obtain some useful data, however, we simulated the $`16^4`$ lattice for four estimates of $`\mu _c`$ each over a range of $`\lambda `$ from $`0.002`$ to $`0.010`$. First we note that the chosen $`dt`$ for the $`\lambda =0.002`$ simulations was too large and although plotted on our graphs, these points were excluded from our fits.

We considered $`0.2855`$, $`0.2870`$, $`0.2920`$ and $`0.2970`$ as estimates of $`\mu _c`$ and tried fits of the form $`A\lambda ^{1/\delta }+B`$. The data and fits are shown in the next four figures, which we discuss in turn.

The $`\mu =0.2855`$ data predicted $`\delta =4.0(6)`$ with a confidence level of $`91`$ percent, but with a large, negative constant $`B=0.17(7)`$. The fit considers the data from $`\lambda `$ ranging from $`0.004`$ to $`0.010`$ because the data at $`\lambda =0.002`$ has large $`dt^2`$ errors. These results are shown in figure 8 The fact that $`B`$ is large and negative suggests that this estimate for $`\mu _c`$ is too low.

In figure 9 we consider the same ideas for the estimate $`\mu _c=0.2870`$ and now find $`\delta =2.5(2)`$ and $`B=0.02(2)`$, but with a rather poor confidence level of a few percent. For $`\mu _c=0.2920`$ (figure 10), we find $`\delta =2.1(2)`$, but $`B`$ is distinctly different from zero, $`B=0.05(2)`$. The confidence level is, however, very good, $`96`$ percent. Finally, for $`\mu _c=0.2970`$ we find $`\delta =2.7(6)`$ and $`B=0.02(4)`$ with a confidence level of about one percent. See figure 11.

In summary, $`\delta `$ is in the mean field ‘ballpark’, but the uncertainly in the simulation results is large. The results give us no reason to doubt the application of mean field theory because they are all in the vicinity of $`\delta =3`$ with some fits higher and some lower. The sensitivity of the simulation estimates to $`\mu _c`$ is considerable. We learn from these measurements that only a focused set of measurements that determine $`\mu _c`$ very well will produce a useful estimate of $`\delta `$.

## IV Comparison with Chiral Perturbation Theory

It is interesting to compare the lattice data to the theoretical predictions of Chiral Perturbation Theory. Two-color QCD at non-zero baryon/quark-number chemical potential has been thoroughly studied at the tree level and at next-to-leading order within Chiral Perturbation Theory . At next-to-leading order, it was found that the phase transition between the normal phase and the superfluid diquark phase is second order, that the critical chemical potential is given by $`\mu _c=m_\pi /2`$ and that the critical exponents are given by mean-field theory. Furthermore, based on the structure of the loop corrections, it was conjectured that the critical exponents are given by mean field at any (finite) order in Chiral Perturbation Theory .

Overall, the next-to-leading-order corrections do not dramatically change the leading-order results . The corrections turn out to be small, as long as $`\mu `$ is close to $`\mu _c`$. Therefore it should be enough to compare the lattice data to Chiral Perturbation Theory at leading-order, at least for the lattices sizes, couplings and parameters considered here.

In Chiral Perturbation Theory at leading order, the diquark condensate at non-zero diquark source is given by

$$\chi ^T\tau _2\chi =\overline{\chi }\chi _0\mathrm{sin}\alpha ,$$
(5)
where $`\overline{\chi }\chi _0`$ is the chiral condensate at zero chemical potential, zero temperature, and zero diquark source, and $`\alpha `$ is given implicitly by

$$4\mu ^2\mathrm{cos}\alpha \mathrm{sin}\alpha =m_\pi ^2\mathrm{sin}(\alpha \varphi ).$$
(6)
In the equation above $`\mathrm{tan}\varphi =\lambda /m`$, and $`m_\pi `$ is the mass of the Goldstone modes at chemical potential $`\mu =0`$, temperature $`T=0`$, and diquark source $`\lambda =0`$ . Furthermore the quark-antiquark (chiral) condensate is given by

$$\overline{\chi }\chi =\overline{\chi }\chi _0\mathrm{cos}\alpha .$$
(7)
At nonzero diquark source there is a crossover between the normal phase and the diquark condensation phase. However, Chiral Perturbation Theory at leading order predicts that $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ does not depend on the chemical potential, the diquark source, or the quark mass. This relation does not hold at next-to-leading order in Chiral Perturbation Theory . However, for small chemical potential, quark mass and diquark source, Chiral Perturbation Theory at leading order should give an accurate description of these observables. Therefore, we expect that $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ depends slightly on the quark mass and the diquark source, in the same way as the quark-antiquark condensate depends on the quark mass in three-color QCD at zero temperature and chemical potential . The behavior of $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ is presented in Figure 12. These graphs confirm what was expected from Chiral Perturbation Theory: $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ is close to constant for small enough $`\mu `$ and $`\lambda `$. Notice that it is impossible to see the critical chemical potential, which is around $`0.29`$ as determined in the previous section, by looking at $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$. The points at $`\lambda =0.010`$ show the effects of finite $`dt^2`$ errors. The points which lie higher were generated in simulations with $`dt=0.01`$ while the 4 lower points were generated in simulations with $`dt=0.005`$.

We could compare the diquark condensate obtained on the lattice with the leading-order result (5). But we notice that within Chiral Perturbation Theory at leading order we have that

$`{\displaystyle \frac{\chi ^T\tau _2\chi }{\overline{\chi }\chi }}=\mathrm{tan}\alpha .`$ (8)
This observable can also be used as an order parameter of the diquark condensation phase. In Chiral Perturbation Theory at leading order, for a given quark mass, diquark source and chemical potential, the ratio (8) only depends on the value of the pion mass $`m_\pi `$ through the minimization equation (6). It is therefore a very suitable observable to use in the comparison between the lattice results and Chiral Perturbation Theory.

Hence we fit the lattice result for $`\chi ^T\tau _2\chi /\overline{\chi }\chi `$ with Chiral Perturbation Theory (8,6). This is a one-parameter fit. Since we use Chiral Perturbation Theory at leading order to analyze our results, we can only take into account data corresponding to chemical potentials and diquark sources that are small enough so that $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ is constant for a given diquark source. The previous figure shows that we have to limit our analysis to data corresponding to $`0.28\mu 0.30`$, where $`\overline{\chi }\chi ^2+\chi ^T\tau _2\chi ^2`$ is close to constant for a given diquark source. If we use the 11 points that correspond to $`\lambda =0.004`$, $`0.005`$, and $`0.006`$, and $`0.28\mu 0.30`$, we get that $`\mu _c=0.3027(1)`$. This is an acceptable fit since $`\chi ^2/dof=2.3`$ (remembering that the norm of the condensate is only approximately constant over this range.) This one-parameter fit and the data we used to perform it are shown in Figure 13a. If we now use the 18 points that correspond to $`\lambda =0.0025`$, $`0.004`$, $`0.005`$, $`0.006`$, and $`0.008`$, and $`0.28\mu 0.30`$, we get that $`\mu _c=0.3028(1)`$, and the fit has a $`\chi ^2/dof=3.1`$. If we use the 25 points that correspond to $`\lambda =0.0025`$, $`0.004`$, $`0.005`$, $`0.006`$, $`0.008`$, and $`0.010`$, and $`0.28\mu 0.30`$, we get that $`\mu _c=0.3026(2)`$. The fit is worse than the previous ones with a $`\chi ^2/dof=12`$. This one-parameter fit and the data we used to perform it are shown in Figure 13a. Again we note that the 3 points at $`\lambda =0.010`$ which lie above the curve are those with $`dt=0.01`$ which have larger $`dt^2`$ errors, and are largely responsible for this large $`\chi ^2`$.

Noting that the new estimates for $`\mu _c`$ are significantly higher than our previous estimates, we consider the inclusion of the ‘data’ at $`\mu =0.31`$ and $`\mu =0.32`$ in our fits. If we include the measurements from $`\mu =0.31`$ to the 25 points that correspond to $`\lambda =0.0025`$, $`0.004`$, $`0.005`$, $`0.006`$, $`0.008`$ used in the fits of the previous paragraph, the $`\chi ^2/dof`$ increases from $`3.1`$ to $`4.9`$. Including both the $`\mu =0.31`$ and $`\mu =0.32`$ measurements further increases the $`\chi ^2/dof`$ to $`5.9`$. In both cases the new $`\mu _c`$ is entirely consistent with the estimates of the previous paragraph.

In summary, we find that Chiral Perturbation Theory at leading order describes the data well, with a critical chemical potential given by $`\mu _c=0.3027(3)`$. This is somewhat higher than our earlier estimates in the previous section, and might explain why our $`\delta `$ estimates were less than compelling. It is also understandable in terms of what was learned from the quenched studies performed in . In Chiral Perturbation Theory, the critical exponents are given by mean-field theory even at next-to-leading order. From the quality of the fits of the lattice results with Chiral Perturbation Theory that are presented above, we conclude that the critical exponents measured on the lattice are consistent with mean field.

We notice that Chiral Perturbation Theory at leading order does not describe the data very well at large chemical potential and diquark source. This is presumably a sign that higher order corrections have to be taken into account. The introduction of $`\mu `$ breaks the original symmetry and allows the condensates to vary independently. This can, in part, be implemented by replacing the effective potential based on non-linear sigma models, by a phenomenological effective potential suggested by linear sigma models, where the condensates are less tightly constrained .

## V Scaling function

The static scaling hypothesis states that the data for $`\chi ^T\tau _2\chi /(\mu \mu _c)^{\beta _{mag}}`$ for different diquark sources and chemical potentials should collapse onto a single curve for $`\mu >\mu _c`$ when plotted against $`\lambda /(\mu \mu _c)^\mathrm{\Delta }`$ . The gap exponent $`\mathrm{\Delta }=\beta _{mag}\delta `$; $`\mathrm{\Delta }=3/2`$ in mean-field theory. We can fit the lattice results to the scaling function if we assume some form for the scaling function. We will try several possibilities for the form of the scaling function: linear, of the form $`f(x)=a+bx^s`$, of the form $`f(x)=a+bx+cx^{2/3}`$ and the simplest mean-field form. The number of parameters changes according to the assumption for the form of the scaling function. We present our results in Table 4 using part of the data ($`\mu =0.297`$, $`0.30`$, $`0.31`$, $`0.32`$, $`0.33`$, and $`0.34`$).

The different forms assumed for the scaling function give similar results for the critical chemical potential and for the two critical exponents that we studied. The best fit corresponds to a critical chemical potential that is consistent with the result we found in the previous section, $`\mu _c=0.2948(7)`$, a value of $`\beta _{mag}=0.58(4)`$ which is consistent with mean field, and a value of $`\delta =2.28(3)`$ which is about three quarters of the mean-field result. None of these fits are very good: the best fit has a $`\chi ^2/dof=17`$. It corresponds to the scaling function of the form $`f(x)=a+bx+cx^{2/3}`$, as presented in Table VI. In Figure 14, we show the scaling function together with the $`95\%`$-confidence ellipse for the parameters $`\beta _{mag}`$ and $`\delta `$ using the best fit.

The quality of these fits is rather poor. Furthermore, we find that there is a sizeable correlation between the different fit parameters: $`\mu _c`$ is strongly anti-correlated with $`\beta _{mag}`$ and mildly with $`\delta `$, whereas $`\beta _{mag}`$ and $`\delta `$ are mildly correlated with each other. Therefore, the results from the fits using the scaling function are somewhat suspect. This can be readily seen by using the same data as above plus the four points that correspond to $`\mu =0.2920`$. The best fit we get has a $`\chi ^2/dof=66`$. It gives $`\mu _c=0.269(8)`$, $`\beta _{mag}=0.7(1)`$, and $`\delta =1.4(2)`$. Notice that $`\mu _c`$ is much lower than before, that the increase in $`\beta _{mag}`$ is sizeable, and that $`\delta `$ has changed by almost $`40\%`$. This result illustrates the problems usually encountered in the use of the scaling function. Since the scaling function is in general different above and below the critical point, the choice of data is crucial, and it is difficult to get a clear result from this type of analysis. The comparison with Chiral Perturbation Theory presented in the previous section does not suffer from these problems.

Finally if we impose mean-field exponents, using the same part of the data as above ($`\mu =0.297`$, $`0.30`$, $`0.31`$, $`0.32`$, $`0.33`$, and $`0.34`$), the best fit we get has poor quality: $`\chi ^2/dof=243`$. It gives $`\mu _c=0.297(8)`$, which is compatible with the results obtained in the previous section where we compared the data with Chiral Perturbation Theory. If we try to use scaling functions derived from the mean-field equation of state as in , the quality of the fit is marginally better, with a $`\chi ^2/dof=205`$, and it gives $`\mu _c=0.308(2)`$.

## VI Runs at finite temperature

We now present our $`12^3\times 6`$ simulations to investigate the interior of the phase diagram of Fig. 1 . We are particularly interested in scanning the diagram starting at $`\mu =0`$ and following phase transitions and lines of crossover phenomena at high $`T`$ out to large $`\mu `$. In particular, the low mass four flavor theory should have a first order transition on the $`\mu =0`$ axis and this transition should penetrate into the phase diagram. In addition, the $`T=0`$ transition to a diquark condensate at $`\mu _c`$ should penetrate into the phase diagram and separate the normal phase from the phase with a diquark condensate as shown in Fig. 1. It would be interesting to understand if and how these two lines of transitions merge inside the phase diagram. Finally, we know from past simulations and Effective Lagrangian analyses that there is a line of first order transitions at high $`\mu `$ and high $`T`$ which separates the diquark condensate phase from the quark-gluon plasma phase . We want to confirm the first order character of the high $`\mu `$ part of this line and to find the low $`\mu `$ end of this line where it changes to second order at a tricritical point, labeled 2 in Fig. 1. This point has also been found within Chiral Perturbation at nonzero temperature and chemical potential .

We are hopeful that some of these features of the phase diagram will also occur in other models, so it would be useful to confirm them in this relatively simple setting. Perhaps, QCD at nonzero Baryon Number chemical potential has a $`T`$-$`\mu `$ phase diagram with some of these features.

We begin by considering the $`\mu =0`$ axis. Our simulations are at nonzero fermion mass $`m=0.05`$ and non-zero diquark source $`\lambda =0.005`$ (equivalent to running at $`m=\sqrt{0.05^2+0.005^2}`$ ). $`m=0.05`$ is large enough that we expect the transition between hadronic matter and the quark-gluon plasma to be weakened and perhaps softened into a crossover. This is, in fact, what we find. (At lower quark mass we would expect to find a first order transition.) In Fig. 15 we show data for the chiral condensate and the Wilson Line and confirm a smooth but rapid crossover in the vicinity of $`\beta =1.9`$.

Measurements at $`\mu =0.10`$ gave similar conclusions: there is a clear but smooth crossover in the vicinity of $`\beta =1.9`$, as shown in figure 16.

It is interesting that a plot of the same quantities at $`\mu =0.20`$ displays a noticeably sharper crossover. Figure 17 suggests that $`\mu =0.20`$ is near a critical point, such as the point 1 in our generic phase diagram, Fig. 1. At this rather large quark mass $`m=0.05`$, this appears merely to be a more pronounced crossover but, perhaps, if $`m`$ were chosen smaller than $`0.05`$, an actual line of first order transitions would be found for the region corresponding to $`\mu >0.20`$ in this phase diagram.

Once we increase $`\mu `$ to $`0.25`$, the induced diquark condensate becomes more significant at low $`T`$, as shown in figure 18.

At $`\mu =0.30`$ the diquark condensate becomes the most interesting and rapidly varying quantity. The chiral condensate is relatively smooth and changes from $`0.493(6)`$ at $`\beta =1.7`$ to $`0.217(1)`$ at $`\beta =1.9`$ while the diquark condensate has varied from $`0.352(5)`$ to $`0.0400(3)`$ over the same range of coupling. The Wilson Line varies from $`0.057(3)`$ to $`0.241(2)`$ over the same range. As shown in the figure 19, the diquark condensate experiences a sharp transition between $`\beta =1.80`$ and $`1.90`$.

That transition becomes even sharper as we increase $`\mu `$ to $`0.35`$ as shown in figure 20.

At $`\mu =0.40`$, the simulations indicate a first order transition near $`\beta =1.85(3)`$. In figure 21 we show both the diquark condensates and the Wilson Line and see strong suggestions of discontinuities.

The really solid evidence for a first order transition comes from the time evolution of the observables at $`\beta =1.87`$ which show signs of metastability. For example, in the figure 22 we show the time evolution of the diquark condensate and display several tunnelings between a state having a condensate near $`0.15`$ and another with a condensate near $`0.40`$.

Runs at $`\mu =0.50`$ and $`0.60`$ show clear discontinuities in both the diquark condensate and the Wilson Line near $`\beta =1.87`$. In figure 23 we show the results for $`\mu =0.60`$.

Finally, we scanned the phase diagram at low $`T`$ and variable $`\mu `$. We accumulated data at $`\beta =1.30`$, relatively strong coupling and near vanishing $`T`$, on the $`12^3\times 6`$ lattice and took measurements over a wide range of $`\mu `$, from $`0.15`$ to $`0.60`$. The results for the diquark condensate are shown in the figure 24 which shows a transition near $`\mu =0.25`$. The power law fit shown in the figure produces a critical index $`\beta _{mag}=0.29(4)`$ at a critical chemical potential $`\mu _c=0.2453(2)`$. The quality of the fit, which extends from $`\mu =0.26`$ to $`0.38`$, is rather poor, having a confidence level of only a fraction of one percent. As we have seen in quenched simulations and in simulations of QCD at finite chemical potential for isospin such small estimates of the critical index are often an indication that the scaling is best described in terms of the scaling form for the effective Lagrangian in the non-linear sigma-model class. .

The estimate for the critical chemical potential, $`\mu _c=0.2453(2)`$, meshes well with our measurements at high $`T`$ above, which indicated that the diquark condensate only started to be numerically significant for $`\mu 0.25`$.

The continuous phase transition is also seen in the induced fermion-number density shown in the figure 25. The dotted line fit to that data, which extends all the way from $`\mu =0.26`$ to $`0.56`$, has a confidence level of $`42`$ percent, a critical index $`0.99(2)`$ and an estimate for the critical chemical potential, $`\mu _c=.2282(2)`$. The scaling law obtained here is in good agreement with field theory expectations of unity. Note that such linearity for these relatively strong couplings has a wider range of validity than would be indicated by the scaling window, as is seen in the explicit scaling forms from effective Lagrangians . Since we only have a single $`\lambda `$ value, we have not attempted a fit to such a form.

## VII Conclusion and Outlook

In this article we have presented a study of the phase diagram of 2-color QCD at nonzero baryon/quark-number chemical potential and nonzero temperature. We have found two phases: a “normal” phase where quark number is conserved and a phase with a diquark condensate which spontaneously breaks quark number. At $`\mu =0`$, for smaller quark masses, there should be a first order finite temperature “deconfinement” transition which divides the hadronic from the quark-gluon plasma phase. Since, even in this case we expect the line of first order transitions emerging from this point to terminate at small $`\mu `$, beyond which there should be a rapid crossover but no transition, the hadronic and quark-gluon plasma phases are not distinct. This is the phase we call the “normal” phase. The chiral condensate is non-zero everywhere. However, where it is large the system is more hadronic in nature; where it is small, the system is more plasma-like. In the “diquark” phase, the chiral condensate decreases with increasing diquark condensate, approaching zero for large $`\mu `$ (even before saturation sets in). We have found that the phase transition between the “normal” phase and the diquark phase for small $`T`$ is second order and compatible with mean field theory. These simulations therefore agree with the predictions of Chiral Perturbation Theory through next-to-leading order .

We have also found that the second order phase transition between the “normal” phase and the diquark phase becomes first order when $`\mu `$ and $`T`$ are increased. Therefore the $`\mu `$-$`T`$ phase diagram contains a tricritical point. It corresponds to a chemical potential near two-thirds of the pion mass, and a temperature only slightly lower than the transition temperature between the hadronic phase and quark-gluon plasma phase at $`\mu =0`$. The tricritical point is also found in Chiral Perturbation Theory at nonzero temperature and chemical potential . Therefore, this phase diagram is consistent with the one obtained within Chiral Perturbation Theory .

Now that the phase diagram of this theory has been mapped out, we can turn to more quantitative properties. In particular, we plan on analyzing the gauge configurations for their topological content. There are interesting predictions  based on Effective Lagrangians for the size, density and interactions among the instantons of the model at large $`\mu `$ which can be investigated by cooling methods which have been very successful at vanishing $`\mu `$. In addition, we can simulate the phase diagram at low $`\mu `$ and high T with lighter quarks and search for the critical point 1. The SU(3) analog of this point is thought to be accessible to heavy ion collisions planned for BNL’s RHIC.

Furthermore, the spectroscopy of the model, with an emphasis on Goldstone and pseudo-Goldstone bosons, is also under consideration. Since the theory’s light modes control its critical behavior and thermodynamics, quantitative results on spectroscopy are quite important.

There are also interesting effective Lagrangian predictions for QCD with chemical potentials associated with the light quarks . We are studying QCD with a chemical potential associated with isospin . This model’s phases at nonzero $`\mu `$ and temperature $`T`$ are expected to be very similar to those discussed here . Although we cannot attack the $`SU(3)`$ theory with a large Baryon number chemical potential, these other situations can be studied both analytically and numerically, and interesting new phases of matter have been found there which should be investigated further.

## Acknowledgment

This work was partially supported by NSF under grant NSF-PHY-0102409 and by the U.S. Department of Energy under contract W-31-109-ENG-38. D.T. is supported in part by “Holderbank”-Stiftung. The simulations were done at NPACI and NERSC. B. Klein, K. Splittorff, M.A. Stephanov and J. Verbaarschot are acknowledged for useful discussions.