Observable signature of the Berezinskii–Kosterlitz–Thouless transition in a planar lattice of Bose–Einstein condensates

We investigate the possibility that Bose–Einstein condensates, loaded on a 2D optical lattice, undergo—at finite temperature—a Berezinskii–Kosterlitz–Thouless transition. We show that—in an experimentally attainable range of parameters—a planar lattice of Bose–Einstein condensates is described by the XY model at finite temperature. We demonstrate that the interference pattern of the expanding condensates provides the experimental signature of the Berezinskii–Kosterlitz–Thouless transition by showing that, near the critical temperature, the component of the momentum distribution and the central peak of the atomic density profile sharply decrease. The finite-temperature transition for a 3D optical lattice is also discussed in this paper, and analogies with superconducting Josephson junction networks are stressed throughout the text.


Introduction
Recent advances in the manipulation of cold atomic gases allow nowadays for storing Bose-Einstein condensates (BECs) in 1D [1]- [7], 2D [8] and 3D [9] optical lattices, opening a pathway for the experimental and theoretical analysis of finite-temperature transitions in atomic systems. In a way similar to the spectacular realizations obtained with superconducting Josephson networks [10,11], arrays of BECs can provide new experimentally realizable systems on which to test well-known paradigms of the statistical mechanics: our interest here is to find out how to detect a finite-temperature defect-mediated transition in 2D optical lattice of BECs.
According to the well-known Mermin-Wagner theorem [12], 2D systems with a continuous symmetry cannot sustain long-range order in the thermodynamic limit at finite temperature; however, a phase transition is still possible and it occurs via the unbinding of point defects like dislocations or vortices. Defect-mediated transitions are widely studied [12] since they provide crucial insights for a wide variety of experiments on thin films. The Berezinskii-Kosterlitz-Thouless (BKT) transition is the paradigmatic example of a defect-mediated transition [12] which is exhibited by the 2D XY model [13,14], describing two-component spins on a 2D lattice. In the low-temperature phase, characterized by the presence of bound vortex-antivortex pairs, the spatial correlations exhibit a power-law decay; above a critical temperature T BKT , the decay is exponential and there is a proliferation of free vortices. The BKT transition has been observed in superconducting Josephson arrays [11] and its predictions are well verified in the measurements of the superfluid density in 4 He films [15].
In finite magnetic systems, the BKT transition point is signalled by the decrease to 0 of a suitably defined magnetization [16]; the existence of such a magnetization in finite systems does not contradict the Mermin-Wagner theorem, since the latter is valid only in the thermodynamic limit. In a large class of experimental situations, the finite size XY model predicts magnetization exponents in agreement with the measurements carried out for layered magnets with planar spin symmetry [17]. We shall show in the following that the Fourier transform at k = 0 of the condensate wavefunction provides the analogue of this magnetization for atomic systems, which are naturally finite (since they are trapped by a harmonic potential).
In this paper, we show that-at a critical temperature T BKT lower than the temperature T BEC at which condensation in each well occurs-2D lattices of BECs may undergo a phase transition to a superfluid regime where the phases of the single-well condensates are coherently aligned allowing for the observation of a BKT transition. We point out, in the following, that the recently realized 2D optical lattice of BECs [8] may allow for the observation of a BKT transition in a bosonic system. Indeed, the thermodynamic properties of the bosonic lattice at finite temperature may be well described, over a suitable and experimentally attainable range of parameters, by the Hamiltonian of the XY model. It is easy to convince oneself that the transition we propose to observe is different from the quantum-phase transition reported in [9], where the system is at T = 0 and the insulator phase (signalled by a reversible destruction of the phase coherence across the lattice) is reached by varying the optical lattice parameters: at variance, here, we propose to fix the system parameters in the superfluid phase and increase the temperature T until the thermally induced vortex proliferation determines the BKT transition.
In addition, we shall show that the experimental signature of the BKT transition in bosonic planar lattices is obtained by measuring, as a function of the temperature, the central peak of the interference pattern of the expanding condensates released from the trap when the confining potentials are switched off. In fact, the peak of the momentum distribution at k = 0 may be regarded as the magnetization of a finite size 2D XY magnet and, as shown in the following, it should exhibit a sharp decrease around T BKT . As a consequence, also the central peak of the atomic density profile decreases around T BKT , signalling the occurrence of the BKT transition.
The plan of the paper is as follows: in section 2, we discuss the Hamiltonian describing a 2D lattice of BECs at finite temperature and we show that-in a suitable range of parameters-the system is described by the XY Hamiltonian. In section 3, we show how the signature of the BKT transition may be evidenced in a bosonic planar lattice and we shall also discuss interesting analogies with superconducting Josephson junction networks. Section 4 is devoted to a discussion of the finite-temperature transition for a 3D optical lattice, and to our concluding remarks. In appendix A, we provide the details of a variational estimate of the coefficients of the Hamiltonian describing 1D, 2D and 3D optical lattices, while appendix B contains the analytical computation of the Fourier transform of a vortex on a lattice.

The finite-temperature Hamiltonian for the system
For 2D optical lattices, when the polarization vectors of the two standing wave laser fields are orthogonal, the resulting periodic potential for the atoms is where k = 2π/λ is the wavevector of the laser beams. The potential maximum of a single standing wave, V 0 = sE R , may be controlled by changing the intensity of the optical lattice beams, and it is conveniently measured in terms of the recoil energy E R =h 2 k 2 /2m (m is the atomic mass), while, typically, s can vary from 0 up to 40 (we mention that very recent experiments with 1D optical lattices reported much larger lattice depths [7]). Gravity is assumed to act along the z-axis.
Usually, superimposed to the optical potential, there is a harmonic magnetic potential where ω z (ω r ) is the axial (radial) trapping frequency. The minima of the 2D periodic potential (1) are located at the points j = (j 1 , j 2 ) · λ 2 with j 1 and j 2 being integers, and the potential around provides the frequency of the wells. Whenω r ω r , ω z , the system realizes a square array of tubes, i.e. an array of harmonic traps elongated along the z-axis [8]. Although the axial dynamics is very interesting in its own right [18], in the following we analyse only the situation in which the axial degrees of freedom are frozen out, which is realizable if ω z is sufficiently large. We also assume that the harmonic oscillator length of the magnetic potential √h /mω r is larger than the size L (in lattice units) of the optical lattice in order to decrease density-inhomogeneity effects and to allow for a safe control of finite-size effects. In addition, when all the relevant energy scales are small compared to the excitation energies, one can safely expand the field operator [19] asˆ with j ( r ) being the Wannier wavefunction localized at the jth well (normalized to 1) and N j =ψ † jψ j the bosonic number operator. Substituting the expansion in the full quantum Hamiltonian describing the bosonic system leads to the Bose-Hubbard model [20,19] where i, j denotes a sum over nearest neighbours, U = (4πh 2 a/m) d r 4 j (a is the s-wave scattering length) and (where N 0 is the average number of atoms per site), when N 0 1 and J/N 2 0 U, the Bose-Hubbard model reduces tô which describes the so-called quantum phase model (see e.g. [21,22]): in equation (7), θ j is the phase of the jth condensate and stands for the Hamiltonian of the classical XY model. When J U and at temperatures T U/k B , the pertinent partition function describing the thermodynamic behaviour of the BECs stored in an optical lattice may be computed using the classical XY model (8). The experimental values are taken, respectively, from [2], [8], and [9], with an average number of 87 Rb atoms N 0 = 1000 (1D), 170 (2D) and 1 (3D).
The close relationship between the Bose-Hubbard model (5) and the quantum phase model (7) is, by now, well known [23]- [25] and it has been widely used in the analysis of superconducting Josephson networks; since, in superconducting networks, N 0 is the average number of excess Cooper pairs at a given site, the requirement N 0 1 is of course satisfied and the condition J/N 2 0 U may also be easily realized for reasonable values of J. At variance, in bosonic lattices, N 0 varies usually between ∼1 and ∼10 3 and the validity of the mapping between the Bose-Hubbard model and the quantum-phase Hamiltonian is not always guaranteed; furthermore, in order to get the XY model (8) from the Bose-Hubbard model (5), U should be much smaller than J, but not vanishing, in order to satisfy to the condition J/N 2 0 U. Of course, the Bose-Hubbard model-for U = 0-describes harmonic oscillators and thus cannot sustain any BKT transition; nevertheless, when N 0 1 and J/N 2 0 U J, the Bose-Hubbard Hamiltonian reduces in the XY model (8), which does display the BKT transition.
A simple estimate of the coefficients J and U may be obtained by approximating the Wannier functions in equation (4) with Gaussians, whose widths are determined variationally. The details of the computation are reported in appendixA (see also [26]). In figure 1, we plot J/U as a function of V 0 for a 2D optical lattice and values for experimental parameters form [8]; for comparison, we plot the same ratio for 1D and 3D optical lattices also. Since, in a mean-field approach, the Mott insulator-superfluid transition occurs, for T = 0, approximately at 2zJ/U ∼ 1 [20,22] (where z is the number of nearest neighbours), from figure 1 one sees that it is much easier to detect a quantum-phase transition for the 3D array (and indeed it has been detected in [9]); no quantum-phase transition has yet been reported in the literature for 2D, and recently, signatures for the quantum-phase transition have been observed for 1D optical lattices [6,7].
For realistic values of the experimental parameters, the mapping between the Bose-Hubbard model (5) and the XY model (8) may be regarded as very accurate. In fact, if one takes, as in [8], λ = 852 nm (so that E R = 2.08 × 10 −23 erg), the number of sites is ∼2000 and the average number of particles per site N 0 ∼ 170, and if one chooses s = V 0 /E R between 18 and 25, one sees that the XY model well describes the properties of bosons in a 2D lattice (for definiteness, we take s = 20). Furthermore, if one assumes a magnetic confinement along z with frequency ω z /2π 500 Hz and a magnetic confinement in the x-y plane with frequency ω r /2π ∼1-5 Hz, one sees that the frequencyω r in each well defined in equation (3) is given byω r = 2π √ s · 6.3 kHz = 2π · 28.2 kHz. From the variational estimate (see appendix A and figure 1), one finally gets K/E R ≈ 0.0004, U/E R ≈ 0.018 and J/E R = 2KN 0 ≈ 0.14, and thus the conditions J U J/N 2 0 are rather well satisfied; the BKT critical temperature T BKT ∼ J/k B is, then, between 10 and 30 nK. A comparison between the thermal energy and the energy of the excitations in the z-direction provides that the condition k B T hω z is needed in order to have frozen axial (i.e., along z) excitations: it turns out thathω z /k B ≈ 25 nK so that, below T BKT , the degrees of freedom along z are not excited. Furthermore, as we shall show later, the BEC temperature T BEC of a single well is larger than 500 nK and thus T BKT T BEC . The size of the system is of the order of 25λ/2 ≈ 10 µm: this length should be compared with √h /mω r , which is ≈ 7µm for ω r = 2π · 2.5 Hz (i.e., the effect of inhomogeneity of the confining magnetic trap in the x-y plane is rather negligible).
Since the BKT transition for the XY model occurs at a temperature T BKT ∼ J/k B , it should be possible to observe the transition with measurements performed at different temperatures. Using a coarse-graining approach to determine the finite-temperature phase-boundary line of the Bose-Hubbard model [27] yields, in 2D and for N 0 1 and J/U 1, a critical temperature T BKT ≈ J/k B , in reasonable agreement with the XY estimates of T BKT ; in addition, in the chosen range of parameters, the condition UN 2 0 /J = UN 0 /K 1 satisfies the finite-temperature stability criterion recently derived by Tsuchiya and Griffin [28].
We notice that the XY Hamiltonian, in the limits N 0 1 and J U J/N 2 0 , can also be retrieved by extending the familiar Gross-Pitaevskii Hamiltonian [29] to finite temperatures T U/k B . Indeed, one can use, for the BEC, the tight-binding ansatz ( r, t) = j ψ j (t) j ( r ), where is the condensate wavefunction of the BEC in the whole optical lattice and ψ j is the wavefunction of the condensate in the jth well (the ψ j values are now classical fields and not operators). Substituting the tight-binding ansatz in the Gross-Pitaevskii Hamiltonian one gets the lattice Hamiltonian which is the classical version of the Bose-Hubbard model and it is also known as the discrete nonlinear Schrödinger Hamiltonian (we use the notation N j ≡ |ψ j | 2 ). Writing N j = N 0 + δN j , one may neglect the quadratic terms [i.e., (δN j ) 2 ] in the hopping part of the Hamiltonian (10), which then reduces to H XY (8). One may then conclude that, within the range of our approximations, the discrete nonlinear Schrödinger Hamiltonian (10) and the XY Hamiltonian (8) provide equivalent descriptions of the system at finite temperature.
Accurate Monte Carlo simulations yield-for the XY model-the BKT critical temperature T BKT = 0.898J/k B [30]. When U J, a BKT transition still occurs at a slightly lower critical temperature T BKT (U ). Intuitively speaking, when U increases, the superfluid region in the phase diagram decreases and, thus, one has T BKT (U ) < T BKT [22]. An estimate of the dependence of T BKT (U ) on U obtained using a renormalization group approach has been reported in [31]: when U = 0, but still U/J 1, the only effect of the interacting term amounts to renormalizing J and lowering the BKT critical temperature [32]. For this reason, in the following, we shall present results from Monte Carlo simulations of the classical XY model only.
In numerical simulations of the finite XY lattice, one defines the critical temperature as the temperature T C (L) at which the correlation length equals the size L of the square lattice [17]; of course, as L → ∞, T C (L) → T BKT . For T = T C (L), the BKT transition is signalled in a finite system by the fact that a suitably defined magnetization, defined by equation (12), decreases to zero [17].
The emerging physical picture is then the following: there are two relevant temperatures for the system, the temperature T BEC at which in each well there is a condensate, and the temperature T BKT at which the condensate phases begin to coherently point in the same direction. Of course, in order to have well-defined condensate phases, one should have T BKT T BEC . The critical temperature T BEC is given by T BEC ≈ 0.94hN [29]. With the numerical values given in [8], T BEC turns out to be 500 nK for s 20. When T < T BEC , the atoms in the well j of the 2D optical lattice may be described by a macroscopic wavefunction ψ j . Furthermore, when the fluctuations around the average number of atoms per site N 0 1 are strongly suppressed, one may put, apart from the factor √ N 0 constant across the array, ψ j = e iθ j . The temperature T BKT is of the order of J/k B : with the experimental parameters of [8] and V 0 between 20 and 25E R , one has that T BKT is between 10 and 30 nK, which is sensibly smaller than the condensation temperature T BEC of the single well.
A similar scenario describes also the phases of planar arrays of superconducting Josephson junctions [21,22]: they exhibit a temperature T BCS at which the metallic grains placed on the sites of the array become (separately) superconducting and the Cooper pairs may be described by macroscopic wavefunctions. At a temperature T BKT < T BCS , the array undergoes a BKT transition and the system, as a whole, becomes superconducting.

Observable signature of the BKT transition
In this section, we shall provide evidence that the experimental signature of the BKT transition in bosonic planar lattices is obtained by measuring, as a function of the temperature, the central peak of the interference pattern obtained after turning off the confining potentials. We notice that the experimental signature for the BKT transition for a continuous (i.e., without optical lattice) weakly interacting 2D Bose gas [33] is also given by the central peak of the atomic density of the expanding condensates.
Firstly, we show that the peak of the momentum distribution at k = 0 is the direct analogue of the magnetization of a finite-size 2D XY magnet. In fact, for the XY magnets, the spins can be written as S j = (cos θ j , sin θ j ) and the magnetization is defined as M = (1/N) | j S j | , where · stands for the thermal average; a spin-wave analysis at low temperatures yields M = (2N) −k B T/8πJ [16,17]. With discrete BECs at T = 0, all the phases θ j are equal at the exhibits a peak at k = 0 ( k is in the first Brillouin zone of the 2D square lattice); the magnetization is then The intuitive picture of the BKT transition is the following: at T = 0, all the spins point in the same direction. As one can see from figure 2(A), a single free vortex modifies the phase · · · · · ·, fit near T C ≈ 1.07J/k B as in [17] (in this figure L = 40).
distribution far away from the vortex core and, thus, the modulus square of its lattice Fourier transform |ψ k | 2 has a minimum at k = 0. At variance, a vortex-antivortex pair (see figure 2(B)) modifies the phase distribution only near the centre of the pair (in this sense is a local defect) and its lattice Fourier transform has a maximum at k = 0. Analytical expressions for the lattice Fourier transformψ k of a single vortex and a vortex-antivortex pair may be worked out in detail and are reported in appendix B. Upon increasing the temperature, vortices are thermally induced. For T < T BKT , only bound vortex pairs are present and, on the average, the spins continue to point in the same direction. When the condensates expand, a large peak (i.e., a magnetization) is observed in the central k = 0 momentum component, as is shown in figure 2(C). Increasing the temperature further, due to the increase in the number of vortex pairs, the central peak density decreases. For T ≈ T BKT , the pairs start to unbind and free vortices begin to appear (see figure 2(D)), determining a sharp decrease of the magnetization around T BKT . At high temperatures, only free vortices are present, leading to a randomization of the phases and to a vanishing magnetization. In figure 3, we plot the intensity of the central peak of the momentum distribution (normalized to the value at T = 0) in a 2D lattice as a function of the reduced temperature k B T/J , evidencing the sharp decrease of the magnetization around the BKT critical temperature.
Let us now turn our attention to the interference patterns of the expanding condensates: after the bosonic system reaches the equilibrium, one may switch off both the harmonic trap and the optical lattice. At this time (t = 0) the momentum distributionψ( p, t = 0), is given by the Fourier transform of ψ( r, t = 0) = j ψ j j ( r ): one finds where˜ ( p) is the Fourier transform of the 3D Wannier functions andh k = (p x , p y ) is the momentum projection on the first Brillouin zone of the 2D optical lattice. Using a Gaussian for the Wannier functions, ( r ) x /2h 2 and, similarly, for χ(p y ) and χ(p z ). Since, for s 1, σ/d = 1/πs 1/4 , the momentum distribution exhibits wellpronounced peaks at the centres of each Brillouin zone: these peaks have different heights and the central k = 0 peak is the largest (see figures 2 and 3 of [8]).
At t = 0, the amplitude of the k = 0 peak of the momentum distribution is simply given by the thermal average ofψ 0 . By measuring the k = 0 peak (i.e., |ψ 0 | 2 ) at different temperatures, one obtains the results plotted in figure 2. The figure has been obtained using a Monte Carlo simulation of the XY magnet for a 40 × 40 array: we find k B T C ≈ 1.07J. In figure 2, we also plot the low-temperature spin wave prediction [16] (solid line), as well as a fit, first used in [17], valid near T C and derived from the renormalization group equations (dotted line). At times different from t = 0, the density profiles are well reproduced by the free expansion of the ideal gas: one obtainsψ( p, t) = χ(p z )e −i[(p z +mgt) 3 −p 3 z ]/6m 2 ghφ (p x , p y , t), i.e., a uniformly accelerating motion along z and a free motion in the x-y plane, withφ(p x , p y , t) = χ(p x )χ(p y )ψ k e −i(p 2 x +p 2 y )t/2hm , giving the central and lateral peaks of the momentum distribution as a function of time for different temperatures.
An intense experimental work is now focusing on the BEC in two dimensions: at present, a crossover to two-dimensions behaviour has been observed for Na [34] and Cs atoms [35]. Our analysis relies on two basic assumptions; namely, the validity of the tight-binding approximation for the Bose-Hubbard Hamiltonian, and the requirement that the condensate in the optical lattice may be regarded as planar [36]. It is easy to see that the first assumption is satisfied if, at T = 0, it satisfies the conditions V 0 µ (where µ is the chemical potential) and, at finite temperature,hω r k B T (so that the thermal energy does not excite higher bands). In [8], it is V 0 = s · k B · 0.15 µK and, for s 18, one hashω r 1 µK, so that the conditionhω r k B T is well satisfied in experiments. The condition V 0 µ is also satisfied with s 18. The planarity requirement is much more restrictive since it amounts to freezing the axial excitations; for this to happen, one should constrain the axial trapping frequency ω z . First of all, one should have ω z ω r : assuming ω z /2π ∼ 500 Hz (or more) and ω r /2π ∼1-10 Hz, this condition is satisfied. Furthermore, at T = 0, it should be alsohω z K (so that the energy associated with the hopping of particles between neighbouring sites does not excite the condensates along z) and, at finite temperature,hω z k B T (so that thermal energy does not excite the axial degrees of freedom).
Since ω z ω r , ifhω z k B T is satisfied one has also thathω r k B T . For s 18, one has K/ h 2 Hz; thus,hω z K. Furthermore, if ω z ∼ 2π × 500 Hz andhω z k B T , one may safely regard our finite temperature analysis to be valid at least up to T ∼ 25 nK. It is clear that what is more demanding is to have ω z /2π of order of 500 Hz or possibly more, and this makes it challenging to have a true planar system. Notice that, with ω z /2π of order of several kHz, one could choose s ≈ 15 in order to have higher critical temperatures T BKT (100 nK or more).

Concluding remarks
We analysed the finite-temperature phase transitions in 2D lattices of BECs. Our study provides evidence for the possibility that BECs loaded on a 2D optical lattice may exhibit-at finite temperature-a new coherent behaviour in which all the phases of the condensates, located in each well of the lattice, point in the same direction. The finite-temperature transition, which is due to the thermal atoms in each well, is mediated by vortex defects and may be experimentally detectable by looking at the interference patterns of the expanding condensates. Our analysis relies on the approximation that the effect of the shallow confining harmonic trap in the x-y plane may be neglected: for a tight confinement, we expect that interesting mesoscopic effects will come into play, affecting the BKT transition. We observe that for a 3D optical lattice (resulting from the potential V opt ( r ) = V 0 [sin 2 (kx) + sin 2 (ky) + sin 2 (kz)]) with N 0 ∼ 100, the system maps on the 3D XY model, which is the prototype of the universality class including the λ normal-superfluid transition in 4 He [37]. Also, in 3D, the central peak of the interference pattern tends to zero at the critical temperature (see figure 4), but the 3D transition is a simple order-disorder transition (whose critical temperature T c is given in the thermodynamic limit by 2.202J/k B [38]). A natural question arising from the comparison of figures 3 and 4 is how to clearly characterize the two transitions: to answer, one may observe that the critical exponent β (defined by M ∝ (T c − T ) β ) is β ≈ 0.23 for the 2D XY universality class [17] and β ≈ 0.35 for the 3D XY universality class [39].
In conclusion, our analysis strengthens, and extends at finite temperature, the striking and deep analogy of bosonic systems with superconducting Josephson junction arrays [1].

Substituting (A.3) in (A.4), one gets the following approximate expression for the energy in the
The optimal values of σ and σ z are determined requiring that ∂E i /∂σ = 0 and ∂E i /∂σ z = 0. It should be stressed that σ being the variational width in the directions x and y, it should be less than λ/2, the distance between minima of the periodic potential. We notice also that, for realistic values of the parameters, the variational width σ in the x and y directions has a very weak dependence on the mean-field term (and, therefore, on the number of particles at each site): in fact, putting σ ∼ 0.2 µm and σ z ∼ 5 µm, one findsh 2 /2mσ 2 ∼ 10 −23 erg g 0 N 0 /2(2π) 3/2 σ 2 σ z ∼ 10 −29 erg. Thus, one sees that the last term in (A.5) contributes very little to the determination of σ. Since V 0 k 2 mω 2 r /2, from (A.5), one finds that the dependence of E i on σ is given by and U E R = 2mg 0 (2π) 7/2h 2 σ z χ 2 . (A.8) The variational procedure used in this appendix is rather accurate also for 1D and 3D optical lattices: in a 1D optical lattice in the x direction, one has . Typical values of the parameters relevant for experiments [2] are ω x ≈ 2π · 10 Hz, ω ⊥ ≈ 2π · 100 Hz, λ = 795 nm, for a total number of sites ∼100 and an average number of particles per site ∼1000. The variational ansatz is still given by (A.3), but now σ y = σ z ≡ σ ⊥ and σ x ≡ σ is the variational width in the laser direction. Proceeding as before, one obtains that σ = λ χ, with χ still given by (A.6). The tunnelling rate K is given by and the nonlinear coefficient U reads U/E R = 2mλg 0 /(2π) 7/2h2 σ 2 ⊥ χ, where χ is given by equation (A.6).

Appendix B. Fourier transform of a vortex on a lattice
In this appendix, we report the analytical computation of the lattice Fourier transformψ k of a single vortex on a lattice described by ψ j = exp (iθ j ); namely, one should computẽ In (B.1), j = (j x , j y ) denotes the sites of a square lattice (having N sites) and k = (k x , k y ) the (quasi)momentum, with k x and k y valued between −π and π. The relationship between the lattice Fourier transform (B.1) and the momentum distributionψ( p) in real space is provided by (13).
To simplify matters, in this appendix, we set the origin O of the coordinates in the centre of the central plaquette: setting the lattice length to 1, the lattice sites are labelled by j = (j 1 /2, j 2 /2) with j 1 and j 2 positive and negative odd integer numbers.