text:
# Amorphous silica between confining walls and under shear: a computer simulation study

## I Introduction

The transport properties of atomic fluids that move through nanoscopic slits between solid walls has been an issue of constant interest over the past decades . Many of the simulation studies investigating such fluids in Couette or Poisseuille flow have considered simple monoatomic systems in which the particles interact via Lennard–Jones (LJ) type forces (see, e.g., ). One of the remarkable results of these studies was that velocity profiles as predicted by hydrodynamics can be observed even if the distance between the confining walls is of the order of only ten molecular diameters. But one has to keep in mind that the monoatomic LJ systems for which this was found have typically a relatively small viscosity (of the order of $`10^4`$ to $`10^3`$ Poise in Argon units) and furthermore, these systems crystallize very easily if one decreases the temperature below the melting point. An open question is what happens with more complicated sheared fluids that are confined between walls, i.e. especially liquids that have a structure that is not similar to the local structure of a closed–packed system of hard spheres like the aforementioned LJ systems. One class of more complex materials for which flow properties in confinement have already been studied to some extent are polymer melts . Very recently it was shown  for such a system that its properties under shear can be well understood within a hydrodynamic description.

In this paper we are interested in amorphous silica (SiO2) which is the paradigm of a network forming liquid. We consider a SiO2 melt that is confined between two crystalline walls and on which a gravitational–like field is exerted in order to yield a stationary Poisseuille flow. It is difficult to address this problem in a simulation because model potentials of silica contain long–ranged Coulomb terms  for which the corresponding force calculations become very time–consuming for a three dimensional system in which periodic boundary conditions are considered only in two directions. This stems from the fact that the Fourier part of the Ewald sums can no longer be calculated by a single loop over the number of particles $`N`$ as in the case of periodic boundary conditions in three dimensions, but one has to compute a double loop that scales with $`N^2`$ . Therefore, very recently, many efforts were undertaken to develop alternative methods  for the summation of “two–dimensional” Coulomb interactions. But most of these new approaches are not established yet and extensive tests that demonstrate that their accuracy is comparable to the “brute force” $`N^2`$–methods are still lacking. Moreover, differently from monoatomic LJ systems, silica is a good glassformer that exhibits a slow dynamics even above the melting temperature. Note that at the experimental melting temperature, $`T=2000`$ K, amorphous silica has a viscosity of about $`10^7`$ Poise . Systems with such a high viscosity cannot be equilibrated nowadays by means of a molecular dynamics computer simulation.

The structure of silica is that of a disordered tetrahedral network in which SiO4 tetrahedra are connected to each other in that they share an oxygen atom at their corner. Extensive simulation studies of free silica surfaces, that used both classical molecular dynamics as well as Car–Parrinello simulations  gave evidence for a different packing of the SiO4 tetrahedra at the surface from that in the bulk. This different packing is manifested very clearly in the ring structure (a ring is a closed loop formed by consecutive Si–O elements). At the surface there is a relatively high probability to find short rings such as two– and three–membered rings whereas in the bulk the occurence of such rings is very rare even at relatively high temperatures. Note that the occurrence of short rings on a free silica surface has been confirmed experimentally by means of atomic force microscopy . A detailed understanding of the structure near the wall is important to give insight into the microscopic origin of the wall–fluid boundary conditions if a flow field is introduced. This might be especially important in the case of silica since it does not form a closed–packed structure and thus the behavior could be very different from that of a hard–sphere like system. Indeed we show below that there is a rearrangement in the ring distribution near the walls if the shear is strong enough.

An important issue for the understanding of a confined fluid under shear are the boundary conditions that enter in the hydrodynamic description in order to take into account the interactions of the fluid with the wall. In principle, these boundary conditions can be derived from the microscopic equations of motion for the fluid and the wall particles. But in practice one often introduces as an ad hoc assumption stick boundary conditions. However, it has been found in simulations of simple fluids that there is often a slip motion of the fluid at the wall . So in general one has to consider slip–stick boundary conditions for a hydrodynamic analysis of fluid flow between walls. Bocquet and Barrat  have shown in a microscopic theory, in which they have related quantities like the slip length with equilibrium time correlation functions via Green–Kubo formulas, that slip occurs if the walls do not provide a strong corrugation. E.g., a large slip effect is observed if the fluid partially wets the solid walls. In a recent simulation of Sokhan et al.  it has been even demonstrated for the fluid flow through a carbon tube that slip–stick boundary conditions are the generic case if one uses realistic parameters in the model potential for the wall–fluid interactions. It is one of the issues of this paper to characterize the structure of the silica melt near the wall and to see how the structure is related to the observed boundary conditions.

Another issue that we address in this paper is to discuss to what extent the velocity profiles from our simulations can be used to estimate the shear viscosity $`\eta `$ of a highly viscous fluid such as silica. An alternative method is to calculate $`\eta `$ from the appropriate Green–Kubo formula  which expresses $`\eta `$ as the time integral over a stress–stress autocorrelation function. We applied the latter approach recently in extensive bulk simulations of silica in which we were able to determine $`\eta `$ in the temperature range 6100 K $`T`$ 3000 K whereby at $`T=3000`$ K the visosity is of the order of $`60`$ Poise . In the present paper we use the latter results as a reference to see how reliable the estimates of $`\eta `$ from our non–equilibrium simulations are. Such a reference is important since it is not obvious that the observed velocity profiles justify a “hydrodynamic analysis” keeping in mind that the distance between the confining walls is on a nanoscopic scale.

The rest of the paper is organized as follows: In Sec. II we describe how we model the walls and the silica melt, and we give the details of the simulation. Sec. III is then devoted to the presentation of the results and eventually, in Sec. IV we summarize and discuss them.

## II Model and Details of the Simulation

A realistic simulation of silica confined between walls requires an interaction potential that yields a reliable description of both the bulk properties and the properties at the interface between the walls and the silica fluid. As far as the bulk properties are concerned recent simulations have shown that the potential that was proposed some years ago by van Beest et al.  (the so–called BKS potential) is able to reproduce various properties of amorphous silica very well, such as its structure, its vibrational and relaxational dynamics, the static specific heat below the glass transition temperature, and the conduction of heat . Recent extensive simulations  that used a combination of a classical molecular dynamics simulation with the BKS potential and ab initio simulations (Car–Parrinello molecular dynamics) gave evidence that the BKS potential yields also a fair description of the structure of free amorphous silica surfaces, and thus, it can be expected that it also reproduces the main features of the structure of amorphous silica near a wall.

The functional form of the BKS potential is given by

$$\varphi _{\alpha \beta }(r)=\frac{q_\alpha q_\beta e^2}{r}+A_{\alpha \beta }\mathrm{exp}\left(B_{\alpha \beta }r\right)\frac{C_{\alpha \beta }}{r^6}\alpha ,\beta [\mathrm{Si},\mathrm{O}],$$
(1)
where $`r`$ is the distance between the ions of type $`\alpha `$ and $`\beta `$. The values of the constants $`q_\alpha ,q_\beta ,A_{\alpha \beta },B_{\alpha \beta }`$, and $`C_{\alpha \beta }`$ can be found in Ref. . For the sake of computational efficiency the short range part of the potential was truncated and shifted at 5.5 Å. This truncation has also the benefit to improve the agreement between the density of the amorphous glass at low temperatures as determined from the simulation with the experimental value. For the calculation of the long–ranged Coulomb interactions we used the method first introduced by Parry  which is a straightforward generalization of the three–dimensional Ewald sums with periodic boundary conditions in all three directions. In the latter case the Fourier sum requires only a loop over the $`N`$ particles of the system which leads, together with the real space sum, to an effective computational load that scales as $`N^{1.5}`$ . In contrast to that in the Parry–Ewald sums for systems with a quasi–two–dimensional geometry the Fourier and the real space part have to be computed by a double loop over all pairs of particles. However, we preferred to use Parry–Ewald sums rather than alternative faster methods  since it provides accurate results and it is easy to implement.

The walls were not constructed as to model a particular material but rather to be a generic surface that can be simulated conveniently. Each wall hence consisted of 563 point particles forming a rigid face–centered cubic lattice with a nearest–neighbor distance of $`2.33`$ Å. These point particles interact with the atoms in the fluid according to a 12–10 potential,

$$v(r)=4ϵ\left[(\sigma /r)^{12}(\sigma /r)^{10}\right],$$
(2)
with $`\sigma =2.1`$ Å, $`ϵ=1.25`$ eV, $`r`$ being the distance between a wall particle and a Si or O atom. Note that we choose the same potential between the wall atoms and both Si and O because we wish to simulate the generic effects of confinement, separated from possible additional effects due to surface enrichment of one component caused by the difference in surface forces, which may be a problem in many real systems. The great advantage of simulation is that by idealizing certain aspects of a problem a much better physical understanding of consequences of particular interactions may be achieved. In this spirit, a wall potential has been “invented” that keeps the slip small (when we apply a constant force in $`x`$–direction to the particles causing a corresponding flow) and provides at the same time a weak bonding of the wall particles to those of the melt such that no covalent bonds of wall particles with Si or O atoms are formed. In order fulfil these two “requirements” the $`r^{10}`$–term in Eq. (2) has turned out to be more suitable than a $`r^6`$–term.

The simulation box had linear dimensions $`L_x=L_y=23.066`$ Å in the directions parallel to the walls (in which also periodic boundary conditions were applied), and $`L_z=31.5`$ Å in the direction perpendicular to the walls. Thus, $`N=1152`$ atoms (384 Si atoms and 768 O atoms) were contained in the system to maintain a density around $`2.3`$ g$`/`$cm3 which is close to the experimental one at zero pressure . In the Parry–Ewald sums the parameter $`\alpha `$ and the cutoff wave–vector for the Fourier part were chosen to $`0.265`$ and $`6`$, respectively . The equations of motion were integrated with the velocity form of the Verlet algorithm with a time step of 1.6 fs. All runs were done in the $`NVT`$ ensemble whereby the temperature was kept constant by coupling the fluid to a Nosé–Hoover thermostat . Whenever an external gravitational–like field was switched on, the thermostat was only applied in the two directions perpendicular to the flow, i.e. in $`y`$– and $`z`$–direction, whereas otherwise the thermostat was used in all three directions. Bocquet and Barrat  also embarked on this strategy to provide a correct adjustment of the temperature when the system is sheared.

We investigated the three temperatures $`T=5200`$ K, $`T=4300`$ K, and $`T=3760`$ K at which we first fully equilibrated the system for $`29`$ ps, $`65`$ ps, and 122 ps, respectively. At $`T=5200`$ K and $`T=3760`$ K we continued with additional runs over 164 ps and 490 ps, respectively, from which we analyzed the equilibrium structure. Then we switched on a gravitational acceleration field of strength $`a_\mathrm{e}=9.6`$ Å$`/`$ps2 that was coupled to the mass of the particles. We mention that although the masses of oxygen and silicon differ by about a factor of two ($`m_\mathrm{O}=15.9994`$ u, $`m_{\mathrm{Si}}=28.086`$ u) we obtained identical velocity profiles for both species (see below), and this is due to the fact that the chemical bonding force between oxygen and silicon atoms is much stronger than the imposed acceleration field. With the acceleration field runs were made over $`736`$ ps, $`1.23`$ ns, and $`3.27`$ ns at $`T=5200`$ K, 4300 K, and 3760 K, respectively. In addition at $`T=5200`$ K we did a run over $`1.72`$ ns with field strength $`a_\mathrm{e}=3.8`$ Å$`/`$ps2. During these runs we stored the positions and velocities of the particles every $`0.16`$ ps (i.e. every 100 time steps). For the calculation of the steady state properties we used only those data that were obtained after two times the time that was used before for the equilibration. The total amount of computer time spent for these simulations was 16 years of single processor time on a Cray T3E.

## III Results

In this section we present the simulation results for the silica melt between confining walls. We characterize the structure of the melt and show how it changes if one switches on a gravitational field, i.e., if a Poisseuille flow is present in the melt. By analyzing the velocity profiles we discuss to what extent their behavior can be understood within a hydrodynamic description.

The snapshot, Fig. 1, illustrates on what length scales the walls affect the structure of the melt. It shows a part of a configuration at $`T=3760`$ K, namely the part that is within a distance less than $`8`$ Å away from one of the walls. One can clearly identify a layering of SiO4 tetrahedra which becomes less pronounced if one moves away from the wall. Note that at $`T=3760`$ K $`95`$% of the silicon atoms are four–fold coordinated by oxygen atoms and 96.5% of the oxygen atoms are two–fold coordinated by silicon atoms (a silicon and an oxygen atom are defined as neighbors if their distance is less than $`2.35`$ Å corresponding to the first minimum of the pair correlation function for Si–O). The rest of the atoms form defects, i.e. silicon atoms that are three– or five–fold coordinated by oxygen atoms and oxygen atoms that are one– or three–fold coordinated by silicon atoms. Thus, since less than 5% of the atoms are involved in such defects it makes sense to characterize the structure also at the relatively high temperature $`T=3760`$ K as a disordered tetrahedral network where SiO4 tetrahedra are connected via bridging oxygens (i.e. every oxygen atom belongs to two tetrahedra).

The layering of the SiO4 network near the wall can be quantified by means of the density profile which is plotted in Fig. 2a across half of the film  for all atoms and for the oxygen and silicon atoms only. In contrast to typical density profiles in simple monoatomic liquids the oscillations of the total profile, which indicate the layering near the walls, do not have a regular character. But this irregular behavior can be understood if one looks at the partial profiles for oxygen and silicon. The explanation that emerges is as follows: As we mentioned before confined fluid SiO2 maintains the irregular network SiO4 tetrahedra as in the bulk. At the same time the tetrahedra have to fill the available space and have to satisfy the bonding energy to the wall. To this end, it is necessary to align the tetrahedra adjacent to the wall such that a two–dimensional plane forms which contains three out of the four oxygen atoms of a tetrahedron, as well as the silicon atom at their center (slightly further away from the wall), while the fourth oxygen atom of the tetrahedron has to be further away from the wall for geometrical reasons: this causes the second peak of the oxygen distribution. Thus, the walls have a tendency to “orient” the network of coupled SiO4 tetrahedra in the fluid locally. It is evident that the oscillations in the local density of both silicon and oxygen are rather regular, like a damped cosine function, but the wavelength and phase of both cosine functions are different: their superposition causes then the rather irregular layering structure of the SiO2 total density. We expect that similar effects also occur in many other associating molecular fluids confined between walls if the wall–fluid interaction is weak enough that it does not affect the chemical ordering in the fluid as it is the case in our system.

For $`z>8`$ Å the total density profile shows only small oscillations around a constant value of 2.3 g$`/`$cm3 which is an indication for bulk behavior. Thus, keeping in mind that the density profile is symmetric, the bulk in our system seems to extend in $`z`$ direction from about $`8`$ Å to $`23.5`$ Å, i.e. it has a width of about 15.5 Å. However, Fig. 2a shows that the deviations from a constant density in the latter region are not due to the statistics: The partial density profiles exhibit still oscillations also in the middle of the film that one would not expect in a pure isotropic bulk system. But these oscillations do not affect the structural quantities that we show below, i.e. the pair correlation functions and the distributions of angles and rings, since these quantities essentially do not differ from the corresponding ones of a pure bulk simulation at the same temperature. Therefore, we consider in the following the region defined by $`8`$ Å $`z23.5`$ Å as the bulk.

The behavior of the total density profiles at the temperatures $`T=5200`$ K and $`T=3760`$ K in equilibrium and with an external force with an acceleration of $`9.6`$ Å$`/`$ps2 can be seen in Fig. 2b. In equilibrium the effect of decreasing temperature on the oscillations near the wall is an increase of the peak heights that is accompanied with a smaller value of $`\rho (z)`$ at the minima between the peaks. Thus, the layering becomes more pronounced if one decreases the temperature. If one switches on the gravitational–like field the effect is similar to an increase of the temperature. In the bulk region the density profiles are not very sensitive to a variation of temperature and/or the presence of the external force. Within the accuracy of our data the same value of about $`2.3`$ g$`/`$cm3 is reached for all four cases under consideration. However, this does not exclude of course the possibility of a dramatic change in the local structure, but as we will see in the following, the local structure is not strongly affected by the considered external forces.

Quantitities that provide information on the microscopic interparticle distances are the partial pair correlation functions $`g_{\alpha \beta }(r)`$ which are proportional to the probability of finding an atom of type $`\alpha `$ at a distance $`r`$ from an atom of type $`\beta `$. In an isotropic system the functions $`g_{\alpha \beta }(r)`$ are normalized by a phase factor $`4\pi r^2`$ . However, for the case that boundaries are present this factor has to be modified in that it is then determined only by that part of the surface $`4\pi r^2`$ of a sphere around a particle that fits into the system. Thus, the normalization factor depends on the distance of the particles from the walls in $`z`$–direction. In order to identify the differences of $`g_{\alpha \beta }(r)`$ close to the wall from those of the bulk we define a bulk region for $`8`$ Å$`z23.5`$ Å and wall regions for distances less or equal $`3`$ Å from the wall (containing respectively the first oxygen and silicon layer, see Fig. 2a). The latter wall layer is denoted in the following by WL. The $`g_{\alpha \beta }(r)`$ are calculated only for particles within the bulk or the WL whereby the aforementioned normalization factors are determined such that also the boundaries of the bulk region and the WL that are located within the fluid film are treated as real boundaries. Fig. 3 shows $`g_{\alpha \beta }(r)`$ for the different correlations at $`T=3760`$ K, each at equilibrium in the bulk and in the WL, as well as in the WL under the acceleration field $`a_\mathrm{e}=9.6`$ Å/ps2. We have not included the curves for the bulk with an external field because they essentially coincide with the equilibrium curves. This means that for the field strength $`a_\mathrm{e}=9.6`$ Å/ps2 the structure as reflected in $`g_{\alpha \beta }(r)`$ is not affected in the bulk by the gravitational field. From this we can consider the latter field strength as a small disturbance, and thus we can expect that the system can be treated as a Newtonian fluid in the bulk. In contrast to this at the walls there are significant changes in $`g_{\alpha \beta }(r)`$ if one switches on the external field in that the peaks become broader and their heights decrease slightly. These changes are most pronounced in $`g_{\mathrm{OO}}(r)`$ which is reasonable since the very first layer at the wall is formed by oxygen atoms.

The differences in $`g_{\alpha \beta }(r)`$ at the WL and in the bulk reflect the different packing of the tetrahedral network in both regions. In $`g_{\mathrm{SiO}}(r)`$ the second and third peaks are shifted to smaller distances in the bulk as compared to the WL. Similar effects are even more pronounced in $`g_{\mathrm{OO}}(r)`$ and $`g_{\mathrm{SiSi}}(r)`$. In these functions one can see that the fcc lattice of the wall affects the structure on length scales of next–nearest Si–Si and O–O neighbors. In $`g_{\mathrm{OO}}(r)`$ for the WL the periodicity of the peaks is approximately given by the lattice constant $`a=2.33`$ Å of the fcc lattice that is formed by the wall atoms. This periodicity appears less pronounced in $`g_{\mathrm{SiSi}}(r)`$ since the Si atoms are not so close to the wall as the O atoms and have thus more freedom to arrange themselves. Moreover the first peak in $`g_{\mathrm{SiSi}}(r)`$ is at a slightly smaller distance for the WL than for the bulk and this is due to the denser packing of SiO4 tetrahedra near the wall.

A further step towards the characterization of the local structure is to consider correlations between triples of particles. Simple quantities to study such correlations are the distributions of angles between three neighboring atoms. Thereby, two atoms of type $`\alpha `$ and $`\beta `$ are defined as neighbors if their distance is smaller than the position of the first minimum in the corresponding $`g_{\alpha \beta }(r)`$. The locations of these minima are at 3.64 Å, 2.35 Å, and 3.21 Å for the Si–Si, Si–O, and O–O correlations, respectively. In Fig. 4a we show the angle distribution functions $`P_{\alpha \beta \gamma }(\theta )`$ for the Si–Si–Si, Si–O–Si, and the O–Si–O angles in the bulk at equilibrium and under shear and at the WL at equilibrium. The temperature in all cases is 3760 K. We see in the figure that in each of the $`P_{\alpha \beta \gamma }(\theta )`$ shown there are only small differences between the bulk curves for the equilibrium and the sheared case which again demonstrates that the structure changes only slightly in the bulk region if a field of strength $`a_\mathrm{e}=9.6`$ Å$`/`$ps2 is switched on. There are also only small differences in $`P_{\mathrm{OSiO}}`$ for the WL as compared to the corresponding distributions in the bulk. Since the O–Si–O angle is an intratetrahedral one this shows that the geometry of the SiO4 tetrahedra is essentially the same at the wall and in the bulk. Remarkable differences are in $`P_{\mathrm{SiSiSi}}`$ and $`P_{\mathrm{SiOSi}}`$ for the WL as compared to corresponding functions for the bulk: First, the main peak in $`P_{\mathrm{SiOSi}}`$ is shifted towards smaller angles which is due to the denser packing of the tetrahedra in the WL. Secondly, there is a shoulder in $`P_{\mathrm{SiOSi}}`$ around 100 which is less pronounced in the WL and $`P_{\mathrm{SiSiSi}}`$ exhibits a secondary peak at about 60 which is more pronounced in the curve for the WL. In Fig. 4b we compare the different $`P_{\alpha \beta \gamma }(\theta )`$ for the WL at equilibrium and for the sheared system. As can be inferred from the figure, in contrast to the bulk case there are changes in the angular distributions in that all three distributions broaden significantly for the system under shear. Moreover, the shoulder in $`P_{\mathrm{SiOSi}}`$ around 100 has about a six times larger amplitude in the latter case.

Recent studies of silica melts  have shown that the shoulder around 100 in $`P_{\mathrm{SiOSi}}`$ is due to two–membered rings whereas the peak at about 60 in $`P_{\mathrm{SiSiSi}}`$ corresponds to the presence of three–membered rings. We look now for the latter ring sizes in our system and investigate the distribution of rings in the network. A ring is defined as follows: One starts from any silicon atom and two of its nearest oxygen neighbors. Then one counts the shortest consecutive sequence of Si–O elements that connect the two latter oxygen atoms and this number is then the ring length $`n`$. Thereby, the ring length $`n=1`$ corresponds to a dangling bond, i.e., an O atom is only attached to one Si atom, and in a two–membered ring ($`n=2`$) two silicon atoms share two oxygens. Note that one can extract (indirectly) information on rings in SiO2 by means of NMR experiments .

Fig. 5a shows the ring distribution function $`P(n)`$ for $`T=3760`$ K in the bulk and in two different wall layers denoted by WL1 and WL2. WL1 and WL2 are defined as the regions which are respectively within a distance of $`6.25`$ Å and $`3.0`$ Å away from the wall corresponding to the second minimum of the density profile for silicon and the first minimum of the total density profile, respectively (see Fig. 2a). Note that WL2 is identical with WL. In each region, i.e. bulk, WL1, and WL2, we took only those rings into account that fit completely into it. Thus, in WL2 those rings are counted that are formed at each case by the first and the second O and Si layers (with respect to the distance from the wall), whereas with WL2 only those rings are taken into account that are formed by the first O and the first Si layer. This is justified because the first two oxygen and the first two silicon layers are well–defined in that the minima in the corresponding density profiles are close to zero density in the case of WL2 and around the small value $`\rho =0.5`$ g$`/`$cm3 in the case of WL1. Furthermore, we can infer from Fig. 2a that in contrast to the second layers the first oxygen layer overlaps strongly with the first silicon layer and the overall thickness of both layers is only about $`2`$ Å. Thus, the first oxygen and the first silicon layer form a quasi–two–dimensional plane and $`P(n)`$ for WL2 gives a distribution of rings that have an orientation parallel to the walls whereas in $`P(n)`$ for WL1 also the rings perpendicular to the walls are included.

As has been also found in pure bulk simulations of SiO2 , in the bulk a maximum is observed around $`n=6`$. This is plausible since in silica the high–temperature crystalline phase at zero pressure, $`\beta `$–cristobalite, exhibits only six–membered rings. In WL1 the probability for $`n6`$ is smaller than in the bulk in favor of a relatively high probability of $`n=3`$ and $`n=4`$. In WL2 it is the other way round: $`n=4`$ and even more $`n=5`$ are less frequent than in the bulk in favor of $`n=8,9,10`$. The ring structure near the walls that corresponds to these findings is as follows: Perpendicular to the walls small rings with $`n=3,4`$ are seen such that, e.g., $`n=3`$ is formed by two silicon atoms from the first silicon layer with a third one from the second silicon layer (for an illustration see Fig. 1). In contrast to that, parallel to the walls (considering the first oxygen and the first silicon layer) an open structure with relatively large rings is observed which compensates somewhat the dense packing of SiO4 tetrahedra perpendicular to the walls.

Fig. 5b shows the behavior of $`P(n)`$ in the bulk at the two temperatures $`T=3760`$ K and $`T=5200`$ K in equilibrium and under shear. We can immediately infer from the figure that the considered shear fields have only a small effect on the structure which confirms the findings for all the quantities discussed before. At $`T=5200`$ K one has a relatively large amount of two– and three–membered rings the frequency of which is more than a factor of two smaller at $`T=3760`$ K. In Ref.  it was shown that their frequency of occurrence decreases further with decreasing temperature such that the amount of two–membered rings falls far below 1% for systems that have typical structural relaxation times of the order of $`1`$ ns.

In contrast to the bulk in WL1 significant changes in the ring distribution take place if the system is sheared (see Fig. 5c), and the external force field affects the ring structure such that small and large rings are formed, while at the same time especially the amount of six–membered rings decreases. Only for the smaller field strength $`a_\mathrm{e}=3.8`$ Å$`/`$ps2 there are no significant changes in the ring distribution. We will see below that the latter is accompanied with a very small slip motion at the walls whereas a large slip velocity is correlated with strong rearrangements in the ring structure. One can also infer the remarkable fact from Fig. 5c that the probability to find rings with $`n=3,4`$ does not change very much when an external acceleration field is switched on. This is reasonable because these small–membered rings, as we have seen before, are located perpendicular to the walls and thus they are very stable to shear forces that are imposed parallel to the walls.

The strongest rearrangements in the ring structure due to a shear field are found when we consider the region WL2. The corresponding curves are shown in Fig. 5d. Again, there are only minor changes in $`P(n)`$ at $`T=5200`$ K and $`a_\mathrm{e}=3.8`$ Å$`/`$ps2 as compared to the corresponding equilibrium case which is, as mentioned before, related to the presence of only a very small slip velocity. For the higher acceleration field, $`a_\mathrm{e}=9.6`$ Å$`/`$ps2, the ring structure becomes more heterogeneous and the effect of the external field is if one would locally increase the temperature. The rearrangements in the ring structure can be summarized as follows: Rings mainly of size $`n=6,7,8`$ are broken under the influence of the shear force and instead small rings with $`n<4`$ and very large rings with $`n9`$ are formed. This is illustrated in Fig. 6 where snapshots at $`T=3760`$ K of the WL2 region are shown for the equilibrium case and the non–equilibrium steady state case with the external field of strength $`a_\mathrm{e}=9.6`$ Å$`/`$ps2. One can clearly identify, e.g., rings with $`n=2`$ and $`n=11`$ in the right panel which are absent in the left panel, i.e. at equilibrium.

Up to now we have seen that in the bulk region the structure essentially does not change for acceleration fields $`a_\mathrm{e}9.6`$ Å$`/`$ps2 whereas significant structural rearrangements occur close to the walls for $`a_\mathrm{e}=9.6`$ Å$`/`$ps2. We demonstrate now how this is reflected in the flow properties, and we discuss to what extent the flow can be understood within a hydrodynamic description.

Fig. 7 shows velocity profiles $`v_x(z)`$ across the film for two different choices of the acceleration in $`x`$–direction for $`T=5200`$ K. The solid lines are fits according to the formula

$$v_x(z)=\frac{\rho a_\mathrm{e}}{2\eta }(zz_1)(L_zz_1z)$$
(3)
with $`\rho `$ being the total mass density, and the intercept $`z_1`$ as well as the shear viscosity $`\eta `$ are treated as adjustable parameters. Eq. (3) is the formula for a one–component fluid, but it is justified to use this equation for our two–component system since the velocity profiles for the individual species, i.e. for silicon and oxygen, are identical within the statistical errors. Note that this was also observed by Koplik and Banavar  for a mixture of two simple liquids confined between two walls. The values of $`\eta `$ that result from the fits with Eq. (3) are quoted in the figure. As it should be, $`\eta `$ is independent of the acceleration $`a_\mathrm{e}`$ within statistical errors. This matches well to the aforementioned result that the structure is essentially not affected in the bulk by the considered shear fields. Furthermore, this analysis shows that hydrodynamics holds down to nanoscopic scales, even if we use a realistic model of a molecular glassforming fluid such as SiO2. Close to the wall the fit with Eq. (3) is obviously not valid for $`a_\mathrm{e}=9.6`$ Å$`/`$ps2 and we find a slip motion of the fluid. As we have seen before this slip is accompanied with structural rearrangements in the ring distribution.

It is conspicuous in Fig. 7 that there is a sharp peak in $`v_x(z)`$ at $`z=1.2`$ Å in the curve for $`a_\mathrm{e}=3.8`$ Å$`/`$ps2 (there is exactly the same peak at $`z=30.3`$ Å because the profile is symmetrized). This is just due to the fact that there is a very small probability for the particles to move into the region $`z<1.5`$ Å, i.e. the region very close to the wall. And if a particle enters this region there is a high probability that it has a relatively high velocity. Of course, on average this velocity should be zero, but in this case the statistics is bad since this event is very rare, and we observe sharp peaks in the regions very close to the walls that can have a negative or positive amplitude with equal probability (note that they appear also in the unsheared system). A similar reasoning might be also true at $`z=3`$ Å, i.e. at the location of the first minimum in the total density profile. Indeed we see a small peak at $`z=3`$ Å in the curve for $`a_\mathrm{e}=9.6`$ Å$`/`$ps2 which we also find at lower temperatures (see Fig. 8). But we cannot conclude from our data whether this feature would disappear with a better statistics.

One still could argue that the fits with Eq. (3) are fortuitous, and therefore the resulting shear viscosity $`\eta `$ then simply is a fit parameter with little physical significance, different from the actual shear viscosity of bulk molten SiO2. Fig. 8 shows that this objection is not true: here the velocity profiles are compared at three temperatures, and the viscosities $`\eta `$ from the fits are compared to viscosities $`\eta _{\mathrm{GK}}`$ obtained earlier from simulations of bulk SiO2 at rest (no walls, no flow) applying Green–Kubo formulas . Note that the latter approach typically involves a relative error of about 20%, and so does the present approach. Within these errors, there is a good agreement, although $`\eta _{\mathrm{GK}}`$ changes for the shown range of temperatures by an order of magnitude. However, the estimated $`\eta `$ seems to be systematically smaller than $`\eta _{\mathrm{GK}}`$ at low temperatures and we cannot exclude that these deviations are due to the influence of the Nosé–Hoover thermostat on the flow velocity. Of course, we can also infer from Fig. 8 that we cannot use our simulation of Poisseuille flow to estimate the viscosity at much lower temperatures than $`3760`$ K: If $`\eta `$ decreases the velocity profile becomes flatter, and one cannot compensate this by higher values of $`a_\mathrm{e}`$ since the slip velocity increases with $`a_\mathrm{e}`$ also.

## IV Summary and discussion

We used molecular dynamics computer simulations to study a silica melt confined between two walls. The interactions between the silicon and oxygen atoms were modelled by the so–called BKS potential for which it has been shown that it describes reliably the bulk as well as free surface properties of silica. The walls consist of particles fixed on a fcc lattice, and wall–fluid interactions were modelled by a simple 12–10 potential. The aim of this work was to investigate how the walls affect the structure of the melt in equilibrium and under shear, and whether it is possible to understand the flow properties of the sheared melt within hydrodynamics.

As in the case of simple liquids a pronounced layering of the melt close to the walls is observed. This layering, as seen in the total density profile, has a very irregular character and this is related to the specific chemical ordering of silicon and oxygen atoms in a tetrahedral network. In the layers that are closest to the walls the tetrahedra prefer an orientation such that three oxygen atoms of each tetrahedron form the first layer, the silicon atoms the second one and the remaining fourth oxygen atoms of the tetrahedra the third one (see Fig. 1).

It has been demonstrated in recent simulation studies of free silica surfaces  that the ring distribution function is a quantity that is well suited to characterize the differences of the surface structure on intermediate length scales from that in the bulk. This holds also true in the present case: The relatively high probability of small and large rings makes the difference of the melt structure at the wall from that of the bulk. Perpendicular to the first wall layer there is a high probability of rings of size $`n=3,4`$ whereby, e.g., a ring with $`n=3`$ is typically formed by two silicon atoms from the first silicon layer and one from the second one (of course, these rings are also decorated with oxygen atoms). Parallel to the walls, at least for temperatures $`T3760`$ K a relatively high probability for rings of length $`n9`$ is found in the wall layers which can be seen as the complement to the occurrence of small–membered rings with $`n=3,4`$. Our findings for the structure of the silica melt near the wall can be expected to be generic if the wall–SiO2 interactions are, as in our case, relatively weak, i.e. no covalent bonds of the wall atoms with those of the fluid are formed as it would be the case, e.g., for an interface between crystalline and amorphous silica.

Under shear, the bulk structure of the silica melt does not change for the field strengths considered in this work, $`a_\mathrm{e}9.6`$ Å$`/`$ps2. This is an indication that it is justified to treat the melt as a Newtonian fluid. At the walls, the behavior of the fluid is different: If a slip motion occurs there are strong rearrangements in the ring distribution as compared to the equilibrium case. Rings of medium size break and, as a result of this, small and large rings are formed (the same effect is achieved by a local increase of the temperature). Thus, we see a rather unusual mechanism for slip motion: Slip seems to occur only if the external (gravitational–like) force field is strong enough to break rings of medium length. This slip mechanism can be expected to be generic for liquids with a tetrahedral order (of course, as before, only if walls are present that do not have a (strong) covalent bonding with the particles in the tetrahedral network).

As expected from hydrodynamic theory of a Newtonian fluid we find parabolic velocity profiles in the bulk region, although the latter region extends in the systems considered in this work only over about 15 Å. The values of the shear viscosity $`\eta `$ that we have extracted from the profiles are in fair agreement with the “Green–Kubo” values from a bulk simulation. So the question arises whether the non–equilibrium molecular dynamics technique used in this work offers a method to estimate transport coefficients. It has been demonstrated for a lattice gas model  that non–equilibrium simulation techniques can be very powerful to gain insight into the understanding of transport processes in binary mixtures. We have shown in this work that with the present simulations this is possible for melts with a viscosity up to about 1 Poise. For systems with higher viscosities (at lower temperatures) a different wall model would be helpful with which the tendency for slip motion is smaller and furthermore, larger system sizes would be required. But for larger system sizes one needs for the computation of the quasi–two–dimensional Coulomb forces faster methods than the Parry–Ewald method used in this work. A promising new method  (also based on Ewald sums) that scales with the number of particles like $`N\mathrm{log}N`$ has been proposed very recently and could be applied to systems like the one considered in this work.

Acknowledgments: The authors are grateful to W. Kob for very valuable discussions and ideas. We also thank F. Varnik and M. Müser for constructive hints and C. Mischler for a critical reading of the manuscript. Furthermore we acknowledge partial financial support from the D.I.P. Project 352–101 and we thank the HLRZ Stuttgart for a generous grant of computer time on the Cray T3E.


format: text

text:
# A lattice of Magneto-Optical and Magnetic traps for cold atoms

## I Introduction

An array of individual microtraps for ultracold atoms, where each trap can be addressed separately, provides a large number of new degrees of freedom for coherent manipulation and control of cold atoms. Besides applications in Quantum Computing, where every lattice site corresponds to a qubit, an addressable array will allow the dynamic definition of quantum dots and quantum wires to perform e.g. transport measurements on fully controllable structures filled with quantum gases.
One possible way to reach this goal could be microstructured traps. Since the first proposals for such traps for neutral atoms Weinstein  a lot of progress has been made in this field Reichel1 ; H"ansel ; Cassettari , including the generation of Bose Einstein Condensates (BEC) in microtraps Ott ; H"ansel1 . The main advantages of miniaturized traps are on the one hand the possibility to create nearly arbitrary potentials and on the other hand the possibility to reach already with small currents traps with much tighter confinement than in large macroscopic traps. Reviews on the basic principles of such microtraps can be found in references Schmiedmayer2 ; Reichel .

In this paper, we present basic periodic trapping configurations which can be realized with multilayer microstructures. The trapping fields for the atoms (Ioffe Pritchard type potentials for magnetic trapping of atoms as well as quadrupole fields for magneto optical trapping) are formed by two layers of crossed wires, which can be individually addressed. We present first experiments with a system of macroscopic wires and show, that it is possible to produce multiple magneto optical traps (MOT) next to a surface in a controlled manner. This could be a first step to load ultracold atoms into an array of microstructured magnetic traps in a parallel way. We envisage that experiments with trapped atoms or even the production of BECs in such arrays of traps can be performed.

The paper is organized as follows: In section II.1 a review of the basic concepts of trapping atoms with a single wire, respectively with microstructures is given. In part II.2 we present a scalable implementation of trapping configuration for cold atoms in periodic arrangements. In section III the experimental setup is introduced and finally, in section IV we present first experiments including multiple MOT-systems based on such multilayer structures.

## II Configurations

### II.1 Basic configurations

The simplest configuration to produce magnetic potentials which can be used to catch cold atoms is the superposition of the magnetic field of an infinitely long wire $`B=|\stackrel{}{B}|=\frac{\mu _0}{2\pi }\frac{I}{r}`$ ($`I`$ is the current in the wire, $`r`$ the distance from the wire) with an external homogenous bias field $`B_b=|\stackrel{}{B}_b|`$ perpendicular to the wire Frisch  (see fig. 1 a). This configuration yields a two dimensional quadrupole field parallel to the wire with a field minimum at a distance of $`r_m=\frac{\mu _0}{2\pi }\frac{I}{B_b}`$. Next to this minimum the magnetic field increases linearly in all directions perpendicular to the wire with a gradient of $`\frac{dB}{dr}|_{r_m}=\frac{2\pi }{\mu _0}\frac{B_b^2}{I}`$ Reichel1 . This type of trap can be extended to a harmonic trap by applying an external field $`B_p=|\stackrel{}{B}_p|`$ parallel to the wire. The resulting curvature in the radial direction at the height $`r_m`$ can be expressed as

$$\frac{d^2B}{dr^2}|_{r_m}=\left(\frac{2\pi }{\mu _0}\right)^2\frac{B_b^4}{B_pI^2}$$
(1)
where I is the current in the wire Schmiedmayer1 . For this one wire configuration, the bias field $`B_b`$ can alternatively be produced by two additional wires on both sides of the trapping wire with currents in opposite direction Thywissen . A similar trapping configuration consisting of four wires in parallel is possible too Thywissen .

The 2-dimensional trapping potential of infinitely long wires can be extended to real 3-D trapping configurations by closing the quadrupoles in the third dimension by adding two wires perpendicular to the main wire. This can be done with one wire by bending it in so called ”U ”or ”Z” configurations (for illustration see fig. 1b,c) Reichel1 , resulting in 3-D quadrupoles (U-configuration) or Ioffe-Pritchard type traps (Z-configuration) respectively.

A 3-D magnetic quadrupole has three main axes which go through the magnetic zero and are orthogonal with respect to each other. To satisfy $`\stackrel{}{}\stackrel{}{B}=0`$, there is one axis where the sign of the field gradient must be different than for the other two. For the further explanation we use this axis to define the orientation of the 3-D quadrupole in space. Furthermore we define the bent parts of the wire to be parallel to the x-direction, where the central part is in the y-direction. The z-direction is perpendicular to the plane defined by the x-and y-axes.

In the ”U” configuration Reichel1  the y-components of the magnetic field of the bent wires cancel out at the origin of the linear quadrupole (see fig. 1b), whereas the z-components add. Due to this the origin of the quadrupole is shifted in the +x direction away from central bar. The gradients in both the x-direction and the y-direction are

$`{\displaystyle \frac{dB}{dx}}|_{r_m}`$ $`=`$ $`{\displaystyle \frac{\mu _0}{2\pi }}{\displaystyle \frac{I_w}{r_m^2}}`$ (2)
$`{\displaystyle \frac{dB}{dy}}|_{r_m}`$ $`=`$ $`{\displaystyle \frac{2\mu _0}{\pi }}{\displaystyle \frac{I_cr_mL}{(L^2+4r_m^2)^2}},`$ (3)
where $`r_m`$ is the distance of the quadrupole minimum from the wire, L is the length of the central bar of the U-trap and $`I_w=I_c=I`$ is the current in the wire. In this formula the shift due to the bent parts is not taken into account. This simplification is valid, if $`r_m`$ is small compared to the spacing $`L`$ between the bent parts of the wire. The axis of the quadrupole in the configuration shown in fig. 1b is lying in the x-z-plane above the central bar between the bent parts of the wire under an angle of $`45^{}`$ with respect to the x-axis. The bent parts of the ”U” configuration do not close the quadrupole in a symmetric way. A more symmetric configuration is a ”H” configuration, which can not be bent in a simple way with a single wire (see inset fig. 1b).

In the ”Z”-configuration Reichel1  a Ioffe-Pritchard-type trap is formed (see fig. 1c). In this configuration the y-components of the magnetic field of the bent parts add and the z-components cancel out in the middle of the central bar. For this reason a magnetic field minimum appears in the middle of the central bar at the distance $`r_m`$. Due to the non vanishing y-components the magnetic field at the field minimum is non zero and can be calculated as

$$B_m=\frac{2\mu _0}{\pi }\frac{I_cr_m}{L^2+4r_m^2},$$
(4)
where L is the length of the central bar and $`I`$ is the current in the wire. The curvature of these traps at the minimum $`r_m`$ are in both the x-direction and the y-direction

$`{\displaystyle \frac{d^2B}{dx^2}}|_{r_m}`$ $`=`$ $`{\displaystyle \frac{\mu _0}{16\pi }}{\displaystyle \frac{I_w(L^2+4r_m^2)}{r_m^5}}`$ (5)
$`{\displaystyle \frac{d^2B}{dy^2}}|_{r_m}`$ $`=`$ $`{\displaystyle \frac{2\mu _0}{\pi }}{\displaystyle \frac{I_c(16r_m^416r_m^2L^2L^4)}{r_m(4r_m^2+L^2)^3}}`$ (6)
This kind of trap is more stable against Majorana transitions than the quadrupole trap and is thus favourable for magnetic trapping of atoms, whereas ”U” type traps are favourable for magneto-optical trapping of atoms, which typically requires a quadrupole configuration. This ”Z” trap is again asymmetric and can be extended to an ”H” form trap, which is a more symmetric situation (see inset fig. 1c).

### II.2 Lattice configurations

To extend these basic schemes to an array of traps, let us consider 2 layers of parallel wires, where the 2 layers are perpendicular to each other (see fig. 2 and fig. 3). For simplification let us assume that all wires carry the same constant current I, the wires are infinitely thin and infinitely long. In addition to the magnetic field produced by these wires an external bias field $`B_b`$ can be applied parallel or perpendicular to the 2 layers. Such a configuration is very versatile and allows to create many different lattice geometries. The possible lattice configurations are built up from simple unit cells, which will now be described in detail.

We consider three different configurations. In the first configuration, the currents $`I`$ in each of the layers are unidirectional (see fig. 2a, configuration A). In the second configuration the currents are unidirectional in one layer (y-direction), and alternating in the second layer (x-direction) (see fig. 2b, configuration B). In the last configuration the currents are bidirectional in both layers (see fig. 3a, configuration C).

Configuration A (currents in both layers unidirectional) has already been discussed in references Reichel2 ; Yin . We describe the concept using simple unit cells consisting of 3 wires segments crossing each other as marked in the box in fig. 2a. This unit cell has a central bar in the middle, enclosed by two wire segments in the upper and lower part of the unit cell. The size of the unit cell is given by the distance between the wires in each layer. The lattice is a simple square lattice. By applying an external bias field in the x-direction as indicated in fig. 2a each ”H” shaped unit cell contains a Ioffe Pritchard trap above the central bar. The basic properties of an individual trap can be calculated from the equations 5 and 6 given above for the ”Z”-trap plus a perturbation term, which will be calculated later.

In configuration B, the currents are unidirectional in one layer (y) and alternating in the second layer (x) (see fig 2 b). Here the simple unit cell consists of 4 wire segments as shown in fig. 2 b). The lattice is rectangular. By applying an external bias field in the x-direction two different quadrupole traps $`Q_{}`$,$`Q_+`$, (see fig. 2b) are contained above the surface in the unit cell, which can be described by the equations 2 and 3 given for a quadrupole trap plus a perturbation term. The difference between the two quadrupoles in the unit cell is the orientation. At quadrupole $`Q_{}`$ the axis is lying in the x,z plane under an angle of $`45^{}`$ with respect to the x-axis, and at quadrupole $`Q_+`$ the axis is under an angle of $`45^{}`$ with respect to the x-axis. The quadrupoles can be used for magneto optical trapping to catch atoms in a reflection MOT Reichel2 ; Lee ; Schneble . Due to the fact that the two quadrupoles have different orientations, the circular polarisation of the MOT beams has to be different to either use the upper or the lower of the quadrupoles as will be shown in the experimental section.

In configuration C, the currents are alternating in both layers (fig. 3a). The unit cell consists of 8 wire segments. Even without an external bias field in this cell 10 different quadrupoles can be found. The bias field for all the different quadrupoles is given by the neighboring wires as mentioned before Thywissen , so no external bias field is necessary. Eight of the quadrupoles can be found above the wires ($`R_1`$..$`R_4`$ and $`L_1`$..$`L_4`$, see fig. 3). The quadrupole axes of all of them are tilted by $`45^o`$ against the surface with respect to the wires and can be used for a reflection MOT Reichel2 ; Lee ; Schneble . Two other quadrupoles ($`P_1`$,$`P_2`$) are built up by two concentric wire squares Thywissen . The quadrupole axes of these two are oriented perpendicular to the surface, and thus they can not be used for a reflection MOT. The unit cell itself is symmetric under a mirror operation around the mirror plane as indicated in fig. 3 b). The lattice is a simple square lattice.

In all three configurations, the magnetic field $`\stackrel{}{B}`$ in the unit cell is given by the magnetic field produced by the wires in the unit cell $`\stackrel{}{B}_c`$ plus a perturbation term from the rest of the lattice $`\stackrel{}{B}_p`$

$$\stackrel{}{B}=\stackrel{}{B}_c+\stackrel{}{B}_p$$
(7)
For configuration A the perturbation in one cell can be expressed by the formula

$`\stackrel{}{B}_p`$ $`=`$ $`{\displaystyle \frac{\mu _0I}{2\pi }}{\displaystyle \underset{\begin{array}{c}n=\mathrm{}\\ n0\end{array}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(xnd)^2+z^2}}\left(\begin{array}{c}xnd\\ 0\\ z\end{array}\right)`$ (13)
$`+`$ $`{\displaystyle \frac{\mu _0I}{2\pi }}{\displaystyle \underset{\begin{array}{c}m=\mathrm{}\\ m0\end{array}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(ymd)^2+z^2}}\left(\begin{array}{c}0\\ ymd\\ z\end{array}\right)`$ (19)
and for configuration B

$`\stackrel{}{B}_p`$ $`=`$ $`{\displaystyle \frac{\mu _0I}{2\pi }}{\displaystyle \underset{\begin{array}{c}n=\mathrm{}\\ n0\end{array}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{(xnd)^2+z^2}}\left(\begin{array}{c}xnd\\ 0\\ z\end{array}\right)`$ (25)
$`+`$ $`{\displaystyle \frac{\mu _0I}{2\pi }}{\displaystyle \underset{\begin{array}{c}m=\mathrm{}\\ m0\end{array}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(ymd)^2+z^2}}\left(\begin{array}{c}0\\ ymd\\ z\end{array}\right)`$ (31)
where d is the distance between the wires and I the current flowing in the wires. The sums over n and m runs over all wires in the two horizontal directions. The terms for $`n=m=0`$ are excluded from the sum since only the perturbation from the rest of the lattice on a given unit cell should be considered here. An analogous situation is given in the third configuration. The perturbations to the magnetic field in one lattice cell will be small, if $`d`$ is much larger than the distance $`z=r_m`$ of the potential minimum from the surface.

For finite systems and $`r_md`$ the sum above is finite and one has to consider edge effects, which perturb the symmetry. This is analogous to finite-size magnetic mirrors used in atom optics Lau . The effect can be compensated by adding compensation wires, as proposed in Zabov .

## III Experimental setup

To test the configurations proposed in section II.2 arrays of magnetic quadrupole traps are created by 2 layers of macroscopic wires mounted inside a vacuum chamber (see fig. 4). Each layer of wires consists of a ribbon UHV cable with a pitch of 1.27 mm (Caburn, Kapton Ribble Wire). The individual wires are kapton insulated with a diameter of the copper core of 0.25 mm. The cables are glued together with UHV compatible epoxy and are mounted on top of each other in a two layers configuration, where the upper layer consists of 10 wires in parallel and the lower layer consists of 8 parallel wires below the first one. The distance between the center of the wires in the two layers is 1.27 mm, given by the wire diameter. The whole wire assembly is attached to a Cu holder. The maximum current in the wires was about 5A. At higher currents, the pressure in the vacuum chamber started to increase probably due to outgasing of the insulating kapton layer. On top of the upper layer of wires a Si plate (thickness: 0.4 mm) covered with 10nm of Cr and 100 nm of Au is mounted, which serves as mirror. The whole setup is mounted upside down inside the chamber. The mirror at the surface is necessary, because we use the reflection MOT principle Reichel2 ; Lee ; Schneble  to catch cold $`{}_{}{}^{87}Rb`$ atoms. For this type of MOT two of the laser beams are parallel to the surface and two other lasers beams are reflected from the surface under an angle of $`45^{}`$. The magnetic quadrupole field necessary for the MOT is generated by the two layers of wires plus an external bias field. This bias field can be adjusted between 0 G and 80 G. The circular polarisation of the four beams depends on the orientation of the magnetic quadrupole. The $`{}_{}{}^{87}Rb`$ is provided by a Rb dispenser SAES ; Fortagh  mounted approximately 4 cm away from the assembly. The optical access is provided by a large window below the sample.

To operate the MOT, we use as lasers two home built grating feedback stabilized diode lasers, operating at 780 nm. The cooling laser is locked on the $`5S_{1/2},F=25P_{3/2},F^{}=3`$ cycling transition using polarisation spectroscopy in a vapor cell. The laser power in this laser is about 17 mW. The repumping laser is locked to the $`5S_{1/2},F=15P_{3/2},F^{}=2`$ transition using the same locking scheme as for the MOT laser. This laser has a total power of 7 mW.

To detect the atoms in the MOT, we used a calibrated 8-bit CCD Camera, placed 8 cm under the surface outside the vacuum chamber. The imaging system is set up in a way, that atoms above the whole surface inside the chamber can be detected. One pixel on the camera corresponds to a $`39\mu `$m$`\times `$$`36\mu m`$ area in the trap region.

## IV Experimental results

By using one wire in the upper one of the two layers together with the external bias field and turning on all the MOT beams one has the situation of a 2D MOT with additional optical molasses on the axis parallel to the surface (see fig. 5), upper left corner). This 2D-MOT can be extended to a 3D-MOT by using the wires of the lower layer transform the 2-D quadrupole to an 3-D quadrupole. In the experiments we used in the lower layer the two outermost wires of our experimental setup, which are at a distance of 10.1 mm (see fig. 5).

We found that the position of the MOT can be precisely controlled by choosing the current carrying wire in the upper layer. The currents were set to 4 A and an external bias field of 3 G was applied. With these values the magnetic field gradients be calculated from equations 2 and 3 as $`\frac{dB}{dx}=11G/cm`$ and $`\frac{dB}{dy}=1.9G/cm`$. The calculated distance $`r_m`$ from the MOT to the surface is 1.6 mm.

In this configuration, a MOT was realized above every wire (see fig. 5). The average number of atoms for the different MOT positions was $`310^5`$. A situation where no external bias field was necessary was reached by using as bias field wires the neighboring wires on the left and on the right of the MOT wire in the upper layer with current flowing in the opposite direction. The number of atoms in such a MOT generated without external field was lower, due to the fact that the capture range was smaller. The reason for this is the fact that the bias field produced by the two neighboring wires is inhomogeneous. The quadrupoles resulting from this configuration extend over a smaller volume and so the capture range is reduced.

In the proposed configuration B two quadrupoles which are different in orientation emerge above the middle wire segment of the unit cell. These two quadrupoles can be used to operate the MOT, if the respective polarisations of the MOT laser beams fit to the orientation of the quadrupoles. In our experimental setup, we realized the unit cell configuration consisting of 4 wires (see fig. 6a). The distance between the parallel wires was 5mm (upper - middle wire), respectively 3.8 mm (middle - lower wire). The current in each of the wires was chosen as 4 A and an external bias field of 3 G was applied. This yields a gradient of $`\frac{dB}{dx}=11G/cm`$ in the x- and $`\frac{dB}{dy}=3.3G/cm`$ in the y-direction. The distance $`r_m`$ above the surface is again $`1.6mm`$. By changing the polarisation of the MOT beams parallel to the surface from $`\sigma ^+`$ to $`\sigma ^{}`$ the MOT changes its position due to the two different orientations of the two quadrupoles within the rectangular unit cell (see fig. 6 b).

In configuration C, 10 quadrupoles are created per unit cell without external fields. As explained before 8 of these quadrupoles are in principle suited to operate a reflection MOT, because the angle between the quadrupole axis and surface is $`45^{}`$. In our experimental configuration only four of them can be seen due to the fact that only for four quadrupoles the axes coincide with the direction of a MOT beam. The accessible positions can again be addressed by changing the polarisation of the MOT beams, where the polarisation configuration is determined by the orientation of the magnetic quadrupoles for the different wire configurations (see fig. 7a). For the generation of the magnetic field in the unit cell we used 3 wires in the y-direction with alternating currents (distance wire-wire 2.5mm) and 3 wires in the x-direction (distance between the wires 5 mm, respectively 3.8 mm, as in configuration B). By building up only one such lattice cell, the quadrupoles at the edges of the lattice cells are not accessible, because only in the middle of the cell the bias field needed is produced by the neighboring wires. So we investigated in this configuration the quadrupoles in the middle of the unit cell ($`R_2`$,$`L_4`$). Then we switched off the current in the left wire, and turned on the current in the right wire to be able to observe the quadrupoles $`R_4`$ and $`L_2`$ (see configuration marked in fig. 7b-e). To reach a stable situation the current in the two middle wires was in both cases 3.5 A (y-direction) and in the outermost wires 5A. These values yield a distance of the quadrupole from the surface of 0.2 mm and a gradient $`\frac{dB}{dx}=20G/cm`$ in the x-direction. The currents in the x direction were 4A. This yields a gradient of $`\frac{dB}{dy}=9.7G/cm`$ in the y-direction. As shown in fig. 7b-e the different quadrupoles can be addressed by changing the polarisations of the MOT beams as expected. By changing the circular polarisation of the beams parallel to the surface we can switch between $`L_4`$ and $`R_2`$ (see fig. 7b,c) and $`L_2`$ and $`R_4`$ (see. fig. 7d,e) respectively. The distance change in the y-direction is 4.4 mm respectively 4 mm and in the x-direction 1 mm. By changing the circular polarisation of the $`45^{}`$ beams (see fig. 7b,d, c,e) one can switch between $`L_4`$ and $`L_2`$, respectively between $`R_2`$ and $`R_4`$. Here the observed change in position in x direction is 2.2 mm. The number of atoms in the different MOTs is between $`810^4`$ and $`2.410^5`$ atoms. The main advantage of this configuration compared to configuration B is the fact that no external bias field in necessary.

Up to now only single unit cells were observed. To built up lattices it is necessary to create arrays of unit cells. This can be done in the experiment by using more wires. For an array of MOTs, we arranged 4 unit cells in configuration B. As explained before it is necessary that the distance between the wires has to be large compared to the distance of the field minimum $`r_m`$ to the wire. To reach this, we used two wires in the x-direction with a distance of 10.16 mm from each other. In the y-direction we used two wires in pairs (current flowing in the same direction) as one effective wire with a current of about 3 A in each of the effective wires (see fig. 8a). The complete configuration in this direction is given by 4 of these wire pairs. The current in the wires of the upper layer was also 3 A, the bias field was 4 G. With this configuration, we could generate a stable situation with 4 MOTs operating above the surface (see fig. 8b). The gradient of the quadrupole in the x-direction is again $`\frac{dB}{dx}=20G/cm`$ and in the y-direction $`\frac{dB}{dy}=4.3G/cm`$. The calculated distance $`r_m`$ above the surface is about $`1`$ mm. The number of atoms in each individual MOT is about $`210^5`$. The distance between the different MOTs is 3.3 mm in the x- and 8 mm in the y-direction. To switch the MOTs individually, we applied a laser pulse split off from the MOT laser to turn the different MOTs on and off by blowing the atoms away (see fig. 8c,d). By turning off one MOT no significant change in atom number of the remaining MOTs could be detected.

## V Summary

In this paper, we have proposed different lattice configurations consisting of magnetic quadrupoles which are generated by 2 orthogonal layers of wires. These quadrupoles could be detected in the different unit cells by operating MOTs with the required polarisation configuration. We could furthermore demonstrate the possibility to use this two layer configuration to trap atoms in a 2x2 MOT array above a surface. We were able to turn the MOTs on and off individually by using an extra laser beam. Such a system of multiple MOTs could be useful to load an array of magnetic microtraps from an array of separate MOTs. The current number of $`10^5`$ atoms in the single MOT will be increased in future experiments by using the Rb dispenser in our chamber in a pulsed operation with higher currents.

## Acknowledgement

We acknowledge support of the European Research Training Network ”Cold Quantum Gases” under Contract No. HPRN-CT-2000-00125.
We thank Axel Görlitz for helpful discussions and proof reading.