Coherent driving and freezing of bosonic matter wave in an optical Lieb lattice

See allHide authors and affiliations

Science Advances  20 Nov 2015:
Vol. 1, no. 10, e1500854
DOI: 10.1126/sciadv.1500854


Although kinetic energy of a massive particle generally has quadratic dependence on its momentum, a flat, dispersionless energy band is realized in crystals with specific lattice structures. Such macroscopic degeneracy causes the emergence of localized eigenstates and has been a key concept in the context of itinerant ferromagnetism. We report the realization of a “Lieb lattice” configuration with an optical lattice, which has a flat energy band as the first excited state. Our optical lattice potential has various degrees of freedom in its manipulation, which enables coherent transfer of a Bose-Einstein condensate into the flat band. In addition to measuring lifetime of the flat band population for different tight-binding parameters, we investigate the inter-sublattice dynamics of the system by projecting the sublattice population onto the band population. This measurement clearly shows the formation of the localized state with the specific sublattice decoupled in the flat band, and even detects the presence of flat-band breaking perturbations, resulting in the delocalization. Our results will open up the possibilities of exploring the physics of flat bands with a highly controllable quantum system.

  • optical lattice
  • Bose-Einstein condensation
  • quantum simulation
  • flat band
  • ultracold atomic gases


Many-body properties of a quantum system show drastic changes according to the geometry of an underlying lattice structure. One of the textbook examples is an antiferromagnet on a frustrated lattice geometry (1), where the geometric frustration prevents spins from Néel ordering and the system exhibits more nontrivial, correlated ground states. A dispersionless flat band realized by specific lattice geometry can also induce frustration of kinetic energy, resulting in a degeneracy of a macroscopic number of momentum eigenstates. Flat bands play an important role in the study of itinerant ferromagnetism because the presence of interaction lifts the bulk degeneracy and chooses the ferromagnetic ground state (24).

A special type of lattice structure known as a Lieb lattice, also referred as a decorated square lattice, has a flat band as the second (first excited) band. It consists of two sublattices: one of them forms a standard square lattice (the A sublattice in Fig. 1A), and the other lies on every side of the square. For convenience, we further divide the latter into the B and C sublattices. The single-particle energy spectrum in the tight-binding limit (Fig. 1B) has the characteristic flat band and the Dirac cone on the corner of the Brillouin zone. This Lieb lattice satisfies the criteria for the occurrence of Lieb’s ferrimagnetism, which states that the half-filled spin-1/2 fermions exhibit nonzero magnetization for a positive on-site interaction (2). Also for bosonic systems, a flat band proposes a fascinating question of whether condensation is possible in the presence of kinetic energy frustration. Theoretical investigation predicts supersolid order for a flat band (5). Here, we also note that a Lieb lattice naturally has the lowest three bands close to each other, providing rich interband physics. The structure of the Lieb lattice is identical to the three-band d-p model, which describes the CuO2 plane of high-Tc superconductors (68).

Fig. 1 Optical Lieb lattice.

(A) Lieb lattice. A unit cell is indicated by the green square. (B) Tight-binding energy band structure of the Lieb lattice. (C) Experimental realization of the Lieb lattice. Black arrows indicate polarizations of the lattice beams. (D) Lattice potential for (slong, sshort, sdiag) = (8, 8, 9.5) at φx = φz = 0 and ψ = π/2. (E) Band structures of the optical Lieb lattice at (slong, sshort, sdiag) = (8, 8, 9.5) (red dashed), (34, 34, 37.4) (solid black).

Ultracold atomic gases in optical lattices have had great success in realizing controllable quantum many-body systems described by well-defined theoretical models of interest, such as the Hubbard model (9, 10). Besides the simple cubic configuration, increasing experimental efforts have been made to create and investigate nonstandard optical lattices that have unique geometric features (1119). The Lieb lattice or its one-dimensional (1D) analog (sawtooth lattice) was recently realized in a photonic lattice (2022) and polaritonic systems (23). However, optical lattice realization has definite advantages: simple and strong interactions, dynamical controllability of system parameters, and availability of both bosonic and fermionic systems. Above all, it can be directly connected to essential models containing key physical concepts.

Here, we demonstrate novel manipulation and detection of an atomic Bose-Einstein condensate (BEC) in a flat band by developing a dynamically controllable optical Lieb lattice. In particular, we invent a method for actively engineering the population and phase on each lattice site, which enables us to coherently transfer atoms into the flat band and observe frozen motion of atoms localized on a specific sublattice. In addition, almost arbitrary superposition of band eigenstates can be prepared, which drives coherent oscillation modes and enables us to map out the characteristic band structure. Novel controllability of our system is highlighted by an experiment that controls the localization and delocalization of an atoms and detects the presence of flat-band breaking perturbations. This work paves the way to a new regime of experimental study of flat band physics with cold atoms.


Formation of an optical Lieb lattice

We construct the optical Lieb lattice by superimposing three types of optical lattices (see Fig. 1, C and D), leading to the potentialEmbedded Image(1)where z indicates the direction of gravity. Here, kL = 2π/λ is a wave number of a long lattice (with a depth Vlong), for which we choose λ = 1064 nm. A short lattice (Vshort) is formed by laser beams at 532 nm. A diagonal lattice (Vdiag) with the wave number Embedded Image is realized by interference of the mutually orthogonal laser beams at 532 nm along the x and z directions. Compared to the proposals of Shen et al. and Apaja et al. (24, 25), Eq. 1 lacks a diagonal lattice with x + z spatial dependence. Although this causes a slightly larger discrepancy from the ideal tight-binding energy bands, the sufficiently deep lattices can reproduce the desired energy spectrum including the flat band, as seen in Fig. 1E. In the following, we specify each lattice depth in unit of long lattice recoil energy as (slong, sshort, sdiag) = (Vlong, Vshort, Vdiag)/ER(1064), where ER(1064) = ħ2kL2/(2m) and m is the atomic mass of 174Yb. The motion of free bosons in the x-z plane is governed by tight-binding HamiltonianEmbedded Image(2)where âi is the annihilation operator on site i and Embedded Image. The nearest-neighbor hopping amplitude J is mainly determined by Vshort, whereas the other lattice depths set an energy offset ES of each sublattice. Excepting the contribution from the zero-point energies of each potential well, they are approximately given by EA ~ −Vlong(x)Vlong(z), EB ~ −Vlong(z)Vdiag, and EC ~ −Vlong(x)Vdiag, which can be independently controlled by tuning the lattice depths. In the y direction, which is perpendicular to the Lieb lattice plane, the atoms are weakly confined in a harmonic trap (1D tube configuration), unless otherwise specified. In addition to each lattice depth, three phase parameters should be set to φx = φz = 0 and ψ = π/2 to realize the characteristic three-sublattice structure of the Lieb lattice depicted in Fig. 1D (Materials and Methods).

Loading BEC into a flat band by phase imprinting

One of the fundamental properties of flat bands is the localization of the wave function as a consequence of quantum-mechanical interference of traveling matter waves. The localization is due to a purely geometric effect, as we briefly explain below. The Hilbert space for the Lieb lattice in the tight-binding regime is spanned by the momentum eigenstates of each sublattice |k, S〉 (S = A, B, C). Nearest-neighbor tunneling induces the quasimomentum-dependent coupling KAB = −2J cos(kxd/2) between the A and B sublattices, and similarly KAC = −2J cos(kzd/2) between the A and C sublattices, where d = 532 nm is the lattice periodicity. The flat band states are the zero-energy eigenstates cos θ|k, B〉 − sin θ|k, C〉 with tan θ = KAB/KAC, which have no amplitude on the A sublattice. Consequently, a wave packet composed of the flat band states remains localized, as the tunneling from a B site and the tunneling from a C site destructively interfere on the adjacent A site. We explore this nature in the following experiments. It is worth noting that the flat band in the Lieb lattice is mathematically equivalent to “dark states” of laser-coupled Λ-type three-level systems in atomic physics (26). In this analogy, three sublattices correspond to the basis of three levels, tunneling amplitudes serve as laser-induced coupling, and the energy difference of each sublattice plays a role in the detuning of laser.

In a Lieb lattice, a flat band is realized as the first excited band; hence, a BEC loaded adiabatically into an optical Lieb lattice is not populated in the flat band. However, tunability of our optical Lieb lattice enables us to coherently transfer the population in the lowest band into the flat band by phase imprinting (see Fig. 2A). The scheme is easily understood by considering tight-binding wave functions in each band. At zero quasimomentum and in the equal-offset condition EA = EB = EC, a simple calculation gives Embedded Image, |2nd〉 = |B〉 − |C〉, and Embedded Image from the 1st to the 3rd band, where we omit the momentum indices from the sublattice eigenstates (see also section S2). Taking advantage of rich controllability in our lattice potential, we can smoothly modify these eigenstates. With sufficiently large Vdiag (equivalently with large EAEB,C), the lowest Bloch state has essentially no amplitude in the A sublattice, allowing realization of the |B〉 + |C〉 state. Next, we apply sudden change in one of the long lattice, Vlong(z). This creates the energy difference between the B and C sublattices, and the relative phase of the condensate wave function starts to evolve with a period 2πħ/(ECEB). On the basis of the initial band structure, this time evolution is a coherent oscillation between |1st〉 and |2nd〉.

Fig. 2 Coherent band transfer.

(A) Principle of the transferring method. (B) Absorption images reveal the coherent oscillations between the |B〉 + |C〉 and |B〉 − |C〉 states. In the upper left image, the first three Brillouin zones are displayed by white, green, and red lines, respectively. (C) Oscillating behavior of the band population during phase imprinting in the absence of lattice confinement along the y direction (blue circles) and with lattice confinement −Vycos2(2kLy) (red squares). Solid lines are the fit results using the single-particle solution of the Schrödinger equations (see main text). Error bars denote SD of three independent measurements.

The explicit procedure of loading and detecting a condensate in the flat band is as follows. We adiabatically load a BEC of 2 × 104 174Yb atoms into the Lieb lattice with (slong, sshort, sdiag) = (8, 8, 20) and apply sudden increase of slong(z) to 26.4 for variable duration. At the same time, we ramp sshort up to 20 to prevent tunneling during the band transfer. After this sequence, we return the lattice depths to the initial values and perform adiabatic turning off of the lattice potential to map quasimomentum to free-particle momentum (band mapping) (27, 28). Figure 2B shows the absorption images taken after 14 ms of the ballistic expansion, which reveal the interband dynamics of a condensate. At zero quasimomentum, atoms in the 2nd and 3rd band are mapped to the same point of the Brillouin zone. In addition, the finite spread of the condensate makes it difficult to precisely distinguish the population in the 2nd Brillouin zone from other neighboring zones. Therefore, instead of plotting the population in the 2nd Brillouin zone, here we count atoms in the 1st Brillouin zone and show the fraction of atoms in the other higher bands in Fig. 2C. For momentum space analysis, see section S3. We fit the data with the function in the form a exp(−t/τ) Fth(t) + b[1 − exp(−t/τ)], where Fth(t) is the numerical solution of the single-particle Schrödinger equation. In fitting the data, we adopt slong(z) during the band transfer as a free parameter and obtain the best fit with slong(z) = 25.5, close to the expected value of 26.4. Although the oscillations involve nonnegligible contributions from the higher bands, at the half period of the first cycle, we expect that >75% of atoms are transferred to the 2nd band. This transfer method is also applicable in the presence of the lattice confinement along the y axis, though the decay time of the oscillation τ = 86 ± 7 μs is much shorter than the case of a weak harmonic confinement (1D tube), τ = 260 ± 10 μs.

Relaxation dynamics of a flat band

We measure the lifetime of atoms in the 2nd band of the optical Lieb lattice. After transferring to the 2nd band, we change the depth sdiag of the diagonal lattice to control the energy gap between the 1st and 2nd bands. As well as the band gap (29), the lifetime of a quantum gas in the excited band is strongly affected by the density overlap with the states in the lower bands (30). As we increase sdiag, the average gap between the 1st and 2nd bands becomes smaller and, at the same time, their density profiles become similar to each other. In the opposite limit of shallow sdiag, the band gap increases and two bands have no density overlap, because the lowest band mostly consists of the A sublattice. We take a variable hold time in the lattice, followed by band mapping to count the atom number in the excited bands. Typical absorption images are shown in Fig. 3A. The decay curves displayed in Fig. 3B show expected behavior of increasing lifetime with decreasing sdiag. In addition, increasing the gap makes the dynamics more clearly separate into two processes: decay of the condensate within the 2nd band (middle image of Fig. 3A) and decay of atoms into the lowest band (bottom image). We find that the curve is well fitted by a double exponential with the form a1 exp(−t1) + a2 exp(−t2) + b. The fast component τ1 shows only weak dependence on sdiag, whereas the slow component τ2 shows more than 20-fold changes from the smallest to the largest sdiag. We also extract the lifetime of the condensate in the 2nd band, τc, by counting atoms on the corner of the 2nd Brillouin zone (Fig. 3A), and find similar behavior with τ1. This implies that the initial fast decay is related to the decay of the condensate, which involves the decay to the lower band with faster time constant compared with the noncondensed atoms.

Fig. 3 Lifetime of atoms in the flat band.

(A) Absorption images for the lifetime measurement of the 2nd band with three different hold times, taken after 14-ms time of flight. The diagonal lattice depth is sdiag = 9.5. The first three Brillouin zones are indicated by the white dashed lines. In the top image, the areas used to evaluate the lifetime of a condensate (τc) are also displayed with the red squares. (B) Decay of the flat band at (slong, sshort) = (8, 8) and variable sdiag. Solid lines are the fit results with double-exponential curves. Error bars denote the SD of three independent measurements. (C) Lifetime of the flat band. τ1,2 are the fast and slow decay time obtained from the data shown in (B), respectively. τc is the e−1 lifetime of condensates. Error bars represent fitting error.

Localization of a wave function in a flat band

As described above, the most intriguing property of a flat band is the localization of the wave function at certain sublattice sites. In the case of the Lieb lattice, the wave function of the flat band vanishes on the A sublattice. Here, we reveal this property by observing the tunneling dynamics of a Bose gas initially condensed at the | − 〉 = |B〉 − |C〉 state, and compare it to the dynamics of the state with opposite relative phase, | + 〉 = |B〉 + |C〉. To observe real-space dynamics of the system, we perform projection measurement of the occupation number in each sublattice, which we call sublattice mapping. In this method, we first quickly change the lattice potential to ((slong(x), slong(z)), sshort, sdiag) = ((8,14), 20, 0). In this configuration, all three sublattices are energetically well separated from one another and the lowest three bands consist of the A, B, and C sublattice, respectively. This maps sublattice occupations to band occupations, which can then be measured by band mapping technique. Figure 4A shows the demonstration of this method, in cases where atoms occupy only one of the sublattices. Note that the populations in the B and C sublattices are mapped to the 2nd Brillouin zones for the 1D lattice along the x and z axis, respectively. This is because the turning off of the diagonal lattice decouples these two directions and the fundamental bands are labeled by the combination of band indices of 1D lattices.

Fig. 4 Tunneling dynamics in the Lieb lattice.

(A) Demonstrating the measurement of sublattice occupancy. Here, sublattice mapping technique is applied to atoms loaded into (left) ((slong(x),slong(z)), sshort, sdiag) = ((8,8), 8, 0), (middle) ((2,8), 8, 19), and (right) ((8,2), 8, 19), corresponding to atoms in A, B, and C sites, respectively. (B) Measured tunneling dynamics of |+〉 and |−〉 initial states in the Lieb lattice with (slong, sshort, sdiag) = (8, 8, 9.5). Solid lines are the fits to the experimental data with damped sinusoidal oscillation (for |+〉) and double exponentials (for |−〉). Inset shows dynamics of the |−〉 state for longer hold times. Error bars denote SD. Illustration of tunneling process for each initial state is also shown on the right-hand side. (C) Frequencies of coherent intersite oscillations. Solid lines are the calculated band gap between the 1st and 2nd (red) and the 1st and 3rd bands (blue). Error bars denote fitting error. (D) Bending flat band. Dynamics of the |−〉 state in the presence of imbalance Δslong = slong(x)slong(z) shows restoration of coherent dynamics. Error bars denote SD.

We prepare the initial state |+〉 by simply loading a BEC into the Lieb lattice with deep Vdiag. On the other hand, the |−〉 state is obtained by applying the band transfer method to the |+〉 state. After changing the lattice depths to satisfy the equal-offset condition EA = EB = EC , dynamics of these initial states is measured by the sublattice mapping. As shown in Fig. 4B, we reveal qualitatively different behaviors of these two states: the |−〉 state shows a significant suppression of the A sublattice occupancy, indicating the freezing of the tunneling dynamics to the A sublattice from the |−〉 state with only a slow decay to the A sublattice, whereas the |+〉 state exhibits coherent oscillations between the A and (BC) sublattices. This clearly features the geometric structure of the Lieb lattice mentioned above. Double-exponential fit to the data for the |−〉 initial state yields τ1 = 0.36 ± 0.04 ms and τ2 = 5.5 ± 0.9 ms, indicating that the leakage to the A sublattice is caused by the decay to the lowest band.

In the Bloch basis, the state |+〉 is expressed as |1st〉 − |3rd〉 and its time evolution is driven by the band gap ΔE3,1 which equals Embedded Image in the tight-binding limit. After a half-period πħE3,1, the state evolves to |A〉 = |1st〉 + |3rd〉, leading to coherent tunneling to the A sublattice. Similarly, it is possible to arrange the initial lattice depths so that the lowest Bloch state has the maximum overlap with a certain superposition of |1st〉 and |2nd〉. Sudden potential change to the Lieb lattice drives oscillation between the B and C sublattices, whose frequency gives the band gap ΔE2,1. We fit these data with a damped sinusoidal oscillation and compare the extracted frequency with the result of single-particle band calculations (see Fig. 4C). Qualitative behavior is well reproduced, whereas quantitative discrepancies are found. This is caused by interactions, as we present a systematic study of the density dependence of the oscillation frequency in section S4.

We further investigate the tunneling dynamics of the |−〉 initial state by adding the perturbations that destroy the flatness of the second band. The flatness is robust against the independent change of nearest-neighbor tunneling amplitudes Jx, Jz along the x and z directions, and energy offset EC, just as a dark state in a Λ-type three-level system persists regardless of laser intensities and detuning from the excited state. However, if the energy difference between the B and C sublattices is introduced—the two-photon Raman off-resonant case—the flat band is destroyed. Note that the finite EBEC induces population in the A sublattice even at k = 0. On the other hand, the direct diagonal tunneling between the B and C sublattices, which is another flat-band breaking term existing in our system, keeps a dark state at k = 0 provided Jx = Jz. We create the energy difference by introducing the imbalance of Δslong = slong(x)slong(z). Figure 4D shows the time dependence of the A sublattice population for the |−〉 initial state. It can be clearly seen that the coherent tunneling dynamics starts to grow as the lattice parameters deviate from the flat-band condition Δslong = 0.


Here, we have successfully implemented the Lieb lattice for ultracold atomic gases and observed the characteristic dynamics of a condensate, including the freeze of the motion in the flat band. This work shows an important ability of our optical lattice setup to make a connection between theory and experiment. The highly controllable lattice allows us to study both a nearly complete flat band where prominent theoretical works have been established, and intentionally perturbed, imperfect flat bands that are relevant to real materials. Relatively short lifetime of atoms in the flat band was observed, although it can be made longer by increasing the band gap to the lowest band. Using Fermi gases with the Fermi energy lying at the flat band can avoid the lifetime problem and will provide an ideal playground for investigating flat band ferromagnetism (3134) and topological phases with artificial gauge fields (35).


Preparation of 174Yb BEC

After collecting about 107 atoms with a magneto-optical trap with the intercombination transition, the atoms were transferred to a crossed optical trap. Then, we performed an evaporative cooling, resulting in an almost pure BEC with about 105 atoms with no discernible thermal component.

All of the optical lattice experiments presented in this paper were subject to additional weak confinement due to a crossed optical dipole trap operating at 532 nm. The Gaussian shape of laser beams for the trap and lattices imposed a harmonic confinement on atoms, whose frequencies are (ωx′, ωy′, ωz)/2π = (147, 37, 105) Hz at the lattice depths of (slong, sshort, sdiag) = (8, 8, 9.5). Here, the x′ and y′ axes were tilted from the lattice axes (x and y) by 45° in the same plane.

Construction of optical Lieb lattice

The relative phases between the long and short lattices (φx, φz) can be adjusted by changing the frequency difference between these lattice beams (36). The proper frequencies that realize the Lieb lattice (φx = φz = 0) were determined by analyzing the momentum distribution of a 174Yb BEC released from the lattice, as in the case of the parameter ψ of the diagonal lattice (section S1). The relative phase between the long and short lattices at the position of atoms depends on the optical path lengths from common retro-reflection mirrors, and in general, two phases φx and φz are not equal. We shifted the frequency of the long lattice laser by an acousto-optic modulator (AOM) inserted in the path for the z axis to simultaneously realize φx = φz = 0. Optimal frequency difference was sensitive to the alignment of the lattice beams, and day-by-day calibration of the phases was needed. Typical drift of the required radio frequency for the compensation AOM was within 5 MHz.

To stabilize the phase ψ, we constructed a Michelson interferometer along the optical path of the diagonal lattice with the frequency-stabilized 507-nm laser. The interferometer had two piezoelectric transducer–mounted mirrors: one was shared with the lattice laser beam for phase stabilization, and another was used to shift the phase over the range 10π, with stabilization kept active. The short-term stability of ψ was estimated to be ±0.007π. The last few optics in front of the chamber were outside of the active stabilization, which caused slow drift of ψ due to changes of environment such as temperature. The typical phase drift was 0.05π per hour, and all measurements of sequential data set was finished within 20 min of the last phase calibration.

At the proper phase parameters φx = φz = 0 and ψ = π/2, the potential depth at the center of each site becomes equal when Vlong = Vshort = Vdiag. In this condition, however, the energy offset EA became lower than EB and EC because of the difference in the zero-point energies. We searched optimal Vdiag by single-particle band calculation (see also section S2 for the derivation of Hubbard parameters).

Throughout the experiment, the Hubbard parameters were set to weakly interacting regime. Without a lattice confinement along the y axis, the typical value of the renormalized on-site interaction (37) was U/J ~ 0.02. Even with a lattice confinement Vy = 10ER(532), the system remained superfluid regime U/J = 1.1, well below the critical point for the 2D Mott transition. To approach the critical point (U/J ~ 15), the lattice depths should be as deep as (slong, sshort, sdiag) = (34, 34, 37.4) and Vy = 15ER(532).

Band occupation measurement and sublattice mapping

To measure the quasimomentum distribution of atoms, we turned off all the lattice potentials with an exponential formEmbedded Image(3)and T = 0.6 ms. The dipole trap was kept constant during band mapping to prevent the movement of the trap center due to gravity and suddenly turned off at t = T. Because of the relatively heavy mass of Yb, the existence of harmonic confinement imposed severe restriction on the choice of mapping time T. We find that T > 1 ms causes considerable deformation of the distribution, whereas T > 1.5 ms is desirable to suppress interband transition. Because of this non-adiabaticity, up to 20% of atoms occupying a certain Brillouin zone were detected in its neighboring zones, depending on the shape of the observed quasimomentum distribution.


Supplementary material for this article is available at

Section S1. Calibration of the relative phase.

Section S2. Tight-binding model for the optical Lieb lattice.

Section S3. Momentum distributions in coherent band transfer.

Section S4. Effect of interactions on inter-sublattice oscillations of a BEC.

Fig. S1. Phase dependence of a time-of-flight signal.

Fig. S2. Tunneling parameters in the optical Lieb lattice.

Fig. S3. Wannier functions of the optical Lieb lattice.

Fig. S4. Momentum space observation of coherent band transfer.

Fig. S5. Density dependence of oscillation frequency.

Table S1. Initial conditions for the inter-sublattice oscillations.

Reference (38)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank K. Noda, K. Inaba, M. Yamashita, I. Danshita, S, Tsuchiya, C. Sato, S. Capponi, Z. Wei, and Q. Zhou for valuable discussions. Funding: This work was supported by the Grant-in-Aid for Scientific Research of the Japan Society for the Promotion of Science (25220711 and 26247064) and the Impulsing Paradigm Change through Disruptive Technologies (ImPACT) program. Author contributions: S.T., H.O., and T.I. performed the experiment and analyzed the data. T.N. and S.N. contributed to building up the optical Lieb lattice. Y.T. supervised the whole project. All authors discussed results and contributed to writing the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from S.T. (taie{at}
View Abstract

Navigate This Article