Disordered O(n) Loop Model and Coupled Conformal Field Theories

A family of models for fluctuating loops in a two dimensional random background is analyzed. The models are formulated as O(n) spin models with quenched inhomogeneous interactions. Using the replica method, the models are mapped to the $M\to 0$ limits of $M$-layered O(n) models coupled each other via $\phi_{1,3}$ primary fields. The renormalization group flow is calculated in the vicinity of the decoupled critical point, by an epsilon expansion around the Ising point ($n=1$), varying $n$ as a continuous parameter. The one-loop beta function suggests the existence of a strongly coupled phase ($0<n<n_*$) near the self-avoiding walk point ($n=0$) and a line of infrared fixed points ($n_*<n<1$) near the Ising point. For the fixed points, the effective central charges are calculated. The scaling dimensions of the energy operator and the spin operator are obtained up to two-loop order. The relation to the random-bond $q$-state Potts model is briefly discussed.


Introduction
A variety of random, or disordered systems are known to have non-trivial universalities. These include the universality consisting of the well-known three symmetry classes which appear in the weak localization effects to the conductance [53] or the level correlations [54] in disordered electronic systems, and the universality in the wave function multifractality at the integer quantum Hall transition in which disorder plays an essential role [55]; even the universality in the distribution of the zeros of the Riemann zeta function is expected to have its origin in some chaotic system [45], and it may possibly be closely related to some disordered system. It is amazing that some of these systems in appropriate limits are phenomenologically well described by such simple ansatz as the one in the random matrix theory. It is, however, by no means certain that these phenomena are well understood theoretically. The most fundamental problem to be understood is how the universalities in disordered systems emerge from a microscopic structure of each system.
A lot of effort has been made to investigate this direction, and in some cases we have a partial answer; we can first formulate a system microscopically, then by using a very crude approximation (such as 'dimensional reduction' used in disordered electronic systems [54]), deduce universality directly in some carefully chosen limit. Even under such fortunate circumstances, it is usually the case that we do not have much knowledge on the way along which the system deviates from the universality away from the limit. This is, roughly speaking, because a disordered system formulated on a microscopic ground typically has highly nonlinear interactions induced by disorder and becomes almost intractable without any crude approximations. Thus, it is important to study a simple model which is tractable with a sensible approximation; analysis of such a model and ideas used there may lead to a practical way to obtain qualitative results on the deviation from the universality in more complicated disordered systems and to an insight on the emergence of the universality itself.
With this broad motivation in mind, among the diverse disordered systems, we take a problem in statistical physics, where relations between microscopic properties and macroscopic behavior of models are often quite intelligible. More specifically, we study a model with quenched disorder defined on a lattice focusing on its possible critical behavior. Understanding the behavior of such a system governed by the Hamiltonian which contains quenched inhomogeneous interactions is especially important, both from the practical perspective that there are no translationally invariant perfect crystals in the real systems and from theoretical interests inspired by the other disorder-induced phenomena.
As is well known, the solid notion of the universality class has been established in the field of critical phenomena without disorder. Nonlinear interactions inherent in a model are essential to understand the existence of the non-trivial universality classes [25]; one should confront with the nonlinearities in a more powerful framework than the mean field approximation. The ideas of scaling and renormalization group (RG) give such a framework. In this framework, one introduces a theory space and a vector field on it, which generates a flow that indicates the scale dependence of a theory considered. Assuming critical points in the second order phase transitions are described by scale invariant field theories, they correspond to the fixed points of the flow. The flow can be determined, step by step, by a path integration on the fluctuations belonging to a certain scale; the nonlinear interactions induce the mutual coupling among degrees of freedom at different scales. As a result, the global configuration of the flow can become non-trivial. Then we can read off the universality classes from the sets of RG eigenvalues which characterize the flow linearized around the fixed points. Furthermore, it is known that the ideas of RG are not only powerful to resolve the universality classes in the critical phenomena of pure systems but also flexible enough to deal with some disordered systems.
For the target of this paper, a system governed by an inhomogeneous Hamiltonian, it is natural to consider a configuration of the interactions as one realization from an ensemble that respects given probability distribution. Basic observables to be discussed can be related to the free energy in a large enough system, which is expected to be self-averaging over the disorder distribution. We should therefore study the quenched average: the average over the disorder distribution, of the free energy obtained by tracing over the statistical mechanical degrees of freedom.
Actually, evaluating the quenched average of the free energy is notoriously difficult task, since we should evaluate the average of the logarithm of the partition function in a non-translational invariant realization of the disorder. It is practical to use either one of the two methods: replica [8] or supersymmetry [54,30]. When these methods are applicable, we are left with an effective theory with the translational invariance, but this time, with the additional nonlinear interactions induced by the disorder. It should be noted that the system acquires enhanced symmetry, namely, the replica permutation symmetry or the supergroup symmetry, respectively.
For a weakly disordered system, one can consider the disorder-induced couplings in the effective theory as a perturbation to the corresponding pure theory. A natural question is whether the induced coupling destructs the critical phenomena of pure system, or not. A basic and general result on this direction is known as the "Harris criterion" that tells when the disorder can be neglected at large scales; if the couplings are irrelevant in the RG sense, the disorder can not change the universality class of the system from that of the pure system [29]. However when they are relevant or marginal, one should work harder to analyze the flow, namely, proceed to calculation of loop diagrams formed by the disorder-induced couplings. This procedure corresponds to the calculation of the RG beta function up to the second or higher order in the couping. In the case that the beta function has a zero in addition to the one corresponding to the pure fixed point, the flow can (depending on the sign of the beta function) transfer, at large scale, the theory to another fixed point; one recognizes that the disordered system belongs to a non-trivial universality class. The RG eigenvalues belonging to the universality class can be predicted perturbatively. In general, without a special reason to be integrable, one can not solve given interacting statistical model. Therefore, the importance of such a crossover from one fixed point to the other can not be overemphasized.
We shall study a one-parameter family of disordered models and discuss their universality class in two dimensions. Thus, let us briefly note the special role of two dimensions in the study of the universality class of pure systems first, and then comment on the current status of corresponding study in disordered systems. Working in two dimensions provides us with an ideal circumstance to study the crossover in depth. First, in two dimensions the conformal symmetry is infinite dimensional, and assuming the conformal symmetry (instead of the scale invariance only) enables us to classify the possible fixed points of unitary theories with minimal symmetry [3]. Second, it can be shown that the RG flow is irreversible for unitary theories ("c-theorem") [22]. Third, there are known examples of the crossover in which, at the non-trivial fixed point reachable via the flow from other fixed points, the realized enhanced-symmetries are known and the non-perturbative results can be obtained [4,5,6].
In the disordered systems in two dimensions, on the other hand, fixed points created by disorder are expected to be described by CFT's, but by more general, non-unitary ones. The lack of the unitarity gives rise to the challenging problems of determining the universality class in disordered systems. Major current examples concerned with the restrictive cases where the supersymmetry method can be used. They typically lead to logarithmic CFT's where the disorder-induced coupling is marginally irrelevant [30,31,58]. In these cases, the coupling goes to zero under the RG flow, and the models do not cause any non-trivial crossover. There is also a consideration on the running of the effective central charge defined on the basis of the supersymmetry [31], in analogy with the c-theorem in pure (unitary) systems [22]. At present, however, not much is known about the crossover cases, and hence our understanding of the disordered critical points is quite limited.
One of a few exceptional examples showing the crossover is the random-bond q-state Potts model in two dimension ‡ [8,9,20]. The model at q = 2 is the random-bond Ising model and does not show the crossover; it is equivalent to the random-mass fermion model, or using the replica method, to the multi-color Gross-Neveu model, in all of which the interaction is marginally irrelevant [30]. Now in the Z 2 -invariant scalar field theory (without disorder), a nontrivial fixed point emerges when the dimension d is considered as a continuous parameter in the range d < 4 [25]. In the random-bond q-state Potts model, it is also fundamental to consider the parameter q as a continuous number. One has then a non-trivial fixed point in the region q > 2. The perturbative calculation of the RG eigenvalues belonging to this non-trivial universality class is well under control. In particular, the theoretical prediction for the exponent of spin-spin correlation function is in good agreement with the numerical simulation [10].
In this paper, we introduce the disordered version of the O(n) model § and study the crossover in it. Our model has its own physical importance at certain integral values of n; it corresponds to the random-bond XY (n = 2), the random-bond Ising (n = 1) and the polymers (or, self-avoiding walks) in random environment (n = 0). But more importantly, we consider the models for continuous values of n; this leads to another example of the crossover in disordered system.
For continuous values of n, as explained below, we have a family of models which describe fluctuating loops in a two dimensional random background. As is well known, elementary excitations in a spin system like the O(n) model can be considered as nonlocal geometrical objects, namely, loops in the high-temperature expansion. It has been repeatedly emphasized that these loops are analogous to the closed trajectories of some particles [43,44]. Now the parameter n controls a statistics of particles in the random potential; we expect distinct, non-trivial behavior according to the value of n. In order to study these, we use the replica method as in the studies of the random-bond Potts model [8,9,20]; the method leads us to consider several, say M, layers of two dimensional O(n) models coupled each other via the disorder-induced coupling, and to take the M → 0 limit in the end.
The loops in the pure O(n) model become scale invariant at critical temperature for |n| ≤ 2. We investigate critical behavior in the disordered O(n) model using conformal ‡ Other exceptions includes the disordered Dirac fermion problems [59,60]. § The O(n) model is another natural extension of the Ising model other than the q-state Potts model. perturbation theory around the one-parameter family of CFT's corresponding to a line of the pure O(n) critical points [12]. The existence of the crossover in our model suggests that there is a one-parameter family of non-unitary CFT's. In this respect, we mention recent development of the stochastic Loewner evolution (SLE) [32], and its application to non-unitary critical points [50,51]. The line of the fixed points in our model may serves as a natural target of such study in disordered system.
The organization of the paper is the following. In Section 2, we formulate the model on a lattice, and explain the types of the quenched disorder considered. Then we use the replica method and take the disorder average. As a result, we reach an intriguing picture of particles going up and down across the two dimensional layers, thus forming a whole connected diagrams. Then we discuss the relation between the observable on a lattice and the scaling fields in a continuum limit. We see that the nontrivial cases in which the disorder is relevant occur for n < 1. In Section 3, we perform one loop calculation and discuss the existence of a non-trivial fixed point. One loop beta function suggests the existence of a threshold n * ; the fixed point exists for n * < n < 1, while the strongly coupled phase exists in 0 < n < n * . In Section 4, we perform the two loop calculation for n * < n < 1, using the full information of the four-point function provided by the O(n) CFT. The next-leading order correction for the thermal exponent and the lowest order correction for the spin exponent are then found. In Section 5, we calculate the effective central charge defined in the replica formalism. We find that this increases, along the flow, against the c-theorem which is responsible for unitary theories. We conclude in Section 6, and comment on few further directions. Appendix A provides the formulas on the critical Liouville field theory used in the paper. In Appendices B and C, we describe the calculation of the integral in two loop calculation of the beta function and in the spin scaling dimensions, respectively. Finally, Appendix D is devoted to the derivation of the integral formula and the expansion techniques. The main formula takes the form of the scattering amplitude reflecting the picture that the particle forming the loop can propagate via intermediate states while going across the replica layers.

Formulation of the Model
In this section, we formulate disordered O(n) loop models on a lattice and discuss their continuum limit. In Section 2.1, three types of the disordered lattice models are mapped to homogeneous coupled lattice models by the replica method. We introduce known CFT for the homogeneous O(n) model and discuss the relation between lattice and continuum in Section 2.2. In Section 2.3, the effective action for the continuum disordered model is derived making use of the operator product expansion (OPE).

Disordered Loop and Coupled Loop Model on a Lattice
We start with the partition function of pure O(n) model on a two-dimensional lattice: where s i is a n-component spin on a site i, and µ(s) represents a measure on the isotropic internal space; using the notation Tr s i for the tracing operation i µ(s i )d n s i ·, it satisfies Tr s i 1 = 1, Tr s i s i = 0 and Tr s i (s i · s i ) = n. The interaction is short-ranged, and the notation i, j refers to a link between the nearest-neighbor sites of the lattice.
A basic idea in this paper is to continue the parameter n to non-integral values and to take advantage of known continuum properties of the corresponding O(n) model [12]. On a lattice, the partition function (1) of the spin system with an arbitrary value of n can be interpreted as a model for fluctuating loops. It becomes especially simple on a honeycomb lattice, and then tracing over the spin degrees of freedom yields, where the summation is taken over the configuration of the closed loops [1]. Now the weight per bond of loops is t, while the weight per closed loop is n, which is called a fugacity and can take non-integral values. The model in (1) shows universal critical behavior when |n| 2, and we may use the same continuum description regardless of specific lattice structures . Another interpretation to make sense out of the non-integral values of n is possible for |n| 2. By orienting the loops, one can assign a complex Boltzmann weight e iχ (e −iχ ) to each clockwise (anti-clockwise) loop with a phase angle χ determined from the relation n = e iχ + e −iχ = 2 cos χ. This notion of the oriented loops here is standard in the coulomb gas (CG) methods [18], where the loops are mapped to the level lines of a height model [18,15,33]. One can also consider the loop as a closed trajectory of some particle. The local phase factor exp(iθχ/2π) is then associated, at a site, with each turn of the particle through an angle θ. This is a lattice version of the spin factor discussed in [44]. In a sense, the value of n controls statistics of the particles.
The qualitative behavior of the O(n) model can be summarized as follows. When |n| 2, the model has three different phases separated by a critical point t = t c ¶ . The region t < t c corresponds to the high temperature phase of the spin model, and the length of loop measured by the unit of the lattice spacing is finite. On the other hand, the average length of the loop is divergent either at t = t c or in t > t c . The O(n) model at the critical point t = t c is called in the "dilute" phase, since the fraction of the number of sites visited by some loop is zero. In t > t c , this fraction becomes non-zero and the corresponding phase is called "low temperature" or "dense" phase. It is known This is valid in the dilute phase but not in the dense phase. For instance, the intersections which may occur in the model on a square lattice become relevant in the dense phase and discriminate it from the model on the honeycomb lattice. For an explanation of the dilute and the dense phase, see below.
¶ On a hexagonal lattice, to be concrete, it is given by that the shapes of the loops in the dilute and the dense phase are described by the SLE κ with the Brownian motion amplitude κ < 4 and κ > 4, respectively [32]. We now give the partition function of the random model as where the local interactions t ij between spins are position dependent and the notation {t} refers to some definite configuration of the interaction. We consider the configuration {t} as a realization taken from some ensemble with a probability distribution functional P [{t}]. One might assume a non-local probability distribution functional, but here we restrict ourselves to study the case of short range correlation. This means that the interaction t = t ij on each link independently respects single distribution function P (t). Later, we shall let P (t) a Gaussian-like distribution.
Given an inhomogeneous realization {t}, the tracing over the spin degrees of freedom is still possible: where a l (i) (i ∈ {1, 2, · · · , L l }) denotes a site which belongs to the l-th loop of a length L l . Now the path of the particles forming loops should avoid the links with higher cost; we expect its behavior to change, according to the value of n. The situation is, to some extent, analogous to the problem of the electrons in a random potential which has been studied in connection with the Anderson localization [52]. As the weight t is different from link to link ( Figure 1-(a)), the exact methods applicable in the pure O(n) model, such as the CG method, can not be applied here.
Assuming the short-range correlation between the disorder, we can study the selfaveraging quantities such as the free energy and the translationally averaged correlation functions by calculating the quenched average [16]. For example, neglecting the surface effects, the free energy of a large enough system A can be considered as a sum of the free energies of many macroscopically large subsystems of A each with a different realization. Since by the short-range assumption these realizations are independent each other, the total free energy per site of A in the thermodynamical limit takes, with the probability one, the quenched averaged value defined by Here, we have used the fact that a free energy per site of a subsystem with N sites is given by f = (− ln Z)/N. The overline is used to indicate the averaging over the distribution. It should be noted that, in this argument, we need a macroscopic number of macroscopically large subsystems.
In this paper, we use the replica method to evaluate the various quantities related to the quenched average of the free energy (5). This method is based on the identity: Using this, the problem of calculating the quenched average ln Z is reduced to the task of obtaining another quenched average Z M correctly for M ≈ 0. To evaluate the latter, we first prepare the M ∈ N layers of replicated inhomogeneous O(n) models with the same realization {t}, and then take an average over the distribution. Since, in general, the moments of the interaction t ij are non-zero, the resulting theory is a coupled M layers of O(n) models. We calculate the quantity Z M for finite M ∈ N, take the limit M → 0 in the end. This procedure gives, at least formally, gives the desired average. Although we perform the detailed analysis in continuum theory, let us proceed, for a while, within the discrete lattice formulation. The purpose is two-fold: (i) to get some insight on underlying physics, as we can concretely see the elementary processes existing in the model and (ii) to see that there are possibly many lattice models described by the same continuum field theory. Now, we consider the M copies of the model (3) and take the disorder average. This gives where in the last line, the first two moments are denoted as t 1 = t ij and t 2 = t 2 ij , and we omit the terms with the higher moments for simplicity. The term t 1 s represents a walk on the same replica layer; the generic case with t 1 = 0 is discussed in the following. In addition to this, the theory acquires nonlinear couplings. For instance, when we look at a particular path of a loop, the coupling t 2 s for a < b is considered as a branching process across the replica direction as shown in Figure 1-(b). On a lattice, evaluating directly the contribution of the diagrams with the second or higher orders in such processes should involve some sophisticated enumerative combinatorics. Notice, however, that the translational invariance of the individual O(n) model is now restored. Thus, in the continuum limit, we can use the knowledge of the homogeneous O(n) CFT.
The restriction on the sum to off-diagonal sector (a < b) in the t 2 term follows from the fact that the original partition function (3) is linear in the bonds (t ij s i · s j ). When we talk about the universality, however, this formal linearity and resulting off-diagonal property should be reconsidered; we should rather consider a larger parameter space including more general couplings generated under block spin transformations.
In order to see this, firstly, it is instructive to remind that the pure model (1) corresponds to the Hamiltonian H = − ln(1 + ts i · s j ), and originates from the high temperature expansion of another spin system which has the same O(n) symmetry but has the different Hamiltonian H = −β s i · s j . From the RG perspective, these two models are regarded as two different points in the same theory space determined from the O(n) symmetry. Actually, both systems are believed to belong to the same universality class. In fact, with the knowledge of the scaling dimensions of fields obtained by the CG method [18], we may argue that these two Hamiltonians are equivalent modulo irrelevant operators. In particular, the field (s i · s j ) 2 represents the double occupancy of the bonds in the partition function, and the presence of it enables the loops to overlap. Although such a field is not contained in the partition function (1), it is natural to think that the term (s i · s j ) 2 is generated, under block spin transformation, from the single occupation (s i · s j ). But, as we will see in the end of Section 2.2, the term (s i · s j ) 2 is irrelevant in the dilute phase, or equivalently, near the critical point. Because of this, we can say that the linear expression (1) nicely summarize the characteristics of the dilute phase in generic models with the O(n) symmetry.
Taking these considerations into account, we now briefly discuss the theory space of the coupled models derived from the disordered ones. Our concern is near the decoupled critical point; applying the knowledge of the pure O(n) model, we note the inessential differences that arise from different formulations of disordered models on a lattice. For concreteness, let us consider another disordered model: with a bond distribution P (β). The formal operations lead to, again, a coupled model: where β 1 = β ij and β 2 = 1 2 β 2 ij − β ij 2 are the first and the second cumulants obtained from the distribution P (β ij ), respectively. Comparing this model (9) with the previous one (7), we notice that the differences between these two are simply the presence of the exponential function and that of the diagonal coupling: i.e. the process of the double occupancy of the bonds on the same replica plane. Such a term can also be generated in the model (7) under block spin transformation. However, when the corresponding pure model is in the dilute phase, the process is irrelevant. For this reason, we expect both of the coupled models (7) and (9) are described by a same field theory. A related model, which is more suitable for taking continuum limit, can be formulated with a probability distribution P s (µ) for disorders µ on each site: In this model, the bond strengthβ ij = µ i + µ j respects a distributionP obtained by the convolution of two P s i.e. P s * P s =P . If this distributionP (β) is identical to the distribution P (β) in (8), then this model becomes similar to that of (8). Still, the model (10) is different from the model (8), since there exists the correlation between the nearest-neighbor bonds (sayβ ij andβ ik ) connected at a site. Nevertheless, the correlation is still short range; we expect a similar behavior at large scales. We can argue that this correlation existing in the model (10) is not essential to distinguish it from the model (8) by examining how these two models are mapped into the theory space of the coupled models. Averaging over the disorder µ in (10) as a leading nonlinear coupling. In particular, the j = k part of this coupling is induced by the correlation between the nearest-neighbor bonds in the model (10). If this coupling were, with some special reason, never generated in the model (9) mapped from (8), we would seriously distinguish (10) from (8) in the first place. But, it has already been contained in (9), implicitly. Indeed, although the expression (9) is written with respect to links, this can be recasted in terms of sites thanks to the discrete rotational and the discrete translational invariance of the coupled model on a lattice. Then we recognize, by expanding the exponential, this coupling is also present in the model (9).
To summarize, in the dilute phase, there are two important processes on a lattice: the walk on the same layer s It is plausible to think that these coupled models discussed here are, at large scales, described by a simple continuum field theory.

Critical point of the homogeneous O(n) model
When we approach a critical point of some lattice model defined in the dimension d, the correlation length ξ of the model, which is usually the order of the lattice cut-off a, becomes arbitrary large. The continuum limit (or, scaling limit) of the lattice model can be reached by taking a → 0, while keeping ξ fixed. There is nice class of local lattice observables {φ lat k } which have the finite limit lim for an appropriate choice of scaling dimensions {2∆ φ k }. With the idea of the block spin transformation, the scaling fields φ k (x k ) may be considered as coarse grained versions Lattice Continuum primary fields scaling dimension energy Table 1. Correspondence between lattice observables and continuum fields. The value of p is given by of φ lat k over some region of size L, centered at x k , with a ≪ L ≪ ξ. These scaling dimensions are determined dynamically from the interactions in the system, and related to the RG eigenvalues y φ k via y φ k = d − 2∆ φ k . When y φ k is positive (negative), the RG flow near the critical point is unstable (stable), and the corresponding observable and fields are relevant (irrelevant).
Around the critical point in the theory space of the O(n) model, there are, at least, two directions which are unstable (i.e. relevant) under the RG flow to the infrared. These directions are spanned by the coupling constants associated with the two most important pairs {φ lat k , φ k } of local lattice observables and the corresponding scaling fields. On the lattice side, one is the energy density j s i · s j and the other is the spin vector s α i . In the scaling limit, they renormalize to become the scaling fields called the energy operator E(x) and the spin operator σ(x). The scaling dimension of E(r) and that of σ(r) is related to the leading thermal and the leading magnetic eigenvalue of the spin system, respectively. These two positive eigenvalues correspond to the linearized RG flow near the critical point, and determine the thermodynamic exponents [16].
So far, we have described the qualitative structure of the theory space focusing on the correspondence between the lattice observables and the continuum scaling fields. It is tantalizing that there is no versatile scheme for a given lattice model to establish this type of correspondence and to proceed to quantitative analysis using a scaling theory at hand. However, the O(n) model was particularly successful in this respect; the geometric nature such as the loop representation of the partition function made it possible to access some of the exact results on the scaling dimensions by the CG method [18] before the emergence of the conformal field theory [3]. These exact results by the CG method, or by the other means [14] then led to the conjecture claiming that a certain one-parameter family of CFT's describe the critical points of O(n) model for the continuous value of n (|n| ≤ 2) [12]. This claim is further checked by the Bethe ansatz result [2].
We assume this one-parameter family of CFT's [12] in the rest of the paper. The correlation functions of primary fields are represented in terms of the vertex operators {e iα k φ(x) } with an appropriate choice of charges {α k }, where φ(x) is a bosonic scalar field [12,15,33]. The necessary formulae on the vertex operator representation are given in Appendix A. In this construction, we have two marginal operators (the screening vertex operators). In the O(n) model, the screening charge α − is related to n as The energy operator and the spin operator is identified with the primary fields as follows: where the value of p is given by p = 1/(2 − 2α 2 − ). These identifications are based on the dimensions originally obtained by the CG methods given in Table 1. Using the formulas (A.5) and (A.6), the dimensions of these operators are checked as where The relation to the two dimensional critical statistical models are summarized below. The O(n) model at n = 2 (α 2 − = 1) and n = 1 (α 2 − = 3 4 ) are the XY model and the Ising model, respectively. The Ising model is the only minimal unitary model [3] that belongs to the O(n) family. Note that the Ising model is also regarded as the q-state Potts model at q = 2 [8]. The model at n = 0 (α 2 − = 2 3 ) corresponds to the polymer, or self avoiding walk (SAW), as the qualitative result was obtained by ǫ = 4−d expansion [26]. The model in two dimensions is related, under the SLE duality [32], to the percolation (Potts model at q = 1). Both models are non-unitary, and various scaling dimensions in the models are numerically studied in ref. [28]. Finally, the model at n = −2 (α 2 − = 1 2 ) is the loop-erased random walk (LERW). The one-parameter family of O(n) CFT's are related to the SLE by where κ is the strength of the Brownian motion which drives the evolution SLE. These critical O(n) models covers 2 ≤ κ ≤ 4.
In the rest of the sub-section, we discuss more on the correspondence between the lattice and the continuum. The relation between the loop configuration and the correlation function in the continuum theory is given; the latter is the object of our study henceforth. The reader who is not interested in the lattice may skip the following.
In this paper, we restrict our consideration to the calculation of the RG flow in the subspace spanned by the coupling constant associated with the energy operator E and the spin operator σ. There is, however, another important RG eigenvalue apart from the ones associated with these two; the eigenvalue is responsible for the geometric property of the loops. In other words, the thermal and the magnetic eigenvalues are not sufficient to characterize the local shape of the loops even at a primitive level. To see this, let us first remind that, in the SAW (n = 0), the fractal dimension D F of the loop is given by the thermal eigenvalue y E , and thus D F = 1/ν [26]; it is then given by The generalization of this result to arbitrary n (|n| ≤ 2) was discussed by relating the O(n) loops to the clusters in the percolation [57]. The fractal dimension for n = 0 should then be written as with σ F = 1 (in the percolation theory, σ F is called the Fisher exponent, and ν the correlation length exponent). Further, the fractal dimension D F is expected to be the RG eigenvalue from the polarization operator (or, the two-leg operator): } with the vector indexes α = β, for which the scaling dimension can be obtained by the CG method (see Table 1). The fractal dimension of the O(n) loop in the dilute phase is then given by D F = 2 − 2∆ σ II = 1 + α 2 − /2 + . Diagrammatically, the polarization operator σ II (x), when inserted into a correlation function, creates two curves emanating from a point x, while the spin operator σ(x) creates only one curve. The curves created by σ II (x) repel each other; they have different colors α = β and hence can not connect by themselves. They need to be annihilated by the other operator σ II (y). By contrast, the insertion of the energy operator E(x) require that a point x be passed by some loop; the curve passing the point will connect each end to form the loop.
However, forming a loop becomes more and more difficult as n tends to zero; the loop segments repel each other in the SAW. In this sense, the role of the energy operator E approaches that of the polarization operator σ II . In fact, at the SAW, these two operators has the same scaling dimension 2∆ = 2 3 , and this explain σ F = 1 in (17). In the CFT context, by using the scaling dimensions obtained by the CG method, σ II (x) was argued to correspond to the φ 1,0 primary field. Then, σ F = 1 can be confirmed by the well-known equivalence between the primary field: in the m-th minimal model with the central charge c = 1 − 6/m(m + 1). Indeed, in the SAW, which is the m = 2 minimal model, φ 1,3 (the energy operator E) should have the same dimension as the field φ 1,0 (the polarization operator σ II ), since There are, of course, many other composite fields constructed as products of the several local lattice observables. A quantitative discussion on such fields is given, which is based on the OPE and the correlation inequalities [19]. Their result suggests, in general situation, that the scaling dimension of a composite field exceeds the sum of those of the elementary fields due to the repulsion between them. This implies higherorder composite fields are more likely to be irrelevant in the RG. In our context, we mention the four-leg operator { j (s i · s j ) 2 , σ VI (x)} as an important example [18,33]. The insertion of the field σ VI (x) into a correlation function corresponds to requiring two of the loop segments to overlap; in the particle picture, the trajectories should encounter at the point x. Again, the dimension of the field σ VI (x) is obtained by the CG method, and identified with the primary field φ 2,0 in the CFT. Now, we see that it is relevant in the dense phase (κ > 4) and irrelevant in the dilute phase (κ < 4) as shown in Fig 2. As used in Section 2.1, the latter observation is a key ingredient in understanding the universality of the pure O(n) model in the dilute phase. Further, at the level of coupled model, this is crucial in our expectation that (7) and (9) are, near the decoupled critical points, described by a same field theory.

Scaling limit of the disordered O(n) models
We shall go into the continuum formulation of the disordered models by using the CFT description of the critical O(n) model. In section 2.1, we have considered, on the lattice, three disordered models whose partition functions are given by (3), (8) and (10). As we have seen, after the disorder averaging, they have much in common when mapped to the coupled models. When the couplings between replicas are sufficiently small, one could, by looking at (7) for instance, guess the action of the coupled models in the continuum limit. But we shall take the other way. Instead, we first consider the continuum limit of the disordered model (10), and then take the disorder average in the continuum theory. The advantage is that we can use the OPE to deal with the composite operators.
In the following, we will use the relation: The first equivalence symbol means that these two pure model belong to the same universality class, as discussed in the paragraph above (8). In the right hand side, we formally write the action of the O(n) CFT as S * . The coupling constant m is proportional to the reduced temperature, which is given by ( Using the relation (19), it is natural to assume the continuum limit of the model (10) is described by where S (a) * is the action of the O(n) CFT on the replica plane (a). Here, we formally write elementary fields (in the sense of [7]) as Φ (a) (x) and replace the tracing Tr s (a) i in (10) to DΦ a (x), which denotes the path integration over the fluctuation of Φ (a) (x). The realization {µ} of the local weights µ i on lattice sites in (10) are now replaced to its continuum counterpart, that is, a configuration of a scalar function m(x). Physically, this m(x) serves as a locally fluctuating reduced temperature. We assume this m(x) respects a single distribution functionP (m(x)) which is independent of the position x. Hereafter, we suppress the argument of the partition function Z [{µ}, n]. Then the averaging operation is written as where Dm(x) denotes a path integral measure x dm(x). Letξ k be a k-th cumulant of the distribution functionP (m(x)). In the important case of a Gaussian distribution Then, by the averaging over the disorder, we obtain In this formal expression, we should note that the nonlinear partŜ I (x) contains the composite operators i.e. the products of the several fields at the same point on the same replica plane. We deal with them by using the fusion rules in the O(n) model. According to the identification E → φ 1,3 in (13), the fusion rule in the thermal sector reads or, equivalently, Here, E ′ = φ 15 is the next-leading energy operator. Since this operator is irrelevant as shown in Figure 2, we neglect the E · E → E ′ part of the OPE (25). The identity field I = φ 11 contributes as a trivial shift of the free energy. The emergence of the energy operator E in the right hand side of (25) is a characteristic of O(n) model with n = 1, which will be discussed in the paragraph below (35). Because of this contribution (E · E → E), qualitatively speaking, we should have a hierarchical flow of the coupling constants: · · · →ξ 3 →ξ 2 →ξ 1 . For example, the diagonal (a = b) part of and hence the flowξ 2 →ξ 1 occurs * . Thus, by redefining the coupling constantsξ k , we get where we use the symbol m 0 , g 0 and ξ k (k ≥ 3) for the new coupling constants. Now, the effective interaction S I (x) is defined without composite operators. It should be noted that we restrict our considerations on the theory space which is replica symmetric; we assume that the coupling constant for, say These coupling constants are determined by the distribution functionP (m(x)) in (21), and are non-zero in general. The massless limit (m 0 → 0) of the decoupled model (g 0 = 0, ξ k = 0) remains obviously as a fixed point. We now change the distribution functionP (m(x)), and gradually turn on the effective interaction S I in (27) while keeping m 0 = 0. The large scale behavior is then dominated by the most relevant field: E (a) (x)E (b) (x). The dimension of this operator is given by 4∆ E , which is the twice the dimension of the energy operator.
As we can infer from (14), this field becomes marginal at n = 1 (α 2 − = 3 4 ). Since O(1) model is the Ising model, this serves as a simple check of the known marginality of the disorder in the random-bond Ising model [8,30]. Hence, we use the parametrization: to perform the epsilon expansion in the next section.

Renormalization group flow in the one-loop calculation
In this section, we discuss the scale dependence of the theory (27) by one-loop RG calculation. It will turn out that under certain circumstances, the disordered O(n) model has a pair of the ultraviolet (UV) and the infrared (IR) fixed points, as in the random-bond q-state Potts model [8,20]. We will investigate the properties of the IR fixed point by which the large scale behavior of the theory is dominated.

The epsilon expansion and procedure of RG
We shall calculate the RG flow of the coupled model (27) perturbatively in the epsilon expansion using the parameter in (28). Generally, a flow is determined only by the * In the coupled model at M = 0, this flow causes shift in the critical temperature.The magnitude of this effect is proportional to the number of diagonal elements M , the structure constant (35) andξ 2 . However, since the disordered model corresponds to M → 0, we assume it is negligible.
most fundamental properties of a theory such as the symmetries and the dimension of the space. In order to grasp the topology of a theory space, it is often useful to think that these symmetries are dependent on some continuous parameters [7,8,20,24], or that the dimension itself is a continuous parameter [25]. Although the change is gradual when we look locally at a generic point as varying these parameters, the global structure of the flow, on the other hand, may change drastically at certain critical values of the parameters. A typical topological transition is caused by a branching of one merged fixed point into a pair of the UV fixed point and the IR fixed point. The idea of epsilon expansion provides us with a firm ground to perform a perturbation calculation in the non-trivial region after the branching.
In the parameter space, we can infer the location of the branching by the presence of a marginal field: the hallmark of the merged fixed points. In the renowned example of the Z 2 -invariant scalar field theory [25], the most relevant interaction φ 4 becomes marginal at the dimension d = 4. In d = 4 − ǫ, the Gaussian fixed point and the Wilson-Fisher fixed point emerge; we can do a solid perturbative calculation by setting an appropriate coupling coordinate in which the coupling constant remains O(ǫ).
In the disordered O(n) model, as mentioned in the previous section, the most relevant part of the effective interaction S I in (27) is the term E (a) (x)E (b) (x), which is marginal at n = 1 and is relevant for n < 1 (see Fig. 2). The next-leading relevant term is irrelevant for n > 0 and is marginal at n = 0 (SAW) ‡. We restrict our analysis on n > 0, keep the most relevant part and discard the other higher-order terms in (27). The implementation of the RG is as follows. We introduce a cut-off length scale r, which serves as an effective lattice spacing, and determine a renormalized coupling g(r) by integrating out the short-distance degrees of freedom: where the insertions of the interaction fields S I are restricted onto a disk of radius r, and the coupling is measured by projecting perturbative contributions on S I (∞). The symbol ... 0 represents the unperturbed correlation function. A trivial scaling factor r −8ǫ is introduced in order to make the renormalized coupling dimensionless. At this stage in RG, as is well known, a finite redefinition of the coupling is possible; physical quantities (scaling dimensions, for instance) is invariant under a finite coordinate ‡ It may be noted that the role of this field is analogous to that of the φ 6 term in the Z 2 -invariant scalar field theory. This term is irrelevant if d > 3 and becomes marginal at d = 3. In 3 < d < 4, we have only the two fixed points mentioned above. reparametrization of theory space [24]. To fix this ambiguity, we choose a minimal subtraction scheme. Then we get the renormalized coupling in the form: where G 2 (ǫ, r) and G 3 (ǫ, r) have only poles in ǫ and no regular part §. The beta function is then obtained by differentiating the renormalized coupling (31) with respect to (ln r).
Using the first two terms in the fusion rule (25), we make contractions and list, in Figure  3, the possible diagrams for the beta function.

One-loop beta function, a line of random fixed points and a strong coupling region
We here calculate one-loop beta function. The process represented by the diagram in Figure 3-(a) involves three layers, and the number of the ways for this contraction is 4(M − 2). Then we have where we have used E(x)E(y) = |x − y| −4∆ E and 4∆ E = 2 − 8ǫ. Besides this, up to one-loop order, there is the other contribution represented by Figure 3-(b) which is obtained by using twice the sub-leading part of the OPE (E·E → E): 2 § More precisely, we keep the combination r 8ǫ /8ǫ = 1/8ǫ + log r + O(ǫ) and discard the other parts.
Here, C E EE (α 2 − ) is the structure constant of the CFT which appear in the OPE: where we write just the primary operators to represent their conformal families including the descendants. The structure constants are, in general, determined by the requirement that the operator algebra be associative. In practice, the crossing symmetry of the fourpoint function is strong enough to fix them [3]; the actual values are obtained by using either the connection matrix of the hypergeometric function [11], or the vertex operator representation of the four-point function [13]. This reads, where we have used temporary notation ρ = α 2 − and γ(x) = Γ(x)/Γ(1 − x). The structure constant C E EE (α 2 − ) has information about the important selection rules in the pure O(n) model. First, observe that the square of the structure constant (35) has a zero at α 2 − = 3 4 (the Ising model, ǫ = 0) and a pole at α 2 − = 2 3 (SAW, ǫ = 1 12 ). The former zero is due to the self-duality of the Ising model [17]. In the vicinity of the critical point, the duality transformation changes the sign of the reduced temperature t. Since the energy operator E couples to t, if the model is invariant under the duality, it should be odd under the duality. Then E is not allowed to appear in the right hand side of the fusion rule (25).
The pole at α 2 − = 2 3 emerges because of the strong repulsion between the loop segments in the limit n → 0. Since the normalization of operators are fixed such that C I EE = 1, what really happens is the divergence of the ratio C E EE /C I EE . When the two loop segments approach, they will either (i) form a complete loop (to contribute the free energy) or (ii) form together a joined loop segment. Since, in the limit n → 0 the process (i) is strongly suppressed by the repulsion between the loop segments, the ratio, indeed, diverges.
Taking these interpretations into the account, more intuitive representations of the second-order contributions are possible. The term G 2,1 (r, ǫ) is represented by the diagram in which two open segments lie on two layers and one closed loop on another layer ( Figure 4-(a)), while the term G 2,2 (r, ǫ) is represented by the diagram in which two open segments lie on two layers (Figure 4-(b)). Now we sum up (32) and (33) to get By introducingg = r 8ǫ g 0 , we solve (36) with respect tog: µ´ µ Figure 4. The second order diagrams; see also Figure 3. (a) G 2,1 (r, ǫ): the three-layer process with a closed loop (b) G 2,2 (r, ǫ): the two-layer process.
Using this, the beta function up to one-loop order is then obtained as This beta function of the coupled model permits an IR fixed point when the coefficient of g(r) 2 is negative. As we take the limit M → 0 for the disordered model, there are a competition between the terms from G 2,1 (r, ǫ) and G 2,2 (r, ǫ). As a result, we get two distinct regions for 0 < n < 1: in the first case we have a monotonic increasing beta function; the theory goes to the strong coupling. In the second case we have an IR fixed point. The one-loop beta function suggests that the IR fixed point vanishes at ǫ = 1 12 · 0.770861 · · ·, which corresponds to n = n * = 0.26168 · · ·. The flow of the coupling constant is schematically shown in Fig. 5. It should be noted, however, that the divergence of the coupling g at n = n * is an artifact of the one-loop calculation; if there is a positive g 3 term ¶, the line of the IR fixed point terminates at (n c , g c ) = (n ′ * , g ′ c * ) with a finite coupling g ′ c * . The strong coupling phase of the disordered model (M → 0) near n ≈ 0 is formed by the dominance of the two-layer process G 2,2 (r, ǫ) over the three-layer process with a closed loop G 2,1 (r, ǫ). The qualitative behavior of the absolute ratio R(n) = |G 2,2 /G 2,1 | is independent of M except a special case with M = 2, where the three-layer process is prohibited. Diagrammatically, a large R(n) implies that the ratio of (the number of inter-layer hopping) to (the average number of complete loops per layer) is also large. Near n ≈ 0 region, the particles, which are destined to form a whole connected diagram, favor to sew layers together rather than to stay on one layer and to walk around avoiding their traces.
We shall henceforth proceed to the two loop calculation to investigate the IR fixed points in the region near n = 1 (ǫ ≪ 1 12 ), where regarding the structure constant C E EE in (35) as O(ǫ) is justifiable. In this region, the calculation turns out to be parallel to that in the study of the random-bond Potts model by Dotsenko, Picco and Pujol [20]. At one-loop level, C E EE = O(ǫ) implies G 2,2 = O(ǫ); in the minimal subtraction scheme, A similar form of beta function appears in ǫ = d − 2 expansion of the random-bond Ising model [16]. ¶ Although the ǫ dependence of the third order coefficient is not calculated in this paper, we anticipate a positive g 3 term from a continuation of the result at ǫ ≪ 1 12 in Section 4. we drop this term. The beta function is then We left the two loop calculation of the intermediate region where C E E,E ≈ 1 as a future problem.

The beta function at two-loop
There are four types of third order diagrams for the beta function as listed in Figure  3-(c), (d), (e) and (f). As the diagram G 3,3 (r, ǫ) in Figure 3-(e) have extra ǫ 2 factor from the structure constants and the diagram G 3,4 (r, ǫ) in Figure 3-(f) has no pole in ǫ, we omit these terms in the following. The G 3,1 (r, ǫ) in Figure 3-(c) can be calculated as where we have used the asymptotic behavior of the integral for 1 ≪ R: It should be emphasized that we get no simple pole in (40). Next, we have G 3,2 (r, ǫ) in Figure 3-(d) containing four-point functions: where we define Actually, instead of I(R, ǫ), we shall calculate I(∞, ǫ) by analytic continuation in Appendix B, and see that I(∞, ǫ) itself has no poles in ǫ. Although, in the Appendix, we prove this fact by applying the contiguity relation between generalized hypergeometric series, the fact can also be interpreted as a result of a cancellation between two poles with opposite sign, as explained in the following + .
In the limit |z| → ∞, the integrand in (44) behaves as |z| −4∆ E (4∆ E = 2 − 8ǫ), since the four-point function, using the identity operator I in (34) as the leading intermediate channel, decouples into two two-point functions. Hence, I(∞, ǫ) has the contributions from both regions |z| < R and |z| > R, each of which contains a pole +2πR 8ǫ /8ǫ and −2πR 8ǫ /8ǫ, respectively; they are canceled out each other. Observing that I(R, ǫ) does not have the latter pole, and estimating the corrections, we obtain As the first term corresponds to the decoupling limit y → x of the four-point function in (42), it leads to the same contribution as (40) where N is the normalization factor with the four-point function given in (A.7). Substituting (45) and the result (B.30) for I(∞, ǫ) into (43), we obtain We get, by collecting the third order terms (40) and (47), the renormalized coupling constant in (31) as Then (48) is solved with respect tog = r 8ǫ g 0 as + A similar cancellation of poles has already been observed in the random-bond Potts model [20].
Consequently, we obtain the RG beta function As is expected from the known equivalence of the random-bond Ising model to the random-mass fermion model and the Gross-Neveu model [30,8], at the Ising point ǫ = 0 (n = 1), the expression (50) reduces to the two-loop beta function for the Gross-Neveu model [34]. Solving β(g c ) = 0 for the disordered model (M → 0), we obtain the coupling constant at the random fixed point as

The correction to the scaling dimensions
In our perturbation theory, the two-point correlation function between fields O(0) and O(∞) is calculated as Again, as in (30), we restrict the insertion of the interaction S I around the field O(0) within the disk of radius r. From this, we can determine a renormalization constant Z O for the field O perturbatively as The two-point functions in the theories with different cut-offs r and sr are related as where, in the second line, we have introduced the anomalous dimension Now consider the large s behavior of (54) and let the upper limit g(rs) of the integral tend to g c of the random fixed point; the beta function β(g) tends to zero, while the anomalous dimensions γ O (g), as we will see, remain finite both for the energy operator E and the spin operator σ. This means the integral is dominated by the contribution from the region g ≈ g c , and we obtain where we mean, by 2∆ IR O and 2∆ UV O , the scaling dimension of the field O at the IR (the random) fixed point and the UV fixed point, respectively.

Scaling dimension of the energy operator
The first and second order diagrams for the renormalization constant of the energy operator Z E are listed in Figure 6-(a), (b), (c) and (d). As the integral in the diagram in Figure 6-(c) does not have pole in ǫ, we omit this term from the calculation. Since the integrals are same as those for the two loop beta function, by adapting combinatorial factors, we have where the results (32), (40) and (47) are used in the second line. Again, by using the expression forg in (49), we get Consequently, by using (51) in (56), we obtain

Scaling dimension of the spin operator
We shall here calculate the renormalization constant for spin operator Z σ up to the third order in g 0 , which gives the lowest order correction to the scaling dimension 2∆ σ at the IR fixed point. According to the operator algebra (see the diagrams in Figure 6), there are no contribution at the first order, and one contribution at the second order S 2 (r, ǫ) and the other one at the third order S 3 (r, ǫ): The second order diagram in Figure 6-(e) is given by while the third order diagram in Figure 6-(f) is given by · |y−x|<r, |z−x|<r, |w−x|<r In evaluating the integrals (61) and (62), we keep only the terms with the lowest powers in ǫ. For this reason, it turns out that we can add certain extra regions to these integrals. First, by adding the region |z − x| > r to (61), we get with K 2 (∞, ǫ) defined as Second, if we add the region |w − x| > r and |z − x| > r to (62), we see that the third order contribution S 3 (r, ǫ) has the similar structure as the second order one S 2 (r, ǫ) in (61). In fact, under a trivial change of variables, we obtain the factorized form: with the definition Then, we write the integrals (64) and (66) using the vertex operator representation: and give the results in Appendix C. Now we obtain, by collecting the contributions (63) and (65), the renormalization constant for the spin operator as where we have used (41) to obtain the third order term. Differentiating this with respect to (ln r) and using the result (C.10) for the integrals K 2 (∞, ǫ) and K 3 (∞, ǫ), we obtain the anomalous dimension for the spin operator as where the expression for the bare coupling constant (49) Now, from (56), we obtain where, in the last line, we have used the elliptic integral of the first kind defined by for the purpose of the comparison with the result in the random-bond Potts model [20]; this will be discussed in Section 6.

Effective central charge
The central charge, in general, characterize a system by counting the number of the critical fluctuations in it. If the replica method is used, however, the central charge becomes rather trivial; we are always left with the vanishing central charge c = 0 because of the formal limit M → 0 in the end. Nevertheless, we can define a so-called effective central charge characterizing a random fixed point as follows [9]. Consider the M layers of the coupled model and the relation where c UV and c IR (M) are the central charge of the pure system and that of the nontrivial fixed point, respectively. The effective central charge is then defined as In order to obtain the c IR (M), it is useful to recall the definition of the C-function and the differential equation satisfied by it [22,16]. To this end, we consider the vicinity of some critical theory S * in the theory space spanned by the coupling constants {g i }: and then specialize to our disordered O(n) model. The response of the action under the transformation z µ → z µ + α µ (z) can be expressed by the stress tensor T µν as As usual, a particular component of the stress tensor T zz = 1 4 (T 11 − T 22 − 2iT 12 ) and the trace of it 4T zz = T 11 + T 22 are denoted as T and Θ, respectively. As the RG flow dg i = β i dt occurs under the dilatation z µ → (1 + dt)z µ , the definition (75) and (76) leads to the relation: From T and Θ, one can consider three-types of the two-point functions: which are measured at the scale τ = ln(zz). The C-function is then defined as The conservation of the stress tensor ∂ µ T µν = 0 yields ∂zT (z,z) + 1 4 ∂ z Θ(z,z) = 0, and one can easily show d dτ After fixing the reference scale τ to τ = 0, the C-function depends only on the coupling constants g i . We introduce the Zamolodchikov metric on the theory space as From (77) and (82), we have For unitary theories, the positivity condition implies that the metric should be positive definite. The C-function is stable at the point with β i = 0 i.e. a fixed point and take the same value as the central charge as we can see from the definition (78)-(81). In our case, the perturbing field is quadratic in energy operator as in (29): and hence the metric (83) is This metric becomes negative in the limit M → 0, thereby violating the c-theorem. Note also that our normalization for the energy operator E is such that C I EE = 1. Substituting the beta function (39) into (84), we obtain From the definition of the effective central charge (74), we have Since we have the IR fixed points in ǫ > 0, this always increase under the RG flow as expected.

Discussion and Conclusions
Before concluding this paper, we discuss the relation to the known results and then comment on possible future directions. The RG flow in our disordered O(n) model near the disordered Ising point, which is inferred from the results (50), (58) and (69), is qualitatively similar to that in the random-bond q-state Potts model [8,20]. The crossover to an IR fixed point occurs only when the most relevant interaction, which is quadratic in energy operators as in (85), becomes relevant; the region correspond to n < 1 for the disordered O(n) model and q > 2 for the random-bond q-state Potts model. This can be understood by recalling how the two families of the pure models, namely, the O(n) models and the q-state Potts models coincide at the Ising point (n = 1, q = 2) [12,17].
The energy operator in the O(n) model is identified with the primary field φ 1,3 as we have stated in (13), while the energy operator in the Potts model is identified with the field φ 2,1 . These primary fields are, in a general minimal CFT, different objects as one expects from the fact that a four-point function of primary fields φ r,s is determined from an ordinary differential equation of order rs. Accordingly, in the vertex operator representation of the four-point functions, we should include two of the screening operators V − in the O(n) model and one of the other screening operator V + in the Potts model, respectively. However, the field φ 1,3 reduces to the other field φ 2,1 at the Ising point (the m = 3 minimal model) because of the equivalence (18). At the level of the scaling dimensions, the reduction can be summarized in Figure 7, in which the line of 2∆ 1,3 and the curve of 2∆ 2,1 intersect at the Ising point (κ = 3).
Up to one-loop order, when the deviation parameter ǫ from the Ising point defined in (28) is small and neglecting the structure constant (35) in the O(n) model can be justified, the beta function is determined only from the scaling dimension of the energy operator as in (39); thus similarity with the random-bond q-state Potts model is natural from Figure 7. By contrast, the reduction of the two-loop beta function (50), which involves four-point functions, to that of the q = 2 random-bond Potts model in the limit n → 1 serves as a non-trivial check of the vertex operator representations.
Remarkably, two lines of the IR fixed points lie the opposite sides of the Ising point; if we introduce a parameter the random fixed points of the disordered O(n) models and those of the random-bond q-state Potts models exist in the region ǫ κ > 0 and ǫ κ < 0, respectively. Another natural parameter ǫ P for the Potts models can be taken as α 2 + = 4 3 −ǫ P , where the IR fixed points emerges in the region ǫ P > 0 * . The definition (28), the relation (16) and α 2 + α 2 − = 1 from (A.3) lead to relations: Then the results for the scaling dimension of the energy operator in the random fixed point of the disordered O(n) model (59) and that of the random-bond q-state Potts models [8,20] can be summarized in one expression as where both results are measured from the Ising value 2∆ E = 1. Naturally, the scaling dimensions at these IR fixed points imply faster decays in the correlation functions: * It should be noted that yet another parameter ǫ defined as α 2 + = 4 3 + ǫ is used in [20].
∆ E < ∆ IR E as one expects with a general RG flow [27]. The ratio between the coefficients in ǫ 2 κ is 9 4 , which is the square of the ratio of the derivative of the line 2∆ 1,3 at the Ising point to that of the curve 2∆ 2,1 in Figure 7.
Contrastingly, the coefficients in the spin scaling dimensions are transcendental numbers. Using the relation (90), our result (71) for the disordered O(n) model is written as Then we quote, from ref. [20], the result for the random-bond q-state Potts model: where, in the second line, we have used the elliptic integral defined in (72). As a first observation, notice that both coefficients are consisting of the numbers K(cos θ) known as "singular values" of the elliptic integrals for which the modular ratios K(cos θ)/K(sin θ) belong to quadratic irrational numbers. The ratios here are √ 1 and √ 3 for the disordered O(n) model and the random-bond q-state Potts model, respectively♯. These modular ratios are exceptional in that the modular angles θ are commensurate to π, namely, θ = π/4 and θ = π/12.
The coefficient in (92) and that in (93) characterize the universality classes of the disordered models as well as of the corresponding coupled two-dimensional CFT's in the zero-layer limits. Now, it is noteworthy that they are simply related to the quantities originating from the specific structures of three dimensional lattices; the coefficient in (92) for the disordered O(n) model and that in (93) for the random-bond q-state Potts model correspond to the Watson's triple integrals [41,42] which represent the special values of Green's functions on the body centered cubic (bcc) lattice and on the face centered cubic (fcc) lattice, respectively. For our disordered O(n) model, to be concrete, the coefficient in (92) can be written as where the quantity in the parenthesis is the Watson integral for the bcc lattice. It would be interesting to see whether this coefficient can be derived, asymptotically at large scales, by a combinatoric analysis of coupled loop model on a lattice such as the one defined from the partition function (7). These remarks are also useful to give representations for the two coefficients in terms of rational integrals in order to see they are "periods" proposed in ref. [38]. Looking back the arguments, both coefficients are obtained not as products but as ratios between ♯ The two singular values have number theoretic origins and belong to the sequence of the values derived by Chowla-Selberg formula [40]; the relevant part of the formula is obtained from the functional equations of the Epstein's zeta functions on the quadratic number fields Q( √ −1) and Q( √ −3).
the complex Selberg integrals [36] and the non-symmetric extensions of them such as the one obtained in Appendix D, and thus it is not a priory clear that they are periods even if each of these Selberg-type integrals independently turns out to be a period. But if we start from the Watson integrals, they can be transformed into rational integrals [41], and hence the coefficients in (92) and that in (93), characterizing the two universality classes interconnected at the disordered Ising CFT, belong to, in fact, periods † †. Now, we list future directions in the following: (i) The fractal dimensions and SLE with the replica symmetry. We have seen that there is a line of IR fixed points for n * < n < 1. An immediate issue to be discussed is the fractal dimension D F of the loops in the IR fixed points. One may figure out the corrections for the Fisher exponent in (17) by perturbative calculation like the one performed in this paper. Related result has recently appeared for the random-bond qstate Potts model [62]. Further, the effective scale invariance on the line of the IR fixed points suggests the possibility of constructing one-parameter family of new SLE's. In order to describe the shape of the disordered loops, stochastic motion over the effective space enhanced by the replica permutation symmetry may be considered in analogy with the SLE driven by a composition of the stochastic motion in the real space and that in the internal space such as an SU(2) [50] or a Z N [51].
(ii) The strong coupling phase. We have seen that the coupled O(n) model has moved to strongly-coupled theory by the two-layer process in Figure 4-(b) in the region 0 < n < n * . The physical picture of the strong coupling region is important for the polymers in disordered environment [56]. Especially for this region, consideration on the theory space with the replica symmetry breaking (RSB) is desirable. Although the RSB-RG flow has been proposed in the random-bond q-state Potts model [21], numerical studies there supports the result of the replica symmetric fixed point rather than that of the RSB fixed point. We anticipate our model has more reason to favor the RSB situation due to formations of replica bound states boosted by the presence of the twolayer process. We should add that a large class of disordered models which includes strongly-coupled layered O(n) models has recently been discussed from the view point of the AdS/CFT duality [61]. Actually, the focus of their analysis is on large-n, while our strong coupling region lies at n ≈ 0; although a direct application seems to be difficult, this direction may be valuable in order to grasp the nature of strong coupling phases induced by disorder.
(iii) Non-local observables. We have studied the scaling dimensions of the fields corresponding to the local segments of loops. In the pure O(n) model, some results are known on non-local observables such as the area of loops [48] or the linking numbers around multiple points [49]. A challenging problem is the study of the off-diagonal pair-correlation between the global shapes of loops in the disordered model. The whole shape of a loop itself is, of course, difficult to handle. But extracting an essential part of such a correlation between non-local objects in a simple disordered system may † † More precisely, as suggested in the expression (94), the coefficient in (92) for the disordered O(n) model belongs to the 1/π-extended period ring P[1/π] defined also in ref. [38]. possibly provide us with the basis of other important problems. In chaotic systems, for instance, the off-diagonal correlation between periodic orbits, which are analogous to the disordered loops studied in this paper, is already a central issue in understanding the quantum level statistics [46,47]. It would be interesting if the possible relation between the disordered loops and periodic orbits in a chaotic system could be clarified.
In conclusion, we have studied the disordered O(n) models in two dimensions. We formulated the inhomogeneous loop model on a lattice. After adopting the replica method in order to take the disorder average, the models were mapped to the coupled layers of the homogeneous loop models. Some of the possible differences in the lattice formulation of the disordered models were argued to be irrelevant using the properties of the pure model in the dilute phase. Then we discussed the continuum limit of the disordered model assuming the identification between the energy operator E and the primary field φ 1,3 in the vicinity of the pure critical O(n) model. We considered the RG flow in the replica symmetric theory space. Since the most relevant interaction E (a) E (b) becomes marginal at the disordered Ising model n = 1, we used the epsilon expansion methods near n = 1 to perform a perturbative calculation. At one loop level, the RG flow suggests that there exists a critical value n = n * , the strong disorder region for 0 < n < n * and the random fixed point for n * < n < 1. These differences are largely controlled by the behavior of the structure constant C E EE which encodes the selection rules in the pure O(n) model. We interpreted the OPE intuitively and explained the picture of the strongly-coupled phase. Then we perform the two loop calculation near n = 1, where C E EE is O(ǫ). The beta function is equivalent to that of the Gross-Neveu model and thus checked the previously known result in the random-bond Ising model. The correction to scaling dimensions of the energy operator and the spin operator was calculated up to two loop. Finally, We calculated the effective central charge and see this increase under the RG flow.

Appendix A. The vertex operators in the critical Liouville field theory
The correlation functions in the O(n) model can be described by the vertex operators in the critical Liouville theory ‡. The basic idea is to deform the gaussian free field theory by putting the background charge (−2α 0 ) at the infinity. This is accomplished by coupling the scalar field ϕ to the scalar curvature R. The action of the theory is formally given by The corresponding central charge is determined from the OPE of the stress-energy tensor and turns out to be The correlation function of the vertex operators V α k = e iα k ϕ in this theory is non-zero only if the overall charge neutrality k α k = 2α 0 is satisfied. The screening vertex operators V ± = e iα ± ϕ are incorporated into the correlation functions in order to ensure the neutrality. The charges of the screening operators are chosen such that they are marginal, because otherwise the conformal invariance of the gaussian theory would be violated. This condition determine the charges: The respectively. In this paper, we fix our convention concerning the sign of the charges such that α 0 > 0 and α 0 < 0 describe the dilute (critical) and the dense (low-temperature) phase of O(n) model, respectively. The charge α r,s and its conjugate α r,s is defined as a linear combination of α + and α − : The scaling dimensions of the primary operators φ r,s is given by In the two-loop calculation in Section 4, we use the vertex operator representations of the four-point functions E(0)E(1)E(z)E(∞) 0 for (46) and σ(0)E(1)E(z)σ(∞) 0 for (67) and (68). In order to satisfy the neutrality, we should include two of the screening operator V − in the vertex correlation functions. The normalization factors in the fourpoint functions are determined in a decoupling limit, such that the structure constant ‡ Although it is often referred as the "Coulomb gas representation" [17,20], we call it the "vertex operator representation" in the critical Liouville theory [15] in distinction with the "Coulomb gas method" [18].
C I EE is fixed to unity. The value can be calculated by the complex Selberg integral [36] as well as by the method described in Appendix D; the results are the same for both types of the four-point functions: where ̟ = Γ 1 4 2 / 2 3/2 π 1/2 = 2.62205755 · · · is the lemniscate constant.
Appendix B. Some details of the integral for the two-loop beta function and the scaling dimension of the energy operator In this Appendix, we thoroughly use the result of Appendix D. We read off the exponents in (D.1) from the vertex operator representation (46). 2a = 4α 13 α 1,3 , 2b = 8α 13 α 1,3 , 2a ′ = 2b ′ = 4α 13 α − , 2f = 4α 1,3 α − , 2g = 4α 2 − For convenience, we interchange the values of a and b, and those of a ′ and b ′ (the latter is trivial, since here a ′ = b ′ ). These interchanges are possible by an obvious symmetry between 0 and 1 in the integral (D.1).
As a result, we have a = 2b = −2f = 4α 13 α 1, . From these parameters, we determine the integral by the "scattering amplitude" formula (D.7). The scattering matrix M are given by (D.8) and read, From (D.20), the initial and the final state basis in (D.7) are given by where, J ± l (l = 1, 2, 3) are real triple integrals defined in (D.9)-(D.14), and the symbol ... is explained in (D. 19). We should discuss the order of these integrals J ± l in ǫ; for this purpose, we can use the series representation derived in Appendix D.3. From (D.22), the prefactors γ ± l are calculated as, Naively, we expect the leading terms in J l = γ l S l is the same order in ǫ as the γ l , since all the triple series S l contains the (i, j, k) = (0, 0, 0), which is unity and hence O(1). Actually, it is not the case for J + 1 ; the series S + 1 is O(ǫ) because of the non-trivial cancellation between the constant terms.
Since this observation makes the crucial point in the calculation, we describe the calculation of J + 1 in detail. First, from (D.23), we have Denoting the summand in (B.4) as S i,j,k , we observe, from the definition of the Pochhammer symbol (x) k = x(x + 1) · · · (x + k − 1), that S i,j,k gets a factor of O(ǫ) whenever each of two conditions {i ≥ 2, i + j ≥ 2} is satisfied. Thus, we let where the leading order of a 1 , a 2 , a 3 and a 4 , a 5 are O(ǫ 0 ) and O(ǫ 1 ), respectively. To be specific, at the accuracy of O(ǫ), we are left with the following: where the values of generalized hypergeometric functions 3 F 2 at unity are used. The arguments, which is always taken at unity, are suppressed in the first equalities in (B.7), (B.8) and henceforth. Now, we see the aforementioned cancellation at O(ǫ 0 ), and get This is an example of identities known as the "contiguity relation" between hypergeometric functions. Note that this type of the leading order cancellation occurs only for S + 1 . We note that the cancellations at this order guarantees a consistency in the RG scheme. Resulting O(ǫ) term can be evaluated numerically, using where ψ(x) is the di-gamma function §. The other terms contributing O(ǫ) in S + 1 are (B.11) § Actually, there have recently been extensive studies on the closed form evaluation of the expansion of the hypergeometric series in algorithmic approach, for example [39]. However, to our knowledge, the desired expansions here have not been published.
To evaluate this double series, it is useful to note the following identity: This holds because (a + b) k has the same binomial expansion as (a + b) k . By comparing the order δ terms of the following expression, one can express the right hand side of (B.11) in terms of a single series: where F δ (...) denotes the coefficients of the order δ terms in the 3 F 2 function at unity. In this case, We use this notation hereafter. Now, collecting the O(ǫ) parts of (B.6)-(B.8) and (B.14), we obtain This completes the calculation of J + 1 . Since the linear relations (D.24) are solved as the knowledge of the other two bases are sufficient in the formula (D.7). Actually, we can calculate J − 1 and J − 3 in similar manners. Let us define constants D and F as follows: where, again, (B.3) has been used for γ − 1 and γ − 3 . Then, D and F are given by where Q 2 denotes well-converging double series which is in good agreement with another value Γ(1/4) 4 /2π = 4̟ 2 = 27.50074327208 · · ·. Assuming the latter value in (B.28) and using (A.7), we obtain parameters a ↔ b and a ′ ↔ b ′ (as in Appendix B), we have a = 2α 2 13 − 2∆ E = 1 2 + 2ǫ, b = 2α p−1,p α 13 = − 1 4 + ǫ, a ′ = f = 2α 13 α − = − 3 2 + 2ǫ, b ′ = 2α p−1,p α − = 1 4 − ǫ, and g = 2α 2 − = 1 2 − 2ǫ for K 2 ; we have almost the same set of the parameters as K 2 but a = 2α 2 13 + 1 − 4∆ E = 1 2 + 6ǫ for K 3 . We substitute these into (D.8) and get, both for K 2 and K 3 , the same matrix: (C.1) where the notation Ω in (C.9) is used to express both K 2 and K 3 in parallel. Actually, only the combination U 2 + W 2 is necessary for the disordered model (M → 0). Substituting the parameters in (C.2) and (C.4) into (D.22) and (D.23), we obtain (C.12) We obtain numerically which is compatible, within error bar (±0.05%), with the value Γ( 1 4 ) 12 8π 6 = 64 ̟ 6 π 3 = 670.78 · · · . (C.14) Unfortunately, the convergence of the triple series U is slow, and thus our numerical accuracy is not good here. Nevertheless, we assume, for U 2 + W 2 , the value in (C.14).

Appendix D. Integrals
In our calculation of the RG functions, we should deal with a multiple integral over C 3 : Since this form of the integral comes from the correlation functions of the vertex operators (or, more plainly, from the interaction between charged particles in a twodimensional plane), it seems to be ubiquitous in physics and mathematics. For example, the integral can be interpreted as a six-particle closed string amplitude [35]. The integral in a special, symmetric case of parameters (a = a ′ , b = b ′ and f = g) is well studied in the context of twisted cohomology, and known as the "complex Selberg integral" [36]. What we need, however, is the formula in a non-symmetric case, when doing perturbation theory around a conformal fixed point. In this respect, a formula for two variables was used in the study of random-bond Potts model by Dotsenko et al. [20]. We extend their results and derive a formula for (D.1) in a systematic way. The formula, obtained in (D.7), takes form of a scattering amplitude.
Appendix D.1. Regularization of the one-dimensional intervals We encounter with strong algebraic singularities that make the multiple integrals divergent. For this reason, we consider an analytic continuation in the parameter of the integrals. In order to keep the discussion clear, it is helpful to make the way of the analytic continuation explicit. The analytic continuation is achieved by the use of a "regularization" of the intervals [36,37], which we now describe using the following one-dimensional simple example.
Consider an integral on a real interval Here, as in figure D1, δ A and δ B are positively oriented circles of radius δ which have centers A and B, and start at A + δ and B − δ, respectively. By replacing the interval [A, B] by the regularized one "reg [A, B]" and taking the limit δ → 0, the value of the integral (D.2) remains same for Re p > −1 and Re q > −1 since the contributions from two additional circles vanish, and is now finite also for Re p < −1 or Re q < −1 unless p ∈ Z or q ∈ Z. In the latter case, adding two circles corresponds to the subtraction of infinite quantities, and the resulting finite value is what is known as the "Hadamard finite part" of the integral . The definition here is natural in the following sense. Starting with the multiple integral (D.1) and introducing an infinitesimal imaginary part, we shall reach an iterated integral in which each integral is on the regularization of the interval rather than on the usual real interval. This regularization comes from the pairing of the paths on the upper and the lower half planes both of which detour the branch points. All the one-dimensional integrals in this Appendix should be regarded as the integral over the regularization of the interval. in (D.1): |z| 2a = (z 2 1 + z 2 2 ) a is now defined on (z 1 , z 2 ) ∈ C × C. Then the path of z 2 is rotated by the angle π/2−2η, without hitting any singularity, for a positive infinitesimal number η. By introducing a new variable z 0 , we rewrite |z| 2a as follows: where the notation z ± = z 1 ± z 0 is introduced. In the following, we call {z + , u + , v + } and {z − , u − , v − }, respectively, holomorphic and antiholomorphic variables. Using a notation X η = η[X + − X − ] for a quantity X, we decompose the integral (D.1) as Here, J + and J − are weakly dependent, through an infinitesimal number η, on the antiholomorphic (−) variables and the holomorphic (+) variables, respectively. If we fix the holomorphic variables first, the dependence of J − on the holomorphic variables determines relative positions of the integration paths and the two branch points 0 and 1 on the complex planes of the antiholomorphic variables. For this relative positions, we can observe there are two very distinct cases. Figure D2. Deformation of the paths on the antiholomorphic planes. Both figures correspond to the sequence of the holomorphic variables z + < v + < u + . (a) Before the deformations.
(b) After the deformations. The five symbols A, C, D, E and F are assigned to each regularization of the interval. The case with v − < u − is omitted here, and thus the regularization B in Figure D3 does not appear.
The first case is trivial, and corresponds to a choice of the holomorphic variables such that at least one of the z + , u + and v + lies outside the interval [0, 1]. Then, according to (D.6), one can find at least one integration path on the antiholomorphic plane such that all of the points 0, 1 and the other paths projected are located on the same side of the plane. Therefore, the path can be contracted to a point by adding a semicircle of infinite radius and the corresponding integral should vanish.
There is, however, the second case in which all of the z + , u + and v + belong to the interval (0, 1). On the antiholomorphic planes, this means all the paths of z − , u − and v − intersect the segment [0, 1] in the same sequence as z + , u + and v + lie on the interval (0, 1). Take a sequence: z + < v + < u + as in Figure D2-(a), for instance. Now, we deform the paths so that the resulting paths encircle the point 1, the real axis from the upper and the lower side ( Figure D2-(b)). Since the path nearest to the point 1 is that of u − , we first deform it, then that of v − , and finally that of z − . Note, in deforming the path of, say v − , the fixed variable u − is projected on the v − -plane as a branch point.
As a result of the deformation, we have a pairing of the paths in the upper and the lower half planes for each variable z − , u − and v − . Each paring of the path can be decomposed into a sum of the integral over the regularizations of intervals (see Appendix D.1 for the definition of the regularizations and ref. [36] for the summation of them).
Taking the presence of the branch cuts into account, we now attach the appropriate factors on these regularizations. Since each factor comes as a difference of two phase factors on the upper and the lower half plane, it takes the form of a sin-function. In the following, we assign the numbers {1, 2, 3} for the sequences of variables {(z ± < v ± < u ± ), (v ± < z ± < u ± ), (v ± < u ± < z ± )}, respectively.
Consequently, from (D.6), we obtain a scattering-type formula ) and M ± in (D.24) in this paper, and the intersection matrix and the monodromy invariant hermitian form in [37] was elaborated. Figure D3. Pairings of the paths (the regularizations of the intervals) that contribute the matrix elements of M. Each branch correspond to the regularization (see also Figure D2 for the branches A, C, D, E and F). The branches B and C, for instance, correspond to the regularization in the v − -plane with the sequence v − < u − and u − < v − , respectively. Each regularization has a certain factor on it due to the presence of the branch cuts. and where we have used the notation The matrix elements of M are conveniently understood if we draw the tree-diagrams which show the sequences of the variables under the process of the deformation ( Figure  D3).
The notation e(x) = exp(πx) and s(x) = sin(πx) are used. The semicircular shape represents the regularization of the interval (see Figure D1). The dotted curves indicate the integrations over the variable z, while the bold curves indicate that of u or v. The vertical lines indicate the relative positions of the branch cuts induced for the other variables.