Scattering approach to two-colour light forces and self-ordering of polarizable particles

Collective coherent scattering of laser light by an ensemble of polarizable point particles creates long range interactions, whose properties can be tailored by choice of injected laser powers, frequencies and polarizations. We use a transfer matrix approach to study the forces induced by non-interfering fields of orthogonal polarization or different frequencies in a 1D geometry and find long range self-ordering of particles without a prescribed order. Adjusting laser frequencies and powers allows to tune inter-particle distances and provides a wide range of possible dynamical couplings not accessible in usual standing light wave geometries with prescribed order. In this work we restrict the examples to two frequencies and polarisations but the framework also allows to treat multicolour light beams with random phases. These dynamical effects should be observable in existing experimental setups with effective 1D geometries such as atoms or nanoparticles coupled to the field of an optical nanofibre or transversely trapped in counterpropagating Gaussian beams.


Introduction
Coherent interference of light scattered from different particles in an extended ensemble of polarizable point like particles leads to important modifications of the forces on the particles as well as to new inter-particle light-forces, even if the light fields are far detuned from any optical resonance [1,2]. While a full 3-D treatment certainly leads to a very rich and complex dynamics [3], key physical effects can already by discussed in effective 1D geometries. One particularly interesting example are atoms in or close to 1D optical micro structures [4,5] as e. g. an optical nanofibre, where even a single atom can strongly modify light propagation and forces [6,7]. In a milestone experiment Rauschenbeutel and coworkers recently managed to trap cold atoms alongside a tapered optical fibre [5] and related setups predict and demonstrate strong back-action and inter-particle interaction [8,9,10] leading to the formation of periodical self-ordered arrays [11,12]. Alternatively, in free space interesting dynamical effects of collective light scattering were recently predicted and studied in standard 1D optical lattices of arXiv:1310.6246v3 [quant-ph] 18 Feb 2014 sufficient optical density [13,14]. One could also consider arrays of optical membranes to study such effects.
In this work we extend an existing model [13,14,15,16] towards light configurations with multiple frequencies and polarizations of the fields illuminating the particles. In particular, this includes a new class of geometries where crystalline order can be dynamically generated and sustained even without prescribing a standing wave lattice geometry. As a generic example the polarizations of two counter-propagating fields can be chosen orthogonally, such that incident and scattered fields do not directly interfere. Light scattering thus occurs for both fields independently and the forces on the particles can simply be added up. However, any structure forming by the scattering of one field component will be seen by all other fields and thus change their scattering properties and the induced forces. On the one hand this mediates nonlinear interaction between the different fields while on the other hand it generates inter particle interactions throughout the sample, inducing a wealth of nonlinear complex dynamical effects. Besides such dynamic self-ordering phenomena, we also study the possibilities to induce tailored long range interactions via multicolour illumination and collective scattering of particles trapped in prescribed optical lattice potentials.
This work is organized as follows: First we introduce the basic definitions and dynamical equations of the well-established generalized multiple scattering model for light forces [13,14,15,16] and extend this framework to support multiple polarizations and frequencies. This formalism is then applied to an orthogonal beam trap consisting of an array of particles modelled as beamsplitters irradiated by two counter propagating beams of orthogonal polarization and possibly different wavenumbers, cf. figure 1a. For two beam splitters we analytically derive conditions on the intensity ratios and wavenumbers to trap or stabilize them at a given separation. These results are then numerically extended to higher particle numbers.
As a generalization and connection to usual experimental setups for optical lattices we then analyze how an additional beam polarized orthogonally to a prescribed standing wave perturbs the trapped particles and induces peculiar interaction patterns, cf. figure 1b in 1D optical lattices.

Multiple scattering approach to multicolour light propagation in linear polarizable media
It is now well established that propagation of far detuned light through a one dimensional atomic lattice or an array of dielectric membranes can be well described in a plane wave approximation with multiple scattering by a corresponding series of beam splitters [13,14,15,16]. A very analogous situation arises when the light is transversely strongly confined by optical structures so that scattering dominantly occurs along a preferred direction.
The spatial dynamics of the electric field E(x, t) = E(x) exp(−iωt) is then described by the 1D Helmholtz-equation where N denotes the total number of beam splitters at positions x 1 , . . . , x N ; ζ := kηα/(2ε 0 ) is a dimensionless parameter proportional to the atomic polarizability α, the wavenumber k = ω/c and the density of particles combined to a single beam splitter, η.
The plane wave solution between two beam splitters, x ∈ (x j , x j+1 ), then reads the amplitudes A j , B j left and C j , D j right of the beam splitter at position x j (cf. figure 1) are related by the linear transformation matrix M BS , with From (2) we read off the propagation matrix The values of the electric fields then are fixed by the incoming beam amplitudes A 1 and D N . The total reflection and transmission coefficients are calculated from the total transfer matrix of the setup and give the remaining amplitudes at the boundaries as B 1 = r tot 1 A 1 + t tot D N and C N = r tot 2 A 1 + t tot D N self-consistently [16]. Using Maxwell's stress tensor [17] yields the time averaged force per unit area on the j th beam splitter as [14] This simple but powerful formalism to calculate the fields and forces on single atoms, atom clouds or other dielectric media such as membranes or elastic dielectrics allows to describe complex dynamics such as self-organization or even laser cooling in any effective 1D geometry [18,19,20,21].
Previous approaches were limited to a single frequency and polarization in a counter propagating geometry. Here we show that it is straightforward to generalize the beam splitter method to allow for multiple frequencies and polarizations. The field propagating in the x-direction shall then be written as where E y (x) [E z (x)] is defined as the component polarized in the direction of e y [e z ] oscillating with frequency ω y = k y c [ω z = k z c]. We want to emphasize that writing the the total field as a sum of linearly polarized fields in (6) is an arbitrary choice. None of the upcoming conclusions would change if we chose another orthogonal basis system (e. g. circular polarizations).
The main assumption of this work is that the particles do not scatter photons from one mode into the other. As long as this is fulfilled, the beam splitter model can be employed for each component independently. This assumption is obviously correct if the beam splitters are made of non-birefringent materials as they are used in many optomechanical experiments.
If the beam splitters are assumed to be single atoms one has to take additional care as these typically have tensor polarizabilities. In this case one would choose counterpropagating circular polarized waves in equation (6) because the two modes then address different atomic transitions. If we additionally assume sufficiently large detuning for each field, we can also neglect mixing due to spontaneous emissions into other Zeeman-levels. The polarizability then loses all spin and polarization dependencies resulting in a scalar quantity This is why the coupling parameter ζ introduced in equation (1) is proportional to a linear atomic polarizability and the wavenumber. For a generalization of the beam splitter method to multi-level atoms with tensor polarizability we refer to a work by Xuereb et al. [22]. In this work we will assume that the atomic polarizability α is the same for k y and k z , hence ζ z = ζ y k z /k y . Of course, a more realistic scenario is easily possible within our framework but it would add unnecessary complexity here. Our central goal is the study of multiple scattering dynamics and not the effect of optical pumping and polarization gradients.
In the following chapters we will study how the introduction of different frequency fields provides new prospects to manipulate arrays of particles, ranging from equidistant lattices to individually tuned inter-particle distances as well as the design and control of motional couplings.

Light forces in counter propagating beams with orthogonal polarization
In this section we explore forces and dynamics of a 1D lattice geometry modelled by a chain of beam splitters at distances d j := x j+1 − x j irradiated from both sides by light with orthogonal polarizations (e y and e z ) and possibly distinct frequencies (ω y and ω z ), cf. figure 1a. In contrast to a standard optical lattice setup as treated before [13,14] no a priori intensity modulation due to wave interference is present and we start with a fully translation invariant field configuration. Hence the light field itself does not prescribe any local ordering and only multiple light scattering from the particles themselves creates local trapping forces. Due to the translation invariance of the setup no stable particle configuration can be expected. However, the coupled particle field dynamics still induce relative order. Hence our central goal is to find conditions, when the light forces induced by two non-interfering beams are nevertheless sufficient to obtain stationary stable particle arrays and how this spontaneous crystal formation arises.

Stability conditions for two beam splitters
To get some first insight, we start with the simplest nontrivial example of two beam splitters at a distance d = x 2 − x 1 . The intensities from the left and right beam are given as I y = cε 0 2 |A y,1 | 2 and I z = cε 0 2 |D z,2 | 2 , respectively. Here we chose the convention that all variables with index y correspond to light polarized in direction of e y which is injected from the left (negative x− axis) and index z corresponds to e z polarized light injected from the right. The individual beam splitters are counted from left to right with integer indices, hence, for example B z,2 is the B-amplitude of the light field polarized parallel to e z at the second beam splitter, cf. (2).
Using (3) and (4) to compute the fields for any given distance d, it is straightforward to obtain the total force on each beam splitter by simply adding the forces generated by the light in each polarization, i. e. F 1 = F y,1 + F z,1 and F 2 = F y,2 + F z,2 . The individual forces F y,1 , F z,1 , F y,2 and F z,2 are obtained from (5).
Despite the simple physical situation the corresponding general analytic solution already is rather unhandy. Thus we first restrict ourselves to real valued ζ neglecting absorption in the beam splitter or equivalently neglecting spontaneous emission in an atom fibre system. Assuming small values of ζ and dropping terms of O(ζ 3 ) and higher, we then find the following approximate formulas for the force on the two particles: Distance d Light force onto the left (7) (solid lines) and right beam splitter (8) (dashed lines) as function of their distance d for ζ = 0.01 and k = k y = k z = 2π/λ. For equal intensity P = I z /I y = 1 and frequency (blue curves) the two forces add to zero and vanish at distances d = λ/8 and d = 3λ/8. For asymmetric intensities P = 0.7 (red curves) we find distances with equal forces F 1 = F 2 but a net center of mass force remains. The black curve shows a similar behaviour occurring for different wavenumbers kz ky = 1.2 of same power P = 1. The red (green) dot marks unstable (stable) stationary points.
For a given set of control parameters, i. e. the intesity ratio P := I z /I y and the wavenumbers k y and k z , the beam splitters will settle at a distance d 0 for which the two forces are equal, i. e. , F 1 | d=d 0 = F 2 | d=d 0 , and the configuration is stable ( In this case, the system can still exhibit centre of mass motion but the particles keep a constant distance. From (7) and (8) we see that a stable configuration in the special case of equal wavenumbers, k = k y = k z requires Independent of the injected laser intensities, which just appear as a multiplicative factor, this corresponds to a pair distance d s 0 = (2n+1)π 4k (n ∈ N). Here the solutions for odd n correspond to a stable configuration, while even n leads to unstable behaviour. As numerical example we plot the full distance dependent forces for three typical sets of parameters in figure 2, where stationary distances of equal force can be read of the intersection points. If these occur at zero force, the centre of mass is stationary as well. For small ζ these distances of zero force on each particle can be approximated by a) b) Figure 3. a) Force on left (orange) and right beam splitter (green) as function of the wavenumber ratio k z /k y and distance for two partly absorbing beam splitters with ζ = 1/12 − i/150 and P = 1 (left figure). b) Stationary distance of two beam splitters with ζ = 0.01 as function of the wavenumber ratio k z /k y and the intesity ratio P obtained by numerical integration of their equation of motion including an effective friction term.
Conditions (10) and (11) imply gives a stable and stationary configuration, where the forces on both beam splitters vanish and small perturbations induce a restoring force, as cf. figure 2. In general, such solutions can only be determined numerically and are not guaranteed to exist for every set of parameters. In Appendix A we show that the line of argument can also be reversed and one can calculate the intensity ratios and wavenumbers needed to obtain a stable configuration for a desired distance d. This allows precise distance control of the particles via intensities and frequencies.
Let us now exhibit some more of the intrinsic complexity of the system in a numerical example. In figure 3a we first plot the forces on the two beam splitters as function of distance and relative wavenumber for fixed equal intensity from both sides. Clearly the intersection of the two force surfaces exhibits a complex pattern with a multitude of stationary distances which can be controlled e.g. via the chosen frequency ratio.
In an alternative approach we can numerically find a stable stationary distance of the two beam splitters as function of intensity and wavenumber ratio by time integration of their motion with some damping added, cf. figure 3b. We see that depending on the parameters for a given initial condition the system can settle to a large range of different stationary distances, exhibiting rather abrupt jumps at certain critical parameter values. Generally a numerical evaluation requires very little effort and can be easily performed for large parameter ranges. Despite the fact that there is no externally prescribed order, the particles mostly tend to arrange at configurations with stationary distance.

Self-ordering dynamics for higher numbers of beam splitters
In principle, determining stationary states for a larger number of beam splitters is straightforward by first solving (3) and (4) for the fields and using these to calculate the forces. However, to determine a completely stationary configuration of N beam splitters for a given input field configuration, we have to solve N nonlinear equations to guarantee a vanishing force at each particle as function of the N − 1 relative distances. This problem can have no or infinitely many solutions. Often one does not get an exact solution, but solutions with vanishingly small centre of mass force.
As a rather tractable example we plot the zero force lines as function of the two relative distances for the case of three beam splitters illuminated by light of equal power, P = 1, but different colour, k z /k y = 1.1, in figure 4. One finds many intersections of these lines, where two forces vanish, but only for a few distances we get triple intersections where the forces on all three particles vanish and stationary order can be achieved. These solutions then still have to be checked for stability against small perturbations to find a stable steady state.
To investigate the dynamics of a higher number of beam splitters it is more instructive to solve the dynamical equations of motion for various initial conditions until an equilibrium configuration is reached. To arrive at a stationary solution we assign a mass to the beam splitters and add an effective friction coefficient µ in the classical Newtonian equations of motion, In the following simulations we assume that the system is in the so called over damped regime, meaning that the characteristic time scale of undamped cloud motion, i. e. the oscillation period, is much smaller than the relaxation time of the cloud's velocity towards a constant value due to viscous friction. Under this assumption the equations of motion (12) reduce to a set of differential equations of first order [14], In figure 5 we show the solutions of (13) for ten beam splitters in a simple orthogonal beam trap with I z = I y and k z = k y = k. In a traditional standing wave trap, the beam splitters would settle at the chosen initial equidistant spacing d OL ≈ λ/2, cf. (23), which can be determined self consistently [14]. However, for two trap beams of orthogonal polarization no prescribed periodicity is present and the particles themselves create field configurations which confine their motion through multiple scattering. Our simulations show that for a large range of operating conditions the light forces generated by two counter propagating beams with orthogonal polarizations will indeed induce an ordering of the particles, i. e. multiple scattering between the beam splitters is sufficiently strong to generate a stable configuration.
Interestingly, the final distances d 1 = d 2 = . . . = d N converge to the same result as obtained for the standing wave optical lattice as the number of beam splitters N is increased, cf. figure 6. Thus orthogonally polarized trap beams have the same trapping properties as a standing wave setup, as N → ∞. The beam splitters themselves then form an effective Bragg reflector to synthesize a standing wave configuration, which traps the particles.
A substantially more complex behaviour is found for the case of two colour illumination with different intensities, I y = I z . The trajectories for some representative cases can be found in figure 7. Interestingly, there is still a wide range of parameters where one obtains stationary patterns, but we generally get a non-equidistant spacing, cf. figure 8, and a finite centre of mass force. As above for two beam splitters, this . For large N we observe asymptotic convergence towards the effective lattice constant as found for a standing wave configuration, cf. equation (23) or [14]. For the red dots we use ζ = 0.01. A small imaginary part of zeta (green dots) ζ = 0.01 + i0.001 decreases the distances but still yields stable configurations. . The colour coding in the background shows the total field intensity I tot := I y + I z during the reorganization process of the system. On the right hand side the trajectories (N = 10) for P = 1.3, k y = k z are shown. In both cases we observe a finite centre of mass force in the long time limit, while the pattern formed is stable but is no longer equidistant.
force can be controlled via the intensity ratio to stabilize the centre of mass or induce controlled motion. Of course, the configuration not only depends on the operating conditions, but also on the initial conditions allowing for a multitude of different stationary configurations.
In summary we conclude that the particles prefer to form crystalline structures held in place by collective multiple scattering. The more particles we have, the more complex these patterns can get and the more different solutions can exist. The complexity of the problem increases further, if one allows for a variation of the individual coupling parameters ζ, e. g. to represent number fluctuations of the atoms trapped at each lattice site or size variation of trapped beads. Note that although appearing similar at first sight, the mechanism is different from standard optical binding of polarizable beads, which works on transverse shaping of the incoming light with the particles acting as small lenses [23], which we neglected in our model. Let us finally note that analogous results should be obtained for a setup using two counter propagating beams of equal polarization, but sufficiently different frequencies, such that scattering between the different colours is suppressed. From the particles point of view, the interference pattern of the combined fields then oscillates so rapidly that they cannot follow and the two forces stemming from the two fields can be calculated independently. Such frequency shifts are a common method to generate 3D optical lattices by using a different frequency in each dimension. But in contrast to those cases, here we get a mutual interaction between the light intensities of the different frequency components. During the evolution the spatial shifts of the beam splitters induced by one field are seen by all other fields and influences their propagation.

Tailored long range interactions in a bichromatic optical lattice
Optical lattices for ultracold atoms are of course an extremely well established and controllable technology. In general, parameters are chosen in a way to avoid back-action of the particles on the fields. The underlying physics helps here to achieve this goal as particles tend to accumulate in zero force regions, where their influence on the lattice light is minimized [13,14]. This is radically changed in the orthogonally polarized beam setup described above, where trapping forces are only created by the back-action of the particles on the two beams and interactions in the lattice occur via multiple collective scattering.
In the following chapter we will consider a second generic example to generate tailored long range interactions in an optical lattice. In particular we study the extra forces introduced by a second perturbation field of different wavelength in a given optical lattice formed by to two strong counter propagating beams of equal wavenumber k and polarization e y , cf. figure 1b. By adding an extra beam of different wavenumber k p and polarization e z we can introduce tailored perturbations and couplings, as its gradient is generally non-zero at the positions of the original lattices sites.
For generality we allow intesity asymmetries for the dominant standing wave field P := I r /I l where the first indices l and r stand for left and right suggesting the direction of incidence. The intensity of the additional perturbating field is called I p . The same index notation will also be used for the corresponding field amplitudes.

Two beam splitters in an bichromatic optical lattice
The first relevant system to study interactions and couplings introduced by an additional field of different frequency are two beam splitters trapped at a distance d sw , cf. (23), in a far detuned optical lattice at stable positions x 0 1 = x 0 − d sw /2 and x 0 2 = x 0 + d sw /2. Here x 0 denotes the centre of mass coordinate calculated following [16], via with u = sgn[Im(rt * )]. The reflection and transmission coefficients r and t of the total system derived from the total transfer matrix are The incident fields of the standing wave component are assumed as A l = 2I l /(ε 0 c) exp(ikx 1 ) and D r = 2I r /(ε 0 c) exp(−ikx 2 ) so that the remaining amplitudes at the boundaries are B l = rA l + tD r , C r = tA l + rD r . This allows to calculate the lattice forces on the first and second particle, The additional perturbation field is described by the amplitudes A p = 2I p /(ε 0 c) exp(ik p x 1 ), B p = rA p and C p = tA p and generates the additional forces (ζ p = k/k p ζ)  Figure 9. Dependence of the coupling constants κ 1 (red) and κ 2 (blue) on the perturbation field intensity I p for ζ = ζ p = 0.1, k = k p and P = 1. As soon as the perturbation field is switched on, the constants differ.
Here we restrict the corresponding added dynamics of the two beam splitters to small, time dependent perturbations ∆x 1 (t), ∆x 2 (t) << d sw from the equilibrium positions x 0 given in (14). Using x 1 (t) = x 0 − d sw /2 + ∆x 1 (t) and x 2 (t) = x 0 + d sw /2 + ∆x 2 (t) and linearising the forces for small ∆x 1 (t) and ∆x 2 (t) yields to the following coupled equations of motion A detailed calculation of the coefficients K, κ 1 , κ 2 and F ext is shown in Appendix B. The equations above correspond to two coupled harmonic oscillators driven by an external force F ext .
Note that the coupling constants κ 1 and κ 2 here are not necessarily equal, cf. also figure 9, as there is no energy conservation enforced for the motion of the two beam splitters. Since the parameters can be chosen in a way so that the coupling constant κ 1 is equal to zero, one-sided couplings can be achieved. This means that only the motion of beam splitter number two is coupled to beam splitter number one, which does not couple to the rest of the system. The direction of this effect is governed by the direction of incidence of the perturbation beam. Besides, κ 2 > 0 holds for all values of I p , meaning that no antisymmetric modes can be obtained if κ 1 < 1, cf. (22). Generally we find that tuning the perturbation field intensity offers a variety of different dynamics not accesible with traditional standing wave setups. This motivates a more detailed treatment of this system, using numerical methods.

Long range coupling of beam splitters in an optical lattice
As shown in [14] the effective self-consistent lattice constant in a standing wave with asymmetry A := I l −Ir √ I l Ir , with I l := 0 c |A l | 2 , I r := 0 c |D r | 2 adjusts to In an optical lattice, multiple scattering induces long range interactions between the particles in the form of collective oscillation modes. In the selfconsistent configuration the particles arrange at intensity maxima at minimal field gradients, so that this interaction is strongly suppressed for small perturbations. Adding, however, a second, perturbative field by a single running wave beam of wavenumber k p injected from one side induces an additional force on each particle perturbing the regular periodic order. This perturbation then acts back on the original standing wave field. Note that a single plane wave by it self would only add a constant force, but this force is modified by multiple scattering depending on the particle distances. As an instructive example we show these perturbing force acting on N = 4 beam splitters in figure 10. We see that the additional force is different for all the particles and changes as a function of the lattice constant relative to the wavelength of the perturbing light. Hence, by a proper choice of parameters almost any combination of magnitudes and signs of forces on the different beam splitters can be achieved.
This can be exploited for different purposes to control and study lattice dynamics. As a first and direct application it is possible to tailor a specific field to induce oscillations of selected particles in the optical lattice by deflecting them from their equilibrium position. As shown in figure 10, the force on the individual beam splitters depends strongly on the prescribed lattice constant. This means that there is a wide range of realizable dynamics as long as the lattice constant can be tuned, cf. (23). This can be potentially refined by simultaneous use of several perturbation frequencies.
For example, using the parameters from figure 10 we anticipate interesting behaviour for a lattice with spacing d sw ≈ 0.23λ p + mπ, m ∈ N as in that case F 1p = F 3p = −F 2p = −F 4p (cf. dash-dotted line in figure 10). This behaviour is verified by calculating the trajectories via (12), the results are shown in figure 11. Obviously it is possible to correlate the motion of distant beam splitters in an optical lattice via the additional beam. After switching on the perturbation at t = 0, particles number one and three show amplified oscillations, while the other's oscillations are damped.
In a second approach the additional field is designed to enhance interactions between selected distant areas in the lattice. As shown in figure 12, exciting an oscillation of one particle weakly coupled to the standing wave field will usually have little effect on the other trapped particles. But after adding a perturbative field with carefully chosen parameters, this oscillation can be transferred to the other particles, forcing them to move along. Note that due to the fact that the additional perturbation (coupling) field is imposed from only one side, this coupling effect is not symmetric and excitations can flow in a desired direction. For example, a perturbing field entering from the left hand side will maximally transfer the motion of the rightmost beam splitter. This setup allows one to correlate the motion of particles and can even be used as a channel to transfer e. g. quantum information through the lattice, very much like phonons in ion crystals.  Figure 12. Example plot for the resonant coupling of three beam splitters, trapped in a standing wave configuration with d 0 = λ p /2. We used ζ = 0.01, ζ p = 0.1, I l = I r = 20I p and k/k p = 0.99. The rightmost beam splitter is displaced from the equilibrium position at t = 0, resulting in a damped oszillation (damping parameter µ = 0.01). The black curves show the resulting dynamics for I p = 0. The green (x 3 (t)), blue (x 2 (t)) and red (x 1 (t)) curves show the dynamics if the coupling field is switched on. Note the resonant coupling between x 1 and x 3 .

Conclusions
We have shown that even in the case of non-interfering counter propagating light fields of different polarization and frequency, stable lattice configurations of particles held in space by multiple coherent scattering are possible. In contrast to conventional optical lattices the light here plays a decisive dynamical role as multiple scattering is essential to form and stabilize the structure. Compared to prescribed optical lattices the physics is much closer to the case of solids, where lattice dynamics in form of phonons not only keeps the atoms in place but also mediates long range interactions. Interestingly in conventional optical lattices such interactions can be tailored by adding additional coupling fields of suitable frequency and polarization. While we have performed our calculations only for 1D geometries, where a semi-analytic scattering approach can be used, similar effects should be present in 3D geometries as well.
In general for very far detuned optical fields these effects will be rather small but their importance will grow with the size of the lattice as well as in transversally confined fields. Particularly strong effects can be expected in fields guided by nano optical devices such as nano fibres or hollow core fibres. Here even for a few particles strong interactions can be expected.
In this work we have restricted ourselves to the bichromatic case for sake of simplicity. Nevertheless one can expect even more complex dynamics for an increasing number of input fields as the forces show a more complex distance dependence. Note that here we have ignored any internal optical resonances of the particles. Working close to such resonances certainly should strongly increase the effects but also will complicate the analysis.
Let us finally mention here that the system not necessarily requires a fixed set of beam splitters as a starting point. As an alternative we can consider each beam splitter to be formed by a small sub ensemble of atoms in a 1D beam configuration, as it has been proposed before [13,14]. In our case of noninterfering counterpropagaing beams, one can expect that under suitable conditions the cold atoms arrange in small groups forming at local field intensity maxima [12]. Groups of atoms at certain spatial sites then commonly form beam splitters shaping a self consistent lattice structure.
In contrast to conventional lattices, backaction of the particles onto the fields is an essential part of the dynamics and the field thus strongly mediates collective interactions. Light scattering on one end of the lattice influences the lattice depth at the other end, which opens a completely new branch of ultracold atom optical lattice physics. Note that also atoms trapped in optical resonator fields [24] exhibit similar dynamical coupling effects but in that case the backaction is strongly restricted by the resonator geometry limiting the available interaction wavenumbers. 2) (dashed line) for a prescribed distance d for ζ = 0.01 and k z /k y = 1.4. The grey regions mark the physically allowed regions P 1 > 0 and P 2 > 0. It impossible to find intensity configurations so that particles can be trapped at distances outside of these regions.
total force F 1 + F 2 vanishes, we solve P 1 = P 2 , finding k ± z = 1 d arccos ± cos(2dk y ) 2(1 + 2 cos(2dk y ) + 2πn, n ∈ Z. Obviously there exists a wide range of parameters which allow stable and trapped configurations of beam splitters in multicolour light beams with orthogonal polarizations (or sufficiently different wavenumbers).

Appendix B. Linearisation of the forces on two beam splitters in a bichromatic optical lattice
In section 3.1 we calculate the equations of motion for two beam splitters in a standing wave geometry perturbed by an additional field with orthogonal polarization. For that purpose we use linearized forces (17), (18), (19) and (20). Here we want to show how this linearization is done and how the constants K, κ 1 , κ 2 and F ext can be calculated.
The force F 1 depends only on the positions x 1 (t) and x 2 (t) of the two beam splitters. Replacing these variables via x 1 (t) = x 0 −d sw /2+∆x 1 (t) and x 2 (t) = x 0 +d sw /2+∆x 2 (t) results in a force dependent on ∆x 1 (t) and ∆(t) := ∆x 2 (t) − ∆x 1 (t). Assuming small ∆x 1 (t) and ∆(t) we perform a 2D Taylor approximation to first order resulting in F 1 = a + b∆x 1 + c(∆x 2 (t) − ∆x 1 (t)) (B.1) where we defined real constants a,b and c which are lengthy expressions depending on the system's parameters.
The same method works for the remaining forces F 2 , F 1p and F 2p , where the latter two only depend on the relative distance ∆(t).
Performing these tedious calculations we find that some of the obtained constants are zero (a = u = 0) while others have the same values. We define −K := b = v, κ 1 := K 2p + c, −κ 2 := K 4p + w with c = w and F ext := K 1p = K 3p . With this results result we get the forces used in the equations of motion (21).