Phase transitions in 5D super Yang-Mills theory

In this paper we study a phase structure of $5D$ ${\cal N}=1$ super Yang-Mills theory with massive matter multiplets and $SU(N)$ gauge group. In particular, we are interested in two cases: theory with $N_f$ massive hypermultiplets in the fundamental representation and theory with one adjoint massive hypermultiplet. If these theories are considered on $S^5$ their partition functions can be localized to matrix integrals, which can be approximated by their values at saddle points in the large-$N$ limit. We solve saddle point equations corresponding to the decompactification limit of both theories. We find that in the case of the fundamental hypermultiplets theory experiences third-order phase transition when coupling is varied. We also show that in the case of one adjoint hypermultiplet theory experiences infinite chain of third-order phase transitions, while interpolating between weak and strong coupling regimes.

that the free energy of this theory behaves as N 3 in the strong coupling limit. This result is consistent with the well-known N 3 behavior of the free energy of 6D theory obtained from the supergravity considerations [8]. Localization results were later generalized to theories living on the manifolds with more complicated geometries in [9][10][11], revealing the same N 3 behavior of the free energy with only difference in the prefactor.
Another 5D theory of interest is the super Yang-Mills with the U Sp(N ) gauge group and infinite coupling constant, which is conformal fixed point in five dimensions. This theory is especially interesting because it has known holographic dual in AdS 6 space [12]. This gauge theory was studied using localization as well in [13] and [14]. It was found that the free energy of this SYM theory has N 5/2 behavior in the planar limit. This result matches the result obtained from the AdS/CF T correspondence. Finally this 5D CFT was also studied on the backgrounds with more general geometries, e.g. on the squashed spheres, in [15,16]. These results allow one to evaluate supersymmetric Rényi in five dimensions [17].
One can also introduce a Chern-Simons term into SU (N ) SYM theory and consider the limit of the infinite Yang-Mills coupling. Then this theory is also 5D superconformal fixed point, which was studied using localization in [18]. It was found that in this case the free energy of theory also behaves as N 5/2 when N is large. This suggests the existence of the holographic dual of 5D supersymmetric Chern-Simons theory, though the dual is not known yet.
Using localization, it was recently found that different supersymmetric theories have a very interesting phase structure. The matrix models obtained from the localization of supersymmetric filed theories with the massive hypermultiplets were shown to experience a phase transition at some critical values of the couplings in the decompactification limit. First evidences of these phase transitions were obtained for 4D N = 2 * SYM on S 4 . It was found that in the limit of infinite radius this theory experiences infinite chain of phase transitions as coupling is varied from weak to strong. In particular, the authors of [19] solved the matrix model obtained from the localization of N = 2 * SYM. They have found that the support of the eigenvalue density in this case consists of many intervals of the length equal to the mass m of the adjoint hypermultiplet. On the boundaries of these intervals density distribution has cusps. When the 't Hooft coupling is increased the length of the support grows as well and at some points new cusps appears in the distribution. Emergence of new cusps signals about the phase transition at the corresponding point. Finally, at very strong coupling cusps of the distribution smooth out and the solution approaches the one obtained in [20]. This interesting phase structure of N = 2 * SYM on S 4 was investigated in details in the series of papers [19,21,22]. Later these results were generalized to the decompactification limit of N = 2 * SYM on the ellipsoids in [23].
Another 4D theory experiencing large-N phase transition is N = 2 SYM with 2N massive hypermultiplets in the fundamental representation considered on S 4 [21]. However, phase structure of this theory appears to be much simpler. It experiences only one phase transition, while interpolating between weak and strong coupling in the decompactification limit. From the matrix model point of view this phase transition takes place when the length of the support of the eigenvalue distribution becomes larger than two times the mass of the hypermultiplet.
Similar phase structures were also obtained in 3D supersymmetric theories. In [24,25] authors considered matrix model obtained from the localization of N = 2 supersymmetric Chern-Simons theory on S 3 coupled to 2N f massive hypermultiplets in the fundamental representation. It was shown that in the decompactification limit this theory experiences third-order phase transition similar to the one in N = 2 SYM with 2N massive fundamental hypermultiplets on S 4 . Finally these results were generalized to the case of massive deformations of the ABJM theory on S 3 [26,27], for which the phase structure of the theory appears to be similar to the case of N = 2 * SYM on S 4 . Namely, it was found that the theory experiences an infinite chain of third-order phase transitions, while interpolating between weak and strong coupling regimes. As in 4D case, these phase transitions are related to the emergence of new cusps in the distribution of the eigenvalue density.
In this paper we try to generalize the results of 3D and 4D to the case of 5D gauge theories 1 . We are particularly interested in the SU (N ) N = 1 SYM coupled to either one adjoint or N f fundamental massive hypermultiplets. To work with these theories we use the matrix model obtained by localizing N = 1 SYM on S 5 with arbitrary hypermultiplet content [3,4]. In general it is not possible to evaluate these matrix integrals. However in this paper we focus only on the large-N limit where the matrix integrals are dominated by their saddle points. In the case of adjoint hypermultiplet this saddle point equations have been solved in the weak and strong coupling limits [6,7]. Here we send the radius of the five-sphere to infinity, which will simplify the saddle point equations and allow us to solve them exactly at any coupling. When the solution is observed we can evaluate supersymmetric observables in N = 1 SYM. In particular, we consider the free energy and the expectation value of the circular Wilson loop. Studying both these observables we can draw conclusions about details of the phase structure of the theory.
This analysis leads to the following picture of the phase structure of considered theories. In the case of the theory coupled to N f fundamental hypermultiplets there is only one thirdorder phase transition taking place at where t = g 2 Y M N is the definition of the 't Hooft coupling used throughout in this paper and m is the mass of the hypermultiplets. Notice that in order to reach this critical point one should consider negative g 2 Y M . As explained in section 2 this does not contradict anything. Negative coupling arises in the theory due to the renormalization and doesn't spoil the convergence of the matrix integral.
The phase structure of the SU (N ) N = 1 SYM with the adjoint hypermultiplet is even more interesting. For this theory we observe an infinite chain of third-order phase transitions 1 Some phase transitions were already found before in 5D N = 1 Chern-Simons theory [18] and in 5D N = 1 SYM theory with adjoint hypermultiplet [28]. However these phase transitions are different from the transitions described in this work. at the following critical points t (n) c = 8π 2 m (n + 1) . (1.2) In the strong coupling limit these phase transitions smooth out and approach the strong coupling solution found in [7]. A nature of these phase transitions is the same as in the mass-deformed ABJM theory and in 4D N = 2 * SYM. As we show, the eigenvalue density at the saddle point of the matrix model has cusps. The number of these cusps grows as we increase the 't Hooft coupling, and each time a new cusps appear the system experiences a phase transition. This paper is organized as follows. In section 2 we review some properties of the matrix model obtained by localizing 5D SYM on S 5 . In particular we consider issue of renormalization of the coupling and find the decompactification limit of this model. In section 3 we solve the saddle point equations of the matrix model corresponding to the theory coupled to N f fundamental hypermultiplets in the decompactification limit. Then calculating the free energy of the matrix model we show that the theory experiences the third-order phase transition at the critical value of the coupling. In section 4 we solve the matrix model with adjoint massive hypermultiplet in the decompactification limit and show how the chain of phase transitions described above arise in this case. Finally, in section 5 we discuss the effect of the squashing of five-sphere on the phase transition. We find that in contrast to the 4D case discussed in [23], phase transitions are not affected by the deformations of the sphere. Some technical details of calculations are included in the appendicies of the paper.

Matrix Model
To study the phase structure of 5d SYM we use the result of supersymmetric localization [3,4], that reduces full field theory path integral to the finite-dimensional matrix integral given by Here g Y M is Yang-Mills coupling, k is Chern-Simons (CS) level, r is the radius of S 5 and integration variableφ is related to the expectation value of the scalar field σ from vector multiplet as the followingφ = −irσ. 2 . One-loop contributions of the vector-and hypermultiplets are given by and where β are the roots, µ are the weights of the representation R and m = −iM r with M being the mass of the the hypermltiplet. The matrix model described above was well studied in some regimes. In particular, planar limit of the SU (N ) SYM was studied in [6,7,29]. In [13,14] this model was studied for the case of 5d superconformal theory. And, finally, planar limit of the CS theory was studied in [18]. In general (2.1) should also include instanton contributions. However we, in this paper will be interested in the large-N limit of the theory, which suppresses instanton contributions. Hence we omit all instanton contributions and consider only classical action together with one-loop determinants in the partition function.
In this paper we put CS term to zero and concentrate on the decompactification limit of the SU (N ) SYM with different hypermultiplet content. In all cases hypermultiplets are massive which leads to the phase transitions similar to the ones obtained in [19,26,21]. Before moving to the solutions of the matrix model we should discuss general properties of this model. More detailed analysis of this matrix model can be found in [7].

Renormalization of the Coupling Constant
As we can see the one-loop contributions for both vector-and hypermultiplets can be written in the following form where P(x) is the infinite product This product is divergent with the divergent part given by where we have also omitted x-independent divergent part. To regularize the product (2.6) we introduce the UV cut-off at t 0 = πΛ 0 r leading to log P = −πΛ 0 r x 2 2 + 3x . (2.8) Using this expression we extract divergent terms in the one-loop determinants Here C 2 (R) is the quadratic Casimir element for representation R, defined by Tr( . Note that we have omitted terms linear inφ, because the gauge group is SU (N ). In this case the tracelessness condition implied forφ forces all linear terms to be zero. The divergent terms in the one-loop contributions (2.9) and (2.10) are proportional to Tr(φ 2 ). So in the full matrix integral (2.1) these divergent terms can be absorbed into the renormalization of the Yang-Mills coupling where the sum in the second term is over all hypermultiplets of the theory. Equation (2.11) reproduces the renormalization obtained in flat space using perturbation theory in [30]. Further we will consider theories with particular matter contents, leading to different effective coupling constants. In each case we will come back to (2.11) and consider particular examples of this renormalization.
After the regularization described above finite parts of the one-loop contributions (2.9) and (2.10) can be written in the form where F V (φ) and F H (φ) are functions defined as follows Here functions f (x) and l(x) are the ones introduced in [3] and [31] correspondingly. For the definitions and properties of these functions see appendix A. Using results of the regularization it is convenient to write down the matrix integral (2.1) in the following form (2. 16) Another observable that can be evaluated using localization technique is the expectation value of the circular Wilson loops where C is wrapping the equator of S 5 . Such loops were considered in [32] and [7] for 5D SU (N ) SYM theory and in [14] for 5D superconformal theories. After the localization expectation value (2.17)

Decompactification limit
The matrix integral (2.16) is not solvable even in the planar limit. Some particular limits of this matrix model were studied in [6,7] for the case of one adjoint hypermultiplet and N f fundamental hypermultiplets. In this paper we study the behavior of the theory with SU (N ) gauge group and different hypermultiplets content in the decompactification (infinite radius) limit. In all equations above we have used the dimensionless variablesφ andm, or equivalently we can assume that the radius of the sphere was set to one. To reproduce the r-dependence of the matrix integral we should rescale our variables as follows 3 : φ → φr ,m → mr . (2.19) The decompactification limit then corresponds to the limit r → ∞. Using asymptotic expressions (A.6) we can find that in this limit F (φ) in (2.16) can be written as (2.20) 3 From now on we will denote all dimensionless parameters with tilde and corresponding dimensionfull ones -without tilde Here we also included 1/r 2 terms linear in φ. These terms will be important when we discuss the effects of finite r in the large radius limit.
We see that in the decompactification limit the matrix model simplifies a lot, as (2.20) does not contain any complicated functions and, as we show further, can be solved explicitly in the planar limit. However, this is not the only reason to consider this limit. As mentioned in the introduction, studies of the decompactification limit in different 3D and 4D theories revealed interesting phase structure of these theories. For instance, phase transitions were found in 3D Chern-Simons coupled to the massive fundamental hypermultiplet [24], in 4D N = 2 super-QCD with massive quarks in the Veneziano limit [21], in mass-deformed ABJM theory [26,27] and finally in 4D N = 2 * SYM theory [19]. In each case theory was considered in the decompactification limit and it was shown that in the finite volume phase transitions disappear. We expect to observe similar picture in 5D SYM theory and hence we should also concentrate on the decompactification limit of the corresponding matrix integral (2.16).
All considerations above are general. They can be applied to any gauge group and hypermultiplet content. From now on we will concentrate on the SU (N ) gauge group, meaning that C 2 (adj) = N and C 2 (f und) = 1 2 . Inspired by the previous results on the phase structure of the supersymmetric gauge theories we consider two type of theories: N = 1 SYM coupled to N f massive fundamental hypermultiplets and N = 1 SYM coupled to one massive adjoint hypermultiplet.

N = 1 SYM with N f fundamental hypermultiplets
In this section we consider N = 1 SYM containing vector multiplet and N f hypermultiplets with the mass m. Half of the hypermultiplets is in fundamental, while another half is in antifundamental representation. 4 Before working out matrix model equations of motion and its solutions we should discuss the renormalization issues, as they play very important role in this case. According to (2.11) effective coupling in the case of N f fundamental hypermultiplets have the form where ζ ≡ N f /N is the Veneziano parameter. Note that the theory has constraint on ζ that can be found from the decompactification limit expression (2.20). We want F (φ) to be positive at large φ in oder for the matrix integral (2.16) to be convergent. In the case of the SU (N ) gauge group with N f /2 fundamental and N f /2 antifundamental hypermultiplets we obtain: where by the ellipses we mean all terms subleading in 1/|φ|. From this expression we see that in order to have well defined matrix integral, parameters of the theory should satisfy N f ≤ 2N or ζ ≤ 2 in terms of the Veneziano parameter. The same restriction was obtained from the flat space prepotantial in [33]. According to (3.1) YM coupling g 2 ef f can be negative, except in the case of ζ = 2. Further in our considerations we will always assume ζ < 2. Then we can redefine UV cut-off Λ 0 so that and use it in the matrix model further. Now we are ready to solve the matrix model with described matter content and coupling renormalization. In the planar limit the matrix integral (2.16) is dominated by the saddle point satisfying Using expressions (A.5) for the derivatives of F V and F H we obtain the saddle point equations (3.5) After taking the decompactification limit we arrive to This equation can be also observed from (2.20). Introducing usual eigenvalue density we rewrite the saddle point equation (3.6) as the integral equation (3.8) We assume that the eigenvalue density ρ(φ) has finite support [−a, a]. Then we should consider two separate cases. In the first case, which corresponds to the phase of the theory below the transition, a < m is satisfied. In the second case we have a > m, which corresponds to the phase above the transition point.
To solve equation (3.8) we differentiate it three times: where (3.9),(3.10) and (3.11) correspond to the first, second and third derivatives of the equation (3.8) respectively. Now we consider these equations for the cases a < m and a > m separately. This equation should be satisfied for any φ ∈ [−a, a], hence we conclude that ρ(φ) is sharply peaked near the endpoints of the distribution at φ = ±a. Then the solution to (3.12) is given by

Solution below the transition point
where the overall coefficient 1/2 comes from the standard normalization condition To be more precise the arguments of the δ-functions in (3.13) should look like (φ ± a ∓ ) with → 0. This small shift should be done in order to include the whole delta function into the density support [−a, a]. Everywhere further we imply these shifts of the δ-function arguments at the distribution endpoints. Note that (3.13) does not satisfy (3.11) at the endpoints of distribution. The best way to make the solution (3.13) more rigorous is to consider equations of motion (3.5) in the limit of large but finite r, taking into account terms subleading in 1/r. These terms correspond to the repulsion of the eigenvalues at small separations, which washes out δ-functions and turn them into the peaks of 1/r width. The form of these peaks can be found analytically. Then the solution (3.13) can be obtained in the limit r → ∞. Details of these calculations can be found in appendix C.
To determine the position of the endpoint a, we substitute solution (3.13) back into the original equation (3.8) which leads to By construction the solution (3.13) works when a < m or Λ < m 1 − 1 2 ζ . Notice that the r.h.s of (3.15) is always positive due to the condition ζ < 2, we have discussed previously.

Solution above the transition point
Now we increase Λ, pass the point Λ c = m 1 − 1 2 ζ and arrive to the phase above the transition point (Phase II), where we have a > m. To solve (3.8) we again use its three derivatives (3.9)-(3.11). First we address (3.11). In the case a < m terms δ(φ ± m) do not contribute into this equation, because arguments of these δ-functions are never zero. However, if a > m δ-functions start contributing into the eigenvalue density. Then it is natural to assume that solution contains these δ-functions on top of the solution (3.13) Here the coefficients in front of δ(φ ± m) are chosen to satisfy (3.11) and the coefficients in front of δ(φ ± a) are found from the normalization condition (3.14) .
The solution found here suffers the same problems as (3.13). To justify it we again consider equations of motion with large but finite r and take the decompactification limit in the very end. This calculation is done in appendix C, where we reproduce the eigenvalue density (3.16) in the limit r → ∞.

Free energy and the order of transition
We now determine if the transition between phases I and II is phase transition and , if so, what is the order of this transition. To answer these questions we find the behavior of the free energy near the critical point Λ c . To calculate the free energy in the decompactification limit we directly use the asymptotic formula (2.20) with an appropriate matter content. In the large-N limit we express the free energy in the integral form These integrals can be evaluated both below and above the phase transition point using corresponding expressions (3.13) and (3.16) for the eigenvalue density. After some calculations we obtain (3.20) where F (I) and F (II) are the free energies of phase I (Λ < Λ c ) and phase II (Λ > Λ c ) respectively. Using (3.20) we find that the free energy has a discontinuity in its third derivative at the critical point Thus at the critical point Λ c = 1 − 1 2 ζ m the theory experience third-order phase transition. Third-order phase transitions are typical for matrix models in the planar limit and has been obtained for many different systems [34,26]. As we have mentioned before there are examples of the phase transitions similar to the one we have described in this section. These are phase transitions in 3D Chern-Simons-Matter theory [24,25] and 4D super-QCD with 2N f massive hypermultiplets [21]. In both examples the phase transition of the theory in the decompactification limit at large-N is of third order. . The parameters of the theory are taken to be r = 40, m = 6, ζ = 3 4 , N = 100 corresponding to Λ c ≈ 3.77.
We have also checked results for the behavior of the free energy for different parameters numerically. In Fig. 2 we show the numerical and analytical results for the certain set of the parameters. The orange dots show the results of the numerical simulation, while the dashed blue lines represent the analytical result (3.20). As we see the analytical results are consistent with the numerical ones.
However from the Fig.2(b) we see that the third derivative, obtained numerically, is not exactly discontinuous but rather interpolates smoothly between two constant values, thus leading to the absence of the third-order phase transition. This happens because numerical calculations were performed for the large, but finite r. The finite volume of the system leads, as usually, to the absence of phase transitions, that are substituted by the crossover transition. However numerical study of the curve of ∂ 3 Λ F for different values of r have shown that it converges to the step function in the decompactification limit r → ∞, signaling about third order phase transition.

Wilson Loops
We also study the behavior of the circular Wilson loops near the critical point. In the planar limit the exponent term in the matrix model expectation value (2.18) does not affect the position of the saddle point, which leads to the following integral expression for the Wilson loop vev: Calculating this integral using the eigenvalue density ρ(φ) from (3.13) and (3.16) we arrive to for the phase I and for the phase II. Using these expressions we find the discontinuity in the second derivative of the Wilson loop expectation value with respect to Λ at the critical point Λ c which is consistent with the third-order phase transition.

N = 1 SYM with adjoint hypermultiplet
In this section we consider phase structure of N = 1 SYM with the massive adjoint hypermultiplet of the mass m. As we will see this phase structure consists of infinite number of third-order phase transitions on the way from weak to strong coupling.
In the case of the adjoint hypermultiplet using (2.11) we obtain i.e. in this case the coupling does not get shifts from the UV cut-off Λ and from now on we use the bare coupling constant g Y M .
In the planar limit matrix integral (2.16) is dominated by the saddle point Using expressions (A.5) for the derivatives of F V and F H we obtain the following equations of motion Restoring the radius dependence (2.19) and taking the decompactification limit, we obtain These equations can be directly observed from (2.20). Their continuous limit has the following form where we have introduced 't Hooft coupling constant t = g 2 Y M N 5 . To solve (4.5) we proceed in the same way as while solving (3.8) for the N = 1 SYM with the fundamental hypermultiplets. We take three derivatives of (4.5) , thus reducing integral equation to the algebraic one. Then we find general solutions to this algebraic equation for different phases of the theory, corresponding to different values of the parameters. Finally, we fix all integration constants and free parameters by substituting general solution back into original equation. Three derivatives of (4.5) are given by Before we start solving these equations let's understand what do we expect to obtain. In analogy with 3D and 4D results [26,19,21] • We start with the theory with 2a < m, where a is the position of the support endpoint for the solution of (4.5).
• Then we increase the 't Hooft coupling t which leads to the widening of the support (increase of a). At some value of t, where 2a = m, we obtain emergence of two resonances, originating from the nearly massless hypermultiplets, i.e. from the terms |φ − ψ| ≈ m in the saddle point equations (4.5). Appearance of these resonances leads to the change of the form of the eigenvalue distribution, signaling about the emergence of a new phase.
• If we further increase coupling new secondary resonances appear every time when 2a = nm with n ∈ Z. Each time this happens we expect to obtain phase transition in the theory.
• In the limit of very large 't Hooft coupling t → ∞ we expect to obtain an infinite number of resonances that will be averaged to the strong coupling solution obtained in [6,7].
Our final goal here is to obtain solutions of (4.5) for arbitrary number of resonances n. However in order to understand how to construct this general solution we start with finding explicit solutions for n = 0 and n = 1 and only then generalize solution to general n.

Simple examples
The phase with no resonances corresponds to the case of 2a < m. To solve (4.5) with these values of the parameters we address (4.7) corresponding to the second derivative of the original equation. This equation then is equivalent to (3.12) and hence can be solved by Here and further below we treat the δ-functions at the endpoints of the distribution in the same manner as we were doing it in the case of the fundamental hypermultiplets, i.e. we always assume the small shift in the argument of the δ-function that brings whole δ-function inside the eigenvalue support. Substituting ansatz (4.9) into the original equation (4.5), we obtain the following relation for the position of the support endpoint (4.10) By construction this solution is valid when 2a (0) < m or, equivalently t < 16π 2 m . Thus the first critical point is Similarly to the case of the fundamental hypermultiplets to justify this solution we consider solutions to the full equations of motion (4.3) in the limit of large but finite r, and then obtain (4.9) together with (4.10) in the limit of infinite r. For the details of this calculation see appendix D.
We also found numerical solutions to (4.3). We show these solutions in Fig.3(a) with the orange dots. The dashed blue lines on the same plot represent positions of the δ-functions in our analytical solution (4.9). As we see from this picture analytical solutions reproduce numerical ones very well. All differences between these solutions appear due to the effects of finite r, which are discussed in appendix D.
Now we increase the 't Hooft coupling and pass the critical point t (1) c , so that now m < 2a < 2m. Under this condition we expect two resonances to appear at φ = ±(m − a). To solve equations of motion (4.5) for this case we divide the eigenvalue support [−a, a] into three regions and find general solution to (4.5) in every region. We also use φ → −φ symmetry of (4.5), that leads to where we used ρ(φ + m) = 0 as φ + m > a and also (φ − m) ∈ C. Using the symmetry property (4.13) we obtain This equation can be satisfied only if ρ A (φ) = 0. In analogy with the previously found solutions we assume that this is true up to isolated points on the boundary of the interval at φ = a, m − a. Then the natural ansatz is Region B (−m + a < φ < m − a). In this region (4.7) takes the form The first and the last terms in this expression cancel each other due to the symmetry properties (4.13) of ρ(φ). The remaining part of the equation then reads suggesting, similarly to the case of n = 0, that solution for the eigenvalue density ρ B (φ) is sharply peaked around endpoints of the interval B at φ = ±(m − a) Note that we have already took one of this δ-functions into account in ρ A (φ) in (4.16), while the second δ-function appears in ρ C (φ) as we will see further.
Region C (−a < φ < a − m). Finally, on the last interval C we do not need to solve any equations in order to find the eigenvalue density. Instead we utilize our previous result (4.16) for the density in the region A together with the symmetry properties (4.13) to obtain (4.20) Notice here that, just as we expected, the second δ-function from ρ B (φ) showed up here. Summing up all our results we find the following ansatz for the eigenvalue density in the phase with one pair of the resonances  Finally substituting these values back into (4.23) we obtain the following expression for the endpoint position Corresponding eigenvalue density is given by By construction this solution is valid, when m < 2a (1) < 2m or equivalently 16π 2 m < t < 24π 2 m . Thus the second critical point is given by As in the previous cases we justify this solution by finding the solution of (4.3) for large but finite r and taking r → ∞ limit in the very end of the calculations. As the result we reproduce solution (4.27) with the right positions and the coefficients of the δ-functions. For the details of this calculation see appendix D. In Fig.3

General solution
Now, using examples of the solution with n = 0, 1 number of resonance pairs, we are ready to find a general solution with an arbitrary number n of the resonance pairs. This solution is valid when inequality mn < 2a < m(n + 1), where n ∈ Z, is satisfied. In this case the cut can contain n pairs of the resonances. The mechanism of the emergence of these resonances is exactly the same as it was in the case of n = 1. Primary resonances appear at φ = ±(a − m) due to the discontinuities in the kernel of (4.5), caused by the peaks at the endpoints of the eigenvalue distribution ψ = ±a. These resonances in turn create secondary resonances at φ = ±(a − 2m) by the same mechanism and so on up to the last resonance pair at φ = ±(a − nm) that fits inside the eigenvalue support.
Thus generalizing solutions (4.9) and (4.21) to the arbitrary number of resonances we expect to obtain the eigenvalue density of the following form (4.29) Now we proceed in the same way as for the cases of n = 0, 1. Namely we substitute ansatz (4.29) into equation (4.6). This substitution results in where in the last line we have shifted the summation index for the terms |φ ± a ∓ m(k + 1)| and |φ ± a ∓ m(k − 1)|. Assuming that c −1 = c n+1 ≡ 0 we arrive to the relatively compact equation where for simplicity we introduced Notice now that for the last two terms the following relations are satisfied everywhere on the support [−a, a] |φ − a + m(n + 1)| + |φ + a − m(n + 1)| = 2(m(n + 1) − a) , |φ + a + m| + |φ − a − m| = 2(m + a) . However the terms in the summation in (4.31) depend on φ and, in analogy with the previously considered cases, we expect the coefficients of these terms to be zero. To prove this statement precisely we follow the same algorithm as in the case of n = 1. However, we should consider cases of even and odd number of the resonance pairs separately. Though they are similar and lead, as we will see soon, to the same result, there are some differences that should be briefly discussed.  In the case of the even n the entire support [−a, a] should be separated into different regions with the boundaries corresponding to the resonances as shown in Fig.4(a). There are two types of intervals we should consider. Intervals of different types are shown with different colors in Fig. 4(a). On this picture and everywhere further for shortness we introduce the notation a k = a − mk.

35)
Intervals −a n 2 +j < φ < a n 2 −j , j = 1, .., n 2 (green intervals in Fig.4). In these regions similarly to the previous case we obtain 2φ , 36) and the corresponding equations of motion (4.31) turns into (4.37) Now the algorithm of defining all (n + 1) coefficients c k is the following. First condition comes from the normalization condition (3.14), that reads for the ansatz (4.29) as 2 n k=0 c k = 1 . (4.38) We then need n more conditions to determine all the coefficients. These conditions can be obtained by considering the intervals shown in Fig.4(a). We start from the region closest to the origin, i.e. a n 2 < φ < −a n 2 +1 , corresponding to j = 0 in (4.35). Then the first sum in this equation is given just by one term 2 n 2 +j As before we assume that this φ-dependence should be canceled leading to C n 2 = 0. Then we move away from the origin to the next interval −a n 2 +1 < φ < a n 2 −1 . This interval corresponds to the equation (4.37) with j = 1, which turns the first sum into the sum of two terms 2 n 2 +j thus resulting in the condition C n 2 + C n 2 +1 = 0. Repeating this procedure for all n intervals on the positive φ half-axis we obtain n equations n k=1 C k = 0 . (4.41) For the case of the odd n all considerations are analogous. As in the previous case we consider the division of the whole eigenvalue support into two types of regions, which is shown in Fig.4(b).
On the intervals an+1 2 −j < φ < −an+1 2 +j , j = 1, .., n−1 2 (blue intervals in Fig.4) equation (4.31) can be written in the form while on the intervals −an+1 2 +j < φ < a n−1 2 −j , j = 0, .., n−1 2 (green intervals in Fig.4) the same equations can be rewritten as the following  Now as in the case of even number of the resonance pairs we consider intervals and corresponding equations (4.42) and (4.43) one by one. To obtain n conditions for the coefficients C k we start with the region −an+1 2 < φ < a n−1 2 . This interval corresponds to the equation (4.35) with j = 0 resulting in C n+1 2 = 0. Repeating procedure for every single interval on the positive half-axis we obtain the same equations (4.41) for the coefficients c k . Now employing equations (4.38) and (4.41) let's determine the coefficients c k in our ansatz (4.29). To do this notice that the general solution for the recurrence relation c k+1 − 2c k + c k−1 = 0 , k = 1, .., n . (4.44) is given by The easiest way to see this, is to write the recurrence relation for large k in form of the differential equation d 2 dk 2 c(k) = 0 with the general solution given by c(k) = α + βk. To determine the coefficients α and β we use boundary condition c n+1 = 0. Substituting k = n + 1 into (4.45) we obtain α = −β(n + 1) . (4.46) Now substituting (4.45) into normalization condition (4.38) we find n k=0 c k = α(n + 1) + β n(n + 1) 2 = 1 2 . (4.47) Finally combining this relation with (4.46) we obtain α = 1 n + 2 , β = − 1 (n + 1)(n + 2) . (4.48) The last ingredient we need for the completeness of the solution is the position of the cut endpoint a. To determine it we should come back to equation (4.31). As we have shown C k = 0 for k = 1, .., n, so that corresponding terms do not contribute into equations. Using α and β determined in (4.48) we obtain expression for C 0 Substituting it into (4.31) we find the position of the support endpoint a (n) = (n + 1)m − 4π 2 t (n + 1)(n + 2) . Obtained solution is valid, when mn < 2a (n) < m(n + 1) or equivalently when 8π 2 m (n + 1) < t < 8π 2 m (n + 2) , (4.51)   c k (δ(φ − a + mk) + δ(φ + a − mk)) , c k = n + 1 − k (n + 1)(n + 2) , (4.53) with the position of the eigenvalue support endpoint a given by (4.50). If we take n = 0, 1 in this general solution we will reproduce corresponding formulas (4.9), and (4.27) obtained previously. We have also checked general solution (4.53) numerically. In

Free energy and the order of the phase transition
Now we are ready to address the question of the order of the phase transitions taking place at the critical points t (n) c . For this we evaluate the free energy corresponding to the eigenvalue density (4.53). We read the free energy in the decompactification limit directly from (2.20) choosing appropriate matter content. In the continuous limit it can be written as follows where we have also introducedF , corresponding to the free energy with excluded dependence on the radius r and rank N of the gauge group. Integrals in the expression (4.54) can be evaluated after substituting our solution (4.53). The first integral is relatively easy to evaluate where in the last step we have used explicit expressions for the coefficients c k (4.53) and the endpoint position a (n) (4.50). The second integral coming from the one-loop determinants is more complicated. Let's consider it in more details here c k c 0 m 3 (k + 1) 3 where we have also used inequalities mn < a < m(n + 1). This expression is much easier to work with, because, as we know, C k = 0 for k = 1, .., n. Then substituting values of the coefficients c 0 = 1 n+2 , c n = 1 (n+1)(n+2) , C 0 = − 1 n+1 and the endpoint position a (n) from (4.50) we finally obtain I (n) 2 = − 1 12t 3 512(1 + n) 2 (2 + n) 2 π 6 − 64mtπ 4 (1 + n)(2 + n)(3 + 2n) + m 3 t 3 (3 + 2n) . (4.58) The free energy of the system is then given by the sum of (4.55) and (4.58) (4.59) To find the order of the phase transition we calculate derivatives of the free energy (4.59) with respect to the 't Hooft coupling at the critical points (4.52). After short computation we observe (4.61) Hence, when 't Hooft coupling t is increased theory undergoes through an infinite chain of phase transitions of the third order. Moreover, with the increase of the coupling transition becomes weaker transforming into crossover transitions or infinite 't Hooft coupling t. Similar behavior was obtained in the case of the mass-deformed ABJM theory [26,27].   We have also evaluated the free energy for the full matrix model (2.16) with massive adjoint hypermltiplet numerically. The results of these numerical simulations are shown in Fig.6 with the orange dots. The dashed lines on the same plots represent our analytical result (4.59) for the free energy. As we see numerical and analytical results are in good agreement. The only difference can be obtained in the graph of the third derivative of the free energy. As we see the free energy, obtained numerically, does not have discontinuities at the critical points and hence leads to the absence of the phase transitions. This result is expected because the numerical solution is found for the finite r and, as known, finite-dimensional systems do not experience phase transitions.

Wilson Loops
We also study the behavior of the circular Wilson loop near the phase transition points. In the planar limit contribution of the e 2πφ prefactor in (2.18)  Substituting solution (4.53) into this expression we easily obtain where we have used expression (4.50) for the position of the endpoint. Calculating derivatives of W (n) with respect to the coupling constant around the critical point t

The strong coupling limit
It is interesting to study the behavior of the solutions (4.53) for a very strong coupling constant, such that 2a m. In terms of our solution, this means that the eigenvalue density contains large number of resonance pairs n 1. Previously in [7] we have obtained that for the strong coupling eigenvalue density at the saddle point is where the endpoint position is given byã = (9 + 4m 2 ) λ/64π 2 and λ is the coupling constant defined as λ ≡ g 2 Y M N/r = t/r. Restoring r-dependence (2.19) and taking the decompactification limit r → ∞ we obtain where the cut endpoint position is given by Now let's compare the expressions above with the corresponding limits of our solution (4.53) and (4.50). In order to do this recall that the solution with n pairs of resonances is valid when 8π 2 m (n + 1) < t < 8π 2 m (n + 2). Thus for very large n we can approximate the number of the resonance pairs by n ≈ mt 8π 2 . For the position of the endpoint using (4.50) we can easily obtain which exactly reproduces the expected result (4.67).
In order to reproduce the eigenvalue density (4.66) we consider the case of large n in (4.53). As n 1 the cut contains a large number of the δ-functions which is then regarded as a continuous distribution. In order to understand which distribution is this, let's consider the following integral where we have used (4.53) for ρ (n) (φ) and c k , and approximated the sum over k by the integral assuming that the sum goes to large enough values of k. As we obtained linear dependence of the integral   [7]. Similar behavior was obtained for the solutions containing resonances in the mass-deformed ABJM theory [26,27] and 4D N = 2 * theory [19]. The same conclusion about behavior of (4.53) at very strong coupling can be drawn numerically. On the plots in Fig.7 the orange dots represent numerical solutions to (4.3). On the same plots we show the solution (4.66) with the dashed blue line. As we see distribution for t = 300 still looks like a bunch of peaks, while for t = 750 it already starts approaching (4.66) and for t = 1250 completely coincides with it up to small regions of the size ∼ 1 r around the endpoints of the distribution. Hence, numerical simulations support the conclusions of this section.

Squashing spheres
Another issue to discuss is the effect of the geometry on the phase transitions described in this paper. Localization methods can be generalized to theories on the space with geometry more complicated than the spheres. In particular 5D SYM was considered on squashed spheres [35,36], Y p,q spaces [9,10] and even more general toric Sasaki-Einstein manifolds [11]. Here we focus on the squashed spheres, which can be embedded into C 3 as follows Then the corresponding matrix model obtained after localizing 5D SYM on this space is where ω ≡ (ω 1 , ω 2 , ω 3 ) and ρ is the volume factor of the squashed sphere, which is given by the ratio of the volumes of squashed and usual spheres Finally, S 3 denotes triple-sine function that can be defined through the following infinite product (n 1 ω 1 + n 2 ω 2 + n 3 ω 3 ) ((n 1 + 1) ω 1 + (n 2 + 1) ω 2 + (n 3 + 1) ω 3 − x) .(5.4) The matrix model (5.2) was partially studied for the case of 5D SCFT, i.e. SYM theory with the U Sp(N ) gauge group and infinite coupling constant g Y M [15,16]. One of the reasons to study supersymmetric theories on the squashed spheres is that their partition functions are related to the supersymmetric Rényi entropy. The latter can be evaluated from the partition function of theory on the n-covering of sphere, which is equivalent to the squashed sphere with the squashing parameters ( 1 n , 1, 1) [17]. In general the matrix model (5.2) is very complicated and can not be solved directly. However, we are only interested in the behavior of the matrix model solutions in the decompactification limit r → ∞ which implies that the arguments of the triple-sine function are large,φ 1. Thus, we can use asymptotic expression for these functions [16,36] log where ω tot. ≡ ω 1 + ω 2 + ω 3 . In the case of the matrix model (5.2), x =φ = rφ and we keep only cubic terms so that the partition function can be rewritten in the following form From this expression we can see that the squashing does not affect our conclusions about phase transitions neither for fundamental nor for adjoint hypermultiplets. The only contribution of the squashing comes from the volume prefactor of the free energy. Similar results were obtained in the strong coupling limit of 5D N = 1 SYM with adjoint hypermultiplet on Y p,q spaces [9] and toric Sasaki-Einstein manifolds [11]. In particular, it was shown that for both cases the only difference in the prepotential of the theory between these spaces and five-sphere is in the overall volume factor. Note that effects of the squashing on phase transitions were also studied for the case of 4D N = 2 * SYM in [23]. However results obtained in that paper differs from ours. The main conclusion is that the squashing of the sphere shifts all critical points of the phase transitions chain to the smaller coupling constant. While in our case as we have seen above squashing does not affect positions of the critical points.

Discussion
In this paper we studied the decompactification limit of the matrix model obtained from 5D N = 1 SYM theory with different hypermultiplet content on S 5 . First we considered N f massive fundamental hypermultiplets coupled to N = 1 SYM theory. Solving the planar limit of the corresponding matrix model we found that the eigenvalue density ρ(φ) at the saddle point consists of an isolated peaks. These peaks are located at the boundaries of the distribution at φ = ±a and at the points φ = ±m. Hence, we concluded that there are two phases of the theory. One phase corresponds to the case when the support includes φ = ±m and another one when it does not. Transition between these two phases takes place at the critical point where t = g 2 Y M N is the 't Hooft coupling used in this paper. Evaluating the free energy we find that the transition between these two phases is a third-order phase transition.
Second theory we consider is N = 1 SYM coupled to the massive adjoint hypermultiplet. In this case we also find that solution for the eigenvalue density consists of isolated peaks. These peaks are located at φ = ±(a − mn), where a is the endpoint of the distribution support, m is the mass of the hypermultiplet and n is an integer number. As the coupling is increased the length of the support grows and it can include more and more peaks. As the result, theory goes through an infinite number of phases which differ by the number of peaks in the solution. Transition between phases with n and n + 1 pairs of peaks take place at Free energy calculation showed that these transitions are third-order phase transitions. In the limit of very strong coupling the peaks, forming the eigenvalue density, smooth out and solution becomes equal to the strong coupling solution derived in [18]. It has been also shown that the effect of the five-spheres squashing comes in the form of an overall prefactor of the free energy in the decompactification limit of the matrix model, while saddle point remains the same. Hence, we concluded that the squashing does not affect neither the positions of the critical points nor the order of the phase transitions.
Finally, we note that the phase structure described in this paper seems to be general for the supersymmetric theories with massive hypermultiplets, as they arise in D = 3, 4, 5 dimensions. Hence, in the future it would be interesting to understand the nature of these phase transitions better.

A Properties of the functions f (x) and l(x)
In this appendix we consider properties of the functions f (x) and l(x) that appear in the expressions for the regularized one-loop contributions (2.12) and (2.13). Function f (x),first introduced in [3], can be written in th following form while l(x), first introduced in [31], is defined by Using these expressions we can derive the following properties of l(x) and f (x): 1. l(x) is an odd function and f (x) is an even function 2. The derivatives of the functions are given by 3. The asymptotic behavior of the functions is given by Using these properties we derive similar properties for F V (x) and F H (x) functions in (2.14) and (2.15): 1. Both F V (x) and F H (x) are symmetric functions 2. The derivatives of the functions are given by 1 4 + (m + x) 2 tanh(π(x + m)) . (A.5) 3. Finally the asymptotic behavior of these functions is given by

B Details of numerical analysis
To obtain numerical results presented in this paper we have used the method introduced in [37] in order to solve the ABJM matrix model. It was also adopted for solving 5D SYM matrix model [7] and 5D Chern-Simons matrix model [18]. Instead of solving the the saddle point equation ∂F ∂φ i = 0, where F is the prepotential, we introduce dependence of the eigenvalues φ i on the effective "time" variable t. Then we substitute our original equation with the effective "heat" equation At large time t solution of this differential equation approaches the solution of the original saddle point equation, provided the choice of τ 's sign is made properly. However, solving (3.5) which corresponds to the case of the fundamental hypermultiplet, we faced some problems. Our numerical method appeared to be unstable and very sensitive towards initial conditions due to the repulsive central potential ∼ −Λφ 2 .
To avoid this problem we used the fact that the renormalization (2.11) of the coupling constant g Y M for the case of N = 1 SYM can be reproduced by considering the case of one adjoint hypermultiplet with very large mass. The easiest way to realize this is to consider (4.3). If the mass of the hypermultiplet is large we can rewrite this equation in the following form where we used tanh π(φ i −φ j ±m) = sign (±mπ) = ±1 and also where t = g 2 Y M N is the 't Hooft coupling. So in order to avoid the problems with the repulsive central potential we introduce adjoint hypermultiplet with large enough mass m. In this case the central potential itself is attractive as g 2 Y M > 0. However the effective coupling (B.3) can be made negative as we need it in equation (3.5). All numerical simulations of the theory with the fundamental hypermultiplets presented in this paper were performed in this way.

C Solution for the finite r: fundamental matter
In this appendix we consider the solution of (3.5) for the large enough but finite r. In section 3 we have neglected all terms subleading in 1/r in this equation in order to consider properly the decompactification limit r → ∞. As the result we have obtained solutions (3.13) and (3.16) containing δ-functions. These δ-functions arise in the solutions because of the kernel in (3.8) that looks like −(φ − ψ) 2 sign (φ − ψ) leading to the attractive potential between the eigenvalues. If not the repulsive central potential, the eigenvalues would just all condense at the origin φ = 0.
Here we will be more accurate with the kernel of integral equation. We now do not neglect the subleading term in the kernel of (3.5), which leads to where ρ(φ) is the eigenvalue density defined by (3.7). Due to the symmetry of (C.1) we expect the eigenvalue density to be symmetric. Note that terms subleading in 1/r are repulsive, hence they wash out δ-functions into peaks of the finite width ∼ 1/r. At the same time we have not included terms subleading in 1/r into the central potential, because they add small central repulsion, thus never playing an important role. Assuming coth (πr(φ i − φ j )) ≈ sign(φ i − φ j ) we also neglect terms exponentially suppressed in large r. However, note that this approximation works pretty well even in the region, where |φ i − φ j | ≈ 1/r due to the factor of π in the argument of coth. Now we consider phases I and II separately.

C.1 Below the transition point
We start with the phase below the transition point which corresponds to a < m, where a is the position of the eigenvalue support endpoint. In this case equation (C.1) and its derivatives are given by General symmetric solution to (C.5) is Now we are left to determine constant C and the position of the support endpoint a. To determine these constants we substitute the general solution (C.6) into normalization condition (3.14) This leads to C = r 2 sinh(ra) , (2Λ + ζm) r 2 − 2ar 2 + 2r coth(ra) = 0 . (C.9) The consistency of these conditions with equations (C.2) and (C.4) can be checked by the direct substitution. In general (C.9) can be solved only numerically for different values of the parameters r, m, ζ and Λ. However, for the large radius, when ra 1, we can simplify the last equation which reproduces the support endpoint (3.15) we found in section 3. Finally, the eigenvalue density below the phase transition point Λ c 6 is (C.11) For the large r this distribution is the function with two sharp peaks of width ∼ 1/r at φ ± a. In the decompactification limit r → ∞ this density is equivalent to (3.13), because the width of the peaks goes to zero, while the integral of ρ(φ) stays 1.
We have also found the numerical solutions to (3.5), which are shown in Fig.8(a) with the orange dots. Dashed blue lines shows the solution (C.11) with the coefficient C and the endpoint a found numerically from the equations (C.9). As we see, solutions found above are consistent with the numerical results.

C.2 Above the transition point
Now we find the solution to (C.1) for the case of the system in the phase II. This phase corresponds to a > m, and to solve equations of motion we should divide the whole eigenvalue support [−a, a] into three parts A, B, C with the symmetry constraints Now we write down equation of motion (C.1) in regions A and B and solve them in this regions separately. Region A (m < φ < a). In this interval (C.1) and it's first, second and third derivatives are As we do not have any symmetry restrictions on ρ A , the general solution to (C.17) is given by Note that we choose solution to be centered around φ = (a + m)/2. This choice is made for convenience. The shifts in cosh and sinh can be adsorbed into constants C 2 and C 3 Region B (−m < φ < m). On this interval φ < m, so that everything is the same as in the case of the phase I. Equation of motion and its three derivatives then can be written in the form of (C.2)-(C.5). Similarly to the phase I we obtain Region C (−a < φ < −m). In this region using symmetry properties (C.13) we find ρ C (φ) = C 2 cosh rφ + r a + m 2 − C 3 sinh rφ + r a + m 2 . (C.20) To complete our solution, we need to find all the parameters C 1 , C 2 , C 3 , a in (C.18),(C. 19) and (C.20). To determine these parameters we substitute solutions (C.18)-(C.20) into equations (C.2)-(C.4), (C.14)-(C.16) and normalization condition (3.14). This leads to where (C.21) is obtained from the normalization condition, while (C.22),(C.23) and (C.24) are obtained after substituting the solution for the eigenvalue density into (C. 16), (C. 15) and (C.3) correspondingly. It can be checked by the direct calculation that remained equations (C.2),(C.14),(C.4) will be satisfied automatically if the coefficients C 1 , C 2 , C 3 , a obey equations we have found above.
To summarize we have found that above the phase transition point, when m < a, the eigenvalue density satisfying equation of motion (C.1) can be written in the form (C. 25) with coefficients C 1 , C 2 , C 3 and the endpoint a of the cut satisfying equations (C.21)-(C.24). Equations (C.21)-(C.24) can be solved exactly in the limit of large radius r. In the case, when rm 1 and r a−m 2 1 we approximate these equations with simpler ones C 1 exp (mr) + 2C 2 exp r a − m 2 = r , This system can be solved leading to ζr exp(−mr) , Note that to find these solution we assumed r a−m 2 1, which means that solutions only makes sense far enough from the critical point where a ≈ m. However this approximation will work very well in the decompactification limit even for the points close to the critical ones.
As we see from (C.27) the position of the cut endpoint is consistent with the one obtained in the section 3. If we now come back to the eigenvalue density (3.13) and take into account coefficients (C.27) we can see that in the decompactification limit r → ∞ solutions are nonzero only at the points φ = ±a and φ = ±m, where they take infinite value. At the same time from the normalization condition (3.14) we know that ρ(φ) integrates to one. Then the natural function to describe this limit is (3.16) and the coefficients in front of the δ-functions should be which leads exactly to (3.16).
We also can find solutions to (C.21)-(C.24) numerically for the particular combinations of the matrix model parameters. In Fig.8(b) the orange dots shows numerical solution to (3.5), while the dashed blue lines shows the solution (C.25) with all coefficients found numerically from (C.21)-(C.24). As we can see numerical solution perfectly coincides with the solution found in this section.

D Solution for the finite r: adjoint matter
In this appendix we consider solution to the saddle-point equation of the matrix model (2.16) for the large but finite r. We have solved these equations in the decompactification limit (see section 4) neglecting all terms subleading in 1/r. As the result we have obtained the general solution (4.53) containing the δ-functions. In order to show more precise how these δ-functions arise we should take into account first subleading terms in the kernels of (3.5), which leads us to where we restore the radius dependence (2.19) and assume that coth(πr(φ−ψ)) ≈ sign(φ−ψ) as well as tanh(πr(φ − ψ ± m)) ≈ sign(φ − ψ ± m) when we consider the limit r 1. The same equation can be obtained by a minimization of the free energy (2.20).
Note that subleading terms correspond to the repulsive force between the eigenvalues at ∼ 1/r distances. Hence these terms wash out δ-functions into the peaks of 1/r width. As we will see further, it is possible to obtain the analytic form of these peaks and reproduce (4.53) in the decompactification limit r → ∞.
To obtain solutions of (D.1) we first reduce this integral equation to the differential one taking three derivatives where equations (D.2),(D.3) and (D.4) corresponds to the first, second and third derivatives w.r.t. φ. Now we are ready to consider solutions for different phases of the theory. In particular, we consider only solutions with no resonances and one pair of resonances. In the decompactification limit r → ∞ we want to observe solutions (4.9) and (4.27). Unfortunately it is very hard to observe solution for the case of arbitrary number of resonance pairs n. However obtaining particular cases of n = 0 and n = 1 will strongly suggest that solutions we have found in section 4 are right.

D.1 Solution with no resonances (n = 0)
We start with the case for which we do not expect any resonances. This happens when 2a < m, where a is the position of the support endpoint as usually. In this case ρ(φ ± m) = 0 as φ ± m appears to be outside the support of the eigenvalue density. Hence (D.4) can be written in the form General solution to this differential equation is given by where C is an integration constant. Notice that we omit sinh(rφ) due to the symmetry ρ(φ) = ρ(−φ).
Finally, substituting this solution into (D.2) we obtain the following condition for the position of the support endpoint a m − 8π 2 t r − ar + coth(ra) = 0 , (D.8) To summarize, we have found that solutions for the eigenvalue density of N = 1 SYM with massive adjoint hypermultiplet and no resonances is given by which exactly reproduces (4.10) found after neglecting the subleading terms from the very beginning. The eigenvalue density (D.9) becomes highly peaked around φ = ±a with the width of the peak ∼ 1 r . In the decompactification limit r → ∞ (D.9) can be considered as the sum of two δ-functions 1 2 (δ(φ + a) + δ(φ − a)), because the peaks becomes infinitely high and narrow and integrate to one. So we see that (D.9) in the decompactification limit coincides with the solution (4.9) found in section 4.
We can also compare the solution found here with the results for the numerical simulations of the full equation of motion (4.3). In particular in Fig.9(a)   Interval A (m − a < φ < a). Note that on this interval φ + m > a so that ρ(φ + m) = 0 and also −a < φ − m < a − m meaning φ − m ∈ C. Using these remarks and symmetry properties (D.12) we can rewrite (4.8) as the following delay differential equation where coefficients α and β where found by substituting these solutions back into (D.13). This procedure is similar to writing out characteristic equation for ODE, however still a bit different due to delays we have in our equation. In particular single exponents can not satisfy (D.14) and we should instead search for the solution in form of cosh αr φ − m 2 and sinh βr φ − m 2 . Shifts in the arguments can be seen as centering of the solution at the middle of the interval A which is made for the convenience and in principle can be adsorbed in constants C 2 and C 3 .
Interval B (a − m < φ < m − a). On this interval situation is simpler then on A, because φ − m < −a and φ + m > a so that both arguments are outside the support and hence ρ(φ ± m) = 0 leading to the following form of (D. where symmetry requirements (D.12) were taken into account.
To summarize we have found = C 2 cosh αr φ − m 2 + C 3 sinh βr φ − m 2 , m − a < φ < a , ρ(φ) = C 1 cosh(rφ) , a − m < φ < m − a , = C 2 cosh αr φ + m 2 − C 3 sinh βr φ + m 2 , − a < φ < a − m , = 0 , |φ| > a , α 2 = 8 17 , β 2 = 8 5 . (D.17) Now we should determine integration constants C 1 , C 2 , C 3 and the position of the support endpoint a using normalization condition (3.14) together with equations (D.1)-(D.3), written out for different intervals. After some algebra we obtain the following system of algebraic equations Assuming that r(a − m 2 ) 1 and r(m − a) 1, or equivalently assuming that the radius of S 5 is large and we are far enough from the critical points, we can rewrite these equations in the following form Solving these equations we find the following expressions for the integration constants C 1 = re −r(m−a) 1 − 2 1 + β −1 α −1 + 3β −1 + 4 , C 2 = 2rαe −αr(a− m 2 ) 1 + β −1 α −1 + 3β −1 + 4 , The last unknown parameter of the solution is the support endpoint position a. To find it we substitute solutions for C 2 and C 3 into (D. 25). Note that the last parenthesis containing C 2 and C 3 is of order r −1 . Hence in the decompactification limit this term is irrelevant and we are left with the relation a = 2m − Indeed the function f (r) is peaked at φ = c for large r and in the limit r → ∞ this peak becomes infinitely narrow and infinitely high, while the integral of f (r) is normalized to 1. Thus in the decompactification limit this function can be considered as the δ-function. We have already used this relation with α = 1 while considering decompactification limit of the solution (D.6) with no resonances. The distribution (D.17) has similar exponential peaks at φ = ±a and φ = ±(m − a) and thus we expect to obtain the following eigenvalue density in the decompactification limit where we have also used solutions (D.26) for the coefficients C 1 , C 2 , C 3 . As we see in the decompactification limit we completely reproduce solution obtained in section 4. Finally, we compare analytical solution (D.17) with the numerical results. We show this comparison in Fig.9(b). The orange dots represent the results of the numerical solution, while dashed blue lines on the same picture show analytical solution (D.17) with the parameters C 1 , C 2 , C 3 and a evaluated numerically from the system of equations (D.18)-(D.21). As we see from this picture analytical results coincide with the numerical ones perfectly. Small discontinuities at the resonances positions φ = ±(m − a) are related to high sensitivity of the exponential factors in this solution as well as in the system of equations (D.18)-(D.21) which we solve numerically in the end. Hence these discontinuities are numerical artifacts and can be neglected.