# Detecting and Focusing on a Nonlinear Target in a Complex Medium

Antton Goïcoechea, Jakob Hüpfl, Stefan Rotter, François Sarrazin, Matthieu Davy

► To cite this version:

Antton Goïcoechea, Jakob Hüpfl, Stefan Rotter, François Sarrazin, Matthieu Davy. Detecting and Focusing on a Nonlinear Target in a Complex Medium. Physical Review Letters, 2025, 134 (18), pp.183802. 10.1103/PhysRevLett.134.183802. hal-05080222

HAL Id: hal-05080222

https://hal.science/hal-05080222v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Detecting and Focusing on a Nonlinear Target in a Complex Medium

Antton Goïcoechea,¹,* Jakob Hüpfl,² Stefan Rotter,² François Sarrazin,¹ and Matthieu Davy¹,†

¹Université de Rennes, CNRS, IETR - UMR 6164; F-35000 Rennes, France

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

(Dated: April 14, 2025)

Wavefront shaping techniques allow waves to be focused on a diffraction-limited target deep inside disordered media. To identify the target position, a guidestar is required that typically emits a frequency-shifted signal. Here we present a noninvasive matrix approach operating at a single frequency only, based on the variation of the field scattered by a nonlinear target illuminated at two different incident powers. The local perturbation induced by the nonlinearity serves as a guide for identifying optimal incident wavefronts. We demonstrate maximal focusing on electronic devices embedded in chaotic microwave cavities and extend our approach to temporal signals. Finally, we exploit the programmability offered by reconfigurable smart surfaces to enhance the intensity delivered to a nonlinear target. Our results pave the way for deep imaging protocols that use any type of nonlinearity as feedback, requiring only the measurement of a monochromatic scattering matrix.

*Introduction.*—Wavefront shaping techniques can par-
tially counteract the effect of disorder by coherently con-
trolling wave-matter interaction [1–3]. Of particular in-
terest is the possibility to focus waves on a diffraction-
limited focal spot inside or behind a strongly scattering
medium [4] or to deposit energy to a target region [5, 6].
In the linear regime, when the field at the target location
is directly accessible, the incident wavefront for focus-
ing can be optimally tailored in space and/or time using
techniques such as phase-conjugation for monochromatic
waves [7–10], time reversal for broadband signals [11] or
the eigenstates of an operator constructed from the scat-
tering matrix [6, 12, 13]. However, embedding a detector
within a scattering medium is an invasive procedure that
rules out many applications in deep optical imaging, wire-
less communications, wireless power transfer, or sensing.
Noninvasive approaches therefore rely on the presence of
a guidestar within the medium [3, 9, 14–16].

The nonlinearity of wave-matter interaction has
emerged in this context as an efficient approach for deep
imaging. Nonlinear techniques rely on the localized feed-
back generated by a nonlinear target. In acoustics, mi-
crobubbles serve as contrast agents for ultrasound imag-
ing [17, 18] while in optics, Raman microscopy [19],
two-photon fluorescence [20] or second-harmonic genera-
tion [21, 22] have been exploited to obtain a diffraction-
limited focal spot [16, 23]. In the microwave regime, most
electronic devices, even as simple as a diode [24–26], ex-
hibit a nonlinear behavior and can be detected by har-
monic radars in cluttered environments [27, 28]. All these
techniques nevertheless require complex experimental se-
tups to detect and/or filter the frequency-shifted nonlin-
ear signal.

A local perturbation of a linear scattering medium can
also serve as a guidestar [29–35]. Any change within a
disordered sample is encoded in the random speckle pat-
tern resulting from the complex interaction of the inci-
dent wave and the sample [36]. Therefore, the derivative

of the scattering matrix $S(\omega)$ with respect to a parame-
ter $\theta$, i.e. $\partial_{\theta}S$, contains information about $\theta$. For unitary
scattering matrices, the eigenvalues of the generalized
Wigner-Smith (WS) operator $Q = -iS^{-1}\partial_{\theta}S$ indicate
how strongly the conjugate quantity to $\theta$ is affected by a
variation [29, 30, 37]. The operator $(\partial_{\theta}S)^{\dagger}\partial_{\theta}S$ also turns
out to measure the content of Fisher information carried
by the scattered wave on the parameter $\theta$ [32, 38, 39].
The eigenstates of these operators thus provide the so-
lution for maximal focusing, micromanipulation or for
optimal sensitivity with respect to $\theta$. However, setting
up these operators requires a variation of the target pa-
rameter(s) in the first place, which is hard to accomplish
in linear static systems without invasive external inter-
vention inside the scattering medium.

Here, we present a noninvasive approach for optimal fo-
cusing on a nonlinear target in static scattering systems.
Most importantly, our approach does not require a mea-
surement of a higher-harmonic response of the target de-
vice or any other frequency-shifted signal. This allows for
the detection of all types of nonlinearities (whether they
induce n-th order harmonic generation, or a Kerr effect
for example), and it does not rely on prior knowledge on
the medium or on the target. All these requirements are
satisfied by leveraging the nonlinear scattering response,
which we probe by tuning the incident power [22, 40].

*Theoretical analysis.*—We consider a system made up of dielectric obstacles described by the dielectric function $\epsilon(\mathbf{r})$, which we probe by coherent monochromatic electromagnetic waves at frequency $\omega$. Our objective is to find an incident field $E_{\omega}^{\text{in}}$ that creates an optimal focus on a target described by the polarization $\mathbf{P}^{\text{NL}}(\mathbf{E})$, which is nonlinear in the electric field $\mathbf{E}$, embedded inside this scattering environment.

To disentangle the linear (L) from the nonlinear (NL)
scattering response, we write the far-field solution $E_{\omega}^{\alpha}(\mathbf{r})$
at frequency $\omega$ for an incident field strength $\alpha$ using the
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Photograph of the experimental setup. A resonator coupled in the near field to a wire antenna connected to a (nonlinear) low-noise amplifier (LNA) is located within a two-dimensional multimode waveguide (cover plate not shown). The scattering matrix is measured using seven antennas at the left interface. An absorbing foam is placed at the right interface to mimic open boundary conditions. (b) Reflection parameter $r = |S_{11}|^2$ of the LNA as a function of incident power $P$ for an empty single-mode waveguide at $f_0 = 12.68$ GHz with the second port connected to the LNA (powered ‘on’ – or not – ‘off’) terminated by an open-circuit condition. (c,d) Spectra of the singular values $\lambda_n$ of $\Delta S$ (c) and asymmetry factor $\mathcal{A}$ (d). The nonlinearity is maximal at the resonance of the nonlinear target $f = f_0$ indicated by a black-dashed line. (e,f) Intensity map within the system for a random illumination (e) or with $c^{\text{in}} = d^*$ (f). Both maps are normalized by the maximum value obtained for the optimized wavefront.</figcaption></figure>

following exact integral equation:

$$ \mathbf{E}_{\omega}^{\alpha}(\mathbf{r}) = \alpha \mathbf{E}_{\omega}^{\text{L}}(\mathbf{r}) + \frac{k^2}{\epsilon_0} \int_{\mathbb{R}^3} d\mathbf{r}' G_{\omega}(\mathbf{r}, \mathbf{r}') \mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha}(\mathbf{r}')] \quad (1) $$

Here, $\mathbf{E}_{\omega}^{\text{L}}(\mathbf{r})$ describes the linear component of the field for normalized amplitude $\alpha = 1$ and $G_{\omega}(\mathbf{r}, \mathbf{r}')$ is the Green's tensor of the system without nonlinearities. The target's nonlinear polarization field component at frequency $\omega$ is given by $\mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha}(\mathbf{r}')]$ and depends in general on all frequency components of the electric field $\mathbf{E}_{\omega_1}^{\alpha}(\mathbf{r}'), \mathbf{E}_{\omega_2}^{\alpha}(\mathbf{r}'), \dots$ at $\mathbf{r}'$. However, for easier notation we write this in Eq. (1) as a dependence on the full field $\mathbf{E}^{\alpha}(\mathbf{r}) = \int_{0}^{\infty} \text{Re}(\mathbf{E}_{\omega}^{\alpha}(\mathbf{r})e^{-i\omega t})d\omega$.

The idea is now to eliminate the linear term and to time-reverse the nonlinear signal to create a focus at the target (see Supplementary Material [41] for details). We achieve this by varying the amplitude of the incident field by $\Delta\alpha$, then the difference of the amplitude-normalized fields $\delta\mathbf{E}_{\omega}(\mathbf{r}) = (\alpha + \Delta\alpha)^{-1}\mathbf{E}_{\omega}^{\alpha+\Delta\alpha}(\mathbf{r}) - \alpha^{-1}\mathbf{E}_{\omega}^{\alpha}(\mathbf{r})$ satisfies

$$ \delta\mathbf{E}_{\omega}(\mathbf{r}) = \frac{k^2}{\epsilon_0} \int_{\mathbb{R}^3} d\mathbf{r}' G_{\omega}(\mathbf{r}, \mathbf{r}') \delta\mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}'). \quad (2) $$

We see that $\delta\mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}') = (\alpha + \Delta\alpha)^{-1}\mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha+\Delta\alpha}(\mathbf{r}')] - \alpha^{-1}\mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha}(\mathbf{r}')]$ acts as a source term for $\delta\mathbf{E}_{\omega}(\mathbf{r})$ with the coupling from the polarization field to the far field being given by the Green's tensor $G_{\omega}(\mathbf{r}, \mathbf{r}')$. Thus, as

long as the incident field interacts with the nonlinearity, the time-reversed of $\delta\mathbf{E}(\mathbf{r})$ provides a focus from the far field onto the nonlinear target.

This becomes especially apparent for a point-like nonlinearity located at $\mathbf{r}_0$, where the field difference provides the exact coupling between the location of the nonlinearity to the far field, i.e. $\delta\mathbf{E}_{\omega}(\mathbf{r}) \propto G_{\omega}(\mathbf{r}, \mathbf{r}_0)$. If we further assume that only a single electromagnetic mode couples to the nonlinearity, then the focus that is created at $\mathbf{r}_0$ with this method is necessarily optimal for a given incident strength $\alpha$. Note that throughout the article, ‘optimal’ and ‘maximal’ refer to states obtainable in the linear regime because the nonlinearity is assumed to be weak (perturbative regime).

From now on, we focus on scalar waves. For linear systems, the linear scattering matrix $S^{\text{L}}$ connects all input and output channels $c^{\text{in}}$ and $c^{\text{out}}$ through $c^{\text{out}} = S^{\text{L}}c^{\text{in}}$. For nonlinear systems, the superposition principle and therefore this linear relation no longer hold since the measured output field in $c^{\text{out}}$ must encapsulate the coupling between the incident field emitted in $c^{\text{in}}$ and the nonlinearity. This nonlinear response emanates from the position of the nonlinearity $\mathbf{r}_0$ and scatters into the output channels quantified by the coefficient vector $\mathbf{d}$ in accordance with the Green's function $G_{\omega}(\mathbf{r}, \mathbf{r}_0)$. Our goal is to split the incident field into a linear part that does not interact with the nonlinearity and a nonlinear part. The Green's identity tells us that all incident fields $c^{\text{in}}$ that
are orthogonal to the phase conjugate of $d$ (i.e. $d^*$) have vanishing electric field at the location of the nonlinearity $E_{\omega}^{L}(r_0) = 0$. This means that the nonlinear contribution to the output field must be a function of $d \cdot c^{\text{in}}$:

$$c^{\text{out}} = S^{\text{L}}c^{\text{in}} + \alpha^{-1}df(d \cdot c^{\text{in}}), \quad (3)$$

where the nonlinearity of the system is captured by the arbitrary nonlinear scalar function $f$ (see Supplementary Material [41]).

If we now probe the system at a given input power $\alpha$, the resulting input-output relation can be expressed through the following power-dependent scattering matrix

$$S_{m,n}^{\alpha} = S_{m,n}^{\text{L}} + \alpha^{-1}d_m f(\alpha d_n), \quad (4)$$

where $S_{m,n}^{\alpha}$ corresponds to the output in channel $m$ for an input in channel $n$. However, it is important to note that $S^{\alpha}$ is not a conventional scattering matrix. While it describes the scattering for the set of incident fields at a given power, it cannot encompass the full nonlinear nature of Eq. (3).

Nevertheless, $S^{\alpha}$ can now be used to identify the optimal focusing input by using the difference matrix

$$\Delta S_{m,n} = S_{m,n}^{\alpha+\Delta\alpha} - S_{m,n}^{\alpha} \\ = d_m \{(\alpha + \Delta\alpha)^{-1}f[(\alpha + \Delta\alpha)d_n] - \alpha^{-1}f(\alpha d_n)\}. \quad (5)$$

The rank of the matrix is equal to the number of point-like nonlinearities (here one), each coupling to the far field through a unique incident wavefront. $d$ can now be estimated by applying a singular value decomposition (SVD) on $\Delta S = U\lambda V^*$. The left singular vector $U_1$ of the largest singular value $\lambda_1$ corresponds to $d$. By applying a phase-conjugation, the incident wavefront $d^*$ provides maximal focusing onto the target.

The nonlinear scattering coefficient of the target depends on the field $E_n(r_0)$ transmitted by channel $n$. In the absence of global symmetries, we have $E_n(r_0) \neq E_m(r_0)$ in general for $n \neq m$. In contrast to linear systems with modulation of an antenna impedance or dielectric permittivity of a subwavelength object [30, 31, 34, 35], the right and left singular vector of $\Delta S$ are therefore not equal.

**Experimental results.**—In our first experiment, we consider a nonlinear target embedded within a two-dimensional scattering system working in the microwave range. The two-dimensional geometry is chosen so that the field can be measured ‘inside’ the system. The target consists of a high-Q dielectric cylinder with a refraction index of $n \approx 6$ coupled in its near field to a wire antenna connected to a low-noise amplifier (LNA). This nonlinear passive target - the LNA is not powered - is conceptually similar to the coupling of a resonator to a short-circuited diode found in Refs [24–26] for nonlinear coherent perfect absorption or the formation of defect modes. However,

we find that our target exhibits a more pronounced nonlinearity due to the presence of transistors in the LNA. This is crucial since the target is excited from its far field. We characterize the reflectivity of the LNA by connecting it to a single-mode waveguide whose reflection parameter $|r(\omega)|^2 = |S_{11}|^2$ is measured for increasing input power $P$. While $|r(\omega)|^2$ is constant at low powers indicating a linear behavior, we observe in Fig. 1(b) that $|r(\omega)|^2$ rapidly decreases for $P > 0$ dBm before saturating for $P > 10$ dBm. The nonlinearity therefore results from enhanced absorption within the LNA at high power.

The nonlinear target is then placed within a multi-mode waveguide (cavity) supporting a single mode in its vertical direction between 12 GHz and 15 GHz [42, 43]. Seven metallic and 13 dielectric scatterers are randomly placed inside the cavity to randomize the field. Placing the nonlinear device in a disordered cavity highlights the performance of the approach in cases where the Green’s functions within the system are unknown. The flux-normalized matrix $S^{\alpha}(\omega)$ is measured using a vector network analyzer (VNA) between $N = 7$ single-mode waveguides that are connected to the left interface of the cavity by coax-to-waveguide transitions. At the right interface, we place an absorbing foam to mimic open boundary conditions. The difference matrix $\Delta S(\omega)$ is constructed from measurements of $S^{\alpha}(\omega)$ at two incident powers $P_1 = 0$ dBm ($S = S^{\text{L}}$) and $P_2 = 12$ dBm.

Two peaks in the spectrum of the first singular value $\lambda_1(\omega)$ of $\Delta S$ are observed at 12.68 and 13.35 GHz in Fig. 1(c), corresponding to resonances of the high-Q dielectric cylinder. The enhancement of the field intensity within the cylinder at resonance results in a strong coupling with the LNA, and thereby increases nonlinear effects. We then scan the normalized field $t_n(x, y; \omega)$ inside the medium for each source channel $n$ by translating a short wire antenna in holes drilled into the top plate of the cavity. These measurements allow us to reconstruct the intensity map $\mathcal{I}(x, y; \omega)$ for any arbitrary incident wavefront $c^{\text{in}}$, $\mathcal{I}(x, y; \omega) = |c^{\text{in}*}(\omega) \cdot t(x, y; \omega)|^2$. For the phase-conjugate of the first singular vector $c^{\text{in}} = d^*$, a strong enhancement of the intensity at the resonance $\omega_n$ is observed in Fig. 1(f) compared to random illumination in Fig. 1(e), since $d$ gives the vector of Green’s functions between the sources and the target. The intensity at the focus at 12.68 GHz is enhanced on average by a factor $\eta \approx 4.7$ relative to a random incident wavefront.

Whereas we probed the system at two incident powers in order to apply a focus onto a target, we can instead exploit the reciprocity of linear systems for the detection of the nonlinearity. For linear systems, the scattering matrix is symmetric; however, the nonreciprocity induced by the nonlinearity breaks the symmetry of $S^{\alpha}$ [44–46], i.e.

$$(S^{\alpha} - (S^{\alpha})^T)_{m,n} = d_m \alpha^{-1} f(\alpha d_n) - d_n \alpha^{-1} f(\alpha d_m). \quad (6)$$

Thus, detecting the presence of a nonlinearity in practice
<figure><img src="image_3.png" /><figcaption>FIG. 2. Experimental results for a powered LNA connected to a single mode coax-to-waveguide transition located at the right interface of the cavity. (a) Spectrum of the singular values $\lambda_n$ of $\Delta S$. (b) Temporal signal at the nonlinear target corresponding to the back-propagation of the reconstructed signal for maximal focusing in space and time (blue line) and to a time-reversal experiment (orange line). The maximum value of the cross-correlation between the two signals is 0.99. (c) Real parts of the field at four times indicated by black dashed lines in (b) obtained from back-propagating the first left singular vector $d^*$. The field around the position of the target (black dot) could not be probed since the port is located outside the cavity. The input antennas (indicated by black arrows) are located 30 cm toward the left of the map.</figcaption></figure>

requires a single measurement of $S(\omega)$ at high power by measuring the norm of the asymmetric part

$$A = \|S^{\alpha} - (S^{\alpha})^{T}\|_{F}, \quad (7)$$

where $\|\cdot\|_{F}$ represents the Frobenius norm. At low power, $A(\omega)$ is dominated by the noise level as the system is operating in the linear regime (see Fig. 1(d)). However, at high power, $A(\omega)$ exhibits similar resonances as $\lambda_1(\omega)$. Note that this reciprocity condition is not completely able to extract $d$.

We now demonstrate spatio-temporal focusing on a nonlinear target. The absorbing foams shown in Fig. 1 are removed and the LNA is connected (without the resonator) to a single-mode waveguide located at the right interface. Apart from the nonlinear target, the cavity is now closed as the antennas at the right interface are terminated by open circuits (making them reflectors). The LNA is now powered so that the vector of transmission coefficients $t(\omega)$ to the target can also be measured by connecting the LNA to the eighth channel of the VNA. This is done only for the purpose of comparison as the measurement of $t(\omega)$ is not necessary to determine the incident waveform. The nonlinearity is nonresonant as the antennas are matched over a broad frequency range (see Fig. 2(a)). Fluctuations in $\lambda_1(\omega)$ arise from multiple scattering within the cavity.

Although the spatial wavefront for optimal focusing $c^{\text{in}}(\omega) = d^{*}(\omega)$ is obtained at each frequency through

an SVD, the frequency-dependence of the global phase $\phi(\omega)$ necessary for optimal temporal focusing is still unknown. We determine $\phi(\omega)$ using the procedure described in Ref. [34] that aligns the phase of each frequency component at the focal point (see Supplemental Material [41]). The temporal signals found from the inverse Fourier transform $s_{\text{opt}}(\tau) = \text{FT}^{-1}[c^{\text{opt}}(\omega) \cdot t(\omega)]$ and $s(\tau) = \text{FT}^{-1}[t^{*}(\omega) \cdot t(\omega)]$ are in excellent agreement with each other, as shown in Fig. 2(b), demonstrating maximal focusing both in time and space on the nonlinear target. The back-propagated signals within the cavity $\mathcal{I}(x, y; \tau) = \text{FT}^{-1}[c^{\text{in}*} \cdot t(x, y; \omega)]$ are presented in Fig. 2(c). Interestingly, the amplitude of the outgoing field at positive times is strongly reduced relative to the incident field at negative times in Fig. 2(c-d), which indicates strong absorption within the lossy nonlinear target.

In environments for which the energy density illuminating a target is small, wavefront shaping techniques may not be sufficient to detect the nonlinear signal. We thus investigate how the environment can be tuned to enhance the signal at the target [31, 47, 48]. For this purpose, we study a three-dimensional enclosure made programmable using a reconfigurable intelligent surface (RIS) (see Fig. 3(a)). For each of the 304 meta-atoms of the RIS, two states with a phase difference of roughly $\pi$ in reflection can be configured electronically. Seven antennas are used to measure a 7 × 7 scattering matrix in the spectral window between 4.8 GHz and 5.8 GHz.
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Schematic of the experiment. We measure the scattering matrix $S$ of a chaotic cavity using $N = 7$ antennas. (b) Optimization result for the first singular value $\lambda_1(f)$ of $\Delta S$ (blue line) compared with an average of 100 random configurations of the Reconfigurable Intelligent Surface (RIS) (orange line). (c) Same as (b) for the intensity measured at the target's position.</figcaption></figure>

The eighth antenna is connected to the nonlinear powered LNA.

The metasurface is first optimized iteratively to maximize the first singular value $\lambda_1(f)$ of $\Delta S$ at a single frequency $f_{opt} = 5.22$ GHz. This corresponds to a modification of the Green's function inside the system such that for the wavefront $d$ the intensity at the position of the nonlinear device is increased. The result is shown in Fig. 3 (b). At $f_{opt}$, the intensity enhancement is 2.3 fold compared to the average value for random configurations (Fig. 3 (c)).

**Conclusion.**— We have demonstrated that nonlinear elements embedded in a complex medium can be detected and localized using measurements of the scattering matrix at a single frequency for two incident powers (no harmonic generation is required). The monochromatic aspect of this technique enables to detect any type of nonlinearity, and is therefore particularly relevant in cases where the specific nonlinear response is unknown (n-th order harmonic generation, Kerr-like, etc.). We have shown that this noninvasive approach enables spatio-temporal focusing on nonlinear targets. The experiments presented here were limited to confined geometries but our technique is directly applicable to open systems. Since smartphones contain nonlinear elements, this could provide a way to enhance the focusing of Wi-Fi signals on these devices in complex environments and may inspire new detection and localization setups. Our approach is broadly applicable to any kind of wave and can readily be used in any domain in which waves are used to probe a medium (such as acoustics, optics, ...). In particular, it provides an efficient way to obtain the focusing wavefront

for wireless power transfer of bioelectronic devices [49].

This work is supported in part by the European Union through a European Regional Development Fund (ERDF), by the Ministry of Higher Education and Research, CNRS, Brittany region, Conseils Départementaux d’Ille-et-Vilaine and Côtes d’Armor, Rennes Métropole, and Lannion Trégor Communauté, through the CPER Project CyMoCod, in part by the French “Agence Nationale de la Recherche” (ANR) under Grant ANR-22-CPJ1-0070-01, and in part by the Direction Générale de l’Armement through the Creach Labs under the project TRsweep. M.D. acknowledges the Institut Universitaire de France. The metasurface prototypes were purchased from Greenerwave. The authors acknowledge P. E. Davy for the 3D rendering of the experimental setup in Fig. 3(a).

<footnote>* antton.goicoechea@univ-rennes.fr</footnote>

<footnote>† matthieu.davy@univ-rennes.fr</footnote>

<footnote>[1] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Controlling waves in space and time for imaging and focusing in complex media, Nature Photon 6, 283 (2012).</footnote>

<footnote>[2] S. Rotter and S. Gigan, Light fields in complex media: Mesoscopic scattering meets wave control, Rev. Mod. Phys. 89, 015005 (2017).</footnote>

<footnote>[3] H. Cao, A. P. Mosk, and S. Rotter, Shaping the propagation of light in complex media, Nat. Phys. 18, 994 (2022).</footnote>

<footnote>[4] I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, Exploiting disorder for perfect focusing, Nature Photon 4, 320 (2010).</footnote>

<footnote>[5] S. Jeong, Y.-R. Lee, W. Choi, S. Kang, J. H. Hong, J.-S. Park, Y.-S. Lim, H.-G. Park, and W. Choi, Focusing of light energy inside a scattering medium by controlling the time-gated multiple light scattering, Nature Photon 12, 277 (2018).</footnote>

<footnote>[6] N. Bender, A. Yamilov, A. Goetschy, H. Yılmaz, C. W. Hsu, and H. Cao, Depth-targeted energy delivery deep inside scattering media, Nat. Phys. 18, 309 (2022).</footnote>

<footnote>[7] Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, Optical phase conjugation for turbidity suppression in biological samples, Nature Photon 2, 110 (2008).</footnote>

<footnote>[8] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media, Phys. Rev. Lett. 104, 100601 (2010).</footnote>

<footnote>[9] X. Xu, H. Liu, and L. V. Wang, Time-reversed ultrasonically encoded optical focusing into scattering media, Nature Photon 5, 154 (2011).</footnote>

<footnote>[10] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, Correlation-enhanced control of wave focusing in disordered media, Nat. Phys. 10.1038/NPHYS4036 (2016).</footnote>

<footnote>[11] A. Derode, P. Roux, and M. Fink, Robust Acoustic Time Reversal with High-Order Multiple Scattering, Phys. Rev. Lett. 75, 4206 (1995).</footnote>

<footnote>[12] M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park,</footnote>
and W. Choi, Maximal energy transport through disordered media with the implementation of transmission eigenchannels, *Nature Photon* **6**, 581 (2012).

[13] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media, *Proc. Natl. Acad. Sci.* **117**, 14645 (2020).

[14] B. Judkewitz, Y. M. Wang, R. Horstmeyer, A. Mathy, and C. Yang, Speckle-scale focusing in the diffusive regime with time reversal of variance-encoded light (TROVE), *Nature Photon* **7**, 300 (2013).

[15] P. Lai, L. Wang, J. W. Tay, and L. V. Wang, Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media, *Nature Photon* **9**, 126 (2015).

[16] R. Horstmeyer, H. Ruan, and C. Yang, Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue, *Nature Photon* **9**, 563 (2015).

[17] D. Miller, Ultrasonic detection of resonant cavitation bubbles in a flow tube by their second-harmonic emissions, *Ultrasonics* **19**, 217 (1981).

[18] S. Qin, C. F. Caskey, and K. W. Ferrara, Ultrasound contrast microbubbles in imaging and therapy: Physical principles and engineering, *Phys. Med. Biol.* **54**, R27 (2009).

[19] C. L. Evans and X. S. Xie, Coherent Anti-Stokes Raman Scattering Microscopy: Chemical Imaging for Biology and Medicine, *Annual Rev. Anal. Chem.* **1**, 883 (2008).

[20] O. Katz, E. Small, Y. Guan, and Y. Silberberg, Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers, *Optica* **1**, 170 (2014).

[21] H. B. De Aguiar, S. Brasselet, and S. Gigan, Enhanced nonlinear imaging through scattering media using transmission-matrix-based wave-front shaping, *Phys. Rev. A* **94**, 043830 (2016).

[22] J. Moon, Y.-C. Cho, S. Kang, M. Jang, and W. Choi, Measuring the scattering tensor of a disordered nonlinear medium, *Nat. Phys.* **19**, 1709 (2023).

[23] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Deep optical imaging within complex scattering media, *Nat Rev Phys* **2**, 141 (2020).

[24] M. Reisner, D. H. Jeon, C. Schindler, H. Schomerus, F. Mortessagne, U. Kuhl, and T. Kottos, Self-Shielded Topological Receiver Protectors, *Phys. Rev. Applied* **13**, 034067 (2020).

[25] D. H. Jeon, M. Reisner, F. Mortessagne, T. Kottos, and U. Kuhl, Non-Hermitian C T -Symmetric Spectral Protection of Nonlinear Defect Modes, *Phys. Rev. Lett.* **125**, 113901 (2020).

[26] S. Suwunnarat, Y. Tang, M. Reisner, F. Mortessagne, U. Kuhl, and T. Kottos, Non-linear coherent perfect absorption in the proximity of exceptional points, *Commun Phys* **5**, 5 (2022).

[27] G. J. Mazzaro, A. F. Martone, K. I. Ranney, and R. M. Narayanan, Nonlinear Radar for Finding RF Electronics: System Design and Recent Advancements, *IEEE Trans. Microwave Theory Techn.* **65**, 1716 (2017).

[28] B. Perez, G. Mazzaro, T. J. Pierson, and D. Kotz, Detecting the Presence of Electronic Devices in Smart Homes Using Harmonic Radar Technology, *Remote Sensing* **14**, 327 (2022).

[29] P. Ambichl, A. Brandstötter, J. Böhm, M. Kühmayer, U. Kuhl, and S. Rotter, Focusing inside Disordered Media with the Generalized Wigner-Smith Operator, *Phys. Rev. Lett.* **119**, 033903 (2017).

[30] M. Horodynski, M. Kühmayer, A. Brandstötter, K. Pichler, Y. V. Fyodorov, U. Kuhl, and S. Rotter, Optimal wave fields for micromanipulation in complex scattering environments, *Nat. Photonics* **14**, 149 (2020).

[31] P. Del Hougne, K. B. Yeo, P. Besnier, and M. Davy, Coherent Wave Control in Complex Media with Arbitrary Wavefronts, *Phys. Rev. Lett.* **126**, 193903 (2021).

[32] D. Bouchet, S. Rotter, and A. P. Mosk, Maximum information states for coherent scattering measurements, *Nat. Phys.* **17**, 564 (2021).

[33] D. Bouchet and E. Bossy, Temporal shaping of wave fields for optimally precise measurements in scattering environments, *Phys. Rev. Research* **5**, 013144 (2023).

[34] K. B. Yeo, C. Leconte, P. del Hougne, P. Besnier, and M. Davy, Time Reversal Communications With Channel State Information Estimated From Impedance Modulation at the Receiver, *IEEE Access* **10**, 91119 (2022).

[35] J. Sol, L. L. Magoarou, and P. del Hougne, Optimal blind focusing on perturbation-inducing targets in sub-unitary complex media (2024), *arxiv:2401.15415 [physics]*.

[36] R. Berkovits, Sensitivity of the multiple-scattering speckle pattern to the motion of a single scatterer, *Phys. Rev. B* **43**, 8638 (1991).

[37] B. Orazbayev, M. Malléjac, N. Bachelard, S. Rotter, and R. Fleury, Wave-momentum shaping for moving objects in heterogeneous and dynamic media, *Nature Physics*, **1**, (2024).

[38] M. Horodynski, D. Bouchet, M. Kühmayer, and S. Rotter, Invariance Property of the Fisher Information in Scattering Media, *Phys. Rev. Lett.* **127**, 233201 (2021).

[39] J. Hüpfl, F. Russo, L. M. Rachbauer, D. Bouchet, J. Lu, U. Kuhl, and S. Rotter, Continuity equation for the flow of Fisher information in wave scattering, *Nat. Phys.* **20**, 1294 (2024).

[40] A. Fleming, C. Conti, and A. Di Falco, Perturbation of Transmission Matrices in Nonlinear Random Media, *Annual der Physik* **531**, 1900091 (2019).

[41] See Supplemental Material at [URL will be inserted by publisher], which includes Refs. [50–51], for the detailed theoretical analysis, additional simulations, and the phase correction method used for spatio-temporal focusing.

[42] M. Davy, M. Kühmayer, S. Gigan, and S. Rotter, Mean path length invariance in wave-scattering beyond the diffusive regime, *Commun Phys* **4**, 85 (2021).

[43] P. Del Hougne, R. Sobry, O. Legrand, F. Mortessagne, U. Kuhl, and M. Davy, Experimental Realization of Optimal Energy Storage in Resonators Embedded in Scattering Media, *Laser & Photonics Reviews* **15**, 2000335 (2021).

[44] D. L. Sounas and A. Alù, Fundamental bounds on the operation of Fano nonlinear isolators, *Phys. Rev. B* **97**, 115431 (2018).

[45] M. Cotrufo, S. A. Mann, H. Moussa, and A. Alu, Nonlinearity-Induced Nonreciprocity—Part I, *IEEE Trans. Microwave Theory Techn.* **69**, 3569 (2021).

[46] C.-Z. Wang, R. Kononchuk, U. Kuhl, and T. Kotzos, Loss-Induced Violation of the Fundamental Transmittance-Asymmetry Bound in Nonlinear Complex Wave Systems, *Phys. Rev. Lett.* **131**, 123801 (2023).

[47] P. Del Hougne, F. Lemoult, M. Fink, and G. Lerosey,
Spatiotemporal Wave Front Shaping in a Microwave Cav-
ity, Phys. Rev. Lett. 117, 134302 (2016).

[48] L. Chen, T. Kottos, and S. M. Anlage, Perfect absorption in complex scattering systems with or without hidden symmetries, Nat Commun 11, 5826 (2020).

[49] J. S. Ho, A. J. Yeh, E. Neofytou, S. Kim, Y. Tanabe, B. Patlolla, R. E. Beygui, and A. S. Y. Poon, Wireless power transfer to deep-tissue microimplants, Proc. Natl. Acad. Sci. U.S.A. 111, 7974 (2014).

[50] A. Yaghjian, Electric dyadic Green's functions in the

source region, Proceedings of the IEEE 68, 248 (1980).

[51] J. Hüpfl, N. Bachelard, M. Kaczvinszki, M. Horodynski, M. Kühmayer, and S. Rotter, Optimal cooling of multiple levitated particles: Theory of far-field wavefront shaping, Phys. Rev. A 107, 023112 (2023).

[52] Y. D. Chong and A. D. Stone, Hidden Black: Coherent Enhancement of Absorption in Strongly Scattering Media, Phys. Rev. Lett. 107, 163901 (2011).

Supplemental Material for “Detecting and Focusing on a Nonlinear Target in a Complex Medium”

I. THEORY

We consider a system without free charges and no magnetisation so that the wave equation is given by

$$
\nabla \times (\nabla \times \mathbf{E}) = -\mu_0 \partial_t^2 \mathbf{D}. \qquad (S1)
$$

We can now split up the electric displacement field **D** = ε₀ε**E** + **P**<sup>NL</sup> into a linear part given by the relative electric permittivity ε and a polarization part **P**<sup>NL</sup>[**E**], which is nonlinear in **E**. We are now going to make use of the Green's tensor G<sub>ω</sub> of the linear system at frequency ω with outgoing boundary conditions, i.e.

$$
\nabla \times (\nabla \times G_{\omega})(\mathbf{r}, \mathbf{r}') - \frac{\omega^2}{c^2} \epsilon(\mathbf{r}) G_{\omega}(\mathbf{r}, \mathbf{r}') = 1\delta(\mathbf{r} - \mathbf{r}'). \quad (S2)
$$

This definition allows us to rewrite the monochromatic solution of the wave equation at frequency ω using the Green’s tensor in the far field

$$
E_{\omega}(\mathbf{r}) = E_{\omega}^{\text{L}}(\mathbf{r}) + \frac{k^2}{\epsilon_0} \int G_{\omega}(\mathbf{r}, \mathbf{r}') P_{\omega}^{\text{NL}}[E(\mathbf{r}')] d\mathbf{r}', \quad (S3)
$$

where $k = \omega/c$, $E_{\omega}^{L}$ describes the incident field component at frequency $\omega$ scattered in the linear system (i.e. $\nabla \times (\nabla \times E_{\omega}^{L}) + k^2 \epsilon E_{\omega}^{L} = 0$) and $P_{\omega}^{NL}(E)$ the nonlinear polarization field at frequency $\omega$ corresponding to the full electric field $E$. We want to emphasize here that $P_{\omega}^{NL}(E)$ depends in general on the full electromagnetic field $E$ at all frequencies and not just the field component $E_{\omega}$ at frequency $\omega$. This occurs in cases such as second harmonic generation, where the polarization field $P_{2\omega}$ at $2\omega$ depends on the electric field component $E_{\omega}$ at $\omega$.

A. Separating the linear and nonlinear contributions

Our first goal is going to be to separate these two contributions. For this we can make us of the fact that the linear scaling of the output field with respect to the input field is broken due to the presence of the nonlinearity. We proceed by varying the strength (quantified by α) of the incident field $E_{\omega}^{L,in} \rightarrow E_{\omega}^{L,in} \alpha$, which produces a predictable linear response in the linear contributions of the scattered field $E_{\omega}^{L,\alpha} = E_{\omega}^{L}\alpha$. Therefore by varying α we can extract the nonlinear contributions of the field by removing exactly those contributions with linear scaling

$$
E_{\omega}^{\alpha+\Delta\alpha} - \frac{\alpha + \Delta\alpha}{\alpha} E_{\omega}^{\alpha} = \frac{k^2}{\epsilon_0} \int G_{\omega}(\mathbf{r}, \mathbf{r}') \left\{ P_{\omega}^{\text{NL}}[E^{\alpha+\Delta\alpha}(\mathbf{r}')] - \frac{\alpha + \Delta\alpha}{\alpha} P_{\omega}^{\text{NL}}[E^{\alpha}(\mathbf{r}')] \right\} d\mathbf{r}' \quad (S4)
$$

whereby the superscript α, α + Δα quantifies the strength of the incident field. The left hand side of the equation now effectively corresponds to a field that originates at the nonlinearity and propagates out of the system according to the Green’s tensor of the corresponding linear system G<sub>ω</sub>(r, r'). This already shows that a nonlinearity can be detected just by observing the divergence of the outgoing field from the linear field strength scaling property. In order to rewrite this expression to Eq. (2) of the main text we divide by α + Δα and define the difference of the normalized fields δE<sub>ω</sub>(r) = (α + Δα)<sup>-1</sup>E<sub>ω</sub><sup>α+Δα</sup>(r) - α<sup>-1</sup>E<sub>ω</sub><sup>α</sup>(r), which satisfies

$$
\delta \mathbf{E}_{\omega}(\mathbf{r}) = \frac{k^2}{\epsilon_0} \int d\mathbf{r}' G_{\omega}(\mathbf{r}, \mathbf{r}') \delta \mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}'), \quad (S5)
$$
with $\delta \mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}') = (\alpha + \Delta\alpha)^{-1} \mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha+\Delta\alpha}(\mathbf{r}')] - \alpha^{-1} \mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}^{\alpha}(\mathbf{r}')]$. In the next section we are going to see how this can be exploited to create a focusing field on a point-like target.

## B. Focusing on the nonlinearity

Our next objective will be to find an incident wavefront that focuses onto the location of the nonlinearity. We reverse the outgoing finite difference field (assumed non-zero from here on) derived in the previous section so that the incident field is given by $\mathbf{E}_{\omega}^{\text{opt}} = \delta \mathbf{E}_{\omega}^{*}$ in the far field and creates a focus at the nonlinear dielectric. This can be seen for systems, where $\epsilon$ is real valued (i.e. no absorption or gain present) by using Green's identity to identify the linear response $\mathbf{E}_{\omega}^{\text{L, opt}}$ of the system to the incident field $\delta \mathbf{E}_{\omega}^{*}$, i.e. a system with the nonlinearity removed,

$$ \mathbf{E}_{\omega}^{\text{L, opt}}(\mathbf{r}_1) = \int_{S_R^3} [\mathbf{n} \times G_{\omega}(\mathbf{r}, \mathbf{r}_1)] \cdot [\nabla \times \mathbf{E}_{\omega}^{\text{L, opt}}] - [\nabla \times G_{\omega}(\mathbf{r}, \mathbf{r}_1)] \cdot [\mathbf{n} \times \mathbf{E}_{\omega}^{\text{L, opt}}] d\sigma, \quad (S6) $$

where $S_R^3$ is the surface of a sphere with large radius $R$ centered at the origin. Due to the Green's function having outgoing boundary conditions (Silver-Müller radiation condition), the outgoing part of $\mathbf{E}_{\omega}^{\text{L, opt}}$ does not contribute to this equation and we thus get

$$ \mathbf{E}_{\omega}^{\text{L, opt}}(\mathbf{r}_1) = \int_{S_R^3} [\mathbf{n} \times G_{\omega}(\mathbf{r}, \mathbf{r}_1)] \cdot [\nabla \times \delta \mathbf{E}_{\omega}^{*}(\mathbf{r})] - [\nabla \times G_{\omega}(\mathbf{r}, \mathbf{r}_1)] \cdot [\mathbf{n} \times \delta \mathbf{E}_{\omega}^{*}(\mathbf{r})] d\sigma. \quad (S7) $$

We can simplify this expression by using the Green's function identity

$$ \begin{aligned} -2i\operatorname{Im}(G_{\omega}(\mathbf{r}_0, \mathbf{r}_1)) &= -2i \int_{S_R^3} (\nabla \times \operatorname{Im}(G_{\omega})(\mathbf{r}, \mathbf{r}_1)) \cdot (\mathbf{n} \times G_{\omega}(\mathbf{r}, \mathbf{r}_0)) - (\mathbf{n} \times \operatorname{Im}(G_{\omega})(\mathbf{r}, \mathbf{r}_1)) \cdot (\nabla \times G_{\omega}(\mathbf{r}, \mathbf{r}_0)) d\sigma \\ &= \int_{S_R^3} (\nabla \times G_{\omega}^{*}(\mathbf{r}, \mathbf{r}_1)) \cdot (\mathbf{n} \times G_{\omega}(\mathbf{r}, \mathbf{r}_0)) - (\mathbf{n} \times G_{\omega}^{*}(\mathbf{r}, \mathbf{r}_1)) \cdot (\nabla \times G_{\omega}(\mathbf{r}, \mathbf{r}_0)) d\sigma. \end{aligned} \quad (S8) $$

The second equality holds due to the outgoing boundary condition of $G_{\omega}$. Note that we avoided the ambiguity of the dyadic Green's function in the source region [S50] in Eq. (S7)/(S8) because $\mathbf{E}_{\omega}^{\text{L, opt}}$, $\operatorname{Im}(G_{\omega})$ are solutions of the source free wave equation inside the system.

By combining Eq. (S5) with (S7) and (S8) this gives us the linear contribution of the field at the nonlinearities

$$ \mathbf{E}_{\omega}^{\text{L, opt}}(\mathbf{r}) = -\frac{2ik^2}{\epsilon_0} \int \operatorname{Im} G_{\omega}(\mathbf{r}, \mathbf{r}') (\delta \mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}'))^* d\mathbf{r}'. \quad (S9) $$

Since we assumed that the probing field produces a nonlinear response, i.e., $\delta \mathbf{E}_{\omega} \neq 0$, it follows directly from the definition that $\delta \mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}', \mathbf{r})$ must also be non-zero. Eq. (S9) therefore implies that the electric field within the nonlinear region $\mathbf{E}_{\omega}^{\text{L, opt}}(\mathbf{r})$ can be expected to be non-zero.

Overall for general nonlinearities (e.g. multiple nonlinearities, extended nonlinearities in space) our method is able to identify incident states that focus at the spatial region, where these non-linearities are located. However, both the spatial correlations of the electromagnetic field, captured by $\operatorname{Im} G_{\omega}(\mathbf{r}, \mathbf{r}', \mathbf{r})$, and the nonlinear response of the target, $\delta \mathbf{P}_{\omega}^{\text{NL}}$, obfuscate the exact field intensity at the target. This highlights the need for additional system information, such as details of the nonlinearities in the system, in order to find the optimal input wavefront for focusing. In our case we overcome this by considering systems with only a single input mode that couples to the nonlinearity such as in the case of a nonlinearity made up of an antenna measuring only a singular polarization direction connected to a LNA or a point-like nonlinearity in restricted systems with only a singular polarization degree of freedom.

## C. Point-like nonlinearity

We consider a point-like nonlinearity in the Rayleigh regime with center at $\mathbf{r}_0$ volume $V_R$ and diameter $R$, so that the far field is given by

$$ \mathbf{E}_{\omega}(\mathbf{r}) = \mathbf{E}_{\omega}^{\text{L}}(\mathbf{r}) + \frac{k^2 V_R}{\epsilon_0} G_{\omega}(\mathbf{r}, \mathbf{r}_0) \mathbf{P}_{\omega}^{\text{NL}}[\mathbf{E}(\mathbf{r}_0)]. \quad (S10) $$
Thus the outgoing field pattern is only given by the Green's tensor of the linear system $G_\omega(\mathbf{r}, \mathbf{r}_0)$ coupling the location of the nonlinearity to the far field. Based on this we can simplify Eq. (S9) giving us the linear contribution

$$ \mathbf{E}_{\omega}^{\text{L, opt}}(\mathbf{r}_0) = -\operatorname{Im}G_{\omega}(\mathbf{r}_0, \mathbf{r}_0) \frac{2ik^2V_R}{\epsilon_0} (\delta\mathbf{P}_{\omega}^{\text{NL}}(\mathbf{r}_0))^*. \quad (S11) $$

Thus in the linear reference system a focus can be created at the location of the nonlinearity. In the case, where $\mathbf{P}_{\omega}^{\text{NL}}$ is restricted to a singular polarization direction (e.g. antenna) then only a singular input wavefront can focus onto our target and the focus is therefore optimal.

## II. FAR FIELD DESCRIPTION

We consider the case of waves at the frequency $\omega$ and a basis of the wavefronts in the far field $\mathbf{E}_{\omega}^n$ so that the in- and outgoing fields have the corresponding coefficients $\mathbf{c}^{\text{in}}, \mathbf{c}^{\text{out}}$ (see section II A). We will for now focus on the case, where the nonlinearity only couples to the system through a single polarization degree of freedom (e.g. an antenna or a system allowing only a single polarization direction), while the more general case is considered in section II E. In linear systems these coefficients can be connected by the linear scattering matrix $S^L$. However in our case this does no longer hold true due to the nonlinearity violating the superposition principle, where the connection $\mathbf{c}^{\text{out}} = \hat{S}(\mathbf{c}^{\text{in}})$ is given by the non-linear scattering operator $\hat{S}$. Luckily, those incident fields that do not interact with the nonlinearity stay linear, which means that for a point-like nonlinearity an incident field quantified by $\mathbf{d}^*$ can be found, so that all orthogonal incident fields do not interact with the nonlinearity, i.e.

$$ \hat{S}(\mathbf{c}^{\text{in}}) - S^{\text{L}}\mathbf{c}^{\text{in}} = 0 \quad (S12) $$

whereby $S^L$ is the scattering matrix of the linear system without the nonlinearity. This can be seen by considering that the vector $\mathbf{d}$ corresponds to the far field coefficients of $G_\omega(\mathbf{r}, \mathbf{r}_0)$ in our chosen basis, i.e. $d_j = \frac{-i}{2\omega\mu_0} \int_{S_R^2} (\mathbf{n} \times (\mathbf{E}_\omega^j(r))) \cdot (\nabla \times G_\omega(r, r_0)) - (\nabla \times (\mathbf{E}_\omega^j(r))) \cdot (\mathbf{n} \times G_\omega(r, r_0)) d\sigma$ (see section II A). Using Green's identity we can see that $d_j = \frac{i}{2\omega\mu_0}(\mathbf{E}_\omega^{L,j})(r_0)$, where $\mathbf{E}_\omega^{L,j}$ corresponds to the solution of the linear system of the far field mode $j$. Thus if we take an incident field in the input channels $\mathbf{c}^{\text{in}} \cdot \mathbf{d} = 0$ then the linear part of the field is given by $\mathbf{E}_\omega^L(r_0) = \sum_j c_j \mathbf{E}_\omega^{L,j}(r_0) \propto \mathbf{d} \cdot \mathbf{c} = 0$, showing that no field is present at the location of the nonlinearity. Due to the non-linear polarization $P^{\text{NL}}$ being only non-zero for electromagnetic fields that interact with the non-linearity ($E(r_0) \neq 0$), this means that only the components of $\mathbf{c}^{\text{in}}$ parallel to $\mathbf{d}^*$ result in differences between $\hat{S}$ and $S^L$, i.e.

$$ \hat{S}(\mathbf{c}^{\text{in}}) - S^{\text{L}}\mathbf{c}^{\text{in}} = \hat{S}(\hat{P}_{d^*}\mathbf{c}^{\text{in}}) - S^{\text{L}}\hat{P}_{d^*}\mathbf{c}^{\text{in}}, \quad (S13) $$

where $\hat{P}_{d^*} = \mathbf{d}^*\mathbf{d}^T/|\mathbf{d}|^2$ is the projection operator onto the vector $\mathbf{d}^*$. We can simplify this even further by noting that Eq. (S10) tells us that the part of the outgoing field due to the non-linearity results in an outgoing field given by $G_\omega(\mathbf{r}, \mathbf{r}_0)$, which is described by $\mathbf{d}$ so that we get

$$ \hat{S}(\mathbf{c}^{\text{in}}) - S^{\text{L}}\mathbf{c}^{\text{in}} = \hat{P}_{d} \left( \hat{S}(\hat{P}_{d^*}\mathbf{c}^{\text{in}}) - S^{\text{L}}\hat{P}_{d^*}\mathbf{c}^{\text{in}} \right), \quad (S14) $$

where $\hat{P}_d$ is the projection operator onto the vector $\mathbf{d}$. By inserting the definition of the orthogonal projection operators we can now describe the scattering of the nonlinear scattering operator by

$$ \hat{S}(\mathbf{c}^{\text{in}}) = S^{\text{L}}\mathbf{c}^{\text{in}} + \frac{\mathbf{d}\mathbf{d}^{\dagger}}{|\mathbf{d}|^2} \left( \hat{S} \left( \frac{\mathbf{d}^{*}\mathbf{d}^{T}\mathbf{c}^{\text{in}}}{|\mathbf{d}|^2} \right) - S^{\text{L}}\frac{\mathbf{d}^{*}\mathbf{d}^{T}\mathbf{c}^{\text{in}}}{|\mathbf{d}|^2} \right) = S^{\text{L}}\mathbf{c}^{\text{in}} + \mathbf{d}f(\mathbf{d} \cdot \mathbf{c}^{\text{in}}), \quad (S15) $$

whereby $S^L$ is the scattering matrix of the linear system without the nonlinearity and $f$ is the nonlinear scalar function describing the deviation of $\hat{S}$ from a linear relation, i.e. $f(x) = \mathbf{d}^\dagger[\hat{S}(\mathbf{d}^*x/|\mathbf{d}|^2) - S^L\mathbf{d}^*x/|\mathbf{d}|^2]/|\mathbf{d}|^2$. This nonlinear operator can now be probed using a basis of incident channels, which can be conveniently summarized in the matrix $S^\alpha$ so that we have

$$ S_{m,n}^{\alpha} = S_{m,n}^{\text{L}} + d_m \alpha^{-1} f(\alpha d_n), \quad (S16) $$

where we choose the basis of incident channels $\mathbf{c}^{\text{in}} = \alpha \mathbf{e}_n$ for all $n$, quantify the strength of the incident fields by $\alpha$ and normalize the outgoing fields by $\alpha$. This has the advantage that for linear system this reduces to the scattering matrix, while we can probe the nonlinear contributions of the scattering operator by varying $\alpha$. Note however that $S^\alpha$ does not fully describe the operator $\hat{S}$, due to the loss of the superposition principle in nonlinear systems, but only the response of the system for a basis of incident channels.
A. Connecting the far field patterns with the scattering matrix

We connect the scattering matrix with the far field pattern using the hermitian form based on the poynting vector [S51]

$$ \langle E_{\omega}^{1}, E_{\omega}^{2} \rangle = \frac{1}{2} \int_{S_{R}^{3}} [(E_{\omega}^{1})^{*} \times H_{\omega}^{2} - (H_{\omega}^{1})^{*} \times E_{\omega}^{2}] \cdot dn, \qquad (S17) $$

where $S_R^3$ is the surface of a sphere with large radius $R$ centered at the origin. We can now use this product to define a basis of incident states $E_\omega^j$ with power normalization, i.e.

$$ \langle E_{\omega}^{i}, E_{\omega}^{j} \rangle = -\delta_{i,j}. \qquad (S18) $$

We now define the input coefficients $c^{\text{in}}$ corresponding to a field $E_\omega$ by

$$ c_{i}^{\text{in}} = -\langle E_{\omega}^{i}, E_{\omega} \rangle. \qquad (S19) $$

Similarly the outgoing field if given by $(E_\omega^j)^*$, which can be used to define the outgoing coefficients $c^{\text{out}}$ by

$$ c_{i}^{\text{out}} = \langle (E_{\omega}^{i})^{*}, E_{\omega} \rangle. \qquad (S20) $$

The scattering operator is now given by the relation $c^{\text{out}} = \hat{S}[c^{\text{in}}]$, which is linear in the case where no non-linearity is present. Finally, we can define the outgoing coefficients $d$ corresponding to the Green's function (e.g. $E_\omega(r) = G_\omega(r, r_0)$, $H_\omega = \nabla \times E_\omega / (i\mu_0\omega)$) by

$$ d_i = \langle (E_{\omega}^{i})^{*}, G_{\omega}(\cdot, r_0) \rangle = \frac{-i}{2\omega\mu_0} \int_{S_R^3} (\mathbf{n} \times (\mathbf{E}_{\omega}^j)(r)) \cdot (\nabla \times G_{\omega}(r, r_0)) - (\nabla \times (\mathbf{E}_{\omega}^j)(r)) \cdot (\mathbf{n} \times G_{\omega}(r, r_0)) d\sigma. \qquad (S21) $$

B. Extracting the focusing field

We will now use Eq. (S16) at two incident powers $\alpha, \alpha + \Delta\alpha$ to define a matrix version of Eq. (S5). The results are summarized by the matrices $S^\alpha, S^{\alpha+\Delta\alpha}$, which gives us

$$ \Delta S_{m,n} = (S^{\alpha+\Delta\alpha} - S^{\alpha})_{m,n} = d_m ((\alpha + \Delta\alpha)^{-1} f((\alpha + \Delta\alpha)d_n) - \alpha^{-1} f(\alpha d_n)). \qquad (S22) $$

This rank-one matrix contains information of at least one incident channel that interacts with the nonlinearity, allowing us to apply a singular value decomposition on $\Delta S$ to extract $d$. Note that if multiple far field modes couple to the nonlinearity (e.g. non-scalar waves, multiple polarization degrees of freedom) or if multiple nonlinearities are present then the non-linearity acts on the subspace of these modes. In this case $f$ would need to be replaced by a function describing the nonlinear interaction between these modes and the rank of the matrix $\Delta S$ will be the dimension of this subspace (see section II E). In general while we can still use $\Delta S$ to identify modes that focus on the nonlinearity, in order to create an optimal focus more information on the non-linearity is needed.

C. Reciprocity and nonlinearity detection

One important property in most linear systems is the reciprocity condition, which can now break due to nonlinear interactions. In our case this turns out to be

$$ (S^{\alpha} - (S^{\alpha})^{T})_{m,n} = d_{m}\alpha^{-1}f(\alpha d_{n}) - d_{n}\alpha^{-1}f(\alpha d_{m}). \qquad (S23) $$

While in theory it is not guaranteed that we will see reciprocity breaking using $S^\alpha$ (e.g. if $d$ corresponds to a basis vector $d_n \propto \delta_{m,n}$), in practice we saw that this can serve as a useful tool for the detection of the nonlinearity.
## D. Lowest singular value and absorption

We find a strong average correlation between the wavefront for maximal focusing $c^{\text{in}} = d^*$ and the eigenvector $U_N$ corresponding to the smallest eigenvalue $\sigma_N$ of $(S^\alpha)^\dagger S^\alpha$ that gives minimal reflection from the cavity and therefore maximal absorption [S52]. The correlation coefficient averaged over the frequency range is $\langle |U_N^* \cdot c^{\text{in}}| \rangle \sim 0.8$. As the cavity is closed, absorption at the target is indeed the main loss mechanism. The correlation coefficient is however below unity since other dissipative mechanisms such as uniform absorption within the cavity also takes place. Small eigenvalues $\sigma_N \to 0$ indicate that the incident energy is almost completely dissipated within the target.

## E. Extended and multiple nonlinearities

We will now consider extended and multiple nonlinearities $D_i$. Similar to section II the far field modes that couple to the regions of these nonlinearities are spanned by $G_\omega(\mathbf{r}, \mathbf{r}_0)$ for $\mathbf{r}_0 \in D = \cup_i D_i$, which we will describe in the far field by the basis $\mathbf{d}^1, \mathbf{d}^2, \dots, \mathbf{d}^N$ of dimension $N$. With the same arguments as in section II we can now show that for an incident field quantified by $c$ and $c \cdot d^i = 0$ for all $i$, the corresponding electric field in the linear system disappears at the locations of the nonlinearities, i.e. $E_\omega^L(r_0) = 0$ for $\mathbf{r}_0 \in D$. As such the vectors $(d^i)^*$ span the space of input modes that can couple to the nonlinearity and $d^i$ the outgoing modes that the nonlinearity couples to. Using this we now write the scattering operator

$$ \hat{S}[c^{\text{in}}] = S^{\text{L}}c^{\text{in}} + \sum_{i=1}^{N} d^i f_i(\mathbf{d}^1 \cdot c^{\text{in}}, \dots, \mathbf{d}^N \cdot c^{\text{in}}), \quad (S24) $$

for the non-linear functions $f_i$. Next we consider the difference matrix given by

$$ \Delta S_{m,n} = \sum_{i} d_{m}^{i} \left( (\alpha + \Delta\alpha)^{-1} f_{i}((\alpha + \Delta\alpha)d_{n}^{1}, \dots, (\alpha + \Delta\alpha)d_{n}^{N}) - \alpha^{-1} f_{i}(\alpha d_{n}^{1}, \dots, \alpha d_{n}^{N}) \right). \quad (S25) $$

We can see that the dimension of $\Delta S$ is at most the number of in coupling modes $N$ and that the space of left singular values is spanned by $d^i$. This shows that we can use $\Delta S$ to extract the incident wavefront that focus on the nonlinearities. However it is important to note that in general the vectors $d^i$ will not be the set of left singular vectors, due to $((\alpha + \Delta\alpha)^{-1} f_i((\alpha + \Delta\alpha)d_n^1, \dots, (\alpha + \Delta\alpha)d_n^N)) - \alpha^{-1} f_i(\alpha d_n^1, \dots, \alpha d_n^N))$ being in general not orthogonal for different $i$.

## III. SIMULATIONS FOR CONTROLLED NONLINEAR BEHAVIOR

We present the results of numerical simulations obtained with COMSOL. The goal here is to demonstrate the versatility of the method with respect to the type of nonlinearity present in the system. We consider a 2D cavity with 9 antennas on each side at a frequency $f = 13.51$ GHz. The field inside the cavity is randomized through the presence of 14 randomly located metallic scatterers of 1 mm radius. The nonlinear target is defined as a dielectric scatterer with a 1 mm radius and Kerr-like index given by:

$$ n = \begin{cases} n_0 + n_2 |\mathbf{E}_\omega|^2 & \text{case 1} \\ n_0 - in_2 |\mathbf{E}_\omega|^2 & \text{case 2,} \end{cases} \quad (S26) $$

where the case 1 corresponds to nonlinear dispersion, and the case 2 corresponds to nonlinear absorption. The intensity maps obtained from the first left singular vectors of $\Delta S$ in each cases (as described in the main text) are shown in Fig. S1. The same locations are used for the scatterers and the target inside the system, and for the value of the nonlinear response $n_2 = 5 \cdot 10^{-6} \text{ m}^2/\text{V}^2$. In both cases the incident wavefront obtained from the SVD focuses on the target regardless of the chosen nonlinear function, showing that the approach is independent of the type of nonlinearity.

## IV. CORRECTION OF THE PHASE FOR BROADBAND SIGNALS

The wavefront for optimal focusing in space and time corresponds to the time-reversed (or equivalently phase-conjugate) of the transmission coefficient $c^{\text{opt}}(\omega) = t^*(\omega)$. The phases at each frequency are aligned at the focus,
<figure><img src="image_5.png" /><figcaption>FIG. S1. (a, b) Intensity map within the system for $c^{\text{in}} = d^*$ for a nonlinearity which intensity-dependent index is given by Eq. (S26). Panel (a) corresponds to case 1 and panel (b) to case 2. Both maps are normalized by the maximum value obtained for the optimized wavefront. The white dots represent the location of the metallic scatterers and the pink circle the location of the nonlinear target.</figcaption></figure>

meaning that the scattered wavefront $c^{\text{opt}}(\omega) = \Delta S(\omega)c^{\text{opt}}(\omega)$ acquires a phase equal to $\arg[t(\omega)]$. We therefore determine the phase $\phi(\omega)$ from the condition $\arg[c^{\text{in}}e^{i\phi}\Delta S c^{\text{in}}e^{i\phi}] = 0$. Because both $c^{\text{in}}e^{i\phi}$ and $c^{\text{in}}e^{i(\phi+\pi)}$ satisfy this condition, we finally exploit the continuity of $\phi(\omega)$ over the bandwidth to correct $\pi$-phase shifts.