Towards a Holographic Bose-Hubbard Model

We present a holographic construction of the large-N Bose-Hubbard model. The model is based on Maxwell fields coupled to charged scalar fields on the AdS2 hard wall. We realize the lobe-shaped phase structure of the Bose-Hubbard model and find that the model admits Mott insulator ground states in the limit of large Coulomb repulsion. In the Mott insulator phases, the bosons are localized on each site. At zero hopping we find that the transitions between Mott insulating phases with different fillings correspond to first order level-crossing phase transitions. At finite hopping we find a holographic phase transition between the Mott phase and a non-homogeneous phase. We then analyze the perturbations of fields around both the Mott insulator phase and inhomogeneous phase. We find almost zero modes in the non-homogeneous phase.


Introduction & Summary
Strongly-coupled many-body systems such as the QGP [1] have been studied using the gauge/gravity correspondence [2], which includes useful string theory embeddings of systems which share many qualitative properties of QCD at both zero and finite temperature/density. Following these advances, holographic models of condensed matter phenomena have been constructed, such as an analog of the particle-vortex duality [3,4], superfluid-insulator transitions [5], superconductivity [6,7,8], and the field theory with Lifshitz symmetry [9] as well as hyperscaling violation [10]. However, many of these models assume translation invariance. In the presence of finite charge density, the DC conductivity becomes infinite unlike in real materials which have a Drude peak characterized by a finite DC conductivity. Recently, holographic models have also been applied to systems in which translation invariance is broken by a periodic function of the chemical potential [11] or by that of scalars [12,13,26,27] to form a lattice. Interestingly, physics without translation invariance can also be captured by introducing a lattice of impurities dual to probe D5(D5) branes wrapped on an asymptotic AdS 2 in the AdS 5 or AdS 5 black hole background. In Ref. [19,20], a dimerization transition on a lattice, which changes the structure of the Fermi surface, was analyzed by coupling defect fermions on D5-D5 branes with itinerant fermions. Many interesting features of the lattice formulation can already be captured in the probe limit [14,15]. These simple holographic probe models however do not easily allow charge transport of fermions on the probe brane.
In this paper, we are interested in and focus on an approximate (bottom-up) model of interacting bosons on a lattice, including charge transport: the Bose-Hubbard model. The Bose-Hubbard model is e.g. realized in ultra-cold atomic experiments using 87 Rb trapped in an optical lattice [16], as well as in Helium atoms moving on substrates. The Bose-Hubbard Hamiltonian is given by where t hop is the hopping parameter giving the mobility of the bosons between neighboring sites and U is the remnant of the repulsive Coulomb interaction 1 between bosons on a single site. b j and b † j are, respectively, annihilation and creation operators for the bosons at site j. n i = b † i b i (no summation over i) is the boson density on site i, and µ b is the chemical potential. In the Bose-Hubbard model without disorder, there exist only two phases, namely, the Mott insulator phase and the superfluid phase. In the Mott insulator phase, bosons are localized on the lattice due to the repulsive interactions. They do not form a coherent state. In the coherent superfluid phase, bosons are delocalized on the lattice and an off-diagonal long-range order (ODLRO), i.e. long-range correlations b † i b j , exists in the superfluid phase. It is known that this condensate becomes of the same order as the particle density. Our holographic model will show the same pattern. When U/t hop is large, bosons are localized and the ground state is in the Mott insulator phase, while when U/t hop is small, bosons acquire kinetic energy derived from the non-zero hopping parameter t hop and are delocalized.
In the rest of the introduction, we will first review the salient details of the mean-field approach of [35] to the Bose-Hubbard model, and the physics of the zero-temperature phase diagram ( §1.1). Then in §1.2 we will introduce our holographic model and compare our results to the phase diagram of [35].

Review of mean-field approach to the Bose-Hubbard model
In this section we review the field theoretic approach of [35] to the phase diagram of the Bose-Hubbard model in the µ b − t hop plane at zero temperature. We will then attempt to construct a holographic model which exhibits the same physics. First let's consider the Hamiltonian (1.1) with t hop = 0. This describes the phase structure of a lattice of bosons with no hopping as a function of chemical potential. The ground state will be the one which minimizes the potential energy, at each site. Taking U > 0, the solution is that n bosons occupy each site for µ b in the range n − 1 < µ b /U < n. If µ b < 0, then n i = 0 at all sites. If µ b /U is exactly a positive integer, m, then (m) = (m + 1), and configurations with m and m + 1 bosons at a single site are degenerate, leading to 2 N possible degenerate ground states for a system with N sites. We see that as we increase µ b /U with t hop = 0, the system undergoes a series of level-crossing phase transitions at each integer value of µ b /U where the density of bosons on each site jumps from n i = n to n i = n + 1. Now let's consider what happens when one turns on a small hopping t hop > 0 when in a phase with n bosons per site, with µ b /U = n − 1/2 + α for −1/2 < α < 1/2. The energy required to add a particle to the system is δE P ∼ (1/2 − α)U , and to remove a particle (create a hole) is δE h ∼ (1/2 + α)U . For t hop U , the kinetic energy gained from allowing an extra particle or hole to hop around the lattice is not large enough to overcome  [35], illustrates the zero-temperature phase diagram as a function of chemical potential and hopping parameter. In their notation, J c = t hop , V = U , and N is the number density. The area inside the lobes is a Mott Insulator phase and outside a superfluid phase. We will reproduce many aspects of this phase diagram with our holographic model. the cost in potential energy of removing a particle or adding a particle. Therefore, for each value of the density n, there is some finite region of fixed density with an energy gap for the creation of particle-hole excitations, where the size of the gap decreases with increasing t hop . Furthermore, this constant-density state is a Mott insulating phase as the fixed density implies it is incompressible.
At fixed nonzero t hop , as one increases or decreases µ b , eventually it becomes energetically favorable to add a particle or hole which can hop around the lattice. At this point, (at zero temperature) the particles will instantly Bose condense and the system will undergo a phase transition to a superfluid. The boundary of the Mott insulating-superfluid phase transition will extend to the t hop = 0 axis at precisely the points where µ b /U is an integer, since at these points the occupation numbers n and n + 1 are degenerate, and there is no energy cost to adding extra particles. Thus we see that the phase boundaries form a series of lobe shapes extending from the t hop = 0 axis to some finite value of t hop . This can be seen in Figure 1, taken from [35].
Generically, the transition from the Mott insulating phase to the superfluid phase is continuous and second-order, driven by the addition of a small number of particles and/or holes to the system, and a change in density from the integer value in the Mott insulating phase to a non-integer value slightly larger or smaller. However, at the tip of the lobes, the Mott insulating phase and the superfluid phase both have integer density n. In this case, the transition is driven by an increase in t hop which allows the bosons to overcome the on-site repulsion and delocalize despite the fact that the density remains constant. As the density µ b /U increases, the amount of kinetic energy required to overcome the potential barrier to create one particle/hole pair decreases, and thus the location of the tips of the lobes decrease. This transition is also continuous and second-order, but in a different universality class than the density-driven transition which occurs at all points along the phase boundary except the tips of the lobes.
One can analytically reproduce the phase boundaries of Figure 1 via a mean-field analysis as in appendix A of [35]. Taking an infinite-range hopping limit, they derive an action of the form where ψ is an order parameter which distinguishes between the superfluid, ψ = 0, and insulating, ψ = 0, phases, N is the total number of lattice sites, and T is temperature.
, (1.4) and minimization of S ∞ corresponds to the condition r(µ b , t hop ) = 0. This reproduces the phase diagram of Figure 1.

A holographic model
The purpose of this paper is to realize the Bose-Hubbard model via holography. There are many different motivations for studying the Bose-Hubbard model via holography. While exact results are known for the Bose-Hubbard model in low dimensions, the higherdimensional Bose-Hubbard model is difficult to analyze using field theory techniques alone. It is however possible to study the Bose-Hubbard model using numerical simulations. Such numerical simulations give quantitative answers, but often little insight into the underlying mechanisms. In contrast, holographic models are at best toy models and so can not be trusted to give quantitatively correct answers for the real Hubbard model, but often give important insights into non-perturbative mechanisms. Most importantly, implementing the Bose-Hubbard model holographically is a first step towards constructing the dual of the Fermi-Hubbard model. Unlike its bosonic cousin, the Fermi-Hubbard model is not easily accessible with numerical techniques due to the notorious fermion sign problem. It is however of paramount theoretical and practical importance. It is the paradigm for a model of interacting fermions on the lattice and is believed to give a good representation of important materials such as, for example, high T c superconductors. On the holographic side we do not see the fundamental charge carriers but only the neutral composite operators such as the charge density and b a † i b ja with i = j. These composite operators are bosonic irrespective of whether the underlying charge carriers created and annihilated by b a † i and b ia respectively are bosonic or fermionic. In the bulk, the difference between the Bose-Hubbard model and the Fermi-Hubbard model (or more general models which have both bosons and fermions) should all be in the details of the holographic construction, such as e.g. interactions and boundary conditions. By focusing first on the Bose-Hubbard model, where we can use the well known phase structure to tune our bulk model, we can establish that holographic techniques indeed work. An obvious next step is to see whether one can generalize our construction to the Fermi-Hubbard model, where holographic techniques have the potential to add new insights.
In our construction we consider a gauge theory living on an AdS 2 background with a hard-wall cutoff. We construct a lattice using a quiver-like gauge theory-we introduce a U (1) Maxwell field corresponding to each site in our lattice, and we couple bifundamental scalar fields charged under U (1) i × U (1) i+1 to the gauge fields at adjacent sites. The hopping parameter will be identified with the source of the operator dual to the bifundamental charged scalar. Charged scalar fields in the bulk can carry the electric (or baryonic) charge and are needed to show the transport of defect (large N ) bosons between the sites. The hard wall at r = r h is needed to realize the Mott insulating phase with gapped excitations. In the presence of the infrared cutoff (i.e. the hard wall), we can safely consider the effective theory for AdS 2 in the probe limit [17,18]. This holographic model should be understood to be dual to a Hubbard-like model with SU (N ) fundamental bosons localized on each site instead of the spin zero bosons of the original Hubbard model. For this system to have a good holographic dual we are implicitly working in the large-N limit. This large-N limit also helps us theoretically analyze the phase transition in finite volume. The dictionary is given in Table 1, where spin indices are given by a = 1 . . . N .
Our main result is the derivation of a lobe-shaped quantum phase transition structure in the chemical potential -hopping plane, separating a Mott insulating phase from an inhomogeneous phase, similar to the phase diagram shown in Figure 1. We determine the phase structure by calculating the free energy from the holographically renormalized on-shell action for the phases present, and determine the phase which minimizes the free energy. At zero hopping our model exhibits a level-crossing quantum phase transition which changes the occupation numbers of the bosons between different insulator phases upon variation of the chemical potential.
The order parameter for the insulator-inhomogeneous phase transition is the vacuum expectation value of the operator dual to the bifundamental scalar, which breaks the diagonal subgroup of the U (1) i × U (1) i+1 global symmetry at neighboring sites. As discussed in §1.1, in the Bose-Hubbard model this quantum phase transition is of second order. However, we generically find a first order transition between the superfluid and insulating phases, except at zero hopping (at the cusps between the lobes), where we find a second order transition in the t hop direction. Furthermore, at the cusp points the transitions at zero hopping are first order in the direction of the chemical potential but second order only when we switch on infinitesimal hopping, which is in accordance with the mean-field analysis of [35]. One important result is the identification of the ground state of the Mott insulating phase which dominates when U/t hop is large.
The content of this paper is as follows. In section 2, we consider the case of zero hopping (t hop = 0). We derive the homogeneous ground state (i.e. same boson density at every site) of the Mott insulator from our holographic model. We furthermore show that the transitions between Mott insulator states with different boson densities are first-order level-crossing phase transitions. In section 3, we consider finite hopping. We holographically derive the Mott insulator/inhomogeneous phase transition in a model with only two lattice sites, showing the lobe-like phase structure. A non-homogeneous state in which the charge densities are not equal on both sites is able to dominate over the homogeneous state for t hop = 0. In section 4, we analyze perturbations around both Mott insulator phase and inhomogeneous phase. We find an almost zero boson in the inhomogeneous phase and a gap in the Mott insulating phase. In section 5 we discuss the generalization to the n-site model. Finally, in section 6 we discuss a string embedding similar to the Hubbard model and several future directions.

Zero Hopping
In this section we first analyze our holographic Bose-Hubbard model at zero hopping. We will reproduce the level-crossing transitions reviewed in §1.1. The theory is based on a quiver-like lattice of decoupled U (1) gauge theories on a AdS 2 hard wall, with one theory at each site k of the spatial lattice. The action is simply given by The metric of the AdS 2 hard wall is given by where we have set the AdS radius to be 1, and placed the hard wall at r = r h . This means the radial coordinate is restricted to r ≥ r h . Note that for the above metric √ −g = 1. This theory is a simplified bottom-up version of the D5 brane lattice of [19], with the internal S 4 coordinates as well as the induced brane geometry of e.g. an AdS-Soliton background [61] being neglected. In the following we choose the radial gauge at every lattice site. In the background (2.6) the Maxwell equations are given by The Maxwell equations are solved by Here, the same chemical potential µ is chosen at every site, as required by thermal equilibrium. The coefficient ρ (k) is the charge on that site and ρ (k) = δS kin /δA (k)t . As worked out in e.g. [48], only Neumann boundary conditions are possible for a gauge field in AdS 2 , which fix the charge density ρ (k) and hence force us to work in the canonical ensemble. 2 We impose a Dirichlet boundary condition The unit charge is defined by introducing the fundamental charged object as the bulk source term . The quantization of the charge then depends on the factor of the kinetic term. For our model, ρ (k) ∈ Z.
which is the generalized version of the Dirichlet boundary condition A t | r=r h = 0, at the hard wall.
Since the on-shell action (2.5) has power-law divergences in r at the boundary r → ∞, we follow [21,22,23] and regularize the action by adding the following cutoff term at large r = R (c.f. e.g. the discussion in sec. 3 of [48]) where √ −h = R is the induced metric at the boundary. A couple of comments on gauge invariance are in order: Note that this counterterm is not manifestly gauge invariant. The boundary term is only invariant under gauge transformations which vanish after the integration over the time direction. As usual in AdS/CFT, we require that suitably quickly decaying gauge transformations do not change the leading coefficient of the boundary expansion of the gauge fields. Only these gauge transformations are truly a redundancy of the theory. Large gauge transformations, that is those that do not vanish at the boundary, change the coupling constants of the theory and so do change the physics.
What is crucial for our analysis is that the charge on each site is quantized. That is we assume that all charge carriers in the theory have charges which are integer multiples of some basic charge and hence ρ (k) ∈ Z. This is clearly true in the standard Bose-Hubbard model: the only charge carriers are the bosons and we can take their charge as the basic quantum of charge. If we assume that our theory can be consistently coupled to monopoles, quantization of both the monopole charge and the electric charge directly follow from the Dirac quantization condition. Any theory that arises as a low energy effective theory for a system made of neutrons, protons and electrons needs to be consistent when coupled to monopoles, as the underlying microscopic theory, QCD and QED, can be consistently coupled to monopoles. Note that this is a purely theoretic statement and true irrespective of whether monopoles actually exist in nature. This argument however only guarantees quantization of the net charge; what we here assume is the stronger statement that charges are quantized on each site. Again, this is also a condition imposed by fiat in the Hubbard model.
Another, possibly more intuitive argument for charge quantization at each defect can be made from the point of view of probe-brane top-down constructions using the AdS soliton geometry [61]. For more information on such a construction c.f. sec. 6. Let us e.g. introduce a probe D5-brane on the AdS 5 times S 5 soliton geometry [61] without considering its back-reaction. The D5-brane wraps two uncompactified directions and an S 4 ⊂ S 5 . In addition these branes carry a quantized F1 flux. Such D5-brane embeddings had first been introduced as the dual to the anti-symmetric Wilson loop in d = 4 N = 4 SYM [39,40,41] and are also the main ingredient in the holographic probe brane lattice constructions of refs. [19,20]. As shown in [38], in response to the quantized F1 flux the transverse embedding function of the D5 brane inside the S 5 (the azimuthal angle) takes quantized constant values, which leads to a quantization of the energy (tension) the D5 brane can have in different embeddings. The electric components of the worldvolume gauge field flux will then be quantized as well [38], leading directly to charge quantization on the branes, and to the idea that the D-brane is actually a bound state of a (quantized) number of F1 strings [38]. This ties in nicely with the intuitive picture that the charge on the D-brane corresponds to the number of F1 strings ending on it, which must be quantized. In our bottom-up model we neglect the internal directions crucial for the above argument, and hence we have to impose the charge quantization condition by hand. Note however that the correct limit in top-down models would be the opposite one, in the following sense: The parameter which quantizes the azimuthal angle of the probe D-brane embedding is [38] the ratio n/N , where n is the quantized charge. At fixed n in the large N limit the probe D brane hence wraps an infinitesimally small internal sphere, and the DBI description would break down and higher derivative corrections become important. The correct limit to work in for top-down models would hence be the limit of large n and N , with their ratio fixed.
Including the counterterms, the holographically renormalized action then becomes S kin + S cut . To compute the free energy, we analytically continue [7,8,25] to Euclidean signature. Taking into account the usual minus sign, the free energy is then given by where β = 1/T is the radius of the time-circle. The plus sign in front of the ρ 2 (k) term is important, as it indicates the repulsive Coulomb interactions between bosons at the same site. Conversely, if it were negative, the Coulomb interaction would be negative, and the system could become unstable.
We want to compare the above free energy (2.12) with that of the Bose-Hubbard model at zero hopping [35]. Its free energy is given by The ρ 2 (k) terms describe the on-site repulsive interactions between the bosons. Matching the parameters (µ b , U ) with our parameters (µ, r h ), we find (2.14) We conclude that for zero hopping the free energy of our holographic model agrees with that of the Bose-Hubbard model with a repulsive on-site Coulomb force. Moreover, the number operator ρ i commutes with the Hamiltonian at zero hopping parameter. Thus, the particle number eigenstates are simultaneously diagonalized in the phase with zero hopping, in accord with charge quantization at each site k. In the Bose-Hubbard model, the homogeneous phase with ρ (k) = ρ (or A (k) t = A t ) at all sites is a Mott insulator for an integer occupation state. We expect the same to be true here, and indeed will show this below in sec. 2.1 by directly evaluating the free energy of the different phases, and also in section 4 by analyzing the fluctuations around this state. We will in particular find that the gap is given by the Coulomb parameter U . Hence our state is a Mott insulator where the fundamental bosons localize because of the on-site Coulomb energy U , set by the gap scale r h of the AdS soliton [61]. It is a non-compressible state with gapped excitations and vanishing quasi-particle density of states at energies below the gap scale. In fact, we find that the whole excitation spectrum in the Mott phase is discrete and gapped, as expected from a system in a AdS hard wall geometry.

Level-crossing First Order Phase Transitions
In the following we analyze the phase structure of our system at zero hopping. We focus in particular on the two-site model, and comment shortly on the (at zero hopping completely analogous) multi-site case at the end of this section. Note that in order for more fundamental bosons to occupy each site we should, according to (2.12), consider negative chemical potential. In this case of negative chemical potential, we can compare our analysis with the phase structure of the Bose-Hubbard model at strongly interacting region U/t hop 1. In Fig. 2, we plot the free energy as the function of µ b /U , where for convenience we chose U = r h = 40 in all our numerical calculations in the rest of this paper. Although this choice seems to be arbitrary at first sight due to the rescaling invariance of of AdS 2 , it is important not to choose the cutoff too small in order to avoid the appearance of other possible instabilities at energy scales above U . Since the occupation number is fixed to be an integer, there is a first order transition called the level-crossing phase transition between two Mott insulator phases at the critical values Decreasing µ induces a level-crossing transition where the occupation number at each Figure 2: The free energy of the two-site model at zero hopping as the function of the chemical potential µ. At µ = 0 the lines plotted are, from bottom to top, (ρ (1) , ρ (2) ) = (0, 0), (0, 1), (1, 1), (2, 0), (1,2), and (2, 2). The free energy for (ρ (1) , ρ (2) ) = (0, 0) is zero. Note that the free energy of the non-homogeneous phase (ρ (1) = ρ (2) ) is symmetric under ρ (1) ↔ ρ (2) , i.e. these states are degenerate.
site increases by 1. 3 Note that, in the figure, the free energies of the systems that are unequally occupied such as (ρ (1) , ρ (2) ) = (0, 1) and (1,2) have degenerate energies at the critical chemical potential (2.15). At the phase transition point, the free energy at each site is the same for two different occupation numbers and as the sites are all independent each can individually chose between the two options. Unequally occupied states are never preferred at zero hopping; they always have a larger free energy than the preferred state away from the phase transition points. As seen in the parameter matching, this phase structure is the same as the phase structure of the Bose-Hubbard model in the strongly interacting region. We see that the ground state is described by a Mott insulator. Without hopping, we can easily generalize the level-crossing phase transition to the case of many sites. We find a phase structure similar to the two-site model. At the critical chemical potential, the number of bosons in the degenerate state is different by 1 and the whole system has a 2 M degeneracy where M is the number of the lattice site. The large degeneracy can be interpreted as a macroscopic entropy.

Finite Hopping
In the holographic theory we add bi-fundamental scalars to describe the hopping parameter on the field theory side. As indicated in Table 1, this bi-fundamental scalar is dual to the operator b † i b j in the Bose-Hubbard model, which is exactly the term we add to the Hamiltonian when including hopping. We then analyze the motion of fundamental bosons of the Bose-Hubbard model between sites as a function of the hopping parameter. We mostly consider a simple two-site model (k = 1, 2) and only briefly comment, in the end, on the generalization to a multi-site theory. The bi-fundamental scalar φ is charged under U (1) 2 as (q, −q). 4 The action of the bi-fundamental scalar is given by µ . The last term is the IR potential describing the interaction of the bi-fundamentals. w is a constant which gives an IR mass for the fields; such a term has been analyzed in [30,31]. In our case, a proper choice of IR potential will be necessary to reproduce the quantum phase transition structure of the Bose-Hubbard model at small hopping.
Since we are interested in a static configuration, we consider fields depending only on the AdS radial direction r, namely, A (i) t (r), φ(r). The equations of motion (EOMs) following from the total action are then given by where l = 1, 2. Henceforth we choose the charge q to be There are several reasons for this choice: First, we should choose q 2 large enough to ensure as few subleading corrections to the t hop term in (3.24) as possible before the VEV term sets in, making the extraction of the VEV numerically simpler. Secondly, for the lowest possible value δρ = 1 in (3.24), the particular choice of (3.18) leads to a rational dimension of the operator dual to φ, Finally, a large q 2 also ensures that the probe limit for the gauge field and the bifundamental is valid even for small values of r h . When the charge q is much larger than the gravitational coupling constant, we can ignore the back-reaction onto the metric [32,33]. We assume this to be the case and completely neglect the gravitational sector of the bulk theory in this work and treat the background metric as fixed. Note that the diagonal gauge field

Homogeneous Mott Insulator
When the homogeneous phase A is considered, the EOMs of the fields φ and A (l) t become independent. We can then solve the EOMa of the fields (3.17) analytically. The solutions become where l = 1, 2. We identify the coefficient t hop and the coefficient of the normalizable mode ϕ 0 with the hopping parameter and vacuum expectation value (VEV) of the bi-local field b † i b j in Bose-Hubbard model, respectively.
We choose the Dirichlet boundary condition φ| r=r h = K , (3.21) similarly to the gauge field. We set to zero the VEV of the bi-fundamental by tuning the parameter K such that while the hopping parameter t hop is used as a free parameter.
In the homogeneous phase, we do not need to add counter-terms for the action S matter . The free energy becomes Note that at zero hopping the last term is simply a constant shift of the free energy, and hence does not affect the identification (2.14). At finite hopping this additional term becomes crucial for the physics of our model, in particular, for the existence of the cusps in the phase diagram fig. 3.

Non-homogeneous Mixed State
When we consider the non-homogeneous case A t , an analytic solution to the equations of motion (3.17) does not exist in general, so we solve (3.17) numerically. Although the axial U (1) is explicitly broken by the hopping parameter, and hence the system does not admit a Goldstone mode, the non-homogeneous phase still is quite similar to the superfluid phase of the Bose-Hubbard model in the following sense: As we will show in sec. 4.2, the phase of the VEV of the kinetic energy operator dual to the bifundamental φ can be related to the superfluid quantum current flowing between different lattice sites in the Bose-Hubbard model in a straightforward way. Although in the body of this work we are going to choose a boundary condition such that the quantum current vanishes, it is naturally present for more generic boundary conditions such as the ones discussed in App. A. We hence would like to think of the non-homogeneous phase as a superfluid-like phase, similar to the actual superfluid phase in the Bose-Hubbard model. Note also that the finite hopping parameter forces the system to live in a state with unequal charge densities, hence forcing additional bifundamental fields to condense spontaneously. If order parameters transforming as fundamentals under one of the gauge groups are coupled to this system in the right way, similar symmetry breaking patterns may arise as well. We are going to comment on both possibilities further below. Note also that although below we are going to work with a charge non-homogeneity of δρ = 1, at larger t hop , higher and higher δρ will presumably be the dominating phase. Also, if we had chosen boundary conditions which allow for a bifundamental VEV, such as the ones in App. A, the VEV would grow asymptotically large with larger t hop .
The solutions of the EOMs (3.17) satisfy the following UV asymptotics: where δρ = ρ (1) − ρ (2) and α t = (−1 + 1 − 4q 2 δρ 2 )/2. In the above asymptotic expansion, the subleading corrections due to the hopping parameter are included since this term becomes important in both the asymptotic expansion of the gauge fields and of the bifundamental. For example, there are finitely many correction terms to the hopping term in the expansion of φ before the VEV term ∼ ϕ v takes over. The same is true in the gauge field expansion. 5 Note that to stay above the BF bound that is to keep α t real, we can only consider the case |δρ| = 1 for the integer occupations ρ (l) and our choice of charge (3.18). The bi-fundamental scalar was introduced as the holographic dual to the bi-local field b † i b j . In particular, the identifications of t hop and ϕ v as the hopping term and vev of this operator should still hold in the non-homogenous phase. We can interpret α t as quantum corrections to the dimensions of this operator (i.e. its anomalous dimensions) from the interactions in the non-homogeneous phase.
We choose the boundary condition in which the subleading term ϕ v in the AdS boundary expansion vanishes. For real t hop , as we choose throughout, the imaginary part of ϕ v is related to the quantum current of our theory as discussed in section 4. Our choice of boundary conditions hence ensures a vanishing quantum current in the ground state. The requirement that the real part of ϕ v also vanishes is however a stronger constraint. The reason for choosing this boundary condition can be motivated as follows. Though our system does not exhibit purely spontaneous breaking of U (1) 1 − U (1) 2 , the symmetry is always explicitly broken by the presence of the source, t hop , and the VEV generated by it (c.f. fig. 5 and eq. (3.28)). This VEV turns out to be non-zero even when ϕ v vanishes. We take the ground state to be the one with vanishing current, so we can safely set the 5 Note that in this paper we use the notation ϕ v for the normalizable piece in the inhomogeneous phase, ϕ 0 for the normalizable piece in the Mott phase, andφ = (1 − 2∆ φ )ϕ v for the actual value of the VEV associated with the normalizable mode in the inhomogeneous phase (c.f. App. A).
imaginary part of ϕ v to zero. We want to interpret setting all of ϕ v to zero as demanding that the U (1) 1 − U (1) 2 breaking VEV is entirely forced upon us by the source and has no spontaneous component, which would correspond to the real part of ϕ v . Of course there is no rigorous way to break up the VEV this way and so this argument can at best serve as a heuristic motivation for our choice. Note that in this way we specify two UV boundary conditions in solving the model, and no boundary conditions at all at the IR hard wall. This is different from how one usually proceeds in AdS/CFT. We checked that our fields behave regularly at the hard wall, so in principle we can always reinterpret this UV boundary condition as a particular IR boundary condition. Whatever value the field takes on the IR wall could be viewed as an input that the ensures vanishing of ϕ v . In the case of IR Neumann boundary conditions, which we analyze in App. A, one can view vanishing of ϕ v as a particular choice of potential on the hard wall. 6 It should probably be noted at this point that by imposing charge quantization with order one charges on the gauge field, we essentially chose two UV boundary conditions in this case as well. We specify the chemical potential (the source) but then only allow the charge to take a few discrete values. Both of those are UV data. In a top-down AdS/CFT setup classical equations in the bulk are only valid when the quantized charge is large (or order N ). In this case the charge density can be treated as a continuous parameter in the UV, which needs to be fixed (as usual) by an IR boundary condition. By insisting on order one charges for the gauge field we break the standard rules of AdS/CFT. This could potentially be justified by working at large but finite N . Since we are forced into this situation for the gauge field already, we should probably not be surprised that we need to follow a similar strategy for the scalar as well. Without imposing vanishing of ϕ v we were unable to produce the nice lobe structure we present in here. Since the on-shell action is divergent for the asymptotics (3.24), we add the following counter-terms Then, the free energy is evaluated by the following sum:  Note that we use the chemical potential µ b defined in (2.14). Inside the lobes, the charge density on both sites is equal, and by analyzing the spectrum of excitations (c.f. sec. 4) we identify this homogeneous phase with the Mott insulating phase of [35]. As t hop is increased, there are regions (which we call inhomogeneous phases) where the non-homogeneous states are favored. Note that the width of the lobe is fixed in units of U , basically by the free energy at zero hopping. According to our experience, changing the w parameter in the IR potential (3.16), the height of the peaks behaves as t hop ∼ 1/w. The cusps appear due to the degeneracy of the ground state between homogeneous phase and non-homogeneous phase at the special points on the µ-axis.  Note that in our model at the particle-hole symmetric points µ b /U = −(n + 1/2), the three phases (n, n), (n, n − 1) and (n + 1, n), compete, and the transition is first order even in the non-homogeneous phase (c.f. fig. 4).
Because the diagonal gauge field A V = A (1) + A (2) decouples from the remaining fields, we can rewrite F as This implies that the energy, which is the second term in the above equation, does not depend on the chemical potential at zero temperature. The first few orders of the strong coupling expansion in terms of t hop and 1/U for the energy is given by where we used |δρ| = 1 and dots include a constant shift and higher order terms. The second term is the leading order φ 2 contribution to the energy (see also [42]). In Fig. 3 we plot the phase structure of the two-site model for r h = 40, w 2 = 1, and Λ = 1. The chemical potential µ b is defined in (2.14). We see that our model reproduces the lobe-like structure of the Bose-Hubbard model. We have regions where the inhomogeneous state is favored at finite t hop . In Fig. 3, the inhomogeneous state extends to the µ-axis at µ b /U ≡ µ/U + 1/2 = 0, −1, −2. The width of the lobe is fixed in units of U . Note that the cusps appear wherever the homogeneous and inhomogeneous states have degenerate free energy.
Our experience shows that (at least for the parameter ranges we explored) we can change the amplitudes of all lobes by changing this parameter w. In particular, we find that large w decreases the amplitude as t hop ∼ 1/w. In Fig. 3, the amplitudes of the lobes do not change as we decrease µ, while in the actual Bose-Hubbard model [35], the amplitudes decrease as 1/ρ (1) , due to the fact that as the number density of the sites increases, one needs less kinetic energy to overcome the potential barrier of removing a particle from one site. We will expand on how to reproduce this feature of the lobe structure in the next subsection, sec. 3.3.
In Fig. 4, we plot the free energy as a function of µ b for the fixed hopping parameter t hop = 0.4. There, the green lines are the thermodynamically prefered phases, while the red lines are the nonprefered phases with higher free energy, and the solid line corresponds to the Mott phase, while the dashed line corresponds to the inhomogeneous phase. Note that fig. 4 is completely consistent with the lobe structure in Fig. 3, that is the values of µ b where the lines of free energy cross precisely correspond to the boundaries of the lobes in Fig 3. In Fig. 5 we plot F − F M ott as a function of t hop . As seen in Fig. 5, the Mott insulating phase is dominating the thermal ensemble in the small t hop regime, while the non-homogeneous phase is dominating in the large t hop regime, as expected. Finally, in Fig. 5, we see the existence of a second order phase transition near t hop = 0, reproducing the behavior of the Mott/Superfluid transition in the Bose-Hubbard model at the cusps at small hopping. 7 However, our model shows a first order phase transition between the Mott phase and the inhomogeneous phase for any t hop except t hop = 0. First order transitions are more common in holographic large N theories, where they often arise from the free energy competition of several saddle points, so the appearance of a first order transition in this model should not be surprising. For comments on how to achieve a continuous phase transition in this model, c.f. Sec. 6.

Decreasing the Amplitude of the Lobe
In the previous sections, we used the IR potential (3.16) without a gauge potential. To realize the decreasing amplitudes of [35] and shown in Fig. 1, we need to deform our holographic model. In this section, by coupling the gauge fields with the IR potential (3.16), we reproduce the 1/ρ (1) behavior of the height of the lobe [35] via holography.
The IR potential can be generalized to the following gauge-invariant IR potential: 8 where F νµ n ν (n µ is the boundary normal satisfying n µ n µ = 1) and F rt . To connect with the level changing transitions in the t hop = 0 case, we require the IR potential to be a purely additive constant to the free energy at zero hopping. In particular it should not explicitly depend on the chemical potential µ b . It can however depend on the charge density in the diagonal sector, i.e. through the quantized field strength F V = i ρ (i) , which does not depend on the radial direction. We plot the phase structure of the model with the above potential in Fig. 6 by setting Λ (2,0) = 1, Λ (1,1) = −3/2, and other Λ (p,q) = 0. Note that 3.30 with this choice of parameters is similar to 3.16, except that ( i ρ 2 i ) now plays the role of the bifundamental mass w 2 , which controls the location of the lobes' apexes. The height of the apexes behaves as t hop ∼ 1/w, as already stated in the previous section. The difference is that the w 4 term is missing. Hence, the difference in free energies between the homogeneous and nonhomogeneous phases varies with the on-site occupation number like ∼ 1/ρ (l) , as desired. Note also that we chose Λ (1,1) = −3/2 to help achieve this specific functional dependence on ρ. Since the IR potential does not affect the EOM, the perturbations and the excited spectrum in section 4 are not affected by this change due to our choice of free boundary conditions in the IR. In App. A we analyze the mixed Neumann boundary conditions that follow from the variation of the action with the choice of IR potential (3.16), and show that the structure of the phase diagram of our model is qualitatively unaffected by this change.

Perturbations and Excitation Spectrum at Small
Hopping In this section, for small hopping parameter, we compute the perturbation around the background of both the Mott insulator phase and the non-homogeneous phase. We show that for Dirichlet boundary conditions at the IR wall, we find no zero modes in the Mott insulator phase, consistent with the existence of a gap. We show that an almost zero mode (E r h ) appears in the non-homogeneous phase, though we do not find a Goldstone mode consistent with superfluidity.
These EOMs are not independent. It can be shown that the radial derivative of the second equation in (4.34) is equal to a linear combination of the third equation and the fourth equation in (4.34). The number of integration constants hence is 5, due to the first order constraint.
In the remainder of this section, we consider an infinitesimal t hop region since finite t hop corrects the power (conformal dimensions) of the asymptotic expansion of perturbations in general, c.f. eq. (3.24). These corrections at finite t hop also make it difficult to compute two point functions numerically, due to possible correction terms to the nonnormalizable modes which are leading compared to the normalizable piece. The small t hop analysis will suffice to substantiate the statements that the phases inside the lobes of Fig. 3 are charged Mott insulators, while the phase outside does not have a gap of the order of the Coulomb parameter U , and hence can't be identified with a Mott insulator. The non-homogeneous phase will in particular show an almost zero mode, i.e. have a parametrically smaller gap. Note that we employ the same boundary conditions for the fluctuations that we used for the background configurations A cl A and φ cl .

Mott Insulator Phase
In Mott insulator phase, φ cl = t hop and A cl = 0. The solutions δA A , δφ I , and δφ R for finite hopping t hop are then given by where α p = 1 2 16q 2 t 2 hop + 1. δφ I has the asymptotic behavior δφ I ∼ ϕ + C 3 r αp− 3 2 (ϕ is constant) near the boundary. As seen in the t hop = 0 case, ϕ and C 3 can be understood as the source term and VEV, respectively. For small hopping parameter (α p → 1 2 ), the Bessel functions and their integrals are replaced by cos and sin functions. The solutions are given by where ϕ is constant. Note that the number of integration constants is 5 consistent with EOMs (4.34). To derive the holographic two point functions we impose the IR boundary conditions δA A | r=r h = 0, δφ R | r=r h = 0 and δφ I | r=r h = 0, i.e. our boundary conditions allow neither the chemical potential nor the zero bifundamental VEV to vary. This shows that ϕ is negligible at O(t hop ) and δφ (2) I,R = −ωδφ The Dirichlet boundary condition makes the differential operator describing the small fluctuations Hermitian and consequently all eigenfrequencies are real. This is expected in the Mott insulating phase. Recall that the AdS boundary expansion becomes δφ I,R ∼ δφ (1) I,R + δφ (2) I,R /r. The two point function of the operator dual to δφ I,R is then given by Thus, the real part of the two point function is the same as the imaginary part of it. These two point functions have a pole at ω = πr h n o (n o ≥ 1) [34]. However, we do not observe a peak at ω = 0, or low lying modes at very small ω; we do not find zero modes or near-zero modes in the spectrum of the Mott insulator phase at least for small t hop . We conclude that the Mott phase is gapped, with the gap of the order of the Coulomb repulsion, ∆ = πr h = πU . When we take into account finite t hop corrections, the position of the peak is corrected in the imaginary part of the Green's function. We can numerically analyze this correction by finding the zero of the constant part of δφ I under the above IR boundary conditions because this constant part is the non-normalizable mode of δφ I . The gap of the excitations decreases only slightly as t hop is varied from zero to the maximal value at the tip of the lobe, t hop,max = 0.473 (see Fig. 7). Hence, the mass gap of the excitations stays of the order of the Coulomb repulsion throughout the Mott phase.
In the Bose-Hubbard model [16], it is known that the low-lying excitations are described by the motion of a fundamental boson from a site to a neighboring site. To move a fundamental boson from a site to a neighboring lattice site costs energy U because of the repulsive Coulomb force between the fundamental bosons. The mass gap obtained above is consistent in order of magnitude with the mass gap ∆ = U of the Mott insulator phase in the actual Bose-Hubbard model.

Non-homogeneous Phase
For small t hop , the background in the non-homogeneous phase can approximated by  : Green and Orange curves show the ratio G R,I /ω dual to δφ R and δφ I as a function of ω, respectively. In Figure, σ I,R ≡ G I,R /ω, respectively. Both σ R,I behave like 1/ω at small ω larger than the energy of the almost zero mode. The peak near ω ∼ 0 does not almost move as t hop increases.
We consider fluctuations around this approximated background in order to calculate the spectrum close to the cusps at t hop = 0 in the phase diagram. Our system of fluctuations is similar to the one in [43], due to the coupling of scalar and longitudinal gauge modes. Since they do not admit analytic solutions, we solve the fluctuation EOMs numerically by the shooting method. We expand the fields at the hard wall as δF = n=0 δF (n) (r − r h ) n , where the parameters are fixed by 5 integration constants. The IR boundary conditions on the perturbations are chosen again to be Dirichlet boundary conditions δA A | r=r h = 0, δφ R | r=r h = 0, and δφ I | r=r h = 0. Setting r h = 40, the expansion around the hard wall is then specified by 2 parameters a, b, δφ I ∼ a(r − r h ), δφ R ∼ b(r − r h ), and δA A ∼ δF (1) (a)(r − r h ), where δF (1) (a) is given by δF (1) (a) = it hop a/ω. One can fix one of these constants to be 1 (a = 1 for example) by rescaling the perturbations. For small t hop , the asymptotic expansion at the AdS boundary becomes The two point functions are plotted in Fig. 8 numerically. We do not observe a peak at ω = 0 for G R but do observe a peak near ω = 0 for G I . The existence of a zero mode in the bifundamental correlator with itself would imply a Goldstone mode in the non-homogeneous phase, consistent with a possible interpretation as a superfluid phase. We however find that the almost zero mode does not change its position as t hop is varied, suggesting that it is not a Goldstone mode. We will comment on ways to introduce superfluidity and hence generate actual spontaneous symmetry breaking in sec. 6. The ratio G R,I /ω is plotted in Fig. 9 as a function of ω. In Fig. 9, σ I,R ≡ G I,R /ω. Both σ R,I behave like 1/ω at small ω below the respective gaps. In particular, σ I should be related with the conductivity of the Bose-Hubbard model because the current in the Bose-Hubbard model is defined by Note that the creation and annihilation operators are dimensionless, while t has dimension of energy, so J has the usual dimension 1.
The identification (4.43) can easily be derived from a conservation equation. Note that charge conservation, as usual, should imply the following continuity equation: (4.44) The right hand side is (minus) the discretized spatial derivative of the current at unit lattice spacing. Note that no lattice spacing appears here, as the lattice spacing a 1 in our lattice is set by t hop ∼ a −1 1 . We can use this relation to identify J l . For time dependent fields, we can look at the r component of Maxwell's equations which we so far neglected, as it is automatically solved for static configurations. Note that below we set q = 1 by a rescaling of the bifundamental fluctuations. 9 This EOM reads where j r is the bulk current associated to the scalar fields. To meaningfully talk about a spatial current we should consider a multi-side model with bi-fundamentals φ l connecting the l-th and (l + 1) − th site. In this case Keeping only the leading order r terms in the near boundary expansion of the fields as in (3.24) this EOM readsρ In this expression we have different leading behaviors for the scalar fields on the various sites, but the version of the Hubbard model we are considering has all t hop,j be equal to the same t hop , which we further can chose to be real. With this our equation of motion (4.47) can be compared to the continuity equation (4.44) to directly give (4.43).

Generalization to the n-site model
So far, we employed the two-site model mostly for computational simplicity, and found its physics to be rather similar to the Bose-Hubbard model. In this section, we briefly 9 In a more careful treatment the charge of the bifundamental under the axial gauge field combination would appear here, which, if we want to retain the standard normalization of the Maxwell term, will be related to the charge of the bifundamental at each site by q A = √ 2q.
introduce the generalization of our model to the n-site model. The action of a model with n sites is given by S = S kin + S matter as In the summation, if we consider a chain model, n + 1 is identified with 1. Other summations over different spatial lattices (triangle, honeycome, Kagome etc.) are straightforward to introduce. We can take the following linear combination to extract the diagonal gauge field V µ as The Maxwell kinetic term is rewritten as The covariant derivative which appears in S matter can be rewritten as D (l) = ∂ µ − iqA Aµ . The Maxwell term of the diagonal gauge field V µ then decouples from the remaining part of the action as A , . . . , A (n−1) A , φ l , . . . , φ n ]. (5.52) This implies that the free energy is of the form F = µ i ρ (i) + E(ρ (i) , t hop ) like (3.28).
For the chain the physics will be similar to the two-site model; in particular, the levelchanging phase transitions will work in the same way, and the phase diagram will be qualitatively unchanged. It would be interesting to explore the phase structure of this model for different lattice configurations and/or beyond-nearest-neighbor hoppings.

Discussion
In this work we have analyzed a holographic dual of the Bose-Hubbard model based on U(1) gauge fields localized on gapped AdS 2 hard wall space-times, which are connected to each other by bifundamentals charged under the respective gauge groups. We have shown that the model admits a good one-to-one holographic dictionary with the operators and parameters showing up in the Bose-Hubbard model, that a model based on two sites already reproduces the lobe-like phase structure in the chemical potential -hopping parameter (µ b − t hop ) plane, that the Mott insulating states have a natural excitation gap of the order of the Coulomb repulsion parameter, and that the transition to the inhomogeneous phase at the cusp points at zero hopping where the lobes meet is second order.
Our holographic model exhibits several differences from the Bose-Hubbard model: Except at the cusp points, the transition to the inhomogeneous phase is generically first order, and no Goldstone modes associated with the spontaneous breaking of the difference U(1) gauge groups are present in the inhomogeneous phase. Instead we find a near-zero mode at unnaturally small frequency appearing in the inhomogeneous phase near the cusp points. This near zero mode does not become a zero mode at vanishing hopping, and hence does not appear to be a Goldstone mode. Whether this conclusion continues to hold in other parts of the phase diagram will be the topic of a future, more complete investigation of the fluctuation spectrum. Finally, the overall vector U(1) in our model is not broken by the hopping, while it is in the condensed state of the Bose-Hubbard model [35].
In view of these differences to the Bose-Hubbard model, the two most interesting questions for future work will be: how to achieve a continuous phase transition between Mott and inhomogeneous phases everywhere along the phase boundaries, and how to achieve superfluidity in the inhomogeneous phase. A continuous phase transition is generically expected in holographic models with spontaneous breaking of U(1) symmetries in the bulk [51]. A hint to the issue is the AdS 2 hard wall geometry we are using, which is not a solution to Einstein's equations, so the first order nature of the phase transition may be an artifact of this shortcoming. An obvious improvement would be to use an AdS 2 hard wall-like geometry as is typically induced on the worldvolume of effectively two-dimensional probe branes embedded into higher dimensional AdS solitons [61]. In these geometries, the radial direction would cap off smoothly, and hence smoothen out the phase structure. Furthermore, such a bottom-up model would be more easily connected to top-down constructions of the bosonic and fermionic Hubbard models (see below).
Even in the hard wall-like geometry induced on the probe brane, however, the transition may still be first order, due to the transition taking place in the presence of a finite source term switched on, the hopping parameter. In this case we will need to introduce additional superfluid order parameters which should condense spontaneously, in order to achieve actual superfluidity in the inhomogeneous phase. There are basically two options: either we can couple fundamental scalars to the U(1) gauge field at each site, or introduce additional bifundamentals. In the former case the U(1) will break spontaneously if the local charge density at a particular site exceeds a critical value given by the charge and mass of the fundamental at that site, while in the latter the different charge density between the two sites to which the bifundamental is connected will be important. Specializing to the two-site model, we for example can break the vector U(1) spontaneously by introducing an additional fundamental at either of the sites. On the other hand, the two-site construction used in the main part of this paper, where we set the normalizable (i.e. in a sense spontaneously generated) part of the kinetic energy VEV to zero by our UV boundary conditions, could easily be amended by introducing a second bifundamental with the same charge as the first one (but maybe different mass), which again connects both sites. In this case however we would require this second bifundamental to condense with zero source term, i.e. not switch on a hopping parameter for it. The combined dynamics of this extended two-site model would then exhibit spontaneous breaking from the second bifundamental, while an explicit hopping VEV would be generated from the "hopping bifundamental". We are planning to present results on these different possibilities, as well as on other improvements of the model, in a follow-up work [62].
In this work we have mostly focused on a simple bottom-up construction. Here we would like to outline how to construct a top-down version of our model using the AdS 5 soliton [61]. We introduce a probe D5-brane on the AdS 5 soliton times S 5 [61] without considering its back-reaction. 10 Recall that the D5 probe branes can not end at the tip of the soliton (hard wall) and the D5-brane has to come back at a turning point like in holographic QCD [44] (see also [19,20]). If we do this in the internal soliton directions, the brane may smoothly cap off before hitting the soliton. However, we can also introduce a D7-brane filling the 2+1 dimensional uncompactified directions of the cigar and wrapping the whole S 5 . When we have a D7-brane sitting at the tip of the AdS soliton, the D5-brane can end on this D7-brane. Such a D7-brane has been identified as the holographic level-rank dual to Chern-Simons theory in [45]. In summary, these two world-volume theories seem to be good top-down constructions. We can expect the fundamental fermions corresponding to the 3-5 string modes to have the following low 10 Focusing on the 3-5 string modes where the ground state of the massless mode is obtained from R-sectors, such a state is given by the fundamental fermions χ ia of U (N ) gauge symmetry on a site. So this D3-D5 model seems to be a good holographic dual to the Fermi-Hubbard model. However, note that in such top-down constructions one is usually forced to work in the strict large N limit, in which we do not expect many differences between the fundamental bosons and the fundamental fermions since large N numbers of particles can occupy states of the same energy as is the case in boson statistics. In other words, there is no restriction from the Pauli exclusion principle at large N . energy effective action at each site: where φ D3 is the transverse scalar, A 0 is the D3 gauge field, andÃ 0 is the D5 gauge field. Considering the background of the D5-brane gauge field [Ã 0 ] a a = µ a and [Ã 0 ] a b = t ab (a = b), we have the following hopping term and chemical potential from the third term in (6.53): This action is similar to semi-holographic fermions [46,47] in the absence of the second fermionic operator. The holographic dual to this field theory would then be a stack of D5 branes, possibly ending on D7 branes, which are separated from each other to reside on the different sites of our model. The nonabelian D5 brane gauge symmetry is higgsed to the U(1) subgroups in this process, and the off-diagonal components of the nonabelian D5 gauge field will become the bifundamentals in our bottom-up construction [19,20]. The construction with D7 branes at the soliton tip seems to have another advantage: In our hard wall model of the bottom-up construction, we added a potential in the IR boundary. In the top-down construction, the interactions between the different D5-branes will be deformed by the Chern-Simons terms on the hard wall. Moreover, the action of 5-7 string modes will introduce additional interactions at the D5/D7 intersections as pointed out in a slightly different set-up in Ref. [48]. In this reference, the D7-D7 are suspended from the AdS boundary, while the D7-brane sits at the tip of the AdS soliton at the large separation limit of D7-D7 in the above model. Such a top-down construction would be very useful, since we are not required to tune free parameters. We are planning to analyze these different top-down constructions in detail in the near future [62].
Generalizing our model to higher-dimensional Bose and Fermi Hubbard models would be interesting because the higher-dimensional Bose and Fermi Hubbard models are more difficult to analyze using field theoretic or numerical techniques. In higher dimensions, Monte Carlo simulations are typically used to analyze the ground and thermal states, but suffer, in the fermionic case, from a sign problem. Different lattices (tetrahedral, triangular, Kagome, etc.) lead to vastly different physics such as spin-charge fractionalization, frustration, quantum spin ices, spin liquids, etc.. It will be very interesting to investigate the possible phases of higher dimensional lattices in this model in future work. It would also be interesting to apply our model to disordered systems in one or higher dimensions, by e.g. randomizing the chemical potentials, the hopping parameter, or other parameters (such as the bifundamental mass [60] or parameters in the IR potential). In disordered systems in the bosonic case, the Bose glass-phase appears in the phase structure between the insulating and superfluid phases. A Bose glass is characterized by a vanishing gap and finite compressibility, but it is an insulator because localization occurs due to the random potential. In the real Bose-Hubbard model, the phase transition to the superfluid phase is known to occur only from the Bose-glass phase [35]. Disorder was introduced in AdS/CFT in Refs. [36,37,58,59,56,57,55]. For Gaussian disorder, the free energy of the disordered system can in particular be evaluated by introducing replica fields and by averaging over disorder: F = −log Z where log Z = (Z n − 1)/n as n → 0 [37].

A.1 Homogeneous Mott Insulator
The homogeneous phase is defined by the condition A t . The EOM of the fields φ and A (l) t are then diagonalized. The EOM of the fields (3.17) are solved analytically as where l = 1, 2. The coefficient t hop and the coefficient of the normalizable mode ϕ M ott ≡ −ϕ 0 are identified with the hopping parameter and vacuum expectation value (VEV) of the bi-local field b † i b j in Bose-Hubbard model, respectively. Note that the minus sign in front of ϕ 0 appears in the formula of VEV (see [49,50]  The boundary term from varying the action is required to vanish at the IR wall, giving rise to the mixed Neumann boundary condition (see also [30,31]) The above boundary condition is a generalized version (due to the IR potential) of a class of modified boundary conditions which can e.g. describe a metal/insulator phase transition [51,52]. The condition (A.56) can be solved for the VEV ϕ 0 , where Q 1/3 = (9 √ Λt hop + 6 + 72Λ 2 w 2 + 48Λ 3 w 3 + 9Λ(9t 2 hop + 4w 2 )) 1 3 . We can show that the real solution (A.57) of the (generally complex) three solutions of (A.56) minimizes the on-shell action below (free energy). Moreover, the two complex solutions are non-zero at t hop = 0 and are hence not preferred (we already chose a gauge in which φ is real). We plot the VEV ϕ M ott (= −ϕ 0 ) of the real solution (A.57) as the function of t hop in Fig. 10. When w 2 < 0, ϕ 0 changes sign, which is not physically preferred either -we would like the VEV of the kinetic energy operator to be positive for positive hopping parameter. We hence choose a positive mass w 2 in what follows.
In the homogeneous phase, S kin is finite and additional counter-terms are not needed. The free energy is given by

A.2 Non-homogeneous Mixed State
For the non-homogeneous case A t there is no analytic solution to the EOM (3.17) in general. We rely on numerical methods to solve the EOM (3.17).
We obtain the following asymptotic behaviors of the solutions: where δρ = ρ (1) − ρ (2) and α t = (−1 + 1 − 4q 2 δρ 2 )/2. We include the subleading corrections in the above asymptotic expansion because subleading corrections are important even for large r. For example, we find finitely many correction terms to the hopping term in the expansion of φ for a specific choice of q 2 δρ 2 . The situation in the gauge field expansion is the same. The condition to stay above the BF bound, i.e. to keep real α t , is To keep real α t , namely, we restrict to the case |δρ| = 1 for the integer occupations ρ (l) and for our choice of charge (3.18). Since at zero hopping the IR potential only contributes an additive shift to the free energy, we still identify t hop andφ ≡ ϕ 0 (1 − 2∆ φ ) with the hopping term and VEV of this operator in the non-homogeneous phase. The additional factor of 1 − 2∆ φ (= −1/5) appears from the requirement of Ward identities in the field theory side [53,54]. We still relate the imaginary part of VEV with the current of our theory as discussed in section 4. α t encodes the anomalous dimension of this operator due to the interactions in the non-homogeneous phase. Our numerical procedure is as follows: instead of shooting from the IR to the UV we shoot from the UV to the IR, and vary the VEV ϕ 0 until the IR boundary condition (A.56) is satisfied. We plot the VEVφ in the non-homogeneous phase as a function of t hop in Fig. 11. The VEVφ increases as w 2 increases. As expected, the IR boundary condition changes with w 2 , and hence the whole solution and in particular the UV VEV varies. The results in the non-homogenous phase should be compared with those in the Mott insulator phase for w 2 = 1 (Dashed line) -in this case the VEV is much larger. Since we expect a small VEV in the Mott phase, it would hence be preferable from a model building point of view to choose the vanishing VEV boundary conditions from the body of this paper in the Mott phase, and the mixed Neumann boundary conditions in the inhomogeneous phase. Figure 12: Phase structure of the two-site model for the mixed Neumann boundary conditions, for r h = 40, w 2 = 1, and Λ = 1. The chemical potential µ b is defined in (2.14). Inside the lobes, the charge density on both sites is equal, analoguous to the situation in sec. 4). We identify this homogeneous phase with the Mott insulating phase of [35]. As t hop is increased, there are regions (which we call inhomogeneous phases) where the non-homogeneous states are thermodynamically favored. The basic structure of the phase diagram is unchanged from the case discussed in the body of this paper.
The bi-fundamental's action is UV divergent in the non-homogeneous case δρ = ±1. To cancel this divergence, the following counter-terms should be added : The free energy is then given by the holographically renormalized action as F = −(S kin + S matter + S cut + S cut,2 )/β. (A.62) Note that the diagonal gauge field A V = A (1) + A (2) decouples from the remaining parts.  Figure 13: The free energy is plotted as the function of µ b for t hop = 0.6. The color coding is as follows: Green lines correspond to the thermodynamically favored phase, while red lines to unstable phases of higher free energy. Solid lines correspond to the Mott phase, dashed lines to the inhomogeneous phase. The Mott insulator phase is hence the stable ground state for most values of µ b at this value of the hopping parameter. The Mott insulator is unstable in some regions of the chemical potential, between the lobes, where the inhomogeneous phase takes over. The free energy hence reflects the lobe structure in Fig. 12. For finite hopping parameter, first order level-changing phase transitions are hence found between homogeneous phase and non-homogeneous phases. Figure 14: Difference of the free energy F − F M ott between the non-homogeneous phase and the Mott insulator phase (ρ = 2) as the function of t hop . From the left to the right, the free energy is plotted for fixed µ b /U = µ/U + 1/2 = −1, −1.125, −1.25, −1.375, −1.5, respectively. It implies that for the cusp point at µ b /U = −1, the non-homogeneous phase is always favored when t hop = 0. µ b /U = −1.5 is a particle-hole symmetric point. After crossing this point, one attains the curves again in reversed order until one reaches the next cusp point at µ b /U = −2.
Thus, F can be rewritten as The above formula shows that the energy E( i ρ (i) , δρ, t hop ) is independent of the chemical potential at zero temperature at least. The phase structure of the two-site model is plotted numerically in Fig. 12 for r h = 40, w 2 = 1, and Λ = 1. Note that the chemical potential µ b is defined in (2.14). The lobeshaped phase structure of the Bose-Hubbard model [35] is also realized with the mixed Neumann boundary conditions, c.f. Fig. 12. For finite t hop , there are regions where the inhomogeneous state is favored. Furthermore, the non-homogeneous phase extends to the µ-axis at µ b /U ≡ µ/U + 1/2 = 0, −1, −2 as seen in Fig. 12. A small VEV ∼ φ is expected near these critical points on the µ b -axis.
The amplitudes of all lobes can be changed by arranging the parameter w in the IR potential accordingly. The amplitude is decreased as w becomes large. Note that when µ decreases, the tips of lobes are not changed in our model even with the mixed Neumann boundary conditions Fig. 12, while the tips of the lobes decrease as 1/ρ (1) in the actual Bose-Hubbard model [35].This shows that the height function of the lobes solely depends both on the choice of IR potential, as well as on the choice of boundary conditions. The free energy is plotted as a function of µ b for fixed hopping parameter t hop = 0.6 in Fig. 13. There, the green lines are the thermodynamically prefered phases, while the red lines are the nonprefered phases with higher free energy, and the solid line corresponds to the Mott phase, while the dashed line corresponds to the inhomogeneous phase. The level-crossing phase transitions of first order are found between the homogeneous phase and non-homogeneous phase at finite hopping, in complete analogy to the free boundary conditions employed in the body of the paper. F for the inhomogeneous phase is plotted as a function of t hop in Fig. 14. Again, a second order phase transition is found near t hop = 0 in Fig. 14. However, a first order phase transition is found between the Mott phase and the non-homogeneous phase for all values of the hopping parameter, except for t hop = 0. This is also in accordance with the findings in the main part of this paper. We conclude that we find no qualitative difference in the phase structure of our model between the boundary conditions we discuss here and the boundary conditions we use in the body of the paper.