Optimal Wave Fields for Micro-manipulation in Complex Scattering Environments

TL;DR

This paper presents a method to optimize wave-fronts for micro-manipulation of objects in complex environments, achieving maximum force and focus efficiency without iterative optimization, verified through microwave experiments.

Contribution

It introduces a non-iterative eigenvalue-based approach to design wave-fronts for optimal micro-manipulation in complex scattering environments.

Findings

Successful demonstration in microwave experiments confirming theoretical predictions

Achieved maximum force and focus efficiency on targets of arbitrary shape

Method applicable to complex disordered media without iterative optimization

Abstract

The manipulation of small objects with light has become an indispensable tool in many areas of research ranging from physics to biology and medicine. Here we demonstrate how to implement micro-manipulation at the optimal level of efficiency for targets of arbitrary shape and inside complex environments such as disordered media. Our approach is to design wave-fronts in the far-field that have optimal properties in the near-field of the target such as to apply to it the strongest possible force, pressure or torque as well as to achieve the most efficient focus at the target position. Free of any iterative optimization, our approach only relies on a simple eigenvalue problem established from the scattering matrix of the system and its dependence on the target parameters. To illustrate this theoretical concept, we perform a proof-of-principle experiment in the microwave regime, which fully confirms our predictions.

Equations

$$S_{11}$$

$$S_{12}$$

$$Q_{n,11}$$

$$Q_{n,12}$$

$${\theta }_n^1$$

$${\theta }_n^2$$

$$I_1$$

$$I_2$$

$$\frac{I_1}{{\theta }_{\epsilon }^1}=\frac{I_2}{{\theta }_{\epsilon }^2}=\frac{2}{k}.$$

$$k_{\mathrm{in}}=\left(\begin{array}[]{cc}k&0\\ 0&-k\end{array}\right)=-k_{\mathrm{out}}.$$

$$-Q_{x,11}$$

$$-Q_{x,12}$$

$${\theta }_x^{1,2}$$

$${\vec{u}}_x^{1,2}$$

$$F_{1,2}$$

$$=\mp \frac{4\sqrt{2}\sin \left(2kLn\right)}{\sqrt{\Gamma }},$$

$$\frac{{\theta }_x^1}{F_1}=\frac{{\theta }_x^2}{F_2}=-\frac{k\left(n^2-1\right)}{2}.$$

$$Q_{L,11}$$

$$Q_{L,12}$$

$${\theta }_L^1$$

$${\theta }_L^2$$

$$P_1$$

$$P_2$$

$$S=-\text{mathbbm}1+2\mathrm{i}{V}^{†}GV,\mathrm{where}G={\left(\Delta +U(\vec{x})+\mathrm{i}V{V}^{†}\right)}^{-1}\mathrm{and}U(\vec{x})=k^2\epsilon (\vec{x}),$$

$$\frac{\mathrm{d}S}{\mathrm{d}\alpha }$$

$$-2\mathrm{i}{G}^{†}V{V}^{†}G=G-G^{†}.$$

$$Q_{\alpha }$$

$$=2\left(\text{mathbbm}1+2\mathrm{i}{V}^{†}{G}^{†}V\right)V^{†}G\frac{\mathrm{d}U}{\mathrm{d}\alpha }GV=2V^{†}\left(G-G+G^{†}\right)\frac{\mathrm{d}U}{\mathrm{d}\alpha }GV$$

$$Q_{\alpha }$$

$$|{\psi }_{\chi }⟩=2\mathrm{i}GV|\chi ⟩,$$

$$⟨\chi |Q_{\alpha }|\chi ⟩=\frac{1}{2}⟨{\psi }_{\chi }|\frac{\mathrm{d}U}{\mathrm{d}\alpha }|{\psi }_{\chi }⟩.$$

$$Q_k=4k{V}^{†}{G}^{†}\epsilon GV\equiv Q_d,$$

$$I\left(|u_{\epsilon }^i⟩\right)\equiv \int |\psi \left(|u_{\epsilon }^i⟩\right){|}^2\mathrm{d}A=\frac{2{\theta }_{\epsilon }^i}{k^2},$$

$$\frac{1}{2}⟨\psi ⟩\frac{\mathrm{d}U}{\mathrm{d}R}\psi$$

$$=\frac{k^2\left(\epsilon -1\right)R}{2}{\int }_0^{2\pi }{\left|\psi \left(\rho =R,\phi \right)\right|}^2\mathrm{d}\phi$$

$$=\frac{k^2\left(\epsilon -1\right)R}{2}P\left(|u_R^i⟩\right).$$

$$P\left(|u_R^i⟩\right)\equiv {\int }_0^{2\pi }{\left|\psi \left(\rho =R\right)\right|}^2\mathrm{d}\phi =\frac{2{\theta }_R^i}{R{k}^2(\epsilon -1)}.$$

Full text

See MainManuscript.pdf

Optimal Wave Fields for Micro-manipulation

in Complex Scattering Environments: supplementary material

M. Horodynski

Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria

M. Kühmayer

Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria

A. Brandstötter

Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria

K. Pichler

Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria

Y. V. Fyodorov

Department of Mathematics, King’s College London, London WC2R 2LS, UK

U. Kuhl

Université Côte d’Azur, CNRS, LPMC, 06108 Nice, France

S. Rotter

Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria

Corresponding author: stefan.rotter@tuwien.ac.at

1 Generalised Wigner-Smith Operators in One Dimension

To illustrate how our concept based on the generalised Wigner Smith operator works in a very simple context, we consider here the problem of scattering at a one-dimensional refractive index barrier (see Fig. S1) of refractive index $n>1$ ranging from $x=-L$ to $x=L$ (the refractive index in the asymptotic regions is $n_{0}=1$). For this simple problem we can explicitly calculate the scattering matrix $S$ and all the possible GWS-operators $Q_{\alpha}$ fully analytically. Specifically, the elements of the scattering matrix $S$ for this setup read as follows

[TABLE]

where $k$ is the wavenumber of the incident field. Note that in contrast to the main text all formulas here are written in terms of the refractive index $n=\sqrt{\varepsilon}$ instead of the dielectric function $\varepsilon$ in order to simplify the notation.

1.1 $Q_{n}$ for Refractive Index Variation (in 1D)

To better understand $Q_{\varepsilon}$, the GWS-operator associated to the total intensity inside a target, we first consider the GWS-operator we get when taking the derivative of the scattering matrix with respect to the refractive index $n$ of the barrier. At the end of this section we also show the results for $\alpha=\varepsilon$. The components of $Q_{n}=-\mathrm{i}S^{-1}\mathrm{d}S/\mathrm{d}n$ are

[TABLE]

The corresponding eigenvalues $\theta_{n}^{1,2}$ and eigenstates $\vec{u}_{n}^{1,2}$ then follow as

[TABLE]

The integrated intensities of the scattering states inside the barrier read as

[TABLE]

This result shows that there is a strict linear relation between the stored intensity in the designated scattering region and the eigenvalues of $Q_{n}$ which enables us to tune the intensity inside the scatterer based on the choice of $Q_{n}$-eigenstates (or a superposition of such eigenstates). Taking eigenstates with large eigenvalues results in scattering states that focus into the target region.

To calculate $Q_{\varepsilon}$ we make use of $\mathrm{d}/\mathrm{d}n=\mathrm{d}/\mathrm{d}\varepsilon\cdot\mathrm{d}\varepsilon/\mathrm{d}n=2n\mathrm{d}/\mathrm{d}\varepsilon$, which directly leads to the following expression

[TABLE]

1.2 $Q_{x}$ for Target Displacement (in 1D)

To construct $Q_{x}$ we take the derivative of the scattering matrix with respect to the position of the barrier. To obtain the eigenvalues associated to the momentum-transfer onto the target we make use of $Q_{x}=k_{\mathrm{in}}-S^{\dagger}k_{\mathrm{out}}S$ [1]. This formula comes from the fact that for only one scatterer (or if we move all scatterers) $Q_{x}$ measures the momentum difference between incoming and outgoing waves. The matrices $k_{\mathrm{in}}$ and $k_{\mathrm{out}}$ read as follows

[TABLE]

The components of $Q_{x}$ then read

[TABLE]

with $\Gamma\equiv-\left(n^{2}-1\right)^{2}\cos\left(4kLn\right)+n^{2}\left(n^{2}+6\right)+1>0\enskip\forall n,k,L$. The solution of the eigenproblem associated with $Q_{x}$ is

[TABLE]

where $N_{1,2}=|\vec{u}^{1,2}_{x}|$ is the norm of the eigenvectors. We then calculate the quantity $F_{1,2}$, which we will show to be proportional to the momentum transfer onto the target:

[TABLE]

where the first argument in the parentheses tells us which eigenvector of $Q_{x}$ is used in the calculation. We then compare these two quantities with the eigenvalues $\theta_{x}^{1,2}$ of $Q_{x}$ and see that

[TABLE]

This result tells us that there is a strict linear relation between the eigenvalues of $Q_{x}$ and the difference in intensities at the left and right boundary of the barrier, $F_{i}$. The eigenvalues $\theta^{1,2}_{x}$ are equal to the momentum difference between incoming and outgoing waves, i.e., $\Delta k=\theta_{x}$ [1]. Conservation of momentum then tells us that a momentum of $\Delta k$ is transferred to the target. We work in the stationary case, thus this momentum transfer onto the target is the same at all times. The average force, $F$, experienced by the target is the momentum transfer per unit time interval and since the momentum transfer is stationary, it is equal to the force. This implies that the quantities $F_{i}$ we have defined above are proportional to the force the target experiences from the scattering of the wave. This enables us to tune the degree of force exerted onto the target and even control the direction of it.

1.3 $Q_{R}$ for Radial Change (in 1D)

The one-dimensional analogue to the radius $R$ of a circle is the length $L$ of a region with raised refractive index. We therefore investigate the behaviour of $Q_{L}=-\mathrm{i}S^{-1}\mathrm{d}S/\mathrm{d}L$, whose elements are given by

[TABLE]

The solution of the corresponding eigenproblem yields

[TABLE]

The intensities at the boundaries of the barrier are

[TABLE]

for each eigenstate of $Q_{L}$. In contrast to $Q_{n}$ whose eigenvalues represent the integrated intensity over the target region, the above result shows that the eigenvalues of $Q_{L}$ correspond to the intensities at the target’s boundary which can be again tuned by choosing a certain $Q_{L}$-eigenstate (or a superposition of such eigenstates).

2 Derivation of Eq. 2

In the following we will sketch the derivation of Eq. 2 from the main text involving the expectation value of the generalised Wigner-Smith (GWS) operator with respect to some arbitrary input vector $\ket{\chi}$. All conventions are taken from [2] (see also [3, 4]). The scattering matrix $S$ for a system described by the Helmholtz equation can be written as

[TABLE]

where $G$ is the system’s Green’s function, $\Delta$ is the Laplacian, $U$ determines the scattering environment and $V$ describes the coupling to the leads. The derivative of the scattering matrix can then be written as,

[TABLE]

where we have assumed that only the scattering environment $U$ depends on $\alpha$, but not the coupling matrix $V$. Under the further assumption of a unitary $S$-matrix the following useful identity was derived in [2] for the Green’s function $G$,

[TABLE]

With a unitary scattering matrix we can then write $Q_{\alpha}=-\mathrm{i}S^{-1}\frac{\mathrm{d}S}{\mathrm{d}\alpha}$ as,

[TABLE]

The wavefunction inside the scattering system for an arbitrary input vector $\ket{\chi}$ can be written as [2],

[TABLE]

resulting in the following expression for the expectation value of $Q_{\alpha}$,

[TABLE]

This novel result now allows us to extract information about the local wave intensity inside the scattering medium out of the asymptotic information encoded in the scattering matrix $S$.

As a simple test of the above relation, we take the derivative of $U$ with respect to $k$ and we get

[TABLE]

which is exactly the result for the dwell-time operator as given in [2, 4]. Keeping in mind that the dwell time operator and the Wigner-Smith time-delay operator are the same for a vanishing $k$-dependence of $V$ [2, 4], we have thus shown that our new relation is perfectly consistent with this established case.

3 Linear Relations in Two Dimensions

Building up on the insights gained in the previous section, we can now prove the linear relations written down in Eqs. (3) and (4) of the main text. We start by proving the relation between $\theta_{\varepsilon}$ and the intensity integrated over the target’s area. Taking the derivative of the scattering environment with respect to the dielectric constant of a single scatterer results in $\mathrm{d}U/\mathrm{d}\varepsilon=k^{2}A_{\mathrm{scat}}\left(\vec{r}\right)$, where $A_{\mathrm{scat}}$ parameterises the area of the scatterer and is equal to one inside the target and zero outside. This then gives

[TABLE]

where the integral is performed over the area $A$ of the target scatterer. This equation is numerically confirmed in Fig. S2a and equips us with a tool to control the wave intensity inside a target’s area – ranging from no intensity to the theoretical maximum.

Next, we derive in detail the linear relation between the eigenvalues of $Q_{R}$ and the pressure, $P\left(\ket{u_{R}^{i}}\right)\equiv\int_{0}^{2\pi}\left|\psi\left(\rho=R\right)\right|^{2}\mathrm{d}\varphi$, applied to the target. We start by considering a dielectric target, for which $\mathrm{d}U/\mathrm{d}R=k^{2}(\varepsilon-1)\delta\left(\rho-R\right)$. This leads us to

[TABLE]

This derivation provides us with the final result that the pressure is in a linear relation to the eigenvalues of $Q_{R}$

[TABLE]

We confirm this equation numerically in Fig. S2c. To prove a linear relation between the eigenvalues of $Q_{R}$ and the transferred pressure in the case of a metallic target with perfect conductance we first consider the simple example of an infinitely long metallic rod of circular cross-section in vacuum. Effectively this configuration can be reduced to a two-dimensional problem. Due to the simplicity of the geometry, different incoming cylindrical wave modes do not mix, therefore the scattering matrix $S$, which in this case is just a unitary diagonal reflection matrix $r$, reads

[TABLE]

where $\mathrm{H}_{n}^{(1)}$ and $\mathrm{H}_{n}^{(2)}$ are the Hankel functions of the first and second kind, respectively. The knowledge of the reflection matrix then allows us to compute the GWS-operator $Q_{R}$, whose elements on the diagonal read

[TABLE]

where $J_{n}$ and $Y_{n}$ are the Bessel functions of the first and second kind, respectively. The wavefunction around the metallic rod comprised of the incident and scattered wave is

[TABLE]

where $b_{n}$ are the components of the input vector. This leads to the following expression in the case when $b_{n}=1$ for only one $n$ and zero for all other

[TABLE]

We thus find

[TABLE]

which we confirm numerically in Fig. S2c for a circular metallic target inside a waveguide featuring the same disorder as in Fig. 1 of the main manuscript. The reason why our derivation of the proportionality constant also works for an arbitrary environment as wells as an arbitrary cross-section, is that the proof of the GWS-operator tells us that its eigenvalues $\theta_{\alpha}$ are in a linear relation with the local wave intensity around the target, irrespectively of the surrounding scattering environment and whether these are waveguide walls or a disordered medium.

Next we consider the eigenvalues of $Q_{x}$, the GWS-operator we get when the parameter $\alpha$ considered is the longitudinal position $x$ of a target. This problem was originally considered in [1], although a proof was only given for the case of longitudinally moving the entire scattering system. Here we close this gap and prove a linear relation between the eigenvalues and the momentum transferred also for a single target inside a disordered medium, using the same strategy as in the preceding paragraph. For a dielectric circular target with radius $R$, the result reads

[TABLE]

where the direction of the shift and the corresponding momentum transfer is generalised to an arbitrary direction, which is parametrised by the unit vector $\hat{\vec{n}}$ and the integral is performed along the boundary of the circular target (a generalization to arbitrary target shapes is also possible). In the presence of a metallic target the result is

[TABLE]

Numerical results confirming these relations for $\hat{\vec{n}}=(1,0)^{T}$ can be found in Fig. S2b.

The same approach also works when considering the eigenvalues of $Q_{\varphi}$ that are in a linear relation with the torque transferred to the target. In the case of a dielectric target this relation reads

[TABLE]

where the integral is taken along the target’s boundary described by the curve $\mathcal{C}$ and parametrised by $\vec{c}$. The expression $\vec{n}\left(\vec{c}\right)\left|\psi\left(\vec{c}\right)\right|^{2}$ denotes the normal force at every point of the boundary excised by the electric field and $\vec{m}_{\perp}\left(\vec{c}\right)$ is the part of the distance from the boundary to the target’s center of mass that is normal to $\vec{n}$, i.e., the lever. We find that in the presence of a metallic target, the eigenvalues $Q_{\varphi}$ are proportional to

[TABLE]

where $\vec{n}\left(\vec{c}\right)\left|\partial_{\vec{n}}\psi\left(\vec{c}\right)\right|^{2}$ denotes the normal force at every point of the boundary excised by the electric field. These two equations are numerically verified in Fig. S2d for the case of a square target.

4 Applying Pressure with the GWS-Operator

In this section we look at the GWS-operator $Q_{R}=-\mathrm{i}S^{-1}\mathrm{d}S/\mathrm{d}R$, where $R$ is the radius of a metallic circular scatterer (our approach also works for more complicated target shapes). For the corresponding eigenvectors $\ket{u_{R}^{i}}$ we analytically showed in the preceding section that there exists a linear relation between the radiation pressure applied to such a metallic target and the corresponding eigenvalue $\theta_{R}^{i}$ of the following form:

[TABLE]

where the integral is along the boundary of the circular scatterer with radius $R$ and $\psi$ is the electric field distribution of the corresponding eigenstate $\ket{u_{R}^{i}}$. In other words, the eigenstates $\ket{u_{R}^{i}}$ of $Q_{R}$ are characterised by a well-defined radiation pressure $P$ that they apply to the dielectric target scatterer. The pressure applied to a chosen target can thus be maximised by injecting eigenstates $\ket{u_{R}^{\mathrm{max}}}$ corresponding to the largest eigenvalues $\theta_{R}^{\mathrm{max}}$. The experimentally measured scattering state of the eigenstate associated to the largest eigenvalue can be seen in Fig. S3a, where the experiment was carried out in a waveguide featuring a disorder (see Fig. 1). We can clearly see that the state has a strong intensity build-up on all sides of the boundary of the chosen scatterer. Constructing $Q_{R}$ with a dielectric scatterer (rather than a metallic one) as a target yields the following linear relation between the radiation pressure $P$ and the eigenvalues $\theta_{R}^{i}$,

[TABLE]

5 Comparing the GWS-Method and the Field Matrix Method

Here we verify explicitly that our micro-manipulation protocol is indeed optimal by comparing it to a procedure proposed in [5]. There the authors achieve a focus onto a point inside a scattering medium by utilizing the field matrix $e(x)$, which relates the field at a depth $x$ inside the scattering medium to the incident field in the waveguide’s leads. The components of the field matrix $e_{ab}(x)$ connect the field in channel $b$ at depth $x$ with the field in channel $a$ at the waveguide’s leads, i.e., $E_{b}(x)=e_{ba}(x)E_{a}$, where $E_{a}$ and $E_{b}(x)$ are the fields in channels $a$ and $b$. Optimal focus, i.e., maximizing $\lvert E_{\beta}(x)\rvert$, at a target point $\beta$ can now be achieved by shaping the incident wavefront as

[TABLE]

with $I_{\beta}(x)=\sum_{a}\lvert e_{\beta a}(x)\rvert^{2}$ (for a proof see [6]). In Fig. S4a we compare the focus achieved by the field matrix method with the one achieved with the highest $Q_{\varepsilon}$-eigenstate for a circular target scatterer with a very small diameter $D$ (a very small target allows us to converge to the limit of focusing on a point, as considered in the field-matrix method). Our comparison, indeed, shows that the highest $Q_{\varepsilon}$-eigenstate is indistinguishable from the one constructed by the field matrix method, i.e., they both deliver the optimal focus in form of a single peak in the center of the scatterer. In contrast to the field matrix method which, however, can only focus onto a single point, the GWS-concept also enables optimal focusing into an extended area of arbitrary size. Focusing into an extended area requires more than one peak such that the field matrix method would then already have to know where exactly all these peaks have to lie, which is equivalent to already knowing the state that leads to optimal focus in the first place. To demonstrate that the GWS-concept is able to focus in an extended area we show in Fig. S4b the spatial intensity distribution of the highest $Q_{\varepsilon}$-eigenstate around a target scatterer with a larger area. In order to maximise the stored intensity, $Q_{\varepsilon}$ constructs the highest eigenstate such that two peaks fit inside the scatterer, resulting in a two-fold increase of the stored intensity as compared to the field-matrix solution with just one peak in the center of the target.

6 Incomplete Access to Scattering Channels

In the main document we explore different realisations of GWS-operators for the ideal case where we have access to the full scattering matrix, $S$. Since it is experimentally challenging to measure the entire scattering matrix, we numerically show that our approach also works in the regime of a sub-unitary $S$-matrix. For this case we consider the same waveguide ($L=6W$) as in the main document but fill it with $300$ randomly distributed Teflon scatterers ($R=0.025W$, $n=1.44$) and a square metallic target (side length $0.165W$) at $x$-position $L/6$. Instead of the full scattering matrix, we now take only the reflection matrix $r$ for the calculation of the GWS-operator, since in the experiment both transmission and reflection measurements are not always possible. Experiments also suffer from a low numerical aperture (LNA), which we take into account in our simulation by removing the $M$ highest modes from $r$ such that we are left with a $(N-M)\times(N-M)$ matrix $r_{\mathrm{LNA}}$, where $N$ is the number of modes. When we use only the reflection matrix to construct an operator $q_{\varphi}=-\mathrm{i}r^{-1}\mathrm{d}r/\mathrm{d}\varphi$ its eigenvalues turn out to be $\theta_{\varphi}^{i}=\left(M_{z}(\vec{u}_{\varphi}^{i})+\mathrm{i}\vec{u}_{\varphi}^{i,\dagger}t^{\dagger}\mathrm{d}t/\mathrm{d}\varphi\vec{u}_{\varphi}^{i}\right)/(2\langle r^{\dagger}r\rangle)$ [1], where $\langle r^{\dagger}r\rangle$ is the global reflectance of the eigenstate and the loss of information manifests itself in complex eigenvalues caused by the second term containing the transmission matrix $t$. From this expression we may thus conclude that our procedure based on the reflection matrix only will work better for higher total reflection and for the case that the target is placed at the beginning of the waveguide (placing it on the end would result in far less intensity reaching it). In order to still maintain a correlation between the eigenvalues of $q_{\varphi}$ and the total torque acting on the target scatterer, we also have to restrict ourselves to channels that are strongly reflecting to minimise the term with the transmission matrix $t$. We achieve this by using a singular-value decomposition of $r_{\mathrm{LNA}}=U\Sigma V^{\dagger}$, where the matrices $U$ and $V$ contain the column wise the left and right singular vectors (not to be confused with the scattering environment and the coupling matrix, respectively) and the matrix $\Sigma=\mathrm{diag}\left(\left\{\sigma_{n}\right\}\right)$ contains the singular values on its diagonal. In order to project on the highly reflecting channels we pick a certain subset of large singular values $\tilde{\Sigma}=\mathrm{diag}\left(\left\{\tilde{\sigma}_{n}\right\}\right)$ as well as the associated left and right singular vectors $\tilde{U}$ and $\tilde{V}$. Equipped with these matrices we can construct an effective inverse $r_{\mathrm{LNA}}^{-1}=\tilde{V}\left(\tilde{U}^{\dagger}r_{\mathrm{LNA}}\tilde{V}\right)^{-1}\tilde{U}^{\dagger}$ and also project the derivative of $r_{\mathrm{LNA}}$ onto this subspace with the proper projection operators $P_{\tilde{U}}=\tilde{U}\tilde{U}^{\dagger}$ and $P_{\tilde{V}}=\tilde{V}\tilde{V}^{\dagger}$ [1]. As a final result we obtain:

[TABLE]

In Fig. S5a one can see that Eq. (53) is almost perfectly fulfilled, i.e., our protocol works even for the case where only parts of the scattering matrix are available. This, however, comes with the caveat that the system needs to be highly reflecting (transmitting) when we work only with the reflection (transmission) matrix. Fig. S5b shows the spatial intensity distribution of the $\tilde{q}_{\varphi}$-eigenstate which transfers the maximal amount of torque.