Supersymmetry Breaking in a Large N Gauge Theory with Gravity Dual

We study phase structure of mass-deformed ABJM theory which is a three dimensional $\mathcal{N}=6$ superconformal theory deformed by mass parameters and has the gauge group $\text{U}(N)\times \text{U}(N)$ with Chern-Simons levels $(k,-k)$ which may have a gravity dual. We discuss that the mass deformed ABJM theory on $S^3$ breaks supersymmetry in a large-$N$ limit if the mass is larger than a critical value. To see some evidence for this conjecture, we compute the partition function exactly, and numerically by using the Monte Carlo Simulation for small $N$. We discover that the partition function has zeroes as a function of the mass deformation parameters if $N\ge k$, which supports the large-$N$ supersymmetry breaking. We also find a solution to the large-$N$ saddle point equations, where the free energy is consistent with the finite $N$ result.


Introduction
Spontaneous supersymmetry (SUSY) breaking in string/M-theory is one of the most important subjects and has been discussed intensively in the context of phenomenology and cosmology. The SUSY breaking in string/M-theory is the super-Higgs phenomenon in general since the theory has gauged SUSY rather than global one. In contrast, various situations in string/M-theory are expected to have holographic descriptions by SUSY quantum field theories (QFT), whose SUSY are global. Therefore it is interesting to discuss SUSY breaking in QFT with a gravity dual which is typically large-N and strongly coupled. As far as we know, the only explicit examples of such problem are the models discussed in [1][2][3][4]

JHEP03(2019)159
in which the SUSY is kinematical, i.e. the SUSY algebra does not includes the Hamiltonian. 1 Main reason for the existence of the very few examples is that it is technically hard since we typically need non-perturbative analysis in this type of problem.
In this paper we study the SUSY breaking problem in so-called massive ABJM theory [7,8] which is a three dimensional N = 6 superconformal theory known as ABJM theory [9] deformed by two mass parameters and has the gauge group U(N ) × U(N ) with Chern-Simons levels (k, −k). It is expected that the ABJM theory without the masses is the low-energy effective theory of N coincident M2-branes and dual to the M-theory on AdS 4 × S 7 /Z k . The mass deformation in this situation corresponds to the introduction of the background flux. Then the massive ABJM theory in the large-N limit is holographically dual to M-theory on asymptotic AdS 4 geometry [10,11].
We argue that the mass deformed ABJM theory on S 3 breaks supersymmetry in the large-N limit with k fixed if the mass is larger than a critical value. We can address this because the partition function of the theory on S 3 can be exactly computed by the localization technique [12][13][14]. Note that the localization technique can be applied to our theory with finite N on S 3 as in Witten index even if it will break the supersymmetry spontaneously in the large-N or large volume limit. The arguments are based on the existences of the zeroes of the partition function which will be related to the SUSY breaking and the phase transition at the critical mass which is expected from the large-N saddle point solution found in [15]. A summary of our arguments for the SUSY breaking of the theory will be explained in section 3. In the previous work [15] a part of the authors analyzed this theory with the two equal mass parameters, which enjoys the N = 6 supersymmetry. We studied the partition function in the M-theory limit (N → ∞ with k kept finite) by using the saddle point approximation, and found that the saddle point solution which gives the free energy F ∼ N 3/2 disappears as we increase the mass deformation parameter to some critical value. Although this would suggest that a phase transition occurs at that point, the whole phase structure is still unclear.
The rest of this paper is organized as follows. After introducing the mass deformed ABJM theory with two mass parameters ζ 1 , ζ 2 in the next section, in section 3 we argue that this theory with ζ 1 = ζ 2 breaks the supersymmetry for ζ 1 /k > 1/4. In the same section we also argue that the supersymmetry breaking does not occur if ζ 1 = 0 or ζ 2 = 0. In section 4 we first study the latter cases with ζ 1 = 0 in detail, and indeed find the large-N free energy obeys N 3/2 -law for any value of ζ 2 . Then, in section 5 we consider the case with both ζ 1 and ζ 2 are non-zero. In this case the partition function Z(k, N, ζ 1 , ζ 2 ) for N ≥ k can have zeroes at some finite ζ 1 , ζ 2 , this was explicitly shown for N = 2. We provide positive evidence for the existence of zeroes from the Monte Carlo simulation. We also argue the physical interpretation for the zeroes, and estimate how the partition function behaves in the large-N limit. Our proposal on the phase structure in the large-N limit is summarized in figure 1. We expect that the partition function vanishes when ζ 1 ζ 2 k 2 = 1 16 and the theory is in the SUSY breaking phase for ζ 1 ζ 2 k 2 ≥ 1 16 . In section 6 we summarize our analysis and propose future directions. We summarize technical details in appendices. In appendix A, starting from the localization formula (2.5) we rewrite the partition function to a simpler matrix model (2.6), which we use in the subsequent sections. In particular, in the new expression the integration is absolutely convergent, hence we can evaluate the partition function numerically by applying the Monte Carlo method. In appendix B we display the exact computation of the partition function with ζ 1 = 0 with k and N being small integers, and summarize the results in appendix C. As N increases these results match with the saddle point approximation in section 4.1, which support the validity of the saddle point approximation. We also compare the exact partition function with the partition function of the linear quiver theory with single hypermultiplet obtained in [16].

JHEP03(2019)159
Lastly let us clarify a relation between our results and the previous results for the Rcharge deformations which were obtained by a saddle point method in the large-N limit [19] and a semi-classical expansion of Fermi gas formalism supposed to be correct up to nonperturbative effects of large-N expansion [34]. It is known that the theory with the Rcharge deformation is formally identical to the mass deformation with pure imaginary mass parameters ζ 1 , ζ 2 ∈ iR at the level of the localization formula (2.5). Therefore one might wonder whether the results for the real masses were simply obtained by an analytic continuation of the ones for the R-charge deformations or equivalently imaginary masses. Our claim is that this is true for small ζ 1 , ζ 2 otherwise this is not. In this paper, we have obtained the exact expressions 2 for Z(N, k, 0, ζ 2 ) and numerical results for Z(N, k, ζ 1 , ζ 2 ) with finite N ≥ 3, and the large-N expressions. The finite N results are completely new since it has not been computed even for the R-charge deformation case. The large-N expressions for small ζ 1 , ζ 2 is formally the same as the simple analytic continuation of the one obtained for the R-charge deformation. In this regard, our result for sufficiently small ζ 1 , ζ 2 is merely an analytic continuation of the known results although it is nontrivial when the analytic continuation works. However, the partition function behaves in a completely different way as |ζ 1 |, |ζ 2 | increases. In the case of R-charge deformation the partition JHEP03(2019)159 function diverges 3 when one of ζ 1 , ζ 2 reaches ζ i = ± ik 4 . On the other hand, the partition function never diverges for real masses, though it changes the behavior as ζ 1 , ζ 2 reaches some critical values. The shape of the phase boundary ζ 1 ζ 2 k 2 = 1 16 which is suggested from our analysis is also different from that for the R-charge deformation. In particular if ζ 1 = 0 the matrix model does not encounter any phase transition regardless of the value of ζ 2 . Although the understanding is not complete, we suspect this phase boundary particular for the real masses is related to the convergence property of the large-N instanton effects which is completely different between ζ 1 , ζ 2 ∈ iR and ζ 1 , ζ 2 ∈ R. We will elaborate this point in section 5.3.

Review on mass deformed ABJM theory
In this section we review some basic facts on the mass-deformed ABJM theory on S 3 . The field content of the ABJM theory consists of, in the 3d N = 2 SUSY notation, a U(N ) k vector multiplet V = (A µ , σ, χ, D), a U(N ) −k vector multiplet V = ( A µ , σ, χ, D), two chiral multiplets Z α = (A α , φ α , F α ) in ( ,¯ ) representation under U(N ) k × U(N ) −k and two chiral multiplets Wα = (Bα, ψα, Gα) in (¯ , ) representation. 4 Here the vector multiplets obey the Chen-Simons action with level ±k, while the action for the chiral multiplets consists of the superpotential together with the following minimal coupling to the vector multiplets We can introduce a mass by turning on a background vector multiplet V (bgd) = (A (bgd) µ , σ (bgd) , χ (bgd) , D (bgd) ) of a global symmetry in the following supersymmetric configuration 5 [17] A (bgd) where we have set the radius of S 3 to r S 3 = 1. Here we turn on the background multiplets of the flavor symmetries U(1) 1 × U(1) 2 × U(1) 3 commuting with the N = 2 supersymmetry under which the chiral multiplets are charged 6 as in table 1. The background gauge fields so that ζ 1 , ζ 2 are interpreted as the mass parameters for the chiral multiplets as m 1 = 2ζ 1 /k for Z 1 , W2 and m 2 = 2ζ 2 /k for Z 2 , W1.

JHEP03(2019)159
which we derive in appendix A. For k = 1 and ζ 1 = ζ 2 = 0, this latter expression coincides with the partition function of the N = 8 U(N ) Yang-Mills theory coupled with a fundamental chiral multiplet, which is dual to the ABJM theory under the SL(2, Z) transformation in the type IIB brane setup. Because of this reason we simply refer to (2.6) as the S-dual representation even for general (k, ζ 1 , ζ 2 ). Note that the integration in the S-dual representation (2.6) is absolutely convergent in contrast to the representation (2.5) where the convergence is achieved by the rapidly oscillating factors. Because of this fact, it is much easier to apply the Monte Carlo simulation of the partition function to (2.6) than (2.5). With the help of the Monte Carlo simulation of (2.6) we will observe a novel behavior of the partition function: the partition function vanishes at some finite values of ζ 1 , ζ 2 , which was not encountered in the undeformed case or the case of the R-charge deformation (ζ 1 , ζ 2 ∈ iR).

Evidence for SUSY breaking
In this section, we discuss why we expect the SUSY breaking of the mass deformed ABJM theory on S 3 in the large-N limit at some finite (ζ 1 , ζ 2 ) and explain our criterion for the SUSY breaking which we will examine in the following sections.
First, in the case of ζ 1 = ζ 2 = ζ, there is a large-N saddle point solution for the original matrix model (2.5) which exist only for 0 ≤ ζ k < 1 4 [15]. This solution becomes the saddle point solution of the massless ABJM theory [20] in the ζ → 0 limit and gives the N 3/2 -law of the free energy: However, this saddle point solution becomes singular in the ζ k → 1 4 limit. There is another large-N solution for any value of ζ. The free energy of this solution is proportional to N 2 and this solution may correspond to a confinement vacuum. 8 Although it is still possible that there exist other solutions, we could not find any solution other than the two solutions with some numerical methods. Therefore we regard these results as strongly indications of a phase transition at ζ k → 1 4 . We expect that this phase transition comes from SUSY breaking as follows. Let us take the mass very large, i.e. ζ k 1, then, at least naively 9 , the hypermultiplets become heavy and decouple from the vector multiplets. The remaining N = 2 SUSY pure Chern-Simons theory will spontaneously break SUSY as shown in [21][22][23] and become the confinement phase in the large-N limit. This expectation is consistent with the above large-N solutions. 8 This statement is not precise because the Chern-Simons interaction remains and theory may be in a gapped phase. Nevertheless we will call the confinement phase also for such case. Note that here we take k/N → 0 limit, thus the Chern-Simons interaction will be ignored for the leading order in the large-N limit and the Yang-Mills term always induced by the renormalization flow. We also note that the N = 2 SUSY pure Yang-Mills theory does not have SUSY vacua. 9 This is naive since we have to integrate over Coulomb branch in the theory on S 3 as seen from (2.5).
Therefore this naive argument implicitly assumes that the integral (2.5) in the large-N limit is dominated by a point in (λi,λi) satisfying |λi −λj −4πζ1,2/k| 1 so that all the hypermultiplets are infinitely massive. See also footnote 30 for related comments.

JHEP03(2019)159
However, for the mass deformed ABJM theory, the SUSY index was computed to be nonzero in [11,24]. In this theory, there are infinitely many discrete classical vacua which are characterized by the fuzzy S 3 solutions given in [8,25], which represent M5-branes. Although the contribution to the index for the trivial vacuum, where all the scalar fields are zero, vanishes as in the pure SUSY Chern-Simons theory, other vacua give the non-zero contributions to the index if there are no coincident M5-branes. One might think that this result contradicted the above argument of the SUSY breaking. However, this result is for the theory on T 3 , not on S 3 . For the N = 2 SUSY theory on S 3 , there are mass terms for the hypermultiplets proportional to the curvature of S 3 . The mass term will lift all of the vacua except the trivial vacuum at the origin classically. 10 Thus, the result of [11,24] on the SUSY index does not exclude the possibility that the mass deformed ABJM theory on S 3 has the SUSY breaking phase. 11 Here note that we do not take the large volume limit.

Criterion for SUSY breaking
By now, we have not defined what the spontaneous SUSY breaking on S 3 is. Usually, the SUSY breaking means that there is no states with zero energy in the theory. For S 3 , we can not define states with an appropriate Hamiltonian and time, thus it is difficult to use this definition. Instead of this definition for the SUSY breaking, the spontaneous breaking of a symmetryQ can be defined as ∃Ô s.t. 0|[Q,Ô]|0 = 0. In the path-integral formalism, this corresponds to where the condensation is the order parameter. Note that this correspondence is valid for the theory with sufficient number of non-compact space directions, in which the notion of vacuum is meaningful, otherwise, QO corresponds to Tr[Q,Ô], 12 not to 0|[Q, O]|0 . Because Q is a symmetry generator which behaves well, we expect that QO = 0 (Tr[Q, O] = 0) is trivial identity due to the invariance of path integral measure (cyclic invariance of Tr). 13 For example, for SUSY quantum mechanics case, the invariance of the Witten index means Tr[(−1)FQ,Q † ] ∼ Tr(−1)FĤ = 0. Thus, the definition of (3.2) is meaningful for the theory with some space with enough number of non-compact directions. Since S 3 is compact, we need to take the large volume limit or large-N limit which can effectively gives extra dimension. If this happens, then there should exist a massless Goldstone fermion in 10 We expect that the energy of the possible metastable SUSY breaking vacuum is proportional to ζ and the free energy will be proportional to ζr S 3 . The extra contribution by the curvature induced mass term to the free energy for the fuzzy sphere solutions will be also proportional to ζr S 3 because the size of the fuzzy sphere grows as ζ grows. Of course, this is not valid except the weak coupling limit and the phase of the theory can be non-trivial. 11 Here, we assume that the theory is regarded as a deformation of the ABJM theory on S 3 for a small ζ/k case. For a enough large ζ/k case, we think that the curvature effect of S 3 is almost negligible, but still remains. This picture will lead the SUSY breaking scenario explained here. 12 For the SUSY, it corresponds to Tr(−1)F {Q,Ô} = Tr[(−1)FQ,Ô]. 13 In the case of Q =SUSY, QO is such as F -term and D-term. Unfortunately we cannot compute F or D by using the supersymmetry localization. We can compute F and D , but they are trivially zero. This is consistent with the fact that there is no SUSY breaking for the theory on S 3 with N finite.

JHEP03(2019)159
the theory which makes the (SUSY) partition function Z vanished. Note that in the localization computation the Z = 0 is not technically realized by Goldstone fermions because the matrix model (2.5) is bosonic. Instead, Z = 0 is realized by the phase cancellation in the integral over Coulomb branch. This is not inconsistent with the above interpretation in terms of Goldstino because the localization technique uses the large deformations of the theory keeping some physical quantities invariant which are Yang-Mills coupling and matter couplings for our case.
Instead of the large volume limit, we take a large-N limit in which the SUSY breaking is meaningful. Thus, we need a criterion of the SUSY breaking in the large-N limit from a finite N result. For the theory in which we can define the Witten index, the vanishing of it, i.e. Z = 0, is the necessary condition for the SUSY breaking for finite volume. For the other theories also, we expect that the massless Goldstone fermion makes Z = 0. Indeed, for a superconformal theory on S 3 , the theory can break the SUSY even for finite N if Z = 0 because the radius of S 3 is not physical. Such theories were discussed in [5,[26][27][28][29]31]. A well-studied example is 3d N = 2 SU(N ) k Chern-Simons theory which has several evidence of SUSY breaking for k < N including vanishing of sphere partition function. In addition, the theory in this regime has vanishing Witten index [21], SUSY breaking brane construction [22,23] and supergravity dual solution [5]. There are also similar results for the case with fundamental matters [26][27][28]30] and ABJ theory [31][32][33]43].
For our case, the theory is not conformal, but we take the large-N limit. Thus, we regard Z = 0 as a criterion of the SUSY breaking. 14,15 It is worth to note that Z = 0 does not necessarily mean SUSY breaking as in Witten index. However, for our case, interpreting Z = 0 as SUSY breaking is the most natural possibility because the mass deformed ABJM theory on S 3 will be smoothly connected to the pure SUSY CS theory in the large mass limit whose SUSY is broken for k < N In the following sections, we will give further supporting arguments for the above picture of the SUSY breaking phase for the mass deformed ABJM theory on S 3 using the S-dual representation of the matrix model. Here we will summarize these argument for the SUSY breaking shortly. The S-dual representation of the matrix model (for k = 1) is obtained from the U(N ) Yang-Mills theory with an adjoint and fundamental matter fields where ζ 1 and ζ 2 corresponds to the FI term and the mass for the adjoint matter, respectively. Because of the FI term (and the mass term), this theory will break the SUSY at the origin of the Coulomb branch moduli space which will be favored by the mass terms induced from the curvature of S 3 . This picture will be right for a generic large value of ζ 1 and ζ 2 . However, for ζ 1 = 0, the FI term vanishes and the SUSY will not break. 16 For 14 In the gravity dual, the SUSY is gauged and the theory is described by a supergravity. In the supergravity, there are massless fermions, however, there are no zero modes around the SUSY vacuum which is an asymptotic AdS4 background. In a SUSY breaking vacuum, some fermions near the boundary have zero modes. 15 From the analogy to the case with bosonic zero mode, an appropriate analysis would be to add an explicit-susy-breaking deformation to kill the zero mode and see what happens in the limit of zero deformation. In this approach, however, we cannot use the result of the localization. 16 The mass deformed ABJM also will not break the SUSY for this case because a half of the hypermultiplets remain massless and do not decouple.

JHEP03(2019)159
this case, as we will see later, we can construct a large-N solution for any value of ζ 2 , thus there are no critical mass for this case. This is consistent with the above picture.
In order to investigate further, we will compute the partition function Z for finite N exactly and numerically using the Monte Carlo method for various points of (ζ 1 , ζ 2 ). We expect that some values of N for which we computed Z are not very large, but enough large for the large-N expansion. Indeed, the computed values of Z are consistent with the large-N solutions for ζ 1 = ζ 2 < k/4 and ζ 1 = 0. These actual computations of Z for finite N shows that as increasing ζ i , Z is decreasing and oscillating, thus Z = 0 for some values of ζ i . We expect that this zero corresponds to the SUSY breaking in the large-N limit. Furthermore, if we increase N with other parameters fixed, the smallest value of ζ i which gives Z = 0 decreasingly approaches to the critical point of the large-N solution. Therefore, the extrapolation of this to the large-N limit may be consistent with the SUSY breaking picture above. 4 The case with one massless hypermultiplet (ζ 1 = 0) In this section we consider the case with ζ 1 = 0. In this case we find a solution to the saddle point equation for the partition function in the S-dual representation (2.6). We can also compute the exact values of the partition function for finite (N, k) by a slight generalization [34,35] of the technique used in the ABJM theory [36,37]. We will see a good agreement of the both results.

Saddle point analysis in the large-N limit
In this subsection we compute the partition function in the large-N limit N → ∞, with fixed (k, ζ 2 ). (4.1) In this limit, we can evaluate the partition function by the saddle point method. To perform the saddle point analysis, we first introduce the effective action S eff by We rearrange the eigenvalues x i such that x i+1 ≥ x i by the permutation symmetry and regard x i as a function of s = i/N − 1/2, which becomes the continuous variable in the large-N limit:

JHEP03(2019)159
Then the summations over i are replaced by the integral over s We look for saddle point solutions by the approach taken in [20] which has been used to derive O(N 3 2 ) behaviors of free energies, rather than the traditional approach often applied for matrix models in the planar limit. 17 This is achieved by taking the following ansatz with an O(1) real 18 function z(s), and perform large-N expansion of S eff (x) to simplify the saddle point equation. It is easy to write down the leading part for the first and second terms in (4.3): We can also expand the third and fourth terms in (4.3) respectively by using the techniques of [15] to see that the leading part of S eff in the large-N limit is proportional to N 3 2 . First we rewrite these terms as (4.9) 17 The traditional approach was taken in [38,39] for ζ1 = 0 = ζ2 identifying 't Hooft coupling with N/N f where N f denotes an additional power put on the cosh (For our case, N f = 1). 18 In our actual analysis, we have looked for solutions with complex z(s) under the ansatz (4.6) but we have found only a real solution as a result. Because of this, we take z(s) to be real for simplicity of explanations in the main text. Precisely speaking, we should first take the variation δS eff δz(s) with z(s) ∈ C before assuming z(s) ∈ R. This induces a new constraint δS eff δIm(z(s)) = 0 in addition to (4.15) and (4.16); nevertheless the final result (4.18) remains the same. 19 Note that odd functions of z − z do not contribute.

JHEP03(2019)159
where z is the abbreviation for z(s ). Note that the O(N 5/2 ) terms, which are the first terms in (4.8) and (4.9), are canceled and only the second terms remain. Here we use the following formula in the large-N limit: 20 where y(s) = √ N v(s) + w(s) and s 0 is the zero of v(s). These formulas are obtained by changing the integration variable and reflected with the fact that the contribution to the integral in l.h.s. of (4.10) and (4.11) is coming from only s ∼ s 0 region in the large-N limit. Using these formulas the second terms in (4.8) and (4.9) can be evaluated as Putting the above computations together, we find the following large-N expansion for the effective action (4.14) The overall scaling N 3/2 in (4.13) implies that the integration (4.2) is dominated in the large-N limit by the saddle point configuration satisfying the following equation of motion together with the boundary condition (4.16) 20 The condition that this evaluation is valid is following [15]: which is satisfied in our case.

JHEP03(2019)159
First let us consider the case for ζ 1 = 0. First of all the equation of motion (4.15) has the following two local solutions depending on sgn(z) where s b and z b are the integration constants, and we have excluded the other two solutions by the conditionż(s) ≥ 0. The bulk solution would be obtained by connecting these solutions appropriately and determining the integration constants so that z(s) satisfies the boundary conditionż(s) = ±∞ (4.16) at every point where z(s) orż(s) is discontinuous. Notice that both of z (±) (s) satisfiesż (±) (s) = ∞ at only a single point s = s b . Therefore, if we split the support −1/2 < s < 1/2 into segments by the points of discontinuity, z(s) on each segment must be given as a smooth junction of z (−) (s) and z (+) (s). Since z (+) (s) cannot be followed by z (−) (s) due to the assumption that z(s) is monotonically increasing, we conclude that the solution is given by a single junction of z (−) (s) with . In summary we obtain the following unique solution as the saddle point configuration: In the language of the eigenvalue density, this solution corresponds to Substituting this solution to (4.13), we find that the partition function in the large-N limit is given as (4.20) For ζ 1 = 0 we could not solve the saddle point equation with the ansatz we used here because the solution can not satisfy the boundary condition (4.16) due to the existence of the imaginary term 2iζ 1 k in (4.15). However, the partition function with ζ 1 = 0, ζ 2 = 0 and that with ζ 1 = 0, ζ 2 = 0 is the same because the partition function is invariant under exchanging ζ 1 and ζ 2 . This fact suggests that even when ζ 2 = 0, ζ 1 = 0, there exists the solution of the saddle point equation in large-N limit and the free energy can be evaluated by the saddle point approximation.

Exact partition function for finite (N, k)
Next we compute the partition function for some finite (N, k) by the technique used in [34]. We start with the partition function written in the Fermi gas formalism If we consider the generating function of the partition function or equivalently the grand partition function ∞ N =0 z N Z(N ), we can show that it is written as the following Fredholm determinant Comparing the coefficient of z N on the both sides, we find that the partition function Z(N ) is determined by Tr ρ n with n ≤ N , as We can compute Tr ρ n by completely the same way as that in the case of R-charge deformation [34]. First we notice that the matrix element x| ρ|y has the following structure For α = 1, this form is in the range of application of Tracy-Widom's lemma [36] which has been very powerful tool to systematically compute Trρ n in various M2-brane theories without masses [37,[40][41][42][43][44][45][46][47]. We can easily extend it to general α as follows. The structure (4.25) can be expressed as a quasi-commutation relation for ρ For a later convenience we have symmetrized the density matrix by another similarity transformation from (A.9).

JHEP03(2019)159
where |p is momentum eigenstate satisfying (4.28) This relation can be generalized straightforwardly for ρ n as This implies that we can compute the matrix element of ρ n from two sets of functions φ (x) and ψ (x) as We can show that the function φ (x) satisfies the following recursion relation 5 General deformation with ζ 1 , ζ 2 = 0 In this section we consider the case for ζ 1 , ζ 2 = 0. Note that this may affect the sign of the partition function because the integrand of (2.6) for ζ 1 = 0 has the oscillation factor e 2iζ 1 k i x i in contrast to the ζ 1 = 0 case, where the integrand was positive semi-definite. Therefore the partition function may be negative or zero depending on the parameters (N, k, ζ 1 , ζ 2 ). For large ζ 1 , ζ 2 , we can easily see that this actually happens as follows. In this limit, the hypermultiplets become very massive and integrating them out leads us to the N = 2 SUSY U(N ) k × U(N ) −k pure Chern-Simons theory schematically. 22 It is known 22 More precisely, integrating out the matter fields induces level shifts of all the possible mixed CS terms which are among the gauge symmetry U(N ) × U(N ), flavor symmetry U(1)1 × U(1)2 and U(1)R symmetry. In the case of the massive ABJM theory, most of the shifts are canceled and we have only contributions from the gauge-U(1)R and flavor U(1)R CS terms but these terms do not affect the zeroes of the partition function. Here the integration of the matter fields is assumed to be at the origin of the Coulomb moduli space. Thus, the decoupling of the matter fields in the large mass limit is possible.

JHEP03(2019)159
that the sphere partition function of the pure Chern-Simons theory vanishes for k < N . In this section, we will see that the zeroes appear also for finite (ζ 1 , ζ 2 ).
In this case we could not find a solution to the saddle point equation. The technique for small integers k, N in section 4.2 is not applicable either. Nevertheless we can evaluate the partition function exactly for N = 1, 2, which suggest the partition function has zeroes as a function of ζ 1 , ζ 2 only for (N, k) = (2, 1), (2,2). We argue a possible interpretation for this zeroes. We further conjecture the zeroes for general k, N , and provide positive evidence from the numerical computation of the partition function for N ≥ 3.

Exact expression for N = 1, 2
In this subsection we review the exact results for N = 1, 2 obtained in [49]. 23 The relation (4.23) between the partition function and Tr ρ n is correct also for general ζ 1 if we take ρ as For N = 1, the partition function is simply given by Z(1, k, ζ 1 , ζ 2 ) = Tr ρ, which can be exactly computed as For N = 2, we need to compute Tr ρ 2 , which is given by the following two dimensional integration After changing the integration variables to x ± = x ± y, we can easily perform the x +integration, which leads to For k ∈ Z + , this integral can be evaluated by considering an integral with the same integrand along a rectangular whose corners are x − = (−∞, ∞, ∞ + 2πik, −∞ + 2πik) [48,49], and we obtain 5) 23 The notation in [49] is related to ours by

N ≥ 3 from Monte Carlo Simulation
In this subsection we provide numerical evidence that the partition function has zeroes at finite (ζ 1 /k, ζ 2 /k) also for N ≥ 3. For this purpose, we apply (Markov chain) Monte Carlo method to the partition function in the S-dual representation (2.6):

Algorithm
First we explain our algorithm. There are two subtleties in applying the Monte Carlo method to our problem. The first subtlety, which will not be problematic as explained below, is that Monte Carlo simulation can directly calculate only "expectation values" or equivalently ratio of two functions rather than Z itself. The second one is that the Boltzmann weight e −S is not positive semi-definite for ζ 1 = 0 and hence cannot be regarded as probability. This problem appears in many contexts such as finite density QCD, real time systems and theories with CS terms. We take care of these subtleties as follows. Instead of Z itself, we consider the ratio

JHEP03(2019)159
where O(x) ζ 1 =0 denotes the expectation value of O(x) under the action S(N, k, ζ 1 = 0, ζ 2 ). Then we approximate the ratio by Hybrid Monte Carlo simulation 24 by taking samples generated with the probability ∼ e −S ζ 1 =0 . Note that studying only the ratio is sufficient for our purpose since Z(N, k, 0, ζ 2 ) is real positive and we are interested in the sign of the partition function. 25 Since we are taking samples of the oscillating function, whose oscillation is controlled by ζ 1 /k, we typically need more statistics for larger ζ 1 /k to obtain precise approximations. Note also that the S-dual representation (2.6) has much milder oscillation than the original matrix model (2.5). This is why we are using the S-dual representation as in [50].

Results
Now we present numerical results for the ratio Z MC (5.10), which has the same sign as the partition function Z(N, k, ζ 1 , ζ 2 ) itself. Figure 2 plots Z MC for (N, k, ζ 2 ) = (4, 1, 1) as a function of ζ 1 . The statistical errors are estimated by Jackknife method although they are practically almost invisible in the figures. The right panel of figure 2 is the zoomup of the left figure in the range ζ 1 ∈ [0.1, 0.2]. From the right figure, we easily see that the partition function takes negative values when ζ 1 = 0.110, 0.115, · · · , 1.14 even if we take into account the errors. Therefore there must be a zero of the partition function for ζ 1 ≤ 0.110 and the plot indicates that the zero is located at 0.105 < ζ 1 < 0.110. 24 This is so-called reweighting method. The application of a similar algorithm to a similar system is explained in appendix A of [50]. 25 Of course we can also compute Z(N, k, ζ1, ζ2) itself by combining ZMC(N, k, ζ1, ζ2) with Z(N, k, 0, ζ2)  Table 3. Bounds and estimate of first zeroes of the partition function for (k, ζ 2 ) = (1, 2).
We have found similar results for other values of (N, k, ζ 2 ) whose samples are shown in figure 3. These figures indicate that the partition function has the zeroes at finite ζ 1 /k for various (N, k, ζ 2 ). Note also that we sometimes encounter subtle cases. For example, in the case of (N, k, ζ 2 ) = (4, 2, 1) shown in the right-bottom of figure 3, the minimum is consistent with both positive and negative Z within the numerical errors. 26 We expect that this type of behavior appear when the partition function is positive semidefinite but has zeroes as in the case of (N, k) = (2, 2) whose analytic result is given in the second line of (5.7). For this type of cases, any numerical simulation with nonzero errors cannot establish existence of zeroes since numerical values at the zeroes must be consistent with all the possible signs of Z within errors. Therefore, for this type of cases, the best thing we can do by numerical simulation is to check existence of points consistent with Z = 0. For all values of (N ≥ 2, k, ζ 2 ) which we have analyzed, we have checked that there exists at least one value of ζ 1 consistent with Z = 0 within errors. For the cases where we have established existence of first zeroes of Z, we write down bounds on the zeroes in tables 2, 3 and 4 for fixed (k, ζ 2 ) (see table 5 for Z MC at the first negative peaks and their errors). We also estimate their precise locations by constructing 27 interpolating functions of all the data points of Z MC for fixed (N, k, ζ 2 ) and finding zeroes of the interpolating functions. We will discuss implications of these values in section 5.4.   Table 5. Z MC at the first negative peaks and their statistical errors for various (N, ζ 2 ) with k = 1.

Physical origins of the zeroes and Fermi gas formalism
In this subsection we discuss physical origins of the zeroes of the partition function. For this purpose, we apply Fermi gas formalism and identify which effects trigger the change of the sign of Z. Note that some techniques in the Fermi gas formalism are not available for ζ 1 = 0 since the Hamiltonian is not hermitian. However, there is a technique which is still available. This is a formal -expansion of Tr ρ n via Wigner transformation where Tr ρ n is expressed as a phase space integral of a function whose explicit representation can be obtained by acting differential operators on ρ(q, p). In this technique, it does not matter whether or not the Hamiltonian is hermitian since the problem is reduced to compute a perturbative series of the explicit two dimensional integral with respect to . Fortunately, this analysis has been already done in [34] for imaginary (ζ 1 , ζ 2 ) in the context of the R-charge deformation and hence we can obtain the -expansion simply by analytic continuation of the result in [34] up to a subtlety discussed below. 28 Once we find Tr ρ n approximated in this way, one can compute the grand potential J(µ) by the following Mellin-Barnes expression where Z(t) = Tr ρ t and the canonical partition function can be obtained from J(µ) by The -expansion of Z(n) takes the form

JHEP03(2019)159
where the second term denotes non-perturbative effects of the -expansion which we are ignoring. The work [34] computed the first four coefficients Z 0 , Z 2 , Z 4 and Z 6 which are explicitly written down in appendix A of [34]. For example, the leading order coefficient Z 0 is given by where we are keeping (ζ 1 /k, ζ 2 /k) fixed and The large-N behavior of Z(N ) can be easily derived by the large-µ behavior of J(µ) which has the following structure Several comments are in order. First, the -expansions for the coefficients C and B are terminated at leading and sub-leading orders respectively: The coefficient A receives all order corrections and it has been conjectured in [34] that the exact answer for A is given by where [44,50] A ABJM (k) = 2ζ(3) . (5.20) If the approximation by J pert (µ) is reliable, then the canonical partition function is approximated by The large-N limit of this formula exhibits the N 3/2 -law: 29 22) which agrees with (4.20) for ζ 1 = 0 and the result of [15] for ζ 1 = ζ 2 = ζ. 29 In the large-N limit, the µ-integral is dominated by µ = N −B C . Therefore the non-perturbative effects in (5.16)

JHEP03(2019)159
The second term in (5.16) is non-perturbative corrections of the large-µ expansion whose exponents can be explicitly derived by the -expansion (5.13). These corrections for the massless case have been identified with membrane instanton effects whose type IIA picture is D2-branes wrapping (warped) RP 3 in AdS 4 × CP 3 [51]. The third term in (5.16) takes over the non-perturbative correction of the -expansion in (5.13) whose exponents cannot be determined by the above arguments. It has been conjectured in [34] that the exponent for imaginary (ζ 1 , ζ 2 ) is given by O(e − 4µ k(1±4iζ 1 /k)(1±4iζ 2 /k) ) (double signs correspond). These corrections for the massless case have been identified with worldsheet instanton effects coming from fundamental strings wrapping CP 1 [52].
Let us estimate when we can trust the approximation by the perturbative part J pert (µ) in the large-µ expansion (5.17), or equivalently when the canonical partition function is well approximated by (5.21). We easily see that the second term in (5.17) is exponentially suppressed for any (ζ 1 , ζ 2 ) and therefore we can ignore the second term in the large-N limit. Then let us focus on the third term which comes from non-perturbative effects of the -expansion. We have not estimated the exponent of the third term for real (ζ 1 , ζ 2 ) precisely. However, the exponent for real (ζ 1 , ζ 2 ) should be the same as the naive analytic continuation of the one for imaginary (ζ 1 , ζ 2 ) in the domain where the partition function is holomorphic with respect to (ζ 1 , ζ 2 ). Therefore, if we assume the above conjecture on the exponent for imaginary (ζ 1 , ζ 2 ) in [34], then we should have the following correction in (5.17) for real (ζ 1 , ζ 2 ): where the double signs are in the same order. Note that this correction is no longer exponentially suppressed for and we cannot approximate the grand potential J(µ) by J pert (µ) in this region. This also implies that the holomorphy of the partition function with respect to (ζ 1 , ζ 2 ) is broken at ζ 1 ζ 2 k 2 = 1 16 because if we start with imaginary (ζ 1 , ζ 2 ), then the naive analytic continuation to real (ζ 1 , ζ 2 ) does not commute with the large-N limit. Namely, if we take the large-N limit first, then the free energy behaves as ∼ N 3/2 and its continuation to real (ζ 1 , ζ 2 ) is also described by the same formula for any (ζ 1 , ζ 2 ) which is very likely different from the large-N limit after the continuation in the domain ζ 1 ζ 2 k 2 ≥ 1 16 . The above estimate is consistent with our numerical results obtained in section 5.2. In figure 4 we compare the ratio Z MC (N, k, ζ 1 , ζ 2 ) = Z(N,k,ζ 1 ,ζ 2 ) Z(N,k,0,ζ 2 ) obtained by the Monte Carlo simulation with the one computed by the approximation (5.21) for some cases with 16 where we expect (5.21) to be good approximation. The plots show that our numerical results agree with the Airy function formula (5.21) and exhibit the N 3/2 -law. Although we explicitly present only the four cases, we have observed similar behaviors for various other values of (k, ζ 1 , ζ 2 ) satisfying ζ 1 ζ 2 k 2 < 1 16 . Figure 5 shows similar plots for 16 where we expect that we cannot trust (5.21) due to the correction (5.23). In contrast to figure 4, we easily see that the numerical results do not agree with (5.21) and JHEP03(2019)159  no longer exhibit the N 3/2 -law. We have also found similar behaviors for various other values of (k, ζ 1 , ζ 2 ) with ζ 1 ζ 2 k 2 ≥ 1 16 . Thus our numerical results support our expectation that the approximation by (5.21) is valid for ζ 1 ζ 2 k 2 < 1 16 .

Conjecture on phase structure in the large-N limit
We discuss the phase structure of the mass deformed ABJM theory in the large-N limit. Let us first recall the results obtained so far by the various analyzes: • In the case of ζ 1 = ζ 2 = ζ, the partition function in the representation (2.5) has been analyzed by the saddle point method in [15] as reviewed in section 3. The saddle point configuration in [15] realizes the O(N 3/2 ) free energy and becomes singular at ζ/k = 1/4.

JHEP03(2019)159
• In section 4.1, we have constructed the saddle point solution for ζ 1 = 0 in the Sdual representation, which gives the O(N 3/2 ) free energy (4.20). This behavior is consistent with the exact results for finite N obtained in section 4.2. Note also that we can obtain the result for ζ 1 = 0, ζ 2 = 0 by the replacement ζ 2 → ζ 1 in (4.20) since the partition function is symmetric under ζ 1 ↔ ζ 2 .
• In section 5.1, we have written down the exact results for N = 1, 2 and arbitrary (k, ζ 1 , ζ 2 ) obtained in [49]. It has turned out that the partition function for N = 2 has the zeroes at finite (ζ 1 , ζ 2 ) while the one for N = 1 does not.
• In section 5.2, we have performed the Monte Carlo simulation for higher N . We have observed that the partition function has the zeroes at finite (ζ 1 , ζ 2 ) given (N, k).
The bounds and estimates on the zeroes given in tables 2, 3 and 4, imply that the first zeroes do not increase by N . It is natural to expect that the partition function becomes zero at some finite (ζ 1 , ζ 2 ) also in the large-N limit.
• In section 5.3, we have argued when one can trust the approximation in terms of the perturbative grand potential (5.17) in the Fermi gas formalism, which gives the O(N 3/2 ) free energy in the large-N limit. We have found that the approximation is reliable for ζ 1 ζ 2 k 2 < 1 16 while in the other regime ζ 1 ζ 2 k 2 ≥ 1 16 , the expected nonperturbative effects (5.23) of the -expansion are no longer exponentially suppressed. Note that for ζ 1 = ζ 2 = ζ, the approximation starts to be invalid at ζ/k = 1/4, which is the same as the condition that the saddle point in [15] becomes singular.
Based on the above results, we propose the following scenario (see figure 1 for schematic picture): (i) For small (ζ 1 , ζ 2 ), the large-N free energy behaves as F ∼ N 3/2 whose explicit form is given by (5.22). We expect that this formula is valid for ζ 1 ζ 2 k 2 < 1 16 , and becomes invalid for ζ 1 ζ 2 k 2 ≥ 1 16 .
(ii) The partition function vanishes at some finite values of (ζ 1 , ζ 2 ). We expect that this occurs at the boundary of the validity of (5.22): ζ 1 ζ 2 k 2 = 1 16 . We interpret this as the SUSY breaking at this point.
We already have strong evidence of the first point by the saddle point analysis for ζ 1 = ζ 2 in [15] and Fermi gas analysis in section 5.3. Now we provide further evidence for the second point. From tables 2, 3 and 4, we observe that the locations of the first zeroes decrease slowly as N increases. Therefore it is plausible that the first zeroes in the large-N limit are at some finite values of (ζ 1 , ζ 2 ). It would be nontrivial whether or not the first zeroes in the large-N limit coincide with our expectation ζ 1 ζ 2 k 2 = 1 16 . We perform consistency checks of this by fitting analysis of our numerical data. Note that the fitting analysis in the current situation is subtle in the following two reasons. First we do not know asymptotic behaviors of the first zeroes for large-N . In other words, we do not know what appropriate fitting functions are a priori. Second, we do not have sufficient data of the first zeroes since it is only for the five values of N (N = 2, 3,5,7,9). Nevertheless JHEP03(2019)159  the fitting analysis provides quite nontrivial consistency checks as we will see soon. We have constructed fitting functions for the first zeroes ζ 1 (N, k, ζ 2 ) with fixed (k, ζ 2 ) and varied N . As a conclusion, we have found that when we find a fitting function nicely interpolating numerical data, asymptotic value of the fitting function at N → ∞ agrees with our expectation ζ 1 = k 2 16ζ 2 . In figure 6, we plot the estimated zeroes for (k, ζ 2 ) = (1, 2) given in table 3 against 1/ √ N . We construct a fitting function for this data by the ansatz where a(k, ζ 2 ) corresponds to the zero in the large-N limit with respect to ζ 1 . We easily see that the fitting function nicely interpolates the data points. Therefore it is natural to compare the asymptotic value of the fitting function at large-N with our expectation. As a result, we have found a(k = 1, ζ 2 = 2) = 0.0309508 +0.0003776 −0.0003646 , which includes our expectation on the zero at large-N : ζ 1 = k 2 16ζ 2 | (k,ζ 2 )=(1,2) = 0.03125. This strongly supports our expectation on the large-N phase structure. Figure 7 shows results by the same fitting function (5.25) for the other values of ζ 2 . It is clear that the fitting functions (5.25) for (k, ζ 2 ) = (1, 1) and (1,5) do not nicely interpolate the data as much as for the case of (k, ζ 2 ) = (1, 2) although the fitting for (k, ζ 2 ) = (1, 5) is better than the one for (k, ζ 2 ) = (1, 1). From the fitting functions, we have found a(k = 1, ζ 2 = 1) = 0.0369738 +0.00242097 −0.00251291 and a(k = 1, ζ 2 = 5) = 0.0154573 +0.00012051 −0.000114297 which do not include our expectation ζ 1 = k 2 16ζ 2 although the result for (k, ζ 2 ) = (1, 5) is not far from the expectation. We interpret that this does not mean invalidity of our expectation since the fitting ansatz (5.25) does not exhibit very nice interpolations for (k, ζ 2 ) = (1, 1) and (1,5), and we need another fitting functions or more data. In figure 8 we also present similar plots to figure 7 by the different fitting function a(k, ζ 2 ) + b(k,ζ 2 ) √ N for (k, ζ 2 ) = (1, 1) and (1,5). Again the fitting functions do not interpolate the data very nicely and have different intercepts from the ones obtained by the ansatz (5.25) although the intercept for (k, ζ 2 ) = (1, 1) includes our expectation: a(k = 1, ζ 2 = 1) = 0.0628649 +0.000690078 −0.000682495 . To summarize we need to find more appropriate fitting functions or data points for larger N in order to further check our expectation except for (k, ζ 2 ) = (1, 2). We leave this for future work.
Moreover, the correlation between the supersymmetry breaking and the singularity in the saddle point approximation was argued for the pure Chern-Simons theory [29]. Hence it would be more than just a minimal scenario for our theory to relate the singularity in the saddle point approximation with the supersymmetry breaking. It would be interesting to test this conjecture by studying the partition function for larger N in future.

Discussion
In this paper we have studied the mass deformed ABJM theory on the three sphere. Based on the argument in section 3, we expect that this theory exhibits a spontaneous supersymmetry breaking in large-N limit at ζ 1 = ζ 2 = k/4. To gain an evidence for this conjecture we have analyzed the partition function of the mass deformed ABJM theory for finite k

JHEP03(2019)159
and N by using the Monte Carlo simulation. As a result we have found that the partition function vanishes at some finite values of ζ 1 , ζ 2 . The numerical results also indicate that the zeroes exist for general N , and that the locus of the first zero stays finite as N increases. These observations are consistent with the expectation in the end of section 3 from the large-N supersymmetry breaking. We also found a saddle point solution for a new slice ζ 1 = 0 by using the S-dual representation of the matrix model. In contrast to the situation for ζ 1 = ζ 2 this solution exists for an arbitrary value of ζ 2 . This is again consistent with our argument in section 3 where the supersymmetry breaking does not occurr in this case. Our result would shed new light to the phase structure of the mass deformed ABJM theory in the M-theory limit, which was unclear in the previous works [15,53].
Another direction is to improve the algorithm of the numerical simulation. In this paper, we have treated the oscillation factor of (2.6) in the quite naive way where we just regard the factor as the observable in the system with ζ 1 = 0. In this approach, we need much more statistics than simulations without oscillating factors so that the simulation at large-N becomes harder. It is nice if one can find more appropriate algorithm such as complex Langevin method and Lefschetz thimble.
It would be also interesting to compare the exact partition functions for ζ 1 = 0 with those in [16]. In that paper they computed the partition function of the U(N ) k × U(N + M ) −k linear quiver superconformal Chern-Simons theory. Naively, for M = 0 this theory can be obtained by taking the decoupling limit ζ 2 → ∞ in the mass deformed ABJM theory. Indeed we observe for N < k/2 that, if we take the limit ζ 2 → ∞ in the exact expressions for the partition function Z(N, k, 0, ζ 2 ) they precisely coincide with those in [16] up to the contribution of decoupled hypermultiplet e −2πN 2 ζ 2 /k (see appendix C.2). On the other hand we also observe that for some cases with N ≥ k/2 the decay is slower than e −2πN 2 ζ 2 /k . Actually in these cases the corresponding linear quiver theory is a bad theory [54], hence it should not be the right decoupling limit of the mass deformed ABJM theory. Though it is still not clear, the correct description of the decoupling limit might be obtained by expanding the Coulomb branch moduli around a configuration which depends on ζ 2 in a non-trivial way.
We believe that the expected Z = 0 behavior means the SUSY breaking. Indeed, this is an analogue of the Witten index for the theory on S 3 and the partition function on S 3 for the pure SUSY CS theory, which is believed to be SUSY broken for k < N [21], was shown to vanish by using the localization technique. 30 However, it is worth to note that 30 The mass deformed ABJM theory on S 3 with large masses at the origin of the Coulomb moduli will be connected to the pure SUSY CS theory. However, as discussed in section 3, there are many classical vacua, some of them are SUSY breaking including quantum corrections [11], and we can not conclude which one dominate (or no one dominates) the path-integral. There exist several evidences which suggest that the massive hypermultiplets in the limit ζ1, ζ2 → ∞ is not decoupled. If all the hypermultiplets with masses 2πζ1/k and 2πζ2/k decouples, the partition function should behave like Z(N, k, ζ1, ζ2) ∼ e −2πN 2 (ζ 1 +ζ 2 )/k in the limit ζ1, ζ2 → ∞. We can see that this is not satisfied for Z(2, 1, ζ1, ζ2) by using the exact expression (5.7). As we have already commented above, Z(N, k, 0, ζ2) for some N, k with N ≥ k/2 also do not follow the same behavior. This suggests that the even in the large mass limit, the hypermultiplets are not completely decoupled. This implies that the theory will not be described near the origin of the Coulomb moduli. We expect that the SUSY breaking (meta-stable) vacua in the flat space will contribute JHEP03(2019)159 Z = 0 does not necessarily mean SUSY breaking as in Witten index. Nevertheless, we expect that interpreting the first zero in the large N limit as the SUSY breaking would be the most natural possibility as discussed through the paper.
Although our analysis strongly suggests that the first zero of the partition function approaches ζ 1 ζ 2 = k 2 /16 in the large N limit and that a large-N phase transition takes place at these points, currently we do not have sufficient information to conclude whether this phase transition corresponds to the SUSY breaking or not. To prove/disprove our proposal, more detailed analyses such as studying other observables are called for. It is interesting even if the Z = 0 does not come from SUSY breaking since this would imply a new exotic phenomenon.
Lastly it is highly important to check our claims in the gravity dual picture. In [17], the authors found solutions in the effective four dimensional supergravity which will be the gravity duals of the R-charge deformed ABJM theory which are related to our setup by the simple analytical continuation at the level of the localization formula (2.5). 31 They did not report any phase transition and the free energy of the gravity dual coincides with the results for the R-charge deformations [19] and our results for the mass deformations with sufficiently small mass by the analytic continuation. We expect that there exists some subtleties for the validity of the naive analytic continuation of the gravity solutions which make the SUSY breaking phase transition. One natural possibility is that non-perturbative effects on the gravity side changes signs of real parts of their exponents at some values of the masses and invalidates the description by the classical supergravity as discussed in section 5.3. We hope to study these issues in near future.

JHEP03(2019)159
Next we rewrite the 1-loop determinant into pair of determinants of N × N matrices by using the Cauchy determinant formula N i<j 2 sinh and then combine them by using the formula After these manipulations we obtain the following expression for the partition function (A.6) Hence the partition function takes the form of the partition function of 1d N particle noninteracting Fermi gas if we regard f (λ , λ ) as the matrix element of a one-particle density matrix λ| ρ|λ with position eigenstates |· In this appendix we explain details on how to solve the recursion relation (4.32) for integer k.

JHEP03(2019)159
C Exact expressions for Z(N, k, 0, ζ 2 ) The technique introduced in section 4.2 allows us to compute the partition function of the mass deformed ABJM theory Z(N, k, 0, ζ 2 ) with ζ 1 = 0 and for small integers N, k.
We can make a more refined comparison between the exact results and large-N expansion as follows. First we notice that the saddle point approximation (4.20) agree with the following expression in the large-N limit where C = 2 which is obtained from the partition function of the ABJM theory with R-charge deformation by ignoring the large-N non-perturbative effects (e − √ N ) and replacing the real deformation parameters ξ, η formally as ξ → 0, η → 4iζ 2 /k (see eq. (1.4) in [34]). By comparing the numerical values of (C.1)-(C.5) and (C.7) we find good agreement. As an example, in figure 11 we display the comparison of the free energy for ζ 2 = 1.
For the other (N, k) we have found the following asymptotic behavior