Fermionic Rational Conformal Field Theories and Modular Linear Differential Equations

We define Modular Linear Differential Equations (MLDE) for the level-two congruence subgroups $\Gamma_\vartheta$, $\Gamma^0(2)$ and $\Gamma_0(2)$ of $\text{SL}_2(\mathbb Z)$. Each subgroup corresponds to one of the spin structures on the torus. The pole structures of the fermionic MLDEs are investigated by exploiting the valence formula for the level-two congruence subgroups. We focus on the first and second order holomorphic MLDEs without poles and use them to find a large class of `Fermionic Rational Conformal Field Theories', which have non-negative integer coefficients in the $q$-series expansion of their characters. We study the detailed properties of these fermionic RCFTs, some of which are supersymmetric. This work also provides a starting point for the classification of the fermionic Modular Tensor Category.

Abstract: We define Modular Linear Differential Equations (MLDE) for the level-two congruence subgroups Γ ϑ , Γ 0 (2) and Γ 0 (2) of SL 2 (Z). Each subgroup corresponds to one of the spin structures on the torus. The pole structures of the fermionic MLDEs are investigated by exploiting the valence formula for the level-two congruence subgroups. We focus on the first and second order holomorphic MLDEs without poles and use them to find a large class of 'Fermionic Rational Conformal Field Theories', which have non-negative integer coefficients in the q-series expansion of their characters. We study the detailed properties of these fermionic RCFTs, some of which are supersymmetric. This work also provides a starting point for the classification of the fermionic Modular Tensor Category.

Introduction and Concluding Remarks
The classification of all unitary conformal field theories in two dimensions certainly plays a key role in our understanding of the critical phenomena. Utilizing the conformal and modular bootstraps, the unitary conformal field theories with c < 1 in particular can be solved and classified completely [1]. They are called minimal models obeying the so-called ADE classification. One of the important features of minimal models is that they have finitely many conformal primaries. It admits a generalization that leads to a very rich class of two-dimensional conformal field theories, namely rational conformal field theories. A rational conformal field theory (RCFT) refers to a conformal field theory whose torus partition function can be expressed as a finite sum of products of holomorphic and antiholomorphic functions. Such holomorphic functions can be understood as characters with respect to an extended chiral algebra that includes Virasoro algebra. It is well-known that, in a given RCFT, the central charge as well as conformal weights are rational numbers [2,3]. One prominent example of RCFTs is the Wess-Zumino-Witten (WZW) model where the extended chiral algebra is the current algebra.
A few approaches based on chiral algebras and lattices have been proposed to solve a more tractable problem of the classification of RCFTs [4][5][6][7]. Since there are RCFTs that fit into neither of those approaches, the classification however remains incomplete. On the other hand, the authors of [8] have proposed a rather different approach based on the modular invariant linear differential equations (MLDE).
Let us first briefly explain how an MLDE can provide a systematic procedure to classify the RCFTs. Suppose that a given RCFT has N independent characters. From the modular invariance of the torus partition function, one can see that the holomorphic characters transform as a vector-valued modular form under SL 2 (Z). It implies that the N characters can be regarded as independent solutions to an N -th order differential equation invariant under SL 2 (Z). This method is particularly useful to obtain a complete list for RCFTs with small numbers of characters, which was studied extensively in [8]. Readers can also find a recent status on this programme in [9] and the references therein.
Rational conformal field theories are also closely related to the modular tensor categories (MTC) that have extensive applications to the study of anyonic systems and topological quantum computation [10][11][12]. When two different RCFTs are related to two MTCs conjugate to each other, it turns out that they satisfy a certain bilinear relation studied recently in [13][14][15][16]. Given that the MTCs of low rank are only classified in [17,18], we expect that the classification of RCFTs based on the MLDE sheds new light on the problem of the classification of MTCs.
In the present work, we extend the MLDE method to classify the fermionic rational conformal field theories, some of which appear to be supersymmetric. A fermionic CFT refers to a conformal field theory which contains operators of half-integer spin. To define a fermionic CFT on a manifold, one has to choose a spin structure. On a torus, there are four different spin structures, (NS,NS), (R,NS), (NS,R) and (R,R). For later convenience, we use the shorthand notation NS, NS, R and R for those spin structures. A fermionic CFT in our study is further restricted to have a certain extended chiral algebra that includes conserved currents of half-integer weight. In other words, there are half-integer spin descendants of the vacuum in the NS sector.
It is well-known that the Jordan-Wigner transformation maps the critical Ising model to a theory of free Majorana fermion. Later on, the Jordan-Wigner transformation has been revisited to fermionize a given bosonic CFT with a non-anomalous Z 2 to a fermionic CFT [19][20][21]. Along the way, the "Beauty and the Beast" N = 1 superconformal theory [22] can be reinterpreted as a fermionization of the Monster CFT [23], and the fermionic minimal models are constructed in [24]. The goal of this paper is to classify such fermionic RCFTs with small numbers of conformal characters systematically via the MLDE method. The classification of fermionic RCFTs will provide the classification of the fermionic modular tensor category that characterizes the fermionic topological phases of matter, modulo the fact that two different fermionic RCFTs sharing the same fusion rule algebra are related to a single fermionic MTC.
In order to discuss the extension of the MLDE method, let us note that the characters of a given fermionic RCFT in NS, NS, and R sectors become a vector-valued modular functions for the level-two congruence subgroups Γ ϑ , Γ 0 (2) and Γ 0 (2) of the modular group. This is because each of these congruence subgroups is associated to a specific spin structure on the torus, i.e., NS, NS, and R sectors. Hence each of the 'fermionic' MLDEs associated to them, transform into each other under SL 2 (Z) transformations.
The first step consists of understanding the relation between the pole structure of the coefficients of the MLDE and the set of zeros of the characters which are solutions to the MLDE. This was neatly understood in the case of SL 2 (Z), see for instance [8], in terms of the zeros of the Wronskian associated to the MLDE. Having automorphic properties, the Wronskian happens to be subject to the so-called valence formula, a classical result in the theory of modular forms constraining the possible set of zeros (together with their multiplicity) in terms of the weight of the automorphic form. This allows to reduce the choice of the pole structure characterizing the order N MLDE down to the choice of a single integer . In this paper, we introduce the equivalent of the valence formula for various relevant level-two congruence subgroups of SL 2 (Z). Equipped with this tool, we generalize the MLDE to the case of level-two congruence subgroups Γ ϑ , Γ 0 (2) and Γ 0 (2) of the modular group. The level-two valence formula boils down to the following relation between the central charge c, the conformal weights in the NS and R sectors h NS and h R , N and : (1.1) After discussing the possible poles structures for a fermionic MLDE in terms of the level-two congruence subgroups valence formula, we focus on the second order fermionic MLDE in the simplest case, that of a trivial pole structure. In this context, we classify the possible values of the central charge and conformal weights corresponding to some fermionic RCFT for which there exist character-like solutions to the fermionic MLDE. More precisely, we consider the solutions with the property that all the coefficients are non-negative integers in q-series. We classify the solutions by six classes, which can be found in table 4. In particular, we also derive a closed-form expression of the S-matrix of these fermionic RCFTs and investigate their fusion coefficients. We discard the solutions of the second order fermionic MLDE when they do not yield a consistent fusion rule algebra. For the solutions with consistent fusion rule algebra, we find the identifications in terms of the N = 1 supersymmetric minimal models or the WZW models.
One can construct the partition function of individual spin structures using the solutions of the second order fermonic MLDE. The sum of four partition functions defines an SL 2 (Z) invariant partition function of a certain bosonic CFT and this procedure is often referred to as the GSO projection [25], or equivalently bosonization. In some cases we study below, we find that the bosonization works with the assumption that the torus par-tition function for the spin structure R becomes constant, Z R = const. Especially when the following three conditions are satisfied, 1. a vacuum descendant of weight 3/2 is present, 2. the supersymmetic unitary bound h R ≥ c/24 is obeyed, 3. Z R = constant, then we suggest interpreting the corresponding solutions as the characters of unitary supersymmetric RCFT. For instance, we will show that the fermionization of su(2) 6 and (e 6 ) 3 WZW model could potentially be understood as the supersymmetric RCFT, as they satisfy the above three conditions. More examples will be discussed in the main context.
We also notice that some solutions can be combined into known partition functions, e.g., the Conway extremal CFT [26] or the N = 1 extension of the (e 8 ) 1 WZW model. This will be referred to as a bilinear relation in what follows. When such bilinear relation is satisfied for the characters of two different RCFTs, these two theories ought to share the same fusion rule algebra. Therefore, we expect that our classification will provide new insight into fermionic MTCs. On the one hand, some bilinear relations are known to appear as an evidence of deconstruction of the Monster group [15,16]. In a similar way, we test the splitting of the supersymmetric VOA for the Conway group Co 0 with c = 12. Specifically, the solution with c = 11 exhibits moonshine phenomena for the Suzuki group as shown in [27]. Further examples of fermionic deconstructions will be discussed in an upcoming paper [28].
In a separate upcoming paper [29], we generalize the work of the present paper to the case of a third order fermionic MLDE. There, we focus on a subfamily of solutions for which the BPS bound is saturated in the Ramond sector. We provide a closed-form expression of the characters and the S-matrix for these solutions.
This article is organized as follows. In Section 2, we review some mathematical facts concerning the modular group SL 2 (Z) and the well-known classification of second order MLDEs by [8]. In Section 3, we introduce some standard results of three level-two congruence subgroups and propose MLDEs with holomorphic coefficients, to which we also refer as (holomorphic) fermionic MLDEs. Then we move on to finding solutions that can possibly be identified as characters of fermionic RCFTs. In Section 4, we classify all possible solutions of fermionic first order MLDEs, which turn out to consist of products of Majorana-Weyl fermions. In Section 5, we study fermionic second order MLDEs with trivial pole structure and find six families of solutions, listed in Table 4. For all the consistent solutions we are able to express them in terms of the characters of some known RCFT.
We would like to dedicate this article to the memory of Professor Tohru Eguchi who passed away last year. K.L. recalls a personal meeting with him about 20 years ago during his visit to Tokyo University and many more wonderful interactions later on. Many of our works got influenced by Prof. Eguchi's works. In relation to this article, we note that Eguchi and Ooguri derived in [30] the third order MLDE for conformal characters of the Ising model as well as the exact form of the characters. It would be the first place where the MLDE made its appearance in the study of conformal field theories. Anderson and Moore wrote a general MLDE and used it show the rationality of c and h for rational conformal field theories [2]. Afterwards, Mathur, Mukhi and Sen [8] expanded and established the modular linear differential equation and utilized it as a tool for the classification of RCFTs. They found the famous MMS series of theories discussed in Section 2.2.
2 Modular Linear Differential Equation for SL 2 (Z)

SL 2 (Z) Group and its Valence Formula
It is natural to consider RCFTs on a torus. The partition function of the theory is required to be modular invariant. Thus, the modular group SL 2 (Z) plays an essential role in our understanding of the conformal field theory. Here we review its important features relevant to our analysis.
The modular group SL 2 (Z) is the group if invertible 2 × 2 matrices with integer coefficients and unit determinant: The modular group is generated by S = 0 1 −1 0 and T = 1 0 1 1 . The modular parameter τ ∈ H = {τ ∈ C|Im(τ ) > 0} transforms by fractional linear transformations as γτ = (aτ + b)/(cτ + d). The fundamental domain D = SL 2 (Z)\H for the modular group SL 2 (Z) is given as with additional identification of the boundaries by T : τ → τ + 1 and S : τ → −1/τ . The fundamental domain is drawn in Figure 1. The τ = i is an orbifold point of order two and τ = ω ≡ e 2πi/3 is an orbifold point of order three. The topology of the fundamental domain is a sphere with a single puncture, the cusp at i∞.
Given an element γ ∈ SL 2 (Z), one defines the following action on functions f (τ ) from H to C: where ρ(γ) is a possibly non-trivial γ-dependent phase, and k an integer called the weight of f . If the function f is periodic under T : τ → τ + 1 so that ρ(T ) = 1, we can have a Fourier expansion: The function f is meromorphic at i∞ if only a finite number of negative powers of q appears in the above expansion, holomorphic at i∞ if there is no negative power of q, and vanishes at i∞ if only positive powers of q appear.
Let us first focus on the SL 2 (Z) weight k forms which have the following property: In this work we are interested in the function space M k of modular forms of weight k mainly. It is generated by E 4 , E 6 Eisenstein series as follows: The dimension of M k for k ≥ 2 are well-known to satisfy .
The space S k of cusp forms of weight k is related to the space of the modular forms as follows: For the class of functions f : H → C which satisfy f | k γ = f for any γ ∈ SL 2 (Z) of given weight k (2.3), holomorphic in H and possibly meromorphic at i∞, there exists a so-called valence formula relating the number of zeros in the fundamental domain SL 2 (Z)\H to the weight k. Although it is a standard textbook material [31,32], let us review the derivation briefly as we want to extend it to the case of level-two congruent subgroups later.
To get the valence formula, we use the Cauchy formula on the integration of ∂τ f f along the contour given in Figure 1, which is inside the fundamental domain SL 2 (Z)\H. The resulting integral can be evaluated in two ways, where ω = e 2πi/3 . For the first equality, the contour integration is split into integrals surrounding all the zeros in the interior of H and becomes − interior p ν p (f ). The second equality comes from the contour integration along the boundary. One can read off the origin of each terms in the RHS of the above equation from the dashed lines in Figure 1. The contribution from each segment is: ν ω (f ), The order ν p (f ) = n if f (τ ) near τ = p has a Laurent expansion with leading term f (τ ) ∼ (τ − p) n . For f ∼ q n near i∞, log(f ) ∼ 2πin log τ and so ν ∞ (f ) = n. At orbifold points i and ω, the contribution gets only 1 2 and 1 3 fractions, respectively as one integrate over one-half or one-third of 2π integration. On the arc C B C and C C D , a point τ and its S-dual points −τ −1 are matched. One changes the variable τ → −1/τ on C B C and use the fact f (−1/τ ) = ρ(S) τ k f (τ ) to get that the contribution from the contour C B C is minus of that from the contour C C D together with −k/12. The factor 1/12 is due to the fact that the angle between ω = e 2πi/3 and i is 2π/12.
In short the valence formula for weight k form of SL 2 (Z) is (2.10) Table 1 shows the validity of the valence formula for well-known modular forms which have no zeros or poles besides i∞, i and ω.

Modular Linear Differential Equation
A two-dimensional rational conformal field theory of central charge c has a finite number of primary operators. Its modular invariant partition function on a torus can thus be expressed as where χ a (χ a ) denote left-moving (right-moving) characters with respect to an extended chiral algebra including the Virasoro algebra. The modular invariance of (2.11) implies that the characters χ a (τ ) transform as a finite-dimensional representation of SL 2 (Z), i.e, they transform under S and T as follows: where the modular matrices S ab and T ab are symmetric and satisfy the relations below where h a is the conformal weight for χ a . Note also that S ab and T ab should satisfy where C is the charge conjugation matrix.
One can see from (2.12) that the characters of the RCFT constitute a vector-valued modular forms of weight zero and thus satisfy a modular-invariant linear differential equation (MLDE) of order N where N is the number of linearly independent characters. Let us review the argument for the derivation of the modular linear differential equation in [8,33]. We start with an (N + 1)-dimensional square matrix made of χ 0 , χ 1 , ..., χ N −1 , f and their derivatives with the Ramanujan-Serre covariant derivative, acting on a weight r modular form, up to the N -th power. This covariant derivative transforms a weight r modular form to weight r + 2 modular form. (See the Appendix B for the details.) If the function f is a linear combination of N characters, the determinant of this (N + 1)-dimensional matrix vanishes, implying that where each coefficient W k is given by (2.17) One can recast (2.16) into where the coefficients φ k (τ ) = (−1) N −k W k /W N are automorphic forms of weight 2N − 2k for SL 2 (Z). These forms φ k (τ ) could have poles at zeros of the Wronskian W N (τ ).
Let us first apply the valence formula (2.10) to the Wronskian W N (τ ), which transforms under S as a modular form of weight N (N − 1) and is invariant under T up to a constant phase ρ(T ). At τ = i∞, the Wronskian has the asymptotic expansion: because each characters χ a is asymptotic to This implies that The zeros of the Wronskian W N become the poles of the coefficient functions φ k . As the modular forms are covariant under the modular transformation, the coefficient φ k of weight 2N − 2k can be expressed as rational functions of E 4 and E 6 , whose denominator is constrained by the parameter .
The second order MLDE with n = 2 has been studied extensively in [8,33]. For the simple case without poles = 0, the MLDE becomes Note that the coefficient of the first derivative vanishes as there exists no modular form of weight two and E 4 is the unique modular form of weight four up to a constant factor. Two independent solutions to the above equation can be regarded as two characters that can be expanded in powers of q as follows Since the valence formula (2.22) says h = c+2 12 , the free parameter µ of (2.24) can be determined as µ = − c(c+4) 576 = (1+6h)(1−6h)

144
. It has been shown in [8] that there exist only ten values of allowed central charges, such that the solution of (2.25) has all q-expansion coefficients given by non-negative integers. Although the first case with c = 2/5 and h = 1/5 appears to be consistent, it has negative fusion coefficients. To resolve this, one interchanges the two characters and obtain a theory with c = −22/5 and h = −1/5. As discussed in [8,33], the above ten solutions provide the characters for the Lee-Yang edge singularity and level-one WZW models for a 1 , a 2 , g 2 , d 4 , f 4 , e 6 , e 7 , e 7 1 2 , e 8 , respectively. This series of Lie groups is known as the Deligne-Cvitanovic series [34]. The explicit solutions to (2.18) as well as their Smatrix for each central charge c can be found in [35]. One can also show that the modular invariant partition function on the torus is diagonal, i.e., where M is a certain non-negative integer. We present the central charge c, and corresponding values of h, and multiplicity M in Table 2.
As a remark, the WZW model for e 8 at level one appears in Table 2, but in fact it is a single-character RCFT. The vacuum character which is its unique character is given by The characters of the above series satisfy also the bilinear relations This shows multiple ways to divide the e 8 level one WZW model into two pairs of complementary CFTs.
In the following sections, we generalize MLDEs to the fermionic theories including supersymmetric theories by considering the modular forms and the valence formula for the related congruence groups Γ θ , Γ 0 (2), and Γ 0 (2).
3 Congruence Subgroups Γ θ , Γ 0 (2) and Γ 0 (2) 3.1 Γ ϑ , Γ 0 (2), Γ 0 (2), Γ(2) modular subgroups and forms The level-two congruence subgroups Γ ϑ , Γ 0 (2) and Γ 0 (2) are related to the modular symmetry group of the fermionic conformal field theories in the NS-NS, R-NS and NS-R spin structure sectors of the partition functions, respectively. We will also refer to these three spin structures as NS, NS and R sectors respectively. Their study provides a good starting point for building up fermionic modular linear differential equations.
The principal congruence subgroup of SL 2 (Z) of level two, Γ(2), is defined as follows: It is generated by S 2 = −1, T 2 , ST 2 S. Note that S 2 = (ST ) 2 = −1. The weight two Γ(2) modular forms, ϑ 4 a (τ ), a = 2, 3, 4, are built out of the Jacobi theta functions and transform under S and T as We refer the reader to Appendix A for the definition and properties of various modular objects used in the core of the article.
The various modular functions λ, j, K satisfy the following simple identities: (3.14) We will use the λ function to re-express the MLDEs in terms of the new variable τ → λ(τ ). It so happens that in some nice situations, doing so will allow us to explicitly solve the ODE, hence providing a closed form expression of the characters. For instance, it is known that the MLDE (2.24) can be written in terms of λ and take the form: which then can explicitly be solved in term of the hypergeometric function 2 F 1 .

Valence Formula for Γ θ
Let us consider the Γ θ valence formula for definiteness. The valence formula for Γ 0 (2) and Γ 0 (2) could be found similarly. The fundamental domain for Γ θ is depicted above in Figure  3. With the S, T 2 generators leave the fundamental domain invariant, one has to fold along τ = 0 vertical line, identifying Re(τ ) = ±1 vertical line by T 2 transformation and the left and right half of the semicircle by S transformation. The fundamental domain D Γ θ has two cusp points at i∞ and 1, and a single orbifold point τ = i of order two.
Analogous to the SL 2 (Z) case, we consider a class of the functions f : H → C that f are holomorphic inside H, at worst meromorphic at the cusps i∞ and 1, and satisfy where ρ : Γ ϑ → U (1) is a possible phase. To derive a valence formula for Γ θ , we are interested in the integration of df /f along the contour in the fundamental domain of Γ θ as shown in Figure 4. We follow the argument in the derivation of the valence formula for Table 3. Weight k (weak) modular forms of Γ θ SL 2 (Z) in Section 2. The contour integration leads to the equalities, When the contour gets shrunk, it receives contributions only from the zeros of f lying on the interior region, that is, not a orbifold point or cusps. This is the first equality. The second one comes from the contour contributions. One can read off the origin of each terms in the RHS of the above equation from the dashed line in Figure 4. The contributions from contour segments are as follows: For f ∼ q n near i∞, log(f ) ∼ 2πin log τ . The integration length of C EA is two and so its contribution is 2ν ∞ (f ) = 2n. At another cusp 1, the contribution becomes just ν 1 . The structure of the zero and the contour at 1 just comes from the ST transformation of the i∞ contour integration in SL 2 (Z). On orbifold points i the contribution gets only 1 2 fraction of ν i (f ) as one integrate over one half of 2π integration. On the arc C BC and C C C , a point τ and its S-dual points −τ −1 are matched. One change the variable τ → −1/τ on C B C and use the fact f (−1/τ ) = ρ(S) τ k f (τ ) to get the contribution from the contour C B C is minus of that from the contour C C D and −k/4. The factor 1/4 arises because the angle between 1 and i is one fourth of 2π.
In short, we obtain the valence formula for Γ ϑ as follows:  Table 3 is a check for the valence formula for modular forms of Γ θ . As studied in the previous subsection 3. 1 The space S k (Γ θ ) of cusp forms for Γ θ is made of the modular forms of Γ θ which vanish at cusps, i∞ and 1, given by The dimensions of M 2k and S 2k are given as follows:

Modular Linear Differential Equation for Γ θ
Let us now see how to use the above valence formula in the context of fermionic RCFTs. For that purpose, we consider a fermionic RCFT which has finite number of characters χ NS i , i = 0, 1, ...N − 1 with conformal weight h i . We repeat the argument for MLDE of SL 2 (Z) in Section 2.2 for these N characters in NS-sector to get the MLDE for Γ θ : The coefficient functions φ k = (−1) N −k W k /W N are the weight 2N −2k rational functions of modular forms of Γ θ . The Wronskian W k here is given in terms of the NS-sector characters χ NS i . We adopt the same strategy as in i.e., studying or classifying MLDE's based on the order of poles of φ k = (−1) N −k W k /W N . Since the characters χ ns i hence W k are always holomorphic inside the fundamental domain, the only possibility of introducing a pole is through the zeros of the denominator W N . It's easy to check that W N transforms as a (meromorphic) weight N (N − 1) modular form of Γ θ up to a possible overall phase. Therefore, we can safely apply the valence formula (3.19) to W N .
Due to the very nature of characters, we are able to further simplify the formula. Notice that at τ = i∞, by definition χ ns j have the following asymptotic expansion, Taking Serre derivative on χ ns j doesn't change their order at i∞, so we learn that In other words, we have: Concerning the cusp at 1, a priori we need to find a suitable local coordinate to extract the leading order. However, there is a way to bypass this problem by noting the fact that 1 is mapped to i∞ under the (ST −1 ) transform. Then the local behaviour of χ NS j around 1 is equivalent to that of ST −1 (χ NS j ) ∼ χ R j around i∞. Thus the Wronskian of NS characters near τ = 1 becomes that of R characters near τ = i∞: Immediately, we see that Mimicking what was done in [8] and also in (2.23), if we define l/2 to be the"number" of zeros inside the fundamental domain, the valence formula (3.19) gives us then the following relation: where the contribution from zeros is given as For a given fermionic RCFT with known c, h NS j and h R j , the value of the zero determines the structure of MLDE satisfied by the characters in NS sector. In this work, we will focus on the case where = 0. In this case the coefficient functions φ k of the MLDE for Γ θ become the modular functions of M 2N −2K (Γ θ ). We explore the possible values of c, h NS j , h R j whose characters have the non-negative integer q 1 2 -expansion coefficient. For the conformal characters in NS-sector, there always exists the unique vacuum character with zero conformal weight h = 0 and q 1 2 expansion: The MLDE for NS, R sector can be obtained by just making T and S transformation of the MLDE for the NS sector. The modular forms φ k would transform accordingly.
Before exploring the solutions of MLDE, let us compare bosonic and fermionic valence formulas (2.22) and (3.28). If we start with N characters in the NS sector, we can bosonize and get 3N characters by combining all characters in NS, NS, and R characters. By a linear combination of χ NS a and χ NS a we get two characters with conformal weight h a for the sum and h a + s a for the difference, where s a is a positive half integer, say, 1 2 or 3 2 , in our case. So we could use the bosonic MLDE for the 3N characters with h NS a , h NS a + s a and h R a with the number of zeros b . For simplicity, we are assuming that all 3N characters are linearly independent. Then two valence formulas imply a consistency condition, As we consider the solutions of MLDE for the fermionic RCFT in the q-expansion, we can read s a trivially. The above relation then provides the information on the b , that is, the pole structure of the coefficient functions for MLDE of bosonic 3N characters.

Fermionic First Order MLDE
As a warm-up exercise, we consider the fermionic RCFT whose torus partition function for each spin structure can be holomorphically factorized. That is to say, the theory has a single character in each sector. For convenience, let us focus on the NS sector. The vacuum character f NS 0 (τ ) ∼ q − c 24 (1 + · · · ) then becomes a solution to the first order MLDE below where φ 0 (τ ) is the Γ θ modular form of weight two. Restricting our attention to the case l = 0, φ 0 (τ ) can be generated by − ϑ 4 2 (τ ) + ϑ 4 4 (τ ) since the space M 2 (Γ θ ) is onedimensional. Therefore the first order MLDE of our interest takes a form of Note that its T and S transformations lead to the first order MLDE for Γ 0 (2) and Γ 0 (2): To have a conformal character f NS 0 ∼ q − c 24 (1 + · · · ),one can determine the coefficient µ to be µ = c 12 . With the help of (3.14), (4.2) can be rewritten in terms of the λ(τ ) variable as follows, The solution is simply given by (4.5) Demanding the solution to have non-negative integer coefficients for the q-expansion leads to c = N 2 for a positive integer N . It is consistent with the fact that the chiral central charge of a fermionic CFT has to be half-integral. This is because the (2 + 1)-dimensional fermionic gravitational Chern-Simons coupling can cancel certain (anomalous) U (1) phases that arise from the modular transformation of partition functions of fermionic CFTs with c L − c R ∈ Z/2.
When N = 1, one can see that f NS 0 is nothing but the NS partition function of a single free Majorana-Weyl fermion, The other partition function for different spin structures can be obtained by acting T and S on ψ NS , Note also that the characters of a free Majorana-Weyl fermion are related to the characters of the Ising model as follows Therefore, the solution (4.6) can be identified with the NS partition function of N -copies of free Majorana-Weyl fermions, with the identities 5 Fermionic Second Order MLDE

The Second Order MLDE
In this section, we consider the fermionic RCFT with two characters in each sector. Let us first focus on the NS sector. The characters for the vacuum and primary state of conformal weight h NS have a series expansion and they are independent solutions of the second order MLDE of the form We restrict our attention to the zero-pole case = 0 where the coefficients φ 0 and φ 1 are spanned by Γ θ modular forms of weight four and two, respectively. The weight two modular form M 2 (Γ θ ) is one-dimensional and generated by −ϑ 4 2 + ϑ 4 4 . The space of weight four modular form M 4 (Γ θ ) is two-dimensional whose basis is chosen to be ϑ 8 3 and E 4 . Therefore, the second order MLDE for NS sector takes the expression: where µ 1 , µ 2 and µ 3 are three independent coefficients. We will mainly focus on the solutions of (5.3) which have the non-negative integer coefficients in the q-expansion. Note that the first term of the vacuum solution should be given by q −c/24 due to the uniqueness of the vacuum state. On the other hand, there is an ambiguity for the normalization of the primary solution denoted by b 0 in (5.1). We will show that the normalization b 0 for the NS sector solution can be fixed with the help of analytic expression of the S-matrix.
The MLDE for NS and R sector can be obtained from (5.3) by applying T and ST transformation. Explicitly, the differential equations are given by which are consistent with the valence formula. In the case of µ 1 = µ 2 = 0, the above three MLDE reduce to the second-order MLDE for SL 2 (Z) whose solutions were discussed in Section 2.2. Thus we restrict to the cases where at least one of µ 1 , µ 2 is non-zero.
Our strategy of solving the NS sector MLDE is to ask if the characters of the series form (5.1) are solutions of (5.3) in every order of q. In addition, we demand the coefficients of the series expansion to be non-negative integers. Let us choose the input to be central charge c, weight of the primary h NS , and a 1 ≥ 0, where a 1 is the coefficient of q 1 2 term in the vacuum character. For non-zero a 1 , the corresponding theory contains primaries of dimension 1/2 and spin 1/2, i.e., free fermions. See, for instance, [36]. Once we require the NS-sector characters to satisfy the MLDE for Γ θ , the coefficients µ 1 , µ 2 , µ 3 are determined by inputting parameters (c, h, a 1 ) as follows: Now let us move our attention to the R sector. The two characters of the R sector have series expansions of the form where h R ± denote the weights of two primaries in the R sector. Without loss of generality, we assume h R + ≥ h R − , i.e., h R − is the Ramond vacuum weight. With the help of (5.5), one can show that the weights h R ± of the R sector are determined from the MLDE (5.4) as below, and the weights h R ± should be a rational number as far as we consider the RCFT. Note that the valence formula (3.28) is satisfied trivially with (5.7). When we assume the resulting theory has supersymmetry, there is a unitary bound for the R sector weights which leads to the constraint on the initial data (c, h NS , a 1 ), For a specific case where a 1 = 0, or h NS = 1 2 , this bound takes a simpler form: Let us further comment on the RCFT with supersymmetry. A necessary condition for a fermionic RCFT to have supersymmetry is that there is at least one weight-3/2 supersymmetry current as a vacuum descendant. Such supersymmetry currents contribute to a 3 q 3 2 in the NS vacuum character. If a supersymmetric RCFT further saturates a unitary bound (5.8), the RCFT has the supersymmetric ground states. Otherwise one can say that the supersymmetry is spontaneously broken. For this reason, we call h R − = c 24 as the BPS condition. When the BPS condition is satisfied and h NS = 1 2 , one can fix h R + and a 1 in terms of c and h NS as follows.
We will discuss later some examples of solutions that violate the unitary bound, thus cannot be considered as the characters for a supersymmetric RCFT.

Solutions of the Second Order MLDE
Let us reformulate the second order MLDEs in terms of λ(τ ) variable to find a closed-form solution, We obtain the MLDE for f NS (λ) and f R (λ) by acting T and ST transformations on the MLDE for f NS (λ). When h NS = 1 2 , the solutions can be expressed in terms of the hypergeometric function as follows, where b 0 denotes a normalization constant of the non-vacuum character and From the known identities of the hypergeometric functions, one can see that (5.13) transform under the S transformation as follows The derivation of (5.13) and (5.15) was inspired by [37], and we briefly review it in the Appendix B for completeness. Whenever we encounter zeros or poles in the denominator or numerator, a careful limit has to be taken.
For the case of h NS = 1 2 , the MLDE takes a form of and there are two independent solutions for the above differential equation, Thus any NS sector character of two character RCFT with h NS = 1 2 should be written as a combination of g 1 (λ) and g 2 (λ).
The structure of NS sector partition function can be fixed by the S-matrix presented in (5.15). More precisely, the NS sector partition function takes a form of where M denote the degeneracy of the NS primary states. Note that the multiplicity M can be determined from the S-matrix as below.
With the help of the identity of hypergeometric function we find the analytic form of the NS sector solutions that work for h NS = 1 2 .
One can then express the partition function of NS and R sector in terms of the NS sector characters f NS 0 (λ) and f NS 1 (λ). Explicitly, where M 1 and M 2 are related to the degeneracy of two characters in the R sector. The q-expansion of characters are presented in the Appendix C.
Let us now describe how we find the initial values (c, h NS , a 1 ) which generate the solutions f NS 0 (τ ) and f NS 1 (τ ) as the q-series with the non-negative integer coefficients. Our approach is first to focus on the vacuum character f NS 0 (τ ) by taking the series expansion of (5.13). One can survey the solutions of MLDE which have rational 0 < c = p p ≤ 24 and non-negative integers for the coefficients a n in (5.1). For instance, when we impose the BPS condition in the R sector, coefficients a 2 of the vacuum character is given as follow : and it should be non-negative integer for a lot of positive rational numbers c = p p and h N S = r s < 1. Assuming the upper bound on the denominators p ≤ p * , s ≤ s * and the spin half current a 1 ≤ a * 1 , we search for values (c, h NS , a 1 ) that make a 2 a non-negative integer. Alternatively, one can choose the initial data to be (c, a 2 , a 1 ) and check a 3 being integral and non-negative for rational 0 < c = p p < 24 with upper bounds on the denominator i.e., p ≤ 12. Then we find the initial data that furnishes non-negative integer a 3 with a 2 in the range 0 ≤ a 2 ≤ 10 6 . As an illustrative example, let us consider the case of a 1 = 0. Then the vacuum solution takes a form of and the coefficients a n for n ≥ 4 are in general the rational function of two parameters a 2 and c. Here we demand at least one a 2 +1 = 0 for ∈ Z >0 to obtain the fermionic characters. Then we search the values of c and a 2 which give us the series with nonnegative integer coefficients.
Following the process described above, our next goal is to find the solutions of the second order MLDE and classify them. We first divide the possible solutions of MLDE into two types: the BPS solution and the non-BPS solution. For the type of BPS solution, we explore both the cases of a 1 = 0 and a 1 = 0. For the non-BPS type, we work out only the case with a 1 = 0. In total, there are six types of solutions and their characteristic profiles are presented in Table 4. The detailed descriptions for each type will be given below, while the list of solutions can be found in the Appendix C.
With the choice of connection α = c 2 − 6h NS + 1, one can show that (5.3) is equivalent to: The deformation α in the covariant derivative also precisely appears in front of the last term. Therefore this second order MLDE can be viewed in some sense as a direct generalization of the Kaneko-Zagier equation. Indeed, the standard Deligne-Cvitanovic exceptional series is reproduced as a particular case of our equation, for which α is simply set to zero, namely when c = 12h NS − 2.     Here M denotes the degeneracy of the non-vacuum primary. When M = 1, the S-matrix is symmetric otherwise one should find the extended S-matrix for the consistent fusion rule algebra. The bracket means that the corresponding theory is effectively a single character theory, appearing as a solution to the second order MLDE, with artificial second solution.
persymmetry. In this subsection, we first focus on the case of a 1 = 0. Sensible solutions appear for the five central charge values of c = 1, 9 4 , 6, 39 4 , 11 and they are listed in Table 5. For the solutions in this type, the NS sector weight is given by h NS = (2c + 3)/30. On the other hand, the R sector weights read h R − = c/24 and h R + = 3(c + 4)/40. The explicit q-series of the solutions can be found in Table 10 and Table 11 in Appendix C.
Fermionic RCFT with c = 1 The solution with c = 1 can be expressed in terms of the characters of the N = 1 supersymmetric minimal models. We denote the unitary supersymmetric minimal model by SM( + 2, ) where the central charge is given by The partition function of a unitary supersymmetric minimal model can be expressed by the characters χ NS m,n and χ R m,n [38], Moreover, using the known S-matrix of SM(6, 4), one can verify that these two solutions transform into themselves under S transformation, The above modular S-transformation matrix is neither unitary nor symmetric, but leaves the matrix diag(1, 2) invariant. Thus Γ θ invariant NS partition function is given by (5.31) The multiplicity M = 2 in (5.31) suggests that there are two different primaries associated with the same character f 1 (τ ). In other words, the theory of our interest has three NS characters, denoted byf 0 ,f 1 andf 2 withf 0 = f (1) 0 1 , that transform under the S transformation as follows Here the sign of exponents can be flipped as the role off 1 andf 2 exchanges. Note that the modular S-matrix (5.33) is symmetric and satisfies S 2 = 1. Once we have the symmetric S-matrix, the fusion rule algebra coefficients are followed by the Verlinde formula where the index 0 denotes the vacuum. Using the Verlinde formula, we can obtain the well-defined fusion algebra given below. Let us comment on the NS sector and R sector partition function. The transformation rule of the NS sector and R sector characters are given by .

(5.36)
One can read the NS sector partition function by taking T -transformation to (5.31).
The Ramond sector partition function can be obtained from the S-matrix (5.36). Because the S-matrix leaves diag(1, 1) invariant, the R sector partition function takes a form of (5.38) In fact, the above fermionic model with c = 1 can be understood with the N = 2 supersymmetric A 1 minimal model, a.k.a. Kazama-Suzuki model with level one. The N = 2 supersymmetric A k minimal model has the central charge and super-Virasoro characters of weight h NS (a, b) and with a, b ∈ Z+ 1 2 and 0 < a, b, (a+b) ≤ k+1. The explicit form of the N = 2 super-Virasoro character in the NS sector is given by Here y = e 2πiz with the chemical potential z for the U (1) R charge. The Kazama-Suzuki model with k = 1 has c = 1, and three NS characters of (h NS , Q NS ) = (0, 0), ( 1 6 , ± 1 3 ). One can see that, after turning off the chemical potential, f NS 0 (τ ) and f NS 1 (τ ) become the N = 2 NS characters χ NS 0,0 (τ, z = 0) = f NS 0 (τ ), χ NS  , as expected. It was shown recently in [27] that this model indeed has N = 1 supersymmeteric vertex operator algebra.
We find that two solutions (5.44) can be expressed in terms of the characters of the N = 1 supersymmetric minimal models SM(6, 4) and SM(8, 6) as follows.
We also note that c = 9 4 is one of the central charge for the minimal N = 2 model (5.39). Indeed we can show that two NS characters with the SU(2) chemical potential turned on can be expressed in terms of the N = 2 NS characters with U (1) R chemical potential tuned on:      Fermionic RCFT with c = 6 One can identify the solutions to the second order MLDE with c = 6 as the fermionization of the six-fold product of the a 1 WZW model with levelone. This model has been recently studied in [39] toward understanding the origin of the Mathieu moonshine of the K3 CFT. In particular, the quantum hexacode was utilized to construct the N = 1 supersymmetry current out of 2 6 primaries of weight 3/2. The levelone a 1 WZW model, denoted by su(2) 1 for later convenience, has central charge c = 1 and contains two characters, one of which corresponds to the vacuum and the other to the primary of conformal weight h = 1/4. The explicit form of two characters is Under S : τ → −1/τ , ψ 0 and ψ 1 transform as The six-fold product of su(2) 1 has 2 6 primaries, many of which become degenerate. There are eventually seven different characters of conformal weights h = 0, 1/4, 1/2, 3/4, 1, 5/4, 3/2, given by We can express two solutions with c = 6 in Table 5 in terms of χ n 4 (τ ), The modular S-matrix can be read off from (5.54) as which implies that the NS partition function has to be The multiplicity 15 in (5.58) which says there are 15 primaries associated with f NS 1 (τ ) is reflected in the non-symmetric and non-unitary modular matrix in (5.57). It implies that we have to find a 16×16 extended S-matrixŜ. To this end, let us first denote the extended characters as 1,α (τ ) = ψ 0 ⊗ · · · i ψ 1 · · · j ψ 1 · · · ⊗ ψ 0 + (0 ↔ 1), where the index α means α = {i, j} for 1 ≤ i < j ≤ 6. For the 16 component vector-valued modular form f 1,α (τ ) , the below 16 × 16 matrix whose components are S 00 = S αα = + 1 4 , S 0α = S α0 = + 1 4 , acts as the extendedŜ-matrix. Now the fusion rules follow from the Verlinde formula, Performing T and S transformation on the NS characters, one can obtain the NS and R partition functions where the c = 6 solution is self-dual. Note that the NS weights h NS andh NS of the dual pairs satisfy h NS +h NS = 1. We find that the NS sector solutions of BPS type I (f 1 (τ ) . Therefore, the bilinear relation suggests that the two theories in a given pair have identical fusion rule algebra. Explicitly, the fusion rules of the c = 11 putative theory are given by (5.3) while the hypothetical theory of c = 39 4 has fusion rules (5.46). The N = 1 superconformal field theory with c = 12 has been constructed by [40] and its partition function can be factorized into K(τ )K(τ ). Analogous to the deconstruction of the Monster CFT by two bosonic RCFTs discussed in [15,16], (5.65) suggests that one can deconstruct the c = 12 superconformal field theory by two fermionic RCFTs. On the one hand, the c = 12 superconformal theory is known to has symmetry group Co 0 = 2.Co 1 [41]. Therefore, it is natural to expect that the dual pairs exhibit moonshine phenomena for some subgroup of Co 0 . For instance let us consider the putative c = 11 theory. Our solution for c = 11 can be identified to the characters of N = 1 SVOA with c = 11 given in [27]. The automorphism group of N = 1 SVOA with c = 11 is known as the Suzuki group, which is one of the maximal subgroup of Co 0 . It would be interesting to investigate further deconstruction examples of c = 12 superconformal theory, and the details will be discussed in an upcoming paper [28].

BPS Type II : BPS solutions with free fermion currents
Let us then consider a mild generalization to explore the space of BPS solutions that contain free fermions. To do so, we consider the case that each NS sector vacuum character accommodates for the q 1/2−c/24 term, namely, As shown in (5.11), the weights of the primaries in each sector are given by when the unitary bound is saturated. As a consequence, the solutions of the second order MLDE are in general expressed as the rational function of c and a 1 . We search for the values of c and a 1 leading to non-negative integer Fourier coefficients of each solution. In what follows, we discuss BPS solutions by imposing c > 0 and a 1 = 0 which are summarized in Table 6. The explicit solutions in q-series can be found in Table 12 and two R characters of weight h = −3/32 and h = −7/32, The solution of c = 3/4 in Table 6 can then be identified as follows, From the modular S-matrix below, Note that all fusion coefficients become non-negative integers.

Supersymmetric ADE WZW models with level one
We claim that the solutions with c = 3/2, 3, 6, 9, 21/2 can be understood as the N = 1 supersymmetric ADE WZW model with level one.
Let us start with the tensor product of free Majorana fermions and the WZW model with level-one. The number of Majorana fermions is given by the rank of a group on which the WZW model is defined. The theory then has two NS characters which take the forms, · χ g 1 (τ ), (5.73) where ψ NS (τ ) = ϑ 3 (τ )/η(τ ) is the NS character for a single Majorana fermion as given in (4.6). χ g 0 (τ ) and χ g 1 (τ ) denote the vacuum and non-degenerate characters of level-one WZW model on g. When g = a 1 , a 2 , d 4 , e 6 , e 7 , the above characters (5.73) reproduce the BPS solutions for c = 3/2, 3, 6, 9, 21/2, respectively.
We notice that the level-one WZW model on a simply-laced g can be realized as free bosons on the root-lattice of g, which is referred to as the Frenkel-Kac construction [42,43]. Thus, we can describe the fermionic RCFTs with c = 3/2, 3, 6, 9, 21/2 as the N = 1 supersymmetric theories of rank rk(g) pairs of a free boson and a fermion.
By construction, the modular matrices of the above BPS solutions are governed by those of the corresponding WZW model: 1. For instance, the modular S-matrix of the c = 3/2 BPS solution becomes It is straightforward to read the fusion rule algebra coefficients using the Verlinde formula. The non-vanishing coefficients are given by (5.76) 2. The c = 3 fermionic RCFTs has two primaries associated with the same characters f NS 1 (τ ). They rotate as 3 and3 under SU(3). In the basis of (f NS 0 ,f NS 1 ,f NS 2 ) = (f NS 0 , f NS 1 , f NS 1 ), the modular S-matrix is given bŷ  The NS partition function is thus given by (5.79) 3. For the N = 1 supersymmetric level-one WZW model on d 4 , we have three primaries that are associated with the same character f NS 1 (τ ) but transform differently under SO(8) as the vector 8 v , the chiral spinor 8 s , and the anti-chiral spinor 8 c . In the basis of (f NS They do not provide consistent two-character theories, because the other character f NS 1 (τ ) does not admit positive integer coefficients for any n. Nonetheless, the character f NS 0 (τ ) is Γ θ invariant and serves as the NS partition function of a certain putative fermionic RCFT.
We argue that the integer parameter n can be further restricted to n = 0, 2, 4, 8, and 24 in order to make the corresponding theories physically well-behaved. To see this, note that the parameter n of (5.83) is the number of free fermions contained in a given model. (5.83) can thus split into two pieces as follows wheref NS (τ ) should be understood as the NS partition function of a c = (12 − n/2) CFT with no free fermion. We found thatf NS (τ ) fails to have positive integer coefficients in the q-expansion for n > 8 except n = 24. Another simple diagnostic to rule out an unphysical NS partition function is to analyse the degeneracy and weights in its R sector partition function [44]. One can show that, for n = 1, 3, 5, 6, and 7, the Ramond partition function has non-integer degeneracy and thus any of fermionic RCFTs without free fermions cannot have thef N S (τ ) as its NS partition function.
As discussed in the previous section, we can identify f NS (τ ) with n = 24 as the NS partition function of the 24 free Majorana-Weyl fermions (4.10). We also note that a system of 24 free Majorana fermions can be related to the O(24) WZW model with level one via the non-abelian bosonization [45]. More explicitly, one can express f NS 0 (τ ) = K(τ ) + 24 in terms of the WZW characters as follows It is worth to mention that (5.83) for n = 0 is nothing but the NS partition function of the c = 12 extremal SCFT, a.k.a. the Conway extremal CFT [26,41]. There are two different but equivalent constructions of it. One of them is based on the N = 1 supersymmetric e 8 WZW model with level one, followed by a Z 2 orbifold [26]. The other is a Z 2 orbifold of 24 free chiral fermions [41]. In both constructions, Z 2 action is the charge conjugation. The second construction implies that the Conway extremal SCFT is the fermionization of the SO(24) WZW model with level one [27]. To be concrete, the NS and R partition function of the extremal SCFT with c = 12 can be expressed in terms of the WZW characters, where M (c=3/2) = 1, M (c=3) = 2, and M (c=6) = 3. Indeed, (5.89) is a supersymmetric generalization of the bilinear relation (2.29), and thus potentially gives new insight into the supersymmetric modular tensor category.
We finally remark on the BPS solutions for c < 0. A complete list of BPS solutions with c < 0 and a 1 = 0 is presented in Table 7. Here both NS sector and R sector solutions have non-negative integer coefficients in q-series. One can show that every solution in Table 7 actually describes a unitary theory upon exchanging the vacuum and non-vacuum character. The precise mapping between BPS c < 0 solutions and BPS unitary solutions is given in Table 7.

Non-BPS solutions
A typical type of fermionic RCFTs can be generated by tensoring the WZW models for a Lie algebra g with arbitrary number of free fermions. Those models do not saturate the supersymmetric unitary bound h R ≥ c/24 in general except the cases discussed in the previous section. For instance, let us consider a tensor product of the (bosonic) a 1 WZW model with level k and three free Majorana fermions. Although the model is known to preserve the N = 1 supersymmetry, the Ramond vacuum energy is greater than c/24, i.e., the bound is not saturated for any k.
Although the space of solutions which contain free fermions but do not saturate the supersymmetric unitary bound deserves further investigation, we restrict our attention to a  simpler problem of classifying the fermionic RCFTs with no free fermion where h R > c/24 in the present work.

Non-BPS Type I : Non-BPS solutions without free fermions
In this subsection, we explore the solutions of the second order MLDE with a 1 = 0 and h NS = 1 2 . We refer such solutions to as the non-BPS solutions because they do not saturate the supersymmetric unitary bound (5.8). The central charge and weights for the solutions in this class are summarized in Table 8. The q-series of the solutions are also given in Table  14 In terms of the above characters, the torus partition functions for various spin structures become Assuming Z R (τ,τ ) = 0, we can see that the GSO projection of the fermionic CFT of our interest leads to the tricritical Ising model, The tricritical Ising model is known to preserve the N = 1 supersymmetry, but has the Ramond ground state which breaks the supersymmetry spontaneously [46]. It justifies our assumption Z R (τ,τ ) = 0. We also note that the solution with c = 7 10 simply corresponds to the characters of the first unitary N = 1 supersymmetic minimal model SM(5, 3). e 7 WZW model with level two We next move on to the solution of c = 133/10. We claim that the bosonization of the fermionic CFT results in the level-two WZW model for e 7 . To see this, let us first note that the NS, NS and R sector characters can be expressed in terms of the WZW characters: where we can see that the Ramond vacuum energy does not saturate the bound (5.8). The torus partition function for each spin structure is given by Under the assumption that Z R (τ,τ ) = 0, one obtains the modular invariant partition function of the level-two WZW model for e 7 via the GSO projection, This result suggests that the fermionic description of the level-two WZW model for e 7 preserves the supersymmetry but has a non-supersymmetric Ramond ground state. We leave this problem of verifying the emergent supersymmetry as a future work.
As a side remark, it turns out that the solutions of c = 133 10 are paired up with the characters of the first unitary N = 1 supersymmetric minimal model to produce the ( E 7,1 ) ⊗2 WZW model, in the same way that the Ising model pairs up with the e 8 WZW model at level-two to give ( E 8,1 ) ⊗2 . Precisely the bilinear relations read where χ  We propose the vanishing Z R (τ,τ ) so that the bosonization of the fermionic RCFT with c = 91/5 provides the non-diagonal modular invariant partition function of the d 7 WZW model with level three, .
We again suppose that the fermionic model of our interest has vanishing Z R (τ,τ ). As a consequence, the GSO projection leads to We also propose that the fermionic RCFT with c = 39/2 preserves the supersymemtry, but its Ramond vacuum breaks the supersymmetry spontaneously.
The solutions with c = 114 5 and 21 Let us consider the solutions with c = 114/5 and c = 21. We first make the conjecture that under the assumption of vanishing Z R (τ,τ ), the solutions with c = 114/5 describes the fermionization of the bosonic RCFT exhibiting the moonshine phenomena for 2. 2 E 6 (2).2 [16]. More precisely, the NS, NS and R sector solutions can be expressed as follows (5.108) Hereχ i (τ ) are the characters of putative bosonic RCFT with c = 114 5 and their explicit forms are given in the equation (3.55) of [16]. Let us combine the partition functions for the NS, NS and R sector, thus the GSO projection reproduce the modular invariant partition function of the hypothetical bosonic RCFT with c = 114 5 with the assumption of Z R = 0. We finally comment that the solutions of c = 21 admits the linear group 2 9 .L 3 (4) as an automorphism group. This can be shown with the help of a rank 21 lattice that constructed in [47]. We will discuss more details of this theory in an upcoming paper [28].
Non-BPS Type II : Non-BPS pairs with h NS = 1 2 Here we discuss the solutions obtained by imposing the two conditions h NS = 1 2 and a 1 = 0. There are seven solutions in this class which have non-negative integers in q-series. The profiles of the seven solutions are listed in Table 9 and their explicit q-series are presented in Table 16 and 17. The R sector conformal weights are given by h R − = c−4 8 and h R + = c 8 . Note that the N = 1 unitary bound is violated for c < 6, because h R − = c−4 8 < c 24 . Therefore the theories with c < 6 in this class cannot be unitary N = 1 supersymmetric CFTs. The unitary bound is saturated for the solution of c = 6 and it agrees with the self-dual solution in BPS type I. As discussed before, the solution with c = 6 can be understood as a hexic product of the su(2) 1 WZW model and we will not discuss it any further in this article. The NS sector characters take the analytic form of where the central charge c is constrained between 4 ≤ c ≤ 15 2 to have non-negative integer coefficients in the q-series. Because (5.110) reduces to the characters of level-one WZW model for d 4 when c = 4, we only focus on the solutions between 9 2 ≤ c ≤ 15 2 to explore the fermionic RCFT. Now it is straightforward to see that the two characters f NS 0 (τ ) and f NS 1 (τ ) satisfy the following identity, where ψ NS denote the NS partition function of a free Majorana-Weyl fermion which introduced in (4.9). This identity indicates that two independent solutions in non-BPS type II split the partition function of n = 2c tensor product of the free Majorana-Weyl fermions.
The vacuum solution f NS 0 (τ ) does not involves free fermion currents because λ ∼ 16q 1 2 thus it satisfies a condition a 1 = 0.
The S transformation of the vacuum solution ought to be written as a linear combination of f NS 0 (τ ) and f NS 1 (τ ). Therefore, f NS 0 (τ )/(ψ NS ) 2c should be a linear polynomial of λ. With the help of the analytic expression (5.110), one can read the S-transformation rule for f NS 0 (τ ) and f NS 1 (τ ), .

(5.112)
We demand 64c − 4c 2 > 0 to have a symmetric extended matrix, thus c should be smaller than 16. This constraint is automatically satisfied since we focus on the solutions between 9 2 ≤ c ≤ 15 2 . For this case, one can show that the NS sector partition function take the form of, (5.113) using the S-matrix (5.112).
We now discuss the bilinear relation which is satisfied by the NS sector solutions in this class. Let us take below three pairs of the non-BPS type II solutions, Using the analytic expression of the solutions (5.110), we arrive at the following bilinear relation. , the right hand side of (5.115) simply becomes K(τ ). Assuming this constraint is satisfied, the duality relations for three pairs (5.114) are given by (5.116) Let us turn our attention to the NS and R sector solutions. To this end, we apply T and ST transformation to (5.110). First we note that the character of a single free fermion for NS and R sector can be obtained as follows, Since the λ variable is mapped to λ λ−1 and λ−1 λ under the T and ST transformation respectively, the NS and R sector characters can be expressed as and their q-series agree with that of the Table 17 with suitable choice of a 0 and b 0 . From (5.118), the modular transformation rules for the NS and R sector characters read .

(5.119)
Because the NS sector partition function has a form of the R sector partition function is given by Non-BPS Type III : One-parameter family with c = 16 The second order MLDE exhibits infinitely many solutions of c = 16. Here we discuss a relation between these infinite solutions and the characters of the c = 16 bosonic RCFT without Kac-Moody symmetry that discussed in [48]. The c = 16 bosonic theory of our interest involves the vacuum and two primaries of weight h = 1 2 , 1. The three characters for the vacuum and primaries are known to have below q-expansion. f h=0 (q) = q −2/3 1 + 2296q 2 + 65536q 3 + 1085468q 4 + · · · , f h=1 (q) = q 1/3 1 + 136q + 4132q 2 + 67712q 3 + · · · , f h= 3 2 (q) = q 5/6 1 + 52q + 1106q 2 + 14808q 3 + · · · . (5.122) We find that the infinitely many solutions with c = 16 are one-parameter family. More precisely, the solutions of non-BPS type III can be expressed as follows.
Let us see how the characters (5.123) form Γ θ invariant partition functions as well as the consistent fusion rule algebra. To this end, it is necessary to find the S-matrices of (5.123). Note that the characters (5.122) ought to solve the third order differential equation with = 0 2 [35], 16 c 12 f (λ) and α i = h i − c 24 . The parameters ν i are expressed in terms of α i as follows, (5.125) We find that the NS, NS and R sector characters can be expressed as linear combination of the three independent solutions of (5.124) with α 0 = − 2 3 , α 1 = − 1 3 , α 2 = 1 6 . Explicitly, thus the fermionic characters (5.123) can be written in terms of the λ variable.
Using the S transformation rule of the λ variable, i.e., λ → 1 − λ, it is straightforward to read the S-matrices for the characters in the NS, NS and R sector. We find that the NS sector characters are transformed as while the transformation rules for the N S and R sector characters read .

(5.128)
There should be a symmetric extended matrix to have the consistent fusion rule algebra. If n ≥ 240, the off-diagonal component becomes 0 or negative thus one cannot find an extended matrix. In other words, the existence of the consistent fusion rule algebra constrains n < 240. This indicate that only finite number of solutions have the consistent fusion rule algebra and can be considered as the characters of certain RCFT.

Non-BPS Type IV : Single-character theories
We now focus on the class of the single-character solutions. For the solutions in this class, the vacuum solutions are allowed to have the non-negative integer coefficients. In contrast, the other solutions exhibit rational coefficients with unbounded denominators, thus we do not pay much attention to these solutions. We will show that the vacuum solutions of the non-BPS type IV are invariant under S transformation by themselves. 3 In principle, the single-character of our interest is arose from the first order MLDE with higher . In some cases, the first order MLDE can be lifted to the higher order MLDE without poles and the single-character can be identified to the vacuum solution of the higher order MLDE with = 0. To illustrate this point, we revisit the character of (e 8 ) 1 WZW model. It has been argued that the character of this theory appears as a solution with non-negative integer coefficients of the second order MLDE [35] On the other hand, the other independent solution of the above MLDE is characterized by rational coefficients in the q-expansion, hence cannot be interpreted as a character of certain RCFT. Because the vacuum solution is modular invariant by itself, (e 8 ) 1 WZW model can be understood as a single-character theory with f 0 (τ ) = j 1/3 = E 4 /η 8 .
As a single character theory, the character f 0 (τ ) has to be solution to a first order MLDE given in [8], where the Wronskian have a zero of order two. Given the equation (5.130), the following identity is true, for any meromorphic modular form H. For a generic choice of H, the above second order MLDE will be characterized by an arbitrary complicated pole structure in the Wronskian. Note that the choice of H is sensitive to the data of the second independent solution of (5.131). Once we choose H = − E 6 3E 4 , then (5.131) reproduce the second order MLDE with = 0, namely (5.129). In conclusion, the unphysical solution of (5.129) can be considered as a remnant of lifting the first order MLDE to the second order MLDE without poles. For the solutions in the non-BPS type IV class, we conjecture that a pole structure in the first order MLDE can be trivialized by lifting it to a higher order MLDE. Among them, eight solutions with c = 17 2 , 9, · · · , 12 are presented in [49]. The explicit form of the solutions in q-series are summarized in Table 18 and 19. The NS sector vacuum character can be expressed in terms of the NS partition function of a free Majorana-Weyl fermion ψ NS as follow.
Because ψ NS is invariant under the Γ θ , it is straightforward to see that f NS 0 (τ ) is a weightzero modular form of Γ θ . Therefore f NS 0 (τ ) can be interpreted as the left-moving N S sector partition function of single-character theories.
To explore the NS and R sector partition function, we take modular transformation to ψ NS . After all, the analytic structure of NS sector and R sector vacuum solutions are given by (5.134) As recently discussed in [44], the R sector constraint could be stronger than that of the NS sector. To demonstrate this aspect, we present the q-expansion of f R 0 (τ ) for c = 19 2 , c = 10 and c = 11 as the illustrative examples.

(5.135)
For c = 19 2 , the coefficients in q-series cannot be the integer values. Therefore, the partition function of c = 19 2 putative CFT cannot have the integer coefficients as well. Only the R sector solutions with c = 10 and c = 11 are allowed to have non-negative coefficients among the solutions of c < 12. Note that the unitary bound h ≥ c 24 is violated when c < 12, thus the solutions for c = 10 and c = 11 cannot be considered as the vacuum character of an unitary supersymmetric RCFT.

A SL 2 (Z) Modular Group and Forms
In this appendix, we present the definitions and properties of various modular forms that appeared in the main text. In what follows, we denote τ as a point in the Poincaré upper half plane H, and q = exp(2iπτ ).

B Modular Matrix of the Fermionic Second Order MLDE
In this appendix, we provide a closed-from expression for the NS sector solutions of the second order MLDE (5.12) and their S-matrix. Since the differential equation for f NS (λ) takes the form of Riemann's differential equation, the solutions of (B.1) can be written in terms of the Riemann's P -symbol [37] f NS (λ) = P where α ± = 1 12 6 + 2c − 24h NS ± 4 − a 1 + 2c − 32h NS + 2a 1 h NS − 4ch NS + 64(h NS ) 2 .