document: the mott - hubbard model of interacting bosons on a lattice has been used to describe superfluid mott - insulator transitions in a variety of systems , e.g. , josephson arrays and granular superconductors @xcite . 
   the recent suggestion @xcite to experimentally observe this transition in a system of cold bosonic atoms in an optical lattice and its successful experimental demonstration @xcite has rekindled the interest in the mott - insulator transition and triggered a great deal of theoretical @xcite and experimental @xcite activity . the possibility to directly manipulate and test the many - body behavior of a system of 
   trapped bosonic atoms in an optical lattice @xcite is very attractive . 
   possible applications include the use of a mott state of bosonic atoms in an optical lattice as a starting point to create controlled multiparticle entanglement as an essential ingredient for quantum computation @xcite    the mott - insulator quantum phase transition is driven by the interplay of the repulsive interaction of bosons on the same lattice site and the kinetic energy . 
   hence the ratio of the onsite energy and the bandwidth forms the key parameter in the system . in optical lattices , this parameter can be easily controlled and varied by several orders of magnitude , enabling detailed studies of the quantum phase transition . 
   probing the system by taking absorption pictures to image the expansion patterns after a reasonable expansion time yields information about the momentum distribution of the state . 
   this procedure was used to experimentally confirm the mott transition in an optical lattice @xcite .    the essential physics of cold bosonic atoms in an optical lattice is captured by a bosonic mott - hubbard model describing the competition between hopping and on - site interaction . 
   a number of approximation schemes have been used to study this model analytically @xcite as well as numerically , using approaches like the gutzwiller mean - field ansatz @xcite , density - matrix renormalization group ( dmrg ) @xcite , exact diagonalization ( ed)@xcite and quantum monte carlo ( qmc ) @xcite .    in this article , we study the short - range correlations , not included by the gutzwiller ansatz , by using perturbation theory . the main purpose is to find corrections to the short - range behavior of the one - particle density matrix , which is directly relevant to experimentally observed expansion patterns . 
   these patterns are important for determining the location of the insulator - superfluid transition . 
   we note that in the insulating state our perturbative approach is identical to the one used in @xcite ( see also @xcite ) , although there the goal was different , viz . , studying corrections to the phase diagram . 
   the remainder of the article is organized as follows : in section [ modsec ] , we will introduce the model and its mean - field solution . the general perturbative approach is briefly outlined in section [ secpa ] , while details may be found in the appendix . 
   numerical results are presented and discussed in section [ secnr ] , first for local observables ( [ secslso ] ) and then for the density matrix ( [ secrho ] ) . 
   implications for expansion patterns both for bulk systems and a harmonic confining potential are discussed in section [ harmotrap ] . 
   the cold bosonic gas in the optical lattice can be described by a mott - hubbard model @xcite @xmath0 here , @xmath1 is the total number of lattice sites , @xmath2 ( @xmath3 ) creates ( annihilates ) a boson on site @xmath4 , @xmath5 , @xmath6 is the on - site repulsion describing the interaction between bosons on the same lattice site , and @xmath7 denotes the chemical potential . 
   the kinetic term includes only hopping between nearest - neighbor sites , this is denoted by the summation index @xmath8 ; @xmath9 is the hopping matrix element that we will assume to be lattice - site independent . 
   finally , @xmath10 describes an external on - site potential that is commonly present in experiments . 
   the gutzwiller ( gw ) approach is based on an ansatz for the many - body ground state that factorizes into single lattice - site wavefunctions @xmath11    the gutzwiller wavefunction represents the ground state of the following mean - field version of the mott - hubbard hamiltonian , eq . 
   ( [ motthubb ] ) : @xmath12 here @xmath13 is the mean - field potential on the @xmath4-th lattice site , which is self - consistently defined as the expectation value of @xmath3 in terms of the gutzwiller wavefunction , @xmath14 @xcite . 
   using the gutzwiller ansatz to obtain an approximate variational solution for the mott - hubbard hamiltonian ( [ motthubb ] ) corresponds , however , to restricting the hilbert space to the subset of product states . 
   consequently , even in higher dimensions , this ansatz fails to describe the correct behavior of short - range correlations between different lattice sites , which are important for experimentally measurable observables , such as expansion patterns ( momentum distributions ) . nevertheless , in the thermodynamic limit and higher dimensions , the gutzwiller wavefunction provides a good approximation in the limits of @xmath15 and @xmath16 ( i.e. , deep in the mott insulator ( mi ) and superfluid ( sf ) phases ) . 
   to get a satisfactory description of the short - range correlations we will now derive perturbative corrections to the gutzwiller mean - field result . 
   our aim is to start from the gutzwiller approximation and improve it by perturbatively including the short - range correlations between lattice sites . 
   we re - express the mott - hubbard hamiltonian ( [ motthubb ] ) by adding the appropriate perturbation to the mean - field hamiltonian , eq . 
   ( [ mf ] ) : @xmath17 with @xmath18 as the mean - field hamiltonian represents a sum of single lattice - site hamiltonians , the excited states @xmath19 and the excitation spectrum @xmath20 can be obtained numerically for each lattice site @xmath4 separately . hence we can write the excitations of @xmath21 as product states of single lattice - site excitations , @xmath22 and the excitation spectrum as a sum over the single lattice - site excitation energies , @xmath23 where @xmath24 describes the set of quantum numbers characterizing the given many - body energy eigenstate .    having obtained the mean - field solution from the gutzwiller ansatz , we can now proceed to improve our wavefunction performing rayleigh - schrdinger perturbation theory @xcite in @xmath25 : @xmath26 where @xmath27 and @xmath28 are the @xmath29-th order corrections to the wavefunction and grand - canonical energy @xmath30 , respectively .    knowing the excited wavefunctions , eq . 
   ( [ es ] ) , and the excitation spectrum , eq . 
   ( [ esp ] ) , perturbative corrections to observables can be calculated explicitly . 
   the resulting expressions up to second order in @xmath25 , which we use in the following , are derived and illustrated diagrammatically in appendix [ appendixa ] . 
   we start by computing the gutzwiller wavefunction numerically using a conjugate - gradient descent method . 
   propagation steps of the time - dependent gutzwiller equations @xcite in imaginary time are performed to test whether the minimum found by the conjugate - gradient descent method is indeed the ground state . 
   afterwards , we calculate eigenenergies and -vectors of each site , with respect to the mean - field hamiltonian . as explained above , this forms the basis for perturbation theory , which we use to calculate the corrections to the observables up to second order . 
   the results obtained from perturbation theory show important modifications to single - site observables as well as to the correlation function . 
   these observables are composed of operators acting only on single lattice sites . 
   hence , they describe local properties and will be less sensitive to the correlations between lattice sites . 
   thus , values for single - site observables obtained from the gutzwiller wavefunction already provide a good approximation in most cases ( unless fluctuations are concerned ) . 
   ( diamonds ) , @xmath31 ( triangles ) , and @xmath32 ( squares ) dimensions . 
   the dashed - dotted lines in fig . 
   [ slso]a - c are the gutzwiller results . 
   ( a ) number fluctuations @xmath33 calculated for a commensurate filling of one boson per lattice site . 
   the dashed line shows the result from the exact diagonalization for 7 lattice sites and @xmath34 bosons . 
   ( b ) order parameter @xmath35 and ( c ) compressibility @xmath36 both computed at a fixed chemical potential @xmath37 . 
   [ slso],width=302 ]    the leading corrections to the mean values of single lattice - site observables ( slso ) are of second order in @xmath25 . 
   the results for the order parameter @xmath38 , the compressibility @xmath36 and the number fluctuations @xmath39 are shown in fig . 
   [ slso ] , where the ration @xmath40 has been scaled by the dimension , to keep the same mean - field transition point . 
   the solid lines in fig . 
   [ slso ] show the results from perturbation theory , as compared to the gutzwiller result ( shown as the dashed  -  dotted line ) . 
   all three quantities show a vanishing perturbative correction both for small and large @xmath40 , as the gutzwiller wavefunction becomes a good approximation in these regimes ( for lattice dimensions @xmath41 ) 
   . as expected , deviations from the mean - field picture are strongest near the mi - sf transition , where higher - order corrections will become more and more important . 
   the order parameter @xmath42 shown in fig . 
   [ slso]b gets suppressed in the sf . 
   perturbative corrections to the gutzwiller result are particularly large in 1d ( where the order parameter vanishes in reality ) , but get smaller with increasing dimension . 
   this is not surprising , as gutzwiller is a mean - field approach and hence a better approximation for higher - dimensional systems . 
   the critical value @xmath43 is not modified within the present perturbative approach . 
   figure  [ slso]c shows the results for the compressibility    @xmath44    the results of perturbation theory ( pt ) show a decrease of the compressibility , pointing to an increasing stiffness of the sf phase induced by the short - range interaction .    finally , we computed the local particle number fluctuations @xmath45 . 
   exact diagonalization calculations have shown that the number fluctuations are changing smoothly at the mi - sf transition @xcite , whereas the gutzwiller result predicts vanishing fluctuations in the mi phase , @xmath46 ( see dashed - dotted line in fig .  [ 
   our perturbative results reproduce the non - vanishing part of @xmath39 in the mi regime and agree well with our result obtained from exact diagonalization of a one - dimensional system of @xmath47 lattice sites . 
   significant deviations from the exact 1d result are seen in the mi - sf transition regions , starting from the mean - field critical value , @xmath48 , up to values of the order of the critical values of @xmath49 usually obtained from dmrg @xcite and qmc @xcite calculations . 
   nevertheless , based on the good agreement of our perturbative result with the exact diagonalization for a large region in the mi , we conclude that the number fluctuations are mainly produced by next - neighbor particle - hole fluctuations included in perturbation theory . 
   the single - particle density matrix @xmath50 is of particular importance as it describes the correlation between the different lattice sites . 
   the correlation function @xmath51 shows off - diagonal long - range order in the sf state ( in dimensions @xmath41 ) , in contrast to the mi phase where @xmath51 decays exponentially . 
   the experimental observation @xcite of the mi transition relies on the different behavior of the density matrix in the mi and sf regimes , which can be visualized by taking absorption pictures of the freely expanding atomic cloud . assuming that the expansion time is long enough and that the gas is dilute enough ( such that atom - atom interactions can be neglected during the expansion ) , the shape of the cloud reflects the initial momentum distribution @xmath52 , which is directly given by the fourier - transform of the density matrix @xmath51 : @xmath53     as a function of site distance @xmath54 . 
   results for a homogeneous lattice of 7 sites with @xmath34 bosons have been obtained from exact diagonalization ( ed ) , from second - order perturbation theory ( pt ) , and from the gutzwiller mean - field ansatz ( gw ) . 
   ( a ) @xmath55 ( mi regime ) , ( b ) @xmath56 , ( c ) @xmath57 , ( d ) @xmath58 ( deep sf regime ) . 
   all calculations use periodic boundary conditions . 
   the mean - field value @xmath43 for the mi transition in the commensurate case with one boson per lattice site is @xmath48 in 1d ( see fig . 
   [ slso ] ) . 
   ( the mf value differs strongly from @xmath59 derived from qmc calculations @xcite or @xmath60 for dmrg calculations @xcite ) . 
   [ vgled],width=302 ]    here , @xmath61 is the fourier transform of the wannier functions @xmath62 describing the wavefunction of a single lattice site . 
   the presence of the factor @xmath61 in eq . 
   ( [ fourier ] ) provides a cutoff at high momenta . 
   the mean - field results for @xmath51 only describe the different long - range behavior in the mi and sf . for a homogeneous lattice 
   the correlation function calculated from the gutzwiller wavefunction gives @xmath63 for the diagonal elements and then drops instantly to @xmath64 for all off - diagonal elements @xmath65 . 
   short - range correlations are not reproduced by the gutzwiller approach . 
   this deviation is particularly severe in the mi , where the mean - field result predicts a completely flat momentum distribution , whereas the short - range correlations ( i.e. the exponential decay of @xmath51 ) yield smooth bumps in the expansion pattern . these can be distinguished from the @xmath66-peaks of the sf only after a sufficiently long expansion time . 
   applying perturbation theory to the gw wavefunction improves the structureless gw correlation function . in fig . 
   [ vgled ] , we have compared the results of gw mean - field theory , of pt , and of exact diagonalization ( ed ) . 
   the diagonalization has been carried out for a small 1d lattice , where it is easily feasible . 
   although there is no long - range order in the sf phase for the 1d case , where the density matrix exhibits a power - law decay towards zero , it is still reasonable to compare the short - range correlations . 
   indeed , we find a nice agreement between perturbation theory and exact results , not only for the mi ( see fig . 
   [ vgled]a ) , but also for the short - range behavior in the sf . 
   this agreement is made possible by the fact that @xmath51 decays only slowly and higher - order corrections show only negligible corrections for small lattices . 
   an example for the sf case is shown in fig . 
   [ vgled]d . 
   however , there is still a considerable difference for @xmath67 in fig . 
   [ vgled]d . 
   this is not surprising as we do pt up to second order . 
   hence correlations over a distance of three and more lattice sites are only corrected by the global mean - field correction for the infinite lattice ( see eq . 
   ( [ cpsione ] ) , eq . 
   ( [ cpsitwo ] ) , and fig . 
   [ 2ndcpsi ] in appendix [ appendixa ] ) . 
   we expect better agreement for larger lattice sizes as finite - size effects , arising in small lattices , are still considerable for @xmath68 sites used in our ed calculations .    finally , for intermediate values of @xmath40 , shown in fig . 
   [ vgled]b , c , we observe a faster drop in the off - diagonal correlations , such that higher - order contributions in the pt become more important . 
   in any dimension , the perturbation @xmath25 is no longer small at the tip of the mi - sf transition lobe ( fig .  [ vgled]c ) , and the pt breaks down . 
   nevertheless , comparing with exact diagonalization results , fig . 
   [ vgled]b shows still good agreement with ed , in contrast to fig . 
   [ vgled]c , which shows clear disagreement . 
   even though the parameters @xmath57 chosen for fig . 
   [ vgled]c are close to the mi - sf transition for 1d - lattices ( as predicted by dmrg @xcite and qmc @xcite calculations ) , pt reproduces the correct slope for the off - diagonal decay and lacks only the wrong offset from the mean field . 
   thus , even for this case , pt represents a qualitative improvement on the gutzwiller result . 
   the results for both approximations , gutzwiller and pt , are expected to become better in higher dimensions ( with the perturbative corrections diminishing in size ) .     in second - order pt as a function of @xmath69 and @xmath40 . 
   the order parameter @xmath70 vanishes inside the mott - insulating lobes , whose mean - field boundaries are given by the white dashed line . 
   plots ( a ) 
   ( f ) display the resulting momentum distribution without the wannier form factor , @xmath71 , calculated for a 2d lattice with @xmath72 lattice sites . 
   ( a)-(c ) are the gutzwiller mean - field results , and ( d)-(f ) are calculated using second - order pt . 
   the inset in ( c),(f ) shows a cut of one peak taken along @xmath73 ; dashed line for gw and solid line for the pt result . 
   arrows indicate the position of the respective plots in the @xmath74 phase diagram . 
   the parameters used are : @xmath75 , @xmath76 for ( a ) and ( d ) ; @xmath77 , @xmath78 for ( b ) and ( e ) ; and @xmath79 , @xmath80 for ( c ) and ( f ) . 
   the gray - scales of plots belonging to the same parameter set are identical . 
   expansion patterns ( a),(b),(d),and ( e ) are normalized to the peak maximum ; ( c ) and ( f ) are normalized to @xmath81 of the peak maximum.[expans ] , width=302 ]    as discussed before , the perturbative enhancement of the description of short - range correlations is expected to lead to strong consequences for the momentum distributions . as an example we discuss a set of momentum distributions @xmath71 for a homogeneous 2d lattice , fig . 
   [ expans]a - f , and compare the gutzwiller results to those improved by pt . 
   the improved pt versions , figs . 
   [ expans]d - f , show much finer structures than the mean - field results , fig . 
   [ expans]a - c . 
   pt predicts broad peaks in the mi regions down to very small values of @xmath40 , fig . 
   [ expans]e , whereas the gutzwiller result without pt shows a structureless flat distribution for the whole mi region , fig . 
   [ expans]a , b . naturally , the modifications of @xmath71 are strongest near the phase transition , figs . 
   [ expans]a and [ expans]d . 
   going towards larger values @xmath40 into the sf phase , pt gives rise to a suppression of the peaks ( inset in fig . 
   [ expans]c , f ) 
   . this suppression can be larger than @xmath82 of the original peak height and stems from the corrections to the mean field . 
   additionally , fig . 
   [ expans]f shows broad peaks induced by the inclusion of short - range correlations . 
   however , for large lattices , these broad peaks are small compared to the ( finite - size broadened ) sf @xmath66-peaks .       on a plane through the trap center for a 3d lattice with @xmath83 lattice sites . 
   ( a ) @xmath84 , @xmath85 , and @xmath86 . 
   ( b ) @xmath75 , @xmath87 , and @xmath88 . 
   ( c ) @xmath84 , @xmath87 , and @xmath86.[spd],width=302 ]    -sites in the presence of a harmonic potential . 
   ( a ) momentum distribution without the wannier form factor , @xmath71 , calculated along the @xmath89 direction . 
   filled symbols are the pt results and open symbols are the gw results . 
   the graphs for the mi surrounded by a sf shell ( diamonds ) are rescaled by a factor @xmath90 . 
   the results denoted by circles correspond to the occupation distribution of fig . 
   [ spd]a , diamonds to fig . 
   [ spd]c , and squares to fig . 
   ( b ) difference between pt and gw results . 
   [ expans3d],title="fig:",width=264 ] -sites in the presence of a harmonic potential . ( a ) momentum distribution without the wannier form factor , @xmath71 , calculated along the @xmath89 direction . 
   filled symbols are the pt results and open symbols are the gw results . 
   the graphs for the mi surrounded by a sf shell ( diamonds ) are rescaled by a factor @xmath90 . 
   the results denoted by circles correspond to the occupation distribution of fig . 
   [ spd]a , diamonds to fig . 
   [ spd]c , and squares to fig . 
   ( b ) difference between pt and gw results . 
   [ expans3d],title="fig:",width=264 ]    in contrast to what was assumed in the last section , optical lattices used in experiments are not homogeneous . magnetic or optical trapping potentials are used @xcite to confine the atomic gas to a finite volume .    the inhomogeneity caused by the trapping potential leads to slowly varying on - site energies , @xmath10 , in the mott - hubbard model eq . 
   ( [ motthubb ] ) , that can be interpreted as a spatially varying chemical potential @xmath91 . 
   consequently the lattice is in general not in a pure mi or sf phase , but shows alternating shells of sf and mi regions . 
   an example for a sf region surrounded by a mi shell is shown in fig . 
   [ spd]b .    considering the slowly varying on - site energy as 
   a spatially varying chemical potential gives a qualitative understanding of the shell structures . 
   the spatial variation of the chemical potential corresponds to a path parallel to the @xmath69 axis in the @xmath69-@xmath40-diagram . 
   starting with the potential minimum , in the trap center , and then moving off the center , decreases the effective local chemical potential . whenever a mi - sf ( sf - mi ) phase boundary is hit along the path in the @xmath69-@xmath40-diagram , a change from a mi to a sf ( sf to mi ) 
   shell appears . 
   the inclusion of short - range correlations gives rise to considerable modifications also in the presence of an inhomogeneous trapping potential . 
   examples , calculated for a 3d - lattice , are shown in fig . 
   [ expans3d ] , where an underlying harmonic potential @xmath92 was chosen . 
   the three situations considered here correspond to a case with a large mi fraction ( fig .  [ spd]a ) , a sf island surrounded by a mi shell ( fig . 
   [ spd]b ) , and a mi island surrounded by a sf phase ( fig . 
   [ spd]c ) . calculating the expansion patterns for these situations , fig . 
   [ expans3d ] , shows that the perturbative corrections arising from the short - range correlations lead to substantially different behavior in the different cases .    for the almost complete mi state we get a correction to all wavevectors @xmath93 , with a fast drop at values close to the peak center @xmath94 , which leads to a peak broadening in the expansion picture ( see circles in fig .  [ expans3d ] ) . 
   particularly large changes were found for the case of a sf island surrounded by a mi phase , fig . 
   [ expans3d ] ( squares ) . 
   again , corrections arise for all wavevectors , but , in contrast to the almost homogeneous case , the largest increase is now found for @xmath94 , with changes of about @xmath82 of the peak - maximum .    finally , the reversed situation , a mi island surrounded by a sf phase ( diamonds in fig .  [ expans3d ] ) does not show an increase of its maximum peak height but a considerable reduction ( over 5% of the peak maximum ) . 
   this is not surprising , as the majority of the lattice sites are now contributing to the sf phase , and a peak reduction was also observed for the bulk sf phase . 
   we would now like to compare our approach with the results obtained by kashurnikov _ 
   @xcite , who used qmc calculations to calculate expansion patterns for a small 3d lattice with harmonic confinement . 
   there is good qualitative agreement in all cases with high superfluid fraction ( compare for example fig . 
   [ kash]a(b ) and fig . 
   2b(c ) in ref . 
   @xcite ) . however , the features of these expansion patterns are already well reproduced using the gutzwiller mean - field ansatz alone , in particular , the satellite peak , which was discussed as a signature of the mi - sf shell structure in @xcite . 
   corrections arising from pt show a suppression of the sf peak , as was discussed above .    however , there are also considerable discrepancies to the qmc results for situations with a large mi fraction , even after implementing second - order pt . in these cases ( fig . 
   [ kash]c(d ) and fig . 
   2d(e ) in ref . 
   @xcite ) , the influence of the mi phase on the expansion picture broadens the peak and leads to a homogeneous background . the discrepancies to the qmc 
   are clearly visible in fig . 
   [ kash]c showing no peak broadening and a sattelite peak in contrast to the qmc results ( fig.2d in @xcite ) . 
   including the short range correlations perturbatively corrects the expansion pattern in the right direction , giving rise to a suppression of the sf - peak . considering the expansion pattern with the clearest mi features ( fig  [ kash]d and fig . 
   2e in ref . 
   @xcite ) , we obtain the correct peak broadening from gw / pt calculation , but a larger ratio of mi background to sf peak . 
   concerning these discrepancies , we note that the expansion pattern is highly sensitive to the value of the mean field . 
   even small deviations can lead to a change in the sf - peak height sufficient to mask the flat distribution of the mi phase . 
   we checked , however , that the observed discrepancy is not due to a lack in accuracy of our numerical calculations . 
   we therefore believe that the discrepancies between qmc and gw / pt for situations with a large mi fraction can be attributed to the insufficiency of gw and low - order pt in describing the long - range correlations in this inhomogeneous situation . 
   in addition , we note that the lattice employed in @xcite is comparatively small for the given harmonic confinement potential ( with no complete shell of empty sites at the perimeter , see insets of fig . 
   [ kash ] ) , i.e. , the choice of boundary conditions ( periodic in the case of our numerical calculations ) may have non - negligible effects on the outer lattice sites .    ) for a system of @xmath95 lattice sites considered in @xcite . 
   all plots are without the wannier form factor , @xmath96 . 
   the insets show the occupation number @xmath97 for a cut along the @xmath98-direction . 
   ( a )  @xmath99 , @xmath100 , @xmath101 . 
   ( b )  @xmath102 , @xmath103 , @xmath104 . 
   ( c )  @xmath102 , @xmath105 , @xmath106 . 
   ( d )  @xmath102 , @xmath107 , @xmath108 . 
   ( note however , that the hamiltonian used in @xcite differs from eq . 
   ( [ motthubb ] ) and hence the parameters are converted to the corresponding quantities used in our definitions . ) [ kash],width=321 ] 
   we have used perturbation theory to incorporate the effects of short - range correlations on top of the mean field solution of the bosonic mott hubbard model . 
   we derived corrections to local quantities , as well as to observable expansion patterns . 
   we numerically calculated the corrections to the mf - result , using pt up to second order , thus including correlations between next and next - nearest lattice sites . 
   modifications to the particle number fluctuations @xmath39 , arising from the pt , gave rise to the expected smooth transition of @xmath39 at the mi - sf transition . 
   moreover , comparing the pt results to the ed results for @xmath39 in 1d lattices showed good agreement for small values of @xmath40 . 
   of particular importance are the corrections to the correlation function , and thus to the expansion patterns . 
   we compared the correlation function obtained from pt with calculations from ed for small 1d - lattices and found good agreement . studying 
   the expansion patterns showed that the inclusion of the short - range correlations to the mf - ansatz gives rise to distinct modifications . 
   a broad peak can be seen in the pt results for the mi regime , visible even down to small values of @xmath40 . comparing pt and mf expansion patterns obtained for parameters in the sf region displayed a considerable suppression of the sf peak in the pt results . 
   additionally , on approaching the sf - mi transition from the sf side , broad peaks underlying the sf peaks were found in the pt - expansion patterns . including a harmonic confinement potential leads to situations where sf and mi regions coexist . 
   hence the perturbative corrections to the expansion pattern become more complex . 
   we studied lattices with different underlying harmonic traps giving rise to different constellations of sf - mi regions , finding modifications of up to @xmath82 of the peak maximum in the expansion patterns . 
   we would like to thank krishnendu sengupta for discussions . 
   our work was supported by the swiss nsf and the nccr nanoscience , as well as a dfg research fellowship ( f.m . ) . 
   schematic diagram illustrating the term @xmath109 appearing in the first - order correction @xmath110 of the density matrix , where @xmath111 is the perturbation connecting sites @xmath4 and @xmath112 . 
   the diagram shows the lattice sites in the horizontal direction . 
   the different steps needed to obtain the matrix element are shown vertically . 
   open circles denote the ground state ( gs ) of the given lattice site , while filled circles are excited states ( es ) of this site , for the mean - field hamiltonian @xmath21.,width=264 ] 
   we use standard stationary perturbation theory @xcite to calculate the corrections to the mean - field results induced by the perturbation @xmath113 where we introduce the new operators @xmath114 the expectation value @xmath115 is taken with respect to the mf wavefunction @xmath116 . defining the operator for the energy denominator as @xmath117 the expectation value of an observable @xmath118 including all corrections up to 
   second order is @xmath119 the first line in eq . 
   ( [ perta ] ) is the mf - result @xmath120 , followed by two contributions which are the first - order corrections . 
   lines two and three in eq . 
   ( [ perta ] ) are the second - order corrections to the mean value .      as an example we will discuss the corrections to the density matrix @xmath121 . 
   the first - order correction to the density matrix is : @xmath122 for two different lattice sites @xmath123 , we find @xmath124 since the gutzwiller ground state is a product state . as a consequence , 
   the contributions @xmath125 and @xmath126 vanish . 
   the only remaining contributions to @xmath127 stem from @xmath128 and @xmath129 .      the graph shows a decomposition of the matrix element , with each row showing the wavefunction at an intermediate step in the evaluation of the matrix element . 
   as we deal with product states @xmath130 we represent the wavefunction by a row of circles , where each circle denotes the state @xmath131 of a particular lattice site @xmath4 . 
   open circles in fig . 
   [ diagramm ] denote a lattice site in its ground state ( gs ) , filled circles refer to an excited state ( es ) of this particular lattice site , with respect to the local mean - field hamiltonian . 
   note that , in general , this can be an arbitrarily highly excited state ( although higher contributions are suppressed by the energy denominator , and a cutoff is used in practice ) . 
   starting with a row of open circles , denoting the gs , @xmath116 , each following row corresponds to the state after the action of @xmath25 or @xmath132 , as indicated on the left side of the graph . as all matrix elements in eq . 
   ( [ perta ] ) can be expressed in terms of a gs expectation value @xmath115 and a sequence of @xmath25 and @xmath132 operators , the first and last row must always be a line of open circles . 
   let us consider for instance the second term in eq . 
   ( [ firstorder ] ) @xmath133 reading the graph in fig . 
   [ diagramm ] from top to bottom corresponds to reading the matrix element from right to left . 
   starting with the gs , @xmath134 , the first row consists of open circles . 
   the second row shows the state after the action of the perturbation @xmath25 . 
   acting with @xmath135 to the right onto the gs results in a state latexmath:[\[v_{ij } |g_0\rangle = \sum_{\alpha,\beta } f_{\alpha,\beta }    @xmath112 are neighboring lattice sites and @xmath137 denotes the state with lattice site @xmath4 ( @xmath112 ) in the exited state @xmath138 ( @xmath139 ) and all other sites in their gs . 
   thus the second row shows the lattice sites @xmath4 and @xmath140 in an excited state ( filled circle ) , as the perturbation , eq . 
   ( [ pert ] ) , allows only next neighbor interactions . 
   finally , the action of @xmath132 has to bring the excited states back to the gs , in order to get a non - vanishing contribution . 
   therefore , in first order pt , only next neighbor corrections to the correlation function arise , as the final row must represent the ground state @xmath141 again . 
   the graph representing the remaining first term in eq . 
   ( [ firstorder ] ) is obtained by rotating the graph in fig . 
   [ diagramm ] by @xmath142 .     and @xmath143 to the second - order correction of the density matrix . 
   ( a ) terms representing contributions of the form eq . 
   ( [ cpsione ] ) . 
   ( b ) contributions of the form eq . 
   ( [ cpsitwo ] ) . 
   all notations are the same as in fig . 
   [ diagramm ] . 
   [ 2ndcpsi],width=302 ]      rewriting the second - order corrections to the density matrix , @xmath144 , in terms of the operators @xmath145 and @xmath146 gives eq . 
   ( [ atoc ] ) , but with @xmath147 replaced by @xmath148 . 
   in contrast to the first - order corrections , now the terms proportional to @xmath149 and @xmath150 also give non - vanishing contributions . 
   using eq . 
   ( [ perta ] ) we obtain : @xmath151 here , the primed sums run over all neighbors @xmath93 to site @xmath112 . 
   the corresponding subset of graphs for @xmath152 and @xmath153 are given by fig . 
   [ 2ndcpsi]a and fig . 
   [ 2ndcpsi]b for eq . 
   ( [ cpsione ] ) and eq . 
   ( [ cpsitwo ] ) respectively . 
   note that all terms of eq . 
   ( [ cpsione ] ) and eq . 
   ( [ cpsitwo ] ) give a correction to all matrix elements of the density matrix independent of the distance between the lattice sites . 
   we can understand these contributions as a modification to the mf value of the density matrix . 
   graphs showing second - order corrections arising from @xmath154 . diagrams ( a ) and ( b ) are coming from direct - neighbor contributions as given by eq . 
   ( [ 2ndccone ] ) and eq . 
   ( [ 2ndcctwo ] ) respectively . diagrams ( c ) and ( d ) are next - nearest - neighbor contributions : ( c ) corresponds to eq . 
   ( [ 2ndccthree ] ) and eq . 
   ( [ 2ndccfour ] ) ; ( d ) corresponds to eq . 
   ( [ 2ndccfive]).,width=264 ]      * lattice site @xmath4 and @xmath112 being direct neighbors . 
   in this case we get @xmath155 corrections for eq . 
   ( [ 2ndccone ] ) and eq . 
   ( [ 2ndcctwo ] ) are shown in fig .  [ 2ndcc]a and fig . 
   [ 2ndcc]b , respectively . 
   * configurations corresponding to two lattice sites @xmath156 connected by two successive hopping steps via site @xmath93 . 
   this gives rise to six contributions : + @xmath157 + an example for the contributions arising from eq . 
   ( [ 2ndccthree ] ) and eq . 
   ( [ 2ndccfour ] ) is shown in fig . 
   [ 2ndcc]c . the last term , eq . 
   ( [ 2ndccfive ] ) , has the representation shown in fig . 
   [ 2ndcc]d . note that for lattices with dimensions @xmath41 , the sites @xmath4,@xmath112 and @xmath93 need not necessarily form a straight line but can form a chevron .
summary: we study the mott - insulator transition of bosonic atoms in optical lattices . 
   using perturbation theory , we analyze the deviations from the mean - field gutzwiller ansatz , which become appreciable for intermediate values of the ratio between hopping amplitude and interaction energy . 
   we discuss corrections to number fluctuations , order parameter , and compressibility . 
   in particular , we improve the description of the short - range correlations in the one - particle density matrix . 
   these corrections are important for experimentally observed expansion patterns , both for bulk lattices and in a confining trap potential .