Collective scattering and oscillation modes of optically bound point particles trapped in a single mode waveguide field

Collective coherent scattering of laser light induces strong light forces between polarizable point particles. These dipole forces are strongly enhanced in magnitude and distance within the field of an optical waveguide so that at low temperature the particles self-order in strongly bound regular patterns. The stationary configurations typically exhibit super-radiant scattering with strong particle and light confinement. Here we study collective excitations of such self-consistent crystalline particle-light structures as function of particle number and pump strength. Multiple scattering and absorption modify the collective particle-field eigenfrequencies and create eigenmodes of surprisingly complex nature. For larger arrays this often leads to dynamical instabilities and disintegration of the structures even if additional damping is present.


I. INTRODUCTION
Laser light scattering off free particles is accompanied by momentum transfer and thus forces. If the incident light fields are far detuned from any optical resonance, energy absorption and spontaneous re-emission play only a minor role and a coherent redistribution of photon momenta determines the force. As such coherent scattering from different particles largely preserves the laser phase, these fields interfere and thus the resulting force on each particle depends on its distance to nearby particles and acquires a collective nature, which depends on all other particle positions [1,2]. While in a 3D geometry the interaction strongly decays with distance and is averaged out by phase randomization from particle positions and motion, in regular particle arrays or in a spatially confined geometry, this interaction can become very strong [3,4]. Particularly strong effects appear if the light fields are optically confined in resonators [5], guided by optical structures as mirrors [6] or confined in waveguides [7][8][9][10]. Besides cold dilute gases, closely related effects have been studied using various kinds of nano-particles in solution [11][12][13].
Following first predictions by Chang and coworkers for near resonant weak pumping of atomic dipoles [7,10], we have recently shown that cold particles illuminated far of resonance, who can freely move along or within an optical nano-fiber field, tend to form regular but complex optical structures through optical long range coupling. They collectively scatter light into the fiber mode and self-trap in the optical potential generated by interference of the scattered light and the pump field [9]. Surprisingly we see that not only the light confines the particle motion but the particles also confine the light forming Bragg like atom mirrors at their outer edges. Numerical simulations predict that in certain cases even more complex structures as cavity arrays are self-formed by the particles [8]. This potentially generates a self-ordered cavity QED system with strong particle photon coupling via so called atomic mirrors [14]. While such dipole energy minimizing configurations were also studied for conventional 1D optical lattices, they could be shown to be generally unstable [15].
In this work we theoretically study the dynamics and stability of such self-organized cavity QED systems. Already in our previous work we could identify configurations where all forces on the particles vanish and one gets a restoring force towards these equilibrium positions for small shifts for each particle. Nevertheless, displacing one particle immediately changes the fields at all other particle positions and creates perturbation forces on them. Here we cannot expect any energy conservation for the particle motion as we are dealing with an open system where the pump laser forms a non depleting energy and momentum reservoir [16,17].
Hence perturbing one particle induces nonconservative collective motion of all particles.
Eventually these oscillations exhibit exponential growth and lead to a disintegration of the whole structure. A central aim of this work is to calculate and study the frequencies and stability of the corresponding linearized eigenmodes of the coupled atom field systems and their application to induce tailored long range interactions in this system. Of course the system is inherently nonlinear so that the full dynamics can only be captured numerically.
Our work is organized as follows. First we quickly overview the results of our recent paper [9]. Out of these, in order to investigate the stability of the particles, we linearize the forces acting on the particles around their equilibrium positions. Calculating the eigenvalues of the so obtained coupling matrix we find a condition which enables and limits the formation of stable configurations. With this we first study configurations in the negligible coupling limit and then pass over to more realistic complex coupling parameters.

II. SCATTERING MATRIX APPROACH TO LIGHT INDUCED MOTION IN 1D
Let us first briefly recall our previously developed multiple scattering approach to calculate the stationary fields and forces on a linear array of particles with transverse illumination as schematically depicted in Fig. 1 [9] . For a given spatial arrangement of point scatterers in the field of a 1D waveguide the light field can be calculated by a scattering matrix approach [17][18][19][20]. The contribution of each particle is represented by a single 3 × 3 matrix parametrized by an effective coupling constant ζ proportional to its linear polarizabilityα and the field mode amplitude at its position. This determines the interaction between the particle and the fiber mode as well as the scattering amplitude in and out of the fiber η. The transmission and reflection coefficients are related to the coupling constant via t = 1/(1−iζ) and r = iζ/(1 − iζ). For symmetric scattering the coupling between the amplitudes left and right of a particle can then be expressed using the following beam splitter matrix M BS [9] where B j , C j , η are the incoming fields allowing to determine the outgoing fields A j , D j , cf.
). Multiplying all the matrices including the propagation phase shifts then allows to determine the field distribution corresponding to the momentary particle positions. Note that we neglect the propagation time of the light within the structure compared to particle motion.
Using simple arguments based on the Maxwell stress tensor the momentary field then determines the optical force on the j-th particle along the x-axis [18,21]: We will use this optical force then as ingredient to obtain Newton's equations for the particle motion, where we add damping describing additional light or background gas induced friction. In previous work we and others have seen that under quite general conditions stationary equilibrium positions can be found, where all particles are force free and feel a restoring force against small shifts. However, in contrast to a standard 1D optical lattice, the absence of a force does not enforce an unaltered propagation of the light [18] but just guarantees a balance between internal and external light scattering [9]. Following Ref. [18] and including a linear friction term proportional to µ the equations of motion for the point particles of mass m then read: Stationary structures require F j (x 1 , . . . , x N ) = 0 for all j. Here these are not simple equidistant lattices, but rather complex particle arrangements with regular central parts including the formation of Bragg like reflectors at both ends of the equilibrium particle distribution

III. COLLECTIVE EXCITATIONS AROUND EQUILIBRIUM POSITIONS
When a particle is weakly displaced from its equilibrium j . In the following we study the spatial properties and eigenfrequencies of such small amplitude particle oscillations analogous to phonons of a lattice. Note that the existence of a stationary stable equilibrium configuration does not guarantee stability with respect to even very small such collective oscillations, which in general can grow in amplitude and disintegrate the particle array.
A Taylor expansion of the forces gives the following set of coupled differential equations with nonlocal coupling matrix Note that these coupling constants would be an ideal basis for a coupled oscillator model in the ultracold gas limit. For the classical point particle model here we simply use the ansatz ξ = Ae iωt in equations (4) and solve for the eigensystem of D, where the oscillation frequencies can be calculated from the eigenvalues λ ν via yielding the linearized solutions: The eigenfrequencies include a real part describing damping or antidamping depending on its sign. If positive, it leads to an amplification of the oscillation amplitude and the particles can not form a stable stationary configuration. Hence, to ensure stable configurations the eigenvalues have to fulfil the condition [18]: Note that added damping can restore the stability only as long as all real parts of the eigenvalues are negative. Without external friction, µ = 0, Eq. (7) simplifies to and as long as we have no imaginary parts of the eigenvalues λ ν and the real parts are negative the particles are simply oscillating around their equilibrium position.
A. Negligible mode coupling limit ζ = 0 In the special case when particle mode coupling is very weak but we have a strong transverse pump field applied, one can set ζ = 0 while still keeping a nonzero scattering into the fiber, thus η = 0. Here the scattering of the field by the particles within the fiber is neglected and any two particles interact equivalently. As predicted before [7] and also found by us in Ref. [9], in this limit the particles in a stationary state are equidistantly distributed and the transfer matrix can easily be calculated explicitly: Here the collectively scattered field intensities are symmetric to right and left and read and for the force on the j-th of N particles we get: with P η = I η /c the radiation pressure resulting from the transverse pump field.
This force vanishes for different lattice constants d = (2n − 1)λ/(2N ) with n ∈ N giving potential stationary solutions. We will now calculate the linearized coupling matrices D for a small displacement of the l-th particle within such a solution. For the amplitude perturbation at the j-th cloud from displacing the l-th cloud by ξ = k we have to distinguish between the two cases if it is at the right or left side of the displaced l-th cloud: with the perturbation b. For j < l the amplitudes can be calculated as follows Let us first look at the stability of the solution with the maximal possible lattice distance where the force on the j-th particle induced by the perturbation of the l-th particle reads: and so D jl = k F j is: D is a symmetric matrix with all eigenvalues λ ν real and as it is a circulant matrix, they and the eigenvectors z ν can easily be calculated to give: As we have seen in Eq. (9), stably oscillating configurations require real and negative eigenvalues. Hence this configuration is stable and we find N phonon modes with frequencies where we have defined ω 2,0 = 2Pηk m as binding oscillation frequency for two particles at ζ = 0. As expected from translation invariance of the system, the lowest ν = 1-mode corresponding to the center of mass motion has a zero eigenfrequency. Interestingly our result exactly coincides with the prediction based on nonlocal fiber enhanced resonant dipole-dipole interaction obtained before in Ref. [7]. In the limit of large particle numbers we find: Consequently the particle binding frequency grows with √ N for large particle numbers stiffening the structure with particle number. Fig. 2 shows the phonon frequency as a function of mode number for different particle numbers. It demonstrates that the second and the N -th mode posses higher frequencies than all the other modes which are almost degenerate.
It is now interesting to see how the atoms interact with respect to their distance. For this we insert Eq. (5) into Eq. (4) and calculate the perturbation induced forces explicitly: exhibiting a coupling between the oscillators proportional to sin |j−l|π N depending on the distance |j − l| and particle number N . Maximum coupling occurs between center and boundary particles at N = 2|j − l| and it vanishes for |j − l| N .
Let us now insert other zero force solutions at smaller distances d n = 2n−1 2N λ into Eq. (5). Here we can explicitly calculate the eigenvalues in the limit of large particle numbers As according to Eq. (8) fully stable configurations are only possible if all eigenvalues are negative, we see that the first case, n = N , was the only stable configuration with a distance below λ. An example how the particles evolve over time is plotted in Fig. 3, where we compare stable configurations setting n = N (Fig. (a)) with unstable ones setting n = N −1 ( Fig. (b)). Note that starting from a stable configuration we obtain correlated oscillations of particles and fields around the stable position with a slow nonlocal excitation transfer between the sites. At negligible coupling, ζ = 0, one finds neither damping nor amplification of the oscillation modes but we still see a higher light intensity confined at the center of the structure. Starting from an unstable zero-force configuration the particles keep their position only for a short time until two particles collapse and the other particles reorder.
Nevertheless, also this new configuration is not stable and the whole order disintegrates.
B. Collective dynamics for finite particle field interaction parameter ζ Particles within the fiber mode will not only couple pump light in and out of the fiber, but also scatter it to the opposite propagation direction, into free space or absorb it. This is effectively described by a non-zero complex value of ζ [18]. As a first consequence this immediately leads to spatially inhomogeneous fields in the mode and thus non-equidistant stationary configurations, which, in general, do not allow for analytic treatment. Although the equation for the force (2) still looks very simple, its explicit form gets already complicated for more than three particles and even the zero-force points can not be found analytically.
Hence to still get some analytical insights we first study systems of only two and three particles.

Two particles
For two particles with no fields injected along the fiber, the light scattered by one particle only interferes with the light scattered by the second particle. Translation and spatial inversion invariance of the system here implies a vanishing sum of the forces on the two particles as long as we neglect asymmetric chiral coupling [22]. The force as function of the real part ζ r and imaginary part ζ i then depends only on particle distance d [9]: and the zero force distances d 0 = λ 2 1 2 + n , n ∈ N are independent on ζ. Interestingly even for strongly absorbing particles as e.g. gold nano-particles ζ i >> ζ r , we thus obtain the same force free binding distances [4]. Intuitively this can be understood from the π/2 relative phase shift of the two scattered fields at this distance preventing interference. To study the stability and strength of binding we calculate the coupling matrix D: and its eigenvalues. Inserting the zero force distances d = 1 4 + n λ the nonzero eigenvalue reads and for d = 3 4 + n λ we get the same eigenvalue with opposite sign. As the non-zero eigenvalue for the latter case is real and negative, these are stable configurations. The phonon frequency characterizing the optical binding strength then reads: which for large values of ζ stays finite and converges to ω 2,0 √ 2 below the noninteracting value. Nevertheless for real ζ ω 2 reaches a maximum of ω 2,0 3 + √ 5 at ζ 0.618 ( Fig. 5(a)), while we find a minimum coupling frequency for blue detuning at ζ −1.618 of ω 2 = ω 2,0 3 − √ 5 ( Fig. 5(b)). The dynamics in these two extreme cases is shown in figure 5. Note that in both cases we can observe that the particles and fields oscillate and are periodically pushed apart whenever the intensity of the light trapped between them is maximal. Note that even in the blue detuned case for low field seeking particles we find stable trapping but much slower oscillations, i.e. weaker confinement.

Three particle dynamics
Let us now add a third particle. Due to symmetry we will restrict ourselves to equidistant configurations and calculate the total transfer matrix. Again the amplitudes of the electric field to the left and right of the particles yield the forces acting on the particles, where as a consequence of symmetry the force on the middle particle is zero and the remaining two sum up to zero, i.e.: Figure 4: Fig. (a) shows the phonon frequency ω ν as a function of the real part of the coupling constant ζ r . The blue line corresponds to ζ i = 0, the red line to ζ i = 1 9 and the green line to ζ i = 1 2 . Fig. (b) shows the phonon frequency ω ν as a function of the imaginary part of the coupling constant ζ i . The blue line corresponds to ζ r = 0, the red line to ζ r = 1 9 and the green line to ζ r = 1 2 . We can see that the function goes to ω 2,0 √ 2 for large values of ζ. For the real part of ζ we can find a maximum for ζ r > 0. F 1 = −F 3 = P η (cos(kd) + cos(2kd) − ζ(2 sin(2kd) + sin(3kd) + sin(4kd))) + O[ζ] 2 , Here the general expressions are complex so that we only present some special cases to exhibit the stability of the particles. For real ζ = 1 9 only one stable equidistant configuration exists with d 1 = d 2 (0.8276 + n)λ. Adding an imaginary part, ζ = 1+i 9 , the stationary distances are enlarged to d 1 = d 2 (0.8305 + n)λ. As shown in Fig. 6 the trajectories of the particles here show surprising features.
Although the stationary state is symmetric with respect to the middle particle, a small perturbation from this equilibrium induces an intricate superposition of two oscillations, namely fast relative oscillations superimposed on a much slower and large amplitude a center of mass oscillation. As the system possesses translational invariance implying momentum conservation, it could be expected in principle to not allow any center of mass oscillations. This is at first sight also confirmed by the oscillation eigenvalues {0, −3.83P η k, −3.51P η k} calculated for real ζ = 1/9 with eigenvectors The eigenvalue λ 1 corresponding to the center of mass oscillation z 1 is indeed zero. However, we can see that the eigenvector z 2 , corresponding to the case where two outer particles oscillate together against the middle particle includes a center of mass oscillation, as the middle particle moves much more then the sum of the outer two. This is also clearly visible in figure 6. At this point we have to recall that we have a strongly coupled atom field system where the pump laser constitutes an external energy and momentum reservoir, which can provide or accept momentum from the particles and conservation is only true for the combined atom field system. The asymmetric oscillation z 2 thus periodically channels momentum in and out of the light field inducing these unexpected center of mass oscillations.
The dynamics gets even more complicated if we allow for light absorption and investigate the eigensystem for an imaginary coupling constant, ζ = 1+i 9 . We find the eigenvalues {0, −2.98P η k, −2.89P η k} with eigenvectors Again we observe a large amplitude center of mass oscillation due to the non inversion symmetric initial perturbations which lead to even amplified center of mass motion. Hence we see that while we have stationary equilibrium configuration, the system tends to disintegrate due to dynamical instability when no external friction force is provided.

IV. NUMERICAL STUDIES OF THE DYNAMICS OF LARGER ENSEMBLES
While the momentary light field for a given configuration of many particles can still be found by simple matrix multiplication, the explicit expression for the forces soon becomes prohibitively complicated and their common zeros cannot be analytically obtained any more.
Numerically, however, this still can be done efficiently even for hundreds of particles and the solutions of the corresponding dynamical equations of motion are still possible. In order to capture typical features of the resulting many particle dynamics but still obtain readable plots, here we study the case of ten particles for various coupling strengths and damping.   (Fig. 7(a)) with the case of partial absorption, where ζ also possesses a small imaginary part ( Fig. 7(b)). In the first case a small perturbation simply induces small correlated oscillations in the vicinity of local field maxima with some slow energy exchange between the particles. In the second case with absorption, the particles start to oscillate first, but these oscillations grow until one reaches the nonlinear regime leading to a complete destruction of the order. As shown in a recent paper [9] one can still find a stable configuration for this case in the over damped limit. The mathematical origin for this instability is that due to the imaginary part of ζ the eigenvalues also acquire imaginary parts. This instability still can be suppressed up to a limiting value by introducing an additional friction coefficient as in Eq. (8).
We had demonstrated that in the weak-coupling-limit, ζ = 0, one has stable oscillatory solutions without the need of damping. Introducing ζ = 0 this behaviour is often lost and even for rather small perturbations the particle oscillations grow beyond the linearized regime. This leads to the collapse of the particle structure and they are expelled to the sides.
It is now interesting to study the grade of instability in such a system, which turn out to be very parameter and size sensitive. To quantitatively show this we compute the eigenvalues of D and examine the real parts of them characterizing the magnitude of the corresponding instabilities.
The maximum real part of all the nonzero eigenvalues of D is plotted in Fig. 8 as a function of real and imaginary part of ζ for different system sizes N . If this value is negative, in principle one can always add enough extra damping to stabilize the configuration, while a positive value implies dynamical instability in any case. The figures show that for odd particle number we get a boundary line to the left of which at least one eigenvalue is positive and Eq. (8) thus predicts an instability, which can not be removed by damping. For negative real ζ and even particle numbers (Fig. (b) and (d)) we find stable configurations with negative eigenvalues but they are very very close to zero. In this case the system is very sensitive to small perturbations and needs strong extra damping for stabilization. Here in numerical simulations we often find unstable dynamics despite using a nominally stable parameter set, as validity of the linearized equations is very limited.
An analogous phenomenon occurs for large particle numbers (Fig. (d)) and large values of ζ r . In this case it is difficult to distinguish between stable configurations with eigenvalues close to zero and unstable configurations. We can observe that for large ζ positive eigenvalues appear, especially for imaginary ζ, which prevent the formation of stable configurations. In general only for small and mostly real ζ (blue region in Fig. 8 ) robust stable configurations can be expected.

V. CONCLUSIONS
Coherently illuminated particles in a 1D trap order in strongly bound regular arrays induced by the interference of the scattered and pump light. If the scattered light is guided by nano-optical waveguide structures one finds very strong binding and long range order with phonon frequencies growing approximately with the square root of the particle number.
In general the scatterers arrange in configurations, where the scattered light is resonantly confined within the structure to maximize the optical trapping potential for the particles.
This can be seen as a controllable prototype 1D implementation of an atom light crystal.
In this regime absorption and scattering of the fields by the particles within the waveguide create a strongly nonlinear dynamical response of the system, which limits the size of stable arrays even in the presence of additional friction forces. This hints at the possibility to create self-ordered nano-cavity QED systems where the particles simultaneously form the resonator via atom mirrors [23] as well as the optical active system. As a next task towards a quantized field treatment one of course would need to estimate the magnitude of the selfordered single photon coupling strength, which requires a sort of dynamic mode description.
Generalizations to several frequencies, polarizations or higher order transverse modes should generate even more intriguing complex physics. In principle similar mechanisms should also be at work in planar 2D setups [13] leading to even stronger transverse light confinement and limiting the extension of optical matter.