Projecting light beams with 3D waveguide arrays
Free-space light beams with complex intensity patterns, or non-trivial phase structure, are demanded in diverse fields, ranging from classical and quantum optical communications, to manipulation and imaging of microparticles and cells. Static or dynamic spatial light modulators, acting on phase or intensity of an incoming light wave, are the conventional choices to produce beams with such non-trivial characteristics. However, interfacing these devices with optical fibers or integrated optical circuits often requires difficult alignment or cumbersome optical setups. Here we explore theoretically and with numerical simulations the potentialities of directly using the output of engineered three-dimensional waveguide arrays, illuminated with linearly polarized light, to project light beams with peculiar structures. We investigate through a collection of illustrative configurations the far field distribution, showing the possibility to achieve orbital angular momentum, or to produce elaborate intensity or phase patterns with several singularity points. We also simulate the propagation of the projected beam, showing the possibility to concentrate light. We note that these devices should be at reach of current technology, thus perspectives are open for the generation of complex free-space optical beams from integrated waveguide circuits.
The ability to shape the wavefront of a light beam, producing non-trivial intensity or phase patterns, is of interest for several applications. In particular, the possibility to exploit the spatial degrees of freedom of a multi-mode field to encode information has attracted increasing research in the recent years, with special interest in orbital angular momentum, both of classical and quantum light states, for free-space [1, 2, 3] or fiber-based communications [4, 5]. Mode beams with angular momentum  as well as structured light [7, 8] are also exploited to trap and manipulate micro and nanoparticles, or even living cells, through optical forces [9, 10]. In fact, the possibility to draw complex interference patterns in the propagation provides enhanced control on the particle position, rotation or movement, capabilities of former importance when dealing with the analysis and sorting of single cells.
Bulk optics components, such as phase or intensity masks [1, 2, 6], as well as dynamic spatial light modulators [3, 4, 8], are the almost exclusive choice to shape the wavefront of such free-space propagating beams. Therefore, the interfacing with fiber networks, integrated optical circuits, or fluidic lab-on-a-chips can pose stability or alignment problems, and limit the effective exploitation of the advantages of these beam shaping methods.
On the other hand, three-dimensional waveguide lattices [11, 12, 13, 14, 15, 16] have emerged in the last decade as powerful devices to investigate engineered multi-mode fields. Due to the high level of control that can be achieved in such structures, they have been employed to observe fundamental effects such as spatial solitons  or edge states in topological insulators , as well as to simulate a variety of quantum physics phenomena[14, 15, 16].
While waveguide lattices with a simple structure may be induced by external optical fields in photorefractive materials , it was the three-dimensional fabrication capability of femtosecond laser waveguide writing  that fostered much of the work in this field. In fact, this micromachining technique is based on the non-linear absorption of femtosecond pulses, focused in the bulk of a dielectric material: judiciously chosen irradiation parameters result in a permanent refractive index increase, localized in the focus region. Translation of the substrate with respect to the laser beam allows to inscribe with full three-dimensional freedom the desired waveguiding paths. However, up to now, the research was mainly focused on the propagation properties of the light within the waveguide arrays and the far-field emission properties were not specifically considered.
Here, we explore the possibility to use properly engineered waveguide lattices to project free-space light beams with controlled wavefront. We numerically and theoretically investigate, also by using relevant examples, the properties of the emitted light from a pattern of waveguide modes, distributed on a two-dimensional cross-section, as it can be achieved with three-dimensional waveguide arrays. In particular we study the possibility to produce peculiar far-field distributions, and one or multiple phase singularity points. We will also show how to obtain a certain degree of control on the propagation properties of the emitted beam, including focusing.
2 Models and approximations
This work investigates the properties of the output emission from waveguide arrays, in which identical single-mode optical waveguides are placed in arbitrary position on a two-dimensional cross-section. We will assume that it is possible to produce arrays maintaining full control on the intensity of the different modes, as well on their relative phase. While we will not go into details of the design and engineering of the arrays, we will discuss their effective feasibility in general terms in section 6.
We will consider the arrays as illuminated with linearly polarized light, with a fixed polarization direction. Of course, including also the polarization degree of freedom would open even wider possibilities on the achievable output light states, but this would be out of the scope of this preliminary study. We will further restrict our study to the paraxial condition.
In this framework, it is common to treat the electromagnetic field as scalar and the calculation of the field distribution during the propagation be achieved in terms of scalar diffraction theory [19, 20]. Considering our case in more detail, the output of an array of identical single-mode optical waveguides may be conveniently modelled with an array of Gaussian modes. If we drop uniform phase terms, the field of a Gaussian mode at a certain propagation coordinate , and as a function of the transverse position , reads:
where the beam radius and the wavefront curvature radius depends on as follows:
Now, the Gaussian modes in the array are assumed to be all identical with regard to (which is placed at the output facet of the array, where we can put the origin =0), but they may have different intensity and phase terms. If are the positions of the waveguides in the transverse plane, their intensity and their phase terms, the scalar field configuration of the output of an array of waveguides can be written as:
The above expression can be used to calculate the field distribution at an arbitrary coordinate within the paraxial approximation, and all numerical simulations shown in this work are based on that.
3 Polygonal arrays and angular momentum
Optical beams possessing angular momentum have attracted increasing interest in the recent years, for diverse applications: on the one hand orbital angular momentum is a powerful degree of freedom for encoding information, both classical and quantum [1, 2]; on the other hand the possibility to exert torque on trapped microparticles opens novel possibilities in optical micromanipulation . In this scenario it can be important to produce such vortex beams directly from optical fibers and waveguides, with reduced or no use of bulk optical components, which require critical alignment. Different approaches have been proposed, from ring resonators  to circular gratings [23, 24], to hybrid plasmonic waveguides .
Three-dimensional waveguide arrays have been also recently investigated to this purpose. Markin and coauthors proposed to directly excite supermodes of a triangular array of coupled waveguides, having cyclical phase distribution, by means of a phase-matched nonlinear process . Guan et al. showed the possibility to encode 15 different states with cyclical phase distribution in a set of 16 non-coupled waveguides, placed on the vertices of a regular polygon: this was achieved by splitting the light in the different waveguides with an integrated-optics star coupler and by using thermo-optic phase shifters for fine adjustment of the phase terms . However, these previous works were limited to the study of the propagating field within the array or to the near field distribution, while the far-field distribution was not studied or employed.
As a matter of fact, one could show that a single linearly-polarized mode or supermode,or a set of separated linearly-polarized waveguide modes has in each point a transverse component of the Poynting vector with vanishing time average (see A). This prevents light propagating on that mode or those modes to have a physical, non-vanishing, longitudinal angular momentum within the array, even in limited regions. However, this does not hold for overlapped optical modes in free space propagation, and thus stimulates the investigation of the free-space beams emitted at the output of this kind of waveguide arrays.
Namely, we study the far-field distribution produced by the output of an array of waveguides, placed (in the cross-section) at the vertices of a regular polygon, with a cyclical phase distribution, similar to the one demonstrated in Ref. . As discussed in section 2, the near field of the output of such an array can be described with good approximation as a set of Gaussian optical modes, that we place in the positions:
where is the radius of the polygon, and we consider the phases of the different modes to be arranged as:
Figure 1 shows some examples of near-field and far-field distributions for some of the simplest cases in which = 1 (phases make a 2 round-trip around the polygon), and polygons with = 4, 6, 8. We consider mode size ( = 15 m) and wavelength ( = 633 nm) commonly employed in three-dimensional waveguide arrays fabricated by femtosecond laser writing . The far-field distribution is calculated by directly evaluating (3) for = 15 mm. The presence of a phase singularity in the center of the far-field distribution is easily observed in all cases, associated with an intensity distribution that closely resembles the typical doughnut shape of the Laguerre-Gauss modes, i.e. the archetypical examples of beams with non-vanishing angular momentum .
for small and for large (, and are the radial, azimuthal and longitudinal coordinates respectively). Indeed, this same expression can be found as a limit for small in the case of Laguerre-Gauss beams  or Bessel beams  of order , which are well-known examples of beams carrying orbital angular momentum. Thus, this ensures that the beam is carrying orbital angular momentum proportional to in its central region.
In fact, field distributions with different angular momentum values can be produced as shown in figure 2, which reports the far-field distribution of a polygonal waveguide array with = 10 waveguides, and phases arranged cyclically for different values of in (5). Different phase vortices can be realized with both clockwise and counter-clockwise orientation, depending on the sign of . On the other hand, for = 0 the far-field distribution converges to a pattern with a single peak in the middle. The doughnut shape of the intensity distribution is less defined for increasing ; in fact, the doughnut would extend farther from the center, where the approximations of B become less accurate. We can further observe that the periodicity of the phase cycle needs not to be commensurate with the number of waveguides composing the array; for instance it is possible in our case to achieve = 3. Finally, we should of course note that it is not possible to effectively produce vortices with equal or higher than : from the sampling theorem it is clear that in distinct phase values (as available in waveguides) there would not be enough information to produce such vortices.
4 Towards more complex field configurations
We shall now discuss the potentials of our approach in producing more complex far-field distributions.
The first case we consider is illustrated in figure 3. Eight waveguides are placed on the vertices of an octagon, and are divided into two series (say, the odd waveguides and the even waveguides). All waveguide modes carry the same intensity. The odd modes have cyclic phases arranged in counter-clockwise direction, while the even modes have cyclic phases arranged in clockwise direction. A phase shift is further added to the phases of the even modes. This produces the two-lobes far-field distributions, whose orientation depends on , that can be seen in figure 3 (b) and (c).
If only the odd or even waveguides were present, we would observe in the far-field a vortex beam distribution such as those in figure 1, with an asymptotic similarity (in the central part) to the Laguerre-Gauss modes with = 1, as shown in B. However, here we have the linear combination of two vortices with opposite direction. It is known that the equally weighted combination of two Laguerre-Gauss modes with = 1 produces a two-lobe Hermite-Gauss mode, oriented along an axis which depends on the phase delay between the two Laguerre-Gauss modes. Indeed, we observe an analogous behaviour in figure 3.
We find this example interesting because it shows the possibility to change profoundly both the phase and the intensity profile of the far-field just by acting on the phase terms of the modes in the array, starting from the same balanced intensity distribution. In fact, from the same eight waveguides arranged on the vertices of an octagon, one could produce phase vortices with different orientation or vorticity such as the ones in figure 2, or other configurations with zero angular momentum such as the ones in figure 3, just by acting on those phases. In waveguide devices, one may think of accurately controlling such phase terms even in a dynamic fashion, by e.g. thermo-optic phase shifters [27, 30].
To guide the design process of waveguide arrays that produce phase and intensity patterns in the far-field with even richer features, it can be useful to go back to a more general picture of the involved diffraction phenomena. In fact, in the Fraunhofer diffraction condition, i.e. for sufficiently large , it is well known that the far-field can be approximated by the Fourier transform of the near-field . More precisely, by dropping common phase terms, one can write in our case:
where the transform is defined as:
and indeed consists in a properly scaled two-dimensional Fourier transform. Practically, the far-field distribution in our case is given by an interference pattern, which is the Fourier transform of the position of the modes, modulated by the diffraction pattern of a single waveguide mode .
One can think, as a consequence, to devise the structure of a waveguide array by, first, anti-transforming a certain desired far-field distribution, and then trying to approximate with a pattern of Gaussian modes (with identical shape but different phases or intensities) the result of this anti-transformation. Naturally, knowledge of the fundamental properties of the Fourier transform is of great help. In particular, it may be useful to remind that the transform of a regular lattice of peaks is again a regular lattice (the reciprocal lattice). Far-field configurations showing regular lattice structure may thus be easier to reproduce by means of the far-field emission of a waveguide array (with a certain amount of regularity).
Figure 4 shows two peculiar examples, that have been engineered following the above described procedure, i.e. by approximating with the waveguide position, intensity and phases, the anti-transform of a certain desired far-field. The first one (pictures on the top row) was designed looking for an “eyeglass” shape, with two counter-rotating and neighbouring phase vortices. The starting point to devise the second one (pictures on the bottom row) was a 22 matrix of phase singularities, two encircled by a clockwise phase orientation, two encircled by phases rotating in the opposite direction. It can be observed that the achieved far field configurations, that show two or four main points of phase singularity, can be produced by very simple dispositions of few Gaussian modes, provided that the relative phases and intensities are set with sufficient accuracy. It may also worth noting that the spiralling phase structure, which is easily observed around the phase singularity points in figures 1 or 2 here is not apparent, because it is dominated by the radial phase variation due to the curvature of the wavefront.
5 Propagation of the projected beam
Having shown the possibility to engineer the intensity and phase distributions in the far-field, it can be interesting to analyse the intensity pattern of the projected beams during their propagation. In fact, as we will show, by properly designing the waveguide position and phases it is possible to influence relevant properties of the propagating beam, such as the focus position.
We can start our analysis from the case of the polygonal arrays, already discussed in section 3 for what concerns the far-field. As already observed (see figure 2), if the waveguide phases are all equal to each other, the projected beam shows in the far-field an intensity pattern with a single peak in the middle, surrounded by a less intense ring, a distribution that is similar to that of a zero-order Bessel beam. A Bessel beam is usually realized by illuminating an axicon with a Gaussian beam or by placing an annular aperture in the back focal plane of a focusing lens . A similar configuration can be easily achieved through the use of few waveguides positioned on the vertices of a polygon. A field distribution with a spatial profile approaching that of a Bessel beam can indeed be obtained, as shown in figure 5, by exploiting an array with = 8 waveguides placed on the vertices of a regular octagon of radius = 25 m. Bessel beams are ideally supposed to be non-diffracting beams, a peculiar property that can be of interest in many application. However in practical cases they show this behaviour only in a limited region, the so-called Bessel zone, which can be estimated by calculating the propagation length where the intensity of the central peak is higher or equal to half of the maximum intensity.
In figure 5 (a) we report the evolution of the intensity of the central peak along the propagation direction. From the analysis of the curve we can estimate the length of the Bessel zone, which is approximately 2 mm long, i.e. much longer than the Rayleigh range of a Gaussian beam with the same wavelength and waist as the waveguides of the array, which is about = 0.28 mm. In figure 5 (b) we also report the transverse intensity pattern at specific propagation distances, that is, in the near-field at = 0 mm, at the coordinates where the axial intensity is equal to the maximum value, to its half and finally at a large distance, which can be considered in the far-field. It can be observed that in the first part of the propagation a reshaping of the intensity pattern of the array occurs, reaching a Bessel-like distribution with a central peak surrounded by rings where the intensity is maximum. Afterwards, the beam conserves its shape during the propagation, showing only a diffraction of the pattern. We calculated also the spot size parameter of the beam  at various distances, defined as where is evaluated along the or coordinate as:
The result is shown in figure 5, and confirms the beam divergence during the propagation. We can thus conclude that the near-field configuration considered in this example experiences a reshaping of the spatial intensity pattern with a sort of focusing after 1.3 mm of propagation.
It is well known in diffractive optics that, by governing the phase front in the near field, it is possible to control the intensity pattern along the propagation, e.g. to achieve focusing or to engineer the shape of the focus itself. As an example, phase shaping has also been used to make the peak intensity evolution of a Bessel-gauss beam uniform within the Bessel zone .
Here, we can think of moving the focus position by properly adjusting the position and phases of the waveguides in polygonal arrays, in particular considering the propagation of the field emitted by waveguides positioned on annuli with different radii and phase shifts. This configuration can be assimilated to that of a Fresnel lens, which is made of a series of concentric rings with different radii, each designed to diffract light so as to obtain constructive interference in the position where the focus is desired. This condition can be achieved by adjusting the phases of the waveguides in each ring according to the following equation:
where is the desired focal distance and the radial coordinate of the waveguide position. As an example, we consider a near-field pattern where the waveguides are distributed over two concentric rings with radii = 25 m and = 45 m. By adjusting the phases in such a way that the condition of constructive interference is satisfied for a focal distance of 2 mm, we obtain the translation of the propagation distance at which the intensity is maximum from = 1.3 mm, as obtained with only one ring, to = 2 mm. The results are shown in the top row of figure 6, where we report the transverse intensity pattern at different propagation distances.
Finally it is worth noting that, if one can separately act on the phase terms of a large number of waveguides, an efficient focusing can be obtained also when starting from a random distribution of the waveguides in the array, provided that the phases are properly set by the relation (10). An example of the intensity pattern that can be achieved at different propagation distances is reported in figure 6, where, starting from a near-field configuration of 16 waveguides randomly distributed in a 50 m radius circular area, a focal zone is obtained after 2 mm of propagation. To evaluate the focusing efficiency we calculated the spot dimension at 1/e and the ratio between the fraction of power confined in this region for the configuration in analysis, and for a Gaussian beam. The random configuration shown in figure 6 yields a focusing efficiency of 50%, higher than the efficiency of 28% that we achieved in the two-rings case.
6 Feasibility considerations and conclusions
Before concluding, we should consider whether the devices and waveguide configurations presented in the previous sections can be effectively fabricated, or our essay has to be considered a mere speculation. The three-dimensional capabilities of femtosecond laser writing technology certainly enables to produce arrays of waveguides with arbitrary positions in the cross-section [11, 13, 16]. However, it may be more difficult to divide the input light, e.g. from an optical fiber, into the waveguides of the array with the required different amplitude ratios and phases.
To this purpose, a very general scheme may be a generalisation of the one adopted in Ref. , where the integrated circuit is divided into three sections. The first section consists in a power dividing device, such as a star coupler or a cascade of directional couplers, to split the input light from the fiber into different branches with the required amplitude ratios. The second section consists in an array of phase shifters (e.g. thermo-optic devices) that allow to tune the phase terms of each waveguide. These sections may be fabricated in the plane and possibly with lithographic technologies. A third section transforms the waveguides from a planar arrangement to the required two-dimensional cross-section, bringing each waveguide to its final position. This needs to be realized by femtosecond laser writing due to the three-dimensionality of the layout. To minimise unwanted parasitic coupling between the waveguides, in the region where they are brought close together to produce the final setting, one might fabricate neighbouring waveguides with sufficiently different propagation constants, so that evanescent field coupling is quenched. With a higher technological effort, the first two sections may be integrated in a single fully-reconfigurable multiport interferometer that enables realising arbitrary unitary transformations between the input and the output waveguides . This would also permit to associate to each different input waveguide a different output configurations, provided that the chosen alternative output configurations are orthogonal supermodes.
As a matter of fact, the capabilities of femtosecond laser waveguide writing could allow reducing the three sections to a single three-dimensional interferometer, at least for specific applications. Let us consider, as an example, a possible way to produce different vortex configurations, as in figure 2, with the same physical array of waveguides. The different vortex configurations have all the waveguides excited with the same amplitude; the phase terms however are arranged cyclically, with different periodicity. These field configurations may be produced with a -inputs--outputs device that implements the Fourier transform of the input modes : by exciting a different input of the device, a different phase arrangement of the output would be produced. Compact integrated optical circuits performing this operation have been indeed recently demonstrated experimentally .
In the light of these considerations, our approach looks reasonably at reach of current technology. In particular, we have shown that it is possible to design waveguide arrays emitting at the output light beams with one or more phase singularity points, possibly carrying non-zero angular momentum. Focusing beams can also be produced, by adequately controlling the phase terms. In addition, the investigation of the far-field distribution here presented is immediately applicable to the distribution in the focal plane of a lens.
This work paves the way to other investigations that could deepen many aspects here just mentioned. For instance, for free-space communication purposes it will be required to understand how to optimize the array design to minimize cross-talks between the emitted modes . The capability to build multi-mode interferometric circuits that associate to each different input a different emitted supermode (e.g. with a different value of angular momentum), could be conveniently used as angular momentum or space division multiplexer/demultiplexer from a fiber optics networks [24, 22, 27]. In quantum optics applications, this would also allow to transfer entanglement from the waveguide path to the angular momentum or other spatial degree of freedom of the photons .
Benefits are also envisaged in optofluidics where lab-on-a-chips are used for single cell manipulation through optical forces . Indeed the integration of engineered waveguide arrays in optofluidic lab-on-a-chips could allow to implement tweezing or manipulation functionalities, as the control of the particle rotation in the optical trap for imaging and analysis purposes, which nowadays require bulk optics equipment. Unexpected applications might also be found in beam shaping for structured illumination, optical filtering or imaging tasks.
Appendix A Vanishing angular momentum within the array
We will show here that light propagating in a set of linearly-polarized and spatially separated guided optical modes cannot carry longitudinal angular momentum. The supporting dielectric structure and refractive index distribution will be considered as uniform along the axis, so that such axis is a well defined propagation direction and the optical mode profiles extend in the plane.
The propagating electromagnetic field of each mode, in each point, decomposed into its transverse and longitudinal components, takes the form:
The transverse Poynting vector, responsible for the longitudinal angular momentum, can be written as:
and its time average is:
If and are taken as real, it naturally holds that and . It is then easy to observe that, in each point . Since the transverse component of the angular momentum density, with respect to an arbitrary point O, is , in each point of the mode field its time average will be vanishing:
The overall angular momentum of the -th mode field would be , where the integration is performed over all the volume of the space and only the contribution of the -th waveguide field has been considered. It is clear from (14) that also the time average of the overall angular momentum is vanishing.
If we consider a set of well-separated modes, the overall angular momentum of the field is just , which has also vanishing time-average. The field of an array of well separated, linearly polarized, optical guided modes, thus possesses no angular momentum. Excitation of several overlapped modes is a necessary condition to have angular momentum within the array: in fact, in that case, since the modes are overlapped, the Poynting vector of the resulting field is not in general equal to the sum of the Poynting vector of the modes.
Appendix B Phase singularity in the far-field of a polygonal waveguide array
We derive here an asymptotic expression for the center of the far-field distribution of a polygonal waveguide array, as those described in section 3, valid in the limits of large propagation coordinate and small distances from the center.
We consider Gaussian modes, monochromatic and linearly polarized, all having their waist at = 0 coordinate, and centered on the vertices of a regular polygon: . All the Gaussian modes have the same intensity, but different phase terms, arranged in a cyclical configuration . The scalar field, at a certain plane , has the form :
where , , and is the curvature radius of the wavefront. It can be assumed for sufficiently large. is a normalization factor which includes common -dependent phase terms.
We start by expanding the squares in the arguments of the exponentials in (15), as , so that:
In the far field, i.e. for sufficiently large, , thus:
where all the constant terms have been included in .
We shall now focus on the central part of the beam and we seek for a power series limit of (17), valid for small . In particular, we consider:
We note that in the same limit, since we have also (and with much better approximation) .
This allows to substitute the last exponential in the sum with its power series expansion, and to approximate , i.e.:
We should now try to understand which is the first relevant power of , at which we can truncate the power-series expansion. It is known that, in general, where several may be vanishing, but . This means:
It is easy to show that is vanishing except in the case . Therefore:
In conclusion, the first relevant power of in (LABEL:somma), i.e. the smallest one that is multiplied by a nonvanishing coefficient, is . In fact, it is for that power of that first appear a sum of the kind (20) with . Hence, by collecting all constants in , we can write:
-  Wang J et al. 2012 Nature Photon. 6 488–96
-  Vallone G et al. 2014 Phys. Rev. Lett. 113 060503
-  Milione G et al. 2015 Opt. Lett. 40 1980–83
-  Bozinovic N et al. 2013 Science 340 1545–48
-  Van Uden R G H et al. 2014 Nature Photon. 8 865–70
-  Garcés-Chávez V et al. 2003 Phys. Rev. Lett. 91 093602
-  MacDonald M P, Spalding G C, Dholakia K 2003 Nature 426 421–424
-  Taylor M A et al. 2015 Nature Photon. 9 669–673
-  Bowman R W and Padgett M J 2013 Rep. Prog. Phys. 76 026401
-  Kreysing M K et al. 2008 Opt. Express 16 16984–16992
-  Szameit A, Blömer D, Burghoff J, Pertsch T, Nolte S and Tünnermann A 2006 Appl. Phys. B 82 507–12
-  Szameit A, Burghoff J, Pertsch T, Nolte S, Tünnermann A and Lederer F 2006 Opt. Express 14 6055–62
-  Rechtsman M C et al. 2013 Nature 496 196–200
-  Corrielli G et al. 2013 Nat. Commun. 4 1555
-  Crespi A et al. 2015 Phys. Rev. Lett. 114 090201
-  Caruso F, Crespi A, Ciriolo A G, Sciarrino F and Osellame R 2016 Nat. Commun. 7 11682
-  Fleischer J W, Segev M, Efremidis N K and Christodoulides D N 2003 Nature 422 147–50
-  Osellame R, Cerullo G and Ramponi R (Eds.) 2012 Femtosecond Laser Micromachining: Photonic and Microfluidic Devices in Transparent Materials (New York: Springer)
-  Svelto O 1998 Principles of Lasers, 4 edition (New York: Springer) p. 145
-  Goodman J W 2005 Introduction to Fourier optics (Roberts and Company Publishers)
-  Hong Z, Zhang J and Drinkwater B W 2015 Phys. Rev. Lett. 114 214301
-  Cai X et al. 2012 Science 338 363–66
-  Doerr C R and Buhl L L 2011 Opt. Lett. 36 1209–11
-  Su T et al. 2012 Opt. Express 20 9396–402
-  Liang Y, Wu H W, Huang B J and Huang X G 2014 Nanoscale 6 12360–65
-  Markin D M, Solntsev A S and Sukhorukov A A 2013 Phys. Rev. A 87 063814
-  Guan B et al. 2014 Opt. Express 22 145–56
-  Allen L, Beijersbergen M W, Spreeuw R J C and Woerdman J P 1992 Phys. Rev. A 45 8185
-  McGloin D and Dholakia K 2005Contemp. Phys. 46 15-28
-  Flamini F et al. 2015 Light Sci. Appl. 4 e354
-  Carolan J et al. 2015 Science 349 711–16
-  Durfee C G, Gemmer J and Moloney J V 2013 Opt. Express 21 15777–15786.
-  Barak R and Ben-Aryeh Y 2007 J. Opt. Soc. Am. B 24 231–40
-  Crespi A et al. 2016 Nat. Commun. 7 10469
-  Zhao N, Li X, Li G and Kahn J M 2015 Nature Photon. 9 822–26.
-  Fickler R, Lapkiewicz R, Huber M, Lavery M P, Padgett M J and Zeilinger A 2014 Nat. Commun. 5 4502
-  Yang T, Bragheri F and Minzioni P 2016 Micromachines 7 90