Ab initio holography

We apply the quantum renormalization group to construct a holographic dual for the U(N) vector model for complex bosons defined on a lattice. The bulk geometry becomes dynamical as the hopping amplitudes which determine connectivity of space are promoted to quantum variables. In the large N limit, the full bulk equations of motion for the dynamical hopping fields are numerically solved for finite systems. From finite size scaling, we show that different phases exhibit distinct geometric features in the bulk. In the insulating phase, the space gets fragmented into isolated islands deep inside the bulk, exhibiting ultra-locality. In the superfluid phase, the bulk exhibits a horizon beyond which the geometry becomes non-local. Right at the horizon, the hopping fields decay with a universal power-law in coordinate distance between sites, while they decay in slower power-laws with continuously varying exponents inside the horizon. At the critical point, the bulk exhibits a local geometry whose characteristic length scale diverges asymptotically in the IR limit.


Introduction
According to the AdS/CFT correspondence [1][2][3], quantum field theories are dual to gravitational theories in a spacetime with one higher dimension. The extra dimension in the bulk can be interpreted as a length scale in the renormalization group (RG) flow [4][5][6][7][8]. However, the flow generated along the extra dimension in holography is different from the conventional RG flow because the bulk theories are in general quantum theories. In order to make the connection between RG and holography precise, one has to introduce quantum RG (QRG) [9][10][11]. According to QRG, general D-dimensional quantum field theories are equivalent to (D + 1)-dimensional theories that include quantum gravity, provided that the quantum field theories can be regularized covariantly in curved spacetime [10,11].

JHEP08(2015)107
Conventional renormalization group describes the flow of the couplings (sources) generated by coarse graining defined in the space of all sources allowed by symmetry. In QRG, only a subset of the sources is included. The sources in the subset are promoted to quantum mechanical operators, and a quantum Hamiltonian governs the evolution of the dynamical sources under scale transformations. Quantum fluctuations in the RG flow precisely capture the effect of other couplings that are not explicitly included. In the bulk, this amounts to the fact that composite operators are generated after lowering the position of the UV boundary by integrating out bulk degrees of freedom [7,8]. Intuitively, sources become quantum mechanical because high-energy modes act as fluctuating sources from the point of view of low-energy modes [9][10][11][12].
Despite the formal mapping between general quantum field theories and gravitational theories, classical and local geometries are expected to arise only for a special set of quantum field theories [13][14][15]. In the presence of a large number of degrees of freedom, the quantum fluctuations of the RG flow become weak, which allows one to use a classical description. However, it is much more non-trivial to have locality in the bulk. Generally, there exist infinitely many non-local operators, and locality arises only when the sources for the non-local operators are suppressed. In QRG, the condition for the emergence of local geometry translates into stringent constraints on the beta functions of the quantum field theory [16,17]. From the constraints, one can try to find quantum field theories that exhibit a local geometry in the bulk. On the other hand, one can use varying degrees of locality that emerge in the bulk as a diagnostic that differentiates one phase from another. The bulk locality serves as a useful order parameter for characterizing phases of matter using geometry, which is one of the primary goals of applying holography to condensed matter systems [18][19][20][21].
In this paper, we apply QRG to a three-dimensional U(N ) vector model regularized on a Euclidean lattice, which can be viewed as a quantum Bose-Hubbard model of N components in two space dimensions and an imaginary time. It supports an interacting field theory at the critical point between the insulating (gapped) phase and the superfluid (long-range ordered) phase. The full bulk equations of motion derived from QRG are solved in the large N limit on finite size lattices. From a finite size scaling analysis, we show that different phases indeed exhibit different degrees of bulk (non-) locality in the thermodynamic limit.
Here is the outline of the paper. In section 2, we start with the U(N ) vector model for N complex bosons defined on a D-dimensional lattice. In this model, the general single-trace operators can be written as N a=1 φ * ia φ ja , where φ ia is a complex boson with flavour a defined at site i. The single-trace operators involve only one flavour contraction and are bi-local in space [22][23][24]. The set of single-trace operators has a special status in that all other singlet operators can be constructed as composites of single-trace operators. They describe the kinetic term that connects pairs of sites through hopping. In the vector model, all single-trace operators are quadratic, and it is not enough to include only singletrace operators to describe the phase transition from the insulating (gapped) phase to the superfluid (symmetry-broken) phase. Therefore, we also include a quartic doubletrace operator. We treat the single-trace hopping terms as deformations to the on-site Figure 1. A schematic holographic phase diagram in the thermodynamic limit. ∆m 2 is the mass which tunes the system from the superfluid phase to the insulating phase. z is the extra dimension in the bulk, which corresponds to a logarithmic length scale in the renormalization group. In the insulating phase, the space gets fragmented beyond the fragmentation scale, z F . In the superfluid phase, non-locality emerges beyond the horizon scale, z H . At the critical point, the bulk exhibits a local geometry.
action that includes the single-trace mass term and the double-trace interaction term. In the conventional RG, a series of coarse graining generates an infinite tower of multi-trace operators [25]. In QRG, one only keeps track of single-trace deformations at the expense of promoting the sources for the single-trace operators to quantum variables. We derive a (D + 1)-dimensional bulk action whose degrees of freedom are scale dependent bi-local hopping fields t ij (z) and their conjugate variables t * ij (z) , where z is the logarithmic length scale. The bulk theory describes the quantum mechanical RG flow in the space of single-trace operators. In the large N limit, the path integral in the bulk can be replaced by saddle point equations whose solution determines the classical geometry in the bulk. In section 3.1, an analytic solution to the bulk equations of motion is found in the deep insulating phase by treating hoppings as small perturbations compared to the on-site terms. Away from the deep insulating phase, we resort to numerical solutions. Section 3.2 and 3.3 outline the numerical scheme. In section 4, the numerical solutions obtained for three dimensional lattices with linear sizes 3 ≤ L ≤ 13 are presented. Then we extract the behaviour of the solution in the thermodynamic limit from finite size scaling. Figure 1, which is the key result of the paper, summarizes the numerical solution in the thermodynamic limit, which is suggested from the finite size scaling analysis. Here ∆m 2 represents the on-site mass term, which tunes the system away from the critical point either to the insulating phase (∆m 2 > 0) or to the superfluid phase (∆m 2 < 0) depending on its sign. Although only the nearest neighbor hoppings are turned on at the UV boundary at z = 0, further neighbor hoppings are generated inside the bulk with z > 0. In the insulating phase, the hopping fields t ij (z) decay exponentially not only in |i − j| but also in z. As a JHEP08(2015)107 result, the lattice is fragmented into decoupled sites in the large z limit. While the scale z F beyond which fragmentation occurs remains small in the deep insulating phase, it diverges as the critical point is approached in the thermodynamic limit. In the superfluid phase, on the other hand, a new scale emerges. Instead of fragmentation, the system loses locality beyond a critical scale z H , which is called the horizon. At the horizon, the hopping fields decay in a universal power-law in |i − j|. Inside the horizon with z > z H , the hopping fields decay with a slower power-law with continuously varying exponents. The horizon scale z H diverges as the critical point is approached from the superfluid side. At the critical point, the bulk exhibits a local geometry with a characteristic length scale that remains finite for finite z and diverges only in the IR limit. In summary, the insulating phase, the superfluid phase and the critical point exhibit ultra-local, non-local and local bulk geometries, respectively. In this sense, locality serves as a holographic order parameter that characterizes different phases of matter. The goal of section 4 is to justify the holographic phase diagram in figure 1. For related works on the holographic description of the vector models, see refs. [22][23][24][25][26][27][28][29].

Holographic action for vector model on a lattice
We consider a bosonic vector model in D-dimensional Euclidean space, where φ refers to the complex boson field with N flavours, |φ| 2 ≡ φ * · φ, m is the mass of the bosons, and λ is the quartic coupling. We regularize the continuum theory on a D-dimensional hypercubic lattice, 3) Here i, j, p, q run over the D-dimensional Euclidean lattice for the N component bosonic fields φ ia with a = 1, · · · , N . S 0 represents on-site terms, and S 1 includes the hopping term and general quartic interactions. To guarantee the action is real, we imposet ij =t * ji andJ ijpq =J * jiqp . Here we consider the case wheret ij is real. We could maket ij complex by turning on a background gauge field. The background gauge field in general breaks the time-reversal symmetry, if we go to Minkowski space by Wick-rotating one of the Euclidean directions into real time. For generalt ij andJ ijpq , the model has the U(N ) symmetry. The symmetry is enhanced to O(2N ) whent ij 's are real andJ ijpq ∝ δ ij δ pq . The action is proportional to N in the large N limit with fixed m 2 , λ,t ij andJ ijpq .

JHEP08(2015)107
There is a redundancy in the parameters of the action, and the action is invariant under (2.5) By shiftingt ii andJ iiii , S 0 could have been entirely absorbed into S 1 . Here S 0 is explicitly singled out because we will use S 0 as 'the reference theory' and treat S 1 as a deformation. In other words, the renormalization group (RG) flow will be defined in terms of the flow of the parameters in S 1 which is defined with respect to the fixed S 0 with m 2 , λ > 0. It is emphasized that one can describe not only the insulating (gapped) phase but also the superfluid (symmetry broken) phase becauset ij can be arbitrarily large to support a Mexican hat potential for φ. The freedom to choose different m 2 and λ amounts to choosing different RG schemes, which does not affect physical observables. We derive a holographic theory for the lattice action in eq. (2.2) using the QRG scheme. We start by writing the partition function: where the measure is given by ia and V is the number of sites in the lattice. The first step is to remove multi-trace operators in S 1 by promoting the sources for the single-trace operators into dynamical variables. Using the representation of the Dirac-Delta function [9] we introduce a pair of complex conjugate link fields, t ij and t * (0) ij for every ordered pair of sites i, j to rewrite the partition function as with i, j running over all sites and ij . The presence of the double-trace operator is also important because it makes the hopping fields t (0) ij genuinely dynamical variables. Note that φ is coupled to dynamical sources t (0) ij which connect all possible sites. Therefore the new action is defined on a globally connected network rather than a fixed lattice. The globally coupled network with dynamical hoppings does not have a fixed lattice structure, and i, j should be regarded as indices for 'events' rather than fixed coordinates. Since there is no fixed lattice, there is a freedom to relabel each event differently, which corresponds to a discrete local coordinate transformation. So we perform a discrete coordinate transformation, where site i is displaced by a site-dependent D-dimensional shift vector, N (0) i [10,11,26]. Since the full set of events must map to itself, the transformation should form cyclic permutations of the events. This is illustrated in figures 2(a), (b). After the local coordinate transformation, S 0 is manifestly invariant and the hopping term is transformed to where the shifted hopping field is given by The point is that the change introduced by the local coordinate transformation can be compensated by a transformation of the hopping fields because t (0) ij is dynamical. Therefore, the dynamical hopping fields play the role of a dynamical metric in the continuum limit, which JHEP08(2015)107 guarantees invariance under a set of D-dimensional local coordinate transformations in the continuum limit. The set of permutations generated by the shift includes D-dimensional local coordinate transformations, which don't necessarily preserve the coordinate volume locally. An example is given in figure 2(c). The fact that the hopping fields play the role of a metric is expected because the physical distance between two sites is determined from the strength of the hopping between them: the larger the hopping between two sites, the shorter the physical distance. Now we generate RG flow by coarse graining the system in real space. The original field φ is split up into high-and low-energy fields, where the low-energy field has a mass slightly larger than m 2 . The missing fluctuations are carried away by the very massive high-energy fields. Then, the high-energy modes are integrated out to generate quantum corrections for the low-energy modes [30,31]. For this, we introduce an auxiliary field Φ with an arbitrary mass µ, Now we rotate the physical field and the auxiliary field into high-and low-energy fields, i is the 'lapse', which is a site-dependent constant that controls the local speed of coarse graining [10,11,32]. Here φ ′ is the lowenergy mode with mass m 2 e 2α (1) i dz , which is larger than m 2 by an infinitesimal amount O(dz), andφ is the high-energy mode with a large massμ i 2 e 2α (1) i dz . After rescaling the (1) i dzφ i and renaming φ ′ i → φ i , we rewrite the partition function as (2.13)

JHEP08(2015)107
Integrating out the high energy modeφ, we obtain the renormalized action for the low energy mode, 14) where quantum corrections are kept only to the linear order O(dz) (see appendix A for details). In the renormalized action, multi-trace operators have been generated for the low energy field. Under the conventional (classical) Wilsonian RG procedure, one would apply the same coarse graining procedure to eq. (2.14) to generate even higher-trace deformations.
In QRG, we remove the multi-trace deformations before repeating the coarse graining. This way, one can keep only single-trace deformations at each step of coarse graining by making the sources for the single-trace deformations fluctuating. Just as before, we introduce new link fields t (1) ij and t * (1) ij to remove the multi-trace deformations, Now we can repeat the coarse graining procedure as before. The coordinates of the φ i 's can be shifted by a new set of discrete translations, N (1) i . Then φ is split into lowenergy and high-energy modes where the low energy mode at site i has a mass m 2 e α (2) i dz . After rescaling the fields, we integrating out the high energy modes to generate quantum JHEP08(2015)107 corrections which include multi-trace operators for the low energy fields. We introduce another set of auxiliary fields t (2) ij and t * (2) ij to remove the multi-trace operators. From here the pattern is clear. After doing this Γ times, the partition function is written as Here t i } are chosen independently at each step of coarse graining. Finally, we take the continuum limit in the RG direction by introducing z = ldz in the limit dz → 0 with z * ≡ Γdz fixed. In the limit where z is continuous, t The partition function is written in terms of the path integral for the scale dependent hopping fields and their conjugate fields, (2.17) Here S UV , defined in eq. (2.8), is the UV boundary action for the dynamical source at z = 0. S bulk is the bulk action, with (2.19) which is defined at the IR boundary z = z * . z * is the scale at which we stop the coarse graining procedure, which can be taken to be infinite.
The new action is written in terms of the fields t ij (z) and t * ij (z) which are bi-local within the original D-dimensional lattice and local in the z direction. If one interprets z as 'time', the action describes a system of bi-local hopping fields under a Hamiltonian evolution in imaginary time. t ij and t * ij correspond to annihilation and creation operators respectively [9], and the last term in eq. (2.18) describes the process where two hoppings merge to become a longer-range hopping. This is illustrated in figure 3. In the path integral language, all RG paths are summed over for the single-trace operators. The weight for each RG path is determined by the bulk action, S Bulk .
In the bulk, the lapse and shift vector parameterize local coordinate transformations. The partition function does not depend on {α i (z), N i (z)} because different choices of lapse and shift merely correspond to different renormalization group schemes [10,11]. As a result, the local coordinate transformation in the bulk is a gauge symmetry. Unlike in the continuum, the action contains terms that are non-linear in the shift. This is because the discrete shift is a finite transformation. When t ij (z) vary slowly in i, j, one can take a continuum limit along the D-dimensional directions by writing where r 1 = ar i , r 2 = ar j in the limit that the lattice spacing a is small. Then, the shift can be performed infinitesimally, N i (z) → n(r, z)dz, and one can drop the terms that are JHEP08(2015)107 non-linear in dz. In this case, the action is linear in n(r, z) and α(r, z), which generate (D + 1)-dimensional coordinate transformations in the bulk. Also, in the continuum limit, the bi-local fields can be represented by an infinite set of fields with arbitrarily large spin.
Here we don't take the continuum limit and proceed with the lattice action.
In the absence of the quartic term (λ = 0), the theory has a larger symmetry generated by [28,29] where V ij is a unitary matrix in the space of all sites in the lattice. The enlarged symmetry is the source of the higher spin gauge symmetry for the free theory with λ = 0. In the interacting theory with λ = 0, the symmetry is broken down to the local coordinate transformation. One might keep the higher spin symmetry in the bulk action even when λ = 0 by absorbing the quartic term into the UV boundary action S UV . In that case, there is no λ in the bulk, and the quartic term is implemented only through the alternative boundary condition determined by S UV [33]. However, this can be done only in the strict large N limit where the fluctuations of the dynamical hopping field induced by S UV are negligible. What is less satisfying in this description for the interacting theory is the fact that 1/N corrections are singular because the path integral for φ is ill-defined without the quartic term in S 0 for fluctuating hopping fields. Therefore we choose to include the quartic term in the bulk for general N , which breaks the higher spin symmetry.
Because the fields t ij , t * ij are singlets under U(N ) rotations, the action in eq. (2.17) is proportional to N in the large N limit. Therefore, the partition function can be obtained by the saddle point approximation in the large N limit. To that effect, all we need to solve are the equations of motion, which we focus on from now on. Although t ij and t * ij are complex conjugates of each other in eq. (2.17), their saddle points configurations do not have to satisfy that condition. This is because the amplitude and phase of t ij can take complex values at the saddle point. Therefore, we treat t ij and t * ij as independent complex fields in the saddle point equation.
In the gauge with α i (z) = 1 and N i (z) = 0, the equations of motion for the bulk fields are given by

JHEP08(2015)107
subject to two boundary conditions which are derived by extremizing the boundary actions at z = 0 and z = z * , respectively. In the equations of motion, sub-leading terms are dropped in the large N limit. From the IR boundary action, one can readily write the IR boundary condition as where O S IR denotes the expectation value of O evaluated with respect to S IR at z = z * . This relationship must hold anywhere inside the bulk, since the RG process can be stopped at any point. Therefore, the on-shell value of t * ij (z) coincides with the 2-point correlation functions at scale z, which are completely determined by t ij (z).
In the end, one has to solve the first order differential equations for t ij and t * ij with one set of boundary conditions imposed at the UV boundary and the other imposed at the IR boundary. It is the IR boundary condition that imposes the constraints that the expectation values of operators (t * ij ) have to satisfy for a given set of sources (t ij ). In particular, the vacuum expectation value has to satisfy the Ward identity associated with the discrete coordinate transformation, where t ′ is defined in eq. (2.11). Since this equation has to be satisfied at all z, it becomes , which is the discrete version of the energy-momentum conservation imposed by N i (z). This condition is automatically implied by the IR boundary condition even in the fixed gauge with N i (z) = 0. It is noted that t * ij cannot be arbitrary at the IR boundary, but it has to be the actual vacuum expectation value computed from the IR boundary action, S IR .

Solutions to saddle point equations
In this section, we examine the solution to the equations of motion both analytically and numerically. We assume that translational symmetry is present, in which case t ij and t * ij depend only on |i − j| ≡ |r i − r j |. Also, since we are maintaining the full SO(D) and D-dimensional inversion symmetries at each scale z, all t ij (z) and t * ij (z) are real. We focus on the three dimensional cubic lattice of linear size L with periodic boundary conditions.
To uniquely determine the solutions we need to impose the boundary conditions in eqs. (2.22) and (2.24). At the UV boundary, eq. (2.22) imposes the Dirichlet boundary condition for t ij (0) =t ij in the absence ofJ ijpq . To impose the IR boundary condition in eq. (2.24) we compute S IR [t ij (z * )] by introducing a Hubbard-Stratonovich field σ which JHEP08(2015)107 satisfies the saddle-point equation in the large N limit (see appendix B), . Then the IR boundary condition becomes .

Analytic solutions in the deep insulating phase
For m 2 ≫t ij and m 4 ≫ λ, the quadratic on-site action dominates. All correlation functions decay exponentially, resulting in the insulating phase. In this case, one can easily obtain an analytic solution to the equations of motion as a perturbation series int ij m 2 and λ m 4 . To order 1/m 6 , the solutions are Both t ij (z) and t * ij (z) decay exponentially as one moves towards the IR boundary. This is expected since bosons are localized at low energies. It is noted that t * ij for i = j decays as e −2z with increasing z. This feature holds generally even when t ij is not small. This is due to the fact that the transformation in eq. (2.12) followed by the rescaling by e −dz implies φ i (z + dz) = e −dz φ i (z) up to the contribution from the auxiliary field introduced at each step of coarse graining which has zero correlation length. This implies that In the insulating phase a large simplification can be made for the IR boundary condition. Since both t ij (z) and t * ij (z) decay exponentially in z except for i = j, we can approximate them as for a sufficiently large z * . In this case, the IR boundary condition reduces to a single-site problem,

Numerical solutions
Since it is hard to solve the full equations of motion analytically in general, we resort to numerical solutions for the superfluid phase and the critical point. We solve the coupled differential equations using a spectral method. The details of the numerical method can be found in appendix C. Numerical solutions are obtained for the cubic lattice with linear sizes 3 ≤ L ≤ 13. We set λ = 1, m 2 = 25, andt ij = 0 except for the on-site and nearest neighbor hoppings. The nearest neighbor hoppings are set to be 1, and the on-site 'hopping' t ii is tuned to change the physical bare mass, to drive the phase transition from the insulating phase to the superfluid phase. The physical massm 2 can be either positive or negative whereas the mass m 2 of the reference action in eq. (2.3) is fixed to be positive. The reason for choosing a positive m 2 is to make sure that only short-distance modes are integrated out in the coarse graining procedure.
Although there is no genuine symmetry breaking phase transition in finite systems, features of the superfluid phase and the critical point can be still identified from the behaviors at intermediate length scales before the finite size effect takes over. Henceforth, we will keep using the terminology "superfluid" to refer to systems in parameter regimes which become real symmetry broken states in the thermodynamic limit. Although only the nearest neighbor hopping is turned on at the UV boundary, all further neighbor hoppings are generated inside the bulk. From now on, we will use the notation t (l,m,n) , t * (l,m,n) to denote the specific fields t ij , t * ij with r j − r i = (l, m, n) in units of the lattice spacing. For example, the three nearest neighbor hopping fields are represented by t (0,0,1) (z), t (0,1,0) (z), t (1,0,0) (z).

Universal IR boundary condition in finite systems
In figures 4 and 5, we plot t ij (z) and t * ij (z) for L = 3 with z * = 8. For each (l, m, n), we display multiple curves at 25 ≥m 2 ≥ −15 that cover from the deep insulating phase to the deep superfluid phase across the critical point, which is atm c 2 = 5.49454 in the thermodynamic limit. Since the curves look very similar away from critical point and they start to change rapidly only in its vicinity, we choose the spacing ofm 2 to be smaller near the critical point and larger away from it. Here, we notice something interesting. No matter what the values ofm 2 are, all t ij (z) and t * ij (z) with i = j approach zero in the large z limit. In other words, the system flows to the on-site problem in the deep IR not only in the insulating phase but also in the deep superfluid state. The universal insulating (on-site) behavior in the deep IR limit is due to the finite size effect. Form 2 <m 2 c , t ij increases with increasing z for small z, which is consistent with the expected RG flow in the superfluid phase where the hopping terms become more important at larger length scales. However, the RG length scale eventually becomes greater than the system size. Then the finite size effect becomes important, and the system "flows" back to the insulating behavior in the deep IR. This has to do with the fact that all phases are adiabatically connected to the insulating phase in finite systems.  At first, the finite size effect appears to be detrimental to our goal of characterizing different phases of matter in terms of different IR geometries. However, one can still study the superfluid phase and the critical point by focusing on the behavior of t ij , t * ij in the intermediate range of z before the finite size effect becomes dominant. Moreover, we can use the finite size effect to our advantage. The universal insulating behavior in the IR limit allows us to use the single-site IR boundary condition in eq. (3.6), which is much simpler than the full IR boundary condition in eq. (3.2). This has a significant implication. A practical difficulty in applying QRG to strongly coupled field theories is that one has to impose the IR boundary condition dynamically. This is done either by extremizing the bulk action in the z * → ∞ limit or imposing eq. (2.24) at a finite z * . If one uses the latter scheme, one has to know S IR . Although this does not pose any problem for the present JHEP08(2015)107 vector model which is exactly solvable in the large N limit, this is in general a difficult task for strongly coupled theories (such as matrix models). However, if one can always impose the insulating boundary condition for finite systems one only needs to know the solution to the on-site problem, which is much easier. This allows one to find solutions for finite systems without knowing the IR behavior of the infinite system a priori. Then, one can extract the behaviors in the thermodynamic limit through finite size scaling. This is the strategy we will employ in the rest of the paper.
In order to guarantee that systems with finite sizes flow to the insulating phase, one has to choose z * to be sufficiently large. The smallerm 2 is, and the larger L is, z * has to be larger to ensure that the system flows to the deep insulating state driven by the finite size effect at the IR boundary. This is because it takes longer RG 'time' before the finite size effect takes over in deeper superfluid phases and larger lattices. We observed that for all values ofm 2 we took within the range 25 ≥m 2 ≥ −15 and all lattice sizes 3 ≤ L ≤ 13, z * = 8 is sufficient. For larger L, it is expected that z * should grow logarithmically in L JHEP08(2015)107 The correlation function at the UV boundary t * ij (z = 0) = 1 N (φ * i · φ j ) UV between farthest sites (l, m, n) ≡ r j − r i = (4, 4, 4) plotted as a function ofm 2 for L = 9. Although there is no true singularity due to finite size effect, the correlation function exhibits a rather sharp kink at the critical point, which is located atm 2 ≈ 5.5 in the thermodynamic limit. The straight line is , which is the expected value in the thermodynamic limit.
because z corresponds to a logarithmic length scale. Although computational cost increases with increasing z * , it is still much cheaper to use the on-site boundary condition with a logarithmically larger z * than using the full boundary condition whose computational cost increases much faster. In appendix D, we show that the solutions obtained from the onsite boundary condition are indeed identical to the ones obtained with the full boundary condition. From now on, we will use the on-site boundary condition which allows us to obtain full numerical solutions for larger lattices. The plots of the hopping fields for larger lattices can be found in appendix D.

Analysis
In this section, we present the main results of the paper. We first establish the presence of the critical point that divides the insulating phase and the superfluid phase based on the correlation functions measured at the UV boundary. Then we examine how different phases exhibit distinct geometric features in the bulk.

UV boundary
The correlation functions, 1 N (φ * i · φ j ) are given by t * ij (z = 0), which are measured at the UV boundary. In the thermodynamic limit, lim r→∞ t * i,i+r (0) is the order parameter for the phase transition between the insulating (symmetric) phase and the superfluid (symmetry broken) phase. Although there is no real phase transition in finite lattices, there exists a rather sharp crossover which becomes the phase transition in the thermodynamic limit. To identify the 'would-be' phase transition in finite systems, we plot t * ij (0) for the largest possible |i−j| as a function ofm 2 for L = 9 in figure 6. Although it is not strictly zero even in the insulating phase due to finite size effect, there is a sharp crossover aroundm 2 ≈ 5.5. The critical point in the thermodynamic limit is atm 2 c = 5.49454 (see appendix E). From now on, we will use the notation to denote the deviation of the mass away from the thermodynamic critical point.
In figure 7 we plot the full correlation functions t * ij (z = 0) as functions of |i − j| for various values of ∆m 2 . In the insulating phase, the correlation functions decay exponentially as expected. The correlation length tends to diverge as the critical point is approached, as is shown in figure 7(a). In the superfluid phase, the correlation functions approach non-zero values in the large distance limit, exhibiting a long range order. As the critical point is approached from the superfluid side, the long range order disappears continuously as is shown in figure 7(b).
At the critical point, we expect that the correlation function decays in a power-law, with ∆ = 1/2 for D = 3. In figure 8(a), the correlation function is shown at ∆m 2 = 0. The curvature in the log-log plot indicates that the correlation function actually decays slower than an algebraic decay. This is due to the finite size effect which overestimates the correlations due to the periodic boundary condition. The finite size effect effectively pushes the system at ∆m 2 = 0 to the superfluid side at finite L. As will be shown in section 4.2.3, the correlation function decays algebraically at a positive value of ∆m 2 that goes to zero only in the infinite L limit. To take this into account, we include a constant piece and fit the correlation function with as is shown in figure 8(b). In the thermodynamic limit, the exponent ∆(L) approaches 0.55 and G(L) tends to vanish as is shown in figures 8(c) and (d). We expect that the 10% error in the critical exponent can be made smaller by increasing the number of basis points JHEP08(2015)107 that are included to parameterize the radial profile of t ij (z) and t * ij (z) in the numerical calculation. For details, see appendix C.

Bulk
Now we set out to understand the different phases from the behaviors of the bulk fields in the (D + 1)-dimensional space.

Insulating phase
In figure 9, the scale dependent hopping fields are plotted in a deep insulating state at ∆m 2 = 6 for L = 13. Although all further neighbor hoppings are generated in the bulk, they remain exponentially small in |i − j| at all z. Furthermore, all hopping fields decay exponentially in z in the IR region. Due to the exponential decay both in |i − j| and z, t ij (z) is significant only for small |i − j| and z. As the connectivity of space becomes weaker with increasing z, the correlation functions t * ij (z) become smaller as is shown in figure 10. As z increases, the overall magnitude of t * ij (z) decreases exponentially while the rate of decay in |i − j| is unchanged. In the large z limit, the space completely loses its connectivity, and gets fragmented into isolated islands. In the RG language, the IR fixed point is described by the decoupled sites. One can quantify the scale of fragmentation in terms of the susceptibility which is proportional to ξ 2 in the insulating phase, where ξ is the correlation length. As z increases, the susceptibility decreases because the hopping amplitudes become smaller in the bulk. In the large z limit, χ saturates to the value determined by the on-site mass term. This is shown in figure 11. Although fragmentation is not a sharp transition, we can choose a convenient criterion for definiteness: we define the fragmentation scale z F to be the scale   Figure 11. Logarithmic plot of the susceptibility χ as a function of z in the insulating phase with ∆m 2 = 1.5 (a) and 0.5 (b) for L = 13. Fragmentation scale is the crossover scale beyond which χ(z) becomes independent of z and becomes a constant whose value is determined by the on-site action S 0 without hopping. For concreteness, we define z F to be the scale at which the second derivative of the logarithm of χ(z) is maximum.
at which d 2 lnχ(z) dz 2 is maximal. It represents a scale around which a crossover occurs from a connected lattice to fragmented space. In figure 11, z F is indicated in the insulating phase, which shows that z F increases with decreasingm 2 , as expected. Fragmentation occurs at longer distance scales when the system is closer to the critical point. The background color in figure 12 represents the value of the susceptibility in the bulk for different values of ∆m 2 . The curve denoted by z F represents the fragmentation scale. The curve z H and other features in the superfluid phase (∆m 2 < 0) of figure 12 will be discussed in the following section.
Near the critical point, the fragmentation scale z F can be used to extract the critical exponent ν, which dictates how the correlation length depends on the distance to the critical point, ξ ∼ (∆m 2 ) −ν . The susceptibility is given by χ ∼ ξ 2 at the UV boundary. Inside the JHEP08(2015)107  Figure 13. The correlation length ξ plotted as a function of ∆m 2 on a log-log scale. The data is obtained from solving the UV boundary theory with lattice size L = 151. In the region where 1 ≪ ξ ≪ L/2, the scaling is close to the one expected in the thermodynamic limit, ξ ∼ (∆m 2 ) −1 , which is shown as a straight line.
bulk, it falls off exponentially in z as χ(z) ∼ ξ 2 e −2z because of eq. (3.5) before saturating to a fixed value of 1 m 2 . Therefore, the fragmentation scale is given by z F ∼ ln ξ ∼ −ν ln(∆m 2 ). In figure 13 we show ξ as a function of (∆m 2 ), which shows ν ≈ 1, the expected value for the interacting model we are studying in the large N limit. Therefore the fragmentation scale should behave as z F ∼ − ln(∆m 2 ) in the thermodynamic limit.

Superfluid phase
In figure 14, we display the hopping fields inside the superfluid phase at ∆m 2 = −3 for L = 13. Near the UV boundary, t ij (z) decays exponentially in |i − j| as is the case in the insulating phase. This has to be true even in the superfluid phase because only nearest neighbor hoppings are turned on at the UV boundary, and further neighbor hoppings are generated gradually in the bulk. The exponential decay of the hopping fields is well captured by with ψ(z) > 0 near the UV boundary as is shown in figure 15(a). However, the rate of the exponential decay, ψ(z) becomes smaller as z increases, and vanishes at a critical scale z H . This is shown in figure 15(b). For z ≥ z H , the hopping fields no longer decay exponentially but decay algebraically with a constant off-set, where the constant off-set turns on continuously as with an exponent β H ≈ 1 for z ≥ z H . After the off-set is subtracted out the hopping fields exhibit pure power-law decays as is shown in figure 16(a). The z dependence of the off-set near z H is displayed in figure 16(b). The scale z H represents a 'horizon' at and beyond which the system loses locality as sites are connected with each other through non-local hoppings. As a result, the 'coordinate distance' |i − j| is no longer a good measure of the 'physical distance' that dictates the actual correlations. Because of the non-local hopping, the physical distance between two sites effectively shrinks to zero. If we were studying a system with a finite temperature JHEP08(2015)107 with a compact Euclidean time direction in such a way that its length β would remain finite as L → ∞, we would first observe this loss of locality only in the Euclidean time direction. In that case, z H corresponds to the usual horizon at which only the thermal circle shrinks to zero. In the present case, the L → ∞ limit corresponds to the zero temperature limit, where all directions are equivalent. As a result, the size of space shrinks to zero in all directions at the horizon. It will be interesting to understand what this generalized horizon corresponds to in Minkowski space.
If the off-set C was indeed nonzero in the thermodynamic limit, it would imply that any two sites remain coupled with a non-vanishing hopping no matter how far they are. However, figure 17(a) shows that C systematically decreases as the system size increases. In order to understand the behavior of C in the thermodynamic limit, we need to do a JHEP08(2015)107  finite size scaling analysis. Figure 17(b) shows that the maximum value of C at fixed ∆m 2 decays as C max ∼ 1 L r with r = 2.81 ± 0.12 for −3 ≤ ∆m 2 ≤ −0.1 as the system size L increases. This implies that C goes to zero, and the hopping fields decay in a purely algebraic manner for z ≥ z H in the thermodynamic limit.
Although the off-set C vanishes in the thermodynamic limit, the location of the horizon z H remains finite, as is shown in figure 18(a). In the deep superfluid phase, z H is more or less independent of ∆m 2 , but z H sharply increases as the critical point is approached from the superfluid side. This is displayed in figure 12. For L = 13, the horizon scale diverges at ∆m 2 = 0.2, and the horizon no longer exists for ∆m 2   of the horizon in the bulk, we define ∆m 2 c (L) to be the critical mass for finite lattices, e.g., ∆m 2 c (13) = 0.2. Although the critical point is at ∆m 2 = 0 in the thermodynamic limit, ∆m 2 c (L) is different from zero for finite L due to the finite size effect. A finite size scaling in figure 18(b) suggests that ∆m 2 c (L) indeed goes to zero in the thermodynamic limit. In the absence of the off-set in the thermodynamic limit, what distinguishes the horizon at z = z H from the region inside the horizon with z > z H is the exponent κ with which the hopping fields decay. Figure 19(a) shows that κ inside the horizon is smaller than κ H defined at the horizon. The exponent κ tends to approach nonzero values for all z in the thermodynamic limit as is shown in figure 19(b).
Remarkably, the exponent κ H at the horizon is independent of ∆m 2 as is shown in figure 20. This suggests that the horizon is characterized by a universal exponent associated with the decay of the hopping fields. The location where the decay of the hopping fields exhibits the universal exponent precisely coincides with the location where C turns on and ψ vanishes in finite lattices. It is interesting to note that the hopping fields exhibit universal 'critical behavior' that is akin to a continuous phase transition. This should not be confused with the usual phase transition where the correlation length diverges as ∆m 2 is tuned from the insulating phase to the superfluid phase. The 'transition' that happens at z H is associated with the divergence of the length scale for the renormalized hopping fields within the superfluid phase. As the 'short-distance' modes are integrated out, the effective action is renormalized by further neighbor hoppings. In the insulating phase, the further neighbor hopping fields remain exponentially small, and the renormalized action remains local at all scales. In the superfluid phase, the system can not keep the locality as further neighbor hoppings are proliferated beyond the critical scale z H . The presence/absence of the horizon can be used as a 'holographic order parameter' that distinguishes the superfluid/insulating phases.
The non-locality in the bulk could have been avoided if one had allowed for spontaneous symmetry breaking by turning on an infinitesimally small symmetry breaking field before taking the thermodynamic limit. In that description, one starts with a new vacuum with broken symmetry, and fluctuations around the new vacuum are described by the Goldstone modes. Here the spontaneous symmetry breaking is not allowed because the large volume limit is taken in the absence of the symmetry breaking field. The emergence of the non-local geometry is a sign that the system cannot maintain both locality and the full symmetry in the superfluid phase.
As z increases further, t ij (z) eventually decays exponentially in z while remaining flat in |i − j|. This leads to fragmentation in the deep IR region even in the superfluid phase. Figure 21 shows that the susceptibility χ becomes negligible beyond a fragmentation scale z F . The fragmentation scale z F for L = 13 is displayed along with the horizon scale z H in figure 12. However, the fragmentation in the superfluid phase is a finite size effect unlike in the insulating phase. This can be seen from figure 22, where z F increases logarithmically in L in the superfluid phase. This is in contrast to the insulating phase where z F is JHEP08(2015)107  Figure 22. L-dependence of z F in the insulating phase and the superfluid phase. The L axis is shown on a logarithmic scale. The non-zero slopes imply that z F diverges logarithmically in L in the superfluid phase. On the contrary, z F remains finite in the thermodynamic limit in the insulating phase. independent of L. In the thermodynamic limit, z F = ∞ in the superfluid phase while z F stays finite in the insulating phase. In this sense, the superfluid phase is characterized by the algebraically non-local geometry in the IR limit.
In figure 23, the correlation functions t * ij (z) are plotted in the bulk in the superfluid phase. t * ij (z) at each z is fitted to eq. (4.3). The constant piece describes the long-range order, and the power-law decay originates from the Goldstone mode in the superfluid phase. The exponent ∆(z) remains constant as a function of z. In the thermodynamic limit, the exponent approaches ∆(L = ∞) ≈ 0.55 ( figure 8(c)). The correlation function in the bulk is insensitive to the divergence in the length scale of the hopping field across the horizon, which is at z H = 2.8 at ∆m 2 = −3.

Critical point
In figure 24(a), we show the hopping fields t ij (z) as a function of |i − j| for L = 13 at the finite size critical point, ∆m 2 = ∆m 2 c (13) = 0.2. The hopping fields are well fit by the form in eq. (4.5). As z increases, the rate of exponential decay ψ(z) decreases as is shown  in figure 24(b). Because of the invariance under the scale transformation at the critical point, ψ(z) is expected go to zero exponentially in the large z limit (it is noted that z is the logarithmic length scale). However, ψ in the IR region deviates from the expected exponential behavior due to finite size effect for small L ( figure 24(b)). To reduce the finite size effect, we push our system size to L = 21 to read off ψ(z) at the critical point. As is shown in figure 25, ψ indeed decays exponentially in the IR region. At the critical point, the hopping fields retain locality with the characteristic length scale 1/ψ. In other words, the sites are more or less globally connected at length scales smaller than 1/ψ, and the theory is local only at larger length scales. However, we don't expect there to be a sense of flat geometry at the length scale of 1/ψ because the present theory is weakly coupled in the large N limit.
As ∆m 2 is tuned across the horizon at a fixed z deep inside the bulk, t ij shows the same critical behaviour as it does for increasing z with fixed ∆m 2 in the superfluid phase, which is discussed in section 4.2.2. This is shown in figures 26(a) and 26(b). On the insulating side, the hopping fields decays exponentially, where the rate of the exponential decay (ψ) continuously vanishes as the critical point is approached. On the superfluid side, the hopping fields decay algebraically with the off-set C from eq. (4.6) that continuously turns on across the horizon for finite size lattices. The dependence of ψ and C on ∆m 2 across the horizon at a fixed z is displayed in figure 27. Once the constant off-set is subtracted, the hopping fields decay algebraically inside the horizon. This is displayed in figure 26(b). While the hopping fields decay with the universal exponent at the horizon, the exponent inside the horizon changes continuously, as is shown in figure 28.  is the case in the superfluid phase, ∆ is independent of z and approaches ∆(∞) ≈ 0.55 in the thermodynamic limit.

Holographic phase diagram
Now we combine all the information to construct a holographic phase diagram in the bulk. In figure 12, the scale z F marks the crossover beyond which the space is fragmented. In the insulating phase, z F is largely independent of L: the fragmentation of space is a generic feature of the insulating phase. On the other hand, z F diverges in the thermodynamic limit for ∆m 2 ≤ 0. This is confirmed from the finite size scaling in figure 22. Therefore, there is no fragmentation outside the insulating phase in the thermodynamic limit. The superfluid phase is distinguished from the insulating phase by the presence of the horizon in the bulk. The horizon is characterized by the power-law decay of the hopping field, t ij (z H ) ∼ 1 |i−j| κ H with the universal exponent that approaches κ H ≈ 1 in the thermodynamic limit. At the critical point, the hopping fields decay as t ij (z) ∼ e −ψ(z)|i−j| |i−j| κ(z) , where ψ(z), κ(z) > 0 and ψ(z) → 0, κ(z) → κ H as z → ∞. Figure 1 summarizes the holographic phase diagram in the thermodynamic limit.

JHEP08(2015)107 5 Summary and discussion
In summary, we derive the holographic action for the U(N ) vector model regularized on a lattice. The bulk equations of motion are solved numerically for finite lattices in three dimensions. From the finite size scaling, we find that the insulating phase, the superfluid phase and the critical point exhibit distinct geometric features in the bulk. The IR geometry of the insulating phase is characterized by ultra-locality with a decoupled lattice. The superfluid phase exhibits a horizon at and beyond which the geometry becomes nonlocal. The critical point shows a local geometry with a characteristic length scale that asymptotically diverges in the IR limit.
Although the U(N ) vector model is exactly solvable in the large N limit, the present holographic description allows one to understand the concrete model from a holographic perspective. This is a first step toward applying QRG to more non-trivial models whose solutions at large N are not known (e.g. matrix models). An advantage of using finite lattices is that the finite size effect can be used in solving the problem using a simplified IR boundary condition numerically. We finish by listing some open problems for future investigation.
Analytic solution. Although the numerical solutions provide a great deal of information on the behaviour of the bulk, it is desirable to have an analytic solution, especially at the critical point. One approach is to reduce the infinite set of equations of motion for the bi-local fields to a finite set by projecting the solutions to the ones constrained by the numerical solution. In general, it will be of interest to better understand theoretical structures of bi-local (or multi-local) field theories.
Finite temperature. The Euclidean theory in D dimensions can be viewed as the (D − 1)dimensional quantum theory at zero temperature in the imaginary time formalism. One can turn on finite temperatures by making the size of the thermal circle finite. In this case, the hopping fields in the temporal direction will be different from those in the spatial directions. At the horizon, the non-locality is expected to develop only in the temporal direction but not in the spatial directions. In other words, the size of the thermal circle will shrink to zero while the spatial area remains nonzero at the horizon.
Application to Fermi surfaces. Having understood the insulator to superfluid phase transition in the bosonic model, one can try to study the fermionic counterpart. In the fermionic system, the superfluid phase is replaced by an itinerant state with a Fermi surface. It will be of great interest to understand how the Fermi surface manifests itself in the bulk geometry [34][35][36][37][38].
Full (D + 1)-dimensional diffeomorphism invariance in the continuum limit. As is discussed in section 2, the holographic action for the lattice model possesses a subset of the full diffeomorphism invariance of the continuum space. At each D-dimensional slice with a constant z, the local permutation symmetry is present. Given that the lattice provides a UV complete theory which flows to the conformal field theory in the continuum (longdistance) limit, it is expected that the holographic theory for the lattice model recovers JHEP08(2015)107 the continuum theory with the full diffeomorphism invariance in the IR region. However, this needs to be understood more explicitly.
Connection to the higher spin theory. In the continuum limit, it may be possible to relate the equations of motion in the bulk with those of the higher spin gauge theories proposed as the holographic dual for the vector models [33,[39][40][41][42][43][44]. However, the connection is not clear a priori because the form of the bulk theory is sensitive to the specific regularization scheme. In particular, we believe that the 1/N corrections cannot be included in the present formalism of the higher spin theory where the quartic interaction is implemented only through the UV boundary condition [33]. This is because the quadratic theory with fluctuating sources is ill-defined unless one keeps the double-trace operator in the bulk, as is discussed in section 2.
Therefore, to derive the bulk action that is continuous in the holographic direction we need to take into account only terms that are linear in 1/μ 2 i . Keeping this in mind we integrate out the high energy fieldsφ,φ * using only their bare action. For this we write Z as is the normalization factor. We need to calculate the correction to the action ∆S (0) ′′ 1 to first order in dz. To make the calculation more tractable we rewrite and rename the terms from the sum The goal is to calculate everyZ p to order dz, then re-exponentiate to get ∆S (0) ′′ 1 to order dz.Z 1 is given bỹ The first term of eq. (A.1) comes from tracing out the high energy modes from the hopping term, while the second term comes from the renormalization of the low energy mass due to integration of the high-low energy 4-boson interaction. The two contributions are illustrated in figure 30.Z 2 is given bỹ −t The three contributions are illustrated in figure 31. The first one represents fusion of hopping links to generate further hopping links (a). The second represents the fusion of two 4-boson vertices to generate a 6-boson vertex (b), while the final two terms refer to the fusion of the 4-vertex with the hopping term (c). There are no more terms in anyZ p that are linear in dz, which implies that ∆S Here we review the solution of the U(N ) vector lattice model at N → ∞ [45]. Our initial partition function at any given scale is We introduce a Hubbard-Stratanovich fieldσ in the following way Redefining We want to replace σ i with its saddle point value at N → ∞, which is assumed to be a site-independent constant σ. To find σ we integrate out the φ and φ * fields after going to momentum space The effective action for σ is now (ignoring constant terms) and the saddle point equation becomes which is the self-consistency condition for the mean field value σ, eq. (3.1).

C Numerical method
Here we discuss in detail the method we use to solve the equations of motion (EOM) of eq. (2.21). The equations themselves are first order and seemingly simple, however each equation for t ij and t * ij involves a sum over all other fields t qp and/or t * qp . The number of EOM grows as the volume V , and the number of terms in the sum also grows as V , leading to a very large number of terms. The equations are also non-linear. We use a tau spectral collocation method [46][47][48] combined with a Newton-Raphson algorithm [49] to solve this system.

JHEP08(2015)107
Each bilocal field is decomposed into a truncated polynomial series where T α x(z) is the α-th Chebyshev polynomial of the first kind. Here x(z) must live on [−1, 1], which is the domain on which the Chebyshev polynomials are orthonormal.
The domain of the RG scale z is [0, z * ], which implies that x(z) = 2z z * − 1. The number n is called the order of the tau spectral approximation, and n + 1 is also the number of collocation points, i.e. the number of points on the x/z-axis at which the spectral Chebyshev representation agrees exactly with the original function. n is also a measure of the accuracy of the solution; for smooth functions, one expects the error in the solution to decrease exponentially with n. In our case we take the collocation points to be the Chebyshev Gauss-Lobatto (CGL) points x α = cos πα n . For each lattice size we take n = 20. We solve our equations for the coefficients a, b instead of the full functions t, t * . This allows us to solve one large (non-differential) matrix equation rather than many differential equations. For this we rewrite the EOM and boundary conditions in terms of the new "vectors" a ij and b ij , where the {ij} label the lattice link and the vector label is the Chebyshev label (labeled as "α" in eq. (C.1)).
The rewritten EOM have 4 types of terms. The first is a constant, which is represented as a vector with that constant as the first component (α = 0) and all other components zero. The second type is a linear term, which is just c ij t ij (or c ij t * ij ) and it becomes c ij a ij (or c ij b ij ). The third type is the differential operator. We denote it by the matrix L and for the CGL points it is given by where c(γ) = 1 + δ γ,0 . The last type of term is the quadratic term, e.g. t ij t pq . We choose to express it in spectral form via a procedure described in section 3.4.4 of [47]. The first step of the procedure is to expand the sum to order m = 2n, leaving all the coefficients zero that are higher order than n, i.e. a ij α = a pq α = 0 for n < α ≤ 2n. Then we perform an inverse discrete Chebyshev transform to obtain the real space coefficients for every field, which are given byã We multiply the real space coefficients to obtain the coefficients of the product in real spacẽ s ij,pq β =ã ij β ·ã pq β , for β ∈ {0, m}. We now perform an inverse discrete Chebyshev transform back to spectral space, which is given by is the spectral representation of the product of t ij (z) t pq (z). Now it is straightforward to express the EOM and boundary conditions in spectral form as a set of algebraic equations. The number of equations is equal to the total number of fields t ij and t * ij times n + 1 (the number of collocation points), which we enumerate and put into vector form R( a ij , b ij ) = 0. Finally, we combine all of the { a ij , b ij } into one big vector c R( c) = 0 . (C.2) As R is quadratic in c, we cannot solve it directly via a matrix inversion. We employ an approximation method, namely a Newton-Raphson method, to solve eq. (C.2). For this we linearize the equations and solve for c 1 . This process is iterated w times to arrive at a solution c w , which satisfies L 2 w ≡ α ( c w − c w−1 ) α < ǫ. We choose ǫ = 10 −11 , which is small enough to indicate that the Newton-Raphson method has converged up to floating point round-off.
An important question is that of convergence: when does the Newton-Raphson method converge, depending on the initial starting vector c 0 , and how many iterations does it take? The convergence domain that we find is summarized below (the reader should keep in mind that in tuningm 2 = m 2 −t ii we changet ii and keep m 2 constant (in this case m 2 = 25)).
The solution in the deep insulating phase is available from the perturbative analytical solution of eq. (3.4). From this starting point one can obtain the solution at a certain distance away from the critical point on the insulating side, say ∆m 2 = 10. From this solution the Newton-Raphson method converges quickly (within 10 iterations) for any point ∆m 2 > 0. Also, it converges at ∆m 2 0, however more slowly (≈ within 20 iterations). This gives us the solutions at and just beyond the critical point on the superfluid side, which are of the most interest physically. Starting from a solution at 10 ≫ ∆m 2 > 0, we see a faster convergence to these desired points. Accessing points in the deep superfluid state however is more difficult. To definitively converge to a point ∆m 2 1 ≪ 0 one has to start with ∆m 2 0 = ∆m 2 1 + δt ii , where δt ii depends on the lattice size L. For L = 13 we find that δt ii ≈ 1 is a good choice. The slower convergence in the superfluid phase is due to the sensitivity of R( c) in eq. (C.2) on c.
In summary, the numerical code can be parallelized when finding solutions on the insulating side, as well as at and near the critical point on the superfuild side, but finding solutions in the deeper superfluid phase needs to be done in serial.  Figure 32. The same plots as in figure 4 with the difference being that the IR BC used is the on-site one of eq. (3.6), instead of the general BC of eq. (3.2). plotted in figures 32 and 33. We can see that the fields look identical to those of figures 4 and 5, and after comparison we find that they indeed lie on top of each other. However, besides the bulk shape it is also important to compare the exact IR values to those obtained using the general BC, especially when the theory starts out in the superfluid phase at the UV. Tables 1 and 2 summarize this comparison for two values ofm 2 = 25 and −15, corresponding to the deep insulating and deep superfluid phases, respectively. In both phases the t * i =j (z * ) values are the only ones that differ at all between the two BCs. For the on-site BC they are zero always, which is what we set them to (10 −15 is the machine precision), and for the general BC this is not the case. However, t * i =j (z * ) for the general BC are still extremely small when compared to t * ii (z * ) in both phases, which is why it works to approximate them by zero. This explains the validity of this approximation.

D Comparison between the general and on-site boundary conditions
In figure 34, we show the hopping fields t ij (z) for larger lattices L = 5, 7, 9 obtained using the on-site BC. They look qualitatively similar to those in figure 32.  Figure 33. The same plots as in figure 5 with the difference being that the IR BC used is the on-site one of eq. (3.6), instead of the general BC of eq. (3.2). m 2 = 25 t ii (z * ) |t i =j (z * )| t * ii (z * ) |t * i =j (z * )| general BC 0.08 ≤ 1.3 × 10 −7 0.04 2 × 10 −10 on-site BC 0.08 ≤ 1.3 × 10 −7 0.04 < 10 −19 Table 1. Comparison at z = z * of general and on-site BCs form 2 = 25 (deep insulating phase).