1 Introduction↩︎

Understanding the ground state energy of interacting quantum many–body systems is a central problem in mathematical physics. Although a complete description is generally out of reach, rigorous results can be obtained in suitable asymptotic regimes. One particularly tractable regime is the dilute limit, where the particle density \(\rho\) is sufficiently small compared to the interaction scale. In this setting, Bose gases exhibit a remarkable universality: to leading order, the ground state energy depends only on a single effective parameter, the two–body scattering length \(\boldsymbol{a}\).

For three–dimensional continuum Bose gases with repulsive interactions, this is seen in the so-called Lee-Huang-Yang-Wu [1], [2] formula \[e_0(\rho)= 4\pi \boldsymbol{a}\rho^2 \Bigl(1+\frac{128}{15\sqrt{\pi}}(\rho \boldsymbol{a}^3)^{1/2}+8(\frac{4\pi}{3}-\sqrt{3})\rho \boldsymbol{a}^3 \ln(\rho \boldsymbol{a}^3)+\dots\Bigr),\] which captures in the dilute regime \(\rho \boldsymbol{a}^3 \to 0\) the correct ground state energy per unit volume, up to corrections of order \(\boldsymbol{a}\rho^2 (\rho \boldsymbol{a}^3)\), which are expected to no longer be universal in \(\boldsymbol{a}\). The leading-order term was rigorously established by Dyson [3] as an upper bound and, over 40 years later, by Lieb and Yngvason [4] as a lower bound (see also [5]). An upper bound matching the second order term was established by Yau and Yin [6] (see also [7]–[9]), while a lower bound finishing the proof of the Lee-Huang-Yang conjecture was established by Fournais and Solovej in [10] (see also [11]). In [12], [13], a new, simpler proof that also establishes the free energy expansion in the positive temperature case was given (see also [14]–[16]). Finally, recently, Brooks, Oldenbrug, Saint Aubin and Schlein [17] established for the first time an upper bound that includes the third order term (so-called Wu term).

A natural question is whether this universality persists in discrete settings. Bose gases realized in optical lattices are described by lattice Hamiltonians, most prominently by the Bose–Hubbard model [18], [19], which has become a standard effective model for interacting bosons and has been extensively studied in both theory and experiment. Such systems arise on a variety of lattice geometries [20] beyond the simple cubic case, which motivates the consideration of general Bravais lattices. In this setting the single–particle dispersion and the associated low–energy kinematics depend strongly on the geometry and hopping amplitudes of the underlying lattice. In contrast to the continuum, both the interaction and the lattice structure influence the two–body problem, and it is not a priori clear whether the leading-order energy retains a universal form independent of the microscopic details.

In this work we show that such universality indeed survives on lattices as far as the leading order term of the energy is concerned. We consider interacting bosons on an arbitrary three–dimensional Bravais lattice with positive hopping amplitudes and on-site repulsive interactions and prove that, in the dilute limit \[e_0(\rho)=4\pi \boldsymbol{a}\rho^2 (1+O(\rho^{1/6}))\] where \(\boldsymbol{a}\) denotes the lattice scattering length (cf. Appendix 7). Thus all microscopic information — both the interaction strength and the lattice geometry — is absorbed into this effective parameter, and the leading-order energy is independent of other details of the underlying lattice. This provides a discrete analogue of Dyson and Lieb–Yngvason theorems for the Bose–Hubbard gas.

It is worth emphasizing that this universality is specific to the leading-order term. While in the continuum the second and (expectedly) third order terms also exhibit a universal structure depending only on the scattering length, in the lattice setting one expects higher-order terms to depend explicitly on the single-particle dispersion and hence on the geometry of the underlying lattice. In particular, beyond order \(\boldsymbol{a}\rho^2\) the energy is not determined solely by the scattering length. Our result therefore identifies the precise regime in which lattice effects are completely absorbed into this effective parameter. We expect that the same leading-order universality holds for more general short-range lattice potentials.

Related universality results have been obtained for fermionic systems, both in the continuum [21]–[25] and on the cubic lattice (with nearest neighbor hopping) [26], [27], where the leading-order energy is again determined by an appropriate scattering parameter. The available lattice proofs for fermions employ techniques that do not readily transfer to bosons. In particular, arguments based on Dyson-type lemmas have no direct counterparts. Therefore, in order to prove our main result, we rely on techniques that have been developed more recently in the context of continuous bosonic systems. For the lower bound we use a localization method coupled with Bogoliubov theory [28]–[30]. To make it work we need to develop estimates for the eigenvalues of Neumann Laplacians on general Bravais lattices. In fact, these bounds lead to the relative error of order \(O(\rho^{1/6})\) in the lower bound. The upper bound is an adaptation of the argument in [7] that allows to include lattice dispersion relations which are not radial.

The remainder of the paper is organized as follows. In Section 2 we introduce the lattice framework and state the main result. The upper bound is obtained via suitable grand-canonical trial states and the equivalence of ensembles. This is presented in Section 3. The corresponding lower bound is proved in Section 4. Several auxiliary results (including the discussion on the Bravais lattices and the lattice scattering length) are collected in the Appendices.
Data sharing is not applicable to this article as no new data were created or analyzed in this study.
Acknowledgements. The work of NM and MN was supported by the National Science Centre (NCN) grant Sonata Bis 13 (project number 2023/50/E/ST1/00439).

2 The model and the main result↩︎

We start by presenting basic definitions and objects of interest. We refer to Appendix 5 for the discussion concerning those aspects.

Our main object of interest is a three dimensional (monoatomic) crystal lattice. To define it we need to specify the underlying Bravais lattice and the neighborhood relation between the points on the lattice.

To this end we first fix three linearly independent vectors \(a_1\), \(a_2\), \(a_3 \in \mathbb{R}^3\) and a matrix \(A\) composed of those vectors as columns. We consider a Bravais lattice \(\Lambda\) defined as \[\label{def:Bravais95lattice} \Lambda = A\mathbb{Z}^3 = \left\{m_1 a_1 + m_2 a_2 + m_3 a_3 \colon m_1,m_2,m_3 \in \mathbb{Z}\right\}\tag{1}\] In this context vectors \(a_i\) are called the primitive (translation) vectors of the lattice \(\Lambda\).

For a given even number \(L \in 2\mathbb{N}\) we consider a finite version of the Bravais lattice 1 of size \(L\), denoted \(\Lambda_L\) and defined as \[\label{def:finite95lattice} \begin{align} \Lambda_L &= A\left(\mathbb{Z}\cap [-L/2,L/2]\right)^3 \\&= \left\{m_1 a_1 + m_2 a_2 + m_3 a_3 \colon m_1,m_2,m_3 = -\frac{L}{2}, -\frac{L}{2} + 1, \dots, \frac{L}{2} - 1, \frac{L}{2}\right\}. \end{align}\tag{2}\] We equip \(\Lambda_L\) with periodic boundary condition (i.e \(\Lambda_L \simeq A\mathbb{Z}^3/A((L+1)\mathbb{Z})^3\) and with the standard counting measure, hence we can define the (one particle) Hilbert space \(\mathcal{H}_L\) of a particle on the lattice \(\Lambda_L\) as \[\mathcal{H}_L = L^2\left(\Lambda_L\right)\] with the inner product \[\langle\psi, \varphi\rangle_{\mathcal{H}_L} = \sum_{x \in \Lambda_L} \overline{\psi(x)} \varphi(x).\] For a given natural number \(N\) we also define \(N\)-particle bosonic Hilbert space \(\mathcal{H}_L^N\) as \[\mathcal{H}_L^N = \bigotimes_{\text{sym}}^N \mathcal{H}_L,\] i.e. the functions of \(N\) variables \(x_1\), \(x_2\), …, \(x_N \in \Lambda_L\) symmetric under the permutations of those variables. The inner product on this space is defined for simple tensors as \[\left \langle\bigotimes_{j=1}^N \psi_j, \bigotimes_{j=1}^N\varphi_j \right \rangle_{\mathcal{H}_L^N} = \prod_{j=1}^N \langle\psi_j, \varphi_j \rangle_{\mathcal{H}_L},\] which can be extended to the whole \(\mathcal{H}_L^N\) by linearity.

Now we will define the neighbor relation on \(\Lambda\) and \(\Lambda_L\). Let \(D\) be a set of all "positive directions" \[D = \{m_1 a_1 + m_2a_2 + m_3a_3 \in \Lambda \setminus\{0\} \colon \text{ the first non-zero } m_j \text{ is positive}\}.\] Note that \(D \cup (-D) = \Lambda \setminus\{0\}\) and \(D \cap (-D) = \emptyset\). To each direction \(v \in D\) and its reverse direction \((-v)\) we will assign a weight \(t(v) = t(-v) \ge 0\). The neighborhood relation on \(\Lambda\) is defined as \[\label{neighbor95relation95def} x \sim y \Longleftrightarrow t(y - x) > 0.\tag{3}\] This relation is symmetric and equips both the infinite lattice \(\Lambda\) and the finite lattice \(\Lambda_L\) with the weighted graph structure, where in the latter \(y-x\) is understood in the sense of periodic boundary condition, i.e. as the element of the \(A\mathbb{Z}^3/A((L+1)\mathbb{Z})^3\) group.

We will make two additional assumptions. The first one is that \[\label{finite95range95assumption} \#\{v \in D \colon t(v) \ne 0\} < \infty,\tag{4}\] meaning that we only consider a finite distance hopping. This assumption is satisfied in most commonly encountered crystal systems in physics. For the future purposes we will also define a parameter \(R_0(t)\) called the hopping length as \[\label{def:hopping95length} R_0(t) = \min \{L \in 2\mathbb{N}\colon \forall_{x \sim 0} \; x \in \Lambda_L\},\tag{5}\] that is the smallest \(L\) such that all the neighbors of point \(x=0\) in the sense of 3 belong to \(\Lambda_L\). In the upcoming proofs we will consider only \(L \ge R_0(t)\) as this condition will allow us to capture all the possible hoppings within the finite volume.

To state the second assumption we will first denote \[\label{D95195def} D_1 := \{a_1, a_2, a_3\} \subset D.\tag{6}\] We will assume \[\label{t95195assumption} t(v) \ne 0 \text{ for } v \in D_1\tag{7}\] that is hopping along the primitive vectors of the lattice \(\Lambda\) is allowed.

The Bose-Hubbard Hamiltonian \(H_{N,L}\) of the \(N\) particle system, acting on \(\mathcal{H}_L^N\), is given by \[H_{N,L} = -\sum_{i = 1}^N\Delta_{i} + U\sum_{i<j}^N\delta_{x_i,x_j},\] where the first term is the kinetic energy. Here \(\Delta_{x_i}\) denotes the lattice (weighted graph) Laplace operator acting on the \(i\)-th particle: \[\Delta_i = \text{Id}\otimes \dots \otimes \Delta \otimes \dots \otimes \text{Id},\] where \(\text{Id}\) is the identity operator on \(L^2(\Lambda_L)\) and \(\Delta\) is a single particle Laplacian standing on the \(i\)-th position. We can write the action of \(\Delta\) explicitly: for \(u \in L^2(\Lambda_L)\) \[\label{def:discrete95laplacian} -\Delta u (x) = \sum_{y \sim x} t(y-x)\big(u(x) -u(y)\big) = \sum_{v \in D} t(v) \big(2 u(x) - u(x+v) - u(x-v)\big),\tag{8}\] where \(x,y \in \Lambda_L\) and \(x \sim y\) denotes the neighborhood relation. Here \(y - x\) is once again understood in the sense of periodic boundary condition. The second term in the Hamiltonian is the interaction energy with \(U > 0\) (i.e. the interaction is repulsive).

The ground state energy \(E_0(N,L)\) of the system is defined by \[\label{def:GSE} E_0(N,L) = \inf_{\substack{\psi \in \mathcal{H}_L^N \\ \|\psi\|=1}} \langle\psi, H_{N,L} \psi \rangle_{\mathcal{H}_L^N}.\tag{9}\] The inner product above is called the expectation value of the operator \(H_N\) in the state \(\psi\). In general, the expectation value of some operator \(T\) acting on the Hilbert space \(\mathcal{H}\) in the state (i.e. normalized vector) \(\psi \in \mathcal{H}\) is defined as \[\label{def:expectation} \langle T \rangle_{\psi} = \langle\psi, T \psi \rangle_{\mathcal{H}}.\tag{10}\] We will use this notation when there will be no ambiguity on which Hilbert space this expectation is evaluated.

We are interested in the ground state energy per unit volume, i.e. \[\label{def:ThermodynamicLimit} e_0(\rho) = \lim_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}}\frac{E_0(N,L)}{|\Lambda_L|}.\tag{11}\] Existence of this limit (under some more general assumptions and even in some broader setting) is known, we refer e.g. to [31] for the details. It is also known that \(e(\rho)\) is a continuous (up to the boundary \(\rho = 0\)) and convex function of \(\rho\).

In order to state the main theorems we introduce the scattering length of the potential which we will denote \(\boldsymbol{a}\) (see Appendix 7 for more details). For the on-site interaction potential that we are dealing with it is defined as \[\label{def:scattering95len} 8\pi \boldsymbol{a}= \frac{U}{U\gamma + 1},\tag{12}\] with \[\gamma = \frac{1}{2}|\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{dp}{\mathop{\mathrm{\varepsilon}}(p)}\] and \(\mathop{\mathrm{\varepsilon}}(p)\) being the lattice dispersion relation, given by \[\label{def:eps} \mathop{\mathrm{\varepsilon}}(p) = \sum_{v \in D} 2t(v)\big(1 - \cos (v \cdot p)\big) = 4\sum_{v \in D} t(v) \sin^2\left(\frac{v \cdot p}{2}\right), \quad p = (p_1,p_2,p_3) \in \widehat \Lambda.\tag{13}\] Here \(\widehat \Lambda\) is the Brillouin zone of the lattice \(\Lambda\) \[\widehat \Lambda= B\mathbb{T}^3 = \left\{b_1 t_1 + b_2t_2 + b_3t_3 \colon -\frac{1}{2} \le t_i < \frac{1}{2} \right\} \text{ with periodic boundary conditions},\] where \(\mathbb{T}^3 = [-\frac{1}{2},\frac{1}{2})^3\) is a three dimensional unit torus (this identification also allows to identify the Haar measure \(dp\) in the integral as the Lebesgue measure), \(|\widehat \Lambda|\) denotes the measure of this set and \(B\) is a matrix composed of columns \(b_1\), \(b_2\), \(b_3\) satisfying \[a_i \cdot b_j = 2\pi \delta_{i,j},\] that is these are the primitive vectors of the reciprocal lattice \(\Lambda^*\). We refer to Appendix 5 for a more detailed discussion. The expression \(\mathop{\mathrm{\varepsilon}}(p)\) can be seen as the eigenvalue corresponding to the eigenfunction \[\chi_p(x) = e^{ip\cdot x}, \quad x \in \Lambda\] of the discrete Laplacian defined in 8 . The sum in 13 is finite due to the assumption 4 . The formula 12 is derived explicitly in Appendix 7.

Due to the assumption 7 there exist \(p_0 > 0\) such that we have the estimate \[\label{eps95bound} \mathop{\mathrm{\varepsilon}}(p) \ge c|p|^2 \text{ for } |p| < p_0\tag{14}\] with \(|p|\) being the Euclidean norm of the vector \(p \in \mathbb{R}^3\) and with the constant \(c\) independent of \(p\) (one can take \(c = \frac{1}{2}\min \{t(a_1), t(a_2), t(a_3)\} > 0\)). This in particular implies that \(\gamma\) is well defined as the function \(1/\mathop{\mathrm{\varepsilon}}(p)\) is integrable near zero. Moreover, as we will see in the proofs of the following propositions, this inequality will be crucial for obtaining the desired results. The assumption 7 itself can be changed in such a way that 14 still holds true, for example assuming that some certain other hopping constants are non-zero (this however would require some modifications in the proofs). For the purpose of this paper we will stick to 4 as this is the simplest case when 14 holds.

The main theorem that we prove is the following

We will prove Theorem 1 by proving the upper and the lower bound separately. In fact the upper bound (see below) provides a better error estimate that the one coming from the lower bound.

As mentioned in the introduction, the proof of the upper bound will be based on [7], adapted to the lattice setting. The main point of this approach is that instead of constructing a sequence of states on \(\mathcal{H}_L^N\), i.e. states with fixed number of particles, we will work in the grand canonical setting. More precisely, we will construct a sequence of states \(\{\Psi_{L,N}\}_{L,N}\) on the Fock spaces \[\mathcal{F}_L :=\mathcal{F}(\mathcal{H}_L) = \bigoplus_{n=0}^\infty \mathcal{H}_L^n, \quad (\mathcal{H}_L^0 = \mathbb{C})\] with each having fixed average number of particles \(\langle\mathcal{N}\rangle_{\Psi_{N,L}} = N\). We will denote the grand canonical Hamiltonian (i.e. the second quantization of \(H_{N,L}\)) as \(H^\text{GC}_L\). This operator acts on the Fock space \(\mathcal{F}_L\), for \(\Psi = (\Psi^{(n)})_{n \in \mathbb{N}}\in \mathcal{F}_L\) its action is given by \[\label{def:gc95ham} (H^\text{GC}_L \Psi)^{(n)} = H_{n,L} \Psi^{(n)}.\tag{15}\] We are going to prove the following result.

The upper bound (with the same error term as above) follows then from the variational principle and the equivalence of ensembles. The adaptation of this well known argument for the discrete setting will be presented at the end of Section 3, after the proof of Proposition 2.

As for the lower bound we will use the method of dividing the large (thermodynamic) lattice \(\Lambda_L\) into smaller sub-lattices \(\Lambda_\ell\) with a properly chosen length scale \(\ell\). In the standard proof of the corresponding bound for the continuous (i.e. non-discrete) setting presented e.g. in [32] the length scale \(\ell\) is chosen is such a way that \(\rho^{1/3} \ll \ell \ll \rho^{-1/2}\), which allows to effectively use the Dyson lemma and obtain the desired result. Since we were unable to prove the discrete analogue of this lemma that would be of use to us, we propose a different approach and choose \(\ell \sim \rho^{-1/2}\), which is commonly known as the Gross-Pitaevskii length scale. Then we will use the method from [29] (also recently used in [30]) to get the operator inequality bounding \(H_{n,\ell}\) from below for certain values of \(n\) and \(\ell\). In the end we will use the obtained bound and the method from [28] to get the following result.

Let us stress that the worse error term than in the upper bound is a consequence of spectral estimates that we derive for general Bravais lattices with general neighbor relations. For example, this error can be improved to be of order \(O(\rho^{1/2} \ln \rho)\) if one considers cubic lattices with nearest neighbor hopping. The proof of Proposition 3 is given Section 4. This will finish the proof of Theorem 1.

In the rest of the paper we use the convention that \(C\) denotes a generic constant (independent of relevant parameters) which may change from line to line.

3 The upper bound↩︎

This section is devoted the proof of Proposition 2 and the corresponding upper bound. As mentioned before, the idea of the proof will follow the one in [7], but with some adaptation to the discrete setting. In particular, we will use methods that do not rely on the spherical symmetry of the system. The proof will be done in a few steps, many of them being by now standard in the field.

3.1 Momentum representation of the Hamiltonian↩︎

It will be convenient to rewrite the grand canonical Hamiltonian 15 in the formalism of creation and annihilation operators - we refer to [33] and [34] for more details concerning second quantization and Bogoliubov transformations. We will also use notation from Appendix 5.

We will fix \(L \in 2\mathbb{N}\) satisfying \(L \ge R_0(t)\) (recall definition 5 ) and work within the momentum representation. Following Appendix 5 (in particular section 5.5) we denote \[\widehat \Lambda_L = \left\{\sum_{j=1}^d m_j \frac{b_j}{L+1} \colon m_j=-\frac{L}{2},-\frac{L}{2} + 1,\dots, \frac{L}{2} - 1, \frac{L}{2}\right\},\] where \(b_j\) ae primitive vectors of the reciprocal lattice \(\Lambda^*\). For \(p \in \widehat \Lambda_L\) we denote by \(a_p\) and \(a_p^*\) annihilation and creation operators of a particle with a momentum \(p \in \widehat \Lambda_L\), that is \[a_p = a(\chi_p), \quad a^*_p = a^*(\chi_p),\] where \[\chi_p(x) = \frac{1}{|\Lambda_L|^{1/2}}e^{ip \cdot x}, \quad x \in \Lambda_L.\] Direct computation of the matrix elements of the one- and two-body operators in \(H_{N,L}\) in the above basis yields \[\label{Ham95mom} H_L^\text{GC}= \sum_{p} \mathop{\mathrm{\varepsilon}}(p) a^*_p a_p +\frac{U}{2 |\Lambda_L|}\sum_{p,q,k} a_{p+k}^* a_{q-k}^* a_q a_p,\tag{16}\] where indices \(p,q,k\) run over \(\widehat \Lambda_L\) and \(\mathop{\mathrm{\varepsilon}}(p)\) is defined analogously as in 13 but only for discrete values of \(p\): \[\label{def:eps95finite} \mathop{\mathrm{\varepsilon}}(p) = \sum_{v \in D} 2t(v)\big(1 - \cos (v \cdot p)\big), \quad p = (p_1,p_2,p_3) \in \widehat \Lambda_L.\tag{17}\] We also note that creation and annihilation operators satisfy standard canonical commutation relations \[[a_p,a_k] = [a_p^*, a_k^*] = 0, \quad [a_p, a_k^*] = \delta_{p,k}.\]

3.2 Construction of the states↩︎

We proceed to the construction of the trial state \(\Psi_L \in \mathcal{F}_L\) with \(\langle\mathcal{N}\rangle_{\Psi_L} = N\). As for this moment values of \(N\) and \(L\) are fixed, we will simplify notation and omit index \(L\) in some of the objects (i.e. \(\Psi_L = \Psi\) etc.).

Consider \(N_0\) satisfying \(0 \le N_0 \le N\) and for \(p \in \widehat \Lambda_L\setminus \{0\}\) consider a (finite) sequence of real numbers \((c_p) \subset \mathbb{R}\), such that \(|c_p| < 1\) and \(c_p = c_{-p}\). Define the state \(\Psi\) as \[\label{state95def} \Psi =\left(e^{-N_0/2}\prod_{p \ne 0} (1-c_p^2)^{1/4} \right) e^{\frac{1}{2}\sum_{p \ne 0} c_p a_p^* a_{-p}^* + \sqrt{N_0}a_0^*}\mathop{\mathrm{|0\rangle}},\tag{18}\] where \(p\) belongs to \(\Lambda^*\). Here the exponent should be understood as a notation for the proper series expansion. One can recognize that \(\Psi\) is a so-called Bogoliubov trial state, that is it is of the form \[\Psi = W\mathbb{U}^*\mathop{\mathrm{|0\rangle}},\] where \[W = W\left(\sqrt{N_0/|\Lambda|}\right) = e^{\sqrt{N_0}(a_0^* - a_0)}\] is the Weyl operator build upon constant function \(\sqrt{N_0/|\Lambda|}\) and \(\mathbb{U}\) is the Bogoliubov transformation given by \[\label{Bog95transform} \mathbb{U}= \exp \left(\sum_{p\in \widehat \Lambda\setminus \{0\}} -\text{artanh } c_p \left(a^*_p a^*_{-p} - a_p a_{-p}\right)\right).\tag{19}\] The state \(\Psi\) is normalized and conserves momentum, meaning that \[\label{mom95cons1} p \ne q \Rightarrow \langle a_p^*a_q\rangle_\Psi = 0 \qquad \text{and} \qquad p \ne -q \Rightarrow \langle a_pa_q\rangle_\Psi = \langle a_p^*a_q^*\rangle_\Psi= 0.\tag{20}\]

3.3 Computation of the energy↩︎

We will now find the energy \(\langle H_L^\text{GC}\rangle_\Psi\) of the system in the state \(\Psi\). This is a well-known computation, so will only state the main steps.

It follows from the properties of the Weyl and Bogoliubov transformations that \[\label{computation1} \langle a_0^*a_0\rangle_\Psi = N_0, \qquad \langle a^*_0a^*_0a_0a_0\rangle_\Psi= \langle a_0a_0\rangle_\Psi = N_0^2\tag{21}\] and for \(p \ne 0\) \[\label{computation2} \langle a_p^*a_p\rangle_\Psi = \frac{c_p^2}{1 - c_p^2}, \qquad \langle a_p^*a_{-p}^*\rangle_\Psi= \langle a_pa_{-p}\rangle_\Psi = \frac{c_p}{1 - c_p^2}.\tag{22}\] The same computation as in [35] leads, using 21 and 22 , to the following expression \[\label{energy} \begin{align} \langle H_L^\text{GC}\rangle_\Psi &= \sum_{p \ne 0} \mathop{\mathrm{\varepsilon}}(p) \langle a^*_p a_p\rangle +\frac{U}{2 |\Lambda_L|}\sum_{p,q \ne 0}\left[ \langle a_{p}^* a_{-p}^*\rangle \langle a_q a_{-q} \rangle + 2\langle a_{p}^* a_{p}\rangle \langle a^*_q a_{q} \rangle \right] \\&+ \frac{UN_0}{2|\Lambda_L|} \sum_{p \ne 0}\left[ 2\langle a_p a_{-p} \rangle + 4\langle a_p^* a_{p} \rangle \right] + \frac{UN_0^2}{2|\Lambda_L|} \\&= \sum_{p \ne 0} \frac{\mathop{\mathrm{\varepsilon}}(p) c_p^2}{1 - c_p^2} + \frac{U}{2|\Lambda_L|} \sum_{p,q \ne 0} \left[\frac{c_pc_q}{(1 - c_p^2)(1-c_q^2)} + \frac{2c_p^2c_q^2}{(1 - c_p^2)(1-c_q^2)}\right] \\&+ \frac{UN_0}{|\Lambda_L|} \sum_{p \ne 0} \left[\frac{c_p}{1 - c_p^2} + \frac{2c_p^2}{1 - c_p^2}\right] + \frac{UN_0^2}{2|\Lambda_L|}. \end{align}\tag{23}\]

As mentioned in the statement of Proposition 2 we will only consider states \(\Psi \in \mathcal{F}\) with fixed expectation value of particle numbers \(N := \langle\mathcal{N}\rangle_\Psi\) and later consider values of \(N\) and \(L\) such that \[\frac{\langle\mathcal{N}\rangle_\Psi}{|\Lambda_L|} = \frac{N}{L^3} \to \rho\] when \(N \to \infty\) and \(L \to \infty\). By 21 and 22 we have \[\langle\mathcal{N}\rangle_\Psi = N_0 + \sum_{p \ne 0}\frac{c_p^2}{1 - c_p^2}\] hence for considered values of \(N\) and \(L\) we get \[\label{rho95constraint} \rho = \frac{N_0}{|\Lambda_L|} + \frac{1}{|\Lambda_L|}\sum_{p \ne 0}\frac{c_p^2}{1 - c_p^2} + o(1)_{L \to \infty},\tag{24}\] where \(o(1)_{L \to \infty}\) denotes the expression converging to zero as \(L \to \infty\). We also observe that (as \(\rho\) is fixed) \[\label{rho94295constraint} \frac{1}{|\Lambda_L|^2}\left(N_0^2 + 2N_0\sum_{p \ne 0}\frac{c_p^2}{1 - c_p^2} + \left(\sum_{p \ne 0}\frac{c_p^2}{1 - c_p^2}\right)^2\right) = \rho^2 + o(1)_{L \to \infty}.\tag{25}\] From now on we will assume that parameters \(N_0\) and \(c_p\) are chosen in such a way that 24 is satisfied.

In the following part it will be convenient to rewrite the expectation value \(\langle H_L^{\text{GC}} \rangle_\Psi\) in 23 to express it in terms of the total density \(\rho\). The observations 24 , 25 and \[\sum_{p,q \ne 0} \frac{c_p^2c_q^2}{(1 - c_p^2)(1- c_q^2)} = \left(\sum_{p\ne 0} \frac{c_p^2}{1 - c_p^2}\right)^2\] show that \[\begin{align} \langle H_L^\text{GC}\rangle_\Psi &= \sum_{p \ne 0} \frac{\mathop{\mathrm{\varepsilon}}(p) c_p^2}{1 - c_p^2} + \frac{U}{2 |\Lambda_L|} \sum_{p,q \ne 0} \frac{c_pc_q}{(1 - c_p^2)(1-c_q^2)} \\&+ U\left(\rho - \frac{1}{|\Lambda_L|}\sum_{q \ne 0}\frac{c_q^2}{1 - c_q^2} + o(1)_{L \to \infty}\right)\sum_{p \ne 0} \frac{c_p + c_p^2}{1 - c_p^2} + \frac{U}{2}|\Lambda_L|\left(\rho^2 +o(1)_{L \to \infty}\right), \end{align}\] which after rewriting gives \[\label{energy95a} \begin{align} \langle H_L^\text{GC}\rangle_\Psi &= \sum_{p \ne 0} \frac{\mathop{\mathrm{\varepsilon}}(p) c_p^2}{1 - c_p^2} + \frac{U\rho(c_p + c_p^2)}{1 - c_p^2} \\&+ \frac{U}{2|\Lambda|}\sum_{p,q \ne 0} \frac{c_pc_q - 2c_q^2(c_p + c_p^2)}{(1 - c_p^2)(1-c_q^2)} + o(1)_{L \to \infty} \cdot \sum_{p \ne 0} \frac{c_p + c_p^2}{1 - c_p^2} \\&+ \frac{U}{2}|\Lambda|\left(\rho^2 +o(1)_{L \to \infty}\right) \end{align}\tag{26}\] We expect (and prove it in further steps) that with proper selection of coefficients \(c_p\) the value of \[\sum_{p,q \ne 0} \frac{c_p^2c_q^2}{(1 - c_p^2)(1-c_q^2)} = \left(\sum_{p \ne 0} \frac{c_p^2}{1 - c_p^2}\right)^2\] will be negligible in the thermodynamic limit in the dilute regime. To this end we rewrite (using symmetry of summation with respect to indices \(p\) and \(q\)) \[\begin{align} \sum_{p,q \ne 0} \frac{c_pc_q - 2c_q^2(c_p + c_p^2)}{(1 - c_p^2)(1-c_q^2)} &= \sum_{p,q \ne 0}\frac{(c_p - c_p^2)(c_q -c_q^2)}{(1 - c_p^2)(1-c_q^2)} - 3\sum_{p,q \ne 0}\frac{c_p^2c_q^2}{(1 - c_p^2)(1-c_q^2)} \\&= \left(\sum_{p \ne 0} \frac{c_p - c_p^2}{1 - c_p^2}\right)^2 - 3\left(\sum_{p \ne 0} \frac{c_p^2}{1 - c_p^2}\right)^2. \end{align}\] With this result expression 26 , after further simplifications \(\frac{c_p - c_p^2}{1 - c_p^2} = \frac{c_p}{1 + c_p}\) and \(\frac{c_p + c_p^2}{1 - c_p^2} = \frac{c_p}{1 - c_p}\) becomes \[\begin{align} \langle H_L^\text{GC}\rangle_\Psi &= \sum_{p \ne 0} \frac{\mathop{\mathrm{\varepsilon}}(p) c_p^2}{1 - c_p^2} + \frac{U\rho c_p }{1 + c_p} + \frac{U}{2|\Lambda_L|}\left(\sum_{p \ne 0} \frac{c_p}{1 - c_p}\right)^2 \\&- \frac{3U}{2|\Lambda_L|}\left(\sum_{p \ne 0} \frac{c_p^2}{1 - c_p^2}\right)^2 + o(1)_{L \to \infty} \cdot \sum_{p \ne 0} \frac{c_p}{1 - c_p} \\&+ \frac{U}{2}|\Lambda_L|\left(\rho^2 +o(1)_{L \to \infty}\right) \end{align}\] As a final modification of this expression we will replace the squared term in the first line with a linear one at expense of some another negligible term in the low density limit. The main idea is to add and subtract a \(\rho w(0)\) term to every element of the sum (also recall that \(w(0)\) is given in 88 ). We will denote \[\label{s95p} s_p = \frac{c_p}{1+c_p}\tag{27}\] and for convenience we will additionally define \(s_0 = 0\). Next we write \[\begin{align} \left(\sum_{p \ne 0} s_p \right)^2 = \left(\sum_{p \in \widehat \Lambda_L} s_p \right)^2 &= \left(\sum_{p \in \widehat \Lambda_L} \big(s_p + \rho w(0)\big)\right)^2 - 2|\Lambda|\rho w(0)\sum_{p \in \widehat \Lambda_L} s_p - |\Lambda_L|^2\rho^2 w(0)^2 \\&= \left(\sum_{p \in \widehat \Lambda_L} \big(s_p + \rho w(0)\big)\right)^2 - 2|\Lambda_L|\rho w(0)\sum_{p \ne 0} s_p - |\Lambda_L|^2\rho^2 w(0)^2. \end{align}\] Eventually we obtain the following expression for the energy in the state \(\Psi\) \[\label{energy95final} \begin{align} \langle H_L^\text{GC}\rangle_\Psi &= \sum_{p \ne 0}\left( \mathop{\mathrm{\varepsilon}}(p)\frac{c_p^2}{1 - c_p^2} + U\rho\frac{c_p}{1 - c_p} - U\rho w(0)\frac{c_p}{1 + c_p}\right) - \frac{U}{2}|\Lambda_L|\rho^2 w(0)^2 \\& + \frac{U}{2|\Lambda_L|}\left(\sum_{p \in \widehat \Lambda_L} \big(s_p + \rho w(0)\big)\right)^2- \frac{3U}{2|\Lambda_L|}\left(\sum_{p \ne 0}\frac{c_p^2}{1 - c_p^2}\right)^2 \\&+ o(1)_{L \to \infty} \cdot \sum_{p \ne 0} \frac{c_p}{1 - c_p} + \frac{U}{2}|\Lambda_L|\left(\rho^2 + o(1)_{L \to \infty}\right). \end{align}\tag{28}\]

3.4 Minimalization↩︎

We proceed to the minimalization procedure. For a start, we are interested in term-by-term minimization of the sum in the first line of 28 , that is the want to minimize: \[\label{min95target} \mathop{\mathrm{\varepsilon}}(p) \frac{c_p^2}{1 - c_p^2} + U\rho \frac{c_p + c_p^2}{1 - c_p^2} - U\rho w(0)\frac{c_p}{1 + c_p}.\tag{29}\] As mentioned before, it will turn out that the quadratic terms (second line in 28 ) will be negligible with the selection of \(c_p\) minimizing this expression. In a more concrete manner, we will start with proving the following following Lemma.

Having minimized the local part of the energy 28 now we will show that the remaining parts are negligible in the dilute limit \(\rho \to 0\).

From Lemma 1 and Lemma 2 we deduce the following corollary concerning the energy.

3.5 Thermodynamic limit↩︎

Now we pass with the expression ?? divided by the volume \(|\Lambda_L|\) to the thermodynamic limit. We will denote this limit as \[\lim_{L \to \infty} \frac{1}{|\Lambda_L|}\langle H_L^\text{GC}\rangle_\Psi = e_\Psi\] The sums in ?? are Riemann sums of a continuous function on \(\widehat \Lambda\), hence they converge to the integrals of the proper expression. More precisely, we obtain \[\label{energy95thermo} \begin{align} e_\Psi &= |\widehat \Lambda|^{-1}\int_{\widehat \Lambda}\frac{1}{2} \left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right)\right)dp \\&+ \frac{U}{2}\rho^2 - \frac{U\rho^2}{2}w(0)^2 + O(\rho^3) \end{align}\tag{41}\] To see the dependence on the scattering length \(a\) we recall the definition 12 and write \[U = 8\pi \boldsymbol{a}(1 + U\gamma),\] hence \[\frac{U}{2}\rho^2 = 4\pi \boldsymbol{a}\rho^2 + 4\pi \boldsymbol{a}U\gamma \rho^2,\] Next note that (e.g. from 88 ) \[4\pi \boldsymbol{a}= \frac{w(0)}{2\gamma},\] so we can rewrite the last line of 41 (besides the error term) as \[\begin{align} \frac{U}{2}\rho^2 - \frac{U\rho^2}{2}w(0)^2 &= 4\pi \boldsymbol{a}\rho^2 + 4\pi \boldsymbol{a}U\gamma \rho^2 - \frac{U\rho}{2}w(0)^2 \\&= 4\pi \boldsymbol{a}\rho^2 + \frac{U\rho^2}{2}w(0)(1 - w(0)) \\&= 4\pi \boldsymbol{a}\rho^2 + \frac{1}{2} U^2\gamma (1 - w(0))^2\rho^2. \end{align}\] Recalling also that \(\gamma\) is given by the integral of the function \(\frac{1}{2\mathop{\mathrm{\varepsilon}}(p)}\), this additionally allows to rewrite 41 as \[\label{energy95thermo2} \begin{align} e_\Psi &= |\widehat \Lambda|^{-1}\int_{\widehat \Lambda}\frac{1}{2} \left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \right)dp \\&+ 4\pi \boldsymbol{a}\rho^2 + O(\rho^3), \end{align}\tag{42}\] where we have joined the previous integral with the integral defining \(\gamma\).

Let us now focus on evaluating the above integral. We note that the integrand is a positive function, which follows from the computation \[\label{positive95integrand} \begin{align} &\;\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \\&= \frac{4U^2\rho^2w(0) - U^2\rho^2(1+w(0))^2}{\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} + \mathop{\mathrm{\varepsilon}}(p) + U\rho\left(1 + w(0)\right)} + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \\&= -\frac{U^2\rho^2(1-w(0))^2}{\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} + \mathop{\mathrm{\varepsilon}}(p) + U\rho\left(1 + w(0)\right)} + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \\&= \frac{U^2\rho^2(1-w(0))^2\left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) + U\rho\left(1 + w(0)\right)\right)}{2\mathop{\mathrm{\varepsilon}}(p)\left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} + \mathop{\mathrm{\varepsilon}}(p) + U\rho\left(1 + w(0)\right)\right)}>0. \end{align}\tag{43}\] Next, similarly as before, we will split the integration into two regions, but this time into regions \(\mathop{\mathrm{\varepsilon}}(p) \ge \delta\) and \(\mathop{\mathrm{\varepsilon}}(p) < \delta\) where \(\delta = \delta(\rho)\) will be chosen as a certain function of \(\rho\). For the first region if \(\rho\) is sufficiently small (with respect to \(U\) and \(w(0)\)) we can Taylor expand the square root up to the terms of order \(\rho^3\). Using the inequality \[(1+x)^{1/2} \le 1 + \frac{1}{2} x - \frac{1}{8}x^2 + \frac{1}{16}x^3\] we get \[\begin{align} \sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} &= \mathop{\mathrm{\varepsilon}}(p)\sqrt{\left(1 + \frac{2U\rho}{\mathop{\mathrm{\varepsilon}}(p)}\right)\left(1 + \frac{2U\rho w(0)}{\mathop{\mathrm{\varepsilon}}(p)}\right)} \\&= \mathop{\mathrm{\varepsilon}}(p)\sqrt{1 + \frac{2U\rho(1+w(0))}{\mathop{\mathrm{\varepsilon}}(p)} + \frac{4U^2\rho^2w(0)}{\mathop{\mathrm{\varepsilon}}(p)^2}} \\&\le \mathop{\mathrm{\varepsilon}}(p) + U\rho(1+w(0)) + \frac{2U^2\rho^2w(0) - \frac{1}{2}U^2\rho^2(1+w(0))^2}{\mathop{\mathrm{\varepsilon}}(p)} + \frac{C}{\mathop{\mathrm{\varepsilon}}(p)^2}\rho^3 \\&\le \mathop{\mathrm{\varepsilon}}(p) + U\rho(1+w(0)) - \frac{\frac{1}{2}U^2\rho^2(1-w(0))^2}{\mathop{\mathrm{\varepsilon}}(p)} + C(\delta)\rho^3, \end{align}\] where \(C > 0\) is a constant coming from the Taylor expansion, independent of \(\rho\) and \[\label{C95delta} C(\delta) = \frac{C}{\min_{\mathop{\mathrm{\varepsilon}}(p) \ge \delta} \mathop{\mathrm{\varepsilon}}(p)^2},\tag{44}\] which is another constant, dependent only on \(\delta\). Using 14 we also note that if \(\delta\) is sufficiently small then \[\label{C95delta95assmptotics} C(\delta) \lesssim \frac{1}{\delta^4}.\tag{45}\] Now we can estimate the integrand as follows \[\begin{align} &\quad \sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \\&\le \mathop{\mathrm{\varepsilon}}(p) + U\rho(1+w(0)) - \frac{\frac{1}{2}U^2\rho^2(1-w(0))^2}{\mathop{\mathrm{\varepsilon}}(p)} + C(\delta)\rho^3 -\mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)} \\&= C(\delta)\rho^3, \end{align}\] hence \[\label{integral62delta} \begin{align} &\;|\widehat \Lambda|^{-1}\int_{\mathop{\mathrm{\varepsilon}}(p) \ge \delta}\frac{1}{2} \left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)}\right)dp \\&\le |\widehat \Lambda|^{-1}\int_{\mathop{\mathrm{\varepsilon}}(p) \ge \delta}\frac{1}{2} C_\delta \rho^3dp \le |\widehat \Lambda|^{-1}\int_{\mathbb{T}^3}\frac{1}{2} C_\delta \rho^3dp =\frac{1}{2} C(\delta)\rho^3. \end{align}\tag{46}\] We will explicitly choose \(\delta = \delta(\rho)\) after the next step.

Now we proceed to the integral on the domain \(\mathop{\mathrm{\varepsilon}}(p) < \delta\). Using the coarea formula (see e.g. [36]) for some general and sufficiently regular function \(f\) integrable near zero we have \[\int_{\mathop{\mathrm{\varepsilon}}(p) < \delta} f(\mathop{\mathrm{\varepsilon}}(p))dp = \int_0^\delta f(r) \left(\int_{\mathop{\mathrm{\varepsilon}}(p) = r} \frac{d \mathcal{H}^2(\xi)}{|\nabla \mathop{\mathrm{\varepsilon}}(\xi)|}\right)dr,\] where \(\mathcal{H}^2\) is the two dimensional Hausdorff (surface) measure. Using formula 13 we see that \[\mathop{\mathrm{\varepsilon}}(p) \le C|p|^2(1 + |p|^2)\] for some constant \(C\) independent of \(p\), hence \[\mathcal{H}^2\left(\{\mathop{\mathrm{\varepsilon}}(p) = r\}\right) \le Cr.\] Moreover \[|\nabla \mathop{\mathrm{\varepsilon}}(p)| > c|p|\] for some constant \(c\), hence for \(\xi \in \{\mathop{\mathrm{\varepsilon}}(p) = r\}\) we have \[\frac{1}{|\nabla \mathop{\mathrm{\varepsilon}}(\xi)|} \le \frac{C}{\sqrt{r}}\] and therefore \[\int_{\mathop{\mathrm{\varepsilon}}(p) = r} \frac{d \mathcal{H}^2(\xi)}{|\nabla \mathop{\mathrm{\varepsilon}}(\xi)|} \le Cr^{1/2}.\] As a result, if the function \(f\) is non-negative, we get \[\int_{\mathop{\mathrm{\varepsilon}}(p) < \delta} f(\mathop{\mathrm{\varepsilon}}(p))dp \le C \int_0^\delta r^{1/2} f(r) dr.\] We are going to use this observation for the integrand as in 42 , that is \[f(\mathop{\mathrm{\varepsilon}}(p)) = \frac{1}{2} \left(\sqrt{\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho\right)\left(\mathop{\mathrm{\varepsilon}}(p) + 2U\rho w(0)\right)} - \mathop{\mathrm{\varepsilon}}(p) - U\rho\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho^2}{2\mathop{\mathrm{\varepsilon}}(p)}\right).\] In 43 we have already noted that this function is positive. We have \[\begin{align} &\;|\widehat \Lambda|^{-1}\int_{\mathop{\mathrm{\varepsilon}}(p) \le \delta}f(\mathop{\mathrm{\varepsilon}}(p))dp \le C\int_0^{\delta} r^{1/2}f(r)dr \\&= C\rho\int_0^{\delta} r^{1/2}\left(\sqrt{\left(\frac{r}{\rho} + 2U\right)\left(\frac{r}{\rho} + 2U w(0)\right)} - \frac{r}{\rho} - U\left(1 + w(0)\right) + \frac{U^2(1 - w(0))^2\rho}{2r}\right)dr \\&= C\rho^{5/2}\int_0^{\delta/\rho} s^{1/2}\left(\sqrt{\left(s + 2U\right)\left(s + 2U w(0)\right)} - s - U(1+w(0)) + \frac{U^2(1-w(0))^2}{2s}\right)ds \\& \le C\rho^{5/2}\int_0^{+\infty} s^{1/2}\left(\sqrt{\left(s + 2U\right)\left(s + 2U w(0)\right)} - s - U(1+w(0)) + \frac{U^2(1-w(0))^2}{2s}\right)ds, \end{align}\] where in the second to last equality we have changed the variable \(r := \rho s\) and in the last inequality we used the fact the integrand is well defined and positive on \(\mathbb{R}_+\). Performing similar computation as in 43 one can check that this integral is convergent, in particular we can conclude \[|\widehat \Lambda|^{-1}\int_{\mathop{\mathrm{\varepsilon}}(p) \le \delta}f(\mathop{\mathrm{\varepsilon}}(p))dp \le C\rho^{5/2}.\] Combining this result with 46 and 44 for any \(\delta \le \rho^{1/8}\) we eventually get \[|\widehat \Lambda|^{-1}\int_{\widehat \Lambda}f(\mathop{\mathrm{\varepsilon}}(p))dp \le C\rho^{5/2}.\] This allows us to estimate \(e_\Psi\) in 42 as \[e_\Psi \le 4\pi \boldsymbol{a}\rho^2(1 + C\rho^{1/2}).\] This finishes the proof of Proposition 2.

3.6 Equivalence of ensembles↩︎

We will now recall the well known argument which shows how to use Proposition 2 in order to obtain the lower bound in Theorem 1. We follow the proof from [37] with some inspiration from [8]. Here, however, we will not assume that \(\frac{N}{|\Lambda_L|}\) is constant. We start with the following lemma.

Now we can proceed to the main problem. We first observe that trivially \[E_0^{\text{GC}}(N,L) \le E_0(N,L)\] as any canonical trial state with \(N\) particles can be lifted to the grand canonical one, occupying only the \(N\)-particle sector of the Fock space (in particular having \(N\) as the expected number of particles). It follows that \[\limsup_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}} \frac{E_0^{GC}(N,L)}{|\Lambda_L|}\le \limsup_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}} \frac{E_0(N,L)}{|\Lambda_L|} = e_0(\rho).\] It remains to prove \[\liminf_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}} \frac{E_0^{\text{GC}}(N,L)}{|\Lambda_L|}\ge e_0(\rho).\] To this end we introduce a variable \(\mu \in \mathbb{R}\) (that can be interpreted as the chemical potential) and for any normalized \(\Psi \in \mathcal{F}_L\) with \(\langle\mathcal{N}\rangle_{\Psi} = N\) we write \[\begin{align} \frac{\langle H\rangle_\Psi}{|\Lambda_L|} &= \frac{1}{|\Lambda_L|} \left[\mu\langle\mathcal{N}\rangle_{\Psi} + \langle H - \mu \mathcal{N}\rangle_{\Psi}\right] \\&= \frac{1}{|\Lambda_L|}\left[\mu N + \sum_{n=1}^\infty \|\Psi^{(n)}\|^2\left(\langle H_n\rangle_{\Psi^{(n)}} - \mu n\right)\right] \\&= \mu \frac{N}{|\Lambda_L|} + \sum_{n=0}^\infty \|\Psi^{(n)}\|^2\left(\frac{\langle H_n\rangle_\Psi}{|\Lambda_L|} - \mu \frac{n}{|\Lambda_L|} \right) \\&\ge \mu \frac{N}{|\Lambda_L|} + \sum_{n=0}^\infty \|\Psi^{(n)}\|^2\left(\frac{E_0(n,L)}{|\Lambda_L|} - \mu \frac{n}{|\Lambda_L|}\right) \\&\ge \mu \frac{N}{|\Lambda_L|} + \sum_{n=0}^\infty \|\Psi^{(n)}\|^2\left(e_0\left(\frac{n}{|\Lambda_L|}\right) - \mu \frac{n}{|\Lambda_L|}\right) \\&\ge \mu \frac{N}{|\Lambda_L|} + \sum_{n=0}^\infty \|\Psi^{(n)}\|^2 \inf_{\tilde{\rho} \ge 0}\left(e_0\left(\tilde{\rho}\right) - \mu \tilde{\rho}\right) \\&= \mu\frac{N}{|\Lambda_L|} + \inf_{\tilde{\rho} \ge 0}\left(e_0\left(\tilde{\rho}\right) - \mu \tilde{\rho}\right), \end{align}\] where in one of the steps we have used the above Lemma. We also recognize \[\inf_{\tilde{\rho} \ge 0}\left(e_0\left(\tilde{\rho}\right) - \mu \tilde{\rho}\right) = -e_0^*(\mu),\] where \(e_0^*(\mu)\) is the Legendre transform of \(e_0(\tilde{\rho})\) (we also use notation \(\tilde{\rho}\) in order not to confuse it with \(\rho\) fixed in the statement of the Theorem 1). Since \(\Psi \in \mathcal{F}_L\) above was arbitrary, we conclude \[\frac{E_0^{\text{GC}}(N,L)}{|\Lambda_L|} \ge \mu\frac{N}{|\Lambda_L|} - e_0^*(\mu).\] Taking the limes inferior of both sides gives \[\liminf_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}} \frac{E_0^{\text{GC}}(N,L)}{|\Lambda_L|} \ge \mu \rho - e_0^*(\mu).\] Furthermore, as the left hand side is independent of \(\mu\), we additionally get \[\liminf_{\substack{N \to \infty \\ L \to \infty \\ N/|\Lambda_L| \to \rho}} \frac{E_0^{\text{GC}}(N,L)}{|\Lambda_L|} \ge \sup_{\mu \in \mathbb{R}} \left(\mu \rho - e_0^*(\mu)\right) = e^{**}_0(\rho) = e_0(\rho),\] where the last equality follows from the fact \(e_0\) is a convex and continuous (up to a boundary) function and for such functions Legendre transform is an involution (i.e. \(f^{**} = f\), see e.g. [8] for a simple proof). The ends the proof of the upper bound in Theorem 1.

4 The lower bound↩︎

In this section we will prove Proposition 3. We will follow the strategy described at the beginning of the paper.

4.1 Division into sub-lattices↩︎

Similarly as in the proof of the lower bound in the continuous setting (see e.g. [32]) we will divide the large (thermodynamic) lattice \(\Lambda_L\) into smaller ones. The upcoming lemma is a well-known result, here will give a proof based on [38] adapted to the lattice setting. Beforehand, in analogy to the definition 9 , we will denote the ground state energy of the \(N\) particle system in the box of side-length \(L\) with Neumann Laplacian as \[E_0^\text{Neu}(N,L) = \inf_{\substack{\psi \in \mathcal{H}^N_L\\\|\psi\| = 1}} \langle\psi, H_{N,L}^\text{Neu}\psi\rangle\] with \[H_{N,L}^\text{Neu}= -\sum_{i = 1}^N\Delta_{\Lambda_L,i}^\text{Neu}+ U\sum_{i<j}^N\delta_{x_i,x_j}\] and \(\Delta^\text{Neu}\) defined (in accordance to the Appendix 6.2) as \[-\Delta^\text{Neu}_{\Lambda_L} u(x) = \sum_{\substack{y \sim x\\y \in \Lambda_L}} t(y-x)\big(u(x) - u(y)\big), \quad u \in L^2(\Lambda_L)\] In the following proof we will use the fact the ground state of the Hamiltonian \(H_{N,L}\) with either periodic or Neumann Laplacian over the symmetric wave functions is the same as the ground state over all wave functions, in particular the ground state energies are the same with or without imposing the symmetry constraint. This statement for the continuous case is proven e.g. in [39], for discrete systems this proof is also valid with some straightforward modifications.

Sometimes it is useful to express the inequality from the Lemma by the coefficients of relative number of boxes, i.e. we replace \(c_n\) with \((\frac{\ell+1}{L+1})^3c_n\). Under this replacement we get the following corollary.

4.2 Neumann eigenvalues↩︎

In order to use the above corollary successfully we need to understand the spectrum of the Neumann Laplacian \(-\Delta^\text{Neu}_{\Lambda_L}\) on the lattice \(\Lambda_L\). In full generality it won’t be possible to derive explicit formula for the eigenvalues, hence we will need to estimate them in a proper way. Before doing that we will derive the explicit form of the spectrum for a very special case of the neighborhood relation on the lattice \(\Lambda\).

Observe the similarity of the expression \(\varepsilon^{\text{Neu}}_\text{special}(k)\) to the dispersion relation 13 or its version 17 for a finite lattice. Recalling that \(\pi \delta_{i,j} = \frac{1}{2} a_i \cdot b_j\) we can write \[\label{Neu95per95relation} \varepsilon^{\text{Neu}}_\text{special}(k) = \mathop{\mathrm{\varepsilon}}\left(k_1\frac{b_1}{2(L+1)} + k_2\frac{b_2}{2(L+1)} + k_3\frac{b_3}{2(L+1)}\right) = \mathop{\mathrm{\varepsilon}}\left(\frac{1}{2(L+1)}Bk\right).\tag{55}\]

For a more general neighborhood relation in \(\Lambda\) the approach used in the proof of Lemma 5 will be unsuccessful as in general the graph \(\Lambda_L\) with edges corresponding to the Neumann Laplacian will not have a structure of a graph product of one-dimensional graphs. However we can still show that in the general case the Neumann eigenvalues are in some sense close to the periodic eigenvalues and derive some bounds. The following Lemma will give the precise statement of those observations and will be used in the next section to prove an appropriate bound on the Hamiltonian.

The method used in the proof of point \((c)\) above can be used to obtain the following useful corollary.

4.3 Analysis on sub-lattices↩︎

In this section we will work with the eigenbasis of the Neumann Laplacian \(-\Delta^\text{Neu}_{\Lambda_\ell}\) for the box of size \(\ell\), denoted as \(\{\psi_k\}_{k \in \{0,1\dots,\ell\}^3}\). Such basis exists as this is a self-adjoint operator on a finite dimensional space. Moreover it is easy to check that the constant function, here denoted as \(\psi_0\), is an eigenvector with eigenvalue zero. Furthermore, as the hopping constants defining \(\Delta^\text{Neu}_{\Lambda_L}\) are real, one has \(\overline{\Delta^\text{Neu}_{\Lambda_\ell} u(x)} = \Delta^\text{Neu}_{\Lambda_\ell} \overline{u(x)}\) (in other words this operator is a complexification of a symmetric operator over the real vector space \(L^2_\mathbb{R}(\Lambda_\ell))\)) and therefore the eigenfunctions \(\psi_k\) can be chosen as real-valued. Our goal is to prove the following proposition.

4.4 Conclusion↩︎

Now we can conclude the proof of Proposition 3, hence finishing the proof of the main Theorem 1. The following proof is based on [28] (which itself is similar to the proof given in [32]).

5 Fourier analysis on a lattice↩︎

5.1 General theory↩︎

In this subsection we will briefly recall the abstract setting in which the Fourier transform can be defined. We refer to e.g. [42] for more systematic approach.

Let \(G\) be locally compact abelian group, let \(\mu\) be Haar measure on \(G\) (that is \(\mu(gA) = \mu(A)\) for any measurable subset \(A \subset G\) and any element \(g \in G\)). We define the (Pontryagin) dual group \(\widehat G\) as \[\widehat G = \text{Hom}(G, S^1),\] where \(S^1\) is a unit circle. The group \(\widehat G\) is also abelian and locally compact, moreover group \(\widehat {\widehat G}\) is isomorphic to \(G\).

For a function \(f \in L^1(G,\mu)\) we define its Fourier transform \(\widehat f:\widehat G \to \mathbb{C}\) as \[\widehat f(\chi) = \int_G f(g)\overline{\chi(g)}d\mu(g).\] As the Haar measure on \(G\) is fixed, once can choose the normalization of the Haar measure \(\nu\) on \(\widehat G\) such that the following inversion formula (defined for a certain class of functions \(f\)) holds \[f(x) = \int_{\widehat G }\widehat f(\chi) \chi(x) d \nu(\chi).\] The Fourier transform also extends to the unitary operator from \(L^2(G,\mu)\) to \(L^2(\widehat G, \nu)\), which is sometimes referred to as the (generalization of) Plancherel theorem.

5.2 Infinite Bravais lattice↩︎

We will present the theory of the \(d\)-dimensional Bravais lattices in the context of the abstract Fourier analysis on groups. Similarly as in 1 we define Bravais lattice \(\Lambda\) as \[\Lambda = A\mathbb{Z}^d = \left\{\sum_{j=1}^d m_j a_j \; \colon \; m_j \in \mathbb{Z}\right\},\] where \(A\) is a \(d \times d\) invertible real matrix and vectors \(a_j\) are its columns. In this context vectors \(a_j\) are called the primitive (translation) vectors of the lattice \(\Lambda\). We treat \(\Lambda\) as an additive topological group with discrete topology and the standard counting measure as its Haar measure.

The reciprocal lattice \(\Lambda^*\) is then defined as \[\Lambda^* = \left\{y \in \mathbb{R}^d \colon y \cdot x \in 2\pi \mathbb{Z}\text{ for all } x \in \Lambda\right\}.\] One can check that \(\Lambda^*\) is also a Bravais lattice generated by primitive vectors \(b_j\), \(j=1,\dots,d\), defined by the relation \[a_i \cdot b_j = 2\pi\delta_{i,j}\] or equivalently \[b_j = \frac{1}{2\pi}(A^T)^{-1}e_j, \quad e_j \text{ - standard basis vector in } \mathbb{R}^d.\]

The Pontryagin dual group \(\widehat \Lambda\) in this context is called the Brillouin zone of the lattice \(\Lambda\). Since \(\Lambda\) is finitely generated every homomorphism \(\chi \in \text{Hom}(\Lambda,S^1)\) is uniquely determined by values \[\chi(a_i) := e^{i\theta_j}\] for some \(\theta_j \in \mathbb{R}\), hence every such \(\chi\) is of the form \[\chi = \chi_\theta, \quad \theta \in \mathbb{R}^d\] with \[\label{chi95theta95def} \chi_\theta\left(\sum_{j=1}^dm_ja_j\right) = e^{i\sum_j m_j\theta_j}.\tag{66}\] We note that if \(y \in \Lambda^*\) then \(\chi_\theta \equiv \chi_{\theta + y}\) as for any \(x = \sum_{j}m_ja_j\in \Lambda\) we have \[\chi_{\theta+y}(x) = e^{i\sum_j m_j(\theta_j + y_j)} = \chi_\theta(x) e^{ix \cdot y} = \chi_\theta(x)\] as \(x \cdot y \in 2\pi\mathbb{Z}\). It follows that the Brillouin zone \(\widehat \Lambda\) can be isomorphically identified as a quotient group \[\widehat \Lambda= \mathbb{R}^d/ \Lambda^*.\] We will further identify \(\widehat \Lambda\) as \[\label{Lh95def} \widehat \Lambda= B\mathbb{T}^d = \left\{\sum_{j=1}^d b_jt_j \; \colon \; t_j \in \left[-\frac{1}{2},\frac{1}{2}\right)\right\} \text{ with periodic boundary condition}.\tag{67}\] Here \(B\) is a matrix whose columns are the primitive vectors of the reciprocal lattice \(b_j\) and \(\mathbb{T}^d\) is the \(d\)-dimensional unit torus \(\left[-\frac{1}{2},\frac{1}{2}\right)\). This identification allows to identify the correct (in the sense of the Fourier inversion formula) Haar measure on \(\widehat \Lambda\) as the normalized Lebesgue measure on the set \(B\mathbb{T}^d\). With this fact we can write the Fourier transform formula on \(\Lambda\) and the inverse transform formula on \(\widehat \Lambda\) explicitly: for \(f \in L^1(\Lambda)\) and \(g \in L^1(\widehat \Lambda)\) we have \[\label{lattice95Fourier95transform} \widehat f(p) = \sum_{x \in \Lambda}f(x)e^{-ip\cdot x}, \quad p \in \widehat \Lambda\tag{68}\] and \[\label{inverse95transform} \check g(x) = |B\mathbb{T}^3|^{-1}\int_{B\mathbb{T}^3} g(p)e^{i p\cdot x}dp, \quad x \in \Lambda.\tag{69}\] By Plancherel theorem \(f \mapsto \widehat f\) and \(g \mapsto \check g\) extend to unitary maps between \(L^2(\Lambda)\) and \(L^2(\widehat \Lambda)\) and the extensions are each others inverses. We also note that \(|B\mathbb{T}^3| = |\det B|\).

5.3 Physical definition↩︎

In physics textbook one can encounter a different definition of the (first) Brillouin zone \(\widehat \Lambda\), namely that it is the Voronoi cell (in this context also called the Wigner-Seitz cell) of the point \(y = 0\) of the reciprocal lattice \(\Lambda^*\), that is \[V_0 = \{p \in \mathbb{R}^d \colon \|p\| = \min_{y \in \Lambda^*} \|p - y\|\}.\] Here \(\|\cdot\|\) is any norm in \(\mathbb{R}^d\). We will show that definition gives a rise to a certain identification of \(\widehat \Lambda\).

For \(y \in \Lambda^*\) denote by \(V_y\) the translation of \(V_0\) by the vector \(y\). This set is the Voronoi cell in \(\mathbb{R}^d\) based on the point \(y\). Next take a subset \(\tilde{V}_0 \subset V_0\) such that the sets \(\tilde{V}_y\) (translations of \(\tilde{V}_0\)) satisfy \[\label{Voronoi95conditions} \bigcup_{y \in \Lambda^*} \tilde{V}_y = \mathbb{R}^d, \quad \tilde{V}_{y_1} \cap \tilde{V}_{y_2} = \emptyset \text{ for } y_1 \ne y_2.\tag{70}\] Those conditions mean that \(\tilde{V}_0\) is \(V_0\) with some parts of the boundary removed.

For a point \(p \in \mathbb{R}^3\) denote as \(y(p)\) such (unique) point in \(\Lambda^*\) that \(p \in \tilde{V}_{py(p)}\). Define a relation \(\sim\) on \(\mathbb{R}^d\) as \[p_1 \sim p_2 \Longleftrightarrow p_1 - y(p_1) = p_2 - y(p_2).\] It it straightforward to check that it is an equivalence relation and that on the quotient \(\mathbb{R}^d\) it is possible to define the addition as \[[p_1]_\sim + [p_2]_\sim := [p_1 + p_2]_\sim.\] This makes \(\mathbb{R}^d/\sim\) an additive topological group with topology induced from \(\mathbb{R}^d\).

We will show that the groups \(\mathbb{R}^d/\sim\) and \(\widehat \Lambda=\text{Hom}(\Lambda,S^1)\) are isomorphic. The isomorphism is given by \[[p]_\sim \mapsto \chi_p(x):= e^{i p \cdot x}.\] It is well defined as if \(p_1 \sim p_2\) then \[p_1 - p_2 = y(p_1) - y(p_2) \in \Lambda^*\] hence \[\chi_{p_1}(x) = e^{i p_1 \cdot x} = e^{i(p_1 - p_2)\cdot x}e^{ip_2 \cdot x} = e^{ip_2 \cdot x} = \chi_{p_2}(x).\] Surjectivity of this mapping follows from the fact every \(\chi \in \text{Hom}(\Lambda,S^1)\) is of the form 66 , for injectivity we check that if \(\chi_{p_1} \equiv \chi_{p_2}\) then for every \(x \in \Lambda\) \[1 = \chi_{p_1}(x) (\chi_{p_2}(x))^{-1}e^{i(p_1 - p_2)\cdot x} \Longrightarrow (p_1 - p_2) \cdot x \in 2\pi \mathbb{Z}\] which means \(p_1 - p_2 \in \Lambda^*\) and hence \([p_1]_\sim = [p_2]_\sim\).

We have thus proven that \(\widehat \Lambda\) can be identified as \(\mathbb{R}^d/\sim\), in particular this group is isomorphic to torus \(\mathbb{T}^d\). For the purpose of this paper we will stick to the definition given in the previous section as it is much easier to work with.

5.4 Distributional Fourier transform↩︎

As we have already identified \(\widehat \Lambda\simeq \mathbb{R}^d/\Lambda^*\), we can use the fact that the latter has a structure of a compact manifold (diffeomorphic to the torus \(\mathbb{T}^d\)) and extend the definition of the Fourier transform on \(\Lambda\) to a distributional one.

Let \(\psi \colon \Lambda \to \mathbb{C}\) be a function with at most polynomial growth, that is \[|\psi(x)| \le C(1 + |x|)^s\] for some constants \(C > 0\) and \(s \in \mathbb{R}\). We define the distributional Fourier transform of \(\psi\), also denoted by \(\widehat \psi\), as a distribution on \(\widehat \Lambda\) given by \[\label{Fourier95def} \langle \widehat \psi, f \rangle = \sum_{x \in \Lambda} \overline{\psi(x)} \check{f}(x),\tag{71}\] where \(f \in \mathcal{D}(\widehat \Lambda) = C^{\infty}(\widehat \Lambda)\) and \(\check{f}\) denotes the inverse Fourier transform 69 . Note that if \(\psi \in L^1(\Lambda)\) then this definition coincides with the standard one 68 .

5.5 Finite Bravais lattice↩︎

For a Bravais lattice \(\Lambda\) as before and for \(L \in 2\mathbb{N}\) we define a finite Bravais lattice \(\Lambda_L\) as \[\begin{align} \Lambda_L &= (A\mathbb{Z}^d)/(LA\mathbb{Z}^d) \\&= \left\{\sum_{j=1}^d m_ja_j \colon \; m_j = -\frac{L}{2}, -\frac{L}{2} + 1, \dots, \frac{L}{2}, \; j=1,\dots,d\right\} \text{ with periodic boundary condition.} \end{align}\] Once again this is an additive group with discrete topology.

Using similar arguments as for the infinite lattice we can show that every \(\chi \in \text{Hom}(\Lambda_L,S^1)\) is of the form \[\label{momentum95representation} \chi(x) = \chi_p(x) = \frac{1}{|\Lambda_L|^{1/2}}e^{i p \cdot x}\tag{72}\] where \(|\Lambda_L| = (L+1)^3\) is the number of points in \(\Lambda_L\) and \(p\) is the element of \[\label{LhL95def} \widehat \Lambda_L := \left\{\sum_{j=1}^d m_j \frac{b_j}{L+1} \colon m_j=-\frac{L}{2},-\frac{L}{2} + 1,\dots, \frac{L}{2} - 1, \frac{L}{2}\right\},\tag{73}\] where \(b_j\) are primitive vectors of the reciprocal lattice \(\Lambda^*\). We will use this identification of \(\widehat \Lambda_L\) for the entire paper.

Not that in 72 we have introduced an additional normalization factor. The reason for it is that when we consider a standard counting measure on \(\Lambda_L\) as its Haar measure then the system \(\{\chi_p\}_{p \in \widehat \Lambda_L}\) forms an orthonormal basis of \(L^2(\Lambda_L)\). We will refer to this system as the momentum basis of \(L^2(\Lambda_L)\).

The correct choice for the Haar measure on \(\widehat \Lambda_L\) is again the standard counting measure on \(\widehat \Lambda_L\). Once again we can write the formulae for the Fourier transform on \(\Lambda_L\) and its inverse on \(\widehat \Lambda_L\) explicitly as \[\label{Fourier95transform95finite} \widehat f(p) = \frac{1}{|\Lambda_L|^{1/2}}\sum_{x \in \Lambda_L} f(x) e^{-i p \cdot x} = \langle\chi_p, f \rangle_{L^2(\Lambda_L)}, \quad p \in \widehat \Lambda_L\tag{74}\] and \[\check g(x) = \frac{1}{|\Lambda_L|^{1/2}} \sum_{p \in \widehat \Lambda_L} g(p) e^{ip \cdot x} = \langle\overline{\chi}_p, g\rangle_{L^2(\widehat \Lambda_L)} \quad x \in \Lambda_L.\] As the Fourier transform is unitary we also note the Parseval identity \[\sum_{x \in \Lambda_L} \sum_{x \in \Lambda_L}\overline{f(x)}g(x) = \sum_{p \in \widehat \Lambda_L} \overline{\widehat f(p)}\widehat g(p).\]

6 Graph calculus↩︎

6.1 Basic definitions↩︎

It might be useful do work within the graph calculus formalism. A graph \(G\) is a couple \(G = (V,E)\), where \(V\) is a finite³ set of vertices and \(E \subset V \times V\) is the set of edges. Note that the edges are directed, that is \((x,y) \ne (y,x)\) for \(x \ne y\). This approach will alow to define directional derivative. We will consider only non-oriented graphs, which in this setting means \[(x,y) \in E \Rightarrow (y,x) \in E.\] We will also assume that there are no self-loops (that is there are no edges of the form \((x,x)\)). We will say that \(x\) and \(y\) are nearest neighbors if \((x,y) \in E\). This defines a symmetric relation on \(V\) that will be denoted as \(x \sim y\). Moreover to each edge \((x,y) \in E\) we will assign a positive real number \(t(x,y)\), which gives rise to the weighted graph structure. Here we will also assume that \(t(x,y) = t(y,x)\) for every edge \((x,y)\).

Consider a subset \(\Omega \subset V\). We define the boundary of \(\Omega\), denoted \(\partial \Omega\), as \[\partial \Omega = \{x \in \Omega \colon \text{there exists } y \not \in \Omega, y\sim x\}\] We also define the set of interior edges \(E_\Omega\) of the set \(\Omega\) as \[E_\Omega = \{(x,y) \in E \colon x,y \in \Omega\}.\] It will be useful to also define the "nearest neighbors boundary" of the set \(\Omega\) defined as \[\partial_\text{nn}\Omega := \{y \not \in \Omega \colon y \sim x \text{ for some }x \in \partial \Omega\}\] and the "nearest neighbors closure" of \(\Omega\) \[\Omega_{\text{nn}} := \Omega \cup \partial_\text{nn}\Omega.\]

With a graph we can associate two Hilbert spaces: the space of functions on the vertices \(L^2(V)\) and functions on the edges \(L^2(E)\) (both with counting measure). For a function \(f: V \to \mathbb{C}\) we define its (discrete) gradient \(\nabla f: E \to \mathbb{C}\) as \[\nabla f(x,y) = \sqrt{t(x,y)}\left(f(y) - f(x)\right).\] For a given edge \((x,y) \in E\) the value \(\nabla f (x,y))\) may be considered as the directional derivative in direction \(x \to y\). As \(\nabla : L^2(V) \to L^2(E)\) we can consider its dual \(\nabla^*:L^2(E) \to L^2(V)\) which satisfies the property that for any \(f \in L^2(V)\) and \(F \in L^2(E)\) we have \[\langle F, \nabla f\rangle_{L^2(E)} = \langle\nabla^*F, f\rangle_{L^2(V)}.\] We can also define \(\nabla^*\) explicitly by the formula \[\label{nabla95star95def} \nabla^* F (x) = \sum_{y \sim x} \sqrt{t(x,y)}\left(F(y,x) - F(x,y)\right).\tag{75}\] Next we can define the discrete divergence \(\text{div}: L^2(E) \to L^2(V)\) as \(\text{div}=-\frac{1}{2} \nabla^*\) and the discrete Laplacian \(\Delta: L^2(V) \to L^2(V)\) as \(\Delta = \text{div}\circ \nabla\). We can check that the action of \(\Delta\) can be written explicitly \[\label{Laplacian95def} \Delta f (x)= \sum_{y \sim x} t(x,y)(f(y) - f(x)) = \sum_{y \sim x} \nabla f(x,y).\tag{76}\] From the definition it is easy to see that the Laplace operator is self-adjoint on \(L^2(V)\). We will show it for completeness: for \(f,g \in L^2(V)\) we have \[\langle f, \Delta g\rangle_{L^2(\Omega)} = \langle g, -\frac{1}{2}\nabla^*\nabla f\rangle_{L^2(E)} = -\frac{1}{2} \langle\nabla g, \nabla f\rangle_{L^2(E)} = \langle-\frac{1}{2} \nabla^* \nabla g, f\rangle_{L^2(V)} = \langle\Delta g, f\rangle_{L^2(V)}\]

We are interested in deriving some properties of the discrete Laplace operator resembling the Green identities that hold for the standard (continuous) Laplacian. Let \(\Omega \subset V\) be a fixed subset. Then we have \[\label{discrete95by95parts} \begin{align} \sum_{x \in \Omega}\overline{g(x)}\Delta f(x) &= \langle\mathop{\mathrm{\mathbb{1}}}_{\Omega}g, \Delta f\rangle_{L^2(V)} = -\frac{1}{2}\langle\mathop{\mathrm{\mathbb{1}}}_\Omega g, \nabla^*\nabla f\rangle_{L^2(V)} = -\frac{1}{2} \langle\nabla \mathop{\mathrm{\mathbb{1}}}_\Omega g, \nabla f\rangle_{L^2(E)} \\& = -\frac{1}{2}\sum_{(x,y) \in E}\nabla (\mathop{\mathrm{\mathbb{1}}}_\Omega \overline{g})(x,y) \cdot \nabla f(x,y) \\& = -\frac{1}{2}\sum_{(x,y) \in E_\Omega} \overline{\nabla g(x,y)} \nabla f(x,y) + \sum_{x \in \partial \Omega}\sum_{\substack{y \not \in \Omega\\ y \sim x}}\overline{g(x)}\nabla f(x,y). \end{align}\tag{77}\] There is no factor \(\frac{1}{2}\) in the second term as cases \(x \in \Omega\), \(y \not \in \Omega\) and \(x \not \in \Omega\), \(y \in \Omega\) are symmetric and give the same contribution. The factor \(\frac{1}{2}\) in the first term is the side effect of considering the ordered pairs in the definition of the edge \((x,y)\), which essentially means every bond between \(x\) and \(y\) is counted twice. This result is the discrete analogue to the standard integration by parts formula \[\int_{\Omega}\overline{g(x)}\Delta f(x) dx = - \int_{\Omega} \overline{\nabla g(x)}\nabla f(x) dx + \int_{\partial \Omega} \overline{g(x)} \frac{\partial f}{\partial n}(x) d\sigma(x).\]

6.2 Neumann Laplacian↩︎

We will define the Neumann Laplacian on some set \(\Omega \subset V\). To this end we first define a quadratic form \[\label{Neumann95form} Q^{\text{Neu}}(f) = \frac{1}{2}\sum_{(x,y) \in E_\Omega} |\nabla f(x,y)|^2.\tag{78}\] We note that the value of \(Q^{\text{Neu}}(f)\) depends only on the restriction of \(f\) to the set \(\Omega\). The Neumann Laplacian \(-\Delta_{\Omega}^{\text{Neu}}\) is defined as the operator on \(L^2(\Omega)\) associated with this quadratic form, meaning that for every \(f \in L^2(\Omega)\) there holds \[\langle f, -\Delta_\Omega^{\text{Neu}}f\rangle_{L^2(\Omega)} = Q(f).\] We can write the action of \(-\Delta_\Omega^{\text{Neu}}\) explicitly: for \(f \in L^2(\Omega)\) we have \[\label{Neumann95Laplacian95def} -\Delta_\Omega^{\text{Neu}} f(x) = \sum_{\substack{y \in \Omega\\y \sim x}} t(x,y)(f(x) - f(y)) = \sum_{\substack{y \in \Omega\\y \sim x}} \nabla f(y,x).\tag{79}\]

Finally we will verify that the operator \(-\Delta_\Omega^{\text{Neu}}\) is self-adjoint, meaning that for every \(f,g \in L^2(\Omega)\) we have \[\langle f, -\Delta_\Omega^{\text{Neu}} g\rangle= \langle-\Delta_\Omega^{\text{Neu}} f, g\rangle.\] This follows from the fact that \(-\Delta_\Omega^{\text{Neu}}\) is the Laplace operator defined as in 76 in the previous subsection for the graph \((\Omega, E_\Omega)\), so self-adjointness follows from the general consideration of graph Laplace operators.

7 The scattering equation on a lattice↩︎

We will start with deriving the formula for the scattering length 12 . To this end we are interested in a solution to the equation (defined on \(\Lambda = A\mathbb{Z}^3\)) \[\label{def:scattering95eq} -\Delta \varphi(x) + \frac{U}{2}\delta_{x,0} \varphi(x) = 0,\tag{81}\] with the condition \[\label{scattering95boundary} \lim_{|x| \to +\infty} \varphi(x) = 1.\tag{82}\] This equation is called the (zero-energy) scattering equation. We will see that this equation has a unique solution, hence it is possible to define the scattering length in a following way.

In order to solve the scattering equation for the moment we will ignore the condition 82 and take the (distributional) Fourier transform (see Appendix 5) of its both sides. A simple computation leads to \[\label{scattering95transform} \mathop{\mathrm{\varepsilon}}(p) \widehat \varphi + \frac{U}{2}\varphi(0) = 0,\tag{83}\] where \(\mathop{\mathrm{\varepsilon}}(p)\) is the dispersion relation, defined in 13 . This equation is satisfied in the sense of distributions, that is after testing against some smooth function on \(\widehat \Lambda\).

For now we will restrict ourselves to the set \(\widehat \Lambda\setminus \{0\}\) and test the above equation with the test function \(\phi\) with \(\mathop{\mathrm{supp}}\phi\) not including zero. On this set \(\left(2\mathop{\mathrm{\varepsilon}}(p)\right)^{-1}\) is a well-defined smooth function and therefore we can multiply both sides of the equation 83 by it. It follows that \[\widehat \varphi = -\frac{U\varphi(0)}{2\mathop{\mathrm{\varepsilon}}(p)}.\] Thus, on this set, we can identify \(\widehat \varphi\) as a \(L^1(\widehat \Lambda)\) function (note that this function would not be integrable in the dimensions \(d=1\) and \(d=2\)).

By restricting our considerations to the set not containing zero, we might have neglected distributions whose support is the one-point set \(\{0\}\). Since distributions supported on one point are the sums of Dirac deltas and their derivatives, we conclude that \[\widehat \varphi = -\frac{U\varphi(0)}{2\mathop{\mathrm{\varepsilon}}(p)} + \sum_{\alpha \colon |\alpha| \le M} c_\alpha \partial^\alpha \delta_0,\] for some \(M \ge 0\) and \(c_\alpha \in \mathbb{C}\). By equation 83 we need to have \[\mathop{\mathrm{\varepsilon}}(p) \cdot \left(\sum_{\alpha \colon |\alpha| \le M} c_\alpha \partial^\alpha \delta_0\right) = 0.\] It follows that \(M = 1\) as the value of function \(\mathop{\mathrm{\varepsilon}}(p)\) and all of its first order derivatives are zero at \(p=0\), whereas values of second order derivatives at \(p=0\) are non-zero. A consequence of this observation is that \[\widehat \varphi = -\frac{U\varphi(0)}{2\mathop{\mathrm{\varepsilon}}(p)} + C_0 \delta_0 + \sum_{j=1}^3C_j \partial_{p_j}\delta_0.\] Using the inverse Fourier transform (see equation 69 in the Appendix) we get \[\varphi(x) = -\frac{U\varphi(0)}{2} |\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{e^{ip \cdot x}}{\mathop{\mathrm{\varepsilon}}(p)}dp + C_0 + \sum_{j=1}^3 C_j x_j.\] The value \(\varphi(0)\) is not yet specified, we need to make sure that this function is self consistent with its value at \(x=0\). Before that we will simplify this expression by using the boundary condition 82 that so far we have omitted. By the Riemann-Lebesgue lemma we have \[\lim_{|x| \to \infty} |\widehat \Lambda|^{-1} \int_{\widehat \Lambda} \frac{e^{ip \cdot x}}{\mathop{\mathrm{\varepsilon}}(p)}dp = 0,\] so this part of the scattering equation solution vanishes. An easy observation also leads to conclusion that in order to satisfy 82 we need to have \(C_0 = 1\) and \(C_j = 0\) for \(j=1,2,3\). We have thus simplified the formula for \(\varphi\) to \[\varphi(x) = 1 - \frac{U\varphi(0)}{2}|\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{e^{ip \cdot x}}{\mathop{\mathrm{\varepsilon}}(p)}dp.\] Computing the value at \(x=0\) we have \[\varphi(0) = 1 - \frac{U\varphi(0)}{2}|\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{1}{\mathop{\mathrm{\varepsilon}}(p)}dp = 1 - U\varphi(0)\gamma,\] where \[\gamma = \frac{1}{2}|\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{1}{\mathop{\mathrm{\varepsilon}}(p)}dp.\] This leads to \[\label{zero95value} \varphi(0) = \frac{1}{1 + U\gamma}\tag{84}\] and \[\label{def:scattering95sol} \varphi(x) = 1 - \frac{1}{2} \cdot \frac{U}{1 + U\gamma} |\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{e^{ip \cdot x}}{\mathop{\mathrm{\varepsilon}}(p)}dp.\tag{85}\] Using 84 in the Definition 11 we can explicitly write \[\label{scattering95sol95U} 8\pi \boldsymbol{a}= \frac{U}{U\gamma + 1},\tag{86}\] which is the definition used in 12

It will also be useful to introduce function \(w(x) := 1 -\varphi(x)\) or explicitly \[w(x) = \frac{1}{2} \cdot \frac{U}{1 + U\gamma} |\widehat \Lambda|^{-1}\int_{\widehat \Lambda} \frac{e^{ip \cdot x}}{\mathop{\mathrm{\varepsilon}}(p)}dp.\] This function satisfies the equation \[\Delta w(x) + \frac{1}{2}U(1-w(x))\delta_{x,0} = 0\] with a condition \[\lim_{|x| \to \infty} w(x) = 0.\] The main advantage of considering this function instead of \(\varphi(x)\) is that its (once again distributional⁴) Fourier transform \(\widehat w(p)\) can be treated as a \(L^1(\widehat \Lambda)\) function (and not only as a distribution): \[\label{w95transform} \widehat w(p) = \frac{U(1-w(0))}{2\mathop{\mathrm{\varepsilon}}(p)} = \frac{U}{1 + U\gamma} \cdot \frac{1}{2\mathop{\mathrm{\varepsilon}}(p)}, \quad p \ne 0.\tag{87}\] We will also note two useful equalities that are frequently used in various parts of the paper: \[\label{useful} w(0) = \frac{U\gamma}{1 + U\gamma}, \quad 1 - w(0) = \frac{1}{1 + U\gamma}.\tag{88}\]

References↩︎

---

1. Note that the definition of \(L(k)\) is correct as the number of points in \(\Lambda_L\) is \(|\Lambda_L| =(L+1)^3\), hence \(|\Lambda_{L(k)}| = k^3(L+1)^3\).↩︎

2. Note that this is the moment where it is important that we do not consider only symmetric wave functions.↩︎

3. We can also consider infinite, but countable sets of vertices. This however requires adding some technical assumptions on summability of functions on vertices and edges.↩︎

4. This function asymptotically behaves as \(\frac{1}{|x|}\) for large \(|x|\), hence it is not summable. This is the reason why we cannot use the standard Fourier transform.↩︎