Two-colour QCD phases and the topology at low temperature and high density

We delineate equilibrium phase structure and topological charge distribution of dense two-colour QCD at low temperature by using a lattice simulation with two-flavour Wilson fermions that has a chemical potential μ and a diquark source j incorporated. We systematically measure the diquark condensate, the Polyakov loop, the quark number density and the chiral condensate with improved accuracy and j → 0 extrapolation over earlier publications; the known qualitative features of the low temperature phase diagram, which is composed of the hadronic, Bose-Einstein condensed (BEC) and BCS phases, are reproduced. In addition, we newly find that around the boundary between the hadronic and BEC phases, nonzero quark number density occurs even in the hadronic phase in contrast to the prediction of the chiral perturbation theory (ChPT), while the diquark condensate approaches zero in a manner that is consistent with the ChPT prediction. At the highest μ, which is of order the inverse of the lattice spacing, all the above observables change drastically, which implies a lattice artifact. Finally, at temperature of order 0.44Tc, where Tc is the chiral transition temperature at zero chemical potential, the topological susceptibility is calculated from a gradient-flow method and found to be almost constant for all the values of μ ranging from the hadronic to BCS phase. This is a contrast to the case of 0.88Tc in which the topological susceptibility becomes small as the hadronic phase changes into the quark-gluon plasma phase. ar X iv :1 91 0. 07 87 2v 1 [ he pla t] 1 7 O ct 2 01 9


Introduction and summary
It is of great importance to probe properties of matter in extreme conditions as encountered in neutron stars and relativistic heavy ion collisions. In fact, recent observations of the gravitational waves from a double neutron star merger event [1], as well as heavy-ion collision experiments performed at the highest energies that can be achieved at present [2], provide empirical information about the equation of state (EOS) of matter that occurs in such systems and its evolution. Since the constituents of such matter are strongly coupled, however, one has difficulty in accurately evaluating the EOS and transport coefficients from QCD.
Theoretically, all one knows at relatively low temperature are the equilibrium properties of matter near zero baryon chemical potential and of matter at sufficiently high baryon chemical potential for weak coupling calculations to work. At zero baryon chemical potential, the chiral condensate and the Polyakov loop are known to gradually change with increasing temperature from lattice QCD simulations at physical pion mass [3,4]. The system undergoes a gradual crossover from a hadronic to a quark-gluon plasma (QGP) phase around a temperature T c at which the chiral susceptibility is peaked. The topological susceptibility is also measured on a lattice, which indicates that the density of instantons decreases from the value at zero temperature to a vanishingly small one [5]. At nonzero baryon chemical potential, several works based on the multi-parameter reweighting [6,7] and the Taylor expansion [8] searched for a critical end point, but unfortunately it remains to be determined. Thanks to a sign problem inherent to finite density lattice QCD, however, the regime in which one can predict from finite density QCD is restricted all the way to a regime of such high densities as to make perturbative, weak coupling calculations reliable. In this regime, a deconfined superfluid BCS state is predicted to occur at sufficiently low temperatures for Cooper pairs of diquarks to be formed via one-gluon exchange in the colour antisymmetric channel [9].
In this work, we examine the temperature and density dependence of the phase and topological structure of dense two-colour QCD by using a lattice simulation with twoflavour Wilson fermions 1 . In this case, the sign problem can be avoided as follows: The fundamental representation (=quark) of the SU(2) gauge group takes a (pseudo)real representation, so that the determinant of the Dirac operator even at finite quark chemical potential µ is a positive-real or negative-real [11][12][13]. Then, the SU(2) gauge theory coupled to the even number of flavours does not suffer from the sign problem.
Even so, there is another problem: A numerical instability occurs in the low temperature and high density regime. This problem comes from a dynamical pair-creation and/or pair-annihilation of the lightest hadrons, so that it is theoretically expected that the instability appears when µ m PS /2 is satisfied at temperatures close to zero [11][12][13]. Here, µ and m PS denote the quark chemical potential and the pseudoscalar meson mass at zero chemical potential. One can see how this kind of instability occurs and can be avoided by introducing a diquark source in terms of the eigenvalue distribution of the Dirac matrix [14]. In this work, we thus introduce the diquark source term in the action to solve the second problem following refs. [15][16][17][18][19][20][21][22][23][24] (see section 2) and attempt to find out the phase structure in the vanishing limit of the source term.
As for the phase structure of dense two-colour QCD with N f = 2, earlier lattice simulations [18][19][20][21][22][23][24] including the diquark source parameter j in the action provided low temperature data for the diquark condensate, the Polyakov loop, the quark number density and the chiral condensate as a function of quark chemical potential µ. By taking the limit of j → 0, a hadronic phase with vanishing diquark condensate was found to change into a superfluid phase with nonzero diquark condensate around µ ∼ m PS /2. The phase diagram predicted from the earlier analyses can be schematically illustrated in figure 1. We remark in passing that the two-colour system has baryons as bosons that are degenerate with mesons in vacuum, a feature clearly different from the usual three-colour system [25][26][27].
Consistency of the µ and j dependence of various quantities with the mean-field prediction by the chiral perturbation theory (ChPT) [15,25,26,28] was examined in detail. Since the ChPT is applicable at sufficiently small µ, it can at most predict the appearance of a diquark Bose-Einstein condensed (BEC) phase around µ ∼ m PS /2 and the µ dependence of the diquark and chiral condensates. At still higher densities, it was found that the BCS phase, which is characterized by the presence of the Fermi sphere in the underlying normal phase, takes over from the diquark BEC phase. However, it is still controversial whether or not one could see deconfinement at sufficiently low temperature.
It is well-known that the two-colour system at the finite density has a similar property Figure 1. Schematic phase diagram of dense two-colour QCD as predicted from earlier works. The details of the phases will be given in subsection 2.2.
to three-colour QCD with isospin chemical potential [29]. In the three-colour QCD, it is expected that the first order confinement-deconfinement transition occurs at nonzero temperature even at asymptotically high densities [29][30][31]. The intuitive picture of the first order phase transition is as follows: In the intermediate temperature and high density regime, a typical momentum of quarks is higher than the size of the Fermi surface (µ) or zero-temperature diquark gap (∆ 0 ) and hence the quarks behave like free particles, while at extremely low temperature the typical momentum of quarks is smaller than µ or ∆ 0 and hence the quarks are quenched. As is well-known, the quenched QCD at low temperature exhibits the confinement, so that the low temperature and high density regime must be described by the confining theory. By analogy, there could be a similar phase transition by changing temperature in a high density regime [32]. Since two-colour QCD at finite temperature exhibits the second order confinement-deconfinement phase transition, the transition might be weaker than the three-colour case. An important and nonperturbative quantity related to the phase structure is the topology. At low temperatures and intermediate densities, theoretical works are based mostly on low energy effective theories [33]. Some of the low energy effective theories focus on instantons since such topological objects are closely related to the vacuum structure. In fact, the instanton-induced interaction among quarks can promote the chiral phase transition as well as the superfluid phase transition [34], just like the one-gluon exchange interactions at extremely high densities. The key question here is how the density of instantons depends on the baryon chemical potential. To answer this question, however, nonperturbative calculations are indispensable. According to the BCS theory, the zero-temperature diquark gap (∆ 0 ) is related to the critical temperature (T SF c ) for the transition from the superfluid to QGP phase at given chemical potential as where γ E denotes Euler's constant [34]. If there are some topological objects in the superfluid phase and they increase the diquark gap at zero temperature 2 , therefore, the phase diagram is expected to be changed from the perturbative prediction at intermediate densities.
As for lattice approach to the topological properties of dense two-colour QCD, the earliest work was done for N f = 4 and 8 [17,20], which shows that the topological susceptibility decreases with µ in the high density regime, while the result of N f = 2 is available only on a rather coarse lattice [20], which suggests that it might be almost constant in the high density regime 3 . However, the temperature and phase dependences of the topological susceptibility were not considered. Since it is well-known that the property of the topological distribution depends strongly on the temperature T at µ = 0 [5], it is important to study the topology at different temperatures with N f fixed. In this work, therefore, we will investigate the phase structure at two relatively low temperatures, T = 0.45T c (section 3) and T = 0.89T c (section 4), where T c denotes the chiral transition temperature at µ = 0. After that, we will study the density dependence of the topological susceptibility (see section 5). Figure 2. Summary of the phase diagram and the topological susceptibility (χ Q ) of dense twocolour QCD revealed by this work. The BCS (deconfined) phase has yet to be found in this work, and will be investigated in future work.
We summarize the main conclusions of this paper, which are schematically displayed in figure 2. By systematically measuring the diquark condensate, the Polyakov loop, the quark number density and the chiral condensate with improved accuracy and j → 0 extrapolation over earlier investigations [18,21,24], the known qualitative features of the low temperature phase diagram are reproduced at T = 0.45T c . In improving the j → 0 extrapolation, we have developed a new reweighting method. In addition, we newly find that around the boundary between the hadronic and BEC phases, so-called the hadronic matter regime, nonzero quark number density occurs even in the hadronic phase in contrast to the ChPT prediction, while the diquark condensate approaches zero in a manner that is consistent with the ChPT prediction. At the highest µ, which is of order the inverse of the lattice spacing, all the above observables suffer from a strong lattice artifact and change drastically. There is no clear evidence for deconfinement at densities below the regime where a lattice artifact takes effect; this is consistent with the analysis with staggered fermions [24], but at odds with that with Wilson fermions [21]. Finally, at temperature of 0.45T c , the topological susceptibility is calculated from a gradient-flow method [38] and found to be almost constant for all the values of µ ranging from the hadronic to BCS (confined) phase. This is a contrast to the case of 0.89T c in which the topological susceptibility becomes small as the hadronic phase changes into the QGP phase.

Simulation detail
In this section we start with the lattice action. We then define the physical quantities to be measured and show how to determine the phases from the behaviour of the quantities. We finally set the simulation parameters.

Lattice action
In this work, as a lattice gauge action, we utilize the Iwasaki gauge action, which is composed of the plaquette term with W 1×1 µν and the rectangular term with W 1×2 µν , where β = 4/g 2 0 in the two-colour theory, and g 0 denotes the bare gauge coupling constant. The coefficients c 0 and c 1 are set to c 1 = −0.331 and c 0 = 1 − 8c 1 .
As a lattice fermion action, we use the two-flavour Wilson fermion action including the quark number operator and the diquark source term, which is given by Here, the indices 1, 2 denote the flavour label, and µ is the quark chemical potential. The additional parameters J andJ correspond to the anti-diquark and diquark source parameters, respectively. For simplicity, we put J =J and assume that it takes a real value. Note that J = jκ, where j is a source parameter in the corresponding continuum theory. The factor κ comes from the rescaling of the Wilson fermion on the lattice as will be specified below. The C in the last two terms is the charge conjugation operator, and τ 2 acts on the colour index. Finally, ∆(µ) is the Wilson-Dirac operator including the number operator, which is explicitly given by where κ is the hopping parameter. Now, the action (2.2) contains three types of fermion bilinears:ψψ,ψψ and ψψ. To build a single kernel matrix from the fermion action, we introduce an extended fermion matrix (M) as The square of the extended matrix can be diagonal because J(=J) is real. We thus obtain (2.6) Note that det[M † M] corresponds to the fermion action for the four-flavour theory, since a single M in eq. (2.5) represents the fermion kernel of the two-flavour theory. To reduce the number of fermions, we take the square root of the extended matrix in the action and utilize the Rational Hybrid Montecarlo (RHMC) algorithm in our numerical simulations.

Observables and definition of phases
Now, we turn to the notation of the phases as described in figure 1. Following earlier investigations [24], we use the name of each phase as shown in which plays the role of an approximate order parameter for confinement. qq denotes the diquark condensate, which is nothing but the order parameter for superfluidity. n q denotes the quark number density, which is the time-like component of a conserved current. Note that this quantity does not require a renormalization and goes to 2N f N c in the high µ limit on the lattice. Incidentally, n tree q is the quark number density defined by the free field on a finite lattice [18]: (2.10) wherek In the continuum limit, n tree,cont q = µ 3 /(3π 2 ). The nonzero value of the diquark condensate indicates that the system is in a superfluid state. We expect that a low density regime in the superfluid phase could be described in terms of a diquark BEC picture. In this picture, the diquark (q-q) pairs behave as bosons that form a condensate when a typical correlation length of the system, namely, the coherence length of the pair, is sufficiently shorter than the inter-quark spacing ∼ µ −1 . This is the case with a rarefied system. With increasing density, however, the inter-quark spacing becomes shorter and shorter. If the spacing is much shorter than the typical correlation length, then, the number of the pairs significantly decreases since only quarks near the Fermi surface that would occur in the absence of attraction between quarks can undergo Cooper pairing. At the same time, the Cooper pairs spatially overlap with each other, which gives rise to a coherent state. It is nevertheless noteworthy that most of the quarks behave like free particles. For this phase, generally referred to as the BCS phase, the quark number density can thus be estimated from the free particle propagator including the chemical potential µ. We thus expect that even on a lattice, the quark number density can be approximated by n tree q in the BCS phase. Furthermore, we expect that qq is proportional to µ 2 , which corresponds to the area of the Fermi surface.
It might be also tempting to consider the chiral condensate, qq ≡ ψ i ψ i for one flavour, as an order parameter of chiral symmetry breaking. In our numerical calculations, however, we utilize the massive Wilson fermions, which explicitly break the chirality. Although we will discuss numerical results for the chiral condensate in each phase, we will not use the quantity for the definition of the phase.
We conclude this subsection by mentioning that a lattice artifact can induce a peculiar phase in the high density regime. There are two possible origins that can result in such a phase in the numerical simulations. The first one can be found by observing the behaviour of the diquark condensate. In fact, when the value of µ approaches 1 in lattice unit, the propagation of the quarks is quenched, leading to a significant decrease in the amplitude of qq with increasing µ. This kind of behaviour was reported in refs. [16,24] in the N f = 2, 4 simulations using the staggered fermions.
The second origin is the finite volume effect, which manifests itself when all the lattice site are occupied by the quarks. The signal of a peculiar phase coming from the finite volume effect can be found by observing the saturation of n q in lattice unit, which corresponds to a 3 n q = 8 in our notation. In the actual simulation of this work, the maximal value of n q is 0.707(2), which is significantly smaller than 8. We can thus safely avoid such an artifact.

Simulation parameters
We perform the simulations with (β, κ, N s , N τ ) = (0.800, 0.159, 16,16) and (0.800, 0.159, 32,8). According to preliminary results in ref. [39], (β, κ) = (0.800, 0.159) gives m PS /m V = 0.823 (9) and am PS = 0.623 (3), where m V is the vector meson mass in vacuum. Furthermore, the chiral phase transition at µ = 0 occurs at (β, N τ ) = (0.900, 10) with the simulations performed at fixed m PS /m V . Eventually, the scale setting utilizing w 0 [40] in the gradient flow method [38] leads to T = 0.45T c and 0.89T c , respectively, for the former and latter simulations.   We start with the Polyakov loop, whose µ-dependence is depicted in figure 3. The symbols (colours) denote the different values of the diquark source (j = J/κ) in lattice unit; the data with circle (red), square (blue), diamond (green) and triangle (magenta) are generated by the aj = 0.00, 0.01, 0.02 and 0.04 simulations, respectively. Note that the theoretical threshold µ value of the numerical instability is µ/m PS = 1/2, beyond which the lightest hadrons are frequently created and annihilated. In fact, we can implement the HMC simulation without the diquark source for µ/m PS < 0.50, while the Metropolis test in the HMC simulation is not accepted even within a tiny Montecarlo step (∼ 1/1000) for µ/m PS ≥ 0.50. In the latter regime, we thus have to utilize the RHMC algorithm including the diquark source in the action.
As can be seen from the left panel of figure 3, the Polyakov loop takes on a nonzero value for µ/m PS 1.2. Normally, it is a signal of the transition or crossover from the confined to deconfined phase. The right panel of figure 3 shows the susceptibility of the Polyakov loop obtained for aj = 0.01-0.04. The susceptibility is marginally peaked at µ = µ D with µ D /m PS ∼ 1.44. A clear j-dependence of the peak position is not observed in our calculation. Since aµ D is close to unity (aµ D = 0.90), the obtained structure of the Polyakov loop might be associated with a lattice artifact as we shall see.

Diquark condensate
The other quantity that helps determine the phase diagram is the expectation value of the diquark condensate, which have been measured by utilizing the noise method for each configuration. It is however difficult to evaluate this expectation value from the extrapolation of j → 0 limit, as is the case with earlier investigations [18,21,24]. Let us then propose a reweighting method with respect to j at fixed β, µ, κ 4 . In this method, the reweighting factor from the original lattice parameters (β 0 , κ 0 , µ 0 , j 0 ) to the measured parameters (β, κ, µ, j) is given by In our simulation, the value of j alone is changed between the configuration generation and the measurement of the observables. Thus, the reweighting factor explicitly reads In our calculations, the value of J 2 = (j · κ) 2 is typically O(10 −6 ) in lattice unit, so that the Taylor expansion, det[1 + A] = e Tr ln[1+A] ∼ 1 + TrA, is valid. In fact, for µ/m PS < 0.50 (aµ ≤ 0.30), the data for aj > 0 are obtained by the reweighting from the configurations generated with aj 0 = 0.00. On the other hand, for µ/m PS ≥ 0.50 (aµ ≥ 0.31), the data for aj ≤ 0.01 are measured by using the reweighting method from the configurations with aj 0 = 0.01. It is numerically confirmed that the leading correction of the reweighting factor (R j − 1) is less than 10 −3 in our calculations. We firstly investigate the consistency of the data obtained by the different methods around µ/m PS = 0.50. In figure 4, the square (blue) data for µ/m PS = 0.48 are generated by the HMC algorithm without the j term in the action and are measured by using the reweighting method from the configurations with aj = 0.00, while the triangle (red) data for the same µ/m PS = 0.48 are obtained by the configurations generated by using the RHMC algorithm for each value of j and are also measured using the same value of j. Both data for each j are 3σ consistent, and furthermore both extrapolated values of qq in the j → 0 limit vanish. We also show the diamond (green) data for µ/m PS = 0.50 (aµ = 0.31) in the figure, where the data are generated by the RHMC algorithm with aj = 0.01, since as explained, at this aµ it is very hard to generate the configurations utilizing the HMC without the j term. The diquark condensate for µ/m PS ≥ 0.50 are measured by the reweighting method from aj = 0.01.
The left and right panels of figure 5 present how the diquark condensate depends on j in the low µ and high µ regimes, respectively. In the j → 0 limit, the left panel shows that the data go to zero for µ/m PS < 0.50, as denoted by circle (red) and square (blue) symbols, while having a nonzero expectation value for µ/m PS ≥ 0.50, as denoted by diamond (green) and triangle (magenta) symbols. The critical value, µ/m PS = 0.50, corresponds to the onset of the numerical instability, at which superfluidity is known to occur [14]. At this critical µ, the ChPT predicts that the diquark condensate behaves like j 1/3 , and in the case of the study with staggered fermions, the numerical evidence was obtained [24]. Although the chi-square of the linear extrapolation in aj becomes a little worse only around µ/m PS = 0.50, there is no clear evidence for this behaviour in our simulation.
In the high µ regime, as can be seen from the right panel in figure 5, the diquark  condensate has a larger j dependence for higher µ. Normally, the extrapolated value of the diquark condensate increases with µ, namely qq ∝ µ 2 , in the BCS phase, but it starts to decrease around µ/m PS = 1.28, which corresponds to aµ = 0.80. This behaviour might be a signal of the strong lattice artifact, as discussed in subsection 2.2. It might be notable that the similar behaviour occurs at high aµ in refs. [16,24], though the lattice fermion in both works is the staggered one. Now, let us exhibit, in figure 6, the phase diagram at T = 0.45T c by using the data of the Polyakov loop and the diquark condensate. According to the definition of the phases given in table 1, the hadronic and superfluid phases appear at T = 0.45T c . Although it remains to be determined whether or not the superfluid phase is composed of the BEC and BCS phases, it will be addressed in the next subsection. Thanks to the effectiveness of the reweighting method concerning j, we can precisely find the location of the transition between the hadronic and superfluid phases as µ B /m PS 1/2. This result is consistent with the ChPT prediction that for µ ≥ m PS /2 the lightest hadrons are frequently created.
The scaling law around the critical point, can be investigated. Here, we have taken µ B = m PS /2 and β m = 0.50, which are the meanfield predictions by ChPT [26], since the determination of µ B and β m from the numerical data was very hard. The chi-square over d.o.f. of the fitting for the 9 data points in 1/2 ≤ µ/m PS < 0.70 is reasonable (χ 2 /(d.o.f.) = 1.31), and we find that the theoretical curve based on the scaling law is consistent with the data as shown in the right panel in figure 6.

Quark number density and chiral condensate
We turn to the quark number density and chiral condensate, which have been measured by utilizing the noise method. The result for the quark number density is shown in figure 7. We can observe from the left panel a regime where the lattice data for the quark number density n q are almost consistent with the free field behaviour (n tree q ). This regime, which ranges 0.72 µ/m PS 1.28, may well be identified as the weakly coupled BCS phase, rather than the strongly coupled BEC phase, which ranges 0.50 µ/m PS 0.72. Some works based on the Nambu-Jona-Lasinio model have also found a roughly consistent crossover point µ/m PS = 0.8-0.85 in refs. [42,43].
It is interesting to note that even in the hadronic phase, the quark number density in the j → 0 limit takes on a nonzero value in a narrow regime just below µ/m PS = 1/2, as can be clearly seen from the right panel in figure 7. In fact, such a regime ranges 0.42 µ/m PS 0.50. Here, this regime is referred to as "hadronic matter". The presence of the hadronic matter is consistent with a recent Schwinger-Dyson analysis [44]. This is apparently at odds with the ChPT prediction that n q becomes nonzero when µ reaches m PS /2. Note, however, that the ChPT prediction holds for sufficiently low temperature and that the lattice data have been obtained at T 0.1m PS . Since the ChPT predicts that the diquark (antidiquark) mass in the medium is m PS − 2µ (m PS + 2µ), the fact that T m PS − 2µ holds at µ 0.42m PS implies that diquarks are thermally excited. It is thus reasonable for the quark number density to start increasing from nearly zero at µ 0.42m PS . This implication could be confirmed by measuring the diquark mass in the medium as in refs. [12,19,45], which will be addressed in the future. Finally, for completeness, we show the chiral condensate in the left panel of figure 8. Interestingly, a plateau appears in the BCS phase (0.72 µ/m PS 1.28). According to the ChPT prediction up to leading order, the sum of the squared diquark condensate and the squared chiral condensate has to be constant in the finite density regime. As shown in the right panel of figure 8, however, this relation seems violated. Even in the BCS phase, the plateau of the chiral condensate and the increment of the diquark condensate with µ (see figure 6) suggest a slight violation of the relation. More careful analyses will be needed to clarify such a violation.
4 Simulation results: Phase diagram at T = 0.89T c Next, we investigate the phase diagram at T = 0.89T c . This temperature is significantly higher than T = 0.45T c as depicted in section 3 but still below T c , the chiral transition temperature at zero chemical potential. It is noteworthy that we can carry out the HMC simulation without the diquark source term even in the high chemical potential regime, µ/m PS ≤ 1.61 (aµ ≤ 1.00). This suggests that at T = 0.89T c , superfluidity does not appear at any density.
To clarify the absence of superfluidity, we have performed the RHMC simulation with the diquark source term (aj = 0.001, 0.003, 0.005, and 0.01) at aµ = 0.30 and 0.70. In the left panel of figure 9, there are two types of data for each aµ: The one is obtained by the reweighting method for nonzero aj where the configurations are generated by the HMC algorithm without the j term. The other one is generated by the RHMC algorithm for each value of aj. These data are consistent within the statistical error bars. In the right panel of figure 9, we show the data obtained in the former way for several values of µ satisfying aµ ≤ 1.0. By carrying out a linear extrapolation in terms of aj, we confirm that the diquark condensates at j = 0 are consistent with zero. Thus, the superfluidity does not appear at this temperature. By using the configurations obtained from the HMC simulation without the diquark source term in the lattice action, we have also measured the Polyakov loop, the chiral condensate, and the quark number density as shown in figure 10. The result for the Polyakov loop suggests that as µ increases, the system tends to be deconfined, while simultaneously the chiral symmetry tends to be restored. Furthermore, the quark number density has a small but nonzero value already around µ/m PS = 0.16, which is significantly smaller than the critical value of µ B /m PS = 1/2, which is given by the ChPT and the lattice study at T = 0.45T c . Note that as long as µ > 0, nonzero n q can occur via thermal excitations. We can thus conclude that at T = 0.89T c , there is only a transition or crossover from the hadronic to the QGP phase at nonzero µ that is closer to µ/m PS = 0.16 rather than µ/m PS = 1/2, while the superfluid phase does not appear even in the high density regime below T c .

Simulation results: The topology at finite density
We are now in a position to examine the topological charge (or instanton number) in each of the phases as elucidated in sections 3 and 4. The topological charge can be measured by incorporating the gluonic definition, into the gradient flow method. The number of configurations measured for each set of the lattice parameters has been set to 50-100, while the flow time of the gradient flow has been fixed at t/a 2 = 40. Figure 11 depicts the histogram of the topological charge obtained at T = 0.45T c for µ/m PS = 0.00, 0.56 and 1.12 (aµ = 0.00, 0.35 and 0.70), respectively. Note that these three cases are typical of the hadronic, BEC and BCS phases, respectively. Remarkably, no clear  difference can be seen among them. This result could have consequence to the diquark gap because the instanton density is expected to control the strength of the attraction between quarks [34]. This result is qualitatively different from what earlier investigations found for the SU(2) N f = 4 theory on 12 3 × 24 lattice [20] and SU(2) N f = 8 theory on 14 3 × 6 lattice [17]. In such investigations, the topological susceptibility was found to decrease in the high chemical potential regime. To understand such a qualitative difference, we also investigate the temperature and the phase dependences of the topological susceptibility.  Figure 12 shows the histogram of Q obtained at T = 0.89T c , where the phase changes from the hadronic to QGP one with increasing µ. At zero µ, there is still a broad distribution of Q, while in the high chemical potential regime (aµ = 0.90), most configurations localize at the Q = 0 sector. Consequently, the topological susceptibility decreases with µ.
We conclude this section by showing, in figure 13, how the topological susceptibility and the Polyakov loop depend on µ at T = 0.45T c and 0.89T c . Here, to make it easy to see, χ Q multiplied by a constant 25 and 200 are shown in the left and right panels, respectively. As can be already seen from figures 11 and 12, there is a crucial difference between these two cases: At T = 0.89T c the topological susceptibility decreases with µ, which seems just opposite to the behaviour of the Polyakov loop, whereas at T = 0.45T c the topological susceptibility is almost independent of µ even in the regime where the Polyakov loop significantly increases. The qualitative difference comes from the underlying phase structure.

Outlook
In this paper we have laid out the phase and topological structure of dense two-colour QCD on lattice. In addition to accuracy improvement, however, many questions remain. First, scale setting [39] has to be finalized for a continuum extrapolation to remove lattice artifacts. Second, it would be significant to find the critical temperature at which the diquark condensate disappears in high density regime. Especially in the BCS phase, such critical temperature would lead to estimates of the zero-temperature diquark gap given in eq. (1.1). Finally, it would be interesting to extend the present work to the cases of lower temperature and larger lattice. These lattice studies would help us study the hadron masses [12,19,45] and the hadron-hadron interactions [46][47][48][49][50].