Partial Deconfinement at Strong Coupling on the Lattice

We provide evidence for partial deconfinement -- the deconfinement of a SU($M$) subgroup of the SU($N$) gauge group -- by using lattice Monte Carlo simulations. We take matrix models as concrete examples. By appropriately fixing the gauge, we observe that the $M\times M$ submatrices deconfine. This gives direct evidence for partial deconfinement at strong coupling. We discuss the applications to QCD and holography.


Introduction
Partial deconfinement [1][2][3][4][5][6] has been firmly established for weakly coupled theories based on analytic methods. In this paper, we focus on strongly coupled theories and use numerical methods based on Lattice Monte Carlo simulations to investigate signals of partial deconfinement. As a concrete setup, we consider the gauged bosonic matrix model with SU(N ) gauge group. This theory's action is the dimensional reduction of the (d + 1)-dimensional Yang-Mills action to (0 + 1)-dimensions. With the Euclidean signature that will be used in Lattice Monte Carlo simulations, the action is given by (1.1) Here I, J run through 1, 2, · · · , d (in this paper we focus on d = 9), β is the circumference of the temporal circle which is related to temperature T by β = T −1 , and X I 's are N × N hermitian matrices. The covariant derivative D t is defined by where A t is the gauge field. This model is often simply called the Yang-Mills matrix model, or bosonic BFSS, because it is the bosonic part of the Banks-Fischler-Shenker-Susskind matrix model [7,8] if d = 9. In the large-N limit, this model exhibits a confinement/deconfinement transition characterized by the increase of the entropy from O(N 0 ) to O(N 2 ) [9]. Concerning this finite temperature transition, partial deconfinement is the phenomenon where only an SU(M ) subgroup of the SU(N ) gauge group deconfines, as pictorially shown in Fig. 1. To characterize partial deconfinement, it is convenient to define a continuous parameter identified by M N whose value can change from 0 to 1. Partial deconfinement happens between the completely confined phase, with M N = 0, and the completely deconfined phase, with M N = 1. These two phases have entropy and energy of order N 0 and N 2 (up to the zero-point energy), respectively. Now, suppose the energy is of order N 2 , where is an order N 0 number much smaller than 1. This is an intermediate 'state '. In fact, the system cannot be in the confined phase because the energy is much larger than N 0 , but it cannot be in the deconfined phase either because the energy is much smaller than N 2 . In partial deconfinement, where an SU(M ) subgroup deconfines, the energy and entropy are of order M 2 , and hence by taking M ∼ √ N such intermediate values of the energy and entropy can be explained. Other explanations that may seem natural are discussed in Sec. 2.
The numerical study in Ref. [10] appears to be consistent with the existence of this partially-deconfined phase. (See also Refs. [11,12] regarding the phase diagram, and Refs. [13,14] for pioneering numerical studies with limited numerical resources which clarified the qualitative nature of the transition.) 1 In the d = 9 Yang-Mills matrix model, the partially-deconfined phase has negative specific heat. Hence, in the canonical ensemble, it is the maximum of the free energy and not preferred thermodynamically. 2 Still, this phase is stable in the microcanonical ensemble. From the point of view of the gauge/gravity duality, this phase is interpreted [1,3] as the dual of a small black hole [16][17][18].
In order to investigate partial deconfinement directly in the nonperturbative regime, we rely on Lattice Monte Carlo simulations. These simulations utilize the path integral formalism of quantum mechanics, where the sum over paths is replaced with a sum over Figure 1: Pictorial representation of partial deconfinement for gauge and adjoint matters degrees of freedom. Only the M × M -block shown in red is excited. This picture is taken from Ref. [4]. In this figure, one specific embedding of SU(M ) into SU(N ) is shown.
'important' field configurations. On the other hand, partial deconfinement is simpler to analyze in the Hamiltonian formalism, with access to the characteristics of individual states in the Hilbert space. Since the field configurations in the path integral are different from the wave functions describing the states, the signals of partial deconfinement become more intricate. Therefore, we devise the following strategy. In Sec. 3, we consider the analytically solvable case of the gauged Gaussian Matrix Model where I = 1, 2, · · · , d. With guidance from analytical results, we derive a few nontrivial features of partial deconfinement in terms of master field. We show that the features of the master field can be seen in the lattice configurations and exemplify them with numerical evidence. Then, in Sec. 4 we move on to investigate the nonperturbative Yang-Mills matrix model, which is the original target of our study. The goal is to determine wether the features of the master field we discovered in the Gaussian matrix model can be applied to the field configurations in our target theory. If that is the case, we can demonstrate that the partially-deconfined phase is an intermediate phase in the confinement/deconfinement transition even at strong coupling. In Sec. 5 we conclude the paper with some discussions regarding the future applications.

Review of partial deconfinement
Partial deconfinement is a phenomenon characterized by deconfinement of an SU(M ) subgroup within the SU(N ) gauge group of a large-N gauge theory. It has been proposed in order to solve a few puzzles associated with the confinement/deconfinement phase transition of gauge theories in light of the gauge/gravity duality. In particular, the original motivation [1] was to explain how the gauge/gravity duality relates the thermodynamics of gauge theories to the physics of black holes: this a very nontrivial problem studied by several papers in the literature, as we discuss below. The clearest example of the connection between thermal phase transitions in gauge theories black holes can be seen in the duality between the thermodynamics of 4d N = 4 super Yang-Mills (SYM) and type IIB superstring theory on AdS 5 ×S 5 [9]. In this duality, the confined and deconfined phases on the gauge side are dual to the thermal AdS geometry and 'large' black hole in AdS space, respectively. In the canonical ensemble, these two phases are separated by a first order phase transition. This transition is called the Hawking-Page transition.
It was immediately realized that between two phases there must be an intermediate state. From the gravity side of the duality, this is simply a very small black hole, which is approximately the ten-dimensional Schwarzschild black hole [16,17]. The energy of this small black hole scales as E ∼ N 2 T −7 , where N 2 corresponds to the inverse of the Newton constant, and this is a stable physical state in the microcanonical ensemble. From the gravity point of view, Schwarzschild black holes with negative specific heat are more realistic than charged black holes such as the 'large' black holes in AdS, which have positive specific heat. Therefore, it is important to understand this intermediate phase.
One of the "puzzles" is how to interpret this small black hole on the gauge side of the duality. In fact, how can a healthy quantum field theory lead to a stable state with a negative specific heat? Partial deconfinement gives a natural answer to this problem [1], introducing a new phase with negative specific heat in thermal phase transitions. Similar phases with negative specific heat are predicted for other theories too. For example, the D0-brane quantum mechanics [7,8] is expected to describe the Schwarzschild black hole in eleven dimensions at very low temperature [19]: such phase would also be understood as a partially-deconfined phase.
While partial deconfinement seems like an original and natural explanation for this problem, it was conjectured to be a general mechanism at work in many theories [2,3]. For example, it has been analytically proven in several weakly coupled theories [4][5][6], by relying on seminal papers [17,20]. These pioneering papers pointed out that the confinement/deconfinement transition characterized by the jump of the energy and entropy can exist even in the weak-coupling limit -it can be kinematical, rather than dynamicaland the intermediate phase resembling the small black hole can exist in general. Even large-N QCD could have such an intermediate phase [21] and it appears to be similar to the cross-over region of real QCD with N = 3. Partial deconfinement gives the precise physical interpretation to this phase [5].
One way to approach partial deconfinement is to analyze the thermodynamics of large-N gauge theories from the point of view of the microcanonical ensemble [3][4][5]. In the microcanonical ensemble, the energy E is varied as a parameter, and the entropy S is maximized at each fixed energy. In the confining phase, E and S are of order N 0 , up to the zero-point energy, while in the deconfining phase they are of order N 2 . This is simply due to the counting of degrees of freedom. In QCD language, we would refer to them as hadrons/glueballs and quarks/gluons. Now we consider a specific value of the energy E = N 2 , where is small but order N 0 number, such as 0.1 or 10 −100 . This is a perfectly reasonable choice because the energy can be varied continuously. On the other hand, the question arises: what kind of phase is realizing this specific energy? It cannot be the confined phase, because the energy is too large and it cannot be the deconfined phase, because the energy is too small. The answer is the partially-deconfined phase with M ∼ √ N . In many theories, including QCD, the canonical and microcanonical ensemble give the same result. The story becomes slightly intricate for systems exhibiting a first-order phase transition in the canonical ensemble, as we will explain in Sec. 2.1. It is worth noting that unless the volume of ordinary space is sufficiently large, the microcanonical ensemble is a physically more realistic setup than the canonical ensemble. This is because the canonical ensemble is typically derived from the microcanonical ensemble as follows. Firstly, let us consider an isolated system in which the energy is conserved, The microcanonical ensemble gives a reasonable statistical description of such system. If the space is sufficiently large, the system can be divided to a small sub-system and large heat bath in thermal equilibrium. Then the small sub-system is described by the canonical ensemble, with the temperature set by the heat bath. By construction, this derivation of the canonical ensemble assumes sufficiently large spatial volume, and hence, is not applicable at a small volume.
Although partial deconfinement has been discovered only recently, from the discussion above it does not look as anything exotic. Another way to approach it is to recall the firstorder transition in a locally interacting system -i.e. the transition between the liquid and solid phases of water -and generalizing it to a system with nonlocal interactions. In the canonical ensemble, water exhibits a first-order phase transition at the temperature of 0 • C and pressure of 1 atmosphere. In the microcanonical ensemble, depending on the energy E the amount of liquid and solid phases change, because of the latent heat at the transition temperature. When E is small/large the completely solid/liquid phase is observed, while in the intermediate range the mixture of two phases appears. Equivalently, the mixture of two phases is realized when the energy is not sufficiently small to be in the completely solid phase, and not sufficiently large to be in the completely liquid phase. The temperature remains fixed because of the short-range nature of the interactions: as long as the interaction at the interface of two phases can be ignored, the temperature cannot change.
A similar mechanism can be applied to the gauge theory phase transition introduced above. In the space of color degrees of freedom two phases -confined and deconfined -can coexist. However, because the interaction between the color degrees of freedom is all-to-all, the temperature can change in a nontrivial way depending on the details of the theory. Pictorially we can visualize three possible patterns as shown in Fig. 2. The blue, orange, and red lines represent the completely confined phase, partially-deconfined phase (or equivalently, partially-confined phase) and completely deconfined phase. These three phases are the counterparts, in color space, of the solid, mixture and liquid phases, respectively. Let us analyze three patterns individually: • The center panel in Fig. 2 would be the easiest one to understand. The Gaussian matrix model studied in Sec. 3 belongs to this class. The temperature does not change, similarly to the case of the mixture of liquid water and ice.
• The Yang-Mills matrix model discussed in Sec. 4 is similar to the left panel in Fig. 2.
In this case, the partially-deconfined phase has a negative specific heat. In the canonical ensemble such phase is not favored thermodynamically and to emphasize this feature we used a dotted line. Strongly coupled 4d N = 4 Yang-Mills and pure Yang-Mills belong to this class too. Depending on the geometry of the ordinary space, instability can set in even in the microcanonical ensemble (see Sec. 2.1 for the details).
• It is likely that QCD with light quarks corresponds to the right panel in Fig. 2, because the lattice QCD simulation suggests a cross-over-like behavior rather than the first order phase transition [22]. The specific heat is positive. The microcanonical ensemble and the canonical ensemble give the same result in this case. It can also be instructive to consider possible objections to the picture presented by partial deconfinement. For example, it can be objected that the phase separation should only take place in ordinary space and not in the internal space of color degrees of freedom. However, it is true that the confinement/deconfinement transition can take place even in matrix models where the ordinary space does not exist by definition. In these cases the phase separation can happen only in the internal 'space' for these theories. Sec. 2.1 contains more details about this point.
Another objection would be to have all the color degrees of freedom get mildly excited and not separate into distinct groups representing the confined and deconfined phases. In the case of water this possibility is forbidden because of the finite latent heat. In the case of the confinement/deconfinement transition we recall that it can take place even in the weak-coupling limit [17,20]. The simplest example is the gauged Gaussian matrix model where the color degrees of freedom are just quantum harmonic oscillators with quantized excitation levels. They cannot be 'mildly' excited precisely because of quantization: the discreteness of the energy spectrum plays the same role of the latent heat.
A third objection can be raised about the SU(M ) block structure of the deconfined sector of partial deconfinement. We could say that there might be other arrangements of the excited degrees of freedom. However, intuitively, since we are maximizing the entropy at fixed energy, it is natural to expect that the solution of such extremization problem preserves a large symmetry. A more precise argument can be made by using the equivalence between color confinement at large N and Bose-Einstein condensation. We refer the interest reader to Ref. [6] where this relation is discussed in details.

Remarks on negative specific heat
Here we want to summarize some remarks about the partial-deconfined phase in the case where it has negative specific heat. We refer again to the leftmost panel of Fig. 2. If the volume of the ordinary space is large, the phase with a negative specific heat is not stable. Any small perturbation can trigger a decay to the co-existence of a completely confined and a completely deconfined phases. This is not necessarily the case if the volume is small and finite. In the case of matrix models, such instability cannot exist by definition, because there is no ordinary space by definition. Moreover, 4d N = 4 super Yang-Mills on S 3 does not have such instability. The partially-deconfined phase in 4d N = 4 super Yang-Mills on S 3 is dual to the small black hole phase, which does not have such instability. It can also be understood via a simple dimensional analysis as follows. In order for the coexistence of two phases in the ordinary space to appear, the radius of S 3 , which we call r, has to be sufficiently larger than the typical length scale of the system, which is the inverse of the temperature of the partially-deconfined phase, β. However β and r are of the same order, as explained in Ref. [9].
When the specific heat is negative, the partially-deconfined phase sits at the maximum of the free energy in the canonical ensemble [3]. The completely confined and completely deconfined phases are the minima of the free energy. The tunneling between the minima can happen only by going beyond this free energy maximum. However, the differences of the free energy between the minima and maximum is of order N 2 and this tunneling is parametrically suppressed at large N . In this way, the local minima is completely stabilized at large N , even when it is not the global minimum. This is very different from the metastable states in the locally interacting systems, such as the supercooled water: even in the thermodynamic limit (large volume) only a small perturbation can destabilize the supercooled water.

Partial deconfinement: the Gaussian matrix model
Let us consider the gauged Gaussian matrix model. The action is given by Eq. (1.2). This model is analytically solvable [17,20] and partial deconfinement has been introduced for the phase transition in Refs. [2,4]. We start from this solvable case in order to analyze how partial deconfinement manifests itself in the path integral formalism.
In the previous section we already mentioned that the Gaussian matrix model has a confinement/deconfinement phase transition which is of first order without hysteresis in the canonical ensemble (the center panel of Fig. 2). The critical temperature is T = T c = 1 log d . In the canonical ensemble, at T = T c , the energy E and entropy S jump from order N 0 to order N 2 , while the Polyakov loop P jumps from 0 to 1 2 (see Fig. 2 of Ref. [4].) In the discussion that follows we have normalized the Polyakov loop such that P = 1 N TrPe i β 0 dtAt , where P stands for path ordering. We also fix the center symmetry ambiguity in the phase of the Polyakov loop by fixing P = |P | for the remainder of this paper.
As functions of the Polyakov loop, the energy and entropy are expressed as 3 and The distribution of the phases of the Polyakov loop is The free energy F = E − T S does not depend on P : Hence all values of P contribute equally to the canonical partition function. Essentially the same results hold for various weakly-coupled theories [17,20]. So far we have introduced important relations for the Gaussian matrix model quantities and we want to relate them to partial deconfinement, following the findings in Ref. [4]. At T = T c , the size of the deconfined sector M jumps from M = 0 to M = N . As proven analytically in Ref. [4], one of the key relations to understand the partial deconfinement mechanism is This equation connects an observable of the system (P ) to the amount of deconfined degrees of freedom M/N . It follows that the distribution of the phases of the Polyakov loop in Eq. (3.3) can be written as Similarly, for the energy E and the entropy S, we can write Note that these relations are valid only at the phase transition temperature T = T c where partial deconfinement takes place. Later in this paper, when we perform the numerical analysis at large but finite N , we determine M from Eq. In order to understand the nature of the states in the microcanonical ensemble, let us rewrite Eqs. (3.6), (3.7) and (3.8) as follows: This way we can separate clearly two different contributions that get summed together. For each equation, the first term of the sum is interpreted as the contribution of the ground state, while the second term is just the value of each observable for an SU(M ) theory at the GWW-transition point (only M degrees of freedom can be excited). Partial deconfinement shown in Fig. 1 naturally explains this M -dependence. One objection to this idea would be that Fig. 1 does not look gauge invariant. Hence, in order to prove partial deconfinement in a gauge-invariant manner, we first show in Sec. 3.1 the explicit construction of the gauge-singlet states in the Hilbert space from the Hamiltonian formalism. Then, in Sec. 3.2, we consider the path integral formalism, which is used in the lattice Monte Carlo simulation and we discuss the properties of the master field.

The Hamiltonian formalism
Following Ref. [4], we explicitly construct the states governing thermodynamics in the gauged Gaussian matrix model. The Hamiltonian iŝ The creation and annihilation operators are defined asÂ † The ground state is the Fock vacuum |0 which is annihilated by all annihilation operators:Â The physical states have to be gauge singlet, e. g.
LetÂ † I be the truncation ofÂ † I to the SU(M )-part. We can construct the states which are SU(M )-invariant but not SU(N )-invariant as LetÔ be a gauge-invariant operator which is a polynomial of O(N 0 ) matrices. We consider these 'short' operators because they do not change the energy too much and we want to study the properties of the states with energy of order N 2 . 4 Then, 2 SU(M )|Ô|SU(M ) 1 = 0, because to connect |SU(M ) 1 to |SU(M ) 2 it is necessary to act O(N 2 ) creation and annihilation operators. This is essentially a superselection rule: different embeddings of SU(M ) to SU(N ) belong to different super-selection sectors. Whether we use a particular embedding or a superposition of all embeddings, we get the same expectation value forÔ.

The path integral formalism and lattice Monte Carlo
So far we only reviewed how partial deconfinement can be seen in terms of the states in the Hamiltonian formalism, following the results in Ref. [4]. Next, we consider the path integral formalism which is used in the lattice Monte Carlo simulations of this paper.

Ensemble properties and master field
An important remark is that the field configurations in the path integral formalism do not have a simple connection to the quantum states in the Hilbert space, other than the fact that the expectation values of gauge-invariant observables agree. This makes the connection with partial deconfinement a little bit more intricate, in particular when trying to analyze the configurations obtained in lattice Monte Carlo simulations. A typical misunderstanding would be "the lattice configurations are the wave functions describing specific states in the Hilbert space"; the absence of such simple connection would be illuminated by noting that lattice configurations have to be averaged in order to obtain the expectation values, unlike the wave function.
At large N , there is a simplification: statistical fluctuations are suppressed at leading order, and we can expect the master field [23] to appear and dominate the path integral. 5 Note that the master field is not the wave function representing the state in the Hilbert space. Still, we can find characteristic features of the master field describing the partially-deconfined phase, by making a 'mapping' between typical states and the master field configuration (typical lattice configurations). Our strategy is to confirm those features by lattice simulation.
In this paper, we refer to the master field as a lattice configuration in the Euclidean path integral of the theory at large N which gives the correct expectation values for properly normalized quantities such as E/N 2 to the leading order in the expansion with respect to 1/N : Here f can be any properly normalized gauge-invariant quantity following the 't Hooft scaling. We need to understand the features of the master field describing the partially 4 The counterpart of this in the case of a finite-N theory at large volume is to consider only the operators with a compact support, in order to make sense of the boundary conditions. 5 For a review of the master field, the readers can refer to Ref. [24].
deconfined phase in the Gaussian matrix model. Then we can start looking for the master field in nontrivial theories such as the Yang-Mills matrix model. 6 Because lattice Monte Carlo simulations are based on importance sampling, the sampleby-sample fluctuations of the properly normalized quantities (e.g. E/N 2 , which is of order N 0 in the large-N limit) are suppressed as N becomes larger. In the strict large-N limit, configurations appearing in lattice simulations can be identified with the master field. However, in actual simulations, we can only study large but finite N values. To learn about the master field in lattice simulations, we simply study the features of configurations sampled by the Markov Chain Monte Carlo (MCMC) algorithm, at sufficiently large N , identifying them with master fields.
When there are multiple saddle points in the path integral, each saddle point has a corresponding master field. In the case of the Gaussian matrix model at T = T c , any M between 0 and N minimizes the free energy, and hence, we need to treat all values of M separately. There is a master field for each M , and we identify them as the dominant configurations at a given (T, M ) pair. This can be done using numerical lattice simulations.
The master field has an ambiguity due to gauge redundancy. In order to eliminate the redundancy, we perform Monte Carlo simulations in the static diagonal gauge. In this gauge, the gauge symmetry is fixed up to S N permutations. The gauge field takes the form where θ 1 , · · · , θ N are independent of t, and θ i ∈ [−π, +π). By using them, the Polyakov loop P is expressed as P = 1 N N j=1 e iθ j . The Polyakov loop phases can be divided into two groups: 1. M of them (we can take them to be θ 1 , · · · , θ M without loss of generality) distributed following the density 1+cos θ 2π ; 2. N − M of them (θ (N −M ) , · · · , θ N ) distributed following the density 1 2π . In terms of the Hilbert space, this corresponds to the separation to the deconfined block and the confined block pictorially shown in Fig. 1 [3,6]. This subdivision is fixing the residual S N permutation symmetry further to S M ×S N −M , where rearrangements inside the two separate groups are indistinguishable in terms of gauge-invariant properties.
We note that in terms of the Euclidean path integral there is still a small residual symmetry. Namely, the same value of θ can appear in both sectors, 7 and the permutation acting on them 8 does not change the distribution of the Polyakov loop phases. As we will see below, this is a feature rather than a bug. 6 A few comments regarding the master field in the completely confined and completely deconfined phases in four-dimensional pure Yang-Mills theory at N = ∞, in the context of lattice Monte Carlo simulations, can be found in Ref. [25]. 7 Strictly speaking, at finite N , because of the Faddeev-Popov term associated with the gauge fixing , θ's cannot exactly coincide. However in the large-N limit neighboring θ's can come infinitesimally close. 8 When we permute θi and θj, we exchange the i-th and j-th rows and columns of the scalars XI as well.
In the rest of this section, we will discuss a few properties of the master field which are related to partial deconfinement. 9 Having the application to the Yang-Mills matrix model in mind, we will demonstrate that such properties are visible in lattice configurations.

Distribution of X I,ij
As a simple characterization of the master field, let us consider the distribution of √ N X I,jj (t), √ 2N ReX I,jk (t) and √ 2N ImX I,jk (t). These are the diagonal and off-diagonal elements of all the scalar hermitean matrices and represent a standard lattice field configuration. We collectively denote them as a random variable 'x' with distribution ρ (X) (x). At T = T c , we want to identify the contributions to this distribution coming from the confined and deconfined sectors. We denote them by ρ (X) con (x) and ρ (X) dec (x), respectively. In fact, we expect a very specific form, where M is related to P by Eq. The distribution ρ (X) (x) has some interesting properties. The part of the action describing the scalar is From this expression, we can see that the distribution of √ 2N ReX I,jk (t) and √ 2N ImX I,jk (t) depends only on θ j − θ k .
From here on, we focus on T = T c where partial deconfinement takes place, and we consider the SU(M )-partially-deconfined phase. If we take the average over j = M + 1, · · · , N or k = M + 1, · · · , N , then the distribution of θ j − θ k is uniform, just as in the completely confined phase. In other words, the distribution of x is the same as in the completely confined phase. On the other hand, if we take the average over j, k = 1, · · · , M we see the difference from the confined phase, because θ j − θ k is not uniform. In particular, the energy E defined by Eq. (3.7) can be directly related to the second moment of ρ (X) (x) by By referring back to Eq. (3.10), we can see that the variances in the equation above can be computed as and The con can be obtained without fixing the S N permutation symmetry, as long as the static diagonal gauge is used. If we pick up only some specific configurations with a common value of P and make a histogram, we should get Eq. (3.19). Numerically, we can determine ρ (X) con and ρ and We can confirm such separation numerically in lattice Monte Carlo simulations. First, from a numerical simulation at fixed d = 2, N and L at T = T c , we sort through the field configurations with fixed Polyakov loop value P such that it reflects the value of M that we are interested it. Those field configurations become the basis for constructing ρ (X) (x) at a pair of M , M values. By using several pairs of M and M , we can construct ρ (X) con (x) and ρ (X) dec (x) solving Eqs. (3.24)-(3.25). In Fig. 3 we confirm that the extracted distributions are indistinguishable from each other, as expected from the equations above. Moreover, the effects due to having only a finite number of lattice sites are very small: using a lattice with only 4 sites is already enough in this case, as it is shown in Fig. 4, showing a good convergence to the continuum limit. The variances in the confined and deconfined sectors calculated from the histograms of ρ (X) con (x) and ρ (X) dec (x) are shown in Table 1. The values agree with the ones at large N , Eqs. (3.22)-(3.23), within the statistical errors and this indicates that the numerical analysis is robust.   Somewhat interestingly, neither ρ (X) At sufficiently low temperature, deep in the completely confined phase, ρ (X) con (x) (which is equivalent to ρ (X) (x) there) approaches a Gaussian, as shown in Fig. 5.

Correlation between scalars and gauge field
Another relevant quantity is the correlation between the scalars X I and the Polyakov loop phases θ i . We consider the quantity K i defined by Numerically, there is an explicit one-to-one correspondence between θ i and K i , given by the components which are labeled by the same index i. This quantity is related to the energy Eq. (3.7) by E = N N i=1 K i . In the completely confined phase, since the entire configuration is in the ground state, and we can identify this contribution with that in the confined sector at T = T c (the light blue in Fig. 1). Next we compute the contribution from the deconfined sector, K i dec . At the Gross-Witten-Wadia (GWW) point in SU(N ) theory, where all elements are thermally excited (the red in Fig. 1), namely at the point T = T c and M = N ,    In the partially deconfined phase, if θ i belongs to the deconfined sector, we expect This is because N −M and M components in |X ij | 2 behave as in the confined and deconfined phases, respectively. For the same reason, we expect K i con = d 2 if θ i belongs to the confined sector. The average should be (3.30) When we look at the distribution of K i 's, one may expect two peaks corresponding to the confined and deconfined sectors. This naive expectation is wrong. In fact, due to the residual symmetry under the exchange of θ's with the same value in the confined and  Fig. 6 we show the two-dimensional histograms of (θ i , K i ) at each M separately. From the two-dimensional histograms we can see only one peak at each θ (represented by a reddish hue). Furthermore we observe that holds with good accuracy when compared to binned histograms. These histograms are created by taking the average over the samples K i falling into a θ i bin of size ∆θ = 0.02. The fluctuation at each fixed θ i can be understood as the finite-N effect which should be suppressed as N becomes larger. By using the distribution of θ i in the confined and deconfined sectors ( 1 2π and 1+cos θ 2π ), we obtain Eqs. (3.27) and (3.29), and hence, also Eq. (3.30).

Partial deconfinement: the Yang-Mills matrix model
The Yang-Mills matrix model with d = 9 exhibits a first order transition near T = 0.885 [10], as sketched in Fig. 7. Slightly different from the Gaussian matrix model, there is a hysteresis in a very narrow temperature range, which can be read off easily from the two-peak signal in the Polyakov loop distribution. Below, we study the properties of the configurations at fixed temperature. As a concrete example, we study T = 0.885, varying the value of P from 0 to 1 2 , along the green dotted line in Fig. 7. We expect that this fixed  At fixed temperature in the transition region, numerically we can find relations similar to (3.9), (3.10) and (3.11) for the Gaussian matrix model [10]. Namely, ρ (P ) (θ) = 1+2P cos θ 2π holds, 10 and by using the identification 2P = M N , we obtain 11 dtTr[X I , X J ] 2 = (N 2 − M 2 )ε 0 + M 2 ε 1 , ε 0 6.14, ε 1 6.60, These relations can be naturally explained if we assume partial deconfinement. The first relation (4.1) can be interpreted that N − M of the phases is in the confined sector, while the other M are in the deconfined sector. The second relation (4.2) would mean that each degree of freedom in the deconfined sector contributes to the increment of the energy by ε 1 − ε 0 . That it appears to be independent of M would be natural because temperature T is fixed. The third relation can be interpreted in a similar manner: M 2 matrix entries are excited as |X ij | 2 ∼ r 1 , while the rest remain |X ij | 2 ∼ r 0 .
Although (4.1), (4.2) and (4.3) are consistent with partial deconfinement, the separation to the SU(M )-and SU(N − M )-sectors shown in Fig. 1 has not been confirmed explicitly in the previous studies. Unlike the Gaussian model, whether the separation to 'confined' and 'deconfined' sectors can work is highly nontrivial due to the interaction. The explicit confirmation of this separation is the goal of this section. Our strategy is to confirm the properties of the master field compatible with partial deconfinement, analogous to the ones we have seen in Sec. 3.
Note also that we do not find a clear theoretical reason that nontrivial M -dependence is forbidden. Due to the interaction between the confined and deconfined sectors, the average contribution in the confined and deconfined sectors may change depending on M . We will come back to this issue later, in Sec. 4.3. In short, whether the M -dependence exists or not does not affect our argument significantly.

The properties of the ensemble and the master field on lattice
In this section, we use the static diagonal gauge, as we did in Sec. 3. Namely the gauge field is fixed to A t = diag θ 1 β , · · · , θ N β . The relations (4.1), (4.2) and (4.3) suggest the deconfinement of the M ×M -block at the upper-left corner, when θ 1 , · · · , θ M are distributed as 1+cos θ 2π and θ M +1 , · · · , θ N are uniformly distributed. If this is correct, then many of the 10 Note that this specific form of ρ (P ) (θ) is not a requirement, though it is observed in various theories.
Note also that we fixed the center symmetry such that P = |P |. 11 The fits performed in Ref. [10] were slightly different, in that the power was not fixed to 2. The Yang-Mills matrix model [10]. Red, orange and blue lines represent the phases dominant at large N , in the canonical or microcanonical ensemble. The red and blue lines are the completely deconfined and confined phases, which are the minimum of the free energy in the canonical ensemble. The orange dotted line is the partially deconfined phase, which is the maximum of the free energy. In order to emphasize that this phase is not favored in the canonical ensemble, we used the dotted line. (Note however that this phase is stable in the microcanonical ensemble.) In the canonical ensemble, there is a first order transition near T = 0.885. We will study the properties of the configurations at T = 0.885, varying the value of P from 0 to 1 2 , along the green dotted line. Ref. [10] found that E and R change as (4.2) and (4.3).
arguments for the Gaussian matrix model presented in Sec. 3 can be repeated without change.
Before showing the results confirming this expectation, let us remark a technical aspect of the simulation. The size of the deconfined sector M can change from 0 to N . Therefore, in order to estimate the quantities such as ρ (X) (x) or K i at fixed M , we need a very long simulation, so that the samples with that specific value of M appear many times. For the Gaussian matrix model, we took this approach, because the simulation cost was low. For the Yang-Mills matrix model, we take a more efficient approach. The idea is to restrict the value of the Polyakov loop P = 1 N N j=1 e iθ j by adding to the action. 12 It enables us to pick up the configurations at a fixed value of M effectively by choosing p 1 and p 2 appropriately, while leaving the configurations at p 1 < |P | < p 2 untouched.

Distribution of X I,ij
The relation Eq. (4.3) suggests the separation of the distribution of X I,ij 's to ρ (X) dec (x) and ρ (X) con (x) just as in the Gaussian model, by using the same expression Eq. (3.19). Let us confirm that this is the case indeed. Here we are assuming that ρ  Table 2.     The agreement with the values in Table 2 is very good. We will show more numerical results for the distributions of x and argue their properties in Appendix B. Though highly nontrivial, this fact alone does not establish the separation into the SU(M )-and SU(N − M )-sectors; logically, it just means the separation between M 2 and N 2 − M 2 degrees of freedom. Furthermore we have assumed that ρ (X) dec (x) and ρ (X) con (x) are independent of M , although, ideally, we do not want to assume it. We will establish the separation to the SU(M )-and SU(N − M )-sectors more rigorously in Sec. 4.2.

Correlation between scalars and gauge field
The correlation between K i ≡ I,j 1 β dt|X I,ij | 2 and θ i is very similar to the case of the Gaussian matrix model. Firstly, because the fit (4.3) works well, holds with good numerical precision. (In the case of the Gaussian matrix model, we could show it analytically.) From this, if we assume that the contributions from the confined and deconfined sectors are always r 0 and r 1 regardless of the value of M , we obtain when θ i is in the deconfined sector. Because the distribution of the Polyakov loop phases in the confined and deconfined sectors are the same as in the Gaussian matrix model ( 1 2π and 1+cos θ 2π ), we naturally expect essentially the same form as Eq. (3.31): In Fig. 10, we show the correlation between θ i and K i obtained by numerical simulations. The values of r 0 and r 1 obtained by using Eq. (4.9) as the fit ansatz are consistent with the values of the variances in Table 2.

Constrained simulation
In this section, we introduce a practically useful approach which makes the M × M -blockstructure manifest. Again, we adopt the lattice regularization in the static diagonal gauge. We separate θ 1 , · · · , θ N to two groups, and define Polyakov loops separately, as (4.10) Then we add the constraint term to the action, Here, γ is a sufficiently large value. Then |P M | and |P N −M | are constrained to be close to 1 2 and 0. We used γ ∼ 10 5 and δ = 0.002 in our simulations. More details about the simulation method will be explained in Appendix A.
If our scenario regarding partial deconfinement is correct, this constraint should fix the S N permutation symmetry and make the upper-left SU(M ) block deconfined while keeping the rest confined, as in Fig. 1 Let us see N = 64 as an example. The distribution of the phases of the Polyakov loop ρ (P ) (θ) obtained from the constrained simulations at T = 0.885 with 24 lattice sites is plotted in Fig. 11. The agreement with 1 2π 1 + M N cos θ is very good. In Fig. 12, E/N 2 and R/N 2 calculated with and without the constraint are compared. We can see good agreement between them. These observations support the expectation that the constraint term in Eq. (4.11) does not alter the theory.   Table 2. The center symmetry is fixed sample by sample such that P = |P |. (a) The two-dimensional histograms of (θ i , K i ). (b) The averaged K i within the narrow bin ∆θ = 0.02. The magenta lines are fit based on Eq. (4.9), by using the best-fit values r 0 = 2.197 and r 1 = 2.302. The values obtained by the fits are consistent with the values of the variances in Table 2. The error bars in each figure are obtained by Jack-Knife analysis.

Distribution of X I,ij
In Sec. 4.1.1, we determined ρ (X) dec (x) and ρ (X) con (x) by using Eq. (3.19), assuming they are independent of M . Now we can determine those distributions much more easily, without assuming the M -independence: ρ (X) dec (x) can be determined from X I,jk with 1 ≤ j, k ≤ M , and ρ (X) con (x) can be determined from the rest. The results are shown in Fig. 13 and Table 3. The error bars are well under control, and we can see a clear difference between ρ (X) dec (x) and ρ (X) con (x). In Fig. 15, the distributions obtained from the constrained and unconstrained simulations are compared. We can see reasonably good agreement. These observations provide us with an explicit confirmation of the M × M -block structure.
By looking at the values of the variances in Table 3 closely, we can see a weak Mdependence. This M -dependence may or may not survive in the continuum limit. We also studied the off-diagonal blocks in the confined sector separately. The result is shown in the same table. We can see a similar M -dependence. In Sec. 4.3, we will give further discussion regarding this. We will explain that the conclusion does not change, even if such M -dependence actually exists.

Correlation between scalars and gauge field
Next let us study the correlation between θ i and K i . We consider two options, with and without taking into account possible M -dependence of ρ (X) con (x) and ρ (X) dec (x). Both options explain the data rather precisely. Firstly let us ignore a possible M -dependence. We use the ansatz for the θ-dependence of K i given by (4.9) and the values of r 0 and r 1 used in Fig. 10. The results are shown   in Fig. 16, with the magenta lines. Of course, Fig. 16 is essentially the same as Fig. 10, except that the constraint term (4.11) is added in the simulations. However this time we can do more: we can easily separate the confined and deconfined sectors and confirm (4.9) separately in each sector; see Fig. 17. This illuminates the residual symmetry in the master field. Next let us consider possible M -dependence as well. This time we can calculate K i dec and K i con directly as  These values can be obtained from Table 3, as and K i con = 9σ 2 con . dec (x), we can write down a reasonable ansatz: The results are shown in the right panels of Fig. 16 and Fig. 17, with the dotted-blue lines. They explain the data very well, and the difference from (4.9) is very small.

Energy
Let us use X dec and X con to denote the deconfined and confined sectors, respectively. By definition, X = X dec + X con . The energy E = − 3N 4β β 0 dtTr[X I , X J ] 2 consists of the contribution from the purely confined part, 16) and the terms involving X dec , E dec would be interpreted as the contribution from the deconfined part and the interaction between confined and deconfined part. A natural guess is that E con does not know whether the SU(M )-sector is deconfined or not. To check it, we construct the counterpart of X con in the completely confined phase. Namely, as shown in Fig. 18, we replace the SU(M )-sector with the confined configuration, by setting both P M and P N −M to be zero. If we take a generic configuration in the completely confined phase, neither P M nor P N −M is zero, although P is zero. Hence we perform another kind of constrained simulation by adding  We calculate the counterparts of E con and E dec , which we denote by E con and E dec , by using this constraint.
The important observation is that E con and E (0) con are very close, as shown in Fig. 19. Therefore, with good numerical precision, the increment of the energy compared to the completely confined phase comes only from E dec .

Summary of the numerical results
Let us summarize the simulation results and see how they are related to partial deconfinement.
In Sec. 4.1.1, we confirmed the separation of the distribution of X I,ij 's to ρ dt|X I,ij | 2 . We see the same kind of correlation as in the Gaussian matrix model, which naturally fit to the partial-deconfinement scenario. In summary, the numerical results shown in Sec. 4.1 are consistent with partial deconfinement, but they are not yet the rock-solid evidence, due to the lack of the demonstration of the SU(M )×SU(N − M )-structure. This is the reason that we needed the constrained simulation introduced in Sec. 4.2.
In Sec. 4.2.1 we confirmed that the constraint term ∆S does not change the theory, except that it fixes the ordering of θ's such that θ 1 , · · · , θ M (resp. θ M +1 , · · · , θ N ) are distributed as 1+cos θ 2π (resp. 1 2π ), which are the form expected for the deconfined sector (resp. confined sector). This means that, if partial deconfinement is actually taking place, then the specific embedding of SU(M ) shown in Fig. 1 should be realized, although we did not touch the scalar fields. In the following subsections, we provided the evidence supporting this expectation. In Sec. 4.2.2, we studied the distributions of X I,ij at the upper-left M × M sector (the red sector in Fig. 1) and the rest (the blue sector in Fig. 1). The results were consistent with ρ dt|X I,ij | 2 . We looked at the statistical features of i = 1, · · · , M and i = M + 1, · · · , N separately, and confirmed that the results are consistent with the specific embedding of SU(M ) shown in Fig. 1. In Sec. 4.2.4, we separated the energy to two parts: E dec , which involves the SU(M )-sector, and E con , which does not involve the SU(M )-sector. Then we showed that the increment of the energy compared to the ground state comes solely from E dec . Again, this is consistent with the SU(M )-partial-deconfinement with the specific embedding of SU(M ) shown in Fig. 1. Based on these observations in Sec. 4.2, we conclude that partial deconfinement is taking place in the Yang-Mills matrix model.
Let us close this section by discussing a potentially small, additional M -dependence. As we have mentioned before, we do not find a theoretical reason that ρ (X) con and ρ (X) dec have to be completely independent of M . Due to the interaction between the confined and deconfined sectors, they might change depending on M .
Via constrained simulations, it is easier to see the M -dependence, if it exists. Actually, as we have seen in Sec. 4.2.2, ρ (X) con and ρ (X) dec appear to have a small M -dependence. In Table 3, we can see that N = 48, 64 and N = 128 exhibit almost the same dependence on M N , and hence, this M -dependence is unlikely to be a finite-N artifact. Somewhat miraculously, the changes of r 0 and r 1 cancel and TrX 2 I can be fit as (4.3) by using the M -independent values. Note also that the analysis in Sec. 4.2.3 was compatible with this M -dependence. Such intricate M -dependence may be a finite-lattice-size artifact, because the lattice size is large but finite (L = 24).
In constrained simulations, we could confirm the separation between two sectors taking into account a possible M -dependence. We did not assume M -independence, and the difference between the confined and deconfined sector was much larger than a possible M -dependence. Therefore, even in case the small M -dependence observed in Sec. 4.2.2 survives in the continuum limit, it does not invalidate our conclusion.

Conclusion and discussion
In this paper, we introduced numerical evidence for partial deconfinement in the Yang-Mills matrix model at strong coupling. In order to establish the numerical methods, we have studied the Gaussian matrix model as well. We identified a few nontrivial properties of the master field which are consistent with partial deconfinement, and confirmed that those properties are visible in lattice configurations. Because the master field can be unique only up to gauge transformation, we used the static diagonal gauge, which drastically simplified the analysis. We expect that other strongly coupled theories exhibit partial deconfinement in the same manner; as discussed in Refs. [3,5,6], heuristic arguments supporting partial deconfinement assume nothing specific to weak coupling.
In this paper, we considered only one fixed value of temperature. It is important to study various different values of temperature and study the temperature dependence at the maximum of free energy, which describe the states realized in the microcanonical ensemble.
A natural future direction is QCD. Suppose, as usual, N = 3 is not too far from N = ∞. Because the QCD phase transition is not of first order [22], the partially deconfined phase is thermodynamically stable [3]. In the large-N limit, the size of the deconfined sector can be read from the distribution of the Polyakov loop phases [6]. If we (perhaps too naively) adopt the relation of weakly-coupled Yang-Mills on S 3 , 0 < P < 1 2 , to roughly identify the partially deconfined phase, partial deconfinement would persist up to several hundred MeV. (See e.g. Ref. [27] regarding the numerical estimate of the renormalized Polyakov loop.) A reasonable starting point to understand the implication of partial deconfinement to QCD would be to study the response of probe fermions to partial deconfinement in simple models such as the Yang-Mills matrix model.
The master fields should be related to classical geometry, when the theory admits weakly-curved gravity dual. 13 Therefore, the specific properties of the master field for the partial deconfinement discussed in this paper should have some geometric interpretation. Obviously, clarifying it within the framework of AdS/CFT is another promising direction of research. At present, we have not yet succeeded to do so, but just to illustrate what the outcome might be, let us end this outlook with a lot of wild speculations. The natural counterpart of the partially-deconfined phase on the gravity side is the small black hole phase [16][17][18] (the top row of Fig. 20). Actually the original motivation to introduce the partially-deconfined phase [1] was to find the dual of the small black hole phase. Up to the 1/N corrections, all the entropy comes from the deconfined sector (black hole), and it appears to be consistent [1] with the Bekenstein-Hawking entropy [33,34]. Hence it is natural to interpret the deconfined and confined phases as the black hole and its exterior. Note that, in general, the confined and deconfined sectors are interacting with each other. 14 Such interaction, which was discussed in Sec. 4.2 when we considered the energy, 13 Previously, the classical dynamics of the Yang-Mills type matrix model has been studied [28][29][30][31][32] with the expectation that typical configurations in the classical theory capture the aspects of the gravitational geometry.
14 Many examples studied previously were weakly-coupled theories and hence the interaction was not important; see e.g. Ref. [4].
could explain the change of the geometry of the exterior compared to the vacuum, due to the existence of the black hole. According to the BFSS proposal [7], block-diagonal configurations, which are partially Higgsed, describe multi-body state. The same interpretation would make sense by using the multiple partially-deconfined sectors (the middle and bottom rows of Fig. 20). Hawking radiation would be described by ripples on the confined sector, or tiny deconfined blocks (strings which are not so long). Local operators can excite tiny deconfined blocks, which propagate in the bulk (confined sector). It would give a natural generalization of the philosophy of BFSS -everything is embedded in matrices -to gauge/gravity dualityà la Maldacena. Note also that the color degrees of freedom in the confined sector can naturally be entangled, and hence, the scenario that the entanglement is responsible for the emergence of the bulk geometry in holography [35,36] would naturally fit to this point of view [4,37]. It would be fun to imagine that the tensor network representing the bulk geometry [38] is hidden in the space of colors. Partial deconfinement may also be related to other mechanisms of emergent geometry such as the ones in the Eguchi-Kawai model [39,40] or IKKT matrix model [41][42][43][44], because the eigenvalue distribution plays the important roles there as well. Numerically tractable targets useful for quantum gravity are the BMN matrix model [45], and perhaps, the BFSS matrix model [7,8,19]. In BFSS matrix model, only the large black hole phase (type IIA black zero-brane) is expected in the 't Hooft large-N limit (i.e. T and λ = g 2 YM N are fixed) [19]. In the very strongly coupled region (T λ 1/3 N −5/9 ), Mtheory becomes a better description, and the eleven-dimensional Schwarzschild black hole should describe thermodynamics [19]. It would be natural to expect that this parameter region, which has negative specific heat, is partially deconfined. Practically, it is not easy to study such strongly coupled region numerically. The situation is better in the BMN model, which is a deformation of the BFSS matrix model with the flux parameter µ. At finite µ, we expect the first order transition at T ∼ µ [46]. At larger µ, the transition temperature is higher, the lattice size needed to study the transition region is smaller, and hence, simulation near the transition temperature is numerically less demanding (see Refs. [47][48][49] for previous attempts), and it might be possible to study the details of the phase transition in near future.
It would also be interesting if the partially deconfined phase could be seen in classical gauge theories, which serve as starting points for learning about real-time dynamics. Different but analogous situation can be seen in the classical real-time dynamics of two-dimensional Yang-Mills theory with adjoint scalar fields [50]. This theory exhibits the 'non-uniform black string phase' [13] which shares a few essential features with the partially-deconfined phase, including the separation in the space of color degrees of freedom [5].
for example, if the initial condition is taken such that θ 1 < θ 2 < · · · < θ N , the target distribution -1+cos θ 2π for θ 1 , · · · , θ M and 1 2π for θ M +1 , · · · , θ N -cannot be realized. To avoid this problem, we randomly choose 1 ≤ p ≤ M and M + 1 ≤ q ≤ N , exchange p-th and q-th row/column and perform Metropolis test. Between each HMC step, 100 random exchanges are performed. Note that only ∆S matters in this Metropolis test. It is easy to see that this procedure does not violate any condition in the Markov Chain Monte Carlo. Therefore, the correct distribution is obtained. That this procedure works shows the non-uniqueness of the 'gauge fixing', as in the Gaussian matrix model. Namely, when the coefficient of the constraint term γ is very large, the permutation takes place only when θ p θ q , which corresponds to the residual permutation symmetry.

A.2 Gaussian matrix model
We have just replaced the potential term of the Yang-Mills matrix model with In this Appendix, we give a few observations regarding the distribution ρ (X) in the Yang-Mills matrix model, which appear to be very different from the case of the Gaussian matrix model. We compare the distributions ρ (X) con (x) and ρ (X) dec (x) with Gaussian distributions whose variances are given by Eq. (4.5) and Eq. (4.6). The results are shown in Fig. 21. We show only the distributions ρ (X) con (x) and ρ (X) dec (x) computed by a pair of (M, M ) = (16,48) because the M -dependence is small. We can see that ρ (X) con (x) and ρ (X) dec (x) are close to Gaussian distributions, while we can see a small but non-vanishing deviation which may be a finite-N or finite-lattice-spacing effect. Note that, in the Gaussian matrix model, ρ (X) con (x) and ρ (X) dec (x) at T = T c are far from being Gaussian distributions. Next, we compare ρ (X) con in the transition region with ρ (X) (x) obtained at low temperature, which is in the completely confined phase. In the left panel of Fig. 22, we show the results of this comparison. We can see that the temperature-dependence is small. Again, this is different from the Gaussian matrix model; see Fig. 5 regarding a large temperaturedependence in the Gaussian matrix model.