Optomechanical deformation and strain in elastic dielectrics

Light forces induced by scattering and absorption in elastic dielectrics lead to local density modulations and deformations. These perturbations in turn modify light propagation in the medium and generate an intricate nonlinear response. We generalise an analytic approach where light propagation in one-dimensional media of inhomogeneous density is modelled as a result of multiple scattering between polarizable slices. Using the Maxwell stress tensor formalism we compute the local optical forces and iteratively approach self-consistent density distributions where the elastic back-action balances gradient- and scattering forces. For an optically trapped dielectric we derive the nonlinear dependence of trap position, stiffness and total deformation on the object's size and field configuration. Generally trapping is enhanced by deformation, which exhibits a periodic change between stretching and compression. This strongly deviates from qualitative expectations based on the change of photon momentum of light crossing the surface of a dielectric. We conclude that optical forces have to be treated as volumetric forces and that a description using the change of photon momentum at the surface of a medium is inappropriate.


Introduction
As light carries momentum besides energy, its propagation through a polarizable medium is accompanied by forces. Although the momentum of a single light quantum is very small, laser light can generate appreciable forces on the microscopic scale. Optical forces are nowadays routinely used to manipulate and trap particles ranging from single atoms and molecules [1,2,3] to plastic beads, biological cells or microbes up to the size of tens of micrometres [4,5,6]. The mechanical motion of even larger objects such as silica mircodisks or suspended mirrors has been damped and cooled by light forces [7,8]. While most of the existing work targets the overall effect on the centre of mass of the particles, it has been shown by us as well as by other groups that these forces do not act homogeneously but exhibit distinct patterns within the medium [9,10,11]. For any elastic medium this leads to local compression or stretching. Of course the modified density also changes the local refractive index and light propagation, which b a l a n c e s c h a n g e s again alters the forces as displayed schematically in figure 1. The resulting coupled complex evolution thus obviously requires self-consistent models and solutions [12]. In addition, as the light mediated interaction is inherently long range, even a small but periodic variation of the refractive index can have a very large overall collective effect coupling distant areas over a large volume.
This work is organised as follows: In section 2 we first present the basic scattering approach to treat the light propagation in an inhomogeneous refractive medium and use a previously developed formalism based on the free space Maxwell stress tensor to calculate the corresponding local force distribution (section 3). This method is then used in section 4 to develop an iterative scheme to calculate the steady state density and field distribution as a function of geometry and field intensity. In section 5 we discuss essential physical consequences predicted by our model at the hand of numerical examples. Finally in section 6 these results are set against common calculations of the total deformation at hand of the change of photon momentum at an interface between to dielectrics.

Multiple scattering model of light propagation in inhomogeneous media
The effective light propagation in a medium can be seen as the result of multiple individual scattering processes, which in general requires intricate numerical treatments, if one cannot make use of material symmetries. Here we restrict ourselves to the simple but still nontrivial case of two incoming counterpropagating plane waves in a transversely homogeneous and linearly polarizable medium. In this limit only forward and backward scattering add up phase coherently, while all amplitudes for transverse scattering average out. From the viewpoint of the forward and backward propagation directions, transverse scattering thus can just be added to an effective absorption rate in the medium. This is certainly not perfectly fulfilled in an actual setup, but still can be expected to give the correct qualitative behaviour, as long as the transverse extensions are much larger than the wavelength of the light. A more realistic treatment, e. g. in terms of Gaussian transverse beams, is possible, but greatly complicates the model and would obscure many interesting physical phenomena found in this simple approach.
Restricting the dynamics to the forward and backward scattering amplitudes along the propagation directions gives a simple and tractable model for our medium via a one dimensional array of N thin slices at positions x 1 , . . . , x N . Here the spatial behaviour of the electric field E(x, t) = Re[E(x) exp(−iωt)]e y is determined by a 1D Helmholtz equation [13,14,11] (∂ 2 The field-induced polarisation density at each slice then is P (x) = αη A E(x) N j=1 δ(x − x j ), where we introduced the dimensionless coupling parameter ζ = kη A α/(2ε 0 ) proportional to the atomic polarizability α and the areal particle density η A within the slice. ε 0 is the vacuum permittivity and k = ω/c the wave number of the optical field. Note that we assumed here that the dipoles in each slice can simply be added up coherently for scattering along the propagation direction. As illustrated in figure 2, the equation above is satisfied by interconnected plane wave solutions [13], here denoted as for x j < x < x j+1 . The amplitudes left and right of a material slice (beam splitter) at position x j are connected via The amplitudes (A j , B j ) and (C j−1 , D j−1 ) are coupled by a simple propagation matrix, Therefore, the amplitudes to the left of the (j + 1) th slice are obtained by a simple multiplication of the previous transfer matrices, The amplitudes A 1 and D N are determined by the amplitudes and phases of waves coming in from the left (i. e. −∞) and from the right (+∞), respectively, and constitute boundary conditions on the solutions for the Helmholtz equation (1). B 1 and C N are obtained by computing the total reflection and transmission amplitudes via Note that the reflection coefficients for left or right incidence on an inhomogeneous setup usually do not coincide, i. e. r l = r r , but the transmission amplitude t is independent of the direction of propagation. More details on the properties of these generalised transfer matrices are given in Appendix A.  . In this case, the field generated by multiple scattering by the beam splitters (blue curve, cf. (2)) reproduces the solution for a homogeneous medium with refractive index n (red curve), if the coupling is chosen as in (6). On the right hand side we see the displaced medium with irregularly spaced slices of the same coupling ζ and the resulting electric field. The background shading illustrates the change in the distances, i. e. the strain u ′ (x) = −( ρ(x) − ρ)/ρ, with dark colours indicating regions of higher density.
For equally spaced, thin polarizable slices we set x j = (j − 1)d 0 , such that x 0 = 0 and x N = (N − 1)d 0 =: L and (4) simplifies to (A j+1 , B j+1 ) T = T h j (A 1 , B 1 ) T , with T h := P d 0 M BS . In an earlier work [11] we showed that choosing a uniform distance d 0 between the slices and setting the coupling parameter leads to the same optical fields as found inside a medium with refractive index n. A sufficiently dense array of beam splitters with spacing d 0 = L/(N − 1) is then in the limit N → ∞ indistinguishable from a homogeneous medium of refractive index n and length L. A decisive step in this work, which allows us to account for local material density variations, is the introduction of a local displacement variable u(x) x j → x j = x j + u(x j ), j = 1, . . . , N.
As illustrated schematically in figure 2, such shifts alter the local fields as well as the total reflection and transmission properties of the object. The distances between the slices then change as A continuous limit can be consistently defined via u( In analogy with the theory of elastic deformations, we call u ′ strain or deformation [15,16], and the relative change in the initially homogeneous local material density ρ simply reads Let us here comment on the notation we will use for the rest of this work. As defined in the paragraph above (6), our coordinates are chosen such that the unperturbed medium occupies the region [0, L]. Introducing a displacement u then shifts the object to [u(0), L + u(L)], with L + u(L) − u(0) =: L. But to ease notation, all the quantities such as the electric field strength E or force F shall remain defined with respect to the original position such that e. g. E(0) [E(L)] always marks the field at the left [right] edge of the medium. The amplitudes at the boundaries then have to be adjusted with corresponding phases, cf. (22). This, however, is relevant for mathematical formulations only, physical discussions and figures are unaffected by this detail. In (7) we introduced a tilde to distinguish the shifted x j from the original x j . For most other quantities such as the fields or forces, we will omit this tedious notation. Only the changed length L, the inhomogeneous density ρ (10) and refractive index n (19) still have to be distinguished from their original values L, ρ and n.
As mentioned before, defining the coupling ζ as in (6) ensures that the solutions of the wave equation (1) agree with the field inside a homogeneous dielectric at positions x j = (j − 1)d 0 , if the fields are assumed to agree at x 1 . In the continuous limit N → ∞, the latter requirement is always fulfilled [11]. Interestingly we still preserve this feature for a model with displaced slices, if we choose the following approach: Let, as in (2), E j (x) denote the plain wave solution of the Helmholtz equation (1) and denote a field defined in the same region, but with a refractive index n j . To obtain the desired equivalence between a stratified dielectric and a set of irregularly spaced slices, we assume for any given j ∈ {1, . . . , N − 2} and demand that with E(↑y) ≡ lim x↑y E(x), The first line shows the conditions that E j and E j+1 are solutions of the Helmholtz equation (1), cf.
(3) or [13], the second line denotes Fresnel's equations for the transition between two dielectrics. In the third line, finally, we demand that the plane wave solutions of (1) agree with the fields inside the dielectrics at positions x j+1 and x j+2 . This leads to the required, successive coupling between ζ, the distances d j = x j − x j−1 and d j+1 , and some indices n j , n j+1 . Solving (13)-(15) under the assumption (12) and demanding solutions independent of the field amplitudes results in two conditions, for j = 1, . . . , N − 1 One can easily check that these conditions give the known relation (6) in the equidistant case where d j = d j+1 ≡ d 0 and n j = n j+1 ≡ n. Unfortunately, we were not able to find solutions with finite values of d j = d j+1 for both conditions. Inspired from (6) one may try to find that this approach satisfies (17), but not (16). However, choosing ζ as in (6), writing d j = d 0 (1 + ∆ j /d 0 ), and taking the continuous limit N → ∞ with ∆ j /d 0 → u ′ (x) (9) alters (18) to satisfying both (16) and (17). With the given inhomogeneous refractive index we can compute the electric field inside a strained, one dimensional dielectric by solving numerically. A comparison with the field computed via the transfer matrix method described in (4) shows excellent agreement, for sufficiently large N.
Another way to approximate the optical field is to expand the transfer matrices in (4) for small local deformations ∆ j and then perform the continuous limit. This analytical approximation works sufficiently well for the typically small strain u ′ obtained in the scope of parameters used in this work. The lengthy results of this approach are presented in Appendix B, equation (B.6).
Inserting the relation between strain and density modifications (10) we finally obtain where we assumed ( ρ − ρ)/ρ ≪ 1 for the final expansion.

Computing the reflection and transmission amplitudes
To find solutions for the fields inside the medium with refractive index distribution n(x), one needs to specify initial values. As discussed for the discrete system in (5), the medium can be described in terms of a transfer matrix such that B 0 = r l A 0 + tD 0 and C L = tA 0 + r r D L , if the electric fields outside the medium are given as E(x) = A 0 exp(ikx)+B 0 exp(−ikx) for x ≤ 0 and E(x) = C L exp(ik(x−L))+D L exp(−ik(x−L)) for x ≥ L. The amplitudes A 0 and D L are determined by the intensities I l,r and phases φ l,r of the fields incident from the left and right and the displacement u, as Therefore, the initial conditions for solutions of (20) are But obviously, the reflection and transmission coefficients r l , r r and t strongly depend on the refractive index n(x). To calculate those one can either use some approximations, cf. Appendix B, equation (B.5), or solve the field equation (20) for specially chosen boundary values, e. g.
allowing the easy computation of r l , r r and t.
, then a beam entering from the left experiences the same medium as one from the right and hence r l = r r . Note that for the homogeneous case where u ′ = 0 and n = n, we recover the usual [17] t h = 2n 2n cos(nkL) − i(n 2 + 1) sin(nkL) ,

Light forces in an inhomogeneous medium
In general, the total electromagnetic force on an object embedded in vacuum is given by [18] where A denotes the surface of the object, n is the normal to A and T αβ is the Maxwell stress tensor Using two planes orthogonal to the direction of propagation (i. e. the x-axis) as integration surfaces and the plane wave fields defined in (2), the time-averaged optical force per area (pressure) on the j th slice simply reads [19] 51i. Please note that in this case the strain is not chosen such that it balances the optical forces as in (39).
Following the beam splitter relation in (3) we rewrite C j = (1 + iζ)A j + iζB j and Taking the naive limit lim N →∞ F j would give a vanishing force per slice as lim N →∞ ζ = 0, cf. (6) with d 0 = L/(N − 1). But assigning each slice to one N th of the object's total length L we can define a force density F (x) := lim N →∞ NF j /L and use Following the derivation of the inhomogeneous refractive index (19) we can replace is a solution of the wave equation (20). Hence we obtain the local optical force density where we used the algebraic limit theorem for lim N →∞ N|ζ| 2 /L = 0. Again, the results from the formula above are fully consistent with forces computed using a large but finite number of slices (4) and (27) as well as an analytical approximation presented in Appendix B, equation (B.9).
In figure 3 we compare the intensity and optical forces in a medium with homogeneous refractive index n to fields and forces obtained by solving (20) and (30), respectively, for a given strain u ′ .

Identification of radiation pressure and dipole force components inside a homogeneous dielectric
In the case of a medium with uniform refractive index n, i. e. with u ′ ≡ 0, the force density computed from (30) can be identified with established expressions for optical forces on dielectric test particles. There, the time-averaged force on a dipole at position where α is the polarizability of the dipole. The first term proportional to Reα is often referred to as dipole or gradient force, the (dissipative) term proportional to Imα is called radiation pressure or scattering force [20].
Inside a homogeneous dielectric we may write the spatial component of the electric where the amplitudes G and H are chosen such that Fresnel's conditions at the boundary of the object are met. Using this field and rewriting (31) into a force density on particles with volume density η V located at x 0 leads to As the coupling parameter in (1) is defined as ζ = kη A α/(2ε 0 ), with η A denoting the areal particle density in each of the N slices, we may write using (6) and (29) This relation also resembles the Lorentz-Lorenz relation for the case of a thin gas [17], where the individual dipoles do not directly interact with each other. It can easily be checked that the force density in (32) together with the Lorentz-Lorenz relation (33) gives exactly the same result as the force computed from (30), if we insert the same field for a homogeneous dielectric. This demonstrates that our approach to calculate fields and forces from multiple scatterers is consistent with well known results derived from more general assumptions.

Integrated force and trap formation
To compute the total force on an extended dielectric in a standing wave we can use the same derivation as for the force on an infinitesimal slice in (27) and get with A 0 and D L denoting the amplitudes at the object's left and right boundaries, as given in (22). Defining the position of the centre of the object x 0 = (u(0) + L + u(L))/2, we can express the total force in terms of ξ : If 4I l I r v 2 ≥ (I l s l − I r s r ) 2 , the force vanishes at every position ξ 0 ∈ Ξ + ∪ Ξ − where Using some trigonometric properties one can show that stable trapping positions are those defined in the set Ξ + . For ξ 0 ∈ Ξ − we find F ′ tot (ξ 0 ) > 0 and hence Ξ − is a set of unstable trapping position. Linearising the total force in (35) around stable trapping positions ξ 0 ∈ Ξ + leads to a trap stiffness of Since the reflection and transmission coefficients strongly depend on the object's size, one finds that the parameter ψ can change its sign abruptly for certain values of L. This leads to sudden jumps between a low-and high-field-seeking behaviour for a trapped object [21,22,23,11]. Figure 5 shows how trap position and trap stiffness is changed by the strain induced on the object by optical forces.

Self-consistent balancing of optical force and elastic back-action
In the previous chapters 2 and 3 we found expressions for the local fields and the local optical force densities in deformed, dielectric media. But depending on the given elastic properties, the strain will result in stress which typically tries to compensate the external volumetric forces.
In this chapter we will investigate the behaviour of a linear elastic, dielectric object subjected to the optical forces described by (30). More precisely, we will provide a framework to compute the equilibrium configuration between the optical forces and the elastic counter reaction in a self-consistent manner. In our computations we will assume only optical forces and neglect thermal or piezoelectric effects as well as surface tension.
Mechanical equilibrium between some general volume force density f and the resulting stress denoted by the tensor σ is given by Cauchy's equilibrium equation [16] for i, j denoting the coordinates of the system. Since the model discussed here considers only one relevant dimension, this equilibrium equation simplifies to f + ∂ x σ = 0. The constitutive relation for a linear elastic, one dimensional object simply reads σ = Eu ′ , where E is Young's modulus and u ′ is the local strain [16].
Hence we see that an equilibrium between the optical force density and the elastic strain requires at every position x ∈ [0, L], with F being a solution of (30). Note that the electric field computed from (20) also depends on the amplitudes at the edges of the object and therefore also on the displacement u, cf. (22).
Solving (39) for an equilibrium requires boundary conditions on the displacement u and the strain u ′ which are determined by the given setup. The displacement is fixed by the assumed trapping mechanism, e. g. if the object is trapped by a standing wave, we have to fulfil (u(0) + L + u(L))/2 ∈ Ξ + . But note that ξ 0 depends on the reflection and transmission coefficients and hence also on the deformation u ′ , cf. (36).
The strain has to be chosen such that the stress σ = Eu ′ at each surface balances external surface forces [16]. Assuming for the moment an object subjected to volumetric optical forces only, we integrate the equilibrium equation (39) at obtain For an object trapped by light fields, we get F tot = 0 and σ(0) = σ(L) = 0, due to the lack of surface pressure. In chapter 5.2, however, we fix the slab by an external mechanism balancing the total optical force via surface interaction. Hence if the left boundary of the slab is retained at x = 0 (i. e. u(0) = 0), then σ(L) = 0 and σ(0) = F tot .
To solve equation (39) numerically, we use an iterative approach where the equilibrium condition is rewritten in the form with u i and u ′ i denoting the displacement and strain obtained by the i th iteration step and F (x)[u i , u ′ i ] is the force density computed using u i and u ′ i . The updated u ′ i+1 and u i+1 can then be obtained by simple numerical integration, with integration constants chosen in accordance to the boundary conditions of the used setup. With the updated optical force densities F (x)[u i+1 , u ′ i+1 ] one can compute the next step of the iteration. An obvious choice for initial values is a homogeneously shifted distribution (i. e. a constant u 0 ) with a given starting length L and refractive index n.
This iterative scheme proved to be sufficiently exact but significantly faster compared to other methods of solving nonlinear equations, such as Newton's method. We also confirmed our computations with force densities obtained from the transfer matrix approach in (4) and the analytic approximation described in Appendix B, respectively.

Examples and physical interpretation
In our basic considerations above we always assumed the object to be exposed to two counterpropagating laser beams of the same, linear polarisation forming a standing wave. These results can easily be extended to describe situations with only one incident beam or with two counterpropagating beams of different polarisation. In the latter case, one has to calculate the intensities and forces separately for each polarisation direction.
In this section we will present four showcase examples to give insight into the large variety of possible results. The first two examples deal with the case where an object is trapped by two counterpropagating beams. The latter examples will treat the case where the object is illuminated by only one beam and externally fixed at one end. For For all considered setups we will see that the interaction between optical forces and elastic back-action strongly depends on the ratio between the initial length L and the wavelength of the deformation beam in the unperturbed medium, λ/Re[n]. Concerning the intensity and the elastic properties of the medium we find that all results grow linearly in (I l + I r )/(Ec), at least in the scope of parameters where a solution of the equilibrium equation (39) could be obtained within a reasonable error tolerance. That is why the local intensities, I(x) = ε 0 c|E(x)| 2 /2, and force densities are given in units proportional to E in the upcoming figures. Note that the numbers used in the simulations are unrealistic in order to exaggerate the effects, since an intensity of I = 0.1Ec would imply I ≃ 30W/µm 2 for E ∼ 1MPa.

Example: Object trapped by two counterpropagating beams
As argued above, an object with initial length L subjected to optical forces will in general experience local deformations and an overall length change. Figure 4 shows the relative length change ∆L = ( L−L)/L for different initial lengths L in the cases where an object is trapped in a standing wave (blue lines) or by two beams of orthogonal polarisation (red curves) and equal intensity. Obviously, both configurations are symmetric regarding an  inversion of x at the centre of the object and therefore we find r l = r r , as mentioned in section 2.1. Surprisingly, in the standing wave case we observe abrupt switching from strong compressive to stretching behaviour around certain initial lengths. A comparison with chapter 3.2 and earlier discussions in [11] shows that these switches are concurrent with jumps of the stable trapping position ξ 0 . Generally speaking, objects with small values of kL are trapped at local maxima of the intensity in the standing wave. But for larger objects, the term ψ in (36) abruptly changes its sign and the object seeks positions centred around intensity minima. As indicated by the dotted grid lines, these jumps occur at lengths of minimal reflection |r h | 2 and maximal transmission |t h | 2 (24), i. e. at L = mλ/(2Re[n]), m ∈ N.
For larger objects we observe a slight decay of the maximal relative length change, which is found also for computations with Im[n] = 0 and thus is not caused by additional radiation pressure only. However, a glance at the reflectivity and transmission shows that the self consistent deformation prevents configurations with zero reflectivity which would usually result in maximum elongation or compression. For the given standing wave trap we generally observe that deformations computed with the steady state equation (39) tend to increase reflectivity and decrease transitivity compared to a homogeneous medium, for both Im[n] = 0 or Im[n] = 0.
In figure 5 we compare the trap stiffness κ from (37) and the trapping position ξ 0 from (36) between unperturbed and self-consistently deformed objects. We find that the deformation significantly increases the trap stiffness, even if the total size of the object remains unchanged, cf. figure 4. An elastic object in a standing wave therefore assists in enforcing its own trap, just like a rabbit who starts to dig when captured in   Figure 6. Dielectric objects trapped in a standing wave configuration with I l = I r = 0.05Ec. The initial lengths are chosen as L ≃ 0.19λ (a) and L ≃ 0.4λ (b); n = 1.3 + 0.0025i. As mentioned in (10), positive strain u ′ (x) (green curves) denotes larger distances between the particles and hence a reduction of material density. The density is increased at positions where the force density (blue lines) F is zero and F ′ < 0, i. e. where the light intensity (depicted red) has a maximum value. As indicated by the marks in figure 4 and verifiable from the signs of u ′ (x), the object on the left gets squeezed whereas the right hand side example shows a stretched object.
a pit. Note that for a compressible slab in a standing wave there is no critical length of zero trap stiffness, but deformation always leads to a stable trapping position. Two examples for optical force densities and the associated deformation in the standing wave setup are shown in figure 6. There we see that the strain is negative (i. e. the material density is increased) at positions where the force density changes from positive values, denoting a force pushing to the right, to negative values associated with a force pushing to the left. One can clearly see the difference between the compressive situation (a) where the object is trapped at maximal intensity and situation (b), where the trapping occurs at minimal intensity and the object experiences a stretching force.
To trap a dielectric slab with two non-interfering plane waves of orthogonal polarisation, the intensities of said beams have to be equal, i. e. I l = I r . In this case there are certain starting lengths L ≃ (2m + 1)/(4Re[n]), m ∈ N 0 , for which the intensities inside the unperturbed object add up to a constant value. Hence for these specific lengths the dominant gradient forces add up to zero and the object's length remains unchanged. Apart from that we find an expanding behaviour for L > λ/(4Re[n]), as can also be seen in the examples in figure 7.
In the case of non-interfering beams, figure 4 also shows that the reflectivity and transitivity is no longer periodic in the object's length L. For the reflectivity we find that the zeros are shifted from L = mλ/(2Re[n]), m ∈ N, towards smaller lengths and no longer coincide with lengths of maximal elongation.   Figure 7. Dielectric objects trapped by two counterpropagating beams with orthogonal polarisation and I l = I r = 0.05Ec. As in figure 6, the initial lengths are L ≃ 0.19λ (a) and L ≃ 0.4λ (b); n = 1.3 + 0.0025i. As can be seen from figure 4, the length of the first example is chosen such that the intensities (red curves) inside the medium almost add up to a constant value. Hence the force densities (blue curves) and the self-consistent deformation (green lines) practically vanish (please note the changed scale on the axes). For the second example we obtain a stretching of ∆L ≃ 6.6 × 10 −3 . Figure 8 shows the relative length change in the case where the left edge of a dielectric slab is fixed by some external mechanism. Here we observe a striking difference whether the object is illuminated from the left (blue) or right hand side (red curves). In the first case we find only stretching behaviour with minor oscillations of the relative length change ∆L. However, if the object is illuminated from the right hand side (i. e. the beam is incident on the free surface) we again find both compression and elongation, depending on the initial length L. Furthermore we observe a stronger impact of absorption than for the previous case with two counterpropagating beams. Now we see that if Im[n] = 0, then ∆L increases (light incident from the left, radiation pressure pushing the object to the right) or decreases (light incident from the right) for larger initial lengths L. As we see from the dash-dotted lines, the length change continues to oscillate around a constant value also for large values of L, if Im[n] = 0.

Example: Object fixed at the left boundary and illuminated by one beam.
For the reflectivity [transitivity], the self consistent strain again results in a shift of the minimal [maximal] values to the left of L = mλ/(2Re[n]), m ∈ N. Note that for the given, asymmetric setup r l = r r , but the difference in the total reflectivity |r l | 2 − |r r | 2 remains very low for the given parameters.
To explain the different results for left and right incident beams, we take a closer look at the interface where a single beam exits the slab, e. g. at hand of figure 9. Let us assume a beam entering from the left and I r = 0. Then the intensity is not only constant    Figure 9. Two dielectric objects fixed at the left edge and illuminated from the left (a) or from the right hand side only (b), with I l,r = 0.1Ec. As predicted from (42), the intensity (red) is locally maximal and the force densities (blue lines) are (almost) zero at the interface where light exits the slab. For left incidence we observe only positive (i. e. stretching) strain u ′ (green lines) whereas both elongating and compressive deformation is feasible for right illumination.
on the right of the object, but Fresnel's formulae also tell us (for a homogeneous object) d dx I(L) = 0 and d 2 and hence the intensity has a local maximum at the edge where the beam exits the slab. For the usually dominant gradient force F gr (x) ∼ d dx I(x) we thus find F gr (L) = 0 and F gr > 0 left of the surface. Considering the steady state equation (39) and omitting the scattering force shows that and hence the strain u ′ has a local minimum at the right edge x = L. If light enters only from the right hand side we find analogous behaviour for the left edge x = 0. Since minimal strain corresponds to a maximum in local material density (10), we conclude that the gradient force tends to accumulate material at the surface where a single light beam exits the object. A similar statement holds for the aforementioned case with two counterpropagating beams oscillating in orthogonal polarisations, but then one has to add up the forces generated by the two beams.
In total, the constraints on the strain as derived from (42), (43) and chapter 4 are found as left incidence only, I r = 0 : right incidence only, So if I l = 0, the strain is fixed at a minimum with negative value on the left edge, at zero on the right edge, and oscillates proportional to the intensity in between. So in total we find both negative (compressive) and positive (expanding) deformation, depending on the length and refractive index of the object. For I r = 0 we see a positive strain at the left edge and a minimum with u ′ (L) = 0 at the right boundary. Hence the deformations oscillate between zero and some positive value and always lead to a total stretching behaviour.

Identifying length changes by probing the reflectivity
As one can see from figures 4 and 8, the relative length change ∆L obtained for a given ratio of total intensity to Young's modulus can be imperceptibly small, especially if the original length is not chosen in an optimal relation to the vacuum wavelength of the trapping beam, λ. But one possibility to detect minor stretching or squeezing for arbitrary initial lengths can be found in the use of a second, weaker laser probing the change in the reflectivity of the medium: Assuming a non-dispersive medium, a probe laser with a vacuum wavelength matching the Fabry-Pérot condition will travel through a slab of length L without being reflected, i. e. |r h (λ p )| 2 ≃ 0. Turning on a powerful laser will deform the object and hence also change the reflectivity for the weak probe beam. Figure 10 shows the relative change of reflectivity between the unperturbed and the strained medium for different wavelengths of the probe laser beam. We can see that as λ p crosses values defined in (46), δR changes from negative to positive values if the object is compressed (i. e. ∆L < 0) and vice versa if ∆L > 0.

Estimating the deformation by computing the photon momentum transfer on a surface
There exist numerous experimental and theoretical papers reporting optical stretching of deformable objects, such as biological cells [24,25], or light induced outward bending of liquid-gas surfaces [26,27]. In the mentioned publications, the light-induced deformation is estimated by considering an effective photon momentum change at the transition from one medium to another. In this context, the optical forces emerge as surface forces, acting on the interface between two regions of different refractive index. Since the considered materials are incompressible, the refractive index in each region remains constant. Let us try here a similar approach to estimate the deformation of an elastic object and put our results in context with these earlier works. But note that the very different physical properties of a linear elastic medium as compared to incompressible water, plane waves instead of Gaussian beams and a wave description instead of geometric optics, do not allow a straightforward comparison of the results. But nevertheless we can investigate whether our findings based on a volumetric description of optical forces are compatible with a concept of surface forces due to photon momentum exchange.
Following the line of [24,25,26] one estimates the time averaged force per area on an interface separating two regions with different indices of refraction n 1 = n 2 ∈ R and fields E(x) = A exp(in 1 kx) + B exp(−in 1 kx), x ≤ 0, and E(x) = C exp(in 2 kx), x > 0, as Here I inc = n 1 cε 0 |A| 2 /2 is the total energy-flux density entering the system, R = |B| 2 /|A| 2 and T = n 2 |C| 2 /(n 1 |A| 2 ) is the reflected and transmitted fraction the energyflux and p i = kn i describes the momentum of a single photon in a medium with index n i , i = 1, 2. As in [24,25,26] we here used Minkowski's version of the momentum of light in dielectric media. Out of curiosity about simple but puzzling arguments on stretching or compression of media in connection with the Abraham-Minkowski controversy [12,26,28,29] we also included results computed by naively inserting Abraham's result for the photon momentum, p i = k/n i , in figure 11. One must note, however, that Abraham's stress tensor would also require a material tensor component. A thorough calculation should always give the same results, independent of the used version of stress tensor [28,29]. Neglecting internal reflections, the force on an extended object with index n 2 embedded in a medium n 1 and subjected to a single beam is estimated to give F tot = F 1,2 +T F 2,1 . The deformation of such an object is then simply the difference in the two forces on the surfaces, reading F def = T F 2,1 −F 1,2 . Assuming a linear elastic medium with Young's modulus E, the relative length-change can be estimated by ∆L = F def /E. For the values used in the previous examples n 1 = 1, n 2 = 1.3, I 0 = 0.1Ec, we obtain ∆L ≃ 0.059 when using Minkowski's momentum and ∆L ≃ −0.046 for Abraham. These deformations have about the same order of magnitude as our full self-consistent computations, depicted e. g. in figure 8, but do not depend on the length of the object.
Formally one can refine these calculations from (48) and include also light incident from the right such that E(x) = C exp(in 2 kx) + D exp(−in 2 kx), x > 0 to obtain For D = 0 this reduces to (48) and for n 1 = n 2 = 1 and p 1 = p 2 = k we recover the force derived previously with the Maxwell stress tensor (27). This allows to formulate a generic wave optics extension for the deformation estimated above in the scope of geometric optics, now including also size dependent reflection and transmission. The total deformation pressure on an object with length L and homogeneous refractive index n ∈ R surrounded by vacuum then reads outside are connected by the homogeneous reflection and transmission coefficients (24), B = r h A + t h D, C = t h A + r h D, and the incoming A and D are given in (22). As expected, a similar calculation for the total force F tot gives the same result as we obtained previously in (34). The resulting relative length change ∆L = F def /E is presented in figure 11. There we find that the estimations using optical surface forces even qualitatively differ from the results we obtained with the present description using the full, volumetric optical forces, cf. figures 4 or 8, even if the force on the surface is adapted to include interference due to internal reflections.
An intuitive example is the object of length L = λ/(2n) where r h = 0 and t h = 1. In a standing wave trap with I l = I r this object is then trapped at ξ 0 = x 0 + (φ l − φ r )/(2k) = 0, cf. chapter 3.2 or figure 5. Using Minkowski's p n = kn we then find that the force due to photon momentum transfer vanishes at each surface and hence ∆L = 0 in figure 11.
But from the examples in figure 6 we deduce that the intensity at the surface is zero, yet the object will contract due to the intensity maximum at its central position. This is because the dipole force pulls each volume element towards the next local maximum of intensity, regardless of whether this volume element is located at the surface or in the bulk of the medium.
In fact, none of our calculations or simulations showed any distinctive effects suggesting a surface force at the boundaries of a dielectric [11]. This is supported by computations on large but finite stacks of polarizable slices where the forces on the first or last slice qualitatively do not differ from those on the second or next to last, respectively. We therefore conclude that optical forces have to be treated as real volumetric forces [12,30] and that a description using the change of photon momentum at the surface of a medium is inappropriate, regardless of using Abraham's or Minkowski's momentum.

Conclusions
Using an implicit calculation of optical fields and forces allows to self-consistently determine the stationary local deformations of an elastic object, where the local stress balances the local light forces by elastic back action. These solutions show a surprisingly variable and nonlinear dependence on the chosen parameters. Generally we see a length and illumination dependent, spatially quasiperiodic strain pattern, which can lead to length stretching as well as compression. As expected, standing wave configurations yield the strongest forces and effective length changes with a clear resonant structure for special ratios of initial object length L and trap beam wavelength λ. In the standing wave setup, variations in the trap wavelength λ lead to discrete jumps of the stable trapping positions. At L = mλ/(2Re[n]), m ∈ N, the particle switches from a position centred around an intensity maximum to one around a field node, which is associated with changes from compression to elongation of the object. Interestingly, in particular close to these instability points, this generally leads to an increase in trap stiffness. We expect that this indicates possible bistability between high and low field seeking behaviour for certain lengths very close to L = mλ/(2Re[n]).
Although the calculations presented here were performed in the scope of elastic media, we believe that the model can be extended to deformable but incompressible media like water or even dilute gases. Here in particular stability thresholds for the homogeneous solutions should prove physically very interesting, as they could lead to stationary flows, periodic density oscillations or light induced density pattern formation and particle ordering in a gas.
Here we limited our considerations to the case, where a steady state solution can be found. As for other nonlinear dynamical effects [14], there are of course regions in parameter space, were no stationary solutions exist and we find self sustained oscillations or even disintegration of the material. Indications of this behaviour appear e. g. in a non converging iteration procedure. At this point we leave this to future work.
we conclude that (T , ·) is a group.
To motivate definition (A.1), let us consider two plane waves E l (x) = A exp(ik(x))+ B exp(−ik(x)), x ≤ 0, and E r (x) = C exp(ik(x − L)) + D exp(−ik(x − L)), x ≥ L, left and right of a dielectric with a homogeneous refractive index n. Hence their amplitudes are coupled as where S n 1 ,n 2 couples the amplitudes at the intersection from a region with refractive index n 1 to a region with index n 2 and P d denotes the propagation matrix over a distance d, i. e.
S n 1 ,n 2 = 1 2n 2 n 2 + n 1 n 2 − n 1 n 2 − n 1 n 2 + n 1 and P d = e ikd 0 0 e −ikd . (A.4) One can easily show that S n,1 · P nL · S 1,n = T(r h , r h , t h ), with r h and t h as given in (24). But the attempt to construct a total transfer matrix for stacked media such as S n K ,1 ·P n K d K ·S n K−1 ,n K · · · S n 2 ,n 1 ·P n 1 d 1 ·S 1,n 1 with arbitrary n i , d i , i = 1, . . . , K, will show that one reflection coefficient is not enough. However, such a system can be described by a more general T (r 1 , r 2 , t) ∈ T such that B = r 1 A + tD and C = r 2 D + tA.

Appendix B. Analytical approximations for electric fields and forces for small deformations
In section 2 and in an earlier work [11] we showed that for equally spaced slices the amplitudes of the electric fields are related as (A j+1 , B j+1 ) T = T h j (A 1 , B 1 ) T , with T h := P d 0 M BS . Choosing the coupling ζ as in (6) Allowing small local density variations d j − d 0 = ∆ j we may use P d j ≃ P d 0 + ik∆ j σ z P d 0 to expand the relation between the field amplitudes from (4) as where σ z = diag(1, −1). In the limit of infinitely many slices within a finite length L, the sum above can be rewritten to an integral and we obtain A(x) B(x) ≃ T(x) + ik Here r h and t h denote the reflection and transmission amplitudes of a homogeneous medium with length L and refractive index n, as given in (24). As mentioned in section 2.1 we find that symmetric local strain, i. e. u ′ (x) = u ′ (L − x), x ∈ [0, L], results in r l ≃ r r and antisymmetric strain gives t ≃ t h and r l − r h ≃ r h − r r .
Using the reflection and transmission amplitudes from (B.5) to replace B(0) = r l A(0) + tD(L) in (B.4) leads to analytical approximations for the local field amplitudes inside a deformed medium A(x) ≃ a h (x) + A 0 (f (L − y)f (y − x) + g(L − y)g(y − x))+ +D L (f (y)f (y − x) − g(y)g(y − x)) , b x (x, y) = −ik n(f (L) + g(L)) A 0 (f (L − y)f (x − y) − g(L − y)g(x − y))+ +D L (f (y)f (x − y) + g(y)g(x − y)) , a L (x, y) = ik(f (x) − g(x)) n(f (L) + g(L)) 2 A 0 (f (L − y) 2 − g(L − y) 2 )+ +D L (f (y)f (L − y) + g(y)g(L − y)) , b L (x, y) = ik(f (x) + g(x)) n(f (L) + g(L)) 2 A 0 (f (L − y) 2 − g(L − y) 2 )+ +D L (f (y)f (L − y) + g(y)g(L − y)) . (B.7) The functions f and g are given in (B.2), the amplitudes at the boundary A 0 ≡ A(0) and D L ≡ D(L) are determined by the incoming intensities and the displacement and can be read off from (22). To obtain an approximation for the forces we use (28), take the limit F (x) := lim N →∞ NF j /L with lim N →∞ Nζ/L = k(n 2 − 1)/2 and insert the amplitudes from (B.6) to find with the terms a X and b X as given above in (B.7). Just like the field amplitudes in (B.6) and the reflection and transmission coefficients in (B.5), the above result is only an approximation for the case of small deformation u ′ . That is, we neglect products of ∆ l ∆ m in (B.3) and hence also correlations of type u ′ (y 1 )u ′ (y 2 ) . . . dy 1 dy 2 . But a comparison of the approximated results in the continuous limit with solutions of the wave equation (20) or numerical computations for a large but finite number of beam splitters confirmed that this first order expansion is sufficient for the scope of parameters used in this work.