Thermal order in large N conformal gauge theories

In this work we explore the possibility of spontaneous breaking of global symmetries at all nonzero temperatures for conformal field theories (CFTs) in $D = 4$ space-time dimensions. We show that such a symmetry-breaking indeed occurs in certain families of non-supersymmetric large $N$ gauge theories at a planar limit. We also show that this phenomenon is accompanied by the system remaining in a persistent Brout-Englert-Higgs (BEH) phase at any temperature. These analyses are motivated by the work done in arXiv:2005.03676 where symmetry-breaking was observed in all thermal states for certain CFTs in fractional dimensions. In our case, the theories demonstrating the above features have gauge groups which are specific products of $SO(N)$ in one family and $SU(N)$ in the other. Working in a perturbative regime at the $N\rightarrow\infty$ limit, we show that the beta functions in these theories yield circles of fixed points in the space of couplings. We explicitly check this structure up to two loops and then present a proof of its survival under all loop corrections. We show that under certain conditions, an interval on this circle of fixed points demonstrates both the spontaneous breaking of a global symmetry as well as a persistent BEH phase at all nonzero temperatures. The broken global symmetry is $\mathbb{Z}_2$ in one family of theories and $U(1)$ in the other. The corresponding order parameters are expectation values of the determinants of bifundamental scalar fields in these theories. We characterize these symmetries as baryon-like symmetries in the respective models.

In this work we explore the possibility of spontaneous breaking of global symmetries at all nonzero temperatures for conformal field theories (CFTs) in D = 4 space-time dimensions. We show that such a symmetry-breaking indeed occurs in certain families of non-supersymmetric large N gauge theories at a planar limit. We also show that this phenomenon is accompanied by the system remaining in a persistent Brout-Englert-Higgs (BEH) phase at any temperature. These analyses are motivated by the work done in [1,2] where symmetry-breaking was observed in all thermal states for certain CFTs in fractional dimensions.
In our case, the theories demonstrating the above features have gauge groups which are specific products of SO(N ) in one family and SU (N ) in the other. Working in a perturbative regime at the N → ∞ limit, we show that the beta functions in these theories yield circles of fixed points in the space of couplings. We explicitly check this structure up to two loops and then present a proof of its survival under all loop corrections. We show that under certain conditions, an interval on this circle of fixed points demonstrates both the spontaneous breaking of a global symmetry as well as a persistent BEH phase at all nonzero temperatures. The broken global symmetry is Z 2 in one family of theories and U (1) in the other. The corresponding order parameters are expectation values of the determinants of bifundamental scalar fields in these theories. We characterize these symmetries as baryon-like symmetries in the respective models.

Introduction
The study of the patterns of spontaneous symmetry breaking (SSB) and their potential restoration plays a key role in understanding the phase structures of matter. From the onset of such investigations it was noted that SSB is a low temperature property. As the temperature is increased, various symmetries that were spontanously broken get restored [3][4][5][6][7][8][9]. The generality of this phenomenon was first challenged in [10], and later followed up by [11][12][13][14][15][16][17][18] as well as many explorations of interesting phenomenological implications for the physics of the early universe [19][20][21][22][23][24][25]. However, these works typically involved either UV-incomplete theories or models with imaginary or random chemical potentials. The former prevents investigating whether the symmetries are restored beyond a UV scale, while the latter requires the identification of concrete unitary theories which reproduce the results. Some holographic constructions also explored the issue of symmetry restoration at finite temperatures [26][27][28][29][30][31][32][33][34][35][36]. In such setups, non-restoration of global symmetries in QFTs translates to the violation of the no-hair theorem for black holes in the dual theories of gravity. However, in all the examples that were studied, either the relevant symmetries are restored at some critical temperature or a metastable symmetry-broken phase exists at high temperatures. So, the feasibility of a symmetry-broken phase which corresponds to the true thermal vacuum up to arbitrarily high temperatures for a unitary UV-complete theory has remained unresolved for several decades.
Recently, in [1,2] it was shown that there exist UV complete systems consisting of fundamental scalar fields in which global continuous symmetries remain consistently broken at arbitrarily high temperatures. The models examined there are Wilson-Fisher-like conformal field theories (CFTs) in (4− ) dimensions 1 symmetric under the group O(N 1 )×O(N 2 ). In the regime 0 < 1 and N 1 , N 2 1 with N 1 = N 2 , these models exhibit various patterns of spontaneous breaking of the O(N 1 ) and O(N 2 ) symmetries in thermal states. Studying CFTs enables one to extend results obtained for one temperature to all temperatures. The reader is referred to [1] for a detailed exposition of the method and the results.
In this paper we venture to study this issue of persistent symmetry breaking in (3 + 1)dimensional large N gauge theories. This positions these theories nearer to realistic particle physics systems than those considered in [1,2]. Moreover, introducing gauge symmetries allows for both asymptotically free theories as well as Banks-Zaks-like conformal field theories [38][39][40].
In either of these cases, the UV completion of the theory is guaranteed. The addition of the gauge particles also enriches the possible phase structures; in particular, it allows the system to be in the Brout-Englert-Higgs (BEH) phase. We will discuss the persistence of both the BEH phase and the associated spontaneous breaking of a global symmetry in the N → ∞ (planar/ Veneziano [41]) limit of these models. We will show that these two phenomena go hand in hand in these theories. We recall that in the standard model of cosmology the electroweak BEH phase is supposed to be modified above a certain temperature [42]. It will turn out that for obtaining the persistent phase behavior, we need to consider systems which are invariant under a gauge group G of the following form: where the G i 's are either SO(N ci ) or SU (N ci ) with two general ranks N c1 and N c2 . We will call the model where G i = SO(N ci ) the real double bifundamental model, and the one where G i = SU (N ci ) the complex double bifundamental model. As we will show in section 4.1, there is a perturbative equivalence between these two models in the planar limit 2 .
From the structure of the above-mentioned gauge group, we can see that there are two sectors (labeled by i) in these models. All the fields belonging to one of these two sectors are invariant under the gauge transformations in the other sector. The matter content in each sector comprises of N f i flavors of massless fermions transforming in the fundamental representation of G i , and N 2 ci massless scalar fields Φ i which transform in the bifundamental representation of (G i ×G i ). We will work in the Veneziano limit where the ratios N f i N ci and N c2 N c1 approach finite values as N ci , N f i → ∞, and the different couplings in the model also scale with appropriate powers of 1 N c1 and 1 N c2 . A more detailed description of these models will be provided in section 2.2.
We will show that for each sector in these models, there is a global internal symmetry which is Z 2 for G i = SO(N ci ) and U (1) for G i = SU (N ci ). We will call this a baryon symmetry in the model since the corresponding order parameter is the expectation value of a gauge invariant composite operator made of N ci number of scalar fields. It is one of these baryon symmetries that will be spontaneously broken at arbitrary nonzero temperatures in some parameter regimes.
Along with this symmetry breaking, we will find the system to exist in the BEH phase in any thermal state.
We will analyze this persistent symmetry breaking and the associated BEH phase by first determining the fixed points in the RG flow of the 't Hooft couplings in the double bifundamental models. For this, we will compute the 2-loop beta functions of the different couplings in the limit N ci → ∞. Searching for the corresponding fixed points will yield the following result: In addition to discrete fixed points, there is a fixed circle in the space of couplings for these 2-loop beta functions. We will also show that this fixed circle in fact survives in the planar limit even when all higher loop corrections to the beta functions are taken into account. 3 2 In this paper, we use the term 'planar limit' to denote the strictly N = ∞ limit. This limit of QFT has been explored in many other contexts (see e.g. [43][44][45]) which enjoy various interesting physical properties that otherwise cannot be continuously obtained from a finite but arbitrarily large N . 3 The RG flow of the different couplings will depend on the renormalization scheme [46][47][48]. In this paper, we work with dimensional regularization and the minimal subtraction scheme. We expect the existence of the conformal Such fixed circles in large N theories were explored by the authors of [49]. Their proof for the existence of exactly marginal interactions between two CFTs was based on the structure of 1-loop beta functions of the corresponding couplings. Since the double bifundamental models fall exactly in the class of systems that they considered, it is expected that there would be a fixed circle for the 1-loop beta functions of the couplings. Here, we take a step further and prove that the fixed circles in these models survive under all loop corrections to the beta functions at the planar limit. To the best of our knowledge, such a result has not been proved earlier for any four dimensional non-supersymmetric large N gauge theory. Moreover, the proof does not rely on any connection to supersymmetric or holographic setups. 4 To study the possibility of spontaneous breaking of the baryon symmetry and the persistent BEH phase for points lying on this fixed circle, we will consider the thermal effective potential of the scalar fields. We will first show that this effective potential is bounded from below for all points on the fixed circle. Next, we will determine the parameter regime in which the effective potential can have a minimum away from the origin in the field space for a subset of points on the fixed circle. Such a minimum of the potential at a nonzero field configuration implies both that the system is in a BEH phase and that the thermal expectation value of the baryon operator does not vanish. The latter leads to the spontaneous breaking of the baryon symmetry.
Assuming N c2 < N c1 , we will show that there is an upper bound on the ratio r ≡ lim for spontaneous breaking of the baryon symmetry to be even possible at any point on the fixed circle. For each value of r below this bound, we will identify the fixed points which demonstrate a symmetry-broken phase at nonzero temperatures.
Before moving on, let us pause here to take a short detour. Enlarging the scope of the theories under investigation to include gauge theories indeed adds more structure. It also brings along with it a large number of subtleties and a considerable amount of small print. Moreover, although our work deals with the spontaneous breaking of global symmetries in gauge theories without gravity, there is evidence that a theory of quantum gravity cannot have any global symmetry [51]. Thus, it is useful to summarize the nuances that arise in the study of phases due to gauge symmetries.
We will now discuss some of these nuances and explain the terminology we will be using in this work.
Compact gauge symmetries are known to express themselves in several possible phases. The Coulomb phase, the Brout-Englert-Higgs (BEH) phase and the confinement phase manifest themselves in nature and are the building blocks of the present standard model. The oblique confinement phase may have manifested itself in interesting condensed matter systems and the conformal manifolds (with the topology of a circle) in the double bifundamental models to persist under all sensible changes in the renormalization scheme. 4 We refer the reader to the discussion in [50] for the status of the search for conformal manifolds in non-SUSY theories anchored in SUSY ones. phase will also hopefully find its place in nature. There are order parameters which distinguish each of them, but the language used to characterize phase transitions corresponding to spontaneous breaking of global symmetries is not strictly appropriate for the phases associated with gauge symmetries. There is no SSB of a gauge symmetry [52]. Furthermore, it was shown [53,54] that in a class of theories containing bosons or fermions in a representation of the gauge group which is not trivial under the center of the gauge group, there are no phase transitions. The quantitative properties of the system may cross over from one appropriate effective description to another as parameters are changed, but there is no phase transition involved in the process. 5 These cases include the standard model with its bosons and its quarks, each being in the fundamental representation of the relevant group. Moreover, it was shown [53] that this also implies that in the presence of temperature there is no phase transition, only a crossover. So far, indeed, no experimental evidence has been found for such a transition. The predicted crossover behavior could be validated in future experiments [55].
We return now to explain how this pertains to the type of models we discuss in this work which contain bifundamental scalars and fundamental fermions. The constraints on the phase diagram imposed by the presence of matter in the fundamental representation do not prevent a priori a transition characterized by the spontaneous breaking of some global symmetry. As we will see, this is indeed the case in the main examples we construct. An additional feature is that we are considering the fixed points in the RG flow of these models at the planar limit. Thus, at this limit, the theories of our interest are conformal. The lack of any intrinsic scale in such conformal theories guarantees that there will not be any phase transition or crossover in them at any nonzero temperature. However, this does not preclude the possibility of a transition (or a crossover) when one switches on the temperature. Beyond this initial possible phase transition, the system persistently remains in a single phase as the temperature is increased to arbitrarily high values.
Before relating the above general feature of CFTs to the large N fixed points of our models, let us contrast this feature with the non-persistent behaviors of theories without conformal symmetry. For instance, in an asymptotically free theory with adjoint scalars, if the system is in a BEH phase at zero temperature, its non-persistent behavior (under increase in temperature) could be to first pass through a phase transition to a confinement phase and then to deconfine. It is also possible to study a non-persistent behavior by deforming a CFT. For example, consider N = 4 SUSY Yang Mills theory with a non-abelian compact gauge group of rank N c in 4 dimensions.
At the planar limit, one can infer the structure of this CFT from AdS/CFT correspondence. If the world volume is R 4 , then at T = 0 there is an exact flat potential. This makes it possible to spontaneously break the scale invariance by choosing an appropriate vacuum. For all the vacua of the theory at T = 0, we expect a phase transition from the system having O(1) degrees of freedom to O(N 2 c ) degrees of freedom as soon as a temperature is turned on 6 . Now, the above CFT can be deformed by taking the world volume to be S 3 × S 1 . Then the flat direction is lifted and the extra scale allows the dimensionless parameter T R where R is the radius of the sphere. From AdS/CFT correspondence, we know that this allows a transition from O(1) to O(N 2 c ) degrees of freedom at a finite nonzero temperature determined by T R [56].
Let us now come back to the discussion of the phases of the large N fixed points in our models. For all these fixed points, we will see that there is no flat direction of the potential at T = 0. This prevents the possibility of spontaneously breaking the scale invariance at zero temperature.
Moreover, the ground state does not demonstrate the breaking of the afore-mentioned baryon symmetries in the model. However, for a subset of these large N fixed points in certain parameter regimes, the system undergoes a transition to a phase where one of the baryon symmetries is spontaneously broken as soon as a temperature is switched on. This phase then survives up to arbitrarily high temperatures. This persistent symmetry breaking is accompanied by the Higgsing of a subset of the gauge bosons. Thus, the system can be qualitatively characterized to be in a BEH-like phase at all nonzero temperatures. It will be called in this work "being in a persistent BEH phase" or simply "being in the BEH phase". We may lapse into less rigorous terminology but we trust the reader will not be misled.
Let us end this introduction by commenting on how finite N corrections may modify the features of the double bifundamental models discussed above. In appendix I, we will study the fixed points of the beta functions in the real double bifundamental model after taking into account the subleading corrections in 1/N . Our strategy there will be to perform double expansions of the beta functions in powers of 1/N and the 't Hooft couplings. Demanding the vanishing of the first few terms in such double expansions, we will show that the degeneracy in the fixed points is lifted by the subleading corrections in 1/N , and the surviving fixed points do not demonstrate thermal order. However, we will also point out a limitation of this analysis in both appendix I.4 and section 6 which leaves open the possibility of there being isolated fixed points with thermal order even at finite N . Therefore, the results for finite N are inconclusive, and a more detailed analysis is required.

Organization of the paper:
As a precursor to all the analyses mentioned above, in section 2 we begin by discussing some familiar examples of conformal QCDs and show explicitly that they fail to exhibit spontaneous breaking of some global symmetries at nonzero temperatures. These examples help us to introduce the analytic techniques necessary for studying the possibility of SSB in a thermal 6 We expect this transition to happen even for large but finite Nc.
state. Moreover, they create the grounds for introducing the double bifundamental models as their natural extensions later in the same section.
In section 3, we first determine the perturbative fixed points of the planar beta functions in the real double bifundamental model. We demonstrate that up to the 2-loop contributions to these beta functions, there is a conformal manifold in the space of couplings which has the topology of a circle. We then show that the thermal effective potential of the scalar fields is stable at all points on this fixed circle. Later, we go on to identify the conditions under which a subset of points on this fixed circle can show spontaneous breaking of a baryon symmetry at any nonzero temperature. Moreover, we demonstrate that such a symmetry breaking is always associated with the Higgsing of a subset of gauge bosons in the model. This implies a persistent BEH phase for these fixed points.
In section 4, we study the complex double bifundamental model. We demonstrate a planar equivalence between this model and the real double bifundamental model. This allows us to extend most of the results of section 3 to this model. In particular, we find that there is a fixed circle for the 2-loop planar beta functions of this theory. Just as before, under certain conditions, a subset of points on this fixed circle demonstrates both spontaneous breaking of a baryon symmetry as well as a persistent BEH phase at all nonzero temperatures.
In section 5, we present a diagrammatic argument for the survival of the fixed circle under all loop corrections to the beta functions in the planar limit. In the process, we also prove an important feature of these planar beta functions, i.e., they are independent of the ratio r = N c2 N c1 . In section 6, we conclude by summarizing our results and commenting on the possibility of the survival of the fixed points under finite N corrections to the beta functions.
In appendix A, we prove that for all the models considered in this paper, the physical masses of the different fields are zero in the ground state.
In appendix B, we show that the Lagrangian of the bifundamental scalar QCD introduced in section 2 is indeed invariant under the Z 2 transformation that we characterize as the baryon symmetry in the model. We also show that this symmetry can be interpreted as an automorphism of the set of classes of gauge-equivalent field configurations in the theory. The analysis presented in this appendix can be generalized to show similar properties of the baryon symmetries in the double bifundamental models.
In appendix C, we present a detailed review of the two-loop beta functions in a general QCD with special emphasis on the case where the gauge group is semi-simple. We also use these general expressions to compute the beta functions of the theories that serve as precursors to the double bifundamental models.
In appendix D, we derive the two-loop beta functions of the real double bifundamental model using the techniques reviewed in the previous appendix.
In appendix E, we analyze some constraints on the unitary fixed points of the planar beta functions in the real double bifundamental model. During this analysis, we restrict our attention to the 2-loop beta functions of the gauge couplings and the 1-loop beta functions of the quartic scalar couplings.
In appendix F, we determine the minima of the thermal effective potential of the scalar fields in the real double bifundamental model for the fixed points where a baryon symmetry is broken.
In appendix G, we discuss the scaling of different planar diagrams in the double bifundamental models at the large N limit.
In appendix H, we derive some general expressions for the planar beta functions of the double trace couplings in the double bifundamental models. These expressions relate these beta functions to the corresponding wave function and vertex renormlizations.
In appendix I, we study the finite N corrections to the fixed points in the real double bifundamental model by considering double expansions of the beta functions in powers of 1 N and the 't Hooft couplings. At the end of the appendix, we also mention a limitation of this analysis.

Survey of gauge theories
In this section, we will first examine two simple variants of QCDs which contain scalars along with fermions. The gauge groups in these two models are SO(N c ) and SO(N c ) × SO(N c ) respectively. The corresponding scalar fields transform in the fundamental and the bifundamental representations of the gauge groups. 7 We will study the Banks-Zaks fixed points of these models in the following subsection. There we will explicitly show that for each of these fixed points, the minimum of the thermal effective potential of the scalar fields lies at the origin of the field space.
This implies a lack of thermal order in these CFTs. Experience in working out these models will then guide us on how to extend them to get theories with the dual phenomena of spontaneous breaking of a global symmetry and the persistent BEH phase at all temperatures. We will discuss these extensions, viz., the double bifundamental models, in section 2.2.  For similar models with SU (Nc) gauge symmetry, the Banks-Zaks fixed points were qualitatively argued to not exhibit symmetry-breaking at nonzero temperatures in [1]. 8 The actual gauge group in this model is Spin(Nc) which is a universal cover of SO(Nc). We neglect this technicality from now on to simplify the presentation. 9 The superscripts in ψ  To tame the UV divergences in the theory, we need to specify a scheme for regularization and renormalization. For all the models we introduce in this paper including the one that we are currently discussing, we choose to work in the dimensional regularization and the minimal subtraction scheme. In this scheme, we take the renormalized masses of all the fields to be zero.
As shown in appendix A, the physical masses of these fields also vanish as a result of this.
The renormalized Lagrangian 10 of this model in the above-mentioned scheme is (2.1) Here F A µν is the coefficient of the generator T A in the expansion of the field strength. We take these generators to be normalized such that the second Dynkin index of the fundamental representation of SO(N c ) is 1 2 , i.e., Tr T A T B = 1 2 δ AB . The quartic terms in the potential of the above Lagrangian can be expressed in terms of symmetric couplings λ (pqrs) abcd as follows: where λ (pqrs) abcd =4 h δ pq δ rs (δ ac δ bd + δ ad δ bc ) + δ pr δ qs (δ ab δ cd + δ ad δ bc ) + δ ps δ qr (δ ab δ cd + δ ac δ bd ) + 8 f δ pq δ rs δ ab δ cd + δ pr δ qs δ ac δ bd + δ ps δ qr δ ad δ bc .

(2.3)
Note that this model has a global O(N s ) symmetry which mixes the different flavors of scalars.
We will show that this flavor symmetry is unbroken at nonzero temperatures for the fixed points in the RG flow of the model. We will also see that associated with this is the absence of the BEH phase in all thermal states. To demonstrate these, let us work in the Veneziano limit where N c , N f , N s → ∞ , while the ratios N f,s /N c are kept fixed at finite values x f,s . In this limit, the couplings in the model scale with N c and N s as follows: where g is the gauge coupling, and λ, h and f are the 't Hooft couplings. The beta functions of these couplings are derived in appendix C.3. Here we provide the forms of the two-loop beta function of the gauge coupling and the one-loop beta functions of the quartic couplings in the Veneziano limit: (2.5) To obtain the fixed points corresponding to these beta functions, let us first set β λ = 0 which gives us For the validity of the perturbation theory at the fixed point, one must have where is a small positive number.
Solving β h = β f = 0, we get the following values of the quartic couplings at the fixed points: (2.9) Note that these upper bounds on x s are all less than 1. Therefore, for the unitary fixed points in this theory, N s < N c . In such a situation, as shown in several papers (see e.g. [58] 13 ), the 11 The reality of the couplings is necessary for the corresponding theory to be unitary. 12 Here, let us remark that identical bounds on xs were already obtained for a model of QCD with fundamental scalar equipped with the SU (Nc) gauge group in [57]. The similarity of these bounds in the two models is a direct consequence of a planar equivalence between them. This perturbative equivalence between the two models in the planar limit is similar to the one we discuss in section 4.1 13 The theories considered in this work were defined in (2 + 1) dimensions. Since the scalar fields, by definition, are in the trivial representation of the spacetime symmetry, the arguments are equally applicable to (3 + 1)-dimensional theories.
flavor symmetry is guaranteed to be preserved under certain conditions even if the scalar fields get thermal expectation values. The pattern of such expectation values would always result in a BEH phase with a residual SO(N c − N s ) gauge symmetry and no surviving Nambu-Goldstone bosons. These conditions are obeyed in the cases we discuss and anyhow, as we will show below, the scalar fields actually have zero expectation value in any thermal state. Hence, the flavor symmetry is preserved without the Higgsing of any of the gauge fields.
To demonstrate this, we will look at the quadratic terms in the thermal effective potential of the scalar fields at the temperature 1 β th . These terms take the following form: where M 2 is the thermal mass matrix of the scalar fields at the temperature 1 β th . The contribution of 1-loop diagrams to this matrix was evaluated for general (3 + 1)-dimensional gauge theories with scalar fields in [10]. For the model that we are considering currently, the form of this 1-loop thermal mass matrix is cd (S) are the generators of SO(N c ) in the representation of the r th flavor of scalar fields. The components of these generators are Substituting the expression of the symmetric coupling λ (pqrs) abcd given in (2.3) and the values of the generators given above, we get the following expression of the thermal mass matrix: where m 2 th is the thermal mass (squared) of the scalar fields, whose value is given by In the N c , N s → ∞ limit it reduces to If m 2 th is positive, then the effective potential of the scalar fields must have a minimum at the origin of the field space. In that case the expectation of the scalar fields would be zero and the flavor symmetry of the scalar fields would be preserved in a thermal state. This would also mean that the SO(N c ) gauge symmetry remains unbroken in any thermal state, and consequently, the system is never in the BEH phase. To determine whether this is indeed the case, we need to   vs. x s for the different fixed points From these graphs, we can see that indeed m 2 th > 0 at each of the fixed points for any value of x s ≤ x max s . Hence, we can conclude that the flavor symmetry is not spontaneously broken and the system is not in a BEH-like phase for these fixed points at nonzero temperatures.

Bifundamental scalar QCD
In the preceding analysis of the model of QCD with fundamental scalars, we found that unitary fixed points could exist only when N s < x max s N c where x max s < 1. We also mentioned that in such a scenario, even if the scalar fields had nonzero thermal expectation values (which they do not), the flavor symmetry would have been preserved along with the Higgsing of some of the gauge fields. It is then natural to wonder whether this bound on N s can be somehow lifted. Here we will show that it is indeed possible to lift the bound by first gauging the SO(N s ) part of the flavor symmetry and then taking N s = N c . 14 In this case, a Z 2 transformation remains as a global 14 We also add a set of Majorana fermions which transform in the fundamental representation of this SO(Nc).
symmetry. Unfortunately, we will find that that this symmetry remains unbroken in a thermal state for the Banks-Zaks-like fixed points of the model. However, the discussion of this model will serve as a prelude to the double bifundamental scalar QCDs that will be introduced in the following subsection. As we will see, some of the large N fixed points in such double bifundamental models indeed demonstrate symmetry breaking at nonzero temperatures along with a persistent BEH phase.
With this motivation, let us now introduce the bifundamental scalar QCD model as a generalization of the model of QCD with fundamental scalars that we discussed earlier. When the SO(N s ) part of the flavor symmetry in this model is gauged and N s is set equal to N c , we end up with a theory whose gauge group is SO(N c ) × SO(N c ). Now, there is an N c × N c matrix of massless real scalar fields which we will collectively denote by Φ. The components of the matrix Φ will be denoted by φ ai . These scalar fields transform in the bifundamental representation of the gauge group, thus justifying the name given to the model. In addition to the scalar fields, for each of the two SO(N c )'s there are N f flavors of massless Majorana fermions which transform in the fundamental representation of that SO(N c ) and are singlets under the other SO(N c ). We will represent these two sets of Majorana fermions by ψ i . 15 The renormalized Lagrangian of this model is given by (2.16) Here (F 1 ) ab µν and (F 2 ) ij µν are the field strengths of the gauge fields (V 1 ) ab µ and (V 2 ) ij µ which correspond to the two SO(N c )'s. These gauge fields are anti-symmetric under the exchange of color Thus, while summing over the components of (F α ) ab µν , as in the first term of the above Lagrangian, one should introduce a factor of 1 2 to count only the independent components. We take the gauge couplings for the two gauge fields to be the same (say, g) 16 . The quartic terms involving the scalar fields in the potential can be expressed in terms of symmetric couplings λ ai,bj,ck,dl as follows: 15 As before, the superscripts denote the flavor indices and subscripts denote the color indices. 16 The equality of these two couplings is guaranteed to be preserved under RG flow due to the exchange symmetry discussed below.
There is also a Z 2 symmetry which exchanges the gauge fields corresponding to the two SO(N c )'s, along with a ψ (p) ↔ χ (p) exchange and a transformation Φ → Φ T of the scalar fields. In addition, there is another Z 2 symmetry that transforms the fields as follows 17 :  17 The form of these transformations as given in (2.19) would change under a gauge transformation. As shown in appendix B, the proper way to interpret this Z2 symmetry is to think of it as an automorphism of a set of different equivalence classes of field configurations. The configurations in each class are related by gauge transformations, whereas those belonging to different classes are gauge-inequivalent. Note that the above interpretation is consistent with the existence of the gauge-invariant order parameter [det Φ] for this symmetry . 18 One could equivalently consider reflection of the elements in a column of Φ along with similar reflections of the components of χ (p) and (V2)µ. But such transformations are related to the reflection of the elements in a row by the afore-mentioned exchange of the two SO(Nc)'s. 19 Note that det Φ is a composite operator, and hence it has to be regulated and renormalized appropriately. One can work with dimensional regularization and the minimal subtraction scheme to preserve the gauge invariance of the operator. We put square brackets on the two sides of det Φ to indicate that we are talking about such a renormalized operator.
These quantum corrections are suppressed compared to the classical result det Φ due to the smallness of the couplings in the perturbative regime. Now, if Φ = 0, then it can be brought to the following diagonal form by an appropriate gauge transformation 20 : As a consequence, det Φ = 0 when Φ = 0. Hence, a nonzero expectation value of the field Φ indicates a spontaneous breaking of the Z 2 symmetry. Another important consequence of such an expectation value is that the SO(N c ) × SO(N c ) gauge symmetry would be broken down to a smaller subgroup leading to the Higgsing of some of the gauge bosons and the system existing in a BEH phase.
To explore the possibility of this dual phenomena for the fixed points in the model, let us first study the RG flow of the couplings. For this, let us work in the N c → ∞ limit. In this limit, we take the couplings in the model to scale with N c as follows: where λ, h and f are the 't Hooft couplings.
The beta functions of these couplings in the planar limit are given by (see C.3) (2.22) Repeating the steps that were followed in case of the model of QCD with fundamental scalars, first let us set β λ = 0 which gives us To ensure the validity of perturbation theory, we will work in the regime where x f → 21 4 − with 0 < 1. There are two unitary fixed points where the quartic couplings are given by with σ = ±1. 20 The argument for this is analogous to the one presented for the real double bifundamental model in appendix F.
To determine whether the Z 2 symmetry is broken or not, we can again compute the 1-loop thermal mass matrix of the scalar fields at a temperature 1 β th . It is given by where λ ai,bj,ck,dl are the symmetric couplings introduced in equations (2.18), and T 1 A (S) and T 2 A (S) are the generators of the two SO(N c )'s in the representations corresponding to the scalar fields.
The components of these generators are given by Substituting the values of the symmetric couplings and the generators, we get (2. 28) In the N c → ∞ limit, m 2 th reduces to (2.29) From the above expression it is clear that m 2 th > 0. Consequently, the thermal expectation value of the scalar field Φ is zero for both the fixed points. Therefore, the baryon symmetry remains unbroken and the persistent BEH phase is absent at all temperatures for these fixed points.

Double bifundamental models
In the previous subsection we discussed the fixed points of the RG flows of two models, viz., QCDs with additional fundamental and bifundamental scalars, where we saw that the global symmetries remain unbroken and the gauge bosons do not get Higgsed at nonzero temperatures. Now, we will consider two natural extensions 21 of these models which interestingly lead to fixed points in the large N limit where certain global symmetries are spontaneously broken and the system is in a persistent BEH phase at all nonzero temperatures. The first extension will involve taking two copies of the bifundamental scalar QCD model and introducing a quartic coupling between the scalar fields in the two copies. We will call this model the real double bifundamental 21 These extensions are partly motivated by the work [59]. model. The second extension will be to modify the gauge group by promoting the orthogonal groups to unitary groups. Moreover, the Majorana fermions will be promoted to Dirac fermions, and the real scalar fields will be promoted to complex scalar fields. We will call this model the complex double bifundamental model. Let us now discuss these models in more detail.
where the two ranks N c1 and N c2 are possibly unequal. For each of the two sectors (labeled by i), the matter content is the same as before, viz., two sets of Majorana fermions (each with Apart from the interactions between the matter fields within each sector that were already present in the bifundamental scalar model, let us now introduce a double trace interaction which couples the scalar fields in the two sectors. A schematic diagram for the matter content of this theory and the interaction between the two sectors is provided in figure 2. The renormalized Lagrangian of the model is given by where (F iα ) µν is the field strength corresponding to the gauge field (V iα ) µ .
We will be studying this model in the Veneziano limit (N c1 , N c2 → ∞) where N c1 → r, and the different couplings scale as follows: Here g i is the gauge coupling of the i th sector and λ i , h i , f i and ζ are the 't Hooft couplings. The symmetries in the model are essentially the same as those for the bifundamental model, except now we have one set of symmetries for each sector. We will be interested mainly in the Z 2 baryon symmetries for the two sectors. The action of the baryon symmetry in the i th sector on the fields in that sector is as follows 22 : where T i is an N ci × N ci diagonal matrix of the following form: Here, the transformation has a nontrivial action only on the first row of the scalar field Φ i , the first component of the fermion ψ the bifundamental scalar QCD with a = 1. As earlier, we can choose [det Φ i ] to be an order parameter for this symmetry. A nonzero thermal expectation value of Φ i would imply that this order parameter is nonzero and consequently, the baryon symmetry is broken in a thermal state.
As earlier, this would also mean that some of the gauge bosons in the i th sector are Higgsed leading to the system being in a persistent BEH phase.
In section 3, we will undertake a study of the fixed points in this model to see whether these dual phenomena actually occur in a thermal state. There we will find that this model has two very interesting properties in the Veneziano limit: 1. Considering up to two loop contributions to the planar beta functions, there is a conformal manifold. It is topologically a circle in the space of the double trace couplings. Later in section 5, we will show that this conformal manifold survives even when all higher loop corrections to the planar beta functions are taken into account.
2. When the ratio of the ranks N c1 and N c2 is sufficiently away from 1, the baryon symmetry corresponding to the sector with the smaller rank is spontaneously broken in a thermal state for a subset of points on this conformal manifold. Moreover, at these fixed points, the gauge symmetry in this sector is broken down from SO(N ci ) × SO(N ci ) to just SO(N ci ) in any thermal state. This leads to Higgsing of half of the gauge bosons in this sector, and thus the system is in a persistent BEH phase.
We will defer further discussion of these features to section 3, and next introduce a closely related model which shares many of the properties mentioned above.

Introduction to the complex double bifundamental model
In the preceding introduction to the real double bifundamental model we saw that the gauge group in that model has the following structure: where G i = SO(N ci ). In addition, there is a discrete (Z 2 ) global symmetry for each sector which we claimed to be spontaneously broken in a thermal state at a subset of fixed points. This should naturally lead to the following question: Is there any model where a continuous global symmetry is similary broken at arbitrarily high temperatures?
To address this question, we will now extend the above-mentioned model by taking the group G i to be the special unitary group SU (N ci ). We will further promote the Majorana fermions to Dirac fermions transforming in the fundamental representation of SU (N ci ), and the real scalar fields to complex scalar fields transforming in the bifundamtal representation of SU (N ci ) × SU (N ci ). 23 The schematic diagram for the matter content in this model is still the one shown in figure 2. The renormalized Lagrangian of the model is given by The components (F iα ) A µν are the coefficients of the generators of SU (N ci ) in the expansion of the field strength (F iα ) µν corresponding to the gauge fields (V iα ) µ . As before, these generators are normalised such that the corresponding second Dynkin index is 1 2 . Analogous to the previous model, we take the gauge couplings for the two SU (N ci )'s in the i th sector to be the same (say, g i ). Therefore, there is a Z 2 symmetry for each sector (labeled by i) which exchanges the two sets of gauge fields in that sector, along with the exchange of the ψ Moreover, there is a symmetry under charge conjugation of each of these fermions. But the symmetries that we will mainly be interested in are two U (1) baryon symmetries, one for each sector, which transform the fields in the respective sectors as follows 24 : where, for each θ ∈ [0, 2π), (T i ) θ is an N ci × N ci diagonal matrix of the following form: (2.39) 24 Let us note here that the transformations given in (2.38) can be reduced to the following simpler form by an appropriate global gauge transformation: (2.37) Furthermore, the above transformation of the fermionic fields can be undone by a corresponding U (1)V transformation. This would leave us with a U (1)/Nci transformation acting just on the scalar fields. We choose to work with the transformation given in (2.38) rather than this simpler transformation to keep the U (1) symmetry manifest in the Nci → ∞ limit. As earlier, this U (1) symmetry can also be interpreted as an automorphism of the set of classes of gauge-equivalent field configurations.
This transformation leaves the remaining fields in the i th sector, as well as all the fields in the other sector, unchanged. We can again take [det Φ i ] to be an order parameter for the baryon symmetry in the i th sector. As before, a nonzero expectation value of Φ i in a thermal states implies a nonzero value of this order parameter resulting in the spontaneous breaking of the corresponding baryon symmetry and a persistent BEH phase.
To study the possibility of the above phenomena, we will look at the fixed points of this model in the Veneziano limit. As mentioned earlier, in this limit N c1 , N c2 → ∞, N c1 → r, while the quartic couplings scale exactly as given in (2.32).
The analysis of this model will be simplified by the fact that it is dual to the real double bifundamental model in the Veneziano limit. We will discuss this planar equivalence between the two models in section 4.1. As a consequence of this equivalence, there will be a conformal manifold (a fixed circle) for this model in the planar limit. Moreover, when the two ranks N c1 and N c2 are sufficiently different, we will see that the baryon symmetry corresponding to the smaller sector is spontaneously broken in a thermal state for some of the fixed points on the manifold. For these fixed points, the gauge symmetry in the smaller sector is broken down from SU (N ci ) × SU (N ci ) to SU (N ci ) in any thermal state, which leads to the Higgsing of half of the gauge bosons in this sector. Thus the system remains in a persistent BEH phase at all temperatures for these fixed points. In section 4, we will discuss these features of the complex double bifundamental model in detail.

Real double bifundamental model
In this section we will restrict our attention to the real double bifundamental model. First, in the following subsection, we will determine the weakly coupled fixed points in this model by employing perturbation theory in the planar (Veneziano) limit that was mentioned in the previous section. We will find that in this limit, the beta functions of the different couplings yield nontrivial fixed points which are analogous to the Banks-Zaks fixed points in more familiar QCDs.
We will show that the planar 2-loop beta functions of the gauge couplings do not depend on the quartic 't Hooft couplings. Hence, the values of the gauge couplings at the fixed points can be determined independently. We will work in a regime where these gauge couplings are small.
For this, we will see that the ratios x f i ≡ N f i N ci need to be tuned as follows: where 0 < i 1. Choosing the gauge couplings to be small, in turn, will allow us to determine the fixed points of the beta functions of the quartic couplings perturbatively. By studying the 1-loop contributions to these beta functions, we will find that instead of a discrete set of fixed points, there is a conformal manifold (topologically, a circle) in the space of couplings. Moreover, we will demonstrate that this degeneracy in the fixed points survives even after the contributions of 2-loop diagrams to the beta functions are taken into account.
Later, in subsection 3.2, we will analyze the thermal effective potential of the scalar fields at the above-mentioned large N fixed points. We will demonstrate that this effective potential is stable at all points on the fixed circle. We will also show that when the ratio r ≡ N c2 N c1 is below a certain bound, some of the points on the fixed circle demonstrate spontaneous breaking of the baryon symmetry in the second sector at all nonzero temperatures. Along with this, we will see that at these fixed points, the gauge symmetry in the second sector gets broken down to SO(N c2 ) in any thermal state, and the system exists in a persistent BEH phase.

Fixed points in the planar limit
The beta functions of the couplings in the real double bifundamental model are determined up to the contributions of 2-loop diagrams in the appendix D. Here, we provide their expressions at the leading order in the planar limit.
• Beta functions of the gauge couplings (up to 2-loops): • 1-loop beta functions of the quartic couplings: • 2-loop corrections to the beta functions of the quartic couplings: Note that in the expression of β 2-loop f i given above, one of the terms has the coupling λ i .
Here, i denotes the complement of i, i.e., i = 2 when i = 1, and i = 1 when i = 2.

Fixed points of the 2-loop beta functions of the gauge couplings and 1-loop beta functions of the quartic couplings
We will now determine the fixed points of the RG flow of the couplings where the two sectors are coupled to each other, i.e., ζ = 0. To determine these fixed points, let us first set β λ i = 0. A nontrivial solution of this equation is As mentioned earlier, we can see that gauge couplings at the fixed point become small while remaining positive when x f 1 and x f 2 approach the value 21 4 from below. Next, we set β 1-loop Now, for convenience, let us define the following combinations of the quartic couplings: In terms of these couplings, the 1-loop beta functions are as follows In appendix E, we show that a unitary fixed point of the above beta functions with ζ = 0 can exist only if x f 1 = x f 2 ≡ x f which leads to the equality of the two gauge couplings: Moreover, as proved in the same appendix, at such a unitary fixed point, we must have 10) or equivalently, Plugging these solutions into the beta functions we get the following simpler expressions: (3.12) Now it turns out that for the fixed points with ζ = 0, β 1-loop Here we have substituted h p by its value given above to obtain f p in terms of λ. Substituting these values of f p and h p in β 1-loop fp , we get (3.14) Finally setting β 1-loop fp = 0, we get Therefore, we see that there is a circle of fixed points for the 1-loop beta functions of the quartic couplings. In the following subsection, we will show that this fixed circle survives even when the 2-loop contributions to the beta functions of the quartic couplings are taken into account.

2-loop corrections to the fixed points of the quartic couplings
To analyse the corrections to the fixed points after including the 2-loop contributions to the beta functions of the quartic couplings, we split these couplings as follows: where h 0 , f 0i and ζ 0 are solutions to the fixed point criteria arising from the 1-loop beta functions.
These quantites are O(λ), and as derived in the previous subsection, they satisfy the following conditions: The quantities δh i , δf i and δζ are corrections to these fixed point values due to 2-loop contributions. They are O(λ 2 ). Now, retaining terms up to O(λ 3 ) in the 1-loop beta functions, we have To analyze the solutions of equations (3.20) and (3.21), let us introduce the following com-binations of the couplings: (3.23) After substituting h 0 , f 0p , δh and (ζ 2 0 + f 2 0m ) by their values, equations (3.20) and (3.21) lead to the following equations: Note that the last two equations can be satisfied by setting Therefore, we see that the 2-loop contributions to the beta functions do not impose any additional constraint on ζ 0 and f 0m . From this, we can conclude that the fixed circle survives even after accounting for the 2-loop terms in the planar beta functions. Moreover, note that the LHS of (3.24) is proportional to δ(ζ 2 + f 2 m ). This means that this equation gives the 2-loop correction to the radius of the fixed circle. In section 5, we will prove that this fixed circle actually survives under all higher loop corrections in the planar limit.

Analysis of symmetry breaking at finite temperature
In this subsection, we will analyze the pattern of symmetry breaking at finite temperature for the large N fixed points in our model. We will begin by demonstrating that the effective potential at any finite temperature is stable for all points on the fixed circle (at least up to leading order in λ). Then we will compute the thermal masses of the fields Φ 1 and Φ 2 and identify the conditions under which some of the fixed points demonstrate the dual phenomena of spontaneous breaking of a baryon symmetry and a persistent BEH phase at nonzero temperatures. In what follows, we will truncate the values of the different couplings on the fixed circle to their leading order.

Stability of the effective potential
The leading order terms in the thermal effective potential of the fields Φ 1 and Φ 2 are O(λ).
At this order, there are only terms which are quadratic and quartic in the fields. The quadratic terms lead to thermal masses of these fields which will be discussed in the following subsection.
Here we will consider only the quartic terms which would determine whether the potential is stable (i.e. bounded from below) or not. 25 These terms are of the following form: where ρ 1 and ρ 2 are symmetric matrices defined as follows: The first two terms in (3.28) are manifestly non-negative if h 1 > 0 and h 2 > 0. A sufficient condition for the remaining terms to be non-negative is After rescaling the couplings by appropriate powers of N c1 and N c2 , these conditions can be summarized as follows: Let us now check whether these conditions are satisfied for the points on the fixed circle. At these points, we have and . From the values of h 1 and h 2 we can conclude that they are positive since 3− To see that the bounds on f 1 and f 2 are satisfied note that (3.33) imposes the following 25 The coefficients of these quartic terms in the classical potential are O(λ). The quantum corrections to this classical potential are suppressed by higher powers of λ. constraint on the value of f m : From this we can conclude that Finally, note that Therefore, we see that all the conditions sufficient for the non-negativity of the quartic terms in the potential are satisfied at all points on the fixed circle. In fact, we have proved a slightly stronger condition, i.e., these quartic terms are strictly positive-definite for nonzero values of the field configurations. This means that as long as the thermal masses are non-negative, the potential increases steadily when one moves away from the origin along any direction in the field space.
This is true even when either of the thermal masses is zero ruling out the possibility of there being any flat direction 26 in such a scenario.

Thermal masses
Now, let us discuss the quadratic terms in the potential. These terms have the following form: As earlier, we can evaluate the contribution of 1-loop Feynman diagrams to the thermal mass matrix M 2 by using the general expression derived in [10]. For the model that we are presently considering, this expression is given by Here β −1 th is the temperature. λ ai,bj,ck,dl are the symmetric couplings introduced in equations (D.6), (D.7) and (D.8). Using the explicit forms of the symmetric couplings and the generators, one can show that the matrix M 2 takes the following form: where m 2 th,1 and m 2 th,2 are the 1-loop thermal masses (squared) of the 2 fields with the following values: If either m 2 th,1 or m 2 th,2 is negative, then the minimum of the potential will be away from the origin in the field space, and hence at least one of the order parameter will be nonzero. This would indicate a spontaneous breaking of the baryon symmetry and the Higgsing of some of the gauge bosons in the respective sector.
Thus, to determine whether these dual phenomena actually occur at any of the points on the large N fixed circle, we need to compute the thermal masses at these points. The values of the couplings at these fixed points are as follows: In addition, the double trace couplings ζ and f m ≡ f 1 −f 2 2 are constrained to lie on the circle Substituting the values of the above couplings in the expressions of the thermal masses and taking the limit N ci → ∞, we get where the parameter r is defined as follows: This parameter did not enter in the 2-loop planar beta functions of the couplings 27 . However, as we can see now, it leaves its imprint in the effective potential of the scalar fields through the 27 In section 5, we will show that this feature of the planar beta functions holds at all orders of the 't Hooft couplings.
thermal masses. There are constraints on the value of this parameter which need to be satisfied to get a symmetry-broken phase at non-zero temperatures. Moreover, even when the constraints on r are satisfied, only a subset of points on the fixed circle demonstrate a symmetry-broken phase. We will next discuss these conditions on r and the fixed points. In what follows, we will assume that r < 1, i.e., N c1 > N c2 . The results that we will get can be generalized to the case r > 1, i.e., N c1 < N c2 , by a (1 ↔ 2) exchange of indices everywhere.

Conditions for symmetry breaking
To analyze the conditions on the parameter r and the fixed points for breaking of the baryon symmetry and a persistent BEH phase at nonzero temperatures, let us first consider the following relation between the couplings f m and ζ at any point on the fixed circle: This relation imposes the following bounds on the values ζ at these fixed points: To avoid carrying along the factor of λ throughout our analysis, let us define a normalized version (z) of the coupling ζ by stripping off the factor of λ: The above-mentioned bounds on ζ then lead to the following bounds on the normalized coupling z: Now, by using equations (3.42) and (3.44), we can express the thermal masses of the scalar fields in terms of the parameter r and the normalised coupling z as follows: (3.48) Note that the above expression of m 2 th,1 is positive-definite as can be seen from the following inequalities: Here, in the first line, we have used the relation |a+b| ≤ |a|+|b|, and the fact that . In the second line, we have imposed the upperbound on |z| given in (3.47). In the last line, we have used the fact that r < 1.

Now, let us see under what conditions m 2
th,2 can be negative. First note that when z = 0, m 2 th,2 is positive: We can also see that for any particular value of r, m 2 th,2 in either of the branches in (3.48) is a continuous function of z. Therefore, for any fixed r, if m 2 th,2 has to be negative at some value of z in the range (− ), then it must also pass through zero at some point (say, z 0 ) in this interval. At this point we have For the existence of a real solution of this equation, the discriminant of the corresponding quadratic polynomial must be non-negative, i.e., 3 2r This puts an upper bound on the value of r for which a symmetry-broken phase can exist for any of the points on the fixed circle. Now, for every value of r lying below this upper bound, there are two solutions of the quadratic equation given above: 2r . Therefore, both these solutions are negative. Hence, the value of z (or equivalently, ζ) for which m 2 th,2 < 0 is always negative 28 .
Notice that when we went from the first line of equation ( In this case, we have (3.54) Note that both the terms in the numerator of the RHS of the above equation are non-negative.
Therefore, taking the positive square roots of both sides of this equation, we get 28 This can also be seen more directly from (3.42) since fm ≤ √ 18 √ 6−39 16 λ < 3 4 λ for all points on the fixed circle.
From the above equation we get Therefore, we find that this solution always belongs to the branch where f m > 0.
In this case, we have ( 3.57) Now, the second term in the numerator of the RHS of the above equation is negative. So, while taking the square root of both sides of the equation, one has to be careful about the sign of the quantity within the brackets. To check this sign, let us look at the magnitude of the expression within the root in the RHS. This expression can be re-cast as follows: Therefore this expression will have a magnitude less than (or equal to) (12r 2 ) 2 (leading to a positive sign of the quantity in brackets in the RHS of (3.57)) only if Thus, for and we can again conclude that the solution lies in the branch where f m > 0.
For r < and the solution lies in the branch where f m < 0.
Let us summarize the above results: • When , there are two two fixed points on the same branch (f m > 0) for which the thermal mass (squared) m 2 th,2 vanishes. The values of the normalized coupling z at these two fixed points are: , again we have two fixed points for which m 2 th,2 vanishes. The values of z at these fixed points are again given by the expressions in (3.62). However, the fixed Now that we have identified the fixed points where m 2 th,2 =0, it is possible to determine the points at which it is negative. To apply the argument of continuity of m 2 th,2 as we move along the fixed circle, we need to work in a particular branch. If there are two fixed points where m 2 th,2 vanishes in that branch, then in the interval between these two points it will be negative. On the other hand, if there is only one fixed point (say, at z = z 0 ) in that branch where m 2 th,2 vanishes, then it will be negative for all values of z below z 0 .
Before enumerating all the fixed points where m 2 th,2 < 0, let us briefly comment on the points where the afore-mentioned branches meet. At these points, we have f m = 0 and z = σ 39 16 where σ = ±1. From (3.48), we find that at these points: . , and m 2 th,2 < 0 when r < Thus, from the above analysis we can finally conclude that m 2 th,2 < 0 for the fixed points satisfying any one of the following conditions: ), ζ ∈ (z min λ, z (2) 0 λ), and f m < 0.
Here z (1) 0 and z (2) 0 are the values of the normalized coupling z defined in (3.62), and z min is the lower bound on the value of z which was given in (3.47), i.e., .
We provide graphical plots of these fixed points in figure 3. In appendix F we have shown that for all these fixed points, the thermal expectation values of the scalar fields (up to gauge transformations) have the following forms: where 'diag' stands for a diagonal matrix with the corresponding diagonal entries given in the brackets. Note that in the zero temperature limit, Φ 2 = 0 due to the vanishing of m 2 th,2 . This means that both the baryon symmetries are unbroken in the ground state.
For any nonzero temperature, the above forms of the thermal expectation values of the scalar fields lead to the following consequences: This relates O 2 to O 1 as follows: Thus, we see that the residual gauge symmetry is just SO(N c2 ). Consequently, half of the generators in the second sector are broken. This, in turn, leads to half of the gauge bosons in this sector getting Higgsed, and the system being in the Brout-Englert-Higgs (BEH) phase at all nonzero temperatures.
Therefore, we find that the spontaneous breaking of the baryon symmetry and the persistent BEH phase always occur together in this model as mentioned in section 1.

Complex double bifundamental model
In this section we will turn our attention to the complex double bifundamental model. We will begin by demonstrating that this model is actually dual to the real double bifundamental model in the Veneziano limit. Therefore, in this limit, all the conclusions that we arrived at in the last section are equally valid for the complex double bifundamental model. In particular, we will see that for the planar beta functions in this model, there is a fixed circle in the space of couplings. Just as before, by looking at the thermal masses of the scalar fields we will determine the conditions under which a subset of these fixed points demonstrate spontaneous breaking of the baryon symmetry and a persistent BEH phase at nonzero temperatures.

Planar equivalence between the two double bifundamental models
Let us now demonstrate that the real and the complex double bifundamental models are perturbatively equivalent in the Veneziano limt. This equivalence between the two models rests on the procedure of orbifolding a parent theory to obtain an equivalent daughter theory. 29 In this case, a real double bifundamental model symmetric under the gauge group and with 2N f i flavors of Majorana fermions in the i th sector serves as the parent theory.
On the other hand, the daughter theory is a complex double bifundamental model symmetric and with N f i flavors of Dirac fermions in the i th sector.
The idea of orbifolding a parent theory is to first identify a discrete global symmetry in it.
One can then project all its fields to the components that are invariant under this symmetry group. Operators composed of these projected fields are also invariant under the discrete global symmetry and hence, are known as neutral operators. One can next project the Lagrangian of the parent theory to its part containing such neutral operators by simply removing the terms that are not invariant under the discrete symmetry. The form of this truncated Lagrangian of the parent theory should be similar to the Lagrangian of the daughter theory. Based on this similarity between the two Lagrangians, it can be shown that in the planar limit there is a mapping between the couplings in the two theories under which the correlators of neutral operators in the parent theory are equal to the corresponding correlators in the daughter theory. In this precise sense, the two theories are dual to each other in the planar limit. We will next show how this procedure of orbifolding of the real double bifundamental model leads to its planar equivalence with the complex double bifundamental model. Our arguments for this equivalence will closely follow those given for similar dualities in [64,65]. 30 To demonstrate the planar equivalence between the two models, let us first express the Majorana fermions in the real double bifundamental model in terms of Weyl spinors as follows 31 : 29 This is a counterpart of string orbifold equivalence discussed in [62,63]. 30 See also [66][67][68][69][70][71][72][73] for discussions on the various aspects of the planar equivalence. 31 From here onwards, we shall be working in the Weyl representation of the Clifford algebra.
where the spinors with the subscript 'L' are left-handed, while those with the subscript 'R' are right-handed. Substituting the above expressions of the Majorana fermions in equation (2.31), we get the following form of the Lagrangian (after removing some total derivatives): where σ µ = (I, − − → σ ), and the dummy index q runs over the values {1, · · · , N f i } for the fermions in the i th sector. We have put a superscript R on the couplings to distinguish them from the couplings in the dual complex double bifundamental model. Now note that the above Lagrangian is invariant under the following Z 2 transformation of the fields 32 : where To project the fields on to their neutral components, we demand the invariance of the fields under 32 Here, the fields in both the sectors are transformed simultaneously.
the above transformations. This yields the following structure of the projected fields: Here (B A iα ) µ and (C S iα ) µ are real (N ci × N ci )-matrices which are anti-symmetric and symmetric respectively. Similarly, Φ i1 and Φ i2 are real (N ci × N ci )-matrices, but they do not have any symmetry property under transposition.
Let us now combine these projected fields to define a new set of fields: Here, ( V iα ) µ is an N ci × N ci Hermitian matrix. Therefore it can serve as a generator of U (N ci ). Note that this is distinct from a generator of the group SU (N ci ) in the complex double bifundamental model (the daughter theory). However, the difference in the number of generators between the two cases is suppressed by a factor of 1 N 2 ci compared to the total number of generators in either of them. Hence, this difference can be ignored while computing correlators in the planar limit [69]. With this caveat, let us go on treating ( V iα ) µ as a gauge field in the projected theory.
The field strength corresponding to this gauge field is given by where g R i is the gauge coupling in the i th sector of the parent theory. The Weyl fermions introduced above can be combined together to form Dirac fermions as shown below: Now writing the Lagrangian of the parent theory in terms of the gauge fields ( V iα ) µ , the complex scalar fields Φ i , and the Dirac fermions ( Ψ where Finally, to obtain the Lagrangian of the daughter theory, we need to rescale the above Lagrangian by a factor of Γ −1 , where Γ is the order of the discrete global symmetry group in the parent theory. This is necessary to account for contribution of projectors on the fields to the correlators of neutral operators in the parent theory [64,65]. In our case, Γ = 2 since the global symmetry group is Z 2 . To express this rescaled Lagrangian in terms of canonically normalized fields, let us rescale the fields appearing in (4.10) as follows: In terms of these rescaled fields, the Lagrangian of the daughter theory takes the following form: (4.14) Note that if we expand the field strengths ( F iα ) µν in terms of generators T A iα normalized so that the corresponding second Dynkin index is 1 2 , then the Lagrangian given above has the same structure as the Lagrangian of the complex double bifundamental model in (2.36). In fact, the two Lagrangians are exactly identical if we impose the following relations between the couplings: where we have now introduced the superscript 'C' to indicate that the couplings correspond to the complex double bifundamental model. From (4.15), one can derive the relations between the 't Hooft couplings 33 in the two dual theories to be As a consistency check, we provide the forms of the 2-loop beta functions of the gauge couplings and the 1-loop beta functions of the quartic couplings in the complex double bifundamental model at the Veneziano limit: From the planar equivalence between the two double bifundamental models, we can obtain the unitary fixed points (with ζ C = 0) corresponding to the beta functions given above. As earlier, for the existence of such unitary fixed points, the ratios When this condition is satisfied, the unitary fixed points are all the points lying on the following manifold in the coupling space: . As before, this conformal manifold has the topology of a circle. The orbifold equivalence with the real double bifundamental model ensures that this fixed circle would survive even after taking the 2-loop contributions to the planar beta functions into account. In fact, in section 5, we will prove that this fixed circle survives under all loop corrections at the planar limit.

Analysis of thermal effective potential and symmetry breaking
Let us now discuss the thermal effective potential of the scalar fields in the complex double bifundamental model for the fixed points given in (4.18). From here onwards, we will drop the superscript 'C' over the couplings.
First we will check the stability of the potential at these fixed points. The quartic terms in the potential (up to leading order in λ) are As in the case of the real double bifundamental model, the single trace interactions are manifestly positive because To ensure that the other terms are positive, we need to check the positivity of f 1 , f 2 and ( f 1 f 2 − ζ 2 ). But this is guaranteed by the relations in equation (4.16), and the positivity of the corresponding quantities in the dual real double bifundamental model. 34 Therefore, we can conclude that the potential is stable for all points on the fixed circle. 34 See the proof of this in section 3.2.1.
Next, consider the quadratic terms in the potential which have the following form: where the thermal masses (squared) are 35 (4.22) In the Veneziano limit (N ci → ∞), we can drop all the subleading terms in 1 N ci and substitute the couplings by their values at the fixed points (given in (4.18)) to obtain N c1 . Note that these thermal masses are equal to the corresponding thermal masses in the dual real double bifundamental model (see (3.42)) under the mapping between the couplings in the two models given in (4.16). This is actually a consequence of the perturbative orbifold equivalence between the two models being valid in thermal states. Since the proof of the perturbative planar (orbifold) equivalence [64] between the parent and the daughter theories relies only on some combinatorics arising from the projection under a discrete automorphism, it can be extended to the finite temperature case without any major modification. 36 This ensures the equality of the 1-loop planar diagrams [10,5] contributing to the thermal masses in the two theories.
The equality of the thermal masses in the two dual theories allows us to just use the results obtained in section 3.2 to investigate the conditions under which one of the thermal masses (squared) is negative. We quote these conditions below.
First, without any loss of generality, we assume that N c2 < N c1 , i.e., r < 1. Under this assumption, on one hand, m 2 th,1 > 0 for all points on the fixed circle. On the other hand, m 2 th,2 < 0 for the fixed points satisfying any one of the following conditions: 0 λ), and f m > 0, 35 These thermal masses are computed by methods analogous to those discussed in section 3.2.2. 36 It would be interesting to see if the equivalence holds non-perturbatively for the double bifundamental models. See [74][75][76][77] for the criteria for such non-perturbative planar equivalences between various models. are the values defined in (3.62), and z min is given in (3.64). For the reader's convenience, we provide these values once more below: . (4.24) The graphical plots of these fixed points are similar to those given in figure 3. The only difference is that the quantity z in those figures has to be now defined as z = ζ 4λ . For all these fixed points, the thermal expectation values of the scalar fields (up to gauge transformations) are as follows 37 : where θ ∈ [0, 2π). Note that in the zero temperature limit, we again have Φ Thus, only one of the two unitary transformations is independent and the residual gauge symmetry is SU (N c2 ). As a result of this, half of the generators in the second sector are broken, and accordingly, half of the gauge bosons in this sector get Higgsed. Therefore, we find that just as in 37 The arguments for arriving at these expectation values are exactly analogous to those in appendix F for the real double bifundamental model. To demonstrate the above result, we will first prove a lemma in the following subsection that the planar beta functions in these models are independent of the ratio r = N c2 N c1 . This ratio is the only quantity that could have led to an asymmetry in the forms of the planar beta functions of the couplings in the two sectors. 38 The above-mentioned lemma ensures the absence of such an asymmetry in the planar limit.
While proving this lemma, we will illustrate the planar diagrams that contribute to connected correlators in the ground state of these models. To extract the planar beta functions from such correlators we work in dimensional regularization 39 and the minimal subtraction scheme. In this scheme, the contributions to the beta functions come from the coefficients of the O( 1 ε ) terms in such connected correlators, i.e., the residues of these correlators corresponding to the pole at ε = 0.
One advantage of working with dimensional regularization is that all diagrams with a massless tadpole 40 vanish in this scheme [78]. The vanishing of such diagrams is not essential to the arguments that we will propose as such tadpoles can at most contribute to the renormalization of the masses 41 and they do not affect the beta functions. However, it simplifies the analysis slightly since we have taken the renormalized masses of all the fields in the double bifundamental models to be zero. As we have argued in appendix A, this also ensures the vanishing of the N c2 could also have led to an asymmetry between the two sectors. But we have already shown that in such a scenario, there is no unitary fixed point with ζ = 0. Therefore, we keep assuming that x f 1 = x f 2 in the Nc1, Nc2 → ∞ limit. 39 In dimensional regularization, we take the theories to be defined in (4 − ε) dimensions. 40 Here, by tadpole, we mean a subdiagram which has only one external vertex. It doesn't matter how this vertex is connected to the rest of the diagram. The loop integrals in such a subdiagram vanish in dimensional regularization when all the propagators in such loops are massless. 41 In fact, these tadpoles are responsible for the generation of the thermal masses of the scalar fields. In a thermal state, the vanishing of such tadpoles is spoiled by the contributions of the Matsubara modes.
corresponding physical masses in the ground state which is a scheme-independent statement and a necessary condition for the theory to be conformal at the fixed points of the beta functions.
Before proceeding further, let us make a brief comment about our convention for representing different diagrams in this section. The fields in the double bifundamental models are matrices, and hence an appropriate convention for drawing the corresponding Feynman diagrams would be 't Hooft's double line notation. Indeed, we make use of this notation while discussing the large N scaling of different diagrams in appendix G. However, in this section, we will represent the propagators of the scalar fields by single lines to avoid unnecessary clutter. We hope the reader will not be confused by this slightly unconventional notation. Now, without further ado, let us delve into the proof of the afore-mentioned lemma. As we will see in section 5.2, the insights gained while proving this lemma will be essential to the proof of the survival of the fixed circle in the planar limit. In sections 3 and 4, we found that the two-loop planar beta functions in the double bifundamental models are independent of the ratio r = N c2 N c1 . Here, we will show that this feature of the planar beta functions actually persists at all orders in the 't Hooft couplings.
Before proceeding with the argument, we note that the interactions in the double bifundamental models can be divided into two classes: single trace and double trace interactions. The couplings corresponding to the former are the gauge couplings (g 1 and g 2 ) 42 and the quartic couplings h 1 and h 2 , whereas those corresponding to the latter are the quartic couplings f 1 , f 2 and ζ.
Now, let us consider the double bifundamental model as two initially decoupled bifundamental scalar QCDs coupled through the double trace coupling ζ. Note that ζ is the only coupling that mediates an interaction between the two bifundamental scalar QCDs. Hence, if the theory starts with ζ = 0 at any energy scale, it will be preserved along the RG flow. In this case, as we saw already, each bifundamental scalar QCD's large N planar beta function would be determined solely by the corresponding 't Hooft couplings (λ Here, a question naturally arises: what happens to the structure of the beta functions after a non-zero double trace interaction ζ between the two sectors is turned on? Since each decoupled sector could have an independent large N limit (i.e. N c1 = N c2 ), one might wonder whether the planar beta functions of the couplings would depend on the additional parameter r = N c2 N c1 . We will now argue that this is not the case when the 't Hooft coupling ζ is defined with an appropriate normalization as ζ ∝ N c1 N c2 ζ.
To prove the r-independence of the planar beta functions, we will use the following feature of the connected planar diagrams that contribute to such beta functions: Any double trace vertex in a such a diagram links two otherwise disconnected planar subdiagrams. 43   First, consider the single trace couplings λ i and h i . We will show that, just like the wave function 43 See e.g. [79] which discusses general large N scaling involving multi-trace vertices. 44 We note that there is a scheme independent way to understand the reason why tadpole diagrams do not contribute to the wavefunction renormalization from their external momentum independence. renormalizations, the vertex renormalzations of these single trace couplings are independent of the double trace couplings in the planar limit. This will, in turn, imply that the planar beta functions of the single trace couplings are totally independent of the double trace couplings 45 . To prove this, let us work with a specific example, viz., the planar vertex normalization of the coupling h 1 . The structure of this vertex normalization can be understood from connected diagrams with 4 external legs of the scalar field Φ 1 . We take the color indices corresponding to these external legs to be (ai, bj, aj, bi) with a = b and i = j. Now, for these external legs, if there is a connected planar diagram containing a double trace vertex, then by the arguments given above we can split it into two disconnected pieces by removing this vertex and joining the residual legs with and ζ ∝ ζ N c1 N c2 , the replacement of the two vertices leads to an additional factor of ( N c1 N c2 ) 2 in the diagram shown in figure 6b compared to the one in figure  6a. However, this factor is exactly cancelled by the replacement of the B I 1 blob by the B I 2 blob. Therefore, we can conclude that the scaling of the two diagrams in the limit N c1 , N c2 → ∞ limit is exactly similar, and there is no relative factor of r = N c2 N c1 between them. Similar arguments can be applied to all planar diagrams contributing to the vertex renormalizations of the double trace quartic couplings. Consequently, all these planar vertex renormalizations are independent of r.
Coupled with same feature of the wave function renormalizations, this proves the r-independence of the planar beta functions of all the double trace quartic couplings. This completes our proof of the lemma.
Finally, we remark that one can extend this lemma to the case of a general weakly coupled QFT which consists of two originally decoupled 3+1 dimensional QFTs, T 1 and T 2 , deformed by an arbitrary number of quartic double trace scalar operators Tr where the scalar φ i belongs to the sector T i . 46 See appendix G for a derivation of the ways in which different diagrams scale in the large N limit. 47 Here, let us remark that the external blobs belonging to B E 1 in such a diagram cannot be replaced by blobs belonging to B E 2 without changing the external legs as well.

Proof of the survival of the fixed circle in the planar limit
Let us now prove that the fixed circle survives at all orders of the 't Hooft couplings in the planar limit. The fact that the planar beta functions are independent of the ratio r = N c2 N c1 ensures that the beta functions of the couplings λ i , h i and f i in the two sectors are related to each other in the planar limit by the following exchanges: This guarantees that the RG flows of the single trace couplings in one sector do not depend on the couplings in the other sector. 48 Therefore, the planar beta functions of these single trace couplings take the following forms: and B E 2 . Since each blob's loop integral is independent of the other blobs, the total vertex renormalization for a given chain of blobs would be a product of the contributions from the individual blobs. This factorization implies that a building block of the planar vertex renormalizations is an infinite sum over all possible planar diagrams within a given blob B * i . Let us denote this sum by C * i . Since we are considering a region where the single trace 't Hooft couplings λ i and h i are the same in the two sectors, the dependence of the quantity C * i on the sector to which it belongs comes only from its overall scaling with N ci , i.e., we have 50 Here we have used the fact that the external blob B E i consists of connected planar diagrams with two external legs which scale like O(1) while the internal blob B I i consists of closed connected 49 These values will receive corrections from higher loop diagrams to their magnitudes determined from the 2-loop beta functions in the previous sections. 50 Here C E and C I depend on the values of λ and h.
planar diagrams which scale like O(N 2 ci ) as N ci → ∞. We will call C E,I i 'building blobs' for obvious reasons.

Now, consider a geometric sum of linear chains of internal building blobs in a given sector
connected by the double trace coupling f i in that sector (see figure 7). We denote these sums for the two sectors by D 1 and D 2 . Based on this diagrammatic expansion, we can see that the quantity D i has the following value 51 : where α is a normalization associated with the vertex factor of the couplings f i and ζ. Its exact value is not essential to the arguments that will follow and hence we will leave it unspecified.
The only thing that we need to keep in mind is that this normalization is the same for all the double trace couplings which follows from the way they appear in the Lagrangians of the two double bifundamental models. internally and connected to an appropriate external building blob C E i at each end. As we explain below, the sum over the contributions of all such diagrams determines the part of the vertex renormalization that appears in the corresponding beta function.
Let us recall that we are working in dimensional regularization and minimal subtraction scheme. In this scheme, the vertex renormalizations of the double trace couplings take the fol- 51 Here let us make a brief comment about the symmetry factors of the different diagrams in this expansion. It will also hold for the expansions given in (5.7). The symmetry factors in these diagrams mostly arise due to the invariance under exchanges of vertices and propagators within individual blobs. Such symmetry factors are absorbed in the definition of the quantity C I and C E , and we need not worry about them while connecting the different blobs via double trace vertices. The exceptions to the above statement are diagrams where all the four legs emanating from a double trace vertex are identical. But the contributions of such diagrams are suppressed by powers of 1 N ci , and hence they can be ignored in the planar limit.
lowing form As discussed in appendix H, the residues (Z In the minimal subtraction scheme, these residues are chosen such that they cancel similar terms arising from loop integrals in the afore-mentioned diagrams. Thus, summing over the residues for all such diagrams, we get The minus sign in the above expressions is there because the vertex renormalizations have to cancel the residues of the different diagrams. The overall factor (C E ) 2 comes from the two external building blobs at the two ends. The first two terms within the residues in the expressions of Z can be understood similarly. Performing the sum over the geometric series appearing in (5.7), we get the following re-summed expressions: Now, we can use of the expression of D i obtained in (5.5) and switch to the 't Hooft couplings .
Here we have used the fact that the vertex renormalization for a 't Hooft coupling is the same as that for the corresponding original coupling. This is because the two couplings are related by a simple rescaling.
Having obtained the re-summed expressions for the residues of the planar vertex renormalizations, let us now look at the other set of quantities that contribute to the beta functions, viz., the residues corresponding to the wave function renormalizations of the scalar fields. Just like the vertex renormalizations, these wave function renormalizations can also be expanded in powers of 1 ε as shown below: Again, as explained in appendix H, the residues (Z Φ 2 ) of the above quantities at the ε = 0 pole contribute to the beta functions of the double trace couplings. In the minimal subtraction scheme, these residues are chosen such that they cancel the residues arising from planar diagrams of the form shown in figure 8. As we discussed in the previous subsection, these diagrams contain only single trace vertices. Since we are probing a region in the coupling space where the single trace couplings in the two sectors are identical, the sums over such diagrams are also the same for the two sectors. Let us denote the common value of these sums by C wfn . 52 To cancel the order O( 1 ε ) term in this quantity we have to set From the above expressions of the vertex and wave function renormalizations in the planar limit, we can now evaluate the planar beta functions of the different double trace couplings. As derived in appendix H, these beta functions are given by we get Note that these planar beta functions are independent of the ratio r = N c2 N c1 as we argued in the previous subsection. Now, let us use differentiation by parts to rewrite the above expressions of β f 1 and β f 2 as (5.14) Here we have also used the fact that the action of the operator From these expressions, we can obtain the planar beta functions of f p = f 1 +f 2 2 and f m = f 1 −f 2 2 as given below: Now, using differentiation by parts once more, we can rewrite the above expression for β fm as By comparing the expressions of β fm and β ζ , we can see that This means that the equations obtained by setting all the planar beta functions to zero (with the condition ζ = 0) are degenerate. The corresponding degenerate fixed points satisfy the following two equations: Notice that the residues in the above equations are functions of λ, h, f p and (f 2 m + ζ 2 ). Moreover, the operator A κ A ∂ ∂κ A acting on these residues takes the following form: . Thus, we have proved the survival of the fixed circles in the double bifundamental models under all loop corrections in the planar limit. As a result, we can conclude that all the points on this circle which demonstrated thermal order are genuine CFTs in this limit. Therefore, the respective baryon symmetries in these theories remain broken up to arbitrarily high temperatures and the systems exist in a persistent BEH phase.

Conclusion and discussion
In this paper we studied the possibility of spontaneous breaking of global symmetries at nonzero temperatures for large N conformal gauge theories in D = 4 space-time dimensions.
We started with some familiar QCD models, viz., QCDs with SO(N c ) and SO(N c ) × SO(N c ) gauge symmetries. The matter contents of these theories consist of Majorana fermions in the fundamental representation of the SO(N c )'s, and real scalar fields transforming in fundamental/bifundamental representation of the gauge group. By analyzing the beta functions of the different couplings in the N c → ∞ limit, we determined the Banks-Zaks-like fixed points in these models. We showed that the conformal theories at these fixed points fail to demonstrate spontaneous breaking of certain global symmetries (a flavor symmetry for the vector model and a Z 2 baryon symmetry for the bifundamental model) at nonzero temperatures.
We then extended these models by taking two copies of the bifundamental scalar QCD (with ranks N c1 and N c2 ) and coupling them via a double trace interaction. We called this the 'real double bifundamental model'. Both by direct perturbative computations up to 2-loops as well as general diagrammatic arguments at all orders in the 't Hooft couplings, we showed that this model has the following interesting property: Within the domain of validity of perturbation theory, the planar beta functions of the different couplings yield a conformal manifold with the topology of a circle. 53 This plane is determined by λ0, h0 and the value of fp obtained by solving the equations in (5.19). 54 The radius of the fixed circle is determined by the value of (f 2 m + ζ 2 ) obtained by solving the equations in (5.19).
Let us emphasize here that the proof of the existence of the conformal manifolds in the double bifundamental models at the planar limit in section 5 is largely model-independent. So it can be used to generate similar families of models with conformal manifolds. One such family would be models where there are two similar large N QFTs T 1 and T 2 , each having a single matrix-valued scalar Φ i , with only a double trace interaction between them 55 . For example, when the gauge groups in each theory are products of SU (N )'s or SO(N )'s, and the scalar fields are in either the adjoint, the symmetric, the anti-symmetric or the bifundamental representations of the gauge groups, then Tr(Φ † 1 Φ 1 ) Tr(Φ † 2 Φ 2 ) can be such a double trace interaction coupling the two theories. It would be interesting to investigate other variants as well, for example, models with multiple double trace interactions between two similar large N QFTs.
Let us now return back to the double bifundamental models. For the fixed points lying on the above-mentioned manifold, we explored the possibility of spontaneous breaking of global symmetries at nonzero temperatures. We found that such a symmetric breaking indeed occurs in certain parameter regimes. The relevant symmetries here are two Z 2 baryon symmetries, one for each sector. We found that in the zero temperature limit, i.e., in the ground state, both these symmetries remain unbroken. Moreover, in this limit, the effective potential of the scalar fields steadily increases as one moves away from the origin in the field space along any direction. This is in contrast with the models discussed in [1,2] where, in the planar limit, a flat direction of the potential in the field space allowed for nonzero vacuum expectation values of the fields which led to the spontaneous breaking of scale invariance. This also made it possible to spontaneously break global O(N ) symmetries in these models even at zero temperature. Here nothing of that sort happens due to the absence of such flat directions of the effective potential.
When a temperature is turned on, the scalar fields in this model pick up thermal masses.
If the square of any of these masses is negative, then the minimum of the thermal effective potential lies away from the origin in the field space which leads to nonzero expectation values of the scalar fields. This, in turn, means that at least one of the Z 2 baryon symmetries is broken in such a thermal state. Thus, to determine whether these baryon symmetries are spontaneously broken at nonzero temperatures for the points lying on the conformal manifold, one just needs to evaluate the thermal masses at these points. We analysed these thermal masses and made the following observations which constitute the main results of the paper: When N c2 < N c1 , the baryon symmetry in the first sector is always unbroken. On top of this, when the ratio r ≡ N c2 N c1 ≤ 6

A. Vanishing of the physical masses in the ground state
In this appendix we will prove that for all the models considered in this paper, the physical masses of the different fields are zero in the ground state. The vanishing of these physical masses is a necessary and sufficient condition for the fixed points corresponding to the beta functions of these models to be genuine CFTs.
To prove the above statement, let us first recall that we have been working in the dimensional regularization and the minimal subtraction scheme. Moreover, we have taken the renormalized masses of all the fields to be zero in this scheme. An important consequence of these is the vanishing of all diagrams with tadpoles in the perturbative expansions of vacuum correlators [78].
As we will now see, this plays a central role in the proof of the vanishing of the physical masses.
Consider any of the models introduced in section 2. Let us denote its fields by {X i }. Suppose the renormalized mass of the field X i is m X i (µ) at the energy scale µ in the above-mentioned scheme. This renormalized mass is related to the corresponding physical mass m ph,X i by the following relation: where Π X i (k 2 , µ) is the self energy of the field X i evaluated at the scale µ as shown in figure   9. This self energy is a function of the Lorentz invariant quadratic combination of the external momentum − → k due to the fact that the theory is Lorentz invariant. Then from the above equation we get Let us see if this can lead to the vanishing of the physical masses. A necessary and sufficient condition 58 for the vanishing of these physical masses is that To see that the above equations are actually true, note that when we set the external momentum to zero, all the propagators in any diagram contributing to Π X i (0, µ 0 ) are exactly as they would be in a tadpole diagram. These propagators all correspond to fields with vanishing renormalized masses. As we mentioned earlier, all such massless tadpoles vanish in our scheme. Thus, the equations in (A.3) are indeed satisfied, and consequently, the physical masses of all the fields are zero.
Let us remark here that the vanishing of the tadpoles also results in the renormalized masses remaining zero at all energy scales. This can be seen by setting the physical masses to zero in (A.1) and then applying arguments similar to those in the previous paragraph.

B. Baryon symmetry in the bifundamental scalar QCD
In this appendix, we will show that the Z 2 transformation given in ( Now, let us see how each term in the above Lagrangian transforms under the transformations given in (2.19). To be specific, we will set a = 1 in these transformations. These transformations can written more compactly as where T is an N c × N c diagonal matrix of the following form: The terms − 1 2 Tr (F 2 ) µν (F 2 ) µν and i 2 χ (p) ( / Dχ (p) ) are clearly invariant as the Z 2 transformation acts trivially on the fields appearing in these terms. The quartic interaction terms are also manifestly invariant as (Φ T Φ) remains unchanged under the transformation. So we just need to check the invariance of the terms − 1 2 Tr (F 1 ) µν (F 1 ) µν , i 2 ψ (p) ( / Dψ (p) ) and 1 2 Tr For this, we can first obtain the transformations of the field strength (F 1 ) µν and the covariant derivatives of ψ (p) and Φ as follows: From these transformations we can then derive the invariance of the above-mentioned terms as shown below: Here, we have used the fact that T † = T T = T −1 . This conclusively demonstrates the invariance of the Lagrangian under the Z 2 transformation. Now, let us show that this Z 2 transformation is an automporhism of a set of classes of gauge-equivalent field configurations. For this, it is necessary and sufficient to demonstrate that performing the Z 2 transformation over a gauge transformation is equivalent to performing another gauge transformation over the Z 2 transformation. To prove this statement, let us consider a gauge which transforms the fields as follows: Now, if we act T on the above configurations, we get

C. Two-loop beta functions of QCD
In this appendix we will review the two-loop beta functions in a general QCD without any Yukawa interaction. We will mainly follow the presentation of [87,88] where the authors employed 60 The fact that this is a gauge transformation can be verifed by checking that T O1T −1 is an orthogonal matrix and det T O1T −1 = 1.
the dimensional regularization and the modified minimal subtraction (MS) scheme. The results match with many other renormalization schemes (e.g. [89]). In particular, they are consistent with the minimal subtraction scheme that we have been working with this in this paper. We refer the reader also to [90] for the detailed forms of such two-loop beta functions.
Let us consider a general QCD Lagrangian (without Yukawa interaction) of the following form 61 : where κ = 1 2 for Majorana fermions and κ = 1 for Dirac fermions. F A µν is the field strength corresponding to a gauge field V A µ associated with an arbitrary compact semi-simple Lie group G. {φ a } and {ψ j } are sets of scalars and fermions respectively, and we denote the representation in which they transform under the gauge group by S and F respectively. We choose the scalars to be real and take the couplings λ abcd to be fully symmetric under permutation of indices. Note that the indices a and j may contain both gauge and flavor indicies. The matter and the gauge sectors are coupled through the covariant derivatives where T A (S) and T A (F ) are the generators of the gauge group in the representations S and F .
Calculation of the beta functions in such a theory largely reduces to a determination of various group theoretical quantities associated with the representations S, F and the G, and different combinations of the quartic couplings λ abcd . When G is not simple, there are several subtleties in the evaluation of the group theoretical invariants because of diagrammatic reasons. We will discuss these subtleties in the following subsections.

C.1. Gauge couplings
The two-loop beta function for the gauge coupling g in the general QCD Lagrangian C.1 is as follows: 61 Just as in the main text, we suppress the gauge-fixing and ghost terms here as well. At the two-loop level, ghost terms affect the beta functions of quartic couplings through wave function renormalization of the gauge and scalar propagators. On the other hand, the beta functions are independent of the gauge fixing term as they should be. 3) The first and second line in the the above expression correspond to the one-loop and two-loop contributions. The quantities C 2 (R) and S 2 (R) are the quadratic Casimir and the second Dynkin index respectively of the representation R (G corresponds to the adjoint representation). Note that there is no contribution from the quartic couplings at this order. Equation (C.3) is valid for any simple gauge group. For a general semi-simple gauge group G = i G i with independent gauge coupling g i for each simple G i , one has to modify this equation to obtain the beta function of g i by the following substitutions: where C i 2 (R) and S i 2 (R) denote the quadratic Casimir and the Dynkin index of the representation R corresponding to the simple Lie group G i . While employing these rules, one should consider only the matter and the gauge fields that transform nontrivially under the group G i .
The above substitution rules can be understood diagrammatically (see figure 10) when we decompose the semi-simple gauge field into the direct sum of simple gauge fields as

C.2. Quartic couplings
• One-loop The contribution of one-loop diagrams to the beta function of the quartic coupling λ abcd is where the group invariants Λ 2 abcd , Λ S abcd , A abcd are defined as follows: Here, 'perm' denotes all possible permutations of the indicies {a, b, c, d}. C 2 (k) is the quadratic Casimir of the representation in which the scalar k transforms under the gauge group. The diagrammatic interpretation of each term is straightforward: Λ 2 abcd and A abcd come from oneparticle irreducible vertices made of two quartic and two gauge interaction vertices respectively, while Λ S abcd comes from the one-loop anomalous dimension of each external scalar propagators (see the diagrams in figure 11). When the gauge group is semi-simple, we need to do the following two substitutions in the expressions of the 1-loop beta functions: where (Λ S ) i abcd and (Λ) i ab,cd are given by Based on the above substitution rules, we find it convenient to define Then the substitution rule for the last term in (C.5) is as follows: Figure 11: One-loop contributions to the running of quartic couplings with corresponding group factors e.g. in the Landau gauge.
• Two-loop More complicated structures arise at the level of two-loop diagrams where 11 different types of combinations contribute to the beta functions: Here the new group factors are defined as follows: (C.12) f g , and f ABC 's denote the structure constants of the gauge group.
For the semi-simple gauge group G = i G i , the following additional substitution rules are necessary because of the internal gluon loops: , where (A g ) i abcd refers to the invariant A g abcd defined in (C.12) for the simple factor G i . Note that the substitution rules for g 2 Λ 2g abcd , g 4 A λ abcd and g 4 A λ abcd can be obtained from the rule already given in the second line of (C.7).
Next we outline the diagrammatic origin of the different group theoretical contributions to the running of two-loop quartic couplings. For simplicity, we work in the Landau gauge.  rules are g 2 Λ ab,cd → i g 2 i (Λ) i ab,cd and g 2 C 2 (R) → i g 2 i C i 2 (R) to take account all possible gluon propagators. Figure 13 • −g 4 35 This factor comes from the g 4 contribution to the two-loop anomalous dimension of the scalar propagator where the corresponding 8 feynman diagrams are given in the figure 1 of [87]. The required substitution rules are , to cover all possible gluons inside a scalar propagator together with all possible matter contribution to the gluon propagator.
• 5 2 A λ abcd g 4 , 1 2 A λ abcd g 4 The first term comes from the diagrams 14a ∼ 14d in figure 14, and the second term originates from 14e. The required substitution rules are to reflect the fact that each gluon propagator has semi-simple indices corresponding to the simple factors {G i }.
These three terms come solely from the cubic or quartic gauge interaction vertices. The four Feynman diagrams shown in figure 15 generate these terms. The first diagram 15a generates A abcd and one needs to use the substitution rules because of the semi-simple index for the gluon lines and the corresponding one-loop insertions.
The second and the third diagrams (15b and 15c) generate A g abcd and correspond to addition of (a) Figure 14 the gluon propagator in the figure 11c. Hence, all the gluons should belong to the same node G i , which is reflected in the substitution rule g 6 A g abcd → i g 6 i (A g ) i abcd . Finally, the fourth diagram 15d generates all three group invariants. While extending A S abcd to the corresponding invariant for a semi-simple gauge group, one should use the substitution rules for A abcd and C 2 (R) introduced earlier. (C.14) To evaluate the beta functions of the different couplings, we need to identify the symmetric couplings and the generators of the gauge group in the representation of the scalar fields. These are given in (2.3) and (2.12) respectively. In addition, we need the quadratic Casimirs and the Dynkin indices of the different representations which are given below: Now, using these, we can evaluate the 2-loop beta function of the gauge coupling from (C.3), and the 1-loop beta function of the quartic couplings from (C.5). We provide the forms of these beta functions below in terms of the 't Hooft couplings (λ, h and f ) defined in (2.4): • Bifundamental scalar QCD (C.17) As explained in C.1 and C.2, for computing the beta functions of the different couplings, we need to use (C.3), (C.5) with appropriate substitution rules to take into account the different internal gluon lines corresponds to the each simple gauge group. The diagrams responsible for these rules are the figure 10 for the two-loop β λ and the figure 11b, 11c for the one-loop β f,h . The necessary ingredients to compute these beta functions are the symmetric couplings and the generators given in (2.18) and (2.26) respectively, and the quadratic Casimirs and the second Dynkin indices of the different representations. We provide these quadratic Casimirs and Dynkin indices below: where the superscript α distinguishes the two SO(N c )'s. Using these ingredients, we can compute the beta functions for the case where the gauge couplings g 1 and g 2 corresponding to the two SO(N c )'s are equal (say, g). We provide the forms of these beta functions below in terms of the 't Hooft couplings (λ, h and f ) defined in (2.21): Nc .

D. Beta functions in the real double bifundamental model
In this appendix, we will derive the beta functions of the different couplings in the real double bifundamental model up to the contributions of 2-loop diagrams. For this, we will use the formalism discussed in appendix C.

D.1. Beta functions of the gauge couplings (up to 2-loops)
In this model the gauge group has the following structure: where each G iγ is an orthogonal group with rank N ci . We will denote all the scalars and fermions in the i th sector transforming nontrivially under the group G i ≡ (G i1 × G i2 ) by S i and F i re-spectively. We will also find it useful to denote all the scalars and the fermions in the model collectively by S and F respectively .
The beta functions (up to 2-loops) of the gauge couplings are given by In the above expression, C 2 and S 2 denote the quadratic Casimir and the second Dynkin index of the corresponding representation. The super scripts (iγ) indicate the simple Lie group G iγ corresponding to which these quantities are computed. The values of these quadratic Casimirs and Dynkin indices are given by C jβ Substituting these values in the expressions of the beta functions of the gauge couplings, we get From the above expression, one can obtain the beta functions of the rescaled couplings λ i ≡ which are given below:

D.2. 1-loop beta functions of the quartic couplings
Let us now evaluate the beta functions of the quartic couplings. To use the results worked out for the 1-loop beta functions of these couplings in [88], we will introduce the following couplings which are symmetric under permutation of indices: (D. 8) In terms of these couplings, the quartic potential of the scalar fields takes the following form: where the summation over repeated indices is implicitly assumed.
The couplings h 1 and f 1 , or equivalently the rescaled couplings h 1 = N c1 h 1 16π 2 and f 1 = N 2 c1 f 1 16π 2 , appear in the expression of λ a 1 i 1 ,b 1 j 1 ,c 1 k 1 ,d 1 l 1 . Therefore, to determine the beta functions of h 1 and f 1 , we will first evaluate the beta function of λ a 1 i 1 ,b 1 j 1 ,c 1 k 1 ,d 1 l 1 . The beta functions of h 2 and f 2 can be then obtained by a (1 ↔ 2) exchange in the indices. Similarly, we will also evaluate the beta function of the rescaled coupling ζ = N c1 N c2 ζ 16π 2 by evaluating the same for λ a 1 i 1 ,b 1 j 1 ,c 2 k 2 ,d 2 l 2 .
The 1-loop beta function of the coupling λ a 1 i 1 ,b 1 j 1 ,c 1 k 1 ,d 1 l 1 , as derived in [88], is given by with the quantity (Λ) iβ a 1 i 1 ,c 1 k 1 ;eu,f v defined as follows: Let us now briefly explain the notations used in the definition of the above objects. In the first and third lines of (D.11), the sums are over all permutations of the indices (a 1 i 1 , b 1 j 1 , c 1 k 1 , d 1 l 1 ). In the second line, S 1 denotes the scalar fields in the first sector. They transform in the bifundemental representation of SO(N c1 ) × SO(N c1 ), and are invariant under the orthogonal transformations in the other sector. Therefore, the quadratic Casimir C iβ 2 (S 1 ) has the following value: The quantities T iβ A (S) in (D.12) are the generators of the representation in which the scalar fields transform under the group G iβ . For example, the generators T 1β can be chosen to take the following values: (D.14) Similarly, the generators T 2β can be chosen by (1 ↔ 2) exchange in the indices of the above expressions.
Note that the generators T 1β   Adding all the contributions we get Similarly, we can obtain the 1-loop beta functions of h 2 and f 2 which are given below: The 1-loop beta function of the coupling λ a 1 i 1 ,b 1 j 1 ,c 2 k 2 ,d 2 l 2 , as derived in [88], is given by with the quantity (Λ) iβ a 1 i 1 ,c 2 k 2 ;eu,f v defined as follows: Here S 2 denotes the scalar fields in the second sector. The rest of the notations are similar to the ones introduced earlier.
The contributions of these terms to the 1-loop beta function of ζ are given below in table 2.
(D. 24) The sum in the first term of the above expression runs over the values {a 1 i 1 , b 1 j 1 , c 1 k 1 , d 1 l 1 }. The different quantities appearing in this expression are defined below: The other structure constants f iβ also have analogous forms.
With these ingredients, we can evaluate the contribution of each of the terms in (D.    Summing up these contributions, we get (D.28) Here Λ 2 (k) = 1 6 λ k,eu,f v,gw λ k,eu,f v,gw , (D.29) The contributions of the terms in the above expressions to the 2-loop corrections in the beta function of ζ are given in table 5.  Adding these contributions we get where i is the complement of i, i.e., for i = 1, i = 2, and for i = 2, i = 1.

E. Constraints on the fixed points in the large N limit
In this appendix, we will discuss some constraints on the fixed points of the RG flow of the couplings in the real double bifundamental model. For this, we will restrict our attention to just the 2-loop planar beta functions of the gauge couplings and the 1-loop planar beta functions of the quartic couplings. In the planar limit (N c1 , N c2 → ∞), the 1-loop beta functions of the quartic couplings have the following forms: (E.1) Here λ 1 and λ 2 are fixed by demanding β λ 1 = β λ 2 = 0 which leads to the following nontrivial solutions: In what follows, we will demonstrate the following constraints on unitary fixed points 62 of the above beta functions.
• Constraint 1: When x f 1 = x f 2 , there is no unitary fixed point where the two sectors are coupled, i.e., ζ = 0.
• Constraint 2: When x f 1 = x f 2 , at any unitary fixed point with ζ = 0, we must have 16 λ, where λ is the common value of the gauge couplings in the two sectors.

(E.4)
One can define similar linear combinations for the double trace couplings to simplify the analyisis: The 1-loop beta functions of these couplings (along with ζ) are as follows: .
By setting β 1-loop ζ = 0 and demanding that ζ = 0, we get Similarly, setting β 1-loop fm = 0, we get Here, we have assumed that From the value of h m given in equation (E.4), we can see that this is equivalent to demanding When σ 1 = σ 2 , this holds true trivially because x f 1 = x f 2 . When σ 1 = −σ 2 , this is true because otherwise λ 1 and λ 2 would have opposite signs. This is not admissible because λ i is related to the gauge coupling by the relation and hence both λ 1 and λ 2 must be positive. 63 Therefore, we can trust equation (E.8). 63 Here the reality of gi is the crucial assumption which is a necessary condition for the unitarity of the theory. Now, combining (E.7) with (E.8), we get Substituting the values of h p and h m (given in (E.4)) into the above equation, we get When σ 1 = σ 2 , we get λ 2 1 = λ 2 2 =⇒ λ 1 = λ 2 since we have already shown that both λ 1 and λ 2 must be positive. However, this cannot be true for The other possibility is that σ 1 = −σ 2 . For instance, consider the case σ 1 = −σ 2 = 1. Then we get (E.14) Now the coefficient (13 + 6 √ 6) > 0, whereas (13 − 6 √ 6) ≈ −1.69694 < 0. Then the signs of λ 2 1 and λ 2 2 are opposite. However, this is not consistent with the reality of λ 1 and λ 2 . Therefore, equation (E.14) cannot be satisfied. Similarly, we can rule out the existence of any unitary fixed point with σ 1 = −σ 2 = −1. Thus, from the above analysis, we can conclude that when there is no unitary fixed point of the 1-loop beta functions with ζ = 0.

E.2. Proof of constraint 2
As before the couplings for the single trace interactions are given by The 1-loop beta functions of the couplings corresponding to the double trace interactions simplify in this case as follows: Let us focus on the last two beta functions given above. Setting them equal to zero and searching for fixed points with ζ = 0, we get Substituting the value of f p obtained from the first equation into the second one, we get The only way in which the above equation can be satisfied is if σ 1 = σ 2 = σ. In this case, we λ. λ.
(E. 22) F. Minima of the thermal effective potential In this appendix, we will investigate the minima of the thermal effective potential of the scalar fields in the the real double bifundamental model. When both the thermal masses (squared), m 2 th,1 and m 2 th,2 , are positive, one can trivially conclude that the minimum lies at the origin of the field space. On the other hand, if either m 2 th,1 or m 2 th,2 is negative, the minima would lie away from the origin, and the baryon symmetry would be broken. To analyze the location of the minima in such a situation, we would restrict our attention to the fixed points discussed in section 3.2 for which r ≡ N c2 N c1 < 1, and m 2 th,1 > 0 while m 2 th,2 < 0. First, let us employ gauge transformations to bring the matrices of the scalar fields to the following diagonal forms: The thermal effective potential (up to leading order in λ) for such a configuration is One can determine the saddle points of this potential by setting its partial derivatives with respect to all the scalar fields equal to zero as shown below: where i denotes the complement of i. The above equation has the following possible solutions: Let us consider a saddle where n i of the N ci diagonal entries of Φ i are nonzero. In F.1, we will show that for n 1 = 0, the solution corresponds to negative values of (φ 1a 1 ) 2 which is in conflict with the reality of φ 1a 1 . So such a saddle point cannot correspond to a minimum of the potential.
Later in F.2, we will argue that the minima actually correspond to n 1 = 0, n 2 = N c2 .
F.1. Saddle points with n 1 = 0 correspond to imaginary field configurations In this subsection we will prove that all the saddle points with n 1 = 0 correspond to imaginary values of φ 1a 1 . For this, let us consider such a saddle where φ a i = 0 for the first n i values of a i .
The nonzero components of the scalar fields satisfy the following equations: Solving these equations, we get the following values of the squares of the fields after substituting the original couplings by the corresponding 't Hooft couplings: To prove that the nonzero values of (φ 1a 1 ) 2 are negative, we will show that both the numerator and the denominator of the quantity within the brackets in the first line of (F.6) are positive.
First, let us consider the denominator: (F.7) Note that in section 3.
where f m ≡ f 1 −f 2 2 . The equation of the fixed circle and the fact that r < 1 further impose the condition (f m + rζ) < |f m | + |ζ| < λ. This inequality, together with the positivity of h and the fact that n 2 ≤ N c2 , ensures that the numerator is positive as shown below: (F.10) Therefore, we can conclude that (φ 1a 1 ) 2 < 0 ∀ a 1 ∈ {1, · · · , n 1 }, (F.11) which implies that φ 1a 1 is imaginary for these values of a 1 . This rules out the possibility of any saddle point with n 1 = 0 being a minimum of the potential.
F.2. The minima correspond to the saddle points with n 1 = 0, n 2 = N c2 Let us now consider the saddle points with n 1 = 0. For these saddle points, φ 1a 1 = 0 for all a 1 ∈ {1, · · · , N c1 }, whereas the nonzero values of φ 2a 2 are as follows: For these saddle points to be candidates for the minima of the potential, the above values of (φ 2a 2 ) 2 have to be positive. This can be verified by noting that for the fixed points under consideration, we have m 2 th,2 < 0, h 2 > 0 and f 2 > 0. The values of potential at these saddle points are Due to the positivity of h 2 at the fixed points, we can see from the above expression that the minima of the potential correspond to n 2 = N c2 . Substituting n 2 by N c2 in (F.12), we get Let us note here that there are 2 N c2 different field configurations that satisfy the above equation.
These configurations are as follows: where σ 1 , · · · , σ N c2 can be 1 or −1. However, many of these configurations are related to each other by SO(N c2 ) × SO(N c2 ) gauge transformations. It can be shown that all these configurations can be categorized into two equivalence classes. The configurations in each class are related to each other by gauge transformations. One can go from one class to the other by flipping the sign of just one of the diagonal entries. Thus, we can choose one representative from each of these classes as follows: h 2 + f 2 1 2 diag{±1, 1, · · · , 1}, (F. 16) and treat these as two distinct minima unrelated by gauge transformations. Since we are working in a perturbative regime where the 't Hooft couplings are small, the thermal expectation values of the field Φ 1 and Φ 2 would correspond to these minima at leading order in perturbation theory.

G. Large N scaling of planar diagrams
In this appendix we will discuss the scaling of different planar diagrams in the double bifundamental models. For this, we will take the two ranks N c1 and N c2 to be of comparable magnitudes, say O(N ), and then work in the limit N → ∞. In this limit, we will look at how different planar diagrams scale with N . In the process, we will demonstrate an important feature of these diagrams, viz., any double trace vertex in such a diagram links two otherwise disconnected subdiagrams. As discussed in section 5, this feature plays a crucial role in proving the survival of the fixed circles in these models under all loop corrections at the planar limit. Throughout this appendix, we will be using 't Hooft's double line notation to represent the diagrams. As we are interested only in how different diagrams scale with N where N is the order of magnitude of the ranks in both the sectors, we will not distinguish between the fields in the two sectors. Thus, unlike section 5, we will not use different colors to show propagators and vertices belonging to different sectors.

G.1. Scaling of bubble diagrams
Let us begin our analysis by considering bubble diagrams, i.e., diagrams which have no external legs. For convenience, we will temporarily rescale all the fields such that their propagators To be specific, let us consider a diagram with m double trace vertices. Each such double trace vertex has the structure shown in figure 16. 64 Departing from the usual convention, we put a tilde over g to avoid any confusion with the gauge couplings. comprising only of single trace vertices and propagators, scales as N 2−2 g i where g i is the genus number corresponding to that piece. Therefore, the overall scaling of the disconnected diagram obtained by removing the double trace vertices is From this we can now estimate the scaling of the original diagram with the double trace vertices.
For this, we just need to recall that in our present convention, each double trace coupling scales as O(1). Therefore, the removal of such vertices from the diagram, by itself, does not lead to any change in the scaling of the diagram. However, while removing a double trace vertex, we are also joining two pairs of propagators in the diagram. This leads to a reduction in the overall number of propagators by 2m. Since each such propagator scales as 1/N in our convention, the above reduction in their number leads to an enhancement by a factor of N 2m . Therefore, comparing with the overall scaling of the disconnected diagram as given in (G.2), we find that the original diagram scales as From the above expression, we can see that the leading contributions come from diagrams where m = m and g i = 0 for all the disconnected pieces. These are precisely the planar diagrams with double trace vertices. Each double trace vertex in such a diagram connects two disconnected planar subdiagrams. As we can see from (G.3), these planar diagrams scale as N 2 , which is identical to the scaling of the planar diagrams with only single trace vertices. Thus, our choice of scaling of the different couplings leads to a consistent large N scaling of all planar bubble diagrams.
Next, we will determine the scaling of different diagrams that contribute to the wave function and vertex renormalizations considered in section 5. Henceforth, we will revert back to our previous convention of unrescaled fields where the corresponding propagators are O(1).

G.2. Scaling of diagrams with 2 external legs
Let us first consider planar connected diagrams that contribute to the wave function renormalizations of the scalar fields. Such a diagram has 2 external legs as show in figure 17. Figure 17: Planar connected diagram with 2 external legs: Here, as well as in all subsequent diagrams, a shaded blob represents a planar connected component.
We can estimate the scaling of this diagram by first joining the two external legs and then summing over the color indices in this propagator. This leads to an enhancement by a factor of N 2 due to the introduction of two new color loops. The resulting diagram is a planar bubble as shown in figure 18. As we have already argued, such a planar bubble scales as N 2 . Therefore, taking into account the above-mentioned enhancement, we can conclude that the original diagram in figure 17 scales as N 0 .

G.3. Scaling of diagrams with 4 external legs
Now, let us look at planar connected diagrams with 4 external legs of the scalar fields. We can have two different types of such diagrams: 1. Diagrams which contribute to the renormalization of the single trace vertices, 2. Diagrams which contribute to the renormalization of the double trace vertices.
We will derive the scaling of both these types of diagrams below.

G.3.1. Diagrams corresponding to vertex renormalizations of single trace couplings
Any planar diagram which contributes to the renormalization of a single trace vertex has the structure shown in figure 19. Figure 19: Planar connected diagram contributing to renormalization of a single trace vertex.
Suppose such a diagram scales as N x . We can estimate the value of x by constructing a bubble diagram via the following procedure: Take two identical copies of the diagram in figure 19 and join the external legs with identical colors in these two copies as shown in figure 20. While joining these legs, sum over corresponding color indices.    The joining of the external legs introduces 4 new color loops leading to an enhancement by a factor of N 4 . As we have already argued, the bubble diagram in figure 22 scales as N 2 . Therefore, taking into account the above-mentioned enhancement, we can conclude that the original diagram in figure 21 scales as N −2 .

G.4. Comments
From the above analysis, we can draw the following two conclusions about connected planar diagrams with 2 and 4 external legs: 1. All such diagrams with a particular set of external legs scale identically with N .
2. Since we have shown that these diagrams can be augmented by definite procedures to construct planar bubble diagrams, the result that double trace vertices link otherwise disconnected planar subdiagrams can be extended to these diagrams as well.
In section 5, we have repeatedly used these two features of the connected planar diagrams to determine the structure of the planar beta functions of the different couplings.
H. General expressions for the planar beta functions of the double trace couplings in the double bifundamental models In this appendix, we will derive some general expressions for the planar beta functions of the double trace couplings (f i and ζ) in terms of the corresponding wave function and vertex renormalizations. For this, we will restrict our attention to the subspace where the planar beta functions of the single trace couplings (λ i and h i ) vanish. As we have shown in section 5, these beta functions (β λ i and β h i ) are independent of the double trace couplings. Hence, their roots can be determined independently. There is a discrete set of such roots out of which only one corresponds to unitary fixed points where the two sectors are coupled 65 . At this root, the single trace couplings in the two sectors are equal, i.e., we have We will now derive the expressions for the beta functions of the double trace couplings on the subspace defined by the above equation.
To derive these expressions, we will work in dimensional regularization and the minimal subtraction scheme. In this scheme, the renormalized double trace couplings are related to the corresponding bare couplings as follows: Here µ denotes the renormalization scale, Z Φ i is the wave function renormalization of the field Φ i , and Z f i and Z ζ are the vertex renormalizations of the respective couplings. The superscript B indicates the bare couplings. Since these bare couplings are independent of µ, they drop out once we differentiate the above equations with respect to µ, and we get where κ A runs over the double trace couplings in the model. Recall that we are working in the subspace where the planar beta functions β λ i and β h i vanish in the ε → 0 limit. Therefore, for 65 Here, let us remind the reader that we need to have x f 1 = x f 2 for the existence of this root.
small nonzero values of ε, these beta functions are β λ i = −ελ 0 , β h i = −εh 0 on the subspace of our interest. In the neighborhood of this subspace, it is convenient to switch to the following basis for the single trace couplings: In terms of this basis, the operator 2 j=1 β λ j ∂ ∂λ j + 2 j=1 β h j ∂ ∂h j on the above-mentioned subspace is as follows Thus, the equations in (H.3) can be rewritten as (H.6) Now, we can expand the wave function and vertex renormalizations in powers of 1 ε as follows: To satisfy the equations in (H.6) at each order in the 1 ε -expansion, we need the beta functions to have the following forms: (H.10) Finally, taking the ε → 0 limit, we get the following expressions of the planar beta functions of the double trace couplings on the subspace λ 1 = λ 2 = λ 0 , h 1 = h 2 = h 0 :

I. Finite N corrections
In this appendix we will discuss the fate of the fixed circle in the real double bifundamental model when the corrections due to finiteness of N ci are taken into account. We will assume that N c1 and N c2 are of comparable magnitudes (say, O(N ), where N is a large number). As we have already seen in appendix E , the values of x f 1 and x f 2 must be the same (say, x 0f ) at the leading order in 1 N for the existence of unitary fixed points where the two sectors are coupled to each other. While considering the finite N corrections to such a fixed point, we will allow for a difference between x f 1 and x f 2 at subleading orders in the 1 N expansion. So, we will take We see that there are two small parameters in the problem: 1 N and . If the degeneracy in the fixed points has to survive for arbitrary small values of these parameters, then each individual term in a double expansion of the beta functions in powers of these parameters must vanish for all points on the fixed circle. We will check whether this is the case for the first few terms in the expansion. Our analysis has a crucial limitation which we will discuss near the end of I.4.
With the above motivation, let us now expand all the couplings in powers of 1/N . We will denote the n th order term in the expansion of a coupling g which lies at a fixed point of the 1-loop beta functions by g (n) . The corrections to these terms due to 2-loop contributions will be denoted by δg (n) . In addition, we find it useful to define the following quantites: (I.4) We will now check whether the circle of fixed points satisfying the above constraints survives when the finite N corrections are taken into account. In particular, we will focus on possible fixed points where the two sectors are coupled, i.e., ζ = 0.
(I. 15) Note that when this condition is satisfied, the line of fixed points survives at this order.
A particular case where the above equation is satisfied is when N c1 = N c2 and x (1) (1) f 2 . As we will see in the following subsection, the 2-loop corrections ensure that this is the only case where one can get a fixed point with ζ = 0.

I.2.2. 2-loop:
Setting the O( 1 N ) in the 2-loop beta functions of the quartic couplings to zero, we get the following equations: 16(f 0m δf (1) x (1) (I.21) The above equation together with (I. 15) can only be satisfied if N c1 = N c2 = N c , and x (1) Setting the O( 1 N 2 ) in the 1-loop beta functions of the quartic couplings to zero, we get Note that this breaks the degeneracy in the fixed points. In particular, when x (2) f , a fixed point that survives is where f 0m = 0, i.e., f 01 = f 02 . Based on our experience with the O( 1 N ) terms, we will not be surprised if the fixed points for x (2) f 1 = x (2) f 2 in (I.27) disappear under 2-loop corrections.

I.4. Summary and comments
From the above analysis, we can come to the following conclusions: 1. The finite N corrections lift the degeneracy in the fixed points. Here, we would like to emphasize a limitation of this analysis. Our requirement for the validity of a fixed point was a bit too stringent as we demanded that each individual term in the double expansions of the beta functions in powers of 1 N and vanishes. There can be isolated fixed points which do not meet this criterion. For instance, when 1 N , the degeneracy in the fixed points of the planar beta functions can be lifted by subleading corrections in 1 N . However, the residual fixed points may not satisfy the above criterion of vanishing of each individual term in the double expansion of the beta functions. One would have to actually re-sum the expansion of the beta functions in the 't Hooft couplings (as we did for the planar diagrams in section 5) at each order in the 1 N expansion to correctly identify all these isolated fixed points. We have discussed this issue in section 6 as well.